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

Abstract¶

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¶

1Introduction¶

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 Mg $^{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 O $_2^{\cdot -}$ and H $_2$ O $_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 Variable	Definition	Units	Variable name in code
$\Delta\Psi$	Mitochondrial membrane potential	mV	DPsi_im_to_matrix
Mitochondrial Matrix State Variables
[ATP] $_x$	Total matrix ATP concentration	M	ATP_matrix
[ADP] $_x$	Total matrix ADP concentration	M	ADP_matrix
[Pi] $_x$	Total matrix Pi concentration	M	Pi_matrix
[NADH] $_x$	Total matrix NADH concentration	M	NADH_matrix
[NAD] $_x$	Total matrix NAD concentration	M	NAD_matrix
[UQH2] $_x$	Total matrix ubiquinol concentration	M	coQH2_matrix
[UQ] $_x$	Total matrix ubiquinone concentration	M	coQ_matrix
[H $^+$ ] $_x$	Matrix free proton concentration	M	h_matrix
[K $^+$ ] $_x$	Matrix free potassium concentration	M	k_matrix
[Mg $^{2+}$ ] $_x$	Matrix free magnesium concentration	M	m_matrix
Intermembrane Space (IMS) State Variables
[c $^{2+}$ ] $_i$	Total cytochrome c $^{2+}$ (reduced) concentration	M	cytocred_im
[c $^{3+}$ ] $_i$	Total cytochrome c $^{3+}$ (oxidized) concentration	M	cytocox_im
[ATP] $_i$	Total IMS ATP concentration	M	ATP_matrix
[ADP] $_i$	Total IMS ADP concentration	M	ADP_matrix
[AMP] $_i$	Total IMS AMP concentration	M	AMP_matrix
[Pi] $_i$	Total matrix Pi concentration	M	Pi_im
[H $^+$ ] $_i$	IMS free proton concentration	M	h_im
[K $^+$ ] $_i$	IMS free potassium concentration	M	k_im
[Mg $^{2+}$ ] $_i$	IMS free magnesium concentration	M	m_im
Cytosolic Space State Variables
[ATP] $_c$	Total cytosolic ATP concentration	M	ATP_c
[ADP] $_c$	Total cytosolic ADP concentration	M	ADP_c
[AMP] $_c$	Total cytosolic AMP concentration	M	AMP_c
[Pi] $_c$	Total cytosolic Pi concentration	M	Pi_c
[CrP] $_c$	Total cytosolic phosphocreatine concentration	M	AMP_c
[Cr] $_c$	Total cytosolic creatine concentration	M	Pi_c
[H $^+$ ] $_c$	cytosolic free proton concentration	M	h_c
[K $^+$ ] $_c$	cytosolic free potassium concentration	M	k_c
[Mg $^{2+}$ ] $_c$	cytosolic free magnesium concentration	M	m_c

2.2Mitochondrial Membrane Potential:¶

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

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

(1)

where $J_{C1}$ , $J_{C3}$ , and $J_{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 F $_1$ F $_0$ ATPase turnover rate is $J_{F1F0}$ and $nH_{F1F0}$ (= 8/3) is the proton flux stoichiometric number associated with the synthesis of one ATP. The fluxes $J_{ANT}$ and $J_{Hleak}$ are the adenine nucleotide translocator and proton leak fluxes.

2.3Mitochondrial Matrix Metabolite State Variables:¶

Metabolite concentrations in the matrix are governed by:

\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*}

(2)

where $Vol_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 $\alpha_{C2}$ ( $=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:

\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*}

(3)

where $Vol_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:

\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*}

(4)

