Location: 12 L Platform 1 model codes @ 4c1ab73f48d7 / USMC / DM+Mijaailovich_huxley_copy / HHM1.m

Author:
aram148 <a.rampadarath@auckland.ac.nz>
Date:
2021-11-04 16:02:42+13:00
Desc:
Updated USMC-Bursztyn model
Permanent Source URI:
http://models.cellml.org/workspace/6b0/rawfile/4c1ab73f48d7150c2094cf8017425482a1594ccf/USMC/DM+Mijaailovich_huxley_copy/HHM1.m
Compare... Download Source 

function [X0,T,R,Force,Stiffness,X,N]=HHM1(trange,N,M,xrange,kappa)
%X0,T,R,Force,Stiffness,Phosphorylation,UnattachedMyosin,UnattachedPhosMyosin,AttachedPhosMyosin,AttachedMyosin
%N number of time points
%trange specifies the interval of time steps
%T time grid
%X0 space grid
T=linspace(trange(1),trange(2),N);
%X=linspace(xrange(2),xrange(1),M);
X0=linspace(xrange(1),xrange(2),M);
%dt=timestep
dt=(trange(2)-trange(1))/(N);


% for i=1
%     t=trange(1)+1*dt;
%     [alpha, ~]=lengthvar(t,f);
%     X=X0-alpha;
%     [k1 k6]=bindingconst(t);%update k1 and k6 for each T
%     
%     [gx  gpx fpx] = bindingfuncs(X,xrange);
    %             

    
%     n11=(fpx .* gx * k5 + gx .* gpx * k2 + gpx * k2 * k6 + gx .* k2 * k5) ./ (fpx .* gx * k1 + fpx .* gx * k5 + k5 * fpx * k1 +...
%         fpx * k1 * k6 + gpx .* gx * k1 + gx .* gpx * k2 + gpx * k1 * k6 + gpx * k2 * k6 + gx * k1 * k5 + gx * k2 * k5);
%     
%     n22=(gx .* gpx + gpx * k6 + gx * k5) * k1 ./ (fpx .* gx * k1 + fpx .* gx * k5 + k5 * fpx * k1 + fpx * k1 * k6 + gpx .* gx * k1 + gx .* gpx * k2 + gpx * k1 * k6 + gpx * k2 * k6 + gx * k1 * k5 + gx * k2 * k5);
%     
%     
%     n33=((k6 + gx) .* fpx * k1) ./ (fpx .* gx * k1 + fpx .* gx * k5 + k5 * fpx * k1 + fpx * k1 * k6 + gpx .* gx * k1 + gx .* gpx * k2 + gpx * k1 * k6 + gpx * k2 * k6 + gx * k1 * k5 + gx * k2 * k5);
% 
%     n44=(k5 * fpx * k1) ./ (fpx .* gx * k1 + fpx .* gx * k5 + k5 * fpx * k1 + fpx * k1 * k6 + gpx .* gx * k1 + gx .* gpx * k2 + gpx * k1 * k6 + gpx * k2 * k6 + gx * k1 * k5 + gx * k2 * k5);

   
    %  Constraint eqn substituted for N_am eqn
%     n11= gx .* (fpx .* gpx + fpx * k5 + gpx * k2 + k2 * k5) ./ (fpx .* gpx .* gx + fpx .* gpx * k1 + fpx .* gx * k1 + fpx .* gx * k5 + fpx * k1 * k5 + fpx * k1 * k6 + gpx .* gx * k1 +...
%         gpx .* gx * k2 + gx * k1 * k5 + gx * k2 * k5);
%      n22=gx .* k1 .* (k5 + gpx) ./ (fpx .* gpx .* gx + fpx .* gpx * k1 + fpx .* gx * k1 + fpx .* gx * k5 + fpx * k1 * k5 + fpx * k1 * k6 + gpx .* gx * k1 + gpx .* gx * k2 + gx * k1 * k5 + gx * k2 * k5);
%      n33=(k6 + gx) .* fpx * k1 ./ (fpx .* gpx .* gx + fpx .* gpx * k1 + fpx .* gx * k1 + fpx .* gx * k5 + fpx * k1 * k5 + fpx * k1 * k6 + gpx .* gx * k1 + gpx .* gx * k2 + gx * k1 * k5 + gx * k2 * k5);
%      n44=fpx .* k1 .* (k5 + gpx) ./ (fpx .* gpx .* gx + fpx .* gpx * k1 + fpx .* gx * k1 + fpx .* gx * k5 + fpx * k1 * k5 + fpx * k1 * k6 + gpx .* gx * k1 + gpx .* gx * k2 + gx * k1 * k5 + gx * k2 * k5);
    %
    %
    %  Constraint substituted for the N_amp equation
