p2p_Amplitude_test.m

close all
clear all
clc

Fs = 1000;
t = 0:1/Fs:5;
t_stim = 0:1/Fs:3;
N = 300; %number of motor unit
i_MU = 1:N; %motor unit identification index
RP = 100; %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_MU); %force generated by a motor unit as a function of its recruitment threshold
T_L = 90; %the longest duration contraction time desired for the pool in unit of ms
RT = 3; % 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
P = exp(b*i_MU); %force generated by a motor unit as a function of its recruitment threshold

RR = 30; %range of recruitment in unit of fold
a = log(RR)/N; %coefficient to establish a range of threshold values
RTE = exp(a*i_MU); %recruitment threshold excitation
MFR = 8; %minimum firing rate constant for all motoneurons
g_e = 1; %missing parameter from the paper
PFR1 = 35; %the peak firing rate of the first recruited motoneuron in unit of impulse/sec
PFRD = 10; %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


testedUnit = 1;
PFR = PFR(testedUnit);
%FR = [2 5 10 15 20 25 30 35 40 45 50];
FR = 1:1:70;
%FR = [1 5 10 20 40];
%FR = [2 15 25 35 70];
t_twitch = 0:1/Fs:1;
twitch = zeros(N,length(t_twitch));
force = zeros(N,length(t));
for j = 1:N
    twitch(j,:) =  P(j).*t_twitch./T(j).*exp(1-t_twitch./T(j));
end

for i = 1:length(FR)
    CT = T(testedUnit);
    ISI = 1/FR(i);
    SR = CT/ISI;
    
    spikeTrain_temp = spikeTrainGenerator(t_stim,Fs,FR(i));
    spikeTrain = [zeros(1,Fs) spikeTrain_temp zeros(1,Fs)];
    if SR <= 0.4
        g = 1;
    else
        S = 1 - exp(-2*SR^3);
        g = (S/SR)./0.3;
    end
    force_temp = conv(spikeTrain,g*twitch(testedUnit,:));
    force = force_temp(1:length(t));
    
    amp(i) = mean(force(2.5*Fs:3.5*Fs));
    p2p_amp(i) = (max(force(2.5*Fs:3.5*Fs))-min(force(2.5*Fs:3.5*Fs)));  
    
    figure(1)
    plot(t,force)
    hold on
end

%%
figure(1)
%xticks([0 1 2 3 4 5])
%yticks(0:2:10)
xlabel('Time (s)')
ylabel('Force (AU)')
set(gca,'TickDir','out');
set(gca,'box','off')

p2p_amp = 1 - p2p_amp./p2p_amp(1);
figure(2)
plot(FR,amp/amp(end)*100,'LineWidth',2)
xlabel('Frequency (Hz)','FontSize',14)
ylabel('Mean Force (%)','FontSize',14)
h1 = line([MFR MFR],[0 120]);
h2 = line([PFR PFR],[0 120]);
set(gca,'TickDir','out');
set(gca,'box','off')

%patch([8 35 35 8],[0 0 10 10],'r')
%set(gca,'children',flipud(get(gca,'children')))

figure(3)
plot(FR,p2p_amp*100,'LineWidth',2)
xlabel('Frequency (Hz)','FontSize',14)
ylabel('Degree of Fusion (%)','FontSize',14)
h3 = line([MFR MFR],[-20 120]);
h4 = line([PFR PFR],[-20 120]);
set(gca,'TickDir','out');
set(gca,'box','off')

figure(4)
plot(amp/amp(end)*100,p2p_amp*100,'LineWidth',2,'color','b')
xlabel('Mean Force','FontSize',14)
ylabel('Degree of Fusion (%)','FontSize',14)
set(gca,'TickDir','out');
set(gca,'box','off')
hold on
plot(0:1:100,0:1:100,'--','color','k')