Computational Modeling of Coupled Energetics and Mechanics in the Rat Ventricular Myocardium

Bahador MarzbanRachel LopezDaniel A. Beard

This paper details a multi-scale model computational model of myocardial energetics—oxidative ATP synthesis, ATP hydrolysis, and phosphate metabolite kinetics—and myocardial mechanics used to analyze data from a rat model of cardiac decompensation and failure. Combined, these two models simulate cardiac mechano-energetics: the coupling between metabolic production of ATP and hydrolysis of ATP to generate mechanical work. The model is used to predict how differences in energetic metabolic state found in failing versus control hearts causally contribute to systolic mechanical dysfunction in heart failure. This Physiome paper describes how to access, run, and manipulate these computer models, how to parameterize the models to match data, and how to compare model predictions to data.

Keywords:cardiac energeticscardiac mechanicsheart failuremechano-energetic couplingmulti-scale modeling


The multi-scale modeling approach is summarized in Figure 1. Model components representing myocardial metabolism Bazil et al., 2016, myocardial cell mechanics Tewari et al., 2016Tewari et al., 2016, myocardial whole-organ pumping Lumens et al., 2009, and a simple lumped circulatory model are integrated together to simulate whole-body cardiovascular function.

All computer codes are available at https://​​github​​.com​​/beards​​-lab​​/Rat​​-Cardiac​​-Energetic.

  • The cardiac metabolic energetic model component is parameterized to match data from individual animals based on the oxidative capacity and cytoplasmic metabolite pools obtained from Lopez et al. (2020).

  • Certain parameters from the cross-bridge and calcium-activation models of Tewari et al. (2016)Tewari et al. (2016) and Campbell et al. (2018) are re-estimated to match data from Janssen et al. (2002) on calcium transients and force-generation in isolated rat cardiac trabeculae, as detailed below.

  • Wall volumes and anatomic parameters associated with the Lumens et al. (2009) heart model are identified based on anatomical data obtained from echocardiography and ex-vivo gross morphological measurements on individual animals from Lopez et al. (2020).

  • The simple lumped parameter circulatory model is identified based on cardiovascular state variables measured under resting conditions.

The resulting integrated model is used to predict the in vivo mechanical function and energetic state of the myocardium under resting conditions in each animal.

Multi-scale modeling of myocardial mechano-energetic function. The model integrates previously developed and validated models of cardiomyocyte dynamics , myocardial energetics , whole-organ cardiac mechanics  and a simple lumped parameter closed-loop circulatory system representing the systemic and pulmonary circuits. Data from multiple experimental modalities are used to identify model components for each individual animal in this study. The model predicts variables representing the in vivo myocardial energetic state, including ATP hydrolysis rate, [ATP], [ADP], [Pi], and the free energy of ATP hydrolysis DGATP in the LV myocardium for each individual animal. Figure reproduced from .

Figure 1:Multi-scale modeling of myocardial mechano-energetic function. The model integrates previously developed and validated models of cardiomyocyte dynamics Tewari et al., 2016Tewari et al., 2016, myocardial energetics Bazil et al., 2016Gao et al., 2019, whole-organ cardiac mechanics Lumens et al., 2009 and a simple lumped parameter closed-loop circulatory system representing the systemic and pulmonary circuits. Data from multiple experimental modalities are used to identify model components for each individual animal in this study. The model predicts variables representing the in vivo myocardial energetic state, including ATP hydrolysis rate, [ATP], [ADP], [Pi], and the free energy of ATP hydrolysis DGATP in the LV myocardium for each individual animal. Figure reproduced from Lopez et al. (2020).

2Model of Cardiac Energy Metabolism

2.1Model Variables:

The cellular energy metabolism model is based on the mitochondrial oxidative phosphorylation model of Bazil et al. (2016). The model is governed by 29 ordinary differential equations governing mitochondrial membrane potential, metabolite species concentrations, and cation (H+^+, K+^+, and Mg2+^{2+}) concentrations in the mitochondrial matrix, inter-membrane space, and cytosol. Table 1 lists the state variables of the model, with a brief description, units used in the model, and the variable name used in the model codes. The original formulation of the model accounted for reactive oxygen species O2_2^{\cdot -} and H2_2O2_2, which are ignored here, and thus the model is modified accordingly from Bazil et al. (2016).

The governing equations for these variables are delineated in Table 1.

Table 1:Energetics Model State Variables

State VariableDefinitionUnitsVariable name in code
ΔΨ\Delta\PsiMitochondrial membrane potentialmVDPsi_im_to_matrix
Mitochondrial Matrix State Variables
[ATP]x_xTotal matrix ATP concentrationMATP_matrix
[ADP]x_xTotal matrix ADP concentrationMADP_matrix
[Pi]x_xTotal matrix Pi concentrationMPi_matrix
[NADH]x_xTotal matrix NADH concentrationMNADH_matrix
[NAD]x_xTotal matrix NAD concentrationMNAD_matrix
[UQH2]x_xTotal matrix ubiquinol concentrationMcoQH2_matrix
[UQ]x_xTotal matrix ubiquinone concentrationMcoQ_matrix
[H+^+]x_xMatrix free proton concentrationMh_matrix
[K+^+]x_xMatrix free potassium concentrationMk_matrix
[Mg2+^{2+}]x_xMatrix free magnesium concentrationMm_matrix
Intermembrane Space (IMS) State Variables
[c2+^{2+}]i_iTotal cytochrome c2+^{2+} (reduced) concentrationMcytocred_im
[c3+^{3+}]i_iTotal cytochrome c3+^{3+} (oxidized) concentrationMcytocox_im
[ATP]i_iTotal IMS ATP concentrationMATP_matrix
[ADP]i_iTotal IMS ADP concentrationMADP_matrix
[AMP]i_iTotal IMS AMP concentrationMAMP_matrix
[Pi]i_iTotal matrix Pi concentrationMPi_im
[H+^+]i_iIMS free proton concentrationMh_im
[K+^+]i_iIMS free potassium concentrationMk_im
[Mg2+^{2+}]i_iIMS free magnesium concentrationMm_im
Cytosolic Space State Variables
[ATP]c_cTotal cytosolic ATP concentrationMATP_c
[ADP]c_cTotal cytosolic ADP concentrationMADP_c
[AMP]c_cTotal cytosolic AMP concentrationMAMP_c
[Pi]c_cTotal cytosolic Pi concentrationMPi_c
[CrP]c_cTotal cytosolic phosphocreatine concentrationMAMP_c
[Cr]c_cTotal cytosolic creatine concentrationMPi_c
[H+^+]c_ccytosolic free proton concentrationMh_c
[K+^+]c_ccytosolic free potassium concentrationMk_c
[Mg2+^{2+}]c_ccytosolic free magnesium concentrationMm_c

2.2Mitochondrial Membrane Potential:

The potential difference across the mitochondrial inner membrane is governed by currents across the membrane:

dΔΨ/dt=(4JC1+2JC3+4JC4nHF1F0JF1F0JANTJHleak)/Cmitod\Delta\Psi / dt = \left( 4J_{C1} + 2J_{C3} + 4J_{C4} - nH_{F1F0} J_{F1F0} - J_{ANT} - J_{Hleak} \right) / C_{mito}