%          n11= gx .* k5 .* (k2 + fpx) ./ (fpx .* gx * k1 + fpx .* gx * k5 + fpx * k1 * k5 + fpx * k1 * k6 + gx * k1 * k5 + gx * k2 * k5);
%      n22=gx * k5 * k1 ./ (fpx .* gx * k1 + fpx .* gx * k5 + fpx * k1 * k5 + fpx * k1 * k6 + gx * k1 * k5 + gx * k2 * k5);
%      n33=(k6 + gx) .* fpx * k1 ./ (fpx .* gx * k1 + fpx .* gx * k5 + fpx * k1 * k5 + fpx * k1 * k6 + gx * k1 * k5 + gx * k2 * k5);
%      n44=fpx * k1 * k5 ./ (fpx .* gx * k1 + fpx .* gx * k5 + fpx * k1 * k5 + fpx * k1 * k6 + gx * k1 * k5 + gx * k2 * k5);
    
    
%     %  Constraint substituted fpr the N_mp equation
%     n11= (fpx .* gx * k5 + gpx .* gx * k2 + gpx * k2 * k6 + gx * k2 * k5 + k2 * k5 * k6) ./ (fpx .* gx * k1 + fpx .* gx * k5 + fpx * k1 * k5 + fpx * k1 * k6 + gpx .* gx * k1 + gpx .* gx * k2 + gpx * k1 * k6 + gpx * k2 * k6 +...
%         gx * k1 * k5 + gx * k2 * k5 + k1 * k5 * k6 + k2 * k5 * k6);
%     n22=(gpx .* gx + gpx * k6 + k5 * gx + k5 * k6) * k1 ./ (fpx .* gx * k1 + fpx .* gx * k5 + fpx * k1 * k5 + fpx * k1 * k6 + gpx .* gx * k1 + gpx .* gx * k2 + gpx * k1 * k6 + gpx * k2 * k6 +...
%         gx * k1 * k5 + gx * k2 * k5 + k1 * k5 * k6 + k2 * k5 * k6);
%     n33=(k6 + gx) .* fpx * k1 ./ (fpx .* gx * k1 + fpx .* gx * k5 + fpx * k1 * k5 + fpx * k1 * k6 + gpx .* gx * k1 + gpx .* gx * k2 + gpx * k1 * k6 + gpx * k2 * k6 + gx * k1 * k5 + gx * k2 * k5 + k1 * k5 * k6 + k2 * k5 * k6);
%     n44=fpx * k1 * k5 ./ (fpx .* gx * k1 + fpx .* gx * k5 + fpx * k1 * k5 + fpx * k1 * k6 + gpx .* gx * k1 + gpx .* gx * k2 + gpx * k1 * k6 + gpx * k2 * k6 + gx * k1 * k5 + gx * k2 * k5 + k1 * k5 * k6 + k2 * k5 * k6);
    
    
    % Constraint Substituded for the N_m equation
