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

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.

## 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.

## 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:

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:

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:

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:

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

| Definition | Units used in code |
| ||||
---|---|---|---|---|---|---|---|

$p_1^0$ |
| unitless | P1_0 | ||||

$p_1^1$ |
| $\mu$m | P1_1 | ||||

$p_1^2$ |
| $\mu$m$^2$ | P1_2 | ||||

$p_2^0$ |
| unitless | P2_0 | ||||

$p_2^1$ |
| $\mu$m | P2_1 | ||||

$p_2^2$ |
| $\mu$m$^2$ | P2_2 | ||||

$p_3^0$ |
| unitless | P3_0 | ||||

$p_3^1$ |
| $\mu$m | P3_1 | ||||

$p_3^2$ |
| $\mu$m$^2$ | P3_2 | ||||

$N$ |
| unitless | N | ||||

$U_{NR}$ |
| unitless | U_NR | ||||

The equations used to simulate the cross-bridge model are

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:

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

### 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:

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

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.

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

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:

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

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

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

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:

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

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

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

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).

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.

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

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:

where

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 | ||
---|---|---|---|---|---|

x$_m,RV$ |
| cm | Xm_RV | ||

x$_m,LV$ |
| cm | Xm_LV | ||

x$_m,SEP$ |
| cm | Xm_SEP | ||

y$_m$ | Radius of midwall junction circle | cm | ym |

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:

Axial and radial components of the tension are computed

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

Transmural pressures are computed

and the pressures in the cavities are computed

### 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 |
| Vw_LV | ||

$V_{w,SEP}$ | Septal wall volume |
| Vw_SEP | ||

$V_{w,RV}$ | RV wall volume |
| Vw_RV | ||

$A_{m,ref}^{LW}$ |
| adjustable, cm$^{2}$ | Amref_LV | ||

$A_{m,ref}^{SEP}$ |
| adjustable, cm$^{2}$ | Amref_SEP | ||

$A_{m,ref}^{RW}$ |
| adjustable, cm$^{2}$ | Amref_RV | ||

$K_{SE}$ | Stiffness of series element | 50000 mmHg/$\mu$m | Kse | ||

$\eta$ | Viscosity coefficient of myofibers | 1 mmHg$\cdot$sec$\cdot \mu$m$^{-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}$.

Flow through a resistive element is calculated

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

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 |
| ||
---|---|---|---|---|---|

$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$ mL$^{-1}$ | R_Ao | ||

$R_{sys}$ | Systemic vasculature resistance | adjustable, mmHg$\cdot$sec$\cdot$ mL$^{-1}$ | R_sys | ||

$R_{SV}$ | Systemic veins resistance | 0.25 mmHg$\cdot$sec$\cdot$ mL$^{-1}$ | R_SV | ||

$R_{TAo}$ | Transmural aortic resistance | 0.5 mmHg$\cdot$sec$\cdot$ mL$^{-1}$ | R_tAo | ||

$R_{TSA}$ | Transmural systemic artery resistance | 4 mmHg$\cdot$sec$\cdot$ mL$^{-1}$ | R_tSA | ||

$R_{PA}$ | Pulmonary vasculature resistance | 7.58 mmHg$\cdot$sec$\cdot$ mL$^{-1}$ | R_PA | ||

$R_{PV}$ | Pulmonary veins resistance | 0.25 mmHg$\cdot$sec$\cdot$ mL$^{-1}$ | R_PV | ||

$R_{VLV}$ | Valve resistance | 0.05 mmHg$\cdot$sec$\cdot$ mL$^{-1}$ | R_VLV | ||

$R_{TAC}$ | Resistance of TAC | adjustable, mmHg$\cdot$sec$\cdot$ mL$^{-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