where JC1J_{C1}, JC3J_{C3}, and JC4J_{C4} are the complex I, III, and IV fluxes, which are associated with pumping 4, 2, and 4 positive charges out of the matrix. The F1_1F0_0 ATPase turnover rate is JF1F0J_{F1F0} and nHF1F0nH_{F1F0} (= 8/3) is the proton flux stoichiometric number associated with the synthesis of one ATP. The fluxes JANTJ_{ANT} and JHleakJ_{Hleak} are the adenine nucleotide translocator and proton leak fluxes.

2.3Mitochondrial Matrix Metabolite State Variables:

Metabolite concentrations in the matrix are governed by:

d[ATP]x/dt=(JF1F0JANT)/Volxd[ADP]x/dt=(JANTJF1F0)/Volxd[Pi]x/dt=(JPICJF1F0)/Volxd[NAD]x/dt=(JC1JDH)/Volxd[NADH]x/dt=(JDHJC1)/Volxd[UQ]x/dt=(JC3αC2JDHJC1)/Volxd[UQH2]x/dt=(JC3+αC2JDH+JC1)/Volx\begin{align*} d[ATP]_x/dt &= \left(J_{F1F0} - J_{ANT} \right) / Vol_x \nonumber \\ d[ADP]_x/dt &= \left(J_{ANT} - J_{F1F0} \right) / Vol_x \nonumber\\ d[Pi]_x/dt &= \left(J_{PIC} - J_{F1F0} \right) / Vol_x \nonumber\\ d[NAD]_x/dt &= \left(J_{C1} - J_{DH} \right) / Vol_x \nonumber\\ d[NADH]_x/dt &= \left(J_{DH} - J_{C1} \right) / Vol_x \nonumber\\ d[UQ]_x/dt &= \left(J_{C3} - \alpha_{C2} J_{DH} - J_{C1} \right) / Vol_x \nonumber\\ d[UQH_2]_x/dt &= \left(-J_{C3} + \alpha_{C2} J_{DH} + J_{C1} \right) / Vol_x \end{align*}

where VolxVol_x is the water volume of the mitochondrial matrix in units of volume of matrix water space per unit mitochondrial volume, the fluxes in the right-hand sides of these expressions are in units of moles per unit liter of mitochondrial volume per unit time, and are defined below. The coefficient αC2\alpha_{C2} (=1/4=1/4) accounts for ubiquinone reduction via complex II.

2.4Inter-Membrane Space (IMS) Metabolite State Variables:

Metabolite concentrations in the intermembrane space are governed by:

d[c2+]i/dt=(2JC4+2JC3)/Volid[c3+]i/dt=(2JC42JC3)/Volid[ATP]i/dt=(JANT+JATPPERM)/Volid[ADP]i/dt=(JANT+JADPPERM)/Volid[AMP]i/dt=(JAMPPERM)/Volid[Pi]i/dt=(JPIC+JPIPERM)/Voli\begin{align*} d[c^{2+} ]_i/dt &= \left(-2J_{C4} + 2J_{C3} \right) / Vol_i \nonumber\\ d[c^{3+} ]_i/dt &= \left(2J_{C4} - 2J_{C3} \right) / Vol_i \nonumber\\ d[ATP]_i/dt &= \left(J_{ANT} + J_{ATPPERM} \right) / Vol_i \nonumber\\ d[ADP]_i/dt &= \left(-J_{ANT} + J_{ADPPERM} \right) / Vol_i \nonumber\\ d[AMP]_i/dt &= \left(J_{AMPPERM} \right) / Vol_i \nonumber\\ d[Pi]_i/dt &= \left(-J_{PIC} + J_{PIPERM} \right) / Vol_i \end{align*}

where VoliVol_i is the water volume of the mitochondrial inter-membrane space in units of volume of IMS water space per unit mitochondrial volume, the fluxes in the right-hand sides of these expressions are in units of moles per unit liter of mitochondrial volume per unit time, and are defined below.

2.5Cytosolic Metabolite State Variables:

Metabolite concentrations in the cytosolic space are governed by:

d[ATP]c/dt=(JATPase+JCK+JAKJATPPERMVRm/VRc)/Volcd[ADP]c/dt=(+JATPaseJCK2JAKJADPPERMVRm/VRc)/Volcd[AMP]c/dt=(+JAKJAMPPERMVRm/VRc)/Volcd[Pi]c/dt=(+JATPaseJPIPERMVRm/VRc)/Volcd[CrP]c/dt=(JCK)/Volcd[Cr]c/dt=(+JAK)/Volc\begin{align*} d[ATP]_c/dt &= \left(-J_{ATPase} + J_{CK} + J_{AK} - J_{ATPPERM} V_Rm/V_Rc \right) / Vol_c \nonumber\\ d[ADP]_c/dt &= \left(+J_{ATPase} - J_{CK} - 2J_{AK} - J_{ADPPERM} V_Rm/V_Rc \right) / Vol_c \nonumber\\ d[AMP]_c/dt &= \left(+J_{AK} - J_{AMPPERM} V_Rm/V_Rc \right) / Vol_c \nonumber\\ d[Pi]_c/dt &= \left(+J_{ATPase} - J_{PIPERM} V_Rm/V_Rc \right) / Vol_c \nonumber\\ d[CrP]_c/dt &= \left(-J_{CK} \right) / Vol_c \nonumber\\ d[Cr]_c/dt &= \left(+J_{AK} \right) / Vol_c \end{align*}

where VolcVol_c is the water volume of the cytosolic space in units of volume of cytosolic water space per unit cell volume. The fluxes in the right-hand sides of these expressions are defined below. The ratio VRm/VRcV_{Rm}/V_{Rc} is ratio of regional volume of the IMS to the cytosolic space. Since the JATPPERMJ_{ATPPERM}, JADPPERMJ_{ADPPERM}, JAMPPERMJ_{AMPPERM}, and JPIPERMJ_{PIPERM} fluxes are in units of mass per unit time per unit mitochondrial volume, the multiplication by VRm/VRcV_{Rm}/V_{Rc} converts the units to mass per unit time per unit cytosolic volume. The units of the other fluxes (cytosolic reaction fluxes) are mass per unit time per unit cytosolic volume.

2.6Cation Concentration State Variables:

The governing equations for the cation (H+^+, K+^+, and Mg2+^{2+}) concentrations in the mitochondrial and extra-mitochondrial compartments are derived using the method outlined in Vinnakota et al. (2009). In brief, the equations account for rapid equilibria between conjugate bases of biochemical weak acid species (e.g., ATP4^{4-}) and cation bound species (e.g., HATP3^{3-}, KATP3^{3-}, and MgATP2^{2-}). The full set of equations is detailed in the supplementary material published with Bazil et al. (2016).

2.7Energetic Model Fluxes:

The fluxes on the right-hand sides of Equations (1,2,3, and 4) are defined in Table 2.

Table 2:Energetics Model Reaction and Transport Fluxes

