- Author:
- GitHub <noreply@github.com>
- Date:
- 2025-03-01 10:02:00+13:00
- Desc:
- Update Mmodel_2025_Human.m
Fixed bug to allow user-inputted metabolites to be registered
- Permanent Source URI:
- https://staging.physiomeproject.org/workspace/d52/rawfile/c86fb8da1200cb2e9d11cd4513a6d3e6dbb19ead/Mmodel_2025_Human.m
%
% Function describing an ODE muscle model for human atria
% published by Musgrave et al in 2025.
%
% Ca2+ dependencies are based closely on those of Land et al. (2017)
%
% XB strain equations based on Razumova et al. (1999) 3 state stiffeness
% distortion model. Sarcomere length dependence on available cross-bridge
% proportion. Strain dependence on cross-bridge detachment (k3 rate) and
% attachment (k1 rate).
%
% Set up to be solved numerically under a range of sarcomere length inputs,
% including perturbation protocols.
%
% Inputs:
% - t: time (seconds)
% - y: states in the model
% - s: sarcomere length (um) input. Either a constant or 2 x n matrix
% with matching time and length points (t in 1st row and L in 2nd)
% - c: intracellular [Ca2+] (uM). Either a constant or 2 x n matrix
% with matching time and Ca points (t in 1st row and Ca in 2nd)
% - params: variable parameters fed into the model. Follows a
% structure including:
% .xb : cross-bridge params (vary for diabetic)
% = [k1 k-1 k2 k-2 k3 phi_x phi_v phi_l K [] phi_s1 phi_s3 kd_ATP []]
% .pass : passive model params (vary for diabetic)
% = [kP1 kP2 eta phi_e]
% .ca : Ca-related parameters (based on Land model, except Ca50)
% = [Ca50 n_Tm ktrpn n_trpn k0 k_0]
% .met : concentration of ATP and Pi (in mM)
% =[ATP Pi]
% .mode : mode for the simulation
% = either 'sarcomere' - for normal simulations with
% sarcomere control or 'loop'- for work loops
% .loop : parameters for the work-loop simulation
% = [Fafter mass visc]
% Outputs:
% - dydt: system of odes to solve
% - F: total (active + passive) stress solution of the model (kPa)
% - Fp: passive stress of the muscle (kPa)
% - k3: ATPase rate during simulation (s^-1)
% - k_3: reverse ATPase rate during simulation (s^-1)
%
% Author: Julia Musgrave
% Date : Sep 2024
function [dydt, F, Fp,k3,k_3] = Mmodel_2025_Human(t,y,s,c,params,met)
% ------------------ initialising state variables -----------------------
if nargin==0
dydt=[1e-4 1e-4 0 0.01 0 0.001 0.79 2.2 0];
return
end
%---------------------------- parameters --------------------------------
try
mode= params.mode; %'sarcomere'; 'loop'
catch
mode = 'sarcomere'; % SL is controlled by default
end
xb_params=params.xb;
p_params=params.passive;
ca_params=params.ca;
% cross-bridge parameters
k1=xb_params(1); % s^-1
k_1=xb_params(2); % s^-1
k2=xb_params(3); % s^-1
k_2=xb_params(4); % s^-1
k3=xb_params(5); % s^-1
phi_x=xb_params(6); % unitless (= 1 for human model)
phi_v=xb_params(7); % unitless
phi_l=xb_params(8); % unitless
K=xb_params(9); % GPa/m (kPa/um)
%Ks=xb_params(10); % zero in human model
phi_s1=xb_params(11); % um^-1
phi_s3=xb_params(12); % um^-1
kd_ATP=xb_params(13); % mM
%kd_Pi=xb_params(14); % not present in human model
% passive force parameters
kP1=p_params(1);
kP2=p_params(2);
eta=p_params(3);
phi_e=p_params(4);
% Ca regulation parameters
Ca50=ca_params(1);
n_Tm=ca_params(2);
ktrpn=ca_params(3);
n_trpn=ca_params(4);
k0=ca_params(5);
k_0=ca_params(6);
% model constants
xC0=0.01; % powerstroke size (um)
Lmax=2.3; % (um)
Lr=2.2*0.85;
ADP=0.036; % mM (based on Tran et al. 2010)
% input metabolite concentrations
if nargin>4
ATP=met(1); % (mM)
Pi=met(2); % (mM)
else
ATP=params.met(1);
Pi=params.met(2);
end
% ------------ getting sarcomere length, velocity and Ca -----------------
if size(y,2)==1; y=y'; end
[L, L0, dL, Ca] = get_sim_inputs(t,s,c,y(:,8));
if mode(1)=='s' && length(s)==1 % isosarcometric contraction
L=L0;
end
%------------------ state variables -------------------------------
% original 5
B=y(:,1);
C=y(:,2);
xB=y(:,3);
xC=y(:,4);
F1=y(:,5);
CaT=y(:,6);
N=y(:,7);
intF=y(:,9);
%XB_frac=(B+C)/(0.0137+0.0325);
% ------------------ cross bridge model ----------------------------------
% algebraic equations (depending on state variables)
k1=k1*exp(-phi_s1*xB);
k3=k3*exp(phi_s3*(xC-xC0));
Z=1 + (L<Lmax).*(L/Lmax-1)*phi_l;
A=max(0,Z-B-C-N);
% incorporating metabolite dependence into rates
k_1=k_1*Pi;
k3=k3*ATP/(ATP+kd_ATP);
k_2=k_2*kd_ATP/(ATP+kd_ATP);
% thermodynamic constraint for k-3 (as in Tran et al. (2010))
G_ATP0 = -30; % kJ/mol
R = 8.314/1000; % kJ/(mol.K)
T=37+273; % K (experiments performed at 37 degC)
G_ATP=G_ATP0+R*T*log((ADP/1000*Pi/1000)/(ATP/1000)); % convert mets to mol/L
k_3 = k1 .* k2 .* k3 ./ (k_1 * k_2 * exp(-G_ATP/(R*T)));
% ------------------ force calculations ----------------------------------
% Active stress (force) equation
Fa=K*(B.*xB+C.*xC); % (kPa)
% Passive stress
F2= kP2*(exp(phi_e*(L-Lr))-1);
Fp=F1+F2;
% Total stress (force)
F=Fa+Fp;
Fafter=0;
% finding dL when length perturbation not imposed
if dL(1) == 'i' % 'iso'
dL=0; % for isosarcometric case (including start of workloop)
if mode(1)=='l'
% funky stuff for loops
Fafter=params.loop(1); % kPa
mass=params.loop(2); % kPa.s^2/(um) = GPa.s^2/m
visc=params.loop(3); % kPa.s/(um) = GPa.s/m
global isotonic ;
% Check isotonic switch
if isotonic
dL = (intF + visc*(Lr - L))/mass;
if dL > 0 % need to switch back to isometric
dL=0;
end
end
if t>0.8 && L<L0% restretch
% start with a gradual acceleration
dL = sin(pi*(t-0.85)/.1)+1;
if t>0.85 % pick a constant velocity to reach starting SL by end of loop
dL = (L0 - L)/(0.95-t);
if t>0.95
dL=0;
end
end
end
% Prevent over extension in restretch
if t>0.8 && isotonic && (L > (L0*1.0001))
dL = 0;
end
if ~isotonic && F > Fafter
isotonic = true;
end
end
end
% ----------------------- ODEs -----------------------------------------
dydt=zeros(size(y));
%troponin
dydt(:,6) = ktrpn*((Ca./Ca50).^n_trpn.*(1-CaT)-CaT); % dCaT/dt
% XB states
dydt(:,1) = k1.*A + k_2.*C -(k_1+k2).*B; % dB/dt
dydt(:,2) = k2*B -(k_2+k3).*C + k_3.*A; % dC/dt
dydt(:,7) = k_0*CaT.^(-n_Tm/2).*A-k0*CaT.^(n_Tm/2).*N; %dN/dt
% XB strain
dydt(:,3) = -phi_x./B.*xB.*(A.*k1 + C.*k_2) + dL*phi_v; % dxB/dt
dydt(:,4) = -phi_x./C.*(xC-xC0).*(B.*k2+A.*k_3) + dL*phi_v; % dxC/dt
% passive stress
dydt(:,5) = kP1*(dL-F1/eta); % dF1/dt
% whole muscle mechanics
dydt(:,8) = dL; % dL/dt
dydt(:,9) = -Fa-Fp+Fafter;
if mode(1)=='l' && ~isotonic
dydt(:,9) = 0;
end
dydt=dydt(:);
end
