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 diﬀerences 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


Introduction
The multi-scale modeling approach is summarized in Fig. 1.Model components representing myocardial metabolism Bazil et al. (2016), myocardial cell mechanics Tewari et al. (2016a,b), myocardial whole-organ pumping Lumens et al. (2009), and a simple lumped circulatory model are integrated together to simulate whole-body cardiovascular function.
• 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. (2016a,b) 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.

Model 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 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 below.

Mitochondrial Membrane Potential:
The potential difference across the mitochondrial inner membrane is governed by currents across the membrane: where J C 1 , J C 3 , and J C 4 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 F 1F 0 and nH F 1F 0 (= 8/3) is the proton flux stoichiometric number associated with the synthesis of one ATP.The fluxes J ANT and J H l e ak are the adenine nucleotide translocator and proton leak fluxes. 3/21

Mitochondrial Matrix Metabolite State Variables:
Metabolite concentrations in the matrix are governed by: where V ol 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 α C 2 (= 1/4) accounts for ubiquinone reduction via complex II.

Inter-Membrane Space (IMS) Metabolite State Variables:
Metabolite concentrations in the intermembrane space are governed by: where V ol 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.

Cytosolic Metabolite State Variables:
Metabolite concentrations in the cytosolic space are governed by: where V ol 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 R m /V R c is ratio of regional volume of the IMS to the cytosolic space.Since the J AT P P E R M , J AD P P E R M , J AM P P E R M , and J P I P E R M fluxes are in units of mass per unit time per unit mitochondrial volume, the multiplication by V R m /V R c 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.
The full set of equations is detailed in the supplementary material published with Bazil et al. (2016).

Energetic Model Fluxes:
The fluxes on the right-hand sides of Equations (1,2,3, and 4) are defined in Table 2. J_AMPPERM_cytoplasm_to_im The mathematical expressions for these fluxes are detailed in Bazil et al. (2016).

Implementation in Multiscale Model:
The cellular energetics model is implemented in a MATLAB script called EnergeticsModelScript.m.This script is used to predict cytosolic 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.

Model Variables and Equations:
A cardiomyocyte mechanics model based on the models of Tewari et al. (2016a,b) and Campbell et al. (2018) is used to simulate the active and passive components of myocardial wall tension 5/21   The equations used to simulate the cross-bridge model are where v is the velocity of sliding (v = d S L/d t where S L 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 7/21 relations: A detailed description of the moment expansion and associated equations is given in Tewari et al. (2016a).

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

Super-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 S R ) to non-relaxed (U N R ) state is force-dependent: where σ X B , is active contractile force and formulated in equation 12.

Overlap function:
Following Rice et al. (2008) the fractional overlap between thin and thick filament is represented as follows: where OV Z −axi s is the overlap region closest to the Z-axis, where OV M −l i ne 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 S L is the length of sarcomere, L t hi ck , L t hi n , L bar e are the length of thick filament, bare region of the thick filament and, the length of the thin filament.

Active and passive force:
The active force generated by cross-bridges is computed from contributions from pre-and postratcheted states: where k st i f f ,1 and k st i f f ,2) are stiffness constants, ∆r is the cross bridge strain associated with ratcheting deformation.The full muscle model (Fig. 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 S E .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 σ 2 is a function of sarcomere length and is calculated σ 2 (S L) = k p assi v e (S L − S L r est ) + σ P assi v e c ol l ag en (S L) (15)

9/21
where k p assi v e is the stiffness parameters for the passive force, S L r est is the sarcomere rest length, and the passive force exerted by collagen is adopted from Rice et al. (2008).
σ P assi v e c ol l ag en (S L) = P con col l ag en e P E x p col l ag en (S L−S L col l ag en ) − 1 i f S L > S L col l ag en 0 otherwise .
where P con col l ag en = 0.01 (unitless) is the scale factor for passive force contributed by collagen.P E x p col l ag en = 70 µm −1 is a model parameter and S L col l ag en = 2.25 µm is the minimum length threshold at which collagen starts to exert force on sarcomere.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.

Model Variables and Equations:
A modified version of the Lumens et al. (2009) TriSeg model is used to simulate left-and rightventricular mechanics, based on the implementation of Tewari et al. (2016b).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. ( 14) the rates of change of sarcomere length in these three segments is given by 10/21 where the parameter λ X B is a scalar used to account for differences in force generation in vivo versus in vitro.(The value value λ X B = 1.4,determined by Tewari et al. (2016b) accounts for slightly lower force generation in vitro versus in vivo.) The S L and d S L d t computed from Eq. ( 17) 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. ( 17), which governs S L (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,r ef ,# 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,S E P , and y m that determine the geometry of the ventricular cavities, are listed in Table 8.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,S E P , and y m -are determined by satisfying the relations ,S E P T x,RV + T x,LV + T x,S E P = 0 T y ,RV + T y ,LV + T y ,S E P = 0 .
Transmural pressures are computed and the pressures in the cavities are computed P LV = −P t r ans,LV P RV = +P t r ans,RV .

Model Parameters
Parameters defining the mass and geometry of the heart are identified from data on individual animals are defined in Table 9.  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 i n − F out is the rate of change of blood volume in the element.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 Fig. 5. 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.

Reproduction 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

Figure 1 .
Figure 1.Multi-scale modeling of myocardial mechano-energetic function.The model integrates previously developed and validated models of cardiomyocyte dynamics Tewari et al. (2016a,b), myocardial energetics Bazil et al. (2016); Gao 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).

Figure 2 .
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 (σ X B ) into a model of muscle mechanics is illustrated in panel B.
Parameters:Certain parameters from the cross-bridge and calcium-activation models ofTewari et al. (2016a,b)   andCampbell et al. (2018) were re-estimated to match data fromJanssen et al. (2002) on calcium transients and force-generation in isolated rat cardiac trabeculae.In brief, experiments were conducted at 37 • C. Calcium transients (Fig.3) were measured at different stimulation frequencies at fixed sarcomere length SL = 2.2 µm.Isometric tension time courses were measured at different stimulation frequencies and sarcomere lengths.Fig.3shows data on peak developed tension (T dev ) as a function of S L at stimulation frequency of 4 Hz; and data on relaxation time from peak to 50% of peak tension (RT 50); peak developed tension (T dev ) and time to peak tension (T T P ) as function of S L.

Figure 3 .
Figure 3. Crossbridge model parameters estimation.For S L = 1.9 µm the C a 50 = 5.89 and Hill coefficient n = 4.63 and for S L = 2.3 µm the C a 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).

Figure 4 .
Figure 4. Cardiovascular system (CVS) diagram.Adopted from Tewari et al.Tewari et al. (2016a,b).C P A , C PV , C S A , C Ao , and C SV represent lumped compliances of pulmonary arteries, pulmonary veins, systemic arteries, aorta, and systemic veins.R pul , R s y s , 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.
panels C and D of the original paper.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.20/21

Table 1 .
Energetics Model State Variables

Table 2 .
Energetics Model Reaction and Transport Fluxes

Table 3 .
Input arguments for cellular energetics model

Table 4 .
Output arguments for cellular energetics 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 Fig.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.6/21

Table 5 .
State variables in cross-bridge model

Table 7 .
Model parameters for cross bridge model

Table 8 .
State variables in TriSeg (Heart) model

Table 9 .
Parameters in TriSeg (Heart) model Table 10 lists the variables of the cardiovascular systems model.These variables and Eqs.21 and 22 are invoked in the MATLAB code Cardiovascularmechancis.m.