FluxDefinitionUnitsVariable name in code
JC1J_{C1}Electron transport chain Complex I fluxmole\cdotsec1^{-1}\cdot(L mito)1^{-1}J_ETC1_im_to_matrix
JC3J_{C3}Electron transport chain Complex III fluxmole\cdotsec1^{-1}\cdot(L mito)1^{-1}J_ETC3_im_to_matrix
JC4J_{C4}Electron transport chain Complex IV fluxmole\cdotsec1^{-1}\cdot(L mito)1^{-1}J_ETC4_im_to_matrix
JF1F0J_{F1F0}Mitochondrial F1_1F0_0 ATPase fluxmole\cdotsec1^{-1}\cdot(L mito)1^{-1}J_F1F0ATPASE_im_to_matrix
JANTJ_{ANT}Adenine nucleotide translocase fluxmole\cdotsec1^{-1}\cdot(L mito)1^{-1}J_ANT_im_to_matrix
JHleakJ_{Hleak}Proton leak fluxmole\cdotsec1^{-1}\cdot(L mito)1^{-1}J_HLEAK_im_to_matrix
JDHJ_{DH}Rate of NADH productionmole\cdotsec1^{-1}\cdot(L mito)1^{-1}J_DH_matrix
JPICJ_{PIC}Mitochondrial phosphate carrier fluxmole\cdotsec1^{-1}\cdot(L mito)1^{-1}J_PIC_im_to_matrix
JATPPERMJ_{ATPPERM}Mitochondrial outer membrane ATP permeationmole\cdotsec1^{-1}\cdot(L mito)1^{-1}J_ATPPERM_cytoplasm_to_im
JADPPERMJ_{ADPPERM}Mitochondrial outer membrane ADP permeationmole\cdotsec1^{-1}\cdot(L mito)1^{-1}J_ADPPERM_cytoplasm_to_im
JAMPPERMJ_{AMPPERM}Mitochondrial outer membrane AMP permeationmole\cdotsec1^{-1}\cdot(L mito)1^{-1}J_AMPPERM_cytoplasm_to_im
JPIPERMJ_{PIPERM}Mitochondrial outer membrane Pi permeationmole\cdotsec1^{-1}\cdot(L mito)1^{-1}J_PIPERM_cytoplasm_to_im
JCKJ_{CK}Cytosolic creatine kinase fluxmole\cdotsec1^{-1}\cdot(L cytosol)1^{-1}J_CK_cytoplasm
JAKJ_{AK}Cytosolic adenylate kinase fluxmole\cdotsec1^{-1}\cdot(L cytosol)1^{-1}J_AK_cytoplasm
JATPaseJ_{ATPase}Cytosolic ATP hydrolysis fluxmole\cdotsec1^{-1}\cdot(L cytosol)1^{-1}J_ATPASE_cytoplasm

The mathematical expressions for these fluxes are detailed in Bazil et al. (2016).

2.8Implementation in Multiscale Model:

The cellular energetics model is implemented in a MATLAB script called EnergeticsModelScript.m. This script is used to predict cytosolic [ATP], [ADP], [AMP], [Pi], [Cr], and [CrP] at a specified input rate of cytosolic ATP hydrolysis. Input and output arguments for the script are listed below in Tables 3 and 4.

Table 3:Input arguments for cellular energetics model

Input variableDefinitionUnits used in codeValues
TANtotal adenine nucleotide poolmole\cdot(l cell)1^{-1}0.0071–0.0086 for sham 0.0052–0.0085 for TAC
CRtottotal creatine poolmole\cdot(l cell)1^{-1}0.0267–0.0330 for sham 0.0146–0.0278 for TAC
TEPtotal exchangeable phosphate poolmole\cdot(l cell)1^{-1}0.0247–0.0298 for sham 0.0181–0.0293 for TAC
Ox_capacityoxidative capacity (relative to control)unitless0.834–1.1526 for sham 0.5287–0.9755 for TAC
x_ATPaseATP hydrolysis ratemmole\cdotsec1^{-1}\cdot(l cell)1^{-1}0.9301–1.392 for sham 0.6758–1.410 for TAC

Table 4:Output arguments for cellular energetics model

Output variableDefinitionUnits used in codeValues
MgATP_cytoplasmcytosolic [MgATP]mmole\cdot(l cytosol water)1^{-1}7.2145–9.4386 for sham 4.9384–9.0634 for TAC
MgADP_cytoplasmcytosolic [MgADP]mmole\cdot(l cytosol water)1^{-1}0.0467–0.0598 for sham 0.0274–0.0549 for TAC
fPi_cytoplasmcytosolic unchelated [Pi]mmole\cdot(l cytosol water)1^{-1}0.4393–1.3073 for sham 1.1437–1.6650 for TAC
MVO2_tissueoxygen consumption rateμ\mumole\cdotmin1^{-1}\cdot(g tissue)1^{-1}7.617–12.380 for sham 6.9104–12.450 for TAC
dGrATPaseATP hydrolysis free energykJ\cdotmole1^{-1}-(66.054–62.71) for sham -(63.574–61.44) for TAC
PCrATPCrP/ATP ratiounitless1.817–2.324 for sham 1.4985–1.897 for TAC
ATP_cytocytosolic total [ATP]mmole\cdot(l cytosol water)1^{-1}8.733–11.438 for sham 5.983–10.959 for TAC
ADP_cytocytosolic total [ADP]mmole\cdot(l cytosol water)1^{-1}0.1164–0.1486 for sham 0.0682–0.1359 for TAC
Pi_cytocytosolic total [Pi]mmole\cdot(l cytosol water)1^{-1}0.7403–2.197 for sham 1.9225–2.798 for TAC

3Cardiomyocyte Mechanics Model

3.1Model Variables and Equations:

A cardiomyocyte mechanics model based on the models of Tewari et al. (2016)Tewari et al. (2016) and Campbell et al. (2018) is used to simulate the active and passive components of myocardial wall tension used in the heart model (§3, below) and to determine the ATP hydrolysis rate used in the energy metabolism model (§1, above). The components of the model are illustrated in Figure 2.

The five states in the cross-bridge model correspond to: the non-permissible (no calcium bound) state NN, the permissible (calcium-bound) state PP, loosely attached state A1A_1, strongly attached state A2A_2, and post-ratcheted state A3A_3. The attached states are distributed over a continuum of cross-bridge strain. To numerically simulate the model a moment-expansion approach is used where ordinary differential equations for the first three moments of the probability distributions of strain of each of the attached states are simulated.

The state variables for the cross-bridge model are tabulated below.

Table 5:State variables in cross-bridge model

DefinitionUnits used in code
Variable name
in code
The 0th^{\rm th} moment of state A1 strain
probability distribution. Equal to the
proportion of cross-bridges in state A1.
The 1st^{\rm st} moment of state A1 strain
probability distribution.
The 2nd^{\rm nd} moment of state A1 strain
probability distribution.
The 0th^{\rm th} moment of state A2 strain
probability distribution. Equal to the
proportion of cross-bridges in state A2.
The 1st^{\rm st} moment of state A2 strain
probability distribution.
The 2nd^{\rm nd} moment of state A2 strain
probability distribution.
The 0th^{\rm th} moment of state A3 strain
probability distribution. Equal to the
proportion of cross-bridges in state A3.
The 1st^{\rm st} moment of state A3 strain
probability distribution.
The 2nd^{\rm nd} moment of state A3 strain
probability distribution.
Non-permissible XB state
Non relaxed state

The equations used to simulate the cross-bridge model are