where $Vol_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 $V_{Rm}/V_{Rc}$ is ratio of regional volume of the IMS to the cytosolic space. Since the $J_{ATPPERM}$ , $J_{ADPPERM}$ , $J_{AMPPERM}$ , and $J_{PIPERM}$ fluxes are in units of mass per unit time per unit mitochondrial volume, the multiplication by $V_{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 Mg $^{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., ATP $^{4-}$ ) and cation bound species (e.g., HATP $^{3-}$ , KATP $^{3-}$ , and MgATP $^{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

Flux	Definition	Units	Variable name in code
$J_{C1}$	Electron transport chain Complex I flux	mole $\cdot$ sec $^{-1}\cdot$ (L mito) $^{-1}$	J_ETC1_im_to_matrix
$J_{C3}$	Electron transport chain Complex III flux	mole $\cdot$ sec $^{-1}\cdot$ (L mito) $^{-1}$	J_ETC3_im_to_matrix
$J_{C4}$	Electron transport chain Complex IV flux	mole $\cdot$ sec $^{-1}\cdot$ (L mito) $^{-1}$	J_ETC4_im_to_matrix
$J_{F1F0}$	Mitochondrial F $_1$ F $_0$ ATPase flux	mole $\cdot$ sec $^{-1}\cdot$ (L mito) $^{-1}$	J_F1F0ATPASE_im_to_matrix
$J_{ANT}$	Adenine nucleotide translocase flux	mole $\cdot$ sec $^{-1}\cdot$ (L mito) $^{-1}$	J_ANT_im_to_matrix
$J_{Hleak}$	Proton leak flux	mole $\cdot$ sec $^{-1}\cdot$ (L mito) $^{-1}$	J_HLEAK_im_to_matrix
$J_{DH}$	Rate of NADH production	mole $\cdot$ sec $^{-1}\cdot$ (L mito) $^{-1}$	J_DH_matrix
$J_{PIC}$	Mitochondrial phosphate carrier flux	mole $\cdot$ sec $^{-1}\cdot$ (L mito) $^{-1}$	J_PIC_im_to_matrix
$J_{ATPPERM}$	Mitochondrial outer membrane ATP permeation	mole $\cdot$ sec $^{-1}\cdot$ (L mito) $^{-1}$	J_ATPPERM_cytoplasm_to_im
$J_{ADPPERM}$	Mitochondrial outer membrane ADP permeation	mole $\cdot$ sec $^{-1}\cdot$ (L mito) $^{-1}$	J_ADPPERM_cytoplasm_to_im
$J_{AMPPERM}$	Mitochondrial outer membrane AMP permeation	mole $\cdot$ sec $^{-1}\cdot$ (L mito) $^{-1}$	J_AMPPERM_cytoplasm_to_im
$J_{PIPERM}$	Mitochondrial outer membrane Pi permeation	mole $\cdot$ sec $^{-1}\cdot$ (L mito) $^{-1}$	J_PIPERM_cytoplasm_to_im
$J_{CK}$	Cytosolic creatine kinase flux	mole $\cdot$ sec $^{-1}\cdot$ (L cytosol) $^{-1}$	J_CK_cytoplasm
$J_{AK}$	Cytosolic adenylate kinase flux	mole $\cdot$ sec $^{-1}\cdot$ (L cytosol) $^{-1}$	J_AK_cytoplasm
$J_{ATPase}$	Cytosolic ATP hydrolysis flux	mole $\cdot$ sec $^{-1}\cdot$ (L cytosol) $^{-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 variable	Definition	Units used in code	Values
TAN	total adenine nucleotide pool	mole $\cdot$ (l cell) $^{-1}$	0.0071–0.0086 for sham 0.0052–0.0085 for TAC
CRtot	total creatine pool	mole $\cdot$ (l cell) $^{-1}$	0.0267–0.0330 for sham 0.0146–0.0278 for TAC
TEP	total exchangeable phosphate pool	mole $\cdot$ (l cell) $^{-1}$	0.0247–0.0298 for sham 0.0181–0.0293 for TAC
Ox_capacity	oxidative capacity (relative to control)	unitless	0.834–1.1526 for sham 0.5287–0.9755 for TAC
x_ATPase	ATP hydrolysis rate	mmole $\cdot$ sec $^{-1}\cdot$ (l cell) $^{-1}$	0.9301–1.392 for sham 0.6758–1.410 for TAC

Table 4:Output arguments for cellular energetics model

Output variable	Definition	Units used in code	Values
MgATP_cytoplasm	cytosolic [MgATP]	mmole $\cdot$ (l cytosol water) $^{-1}$	7.2145–9.4386 for sham 4.9384–9.0634 for TAC
MgADP_cytoplasm	cytosolic [MgADP]	mmole $\cdot$ (l cytosol water) $^{-1}$	0.0467–0.0598 for sham 0.0274–0.0549 for TAC
fPi_cytoplasm	cytosolic unchelated [Pi]	mmole $\cdot$ (l cytosol water) $^{-1}$	0.4393–1.3073 for sham 1.1437–1.6650 for TAC
MVO2_tissue	oxygen consumption rate	$\mu$ mole $\cdot$ min $^{-1}\cdot$ (g tissue) $^{-1}$	7.617–12.380 for sham 6.9104–12.450 for TAC
dGrATPase	ATP hydrolysis free energy	kJ $\cdot$ mole $^{-1}$	-(66.054–62.71) for sham -(63.574–61.44) for TAC
PCrATP	CrP/ATP ratio	unitless	1.817–2.324 for sham 1.4985–1.897 for TAC
ATP_cyto	cytosolic total [ATP]	mmole $\cdot$ (l cytosol water) $^{-1}$	8.733–11.438 for sham 5.983–10.959 for TAC
ADP_cyto	cytosolic total [ADP]	mmole $\cdot$ (l cytosol water) $^{-1}$	0.1164–0.1486 for sham 0.0682–0.1359 for TAC
Pi_cyto	cytosolic total [Pi]	mmole $\cdot$ (l cytosol water) $^{-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 $N$ , the permissible (calcium-bound) state $P$ , loosely attached state $A_1$ , strongly attached state $A_2$ , and post-ratcheted state $A_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

State
Variable

Definition

Units used in code

Variable name
in code

p_1^0

The 0 $^{\rm th}$ moment of state A1 strain
probability distribution. Equal to the
proportion of cross-bridges in state A1.

unitless

P1_0

p_1^1

The 1 $^{\rm st}$ moment of state A1 strain
probability distribution.

\mu

P1_1

p_1^2

The 2 $^{\rm nd}$ moment of state A1 strain
probability distribution.

\mu

^2

P1_2

p_2^0

The 0 $^{\rm th}$ moment of state A2 strain
probability distribution. Equal to the
proportion of cross-bridges in state A2.

unitless

P2_0

p_2^1

The 1 $^{\rm st}$ moment of state A2 strain
probability distribution.

\mu

P2_1

p_2^2

The 2 $^{\rm nd}$ moment of state A2 strain
probability distribution.

\mu

^2

P2_2

p_3^0

The 0 $^{\rm th}$ moment of state A3 strain
probability distribution. Equal to the
proportion of cross-bridges in state A3.

unitless

P3_0

p_3^1

The 1 $^{\rm st}$ moment of state A3 strain
probability distribution.

\mu

P3_1

p_3^2

The 2 $^{\rm nd}$ moment of state A3 strain
probability distribution.

\mu

^2

P3_2

N

Non-permissible XB state

unitless

U_{NR}

Non relaxed state

unitless

U_NR

The equations used to simulate the cross-bridge model are

\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*}

(5)

where v is the velocity of sliding ( $v=dSL/dt$ where $SL$ 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:

\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*}

(6)

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 ( $\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:

\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*}

(7)

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

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

(8)

where $P$ 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.

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

(9)

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

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

(10)

where $\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:

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

(11)

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

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

(12)

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

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

(13)

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

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

(14)

Here $SL$ is the length of sarcomere, $L_{thick}$ , $L_{thin}$ , $L_{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

Parameter	Definition	Value and units	Parameter name in code
$L_thin$	thin filament length	1200 nm	L_thin
$L_thick$	thick filament length	1670 nm	L_thick
$L_bare$	bare length of the thick filament	100 nm	L_hbare
$OV_{Z-axis}$	overlap region closest to the Z-axis	nm	sovr_ze
$OV_{M-line}$	overlap region closest to the M-line	nm	sovr_cle
$LOV$	length of overlap	nm	L_sovr
$OV_{thick}$	fraction of thick filament overlap	unitless	N_overlap

3.5Active and passive force:¶

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

\sigma_{XB} (t) = OV_{thick} ( k_{stiff,1} \left(p_2^2 + p_3^2 \right) + k_{stiff,2} \Delta r p_3^0 )

(15)

where $k_{stiff,1}$ and $k_{stiff,2)}$ are stiffness constants, $\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, $F_1$ and $F_2$ , and a series element force $F_{SE}$ . Overall force balance for the model yields

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

(16)

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

\sigma_1 = \eta \frac{dSL}{dt} \, .

(17)

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

\sigma_2 (SL)= k_{passive} \left(SL - SL_{rest} \right) + \sigma_{Passive_collagen} (SL)

(18)

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

\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*}

(19)

where $Pcon_{collagen} = 0.01$ (unitless) is the scale factor for passive force contributed by collagen. $PExp_{collagen} = 70$ $\mu$ m $^{-1}$ is a model parameter and $SL_{collagen} = 2.25$ $\mu$ m 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 $\mu$ m. Isometric tension time courses were measured at different stimulation frequencies and sarcomere lengths. Figure 3 shows data on peak developed tension ( $T_{dev}$ ) as a function of $SL$ at stimulation frequency of 4 Hz; and data on relaxation time from peak to 50% of peak tension ( $RT50$ ); peak developed tension ( $T_{dev}$ ) and time to peak tension ( $TTP$ ) as function of $SL$ .

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.9$   $\mu$ m the $Ca_{50} = 5.89$ and Hill coefficient $n = 4.63$ and for $SL = 2.3$   $\mu$ m 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 Janssen *et al.*, 2002.

Table 7:Model parameters for cross bridge model

Parameter	Definition	Value and units	Parameter name in code	Reference
$k_{stiff,1}$	Stiffness constant due to myosin–actin interaction	13013 mmHg/ $\mu$ m	kstiff1	Fit to data in Figure 3
$k_{stiff,2}$	Stiffness constant due to working stroke of XBs	341590 mmHg/ $\mu$ m	kstiff2	Fit to data in Figure 3
$k_{passive}$	Passive stiffness constant	25 mmHg/ $\mu$ m	k_passive	Fit to data in Figure 3
$SL_{rest}$	Sarcomere length at which passive force is zero	1.8 $\mu$ m	SL_rest_pas	Campbell et al. (2018)
$\alpha_1$	Strain-dependency parameter	10 $\mu$ m $^{-1}$	alpha1	Tewari et al. (2016)
$\alpha_2$	Strain-dependency parameter	9 $\mu$ m $^{-1}$	alpha2	Tewari et al. (2016)
$\alpha_3$	Strain-dependency parameter	5.93 $\mu$ m $^{-2}$	alpha3	Fit to data in Figure 3
$s_3$	Strain-dependency parameter	$9.9 \times 10^{-3}$ $\mu$ m	s3	Tewari et al. (2016)
$k_{coop}$	Strength of thin filament cooperativity	1.857	K_coop	Fit to data in Figure 3
$k_{on}$	Rate constant of Ca binding to troponin C	101.2  $\mu$ M $^{-1}\cdot$ s $^{-1}$	k_on	Fit to data in Figure 3
$k_{off}$	Rate constant of Ca unbinding from troponin C	101.2 s $^{-1}$	k_off	Fit to data in Figure 3
$k_{force}$	Force-dependent rate constant of super relax transition	$1.169 \times 10^{-3}$ N $^{-1}\cdot$ m $^{-2}$	kforce	Fit to data in Figure 3
$k_{SR}$	Rate constant of force-dependent super relax transition	14.44 s $^{-1}$	ksr	Fit to data in Figure 3
$k_{-SR}$	Rate constant of force-dependent super relax transition	50.03 s $^{-1}$	kmsr	Fit to data in Figure 3
$k_{MgATP}$	[MgATP] dissociation constant	489.7 $\mu$ M	K_T	Tewari et al. (2016)
$k_{MgADP}$	[MgADP] dissociation constant	194.0 $\mu$ M	K_D	Tewari et al. (2016)
$k_{Pi}$	[Pi] dissociation constant	4.0 mM	K_Pi	Tewari et al. (2016)
$k_a$	Myosin–actin rate of attachment	559.6 s $^{-1}$	ka	Fit to data in Figure 3
$k_d$	Myosin–actin rate of un-attachment	304.7 s $^{-1}$	kd	Fit to data in Figure 3
$k_1$	Transition rate constant	112.4 s $^{-1}$	k1	Fit to data in Figure 3
$k_{-1}$	Rate of strongly-bound to weakly-bound transition	21.30 s $^{-1}$	km1	Tewari et al. (2016)
$k_2$	Rate of ratcheting	811.7 s $^{-1}$	k2	Tewari et al. (2016)
$k_{-2}$	Reverse ratcheting rate	43.25 s $^{-1}$	km2	Tewari et al. (2016)
$k_3$	Myosin–actin detachment rate	144.6 s $^{-1}$	k3	Fit 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

\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*}

(20)

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

The $SL$ and $\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)$ 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:

\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*}

