A state-space model to quantify common input to motor neurons (Feeney et al 2017)

 Download zip file 
Help downloading and running models
Accession:238912
"... We introduce a space-state model in which the discharge activity of motor neurons is modeled as inhomogeneous Poisson processes and propose a method to quantify an abstract latent trajectory that represents the common input received by motor neurons. The approach also approximates the variation in synaptic noise in the common input signal. The model is validated with four data sets: a simulation of 120 motor units, a pair of integrate-and-fire neurons with a Renshaw cell providing inhibitory feedback, the discharge activity of 10 integrate-and-fire neurons, and the discharge times of concurrently active motor units during an isometric voluntary contraction. The simulations revealed that a latent state-space model is able to quantify the trajectory and variability of the common input signal across all four conditions. When compared with the cumulative spike train method of characterizing common input, the state-space approach was more sensitive to the details of the common input current and was less influenced by the duration of the signal. The state-space approach appears to be capable of detecting rather modest changes in common input signals across conditions."
Reference:
1 . Feeney DF, Meyer FG, Noone N, Enoka RM (2017) A latent low-dimensional common input drives a pool of motor neurons: a probabilistic latent state-space model. J Neurophysiol 118:2238-2250 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Synapse;
Brain Region(s)/Organism:
Cell Type(s): Spinal cord renshaw cell; Abstract integrate-and-fire leaky neuron;
Channel(s):
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: MATLAB;
Model Concept(s): Action Potentials; Synaptic noise;
Implementer(s):
%% The most realistic set of data below
%As of 3/19, this creates a set of N neurons with a latent
%state . As a result, A is 2x2 identity with some noise and multiplied through
%the 500 iterations. C is Nx2 and encodes the effect of the latent state on
%each neuron (each row). C was estimated to create realistic discharge rates. 
%Xs is how X changes through the trial and includes
%Gaussian noise awgn. 
%z(t) = Cx(t). Discharge rate is exp(z(t)).
%Start by estimating the number of active motor units given a desired force
%Output of the model is stored in seq, which is run through PLDSExample

clear all
%% Make force trace 
GoalF = 5;
force(1:101) = 0:GoalF/100:GoalF;                      %Linear 1s ramp with 4 second hold of a target force
force(101:500) = force(101);

%% Parameters of model: X and A
Xs = zeros(500,2);
A = [1,0; 0, 1];                        %This is the linear approximation of the dynamics and should always be dxd where d is dimensionality of latent state
Xi = [1 ; 1];                           %This is the initial latent state
Xs(1,:) = A * Xi;                       %How the latent state evolves through the trial

%% Estimate the numbers of active motor units for the 5 second trial. Based on Fuglevand, 1993 on the excitation level a unit begins
%Discharging action potentials. The force in this case is a proxy for the
%synaptic current

RTE = zeros(1,120); %Preallocating the vector of recruitment excitation levels. Sorry for the for loop
a = log(60)/120;    %Parameters of exponential relation. Last unit recruited at 60% MVC
for ff = 1:120
    RTE(ff) = exp(a*ff);
end

%Need to calcualte, at each instant, how far above recruitment threshold a
%unit is
AboveRTE = zeros(120,500);                  %Preallocate matrix where each row is a neuron and length of 500 ms.
for tt = 1:500
    AboveRTE(:,tt) = force(tt) - RTE;          %Find the difference between recruitment threshold and current force output
end

ActiveUnits =zeros(1,120);          %Preallocate how many units will eventually be 'turned on' by this force

for yy = 1:500
    nums = find(AboveRTE(:,yy)>0);                %Find all entries above zero
    ActiveUnits(yy) = length(nums);     %Find the index where RTE is greater than the force.This is approx # motor units active 
end
AboveRTE(AboveRTE<0)=0; %If units are not above RTE during the trial, they should be disregarded
neurons = ActiveUnits(500); %Total number of units discharging action potentials during the trial

%% Set the yDim to the number of motor activated (from the previous step)
Trials = 1; yDim = ActiveUnits(500); T = 500;
seq = struct();

%% The following will need to adjust the firing rates for the N active motor units using x, A, and u
for i = 2:500
    Xs(i,:) = A * Xs(i-1,:)';  %Evolution of latent state. Currently a boring identity matrix
end
Xs = awgn(Xs,35);                %Adding Gaussian noise to evolution, which models synaptic fluctations 
 
C = zeros(yDim,2);       %making driving matrix C be n rows (one for each neuron) and two columns (one for each latent input)
for mm = 1:yDim
    C(mm,1) = 2 + 0.00541*exp((0.049)*mm);        %Set the linear equation for how C changes for each neuron for input conductance
    C(mm,2) = 0.5 + 0.001666 * exp((0.048)*mm);     %Set linear equaiton for C for AHP duration
end