dp10dt=kaP(t)UNROVthickkd~p10k1~((p10α1p11+1/2α12p12)+k1(p20+α1p21+12α12p22)dp11dt=vp10kd~p11k1~(p11α1p12)+k1(p21+α1p22)dp12dt=2vp11kd~p12k1~p12+k1p22dp20dt=k1~(p10α1p11+12α12p12)k1(p20+α1p21+12α12p22)k2(p20α2p21+12α22p22)+k2~p30dp21dt=vp20+k1~(p11α1p12)k1(p21+α1p22)k2(p21α2p22)+k2~p31dp22dt=2vp21+k1~p12k1p22k2p22+k2~p32dp30dt=k2(p20α2p21+12α22p22)k2~p30k3~(p30α3s32p30+2α3s3p31+p32)dp31dt=vp30+k2(p21α2p22)k2~p31k3~(p31α3s32p31+2α3s3p32)dp32dt=2vp31+k2p22k2~p32k3~(p32+α3s32p32)\begin{align*} \frac{dp_1^0}{dt} &= k_a P(t) U_{NR} OV_{thick} - \tilde{k_d} p_1^0 - \tilde{k_1}\left((p_1^0-\alpha_1 p_1^1+1/2 \alpha_1^2 p_1^2\right) + k_{-1} \left( p_2^0+\alpha_1 p_2^1+ \frac{1}{2}\alpha_1^2 p_2^2 \right) \nonumber \\ \frac{dp_1^1}{dt} &= v p_1^0 - \tilde{k_d} p_1^1 - \tilde{k_1} \left( p_1^1-\alpha_1 p_1^2 \right) + k_{-1} \left(p_2^1 + \alpha_1 p_2^2 \right) \nonumber \\ \frac{dp_1^2}{dt} &= 2vp_1^1 - \tilde{k_d} p_1^2 - \tilde{k_1} p_1^2 + k_{-1} p_2^2 \nonumber \\ \frac{dp_2^0}{dt} &= \tilde{k_1}\left(p_1^0 - \alpha_1 p_1^1 + \frac{1}{2} \alpha_1^2 p_1^2 \right) - k_{-1} \left(p_2^0 + \alpha_1 p_2^1 + \frac{1}{2} \alpha_1^2 p_2^2\right) - k_2 \left(p_2^0 - \alpha_2 p_2^1 + \frac{1}{2} \alpha_2^2 p_2^2 \right) + \widetilde{k_{-2}}p_3^0 \nonumber \\ \frac{dp_2^1}{dt} &= v p_2^0 +\tilde{k_1}\left(p_1^1-\alpha_1 p_1^2 \right) - k_{-1}\left(p_2^1 + \alpha_1 p_2^2 \right) - k_2 \left(p_2^1-\alpha_2 p_2^2 \right) + \widetilde{k_{-2}}p_3^1 \nonumber \\ \frac{dp_2^2}{dt} &= 2v p_2^1 +\tilde{k_1}p_1^2 - k_{-1} p_2^2 - k_2 p_2^2+ \widetilde{k_{-2}}p_3^2 \nonumber \\ \frac{dp_3^0}{dt} &= k_2 \left(p_2^0 - \alpha_2 p_2^1 + \frac{1}{2}\alpha_2^2 p_2^2 \right) - \widetilde{k_{-2}}p_3^0 - \tilde{k_3}\left(p_3^0 - \alpha_3 s_3^2 p_3^0 + 2\alpha_3 s_3 p_3^1+p_3^2 \right) \nonumber \\ \frac{dp_3^1}{dt} &= v p_3^0 + k_2 \left(p_2^1 - \alpha_2 p_2^2 \right) - \widetilde{k_{-2}}p_3^1 - \tilde{k_3}\left(p_3^1 - \alpha_3 s_3^2 p_3^1 + 2\alpha_3 s_3 p_3^2 \right) \nonumber \\ \frac{dp_3^2}{dt} &= 2v p_3^1+k_2 p_2^2- \widetilde{k_{-2}}p_3^2 - \tilde{k_3}\left(p_3^2+\alpha_3 s_3^2 p_3^2 \right) \end{align*}

where v is the velocity of sliding (v=dSL/dtv=dSL/dt where SLSL is the sarcomere length, used below in the heart model of Section 4).

Metabolite concentrations affect the apparent rate constants in the model via the following relations:

kd~=kd[Pi]/KPi1+[Pi]/KPik1~=k111+[Pi]/KPik2~=k2[MgADP]/KMgADP1+[MgADP]/KMgADP+[MgATP]/KMgATPk3~=k3[MgATP]/KMgATP1+[MgADP]/KMgADP+[MgATP]/KMgATP\begin{align*} \tilde{k_d} &= k_d \frac{ [{\rm Pi}]/K_{Pi} }{1+ [{\rm Pi}]/K_{Pi} } \nonumber \\ \tilde{k_1} &= k_1 \frac{1}{1+ [{\rm Pi}]/K_{Pi} } \nonumber \\ \widetilde{k_{-2} } &= k_{-2} \frac{[{\rm MgADP}]/K_{MgADP} } {1+ [{\rm MgADP}]/K_{MgADP} + [{\rm MgATP}]/K_{MgATP} } \nonumber \\ \tilde{k_3} &= k_3 \frac{[{\rm MgATP}]/K_{MgATP} } {1+ [{\rm MgADP}]/K_{MgADP} + [{\rm MgATP}]/K_{MgATP} } \end{align*}

A detailed description of the moment expansion and associated equations is given in Tewari et al. (2016).

Cardiomyocyte mechanics model. The multi-scale strain-dependent model for the cross-bridge cycle is illustrated in panel A. The integration of the cross-bridge force (\sigma_{XB}) into a model of muscle mechanics is illustrated in panel B.

Figure 2:Cardiomyocyte mechanics model. The multi-scale strain-dependent model for the cross-bridge cycle is illustrated in panel A. The integration of the cross-bridge force (σXB\sigma_{XB}) into a model of muscle mechanics is illustrated in panel B.

3.2Calcium activation:

The calcium activation model is adopted from Campbell et al. (2018) model with minor modifications. The equations for calcium-mediated transition from the N to the P state are:

Jon=kon[Ca2+]N(1+kcoop(1N))Joff=koffP(1+kcoopN)\begin{align*} J_{on} &= k_{on} [{\rm Ca}^{2+} ] N \left(1 + k_{coop} (1-N) \right) \nonumber \\ J_{off} &= k_{off} P \left( 1 + k_{coop} N \right) \end{align*}

The term kcoop(1N)k_{coop} (1-N) is representative of cooperative activation. The variable N represents the non-permissible state:

P=1Np10p20p30dNdt=Jon+Joff\begin{align*} P &= 1 - N - p_1^0 - p_2^0 - p_3^0 \nonumber \\ \frac{dN}{dt} &= -J_{on} + J_{off} \end{align*}

where PP is the permissible (calcium-bound).

3.3Super-relaxed state:

The Campbell et al. model for calcium activation includes a transition between a super-relaxed and not relaxed state.

USRUNRUSR+UNR=1\begin{align*} U_{SR} \leftrightharpoons U_NR \nonumber \\ U_{SR} + U_{NR} = 1 \end{align*}

where the transition from super-relaxed (USRU_{SR}) to non-relaxed (UNRU_{NR}) state is force-dependent:

dUNRdt=kSR(1+kforceσXB)USRkSRUNR.\frac{dU_{NR}}{dt} = k_{SR} \left( 1 + k_{force} \sigma_{XB} \right)U_{SR} - k_{-SR} U_{NR} \, .

where σXB\sigma_{XB}, is active contractile force and formulated in equation 15.

3.4Overlap function:

Following Rice et al. (2008) the fractional overlap between thin and thick filament is represented as follows:

OVZaxis=min(Lthick/2,SL/2),OV_{Z-axis} = min\left(L_{thick}/2,SL/2\right) \, ,

where OVZaxisOV_{Z-axis} is the overlap region closest to the Z-axis,

OVMline=max(SL/2(SLLthin),Lbare/2),OV_{M-line} = max\left(SL/ 2 - (SL - L_{thin} \right), L_{bare}/2) \, ,

where OVMlineOV_{M-line} is the overlap region closest to the M-line. The length of overlap LOVLOV is computed as following:

LOV=OVZaxisOVMline.LOV = OV_{Z-axis} - OV_{M-line} \, .

Using length of overlap LOV, fraction of thick filament overlap is computed as following:

OVthick=2LOV/(LthickLbare).\begin{align*} OV_{thick} = 2 LOV / (L_{thick}- L_{bare} ) \, . \end{align*}

Here SLSL is the length of sarcomere, LthickL_{thick} , LthinL_{thin} , LbareL_{bare} are the length of thick filament, bare region of the thick filament and, the length of the thin filament.

Table 6:Sarcomere overlap function parameters

ParameterDefinitionValue and unitsParameter name in code
LthinL_thinthin filament length1200 nmL_thin
LthickL_thickthick filament length1670 nmL_thick
LbareL_barebare length of the thick filament100 nmL_hbare
OVZaxisOV_{Z-axis}overlap region closest to the Z-axisnmsovr_ze
OVMlineOV_{M-line}overlap region closest to the M-linenmsovr_cle
LOVLOVlength of overlapnmL_sovr
OVthickOV_{thick}fraction of thick filament overlapunitlessN_overlap

3.5Active and passive force:

The active force generated by cross-bridges is computed from contributions from pre- and post-ratcheted states:

σXB(t)=OVthick(kstiff,1(p22+p32)+kstiff,2Δrp30)\sigma_{XB} (t) = OV_{thick} ( k_{stiff,1} \left(p_2^2 + p_3^2 \right) + k_{stiff,2} \Delta r p_3^0 )

where kstiff,1k_{stiff,1} and kstiff,2)k_{stiff,2)} are stiffness constants, Δr\Delta r is the cross bridge strain associated with ratcheting deformation. The full muscle model (Figure 2) includes contributions from the active force generated by the cross-bridge mechanics, the viscous and passive forces associated with the muscle, F1F_1 and F2F_2, and a series element force FSEF_{SE}. Overall force balance for the model yields

