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 12/16, this creates a set of N neurons with a latent
%state that is 2 dimensional. 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

clear all
%% Make force trace 
GoalF = 7;
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. Once it passes a percentage of max contraction (MVC),
%the unit is 'turned on'

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 in this case.
for ff = 1:120
    RTE(ff) = exp(a*ff);
end

%Need to calcualte, at each instant, how far above recruitment threshold a
%each motor 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). Seq is fed into decomposition algorithm.
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 the response to the 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. 0.5 is just a scaling factor
LinMult = tanh(LinMult) * 0.5 + 0.5;             %Set up an approximate sigmoid relation
LinMult = LinMult(1:ActiveUnits(500),:);         %Get to correct size for a given trial

FiringRate = zeros(ActiveUnits(500),500);           %preallocate a discharge rate (intensity) matrix
AboveRTE = AboveRTE(1:ActiveUnits(500),1:500);
NoiseMat = ones(ActiveUnits(500),1); 
NoiseMat = NoiseMat .* (AboveRTE(:,500).* 5);       %In this way, a unit closer to recruitment threshold receives more noise than one further away
NoiseMat(NoiseMat<0) = 0;                           %5 is another scaling factor


%% 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);

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;

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



plot(params.model.Xpca(1,:),'k', 'LineWidth',2)
hold on
%plot(params.model.Xpca(2,:), 'LineWidth',2, 'Color',[0.501960813999176 0.501960813999176 0.501960813999176])
set(gca,'TickDir','out', 'fontsize',16);%, 'XTickLabel',{'0','1','2','3','4','5'}); 
box off
set(gcf,'color','w')
hold off
set(gcf, 'WindowStyle','normal')
set(gcf, 'Units', 'Inches', 'Position', [0, 0, 6, 4], 'PaperUnits', 'Inches', 'PaperSize', [4, 6])

Loading data, please wait...