sizePrinciple_fig.m

close all
clear all
clc

%--------------------------------------------------------------------------
% Motor unit parameters
%--------------------------------------------------------------------------
Fs = 1000;
time = 0:1/Fs:5;

modelParameter.N = 120;
modelParameter.RR = 30;
modelParameter.MFR = 8;
modelParameter.g_e = 1;
modelParameter.PFR1 = 35;
modelParameter.PFRD = 10;
modelParameter.RP = 100;
modelParameter.T_L = 90;
modelParameter.RT = 3;

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


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

%--------------------------------------------------------------------------
% Simulation parameters
time_stim = 0:1/Fs:3;
mean_force = zeros(1,N);
%--------------------------------------------------------------------------
for n = 1:N
    FR = PFR(n);
    spikeTrain_temp = spikeTrainGenerator(time_stim,Fs,FR);
    spikeTrain = [zeros(1,1*Fs) spikeTrain_temp zeros(1,1*Fs)];
    
    StimulusRate = T(n)*FR;    
    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
    
    force_temp = conv(spikeTrain,g*twitch(n,:));
    force = force_temp(1:length(time));
    
    mean_force(n) = mean(force(2*Fs:3*Fs));  
    
end

figure(1)
plot(mean_force./mean_force(end)*100,RTE./Emax*100,'ok')
xlim([0 105])
xlabel('Force Production Capability (%)','FontSize',14)
ylabel('Recruitment Threshold (%)','FontSize',14)