LinMult = AboveRTE ./5;                         %Create column array of force differences to model sigmoid relation between input force and discharge rate
LinMult = tanh(LinMult) * 0.5 + 0.5;             
LinMult = LinMult(1:ActiveUnits(500),:);            %Get to correct size

FiringRate = zeros(ActiveUnits(500),500);           %preallocate
AboveRTE = AboveRTE(1:ActiveUnits(500),1:500);
NoiseMat = ones(ActiveUnits(500),1); 
NoiseMat = NoiseMat .* (AboveRTE(:,500).* 5);
NoiseMat(NoiseMat<0) = 0;


%% Use the noise plut the discharge rates and send through exponential.
for kk = 1:500
    FiringRate(:,kk) = (C * Xs(kk,:)').* LinMult(:,kk); %Average firing rate for each neuron in column array 
end

for ooo = 1:ActiveUnits(500)
    FiringRate(ooo,:) =  awgn(FiringRate(ooo,:), NoiseMat(ooo)); %Add noise to the firing rate
end

FiringRate = exp(FiringRate); %Send through exponential

MultFR = FiringRate;

for rr = 1:ActiveUnits(500)
    for ss = 1:length(AboveRTE)
        if AboveRTE(rr,ss) == 0;
            MultFR(rr,ss) = 0;
        end
    end
end
AvgFiringRate = mean(MultFR(:,300:500),2);
%% Find the lambda value for PoissonRnd.
rng(123);
DummyVariable = zeros(ActiveUnits(500),5); %This is going to be the 5 seconds drawn from the Poisson Distribution
for JA = 1:ActiveUnits(500)
    DummyVariable(JA,:) = poissrnd(AvgFiringRate(JA),1,5);
end

Try = zeros(67,500);
for dan = 1:ActiveUnits(500)
   Try(dan,1:100) = binornd(ones(1,100),DummyVariable(dan,1)/100); 
   Try(dan,101:200) = binornd(ones(1,100),DummyVariable(dan,2)/100);
   Try(dan,201:300) = binornd(ones(1,100),DummyVariable(dan,3)/100);
   Try(dan,301:400) = binornd(ones(1,100),DummyVariable(dan,4)/100);
   Try(dan,401:500) = binornd(ones(1,100),DummyVariable(dan,2)/100);
end


seq.y = Try;
seq.T = T;
seq.yr = MultFR;
seq.x = Xs;
%seq.u = ones(2,500).*0.5; Unsure the best inputs for u

NumPlots = floor(ActiveUnits(500)/10); 
figure(1)
for aa = 1:NumPlots
    subplot(3,5,aa);
    plot(1:500, MultFR(aa*10,:),'*');
    xlabel('bin');
    ylabel('Firing Rate Hz');
    axis([0,500,0,20]);
end

%% print DR and CV values for time 0.5 - 1 s and 3.5 to 4 s
Multiplier = floor(ActiveUnits(500)/10);

ISI = 1./MultFR;  %Calculate the interspike interval at each time instant
ISI2 = ISI;
ISI(ISI==inf)=0;
ISI2(ISI2==inf)=0;
ISI2 = 10^3 .*ISI2;               %Convert to ms from seconds

%%% Time vector voodoo %%%
rounded_ISI = floor(ISI2);
rounded_ISI(isnan(rounded_ISI)) = 0; 
rounded2 = rounded_ISI(:,25:end);
for dd = 1:length(AboveRTE(:,500))
    discharge_index(dd,:) = cumsum(rounded_ISI(dd,:));
end


SDlow = zeros(1,Multiplier);
AVlow = zeros(1,Multiplier);
CVlow = zeros(1,Multiplier);
SDhi = zeros(1,Multiplier);
AVhi = zeros(1,Multiplier);
CVhi = zeros(1,Multiplier);

SDlow(1) = (std(ISI(1,50:100)))^2;          %Calculating SD^2 for ISI
AVlow(1) = (mean(ISI(1,50:100)))^3;         %Calculating avg^3 for ISI
CVlow(1) = sqrt((SDlow(1)/AVlow(1)))*100;
SDhi(1) = (std(ISI(1,350:400)))^2;
AVhi(1) = (mean(ISI(1,350:400)))^3;
CVhi(1) = sqrt((SDhi(1)/AVhi(1)))*100;
DRlow(1) = mean(MultFR(1,50:100));
DRhi(1) = mean(MultFR(1,350:400));
outcomes(1,:) = [AVlow(1), CVlow(1), AVhi(1),CVhi(1)];

for rr = 1:Multiplier
    SDlow(rr) = (std(ISI(rr*10,101:201)))^2;
    AVlow(rr) = (mean(ISI(rr*10,101:201)))^3;
    CVlow(rr) = sqrt((SDlow(rr)/AVlow(rr)))*100;
    SDhi(rr) = (std(ISI(rr*10,350:400)))^2;
    AVhi(rr) = (mean(ISI(rr*10,350:400)))^3;
    CVhi(rr) = sqrt((SDhi(rr)/AVhi(rr)))*100;
    DRlow(rr) = mean(MultFR(rr*10,101:201));
    DRhi(rr) = mean(MultFR(rr*10,350:400));
    outcomes(rr,:) = [DRlow(rr), CVlow(rr), DRhi(rr),CVhi(rr)];
end

Stats = {'DRi' 'CVi' 'DRf' 'CVf'};                  %Stats outputs a cell array with Discharge rate and CV initially (0.5-1 s) and at end (3.5- 4 s)
Stats(2:Multiplier+1,1) = num2cell(DRlow);          %(It's a cell array so each column can be labeled)
Stats(2:Multiplier+1,2) = num2cell(CVlow);
Stats(2:Multiplier+1,3) = num2cell(DRhi);
Stats(2:Multiplier+1,4) = num2cell(CVhi);

%% Isometric force production model from motor unit firing times

P = zeros(1,ActiveUnits(500));   %Preallocate vector of peak twitch force for each motor unit
b = log(100)/120;            %Constant multiplier for setting up force multiplier
for nn = 1:ActiveUnits(500)
    P(nn) = exp(b*nn);
end

TwitchTimes = zeros(1,ActiveUnits(500));  %Preallocate twitch contraction times
c = logb(100,3);
d = 1/c;
for ccc = 1:ActiveUnits(500)
    TwitchTimes(ccc) = 90*((1/P(ccc))^d);
end

NormTi = zeros(ActiveUnits(500),length(Xs));        
for hhh = 1:length(Xs)
    NormTi(:,hhh) = (TwitchTimes')./ISI2(:,hhh);
end
NormTi(NormTi==inf)=NaN;

force_mat = zeros(ActiveUnits(500),length(Xs));  %Preallocate the force for each motor unit at each time point
PoverT = P./TwitchTimes;
for fff = 1:length(Xs)
    force_mat(:,fff) = PoverT.* exp(1-(1/TwitchTimes')); %Gives the force output with no variability
end

S = -2.*((NormTi).^3);        %Makes the sigmoid to determine gain
S = 1 - exp(NormTi);
G = abs(S./NormTi);
G(G==inf)=0;

for rrrr = 1:ActiveUnits(500)
    for cccc = 1:length(Xs)
        if NormTi(rrrr,cccc) > 0.4
            force_mat(rrrr,cccc) = force_mat(rrrr,cccc) * G(rrrr,cccc);
        end
    end
end

for hello = 1:ActiveUnits(500)
    for bye = 1:length(Xs)
        if AboveRTE(hello,bye) == 0;
           force_mat(hello,bye) = 0;
        end
    end
end

Total_Force = sum(force_mat);
figure(2);
plot(Total_Force);          %Scaling up to arbitrary units
% ylim([0,10]);
ylabel('Force (au)');
xlabel('Time (ms)');
mean_force = mean(Total_Force(301:500));
SD_force = std(Total_Force(301:500));
CV_force = (SD_force./mean_force) * 100

%% Find a spectral density estiamte. 

Fout = Total_Force;
FS = 100;
figure(3);
pmtm(Fout,4,length(Fout),FS);
OutputDR = mean(MultFR(:,300:500),2);

%% Interval histograms

ISI2(ISI2 == 0) = NaN; 
ISI2(ISI2 > 200) = NaN;
figure(4)
histogram(ISI2(47,:))

Kurtosis = kurtosis(ISI2(47,300:500));
Skewness = skewness(ISI2(47,300:500));

Kurt_Value = nanmean(Kurtosis);
Skew_Value = nanmean(Skewness);
%dlmwrite('Model_Output2.csv',AvgFiringRate','-append');

figure(5);
L = logical(Try);
plotSpikeRaster(L,'PlotType','Scatter');

%% Raster 
rng(123);

lambda = zeros(neurons, length(force));
for j = 1:length(force)
    for mmmm = 1:neurons
        lambda(mmmm,j) = MultFR(mmmm,j);    
    end 
end
   
time = floor(exprnd(lambda)); 
Spike_vec = zeros(neurons,5000);
for lm = 1:neurons
    for jj = 10:length(force)
        if time(lm,jj) > 0
            Spike_vec(lm,(jj*10)+time(lm,jj)) = 1;
        else
            Spike_vec(lm,(jj*10)+time(lm,jj)) = 0;
        end
    end
end

newL = floor(length(Spike_vec)/10)
for nn = 1:neurons
    for d = 0:(newL-1)
        dsSub(nn,d+1) = sum(Spike_vec(nn,(d*10)+1:(d*10)+10));
    end
end

seq.y = dsSub
seq.T = length(dsSub);
seq.yr = MultFR;
%seq.u = ones(2,500).*0.5; Unsure the best inputs for u

figure;
plotSpikeRaster(logical(dsSub),'PlotType','Scatter');
set(gcf,'color','white')

Loading data, please wait...