(21)

where

\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*}

(22)

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

Table 8:State variables in TriSeg (Heart) model

State Variable

Definition

Units

Variable name in code

_m,RV

Maximal axial distance
from RV midwall surface to origin

Xm_RV

_m,LV

Maximal axial distance
from LV midwall surface to origin

Xm_LV

_m,SEP

Maximal axial distance
from SEP midwall surface to origin

Xm_SEP

_m

Radius of midwall junction circle

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:

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

(23)

Axial and radial components of the tension are computed

\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*}

(24)

The four unknowns of the model— $x_{m,RV}$ , $x_{m,LV}$ , $x_{m,SEP}$ , and $y_{m}$ —are determined by satisfying the relations

\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*}

(25)

Transmural pressures are computed

P_{trans,\#} = \frac{2T_{x,\#}}{y_{m}}

(26)

and the pressures in the cavities are computed

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

(27)

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

Parameter

Definition

Value and Units

Variable name in code

V_{w,LV}

LV wall volume

input
(experimental data), mL

Vw_LV

V_{w,SEP}

Septal wall volume

input
(experimental data), mL

Vw_SEP

V_{w,RV}

RV wall volume

input
(experimental data), mL

Vw_RV

A_{m,ref}^{LW}

LV midwall reference
surface area

adjustable, cm

^{2}

Amref_LV

A_{m,ref}^{SEP}

Septal midwall reference
surface area

adjustable, cm

^{2}

Amref_SEP

A_{m,ref}^{RW}

RV midwall reference
surface area

adjustable, cm

^{2}

Amref_RV

K_{SE}

Stiffness of series element

50000 mmHg/

\mu

Kse

\eta

Viscosity coefficient of myofibers

1 mmHg

\cdot

sec

\cdot \mu

^{-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 $R_{pul}$ , pulmonary arterial and venous compliances $C_{PA}$ and $C_{PV}$ , systemic arterials resistances $R_{Ao}$ and $R_{sys}$ , systemic arterial compliance $C_{SA}$ , and systemic venous compliance $C_{SV}$ , systemic arterial resistance and aortic compliance $C_{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). $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 Lumens *et al.* (2009), described below.

Flow through a resistive element is calculated

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

(28)

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

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

(29)

where $F_{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

Parameter	Definition	Units used in code	Variable name in code
Cardiovascular system model state variables
$V_{LV}$	Left ventricle volume	mL	V_LV
$V_{RV}$	Right ventricle volume	mL	V_RV
$V_{SV}$	Volume of systemic vein	mL	V_SV
$V_{PV}$	Volume of pulmonary vein	mL	V_PV
$V_{SA}$	Volume of systemic artery	mL	V_SA
$V_{PA}$	Volume of pulmonary artery	mL	V_PA
$V_{Ao}$	Volume of aorta	mL	V_Ao
Pressures computed from volume state variables
$P_{SV}$	Systemic venous pressure	mmHg	P_SV
$P_{PV}$	Pulmonary venous pressure	mmHg	P_PV
$P_{PA}$	Pulmonary arterial pressure	mmHg	P_PA
$P_{Ao}$	Aortic pressure (proximal to TAC)	mmHg	P_Ao
$P_{SA}$	Systemic arterial pressure (distal to TAC)	mmHg	P_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

Parameter

Definition

Value and Units

Variable name
in code

C_{Ao}

Proximal aortic compliance

0.0022045 mL

\cdot

mmHg

^{-1}

C_Ao

C_{SA}

Systemic arterial compliance

0.0077157 mL

\cdot

mmHg

^{-1}

C_SA

C_{SV}

Systemic venous compliance

2.5 mL

\cdot

mmHg

^{-1}

C_SV

C_{PA}

Pulmonary arterial compliance

0.013778 mL

\cdot

mmHg

^{-1}

C_PA

C_{PV}

Pulmonary venous compliance

0.25 mL

\cdot

mmHg

^{-1}

C_PV

R_{Ao}

Proximal aortic resistance

2.5 mmHg

\cdot

sec

\cdot

^{-1}

R_Ao

R_{sys}

Systemic vasculature resistance

adjustable, mmHg

\cdot

sec

\cdot

^{-1}

R_sys

R_{SV}

Systemic veins resistance

0.25 mmHg

\cdot

sec

\cdot

^{-1}

R_SV

R_{TAo}

Transmural aortic resistance

0.5 mmHg

\cdot

sec

\cdot

^{-1}

R_tAo

R_{TSA}

Transmural systemic artery resistance

4 mmHg

\cdot

sec

\cdot

^{-1}

R_tSA

R_{PA}

Pulmonary vasculature resistance

7.58 mmHg

\cdot

sec

\cdot

^{-1}

R_PA

R_{PV}

Pulmonary veins resistance

0.25 mmHg

\cdot

sec

\cdot

^{-1}

R_PV

R_{VLV}

Valve resistance

0.05 mmHg

\cdot

sec

\cdot

^{-1}

R_VLV

R_{TAC}

Resistance of TAC

adjustable, mmHg

\cdot

sec

\cdot

^{-1}

R_TAC

The systemic compliances ( $C_{SA}$ , $C_{Ao}$ , $C_{SV}$ ) are fixed to produce a pulse pressure of roughly 33 mmHg for simulations of sham control rats. The pulmonary compliance ( $C_{PA}$ , $C_{PV}$ ) are fixed to roughly have the target value of 12 mmHg for the pulmonary pulse pressure. The resistance $R_{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 $R_{SV}$ is set so that the mean pressure in the systemic veins for sham control rats is 3 mmHg. The pulmonary resistances $R_{PA}$ and $R_{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 $R_{TAC}$ and $R_{sys}$ . The resistance $R_{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 $\Delta P_{TAC} = \frac{1}{2} \rho (V_{2}^{2} - V_{1}^{2})$ , where $V_{1}$ and $V_{2}$ are the velocities on either side of the TAC constriction. Given an estimated pressure drop the resistance is computed $R_{TAC} = \frac{\Delta P_{TAC}}{CO}$ , where CO is the cardiac output. The systemic resistance $R_{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¶

6.1Summary:¶

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

Parameter	Variable name in code	Definition	Units used in code	Values
HR	HR	Heart rate	bpm	318.6–367 for sham 294–349 for TAC. Set to measured value
Vw,LV	Vw_LV	LV wall volume	mL	0.564–0.853 for sham 0.863–1.22 for TAC. Set as 2/3 of measured LV volume
Vw,SEP	Vw_SEP	Septal wall volume	mL	0.2821–0.4266 for sham 0.431–0.61 for TAC. Set as 1/3 of measured LV volume
Vw,RV	Vw_RV	RV wall volume	mL	0.282–0.373 for sham 0.249–0.605 for TAC. Set as measured RV free wall volume
$A_{m,ref}^{LW}$	Amref_LV	LV midwall reference surface area	cm $^{2}$	1.6865–2.465 for SHAM 2.299–2.966 for TAC. Adjusted to fit data
$A_{m,ref}^{SEP}$	Amref_SEP	Septal midwall reference surface area	cm $^{2}$	1.055–1.487 for sham 1.299–1.620 for TAC. Adjusted to fit data
$A_{m,ref}^{RW}$	Amref_RV	RV midwall reference surface area	cm $^{2}$	2.629–3.612 for sham 3.135–5.538 for TAC. Adjusted to fit data
$k_{force}$	k_force	Model parameter for force dependent super relax transition	N $^{-1}$ m $^{2}$	(0.905–3.312) $\times 10^{-3}$ for sham (1.464–3.783) $\times 10^{-3}$ for TAC. Adjusted to fit data
$k_{SR}$	ksr	On rate constant for super relax state	s $^{-1}$	9.021–33.013 for sham 14.598–37.708 for TAC. Adjusted to fit data
$R_{sys}$	R_sys	Systemic vasculature resistance	mmHg $\cdot$ sec $\cdot$ mL $^{-1}$	38.071–71.886 for sham 50.577–125.33 for TAC. Adjusted to fit data
$R_{TAC}$	R_TAC	resistance of TAC	mmHg $\cdot$ sec $\cdot$ mL $^{-1}$	0 for sham 7.11–22.25 for TAC. Adjusted to fit data
MgATP_cyto-plasm	MgATP	cytosolic [MgATP]	mmole $\cdot$ (L cytosol water) $^{-1}$	7.214–9.439 for sham 4.938–9.063 for TAC . Predicted from energetics model
MgADP_cyto-plasm	MgADP	cytosolic [MgADP]	mmole $\cdot$ (L cytosol water) $^{-1}$	0.047–0.06 for sham 0.027–0.055 for TAC. Predicted from energetics model
Pi_cyto-plasm	Pi	cytosolic total [Pi]	mmole $\cdot$ (L cytosol water) $^{-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.31$ mL $\cdot$ (min $\cdot$ g) $^{-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 $\rho_{cell} = 0.694$ mL $\cdot$ g $^{-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.30$ mmol $\cdot$ (sec $\cdot$ L cell) $^{-1}$ . (This value is approximately 2.5 times higher than the estimated resting ATP hydrolysis rate of to $0.547$ mmol $\cdot$ (sec $\cdot$ L cell) $^{-1}$ for human myocardium Gao et al., 2019.)

Simulations of the cardiovascular mechanics model predict an average cross-bridge cycling rate of 5.03 sec $^{-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

\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*}

(30)

If we assume that roughly 3/4 of the ATP consumed by the cardiomyocyte is consumed by the myosin ATPase, we estimate from this constant of proportionality the density of cross-bridge-forming units in a cardiomyocyte to be roughly 0.34 mmol per liter of cell. This density has been independently estimated to be 0.25 mmol per liter of cell for skeletal myocytes Barclay et al., 2010.

6.3Fitting data on individual rats:¶

The full set of adjustable parameters invoked in the cardiovascular systems model is listed in Table 12. Certain adjustable parameters are set based on direct measurements and others are adjusted so that simulation outputs match measured data. Seven parameter values— $A_{m,ref}^{LW"}$ , $A_{m,ref} ^{SEP}$ , $A_{m,ref}^{RW}$ , $k_{force}$ , $k_{SR}$ , $R_{sys}$ , $R_{TAC}$ —are adjusted to fit model predictions to data from individual animals on end-systolic and end-diastolic volumes and estimated pressure drop across the aortic constriction in TAC animals, and to simultaneously maintain a fixed mean arterial pressure of 93.3 mmHg, to maintain end-diastolic sarcomere lengths in the LV, septum, and RV of 2.2 $\mu$ m.

Ranges of estimated values of are listed in Table 12. Values of anatomical/geometric parameters representing heart masses and reference areas are higher in TAC rats than in sham control rats, reflecting hypertrophic remodeling. The parameters $k_{force}$ and $k_{SR}$ govern the transition out of the inaccessible super-relaxed state in the calcium-activation model. Increased in the values of these parameters represent increased levels of phosphorylation of myosin binding protein C. Thus, the higher values of these parameters in TAC compared to sham animals represent a prediction that phosphorylation of this protein is increased in the TAC animals.

Simulations of TAC rats consistently show lower ATP, ADP, and CrP, reflecting reductions in the total adenine nucleotide and creatine pools. (See below.) Lower ADP levels require a compensatory increase in inorganic phosphate to maintain ATP synthesis.

The ranges of values of cross-bridge cycle rate, ATP hydrolysis rate, and end-systolic and end-diastolic volumes are listed in Table 13.

Table 13:Output arguments for cardiovascular mechanic model

Variable	Definition	Units used in code	Values
rate_of_XB_turn-over_ave	Cross bridge cycling rate	sec $^{-1}$	3.5874–5.3658 for sham 2.606–5.4381 for TAC
x_ATPase	ATP hydrolysis rate	mmol $\cdot$ (sec $\cdot$ L $\cdot$ cell) $^{-1}$	0.9301–1.392 for sham 0.6758–1.410 for TAC
EDLV	End diastolic left ventricular volume	mL	0.303–0.547 for sham 0.407–0.650 for TAC
ESLV	End systolic left ventricular volume	mL	0.093 – 0.232 for sham 0.138 – 0.423 for TAC
MAP	Mean arterial pressure	mmHg	93.3 for both TAC and sham

The set of input parameters invoked in the cardiovascular systems model is listed in Table 3. The input parameters TAN, CRtot, and Ox_capacity are all measured for each individual rat. The ATP hydrolysis rate ( $x_{ATPase}$ ) is predicted by the mechanics model (see above). The relationships between metabolite pools from Gao et al. (2019) are used to estimate the total exchangeable phosphate (TEP):

{\rm TEP} = P_{0} - p ((A_{0} - {\rm TAN} )/ a)

(31)

where, $P_{0}= 35.446$ mmol $\cdot$ (L cell) $^{-1}$ , $A_{0} = 10.26$ mmol $\cdot$ (L cell) $^{-1}$ , are reference TEP and TAN values, $p = 0.283$ mmol $\cdot$ (L cell $\cdot$ year) $^{-1}$ and $a = 0.082$ mmol $\cdot$ (L cell $\cdot$ year) $^{-1}$ are the slopes of the relationships between mean TEP and mean TAN and age in humans, respectively.

Output arguments from the energetics model, including the metabolite levels used in the mechanics model, are listed in Table 4.

6.4Predictions associated with changing metabolic/energetic parameterization:¶

The predictions illustrated in Fig. 9 of Lopez et al. (2020) are made by replacing the parameters representing the metabolic state of sham control rats with those representing TAC rats and by replacing the parameters representing the metabolic state of TAC rats with those representing sham control rats.

For predictions of how mechanical function in sham rats changes when the metabolic profile is replaced with that of the average TAC rat (Fig. 9A - 9B Lopez et al. (2020)), the input metabolic parameters are set to:

\begin{align*} {\rm TAN} &= 0.006976 \, {\rm M} \nonumber \\ {\rm CRtot} &= 0.02303 \, {\rm M} \nonumber \\ {\rm TEP} &= 0.02411 \, {\rm M} \nonumber \\ {\rm Ox\_capacity} &= 0.7482 \, {\rm (unitless)} \nonumber \end{align*}

(32)

Note that the above average TAC rat metabolite profile are based on $n = 10$ TAC rats. (TAC #7 is excluded for reasons described in Lopez et al. (2020).)

Given these values specifying the metabolic model, the blood volume, $k_{force}$ , and $k_{SR}$ are increased so that the model-predicted mean arterial pressure was 93.3 mmHg. In other words, it is assumed that baseline cardiac output and system pressure are maintained at the original physiological levels via compensatory increases in preload and myosin binding protein C phosphorylation. To obtain the model predictions in Lopez et al. first the $k_{force}$ and $k_{SR}$ are proportionately increased to compensate for the metabolic dysfunction. If it is not possible to restore mean arterial pressure to the physiological level with that change alone, blood volume is increased until the mean pressure of 93.3 mmHg is reached.

For predictions of how mechanical function in TAC rats changes when the metabolic profile is replaced with that of the average sham rat (Fig. 9C - 9D), the input metabolic parameters are set to:

\begin{align*} {\rm TAN} &= 0.007624 \, {\rm M} \nonumber \\ {\rm CRtot} &= 0.03027 \, {\rm M} \nonumber \\ {\rm TEP} &= 0.02635 \, {\rm M} \nonumber \\ {\rm Ox\_capacity} &= 1 \, {\rm (unitless)} \nonumber \end{align*}

(33)

Given these values specifying the metabolic model, the blood volume, $k_{force}$ , and $k_{SR}$ , and $R_{sys}$ are adjusted so that the model-predicted mean arterial pressure is 93.3 mmHg and the predicted end-diastolic volume with mean sham metabolic parameters is equal to that of the original TAC rat. In other words, it was assumed that baseline system pressure and diastolic filling level are maintained at the original physiological levels via compensatory reduction in preload, myosin binding protein C phosphorylation, and systemic resistance.

7Running the model¶

7.1Summary of codes¶

The simulation package consists of 5 MATLAB files:

CardiovascualarMechanics.m
dXdT_CardiovascularMechanics.m
EnergeticsModelScript.m
dXdT_energetics.m
TrisegEquations.m

Values of all adjustable parameters are stored in spreadsheet files “Adjustable_paramaters_table_rest.xlsx” for the baseline simulations and “Adjustable_parameters_table_SWAP.xlsx” for simulations with replaced metabolic parameters.

The CardiovascularMechanics.m function is the main driver to run the mechanics model for a given animal. For instance, assigning the variable “rat_number” to 1 will execute simulations for SHAM rat number 1. Executing the script will load the parameters associated with this animal, run the cardiovascular systems model for 120 heart beats to attain a periodic steady state, and then plot the predicted left ventricular pressure, aortic pressure, and arterial pressure along with the left-ventricular pressure volume loop for this individual animal. The target end-diastolic and end-systolic volumes and the pres-sure drop across the TAC will be indicated by dashed lines in the figures. The model will compute the predicted cross-bridge cycling rate in the LV free wall and the associated ATP hydrolysis rate.

The energetics model for a given animal is called within the cardiovascular model. Executing this script will read the metabolic/energetic parameters associated with this animal and run the model to calculate the cytosolic metabolite levels. Note the ATPase hydrolysis rate for steady state is pre-identified and is listed in column 9 of the input adjustable_parameters_table_rest.xlsx.

The numbering for rats in the input “data1.xlsx” file is as follows: the first 8 rats are SHAM rats and rat number 9 is the mean sham rat; rat number 10 to 19 are the TAC rats (TAC rat# 7 is excluded). Therefore the first TAC rat is rat number 10 in the simulations. For example, to simulate the model for TAC rat #1 we need assign number 10 to the variable “rat_numbers = 10” in the CardiovascularMehanics.m and run the code.

The model will generate output shown below in Figure 5.

In addition to the plots illustrated above, the running the model for this individual animal results in a predicted ATP hydrolysis rate of 0.675 mmole $\cdot$ sec $^{-1}\cdot$ (L cell) $^{-1}$ . Note that the input ATP hydrolysis rate is the input parameters for the energetic model for this animal (column 9, rat number 10) of the adjustable variables in the input file “adjustable_parameters_table_rest.xlsx”) is equal to this value. Also note that the predicted values for MgATP_cytoplasm, MgADP_cytoplasm, and Pi_cyto from the energetics model are equal to the associated input parameters for the mechanics model for this animal.

7.2Reproduction of simulations highlighted in original paper¶

The procedure detailed above may be followed to reproduce the model fits to data for any of the individual TAC or sham rats analyzed in the original paper. For example, Fig. 7 of the original paper shows model fits to sham rat #3 and TAC rat #2. To simulate the model for sham rat #3 the user assigns the number 3 to the variable “rat_numbers = 3” in the CardiovascularMehanics.m and runs the code. The resulting model output is equivalent to that presented in Fig. 7 panels A and B of the original paper.

To simulate the model for TAC rat #2 the user assigns the number 11 to the variable “rat_numbers = 11” in the CardiovascularMehanics.m and runs the code. The resulting model output is equivalent to that presented in Fig. 7 panels C and D of the original paper.

Fig. 9 of the original paper illustrates the predicted effects of replacing the metabolic profile of sham rat with that of at TAC rat, and the predicted effects of replacing the metabolic profile of a TAC rat with that of a sham rat. The solid lines in Fig. 9 of the original paper represent the resting-state simulations from Fig. 7 of the paper. The dashed lines correspond to simulations associated with the metabolite swaps.

To simulate the model for sham rat #3 with TAC metabolism settings the user assigns the number 3 to the variable “rat_numbers = 3” in the CardiovascularMehanics.m and also assigns the number 1 to the variable “flag_swap_metabolite = 1” and runs the code. The resulting model output is equivalent to that presented in Fig. 9 panels A and B of the original paper. To simulate the model for TAC rat #2 with sham metabolism settings the user assigns the number 11 to the variable “rat_numbers = 11” in the CardiovascularMehanics.m and also assigns the number 1 to the variable “flag_swap_metabolite = 1” and runs the code. The resulting model output is equivalent to that presented in Fig. 9 panels C and D of the original paper.

8Glossary of model codes¶

EnergeticsModelScript.m: This function is used to compute the cellular energetics concentration variables for the myocardium, given input values of mitochondrial oxidative capacity, TAN, TEP, and CRtot metabolite pool values, and the rate of cellular ATP hydrolysis. The function computes the steady state of the cellular energetics model by simulating the model governed by the differential equations in dXdT_energetics.m.
dXdT_energetics.m: This function is an implementation of the Bazil et al. (2016) model of rat myocardial mitochondrial oxidative phosphorylation. The model has 29 state variables, listed in Table 1.
Cardiovascularmechanics.m: This function simulates the pressure and flows in the whole-body cardio-vascular systems model of  1, governed by the five-compartment lumped-parameter cardiovascular systems model coupled to the Lumens et al. TriSeg heart model. The inputs to the model include the cytosolic metabolite concentrations predicted by EnergeticsModelScript.m. The outputs of the model are the myocardial ATP hydrolysis rate (used as an input for the energetics model) and the cardiovascular variables, EDLV, ESLV, MAP, rate of ATP cellular hydrolysis, to be compared to measurements for individual rats.
dXdT_cardiovascular_mechanics.m: This function is an implantation of the whole organ cardiovasucalar mechanics model. The model has 47 state variables listed in Tables 5, 8, and 10.
TriSeg.m: This function runs the TriSeg model equations to obtain estimates for initial value for ODE solver and the associated algebraic equations.

References¶

Bazil, J. N., Beard, D. A., & Vinnakota, K. C. (2016). Catalytic Coupling of Oxidative Phosphorylation, ATP Demand, and Reactive Oxygen Species Generation. Biophys J, 110(4), 962–971. 10.1016/j.bpj.2015.09.036
Tewari, S. G., M., B. S., Palmer, B. M., & Beard, D. A. (2016). Dynamics of cross-bridge cycling, ATP hydrolysis, force generation, and deformation in cardiac muscle. J Mol Cell Cardiol, 96, 11–25. 10.1016/j.yjmcc.2015.02.006
Tewari, S. G., M., B. S., Vinnakota, K. C., Rice, J. M., Janssen, P. M., & Beard, D. A. (2016). Influence of metabolic dysfunction on cardiac mechanics in decompensated hypertrophy and heart failure. J Mol Cell Cardiol, 94, 162–175. 10.1016/j.yjmcc.2016.04.003
Lumens, J., Delhaas, T., Kirn, B., & Arts, T. (2009). Three-wall segment (TriSeg) model describing mechanics and hemodynamics of ventricular interaction. Ann Biomed Eng, 37(11), 2234–2255. 10.1007/s10439-009-9774-2
Lopez, R., Marzban, B., Gao, X., Lauinger, E., Van den Bergh, F., Whitesall, S. E., Converso-Baran, K., Burant, C. F., Michele, D. E., & Beard, D. A. (2020). Impaired Myocardial Energetics Causes Mechanical Dysfunction in Decompensated Failing Hearts. Function. 10.1093/function/zqaa018