Location: 12 L Platform 1 model codes @ 4c1ab73f48d7 / USMC / DM+Mijaailovich_huxley_copy / USM_DM.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/USM_DM.m

%%   Date: 01/03/2021
% 
%    Author: Anand Rampadarath
%    Version 1
%    This model represents the Uterine smooth muscle cell as 
%    developed by Maggio et al 2012. It is a variant of the 
%    Hai-Murphy model including spatian dependence of the 
%    cross-bridges, similar to Mijailovich et al. 
%    We us the DM approximation to reduce the computational complexity.
% 

function [X0,T,R,force]=USM_DM(trange,N,M,xrange)
%force,stiffness,phosphorylation,m10,m20,nmp,nm,C_cai
%N number of time points
%trange specifies the interval of time steps
%T time grid
%X0 space grid
%gamma=25 rho=1 s=airway generation lambda=force activation

   T=linspace(trange(1),trange(2),N);
   X0=linspace(xrange(1),xrange(2),M);
   dt=(trange(2)-trange(1))/(N); % dt must be <=0.1
  
   M10=0.309120293063990;
   M11=0.118325676501984;
   M12=0.066855519431672;
   M20=0.134379171029593;
   M21=0.045132660622671;
   M22=0.101673063419951;
   C=0.6;
    
   R=[M10;M11;M12;M20;M21;M22;C];
   kappa = 5.0859;
   load('k1.mat');
  
%     g=@(r) mis1(r);
%     options = optimoptions('fsolve');
%     [R,~,exitflag]=fsolve(g,R,options);
%     return

 for i=1:N %loop over number of timesteps to update then run euler
  
        t=trange(1)+i*dt;
         X=X0;
      %V=velvar(t,f); 
%       [~,V]=lengthvar(t,f);
      
      [F1,~,~,~,~]=DM_funcs_USM(t,R,k_1(i));
      [F2,~,~,~,~]=DM_funcs_USM(t+(dt/2),R+(dt/2)*F1,k_1(i));
      [F3,~,~,~,~]=DM_funcs_USM(t+(dt/2),R+(dt/2)*F2,k_1(i));
      [F4,q1,q2,p1,p2,C_cai(i)]=DM_funcs_USM(t+dt,R+(dt*F3),k_1(i));

      Rnew = R + (dt/6)*(F1+(2*F2)+(2*F3)+F4);
      R=Rnew;
      
      % A rebuild of the distribution approximation for 
     n1=(R(1)/((sqrt(2*pi))*q1)).*exp(-((X-p1).^(2))/(2*(q1.^(2))));
     n2=(R(4)/((sqrt(2*pi))*q2)).*exp(-((X-p2).^(2))/(2*(q2.^(2))));
     nmp(i)=1-R(7)-R(1);
     nm(i)=1-R(4)-R(1)-nmp(i);
     m10(i)=R(1);
     m11(i)=R(2);
     m21(i)=R(5);
     m20(i)=R(4);
     force(i)=(R(2)+R(5));
     stiffness(i)=R(1)+R(4);
     nMpave(i)=nmp(i);
     nAmpave(i)=R(1);
     phosphorylation(i)=(nMpave(i)+nAmpave(i));

 end
 
 plot(T,force,T,phosphorylation,T,stiffness)
end