Code availability

This code was written in MATLAB, but this file format is not supported by this upload system.

Below is the text of the code.

To receive a copy of the MATLAB code saved in MATLAB format, you can contact the authors directly.

function Rotation_multiple_clutch(T,rt,pd,ksub,rdsd,deg0)
% Input: simulating step dt, total time T, time step for storage storesize, 
%        amplitude ratio rt, period pd, substrate stiffness ksub, 
%        random seeds rdsd, initial angle rdn

% Output: None. Save the simulation data to file below according to given parameters.

% Seed random number generator using input
rng(rdsd);

% Motor Clutch parameters from original Chan Odde Model
nmotor = 75; % number of myosin motors
nclutch = 75; % number of clutches
Fm = -2; % pN, single myosin motor stall force
Fstall = nmotor*Fm; % pN, total motor stall force 
vu = -120; % nm/s, unloaded myosin motor velocity
kon = 1; % s^{-1}, on-rate constant of clutch binding to F-actin
koff = 0.1; % s^{-1}, off-rate constant for clutch unbinding from F-actin
Fb = -2; % pN, characteristic bond rupture force
kc = 0.8; % pN/nm, molecular clutch spring constant
%ksub = ; % pN/nm, compliant substrate spring constant (Input)

% New parameters introduced into model 
l0 = 5e3; % nm, initial length of filament
%deg0 = ; degrees, initial filament orientation (Input)
theta0 = deg2rad(deg0); %radians, initial filament orientation 
vpmax = 120; % nm/s, maximum filament polymerization rate 
lmin = 500; % nm, minimal filament length
lmax = 1e4; % nm, maximal filament length
%rt = ; stretch ratio (Input), if rt = 1.1, substrate is stretched 110% of its
%original size
%pd = ; s, stretch cosine function period (1/frequency) (Input)

% setup
dt = 1e-2; % s, lTime step duration 
%T = ; % s,  total time duration (Input)
Nstep = floor(T/dt)*10; % maximum number of time steps
storesize = dt; % s, time step for storage
Nstore = floor(T/storesize);

% initialization
vflm = vu; % Set filament velocity to unloaded velocity 
l = l0; % Set filament length to initial length
theta = theta0; % Set orientation to initial orientation 
theta_rltv = theta0; % Set relative filament orientation 
init_angle = theta_rltv;
lA = l0; % Initialize anchor point to end of filament 
lsub = 0; % nm, displacement of substrate

clstate = false(1,nclutch); % Set all clutches to unbound 
lc = l0*ones(1,nclutch); % nm, displacement of engaged molecular clutches

koffF = koff*ones(1,nclutch); % initial off rate
Fclutch = zeros(1,nclutch); % Force on each clutch
a = zeros(1,nclutch); % propensities

% Initialize data storage 
l_all = zeros(1,Nstore);
lsub_all = zeros(1,Nstore);
lc_all = zeros(nclutch,Nstore);
lA_all = zeros(1,Nstore);
vflm_all = zeros(1,Nstore);
t_all = zeros(1,Nstore); 
theta_all = zeros(1,Nstore);
theta_rltv_all = zeros(1,Nstore);
nonstate_all = zeros(1,Nstore);

l_avg_all = zeros(1,Nstore);
lsub_avg_all = zeros(1,Nstore);
t_avg_all = zeros(1,Nstore);
theta_avg_all= zeros(1,Nstore);
theta_rltv_avg_all = zeros(1,Nstore);
lA_avg_all = zeros(1,Nstore);
vflm_avg_all = zeros(1,Nstore);
nonstate_avg_all = zeros(1,Nstore);

       
indstore1 = 1;
indstore = 1;

t=0;

l_avg = 0;
lsub_avg = 0;
theta_avg = 0;
theta_rltv_avg = 0;
lA_avg = 0;
vflm_avg = 0;
nonstate_avg = 0;
counter = 0;

for ind = 1:Nstep
   
   if t>T break; end

   % asynchronized update
   
   % Get clutch states
   offst = find(clstate == false);
   onst = find(clstate == true);   
   

   % Gillespie algorithm to determine tau and update t 
   a(offst) = kon; 
   a(onst) = koffF(onst); 
   a0 = sum(a); 
   r1 = rand(); 
   tau = -log(r1)/a0;
   
   % Update polymerization rate according to the current length of
   % filament
   vpoly = vpmax*((lmax-l)/lmax);
   
   % Update filament length according to polymerization rate law
   l = l + (vpoly + vflm)*min(tau, dt);
   if l <= lmin
       disp('Filament too short!');
   end
   
   % Update time step with an upper bound according to gillespie
   % algorithm
   if tau < dt
       r2 = rand(); 
       b = cumsum(a);
       indr = find(r2*a0<b,1); 
       t = t + tau;
       indrstate = clstate(indr); % Find state of the clutch to be changed
       flag = true; % Change state
   else 