%     n11= (fpx .* gpx .* gx + fpx .* gpx * k6 + fpx .* gx * k5 + fpx * k5 * k6 + gpx .* gx * k2 + gpx * k2 * k6 + gx * k2 * k5 + k2 * k5 * k6) ./ (fpx .* gpx .* gx + fpx .* gpx * k6 + fpx .* gx * k1 +...
%         fpx .* gx * k5 + fpx * k1 * k5 + fpx * k1 * k6 + fpx * k5 * k6 + gpx .* gx * k1 + gpx .* gx * k2 + gpx * k1 * k6 + gpx * k2 * k6 + gx * k1 * k5 + gx * k2 * k5 + k1 * k5 * k6 + k2 * k5 * k6);
%     n22=(gpx .* gx + gpx * k6 + k5 * gx + k5 * k6) * k1 ./ (fpx .* gpx .* gx + fpx .* gpx * k6 + fpx .* gx * k1 + fpx .* gx * k5 + fpx * k1 * k5 + fpx * k1 * k6 + fpx * k5 * k6 +...
%         gpx .* gx * k1 + gpx .* gx * k2 + gpx * k1 * k6 + gpx * k2 * k6 + gx * k1 * k5 + gx * k2 * k5 + k1 * k5 * k6 + k2 * k5 * k6);
%     n33=(k6 + gx) .* fpx * k1 ./ (fpx .* gpx .* gx + fpx .* gpx * k6 + fpx .* gx * k1 + fpx .* gx * k5 + fpx * k1 * k5 + fpx * k1 * k6 + fpx * k5 * k6 + gpx .* gx * k1 +...
%         gpx .* gx * k2 + gpx * k1 * k6 + gpx * k2 * k6 + gx * k1 * k5 + gx * k2 * k5 + k1 * k5 * k6 + k2 * k5 * k6);
%     n44=fpx * k1 * k5 ./ (fpx .* gpx .* gx + fpx .* gpx * k6 + fpx .* gx * k1 + fpx .* gx * k5 + fpx * k1 * k5 + fpx * k1 * k6 + fpx * k5 * k6 + gpx .* gx * k1 + gpx .* gx * k2 +...
%         gpx * k1 * k6 + gpx * k2 * k6 + gx * k1 * k5 + gx * k2 * k5 + k1 * k5 * k6 + k2 * k5 * k6);
    

    
%     R = [n11 n22 n33 n44]; % R is a (4M,1) vector

n11=ones(M,1);

n22=zeros(M,1);
n33=zeros(M,1);

n44=zeros(M,1);

R = [n11; n22; n33; n44];
 
% end
% return


for i=1:N %loop over time to update then run euler
    
    t=0+i*dt;
    [F1,~]=funcs1(t,R,X0,M,T);
    [F2,~]=funcs1(t+(dt/2),R+((dt/2)*F1),X0,M,T);
    [F3,~]=funcs1(t+(dt/2),R+((dt/2)*F2),X0,M,T);
    [F4,X]=funcs1(t+dt,R+(dt*F3),X0,M,T);
    
    %main euler loop
    Rnew = R + (dt/6)*(F1+(2*F2)+(2*F3)+F4);
    R=Rnew;
    
    nM=R(1:M);
    nMp=R(M+1:2*M);
    nAMp=R((2*M)+1:3*M);
    nAM=R((3*M)+1:4*M);
    
    %       if (mod(i,500)==0);%&& i>=2000 && i<=2010)
    % % %
    % %             figure(1)
    % %             plot(X,nAMp)
    % %             title('nAMp for Mijailovich distribution')
    % %             drawnow
    % %             hold on
    % %            % pause
    %
    % %             figure(2)
    % %             plot(X,nAM)
    % %             title('nAM for Mijailovich  distribution')
    % %             drawnow
    % %             hold on
    %             %pause
    % % % % %
    %             figure(3)
    %             plot(X,nAM+nAMp)
    %             title('nAm+nAMp for Mijailovich')
    %             drawnow
    %             hold on
    %              %pause
    UnattachedMyosin(i)=trapz(X0,nM.*(X0'>0) .* (X0'<1));
    UnattachedPhosMyosin(i)=trapz(X0,nMp.*(X0'>=0) .* (X0'<=1));
    AttachedPhosMyosin(i)=trapz(X0,nAMp.*(X0'>=0) .* (X0'<=1));
    AttachedMyosin(i)=trapz(X0,nAM);
    Force(i)=kappa*trapz(X,X.'.*(nAMp+nAM));
    
    Stiffness(i)=trapz(X,(nAMp+nAM));
    Phosphorylation(i)= UnattachedPhosMyosin(i)+AttachedPhosMyosin(i);
end