MotorUnitModel.m

function [output] = MotorUnitModel(time,U,modelParameter,Fs)

N = modelParameter.N; %number of motor unit
i = 1:N; %motor unit identification index
RR = modelParameter.RR; %range of recruitment in unit of fold
a = log(RR)/N; %coefficient to establish a range of threshold values
RTE = exp(a*i); %recruitment threshold excitation
MFR = modelParameter.MFR; %minimum firing rate constant for all motoneurons
g_e = modelParameter.g_e; %missing parameter from the paper
PFR1 = modelParameter.PFR1; %the peak firing rate of the first recruited motoneuron in unit of impulse/sec
PFRD = modelParameter.PFRD; %the desired difference in peak firing rates between the first and last units in unit of impulse/sec
RTEn = exp(a*N); %recruitment threshold of the last motor unit
PFR = PFR1 - PFRD * (RTE./RTEn); %peak firing rate
PFRn = PFR1 - PFRD; %peak firing rate of the last motor unit
Emax = RTEn + (PFRn - MFR)/g_e; %maximum excitatory input
cv = modelParameter.cv; %ISI variability as per coefficient of variation (=mean/SD)

RP = modelParameter.RP; %range of twich force across motor untis in unit of fold
b = log(RP)/N; %coefficient to establish a range of twich force values
P = exp(b*i); %force generated by a motor unit as a function of its recruitment threshold
T_L = modelParameter.T_L; %the longest duration contraction time desired for the pool in unit of ms
RT = modelParameter.RT; % range of contraction time in unit of fold
c = log(100)/log(RT); %coefficient to establish a range of contraction time values
T = (T_L.* (1./P).^(1/c))./1000; %contraction time
t_twitch = 0:1/Fs:1;
twitch = zeros(N,length(t_twitch));

for j = 1:N
    twitch(j,:) =  P(j).*t_twitch./T(j).*exp(1-t_twitch./T(j));
end

FR_mat = zeros(N,length(time));
outputG = zeros(N,length(time));
spike_time = zeros(1,N);
spike_train = zeros(N,length(time));
force = zeros(N,length(time));
outputIndex = cell(1,N);

E = U*Emax;

for t = 1:length(time)
    if t > 1
        FR = g_e.*(E(t) - RTE) + MFR;
        FR(FR<MFR) = 0;
        FR_mat(:,t) = FR;
       
        index_1 = i(FR >= MFR & FR_mat(:,t-1)' == 0);        
        index_2 = i(FR >= MFR & spike_time==t);
        index = [index_1 index_2];
        outputIndex{t} = index;
        for j = 1:length(index)
            n = index(j);
            if FR(n) >= PFR(n)
                FR(n) = PFR(n);
                FR_mat(n,t) = FR(n);
            end
            spike_train_temp = zeros(1,length(time));
            if ~any(spike_train(n,:)) % initial time
                spike_train(n,t) = 1;
                spike_train_temp(t) = 1;
                mu = 1/FR(n);
                Z = randn(1);
                Z(Z>3.9) = 3.9;
                Z(Z<-3.9) = -3.9;
                spike_time_temp = (mu + mu*cv*Z)*Fs;
                spike_time(n) = round(spike_time_temp) + t;                                
                force_temp = conv(spike_train_temp,twitch(n,:));
                force(n,:) = force(n,:)+ force_temp(1:length(time));
            else
                if spike_time(n) == t
                    spike_train(n,t) = 1;
                    spike_train_temp(t) = 1;
                    mu = 1/FR(n);
                    Z = randn(1);
                    Z(Z>3.9) = 3.9;
                    Z(Z<-3.9) = -3.9;
                    spike_time_temp = (mu + mu*cv*Z)*Fs;
                    spike_time(n) = round(spike_time_temp) + t;
                    
                    ISI = (spike_time(n) - t)/Fs;
                    StimulusRate = T(n)/ISI;
                    if StimulusRate > 0 && StimulusRate <= 0.4
                        g = 1;
                    elseif StimulusRate > 0.4
                        S_MU = 1 - exp(-2*(StimulusRate)^3);
                        g = (S_MU/StimulusRate)/0.3;
                    end
                    outputG(n,t) = g;
                                      
                    force_temp = conv(spike_train_temp,g*twitch(n,:));
                    force(n,:) = force(n,:)+ force_temp(1:length(time));
                elseif FR_mat(n,t-1) == 0
                    spike_train(n,t) = 1;
                    spike_train_temp(t) = 1;
                    mu = 1/FR(n);
                    Z = randn(1);
                    Z(Z>3.9) = 3.9;
                    Z(Z<-3.9) = -3.9;
                    spike_time_temp = (mu + mu*cv*Z)*Fs;
                    spike_time(n) = round(spike_time_temp) + t;                                        
                    force_temp = conv(spike_train_temp,twitch(n,:));
                    force(n,:) = force(n,:)+ force_temp(1:length(time));
                    
                end
            end
            
        end
    end
    
end

output.Force = force;
output.SpikeTrain = spike_train;
output.TotalForce = sum(force);
output.FR = FR_mat;
output.g = outputG;
output.index = outputIndex;
output.Ur = RTEn/Emax;

end