σSE(t)=σXB(t)+σ1(t)+σ2(t).\sigma_{SE} (t) = \sigma_{XB} (t) + \sigma_1 (t) + \sigma_2 (t) \, .

The stress contributed from the dashpot (viscous) is determined from the rate of change of sarcomere length.

σ1=ηdSLdt.\sigma_1 = \eta \frac{dSL}{dt} \, .

The passive force σ2\sigma_2 is a function of sarcomere length and is calculated

σ2(SL)=kpassive(SLSLrest)+σPassivecollagen(SL)\sigma_2 (SL)= k_{passive} \left(SL - SL_{rest} \right) + \sigma_{Passive_collagen} (SL)

where kpassivek_{passive} is the stiffness parameters for the passive force, SLrestSL_{rest} is the sarcomere rest length, and the passive force exerted by collagen is adopted from Rice et al. (2008).

σPassivecollagen(SL)={Pconcollagen[ePExpcollagen(SLSLcollagen)1]ifSL>SLcollagen0otherwise.\begin{align*} \sigma_{Passive_collagen} (SL) = \left\{ \begin{array}{ll} Pcon_{collagen} \left[ e^{PExp_{collagen}(SL - SL_{collagen})} -1 \right] & if SL > SL_{collagen} \\ 0 & {\rm otherwise} \end{array} \right. \, . \end{align*}

where Pconcollagen=0.01Pcon_{collagen} = 0.01 (unitless) is the scale factor for passive force contributed by collagen. PExpcollagen=70PExp_{collagen} = 70 μ\mum1^{-1} is a model parameter and SLcollagen=2.25SL_{collagen} = 2.25 μ\mum is the minimum length threshold at which collagen starts to exert force on sarcomere.

3.6Model Parameters:

Certain parameters from the cross-bridge and calcium-activation models of Tewari et al. (2016)Tewari et al. (2016) and Campbell et al. (2018) were re-estimated to match data from Janssen et al. (2002) on calcium transients and force-generation in isolated rat cardiac trabeculae. In brief, experiments were conducted at 37^\circ C. Calcium transients (Figure 3) were measured at different stimulation frequencies at fixed sarcomere length SL = 2.2 μ\mum. Isometric tension time courses were measured at different stimulation frequencies and sarcomere lengths. Figure 3 shows data on peak developed tension (TdevT_{dev}) as a function of SLSL at stimulation frequency of 4 Hz; and data on relaxation time from peak to 50% of peak tension (RT50RT50); peak developed tension (TdevT_{dev}) and time to peak tension (TTPTTP) as function of SLSL.

Model simulations were fit to these data to estimate unknown parameters in the calcium activation and cross-bridge kinetics components of the model. Specifically, parameters adjusted to values different from those in the original publication are indicated below in Table 7.

Crossbridge model parameters estimation. For SL =1.9 \mum the Ca_{50} = 5.89 and Hill coefficient n = 4.63 and for SL = 2.3 \mum the Ca_{50} =6.001 and Hill coefficient n = 4.47. Error bars shown in panel (B) represent standard error from the n = 9 data-set .

Figure 3:Crossbridge model parameters estimation. For SL=1.9SL =1.9μ\mum the Ca50=5.89Ca_{50} = 5.89 and Hill coefficient n=4.63n = 4.63 and for SL=2.3SL = 2.3μ\mum the Ca50=6.001Ca_{50} =6.001 and Hill coefficient n=4.47n = 4.47. Error bars shown in panel (B) represent standard error from the n = 9 data-set Janssen et al., 2002.

Table 7:Model parameters for cross bridge model

ParameterDefinitionValue and unitsParameter name in codeReference
kstiff,1k_{stiff,1}Stiffness constant due to myosin–actin interaction13013 mmHg/μ\mumkstiff1Fit to data in Figure 3
kstiff,2k_{stiff,2}Stiffness constant due to working stroke of XBs341590 mmHg/μ\mumkstiff2Fit to data in Figure 3
kpassivek_{passive}Passive stiffness constant25 mmHg/μ\mumk_passiveFit to data in Figure 3
SLrestSL_{rest}Sarcomere length at which passive force is zero1.8 μ\mumSL_rest_pasCampbell et al. (2018)
α1\alpha_1Strain-dependency parameter10 μ\mum1^{-1}alpha1Tewari et al. (2016)
α2\alpha_2Strain-dependency parameter9 μ\mum1^{-1}alpha2Tewari et al. (2016)
α3\alpha_3Strain-dependency parameter5.93 μ\mum2^{-2}alpha3Fit to data in Figure 3
s3s_3Strain-dependency parameter9.9×1039.9 \times 10^{-3} μ\mums3Tewari et al. (2016)
kcoopk_{coop}Strength of thin filament cooperativity1.857K_coopFit to data in Figure 3
konk_{on}Rate constant of Ca binding to troponin C101.2 μ\muM1^{-1}\cdots1^{-1}k_onFit to data in Figure 3
koffk_{off}Rate constant of Ca unbinding from troponin C101.2 s1^{-1}k_offFit to data in Figure 3
kforcek_{force}Force-dependent rate constant of super relax transition1.169×1031.169 \times 10^{-3} N1^{-1}\cdotm2^{-2}kforceFit to data in Figure 3
kSRk_{SR}Rate constant of force-dependent super relax transition14.44 s1^{-1}ksrFit to data in Figure 3
kSRk_{-SR}Rate constant of force-dependent super relax transition50.03 s1^{-1}kmsrFit to data in Figure 3
kMgATPk_{MgATP}[MgATP] dissociation constant489.7 μ\muMK_TTewari et al. (2016)
kMgADPk_{MgADP}[MgADP] dissociation constant194.0 μ\muMK_DTewari et al. (2016)
kPik_{Pi}[Pi] dissociation constant4.0 mMK_PiTewari et al. (2016)
kak_aMyosin–actin rate of attachment559.6 s1^{-1}kaFit to data in Figure 3
kdk_dMyosin–actin rate of un-attachment304.7 s1^{-1}kdFit to data in Figure 3
k1k_1Transition rate constant112.4 s1^{-1}k1Fit to data in Figure 3
k1k_{-1}Rate of strongly-bound to weakly-bound transition21.30 s1^{-1}km1Tewari et al. (2016)
k2k_2Rate of ratcheting811.7 s1^{-1}k2Tewari et al. (2016)
k2k_{-2}Reverse ratcheting rate43.25 s1^{-1}km2Tewari et al. (2016)
k3k_3Myosin–actin detachment rate144.6 s1^{-1}k3Fit to data in Figure 3

4Heart Model

4.1Model Variables and Equations:

A modified version of the Lumens et al. (2009) TriSeg model is used to simulate left- and right-ventricular mechanics, based on the implementation of Tewari et al. (2016). Tension development in each of the left-ventricular free wall, septum, and right-ventricular free wall is simulated using a cell mechanics model to represent each of these segments. From Eq. (17) the rates of change of sarcomere length in these three segments is given by

dSLRVdt=σSE,RVσ2(SLRV)λXBσXB,RV(t)ηdSLLVdt=σSE,LVσ2(SLLV)λXBσXB,LV(t)ηdSLSEPdt=σSE,SEPσ2(SLSEP)λXBσXB,SEP(t)η\begin{align*} \frac{dSL_{RV}}{dt} &= \frac{\sigma_{SE,RV} - \sigma_{2}(SL_{RV}) - \lambda_{XB}\sigma_{XB,RV}(t)}{\eta} \nonumber \\ \frac{dSL_{LV}}{dt} &= \frac{\sigma_{SE,LV} - \sigma_{2}(SL_{LV}) - \lambda_{XB}\sigma_{XB,LV}(t)}{\eta} \nonumber \\ \frac{dSL_{SEP}}{dt} &= \frac{\sigma_{SE,SEP} - \sigma_{2}(SL_{SEP}) - \lambda_{XB}\sigma_{XB,SEP}(t)}{\eta} \end{align*}

where the parameter λXB\lambda_{XB} is a scalar used to account for differences in force generation in vivo versus in vitro. (The value value λXB\lambda_{XB} = 1.4, determined by Tewari et al. (2016) accounts for slightly lower force generation in vitro versus in vivo.)

The SLSL and dSLdt\frac{dSL}{dt} computed from Eq. (20) are used in Eqs. (5) which govern cross-bridge dynamics in each segment. The dynamical state of the cross-bridge model in each segment, in turn, appears in Eq. (20), which governs SL(t)SL(t) for each segment.

The series element elastic force for each segment is computed to be proportional to the difference between the sarcomere length and the sarcomere length calculated from natural myofiber strain:

σSE,RV=KSE(SL0,RVSLRV)σSE,LV=KSE(SL0,LVSLLV)σSE,SEP=KSE(SL0,SEPSLSEP)\begin{align*} \sigma_{SE,RV} &= K_{SE}(SL_{0,RV} - SL_{RV}) \nonumber \\ \sigma_{SE,LV} &= K_{SE}(SL_{0,LV} - SL_{LV}) \nonumber \\ \sigma_{SE,SEP} &= K_{SE}(SL_{0,SEP} - SL_{SEP}) \nonumber \end{align*}


SL0,#=SLrefexp(ϵf,#))ϵf,#=12ln(Am,#Am,ref,#)112z#20.019z#4z#=3Cm,#VW,#2Am,#Vm,#=π6xm,#(xm,#2+3ym2)Am,#=π(xm,#2+ym2)Cm,#=2xm,#xm,#2+ym2.\begin{align*} SL_{0,\#} &= SL_{ref} exp(\epsilon_{f,\#}) ) \nonumber \\ \epsilon_{f,\#} &= \frac{1}{2} ln (\frac{A_{m,\#}}{A_{m,ref,\#}}) - \frac{1}{12} z_{\#}^{2} - 0.019 z_{\#}^{4} \nonumber \\ z_{\#} &= \frac{3C_{m,\#}V_{W,\#}}{2A_{m,\#}} \nonumber \\ V_{m,\#} &= \frac{\pi}{6} x_{m,\#} (x_{m,\#}^2 + 3y_{m}^{2}) \nonumber \\ A_{m,\#} &= \pi(x_{m,\#}^2 + y_{m}^{2}) \nonumber \\ C_{m,\#} &= \frac{2x_{m,\#}}{x_{m,\#}^{2}+y_{m}^{2}} \, . \nonumber \end{align*}

Here Am,#A_{m,\#} the midwall surface area of segment #, Am,ref,#A_{m,ref,\#} is a reference midwall surface area, Cm,#C_{m,\#} is the curvature of the midwall surface, Vw,#V_{w,\#} is the wall volume of wall segment i, Vm,#V_{m,\#} is the midwall volume, and xm,#x_{m,\#} and ymy_{m} determine the geometry of the LV and RV cavity (see Lumens et al. (2009)). The four variables of the TriSeg heart model, xm,RVx_{m,RV} , xm,LVx_{m,LV}, xm,SEPx_{m,SEP}, and ymy_{m} that determine the geometry of the ventricular cavities, are listed in Table 8.

Table 8:State variables in TriSeg (Heart) model

State VariableDefinitionUnitsVariable name in code
Maximal axial distance
from RV midwall surface to origin
Maximal axial distance
from LV midwall surface to origin
Maximal axial distance
from SEP midwall surface to origin
ym_mRadius of midwall junction circlecmym

For given wall volumes and ventricular volumes, the geometry of the heart is solved such that equilibrium of radial and axial tensile forces is achieved at the junction margin (i.e., where the three wall segments meet forming ventricular cavities).

Tension in the midwall of each segment is calculated as a function of stress:

Tm,#=VW,#σSE,#2Am,#(1+z23+z45)T_{m,\#} = \frac {V_{W,\#} \sigma_{SE,\#}}{2A_{m,\#}} \left(1+\frac{z^{2}}{3}+\frac{z^{4}}{5}\right)

Axial and radial components of the tension are computed

Tx,#=Tm,#2xm,#ymxm,#2+ym2Ty,#=Tm,#xm,#2+ym2xm,#2+ym2\begin{align*} T_{x,\#} = T_{m,\#} \frac{2x_{m,\#}y_{m}}{x_{m,\#}^{2}+y_{m}^{2}} \nonumber \\ T_{y,\#} = T_{m,\#} \frac{-x_{m,\#}^{2}+y_{m}^{2}}{x_{m,\#}^{2}+y_{m}^{2}} \end{align*}

The four unknowns of the model—xm,RVx_{m,RV} , xm,LVx_{m,LV}, xm,SEPx_{m,SEP}, and ymy_{m}—are determined by satisfying the relations

Vm,LV=VLV12Vw,LV12Vw,SEP+Vm,SEPVm,RV=+VRV12Vw,RV12Vw,SEPVm,SEPTx,RV+Tx,LV+Tx,SEP=0Ty,RV+Ty,LV+Ty,SEP=0.\begin{align*} & V_{m,LV} = -V_{LV} - \frac{1}{2} V_{w,LV} - \frac{1}{2} V_{w,SEP} + V_{m,SEP} \nonumber \\ & V_{m,RV} = +V_{RV} - \frac{1}{2} V_{w,RV} - \frac{1}{2} V_{w,SEP} - V_{m,SEP} \nonumber \\ & T_{x,RV} + T_{x,LV} + T_{x,SEP} = 0 \nonumber \\ & T_{y,RV} + T_{y,LV} + T_{y,SEP} = 0 \nonumber \, . \end{align*}

Transmural pressures are computed

Ptrans,#=2Tx,#ymP_{trans,\#} = \frac{2T_{x,\#}}{y_{m}}

and the pressures in the cavities are computed

PLV=Ptrans,LVPRV=+Ptrans,RV.\begin{align*} P_{LV} &= - P_{trans,LV} \nonumber \\ P_{RV} &= + P_{trans,RV} \nonumber \, . \end{align*}

4.2Model Parameters

Parameters defining the mass and geometry of the heart are identified from data on individual animals are defined in Table 9.

Table 9:Parameters in TriSeg (Heart) model

ParameterDefinitionValue and UnitsVariable name in code
Vw,LVV_{w,LV}LV wall volume
(experimental data), mL
Vw,SEPV_{w,SEP}Septal wall volume
(experimental data), mL
Vw,RVV_{w,RV}RV wall volume
(experimental data), mL
LV midwall reference
surface area
adjustable, cm2^{2}Amref_LV
Septal midwall reference
surface area
adjustable, cm2^{2}Amref_SEP
RV midwall reference
surface area
adjustable, cm2^{2}Amref_RV
KSEK_{SE}Stiffness of series element50000 mmHg/μ\mumKse
η\etaViscosity coefficient of myofibers1 mmHg\cdotsecμ\cdot \mum1^{-1}eta

5Lumped-Parameter Cardiovascular Systems Model

5.1Model Variables and Equations:

The lumped-parameter model illustrated in Figure 4 is used to simulate pressures and flows in the systemic and pulmonary circuits. This simple lumped model invokes eight parameters representing: pulmonary resistance RpulR_{pul}, pulmonary arterial and venous compliances CPAC_{PA} and CPVC_{PV}, systemic arterials resistances RAoR_{Ao} and RsysR_{sys}, systemic arterial compliance CSAC_{SA}, and systemic venous compliance CSVC_{SV}, systemic arterial resistance and aortic compliance CAoC_{Ao}.

Cardiovascular system (CVS) diagram. Adopted from . C_{PA}, C_{PV}, C_{SA}, C_{Ao}, and C_{SV} represent lumped compliances of pulmonary arteries, pulmonary veins, systemic arteries, aorta, and systemic veins. R_{pul}, R_{sys}, and R_{Ao} represent vascular resistances. The lumped-parameter representations of the systemic and pulmonary circuits are coupled to the TriSeg heart model of , described below.

Figure 4:Cardiovascular system (CVS) diagram. Adopted from Tewari et al. (2016)Tewari et al. (2016). CPAC_{PA}, CPVC_{PV}, CSAC_{SA}, CAoC_{Ao}, and CSVC_{SV} represent lumped compliances of pulmonary arteries, pulmonary veins, systemic arteries, aorta, and systemic veins. RpulR_{pul}, RsysR_{sys}, and RAoR_{Ao} represent vascular resistances. The lumped-parameter representations of the systemic and pulmonary circuits are coupled to the TriSeg heart model of Lumens et al. (2009), described below.

Flow through a resistive element is calculated

q=P1P2Rq =\frac{P_{1} - P_{2}}{R}

where P1P2P_{1} - P_{2} is the pressure drop across the element, and R is the resistance. Pressure in a compliant/capacitive element is governed by

dPdt=FinFoutC\frac{dP}{dt} =\frac{F_{in} - F_{out}}{C}

where FinFoutF_{in} - F_{out} is the rate of change of blood volume in the element. Table 10 lists the variables of the cardiovascular systems model. These variables and Eqs.  28 and  29 are invoked in the MATLAB code Cardiovascularmechancis.m.

Table 10:Variables used in lumped-parameter cardiovascular model

ParameterDefinitionUnits used in codeVariable name in code
Cardiovascular system model state variables
VLVV_{LV}Left ventricle volumemLV_LV
VRVV_{RV}Right ventricle volumemLV_RV
VSVV_{SV}Volume of systemic veinmLV_SV
VPVV_{PV}Volume of pulmonary veinmLV_PV
VSAV_{SA}Volume of systemic arterymLV_SA
VPAV_{PA}Volume of pulmonary arterymLV_PA
VAoV_{Ao}Volume of aortamLV_Ao
Pressures computed from volume state variables
PSVP_{SV}Systemic venous pressuremmHgP_SV
PPVP_{PV}Pulmonary venous pressuremmHgP_PV
PPAP_{PA}Pulmonary arterial pressuremmHgP_PA
PAoP_{Ao}Aortic pressure (proximal to TAC)mmHgP_Ao
PSAP_{SA}Systemic arterial pressure (distal to TAC)mmHgP_SA

5.2Identification of Model Parameters:

Table 11 lists the parameters used in the lumped circulatory model.

Table 11:Parameter used in the CVS model

ParameterDefinitionValue and Units
Variable name
in code
CAoC_{Ao}Proximal aortic compliance0.0022045 mL\cdotmmHg1^{-1}C_Ao
CSAC_{SA}Systemic arterial compliance0.0077157 mL\cdotmmHg1^{-1}C_SA
CSVC_{SV}Systemic venous compliance2.5 mL\cdotmmHg1^{-1}C_SV
CPAC_{PA}Pulmonary arterial compliance0.013778 mL\cdotmmHg1^{-1}C_PA
CPVC_{PV}Pulmonary venous compliance0.25 mL\cdotmmHg1^{-1}C_PV
RAoR_{Ao}Proximal aortic resistance2.5 mmHg\cdotsec\cdot mL1^{-1}R_Ao
RsysR_{sys}Systemic vasculature resistanceadjustable, mmHg\cdotsec\cdot mL1^{-1}R_sys
RSVR_{SV}Systemic veins resistance0.25 mmHg\cdotsec\cdot mL1^{-1}R_SV
RTAoR_{TAo}Transmural aortic resistance0.5 mmHg\cdotsec\cdot mL1^{-1}R_tAo
RTSAR_{TSA}Transmural systemic artery resistance4 mmHg\cdotsec\cdot mL1^{-1}R_tSA
RPAR_{PA}Pulmonary vasculature resistance7.58 mmHg\cdotsec\cdot mL1^{-1}R_PA
RPVR_{PV}Pulmonary veins resistance0.25 mmHg\cdotsec\cdot mL1^{-1}R_PV
RVLVR_{VLV}Valve resistance0.05 mmHg\cdotsec\cdot mL1^{-1}R_VLV
RTACR_{TAC}Resistance of TACadjustable, mmHg\cdotsec\cdot mL1^{-1}R_TAC

The systemic compliances (CSAC_{SA},CAoC_{Ao}, CSVC_{SV}) are fixed to produce a pulse pressure of roughly 33 mmHg for simulations of sham control rats. The pulmonary compliance (CPAC_{PA}, CPVC_{PV}) are fixed to roughly have the target value of 12 mmHg for the pulmonary pulse pressure. The resistance RAoR_{Ao} is arbitrarily set to have a small pressure drop of 4 mmHg between the aorta and systemic arteries for cardiac output of the mean value of 95 mL per min. The systemic venous resistance RSVR_{SV} is set so that the mean pressure in the systemic veins for sham control rats is 3 mmHg. The pulmonary resistances RPAR_{PA} and RPVR_{PV} are set to give a mean pulmonary arterial pressure of 21 mmHg and pulmonary venous pressure of 9 mmHg in the sham control rats.

The two parameters in the circulation model that are adjusted to match measured data on individual TAC and sham rats are RTACR_{TAC} and RsysR_{sys}. The resistance RTACR_{TAC} represents the resistance across the transverse aortic constriction (TAC), and is set to zero in sham-operated rats. In TAC rats, the value of this resistance is obtained based on the pressure gradient across the TAC constriction estimated from ultrasound measurements of the velocity gradient. The pressure drop across the constriction is computed ΔPTAC=12ρ(V22V12)\Delta P_{TAC} = \frac{1}{2} \rho (V_{2}^{2} - V_{1}^{2}), where V1V_{1} and V2V_{2} are the velocities on either side of the TAC constriction. Given an estimated pressure drop the resistance is computed RTAC=ΔPTACCOR_{TAC} = \frac{\Delta P_{TAC}}{CO}, where CO is the cardiac output. The systemic resistance RsysR_{sys} is adjusted so that the mean arterial pressure (MAP) is maintained at 93.3 mmHg. (Model fitting procedures are detailed below in Section 6.)

6Model fits and predictions associated with individual rats


The cardiac energy metabolism and the whole-body cardiovascular mechanics model (which includes the heart model) are implemented as separate modules. These models are matched to data on an individual animal basis. The cardiac energetics model takes as an input the myocardial ATP consumption rate, the measured metabolite pools levels, and the measured mitochondrial ATP synthesis capacity, and outputs the cytoplasm concentrations of phosphate metabolites. The cardiac mechanics code takes as an input the cytoplasmic concentrations of phosphate metabolites (namely, ATP, ADP, and Pi) and computes as an output the ventricular end-systolic and end-diastolic volumes and arterial pressures to compare to measured data, and the myocardial ATP hydrolysis rate to use in the energetics module. The energetics and mechanics models are iteratively run until they simultaneously converge to a steady state at fit the target cardiovascular data.

Table 12:Input arguments and adjustable parameters for cardiovascular mechanics model

ParameterVariable name in codeDefinitionUnits used in codeValues
HRHRHeart ratebpm318.6–367 for sham 294–349 for TAC. Set to measured value
Vw,LVVw_LVLV wall volumemL0.564–0.853 for sham 0.863–1.22 for TAC. Set as 2/3 of measured LV volume
Vw,SEPVw_SEPSeptal wall volumemL0.2821–0.4266 for sham 0.431–0.61 for TAC. Set as 1/3 of measured LV volume
Vw,RVVw_RVRV wall volumemL0.282–0.373 for sham 0.249–0.605 for TAC. Set as measured RV free wall volume
Am,refLWA_{m,ref}^{LW}Amref_LVLV midwall reference surface areacm2^{2}1.6865–2.465 for SHAM 2.299–2.966 for TAC. Adjusted to fit data
Am,refSEPA_{m,ref}^{SEP}Amref_SEPSeptal midwall reference surface areacm2^{2}1.055–1.487 for sham 1.299–1.620 for TAC. Adjusted to fit data
Am,refRWA_{m,ref}^{RW}Amref_RVRV midwall reference surface areacm2^{2}2.629–3.612 for sham 3.135–5.538 for TAC. Adjusted to fit data
kforcek_{force}k_forceModel parameter for force dependent super relax transitionN1^{-1} m2^{2}(0.905–3.312)×103\times 10^{-3} for sham (1.464–3.783)×103\times 10^{-3} for TAC. Adjusted to fit data
kSRk_{SR}ksrOn rate constant for super relax states1^{-1}9.021–33.013 for sham 14.598–37.708 for TAC. Adjusted to fit data
RsysR_{sys}R_sysSystemic vasculature resistancemmHg\cdotsec\cdot mL1^{-1}38.071–71.886 for sham 50.577–125.33 for TAC. Adjusted to fit data
RTACR_{TAC}R_TACresistance of TACmmHg\cdotsec\cdot mL1^{-1}0 for sham 7.11–22.25 for TAC. Adjusted to fit data
MgATP_cyto-plasmMgATPcytosolic [MgATP]mmole\cdot(L cytosol water)1^{-1}7.214–9.439 for sham 4.938–9.063 for TAC . Predicted from energetics model
MgADP_cyto-plasmMgADPcytosolic [MgADP]mmole\cdot(L cytosol water)1^{-1}0.047–0.06 for sham 0.027–0.055 for TAC. Predicted from energetics model
Pi_cyto-plasmPicytosolic total [Pi]mmole\cdot(L cytosol water)1^{-1}0.740–2.197 for sham 1.922–2.798 for TAC. Predicted from energetics model

6.2Relationship between cross-bridge cycle and ATP hydrolysis rates:

The relationship between cross-bridge cycle rate and J_ATP, the rate myocardial oxygen consumption rate, is based on matching the myocardial oxygen consumption rate predicted by the model to that observed for the working rat heart. Duvelleroy et al. (1976) report a mean oxygen consumption rate of approximately MVO2=0.31MVO2 = 0.31 mL\cdot(min\cdotg)1^{-1} for work rates corresponding to resting state in blood perfused working hearts. Using a myocardial cell density (in terms of cardiomyocyte volume per unit mass of myocardium) of ρcell=0.694\rho_{cell} = 0.694 mL\cdotg1^{-1} and assuming a P/O2 ratio (moles of ATP synthesized per mole of O2 consumed) of 4.5 Wu et al., 2008, we estimate a resting ATP hydrolysis rate of 1.301.30 mmol\cdot(sec\cdotL cell)1^{-1}. (This value is approximately 2.5 times higher than the estimated resting ATP hydrolysis rate of to 0.5470.547 mmol\cdot(sec\cdotL cell)1^{-1} for human myocardium Gao et al., 2019.)

Simulations of the cardiovascular mechanics model predict an average cross-bridge cycling rate of 5.03 sec1^{-1} for the mean sham-operated control rat. Assuming a fixed proportionality between cross-bridge cycling rate and myocardial ATP hydrolysis rate, yields a constant of proportionality of

ATPhydrolysisrate(mmol(secLcell)1)=0.258(mmol(Lcell)1)×crossbridgecyclingrate(sec1).\begin{align*} & & {\rm ATP\, hydrolysis\, rate (mmol\cdot(sec \cdot L \, cell)^{-1})} = \nonumber \\ & & {\rm 0.258 \, (mmol \cdot (L \, cell)^{-1}) \times crossbridge \, cycling \,rate \, (sec^{-1}) } \, . \end{align*}