% If tau is larger than dt
       t = t+dt; 
       tau = dt;
       indr = 1e4; % false index, if used, will report error
       flag = false; % will not change state
       
   end
   
   % Update on-state clutch binding position lc
   lc(onst) = lc(onst) + vflm*tau;
   indneg = find(lc<=0); 
   clstate(indneg) = false; 
   non = sum(clstate);
   
   % Clutch binding/unbinding and updating of variables 
   
   % If no clutches are bound, and one clutch binds, initialize new anchor
   % point and set clutch location to tip of filament 
   if non == 0 && flag == true
       
       clstate(indr) = true;
       % Calculate the lagrangian coordinate (X,Y) of the anchor pointp
       X = l*cos(theta)/((1-cos(2*pi*t/pd))/2*(rt-1)+1);
       Y = l*sin(theta);
       
       % Calculate relative angle of filament
       theta_rltv = atan(Y/X);
       
       % Update clutch binding positions and off-rate
       lc(indr) = l; 
       koffF(indr) = koff;

       vflm = vu; 
       
   % Elseif any clutches are bound     
   elseif non>0 
       
       % Update anchor point position and lA 
       xA = X*((1-cos(2*pi*t/pd))/2*(rt-1)+1); % Postion of the clutch at current time point
       yA = Y;
       lA = sqrt(xA^2+yA^2);

       % Update filament angle
       theta = atan(yA/xA);
       
       % Update lsub according to force-balance
       lsub = (kc*sum(lc(clstate))+ksub*lA)/(kc*sum(clstate)+ksub);

       % Update filament velocity
       vflm = vu*(1-ksub*(lsub-lA)/Fstall);

       % Update on-state clutch koffF
       indc1 = find(clstate == true); 
       koffF(indc1) = koff*exp(abs(kc*(lsub-lc(indc1))/Fb));

       % Update new binding/unbinding clutch
       if flag == true
           if indrstate == false && l>=lsub
               clstate(indr) = true;
               lc(indr) = lsub; 
               koffF(indr) = koff;
           elseif indrstate == false && l<lsub
               % clutch stays off, too far to bind
           elseif indrstate == true
               clstate(indr) = false;
               if sum(clstate) == 0 % all clutch unbound
                   vflm = vu;
               end
           end
       end
       

   end
   
   % Save discrete data
   if t>indstore1*storesize
       l_all(indstore1) = l;
       lsub_all(indstore1) = lsub;
       t_all(indstore1) = t;
       theta_all(indstore1) = theta;
       theta_rltv_all(indstore1) = theta_rltv;
       
       lc_all(:,indstore1) = lc.';
       lA_all(indstore1) = lA;
       vflm_all(indstore1) = vflm;
       nonstate_all(indstore1) = sum(clstate);
       
       % increment indstore
       indstore1 = indstore1 +1;
   end
   
   % Averaged data
   l_avg = l_avg+l;
   lsub_avg = lsub_avg+lsub;
   theta_avg = theta_avg+theta;
   theta_rltv_avg = theta_rltv_avg+theta_rltv;
   lA_avg = lA_avg+lA;
   vflm_avg = vflm_avg+vflm;
   nonstate_avg = nonstate_avg+sum(clstate);
   counter = counter + 1;
   if t>indstore*storesize
       l_avg_all(indstore) = l_avg/counter;
       lsub_avg_all(indstore) = lsub_avg/counter;
       t_avg_all(indstore) = t;
       theta_avg_all(indstore) = theta_avg/counter;
       theta_rltv_avg_all(indstore) = theta_rltv_avg/counter;
       lA_avg_all(indstore) = lA_avg/counter;
       vflm_avg_all(indstore) = vflm_avg/counter;
       nonstate_avg_all(indstore) = nonstate_avg/counter;
       
       % increment indstore
       indstore = indstore +1;
       
       % temporary values
       l_avg = 0;
       lsub_avg = 0;
       theta_avg = 0;
       theta_rltv_avg = 0;
       lA_avg = 0;
       vflm_avg = 0;
       nonstate_avg = 0;
       counter = 0;
       
   end
   

   end_angle = theta_rltv;
   
end

% Save data to file
filename = ;
save(filename);