# Reduced left ventricular dynamics modeling based on a cylindrical assumption

Martin Genet, Jérôme Diaz, Dominique Chapelle, Philippe Moireau

► To cite this version:

Martin Genet, Jérôme Diaz, Dominique Chapelle, Philippe Moireau. Reduced left ventricular dynamics modeling based on a cylindrical assumption. International Journal for Numerical Methods in Biomedical Engineering, 2023, 10.1002/cnm.3711. hal-03832575v2

HAL Id: hal-03832575

https://hal.science/hal-03832575v2

Submitted on 22 May 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
SPECIAL ISSUE ARTICLE

# Reduced left ventricular dynamics modeling based on a cylindrical assumption

Martin Genet*<fnref n="1" />,<fnref n="2" /> | Jérôme Diaz<fnref n="2" />,<fnref n="1" /> | Dominique Chapelle<fnref n="2" />,<fnref n="1" /> | Philippe Moireau<fnref n="2" />,<fnref n="1" />

**Correspondence**

*Martin Genet, 1 rue Honoré d'Estienne d'Orves, 91120 Palaiseau, France. Email: martin.genet@polytechnique.edu

**Abstract**

Biomechanical modeling and simulation is expected to play a significant role in the development of the next generation tools in many fields of medicine. However, full-order finite element models of complex organs such as the heart can be computationally very expensive, thus limiting their practical usability. Therefore, reduced models are much valuable to be used, e.g., for pre-calibration of full-order models, fast predictions, real-time applications, etc.. In this work, focused on the left ventricle, we develop a reduced model by defining reduced geometry & kinematics while keeping general motion and behavior laws, allowing to derive a reduced model where all variables & parameters have a strong physical meaning. More specifically, we propose a reduced ventricular model based on cylindrical geometry & kinematics, which allows to describe the myofiber orientation through the ventricular wall and to represent contraction patterns such as ventricular twist, two important features of ventricular mechanics. Our model is based on the original cylindrical model of [Guc-cione, McCulloch, & Waldman 1991; Guccione, Waldman, & McCulloch 1993], albeit with multiple differences: we propose a fully dynamical formulation, integrated into an open-loop lumped circulation model, and based on a material behavior that incorporates a fine description of contraction mechanisms; moreover, the issue of the cylinder closure has been completely reformulated; our numerical approach is novel as well, with consistent spatial (finite element) and time discretizations. Finally, we analyse the sensitivity of the model response to various numerical and physical parameters, and study its physiological response.

**KEYWORDS:**

Reduced-order modeling; Continuum mechanics on manifold; Cardiac modeling; Computational mechanics

## 1 | INTRODUCTION

Biomechanical modeling and simulation is expected to play a significant role in the development of the next generation tools in many fields of medicine such as neurology [Vinje et al. 2019], pulmonology [Morton et al. 2018] and cardiology [Lee et al. 2014]. Indeed, it allows to structure and improve our understanding of the function of tissues and organs in health and disease, and can help to design better treatments and devices for various conditions. Moreover, it is a necessary step toward precision

<footnote n="1">¹LMS, École Polytechnique / CNRS / Institut Polytechnique de Paris, Palaiseau, France</footnote>

<footnote n="2">²Inria, MΞDISIM Team, Inria Saclay-Ile de France, Palaiseau, France</footnote>
medicine, where models are combined with clinical data through model personalization – or data assimilation – procedures in order to generate objective and quantitative diagnosis and/or prognosis information [Chabiniok et al. 2016; Patte et al. 2022].

However, such models can be computationally expensive, thus limiting their practical usability. For instance for the heart, cutting edge finite element models, even when restricted to the left ventricle, can require hours of computation for a single heart beat [Chapelle, Le Tallec, Moireau, & Sorine 2012; Quarteroni, Lassila, Rossi, & Ruiz-Baier 2017], and even more for parameter identification [Chabiniok et al. 2011; Genet et al. 2014]. Reduced models are much valuable to be used, e.g., for pre-calibration of full-order models, fast predictions, real-time applications such as critical care patient monitoring [Le Gall et al. 2020], etc.

Many approaches to derive reduced models have been developed. Organ-scale “lumped parameters” models have been proposed [Keener & Sneyd 2009], especially for whole-system simulations and analysis; they are, however, limited in their ability to integrate knowledge for instance of the tissue complex behavior. Other approaches consist in deriving a reduced model from a full-order model through model order reduction techniques [Brenner et al. 2021]; this is an active field of research, notably to build methods with controlled and interpretable reduction error. An intermediate approach, employed in this paper, consists in defining reduced geometry and kinematics, while keeping general motion and behavior laws, allowing to derive reduced models where all variables and parameters have a strong physical meaning.

This work is focused on the left ventricle, which is the most powerful of the heart’s four chambers, pumping blood toward the entire body. Spherical models have been proposed [Caruel, Chabiniok, Moireau, Lecarpentier, & Chapelle 2014; Hanna 1973], allowing to lump all spatial dimensions and reduce the ventricular dynamics into a simple system of temporal ordinary differential equations. They are, however, intrinsically limited in their ability to describe (i) the myofiber orientation through the thickness of the ventricular wall, and (ii) fine contraction patterns such as ventricular twist.

Twist is an important feature of the ventricular motion, which is induced by the distribution of myofiber orientation through the ventricular wall [Arts, Reneman, & Veenstra 1979; Taber, Yang, & Podszus 1996], and is strongly influenced by the ventricular shape [van Dalen et al. 2010]. It is key to ventricular function, notably in the intricate temporal coupling between, on the one hand, twisting and ejection, and, on the other hand between untwisting and filling [Rademakers et al. 1992; Sengupta, Tajik, Chandrasekaran, & Khandheria 2008]. It may also serve as an early indicator for various cardiac diseases and conditions [Castellanos et al. 2021; Wang, Khoury, Yue, Torre-Amione, & Nagueh 2008]. In practice, twist can be extracted from MRI or US imaging [Sengupta et al. 2008; Young & Cowan 2012], through motion tracking/image registration [Berberoğlu, Stoeck, Kozerke, & Genet 2022; Genet, Stoeck, von Deuster, Lee, & Kozerke 2018].

Reduced models based on cylindrical geometry and kinematics have thus been proposed in the literature [Arts et al. 1979; Guccione et al. 1991 1993]. Compared to spherical models, cylindrical models notably have the ability to better represent the fiber orientation distribution within the ventricular wall, and the ventricular torsion. However, cylindrical models are usually governed by partial differential equations as at least one spatial dimension remains, so that special care needs to be taken in order to perform efficient simulations.

In this paper, we propose a novel reduced ventricular model based on the cylindrical geometry and associated kinematics proposed by [Guccione et al. 1991 1993], while maintaining the constitutive ingredients of the full-order model of [Chapelle et al. 2012] and corresponding spherical reduction [Caruel et al. 2014; Manganotti, Caforio, Kimmig, Moireau, & Imperiale 2021], thus completing a tightly connected hierarchy of ventricular models. There are multiple major differences with the original cylindrical model of [Guccione et al. 1991 1993]: we propose a fully dynamical formulation, integrated into an open-loop lumped circulation model, and based on a material behavior that incorporates a fine description of the contraction mechanisms; moreover, the issue of the cylinder closure has been reformulated completely; our numerical approach is novel as well, with consistent spatial (finite element) and time discretizations.

The paper is organized as follows. In Section 2 we describe the model in details, as well as our proposed numerical solution approach. We start by recalling the kinematics proposed by [Guccione et al. 1991 1993] (Section 2.1), then present the formulation of the law of motion onto the manifold generated by the kinematics (Section 2.2), and recall the material model proposed by [Chapelle et al. 2012] (Section 2.3). Then, we specify the boundary conditions, including the internal ventricular pressure and coupling to general circulation, as well as the closure of the cylindrical ventricle (Section 2.4). We complete the continuous model description by establishing fundamental energy balance principles (Section 2.6), and recall the full variational formulation (Section 2.5). We then describe a consistent temporal (Section 2.7.1) and spatial (Section 2.7.2) discretization. In Section 3 we illustrate, analyze and discuss the model response. We start by studying the impact of numerical parameters (Section 3.1) to find optimal values, and then analyze the model response associated with “normal” physical parameters (Section 3.2). Then, we study the impact of various physical parameters, including the ventricular geometry (Section 3.3.1) and fiber orientation distribution (Section 3.3.2). Finally, in Section 4 we draw some conclusions and provide some perspectives to this work.
## 2 | METHODS: MODEL FORMULATION AND DISCRETIZATION

### 2.1 | Cylindrical geometry and kinematics

Following [Guccione et al. 1991 1993], we aim at representing the left ventricle as a closed cylinder of length L, internal radius $R_i$ and external radius $R_e$, as shown in Figure 1. The reference configuration is defined as $\Omega := \Omega_R \times \Omega_\Theta \times \Omega_Z$ with $\Omega_R := [R_i, R_e]$, $\Omega_\Theta := [0, 2\pi]$ and $\Omega_Z := [0, L]$. Its boundary is denoted by $\Gamma := \partial\Omega = \Gamma_i \cup \Gamma_e \cup \Gamma_+ \cup \Gamma_-$, with $\Gamma_i, \Gamma_e, \Gamma_+$ and $\Gamma_-$ its internal, external, top and bottom boundaries, respectively.

<figure><img src="image_2.png" /><figcaption>**FIGURE 1** Geometry and kinematics of the cylindrical model. Note that the cylinder is actually closed (see Figure 2), but the lids are not represented here. Note also that the mesh represented here is not a computational mesh (we will see that the model leads to a 1D only mathematical problem in the thickness of the ventricle), it only serves to illustrate the considered geometry and kinematics.</figcaption></figure>

The admissible deformation is constrained by the following map [Guccione et al. 1991]

$$ \psi(t) := \begin{cases} \Omega \to \omega(t) \\ (R, \Theta, Z) \mapsto \begin{cases} r = R + \rho(R, t) \\ \theta = \Theta + \beta(t) Z + \varphi(R, t) \\ z = (1 + \varepsilon(t)) Z + \eta(R, t) \end{cases} \end{cases}, \quad (1) $$
which is parametrized by three scalar fields $\rho$, $\varphi$ and $\eta$ belonging to some functional space $\mathcal{H}(\Omega_R)$ assumed to be regular enough and two additional scalars $\beta$ and $\varepsilon$. Basically, $\rho$ drives the circumferential deformation (e.g., the circumferential shortening during systole), while its derivative with respect to $R$ – which we denote by $\rho'$ – drives the radial deformation (e.g., the radial thickening during systole), and $\varepsilon$ the longitudinal deformation (e.g., the longitudinal shortening during systole). $\beta$ describes the global twist of the ventricle, i.e., the global ventricular torsion. $\varphi$ and $\eta$ describe more subtle deformation features, i.e., the in-plane and out-of-plane shear of cylindrical slices. Let us denote by

$$ \zeta := [\rho, \beta, \varphi, \varepsilon, \eta] \qquad (2) $$

the set of quantities characterizing the kinematics, belonging to the tensor product space denoted by

$$ \mathcal{V} := \mathcal{H}(\Omega_R) \times \mathbb{R} \times \mathcal{H}(\Omega_R) \times \mathbb{R} \times \mathcal{H}(\Omega_R), \qquad (3) $$

and henceforth denote the deformation map by $\psi[\zeta](R, \Theta, Z, t)$ to explicitly indicate the dependency in $\zeta$, which we will do for other quantities as well. For instance, for the internal volume of the cavity, we write

$$ V[\zeta](t) = \pi r_i(t)^2 l(t) = \pi (R_i + \rho(R_i, t))^2 (1 + \varepsilon(t)) L. \qquad (4) $$

In the sequel, every spatial entity will be expressed within the cylindrical coordinate system. We denote by $(e_R, e_\Theta, e_Z)$ the cylindrical basis of reference, and $(e_r, e_\theta, e_z)$ the cylindrical basis associated with the map $\psi$ (see Figure 1). Thus, a material point of cylindrical coordinates $(R, \Theta, Z)$ in the reference configuration, whose reference position in the reference basis (i.e., vector components in the $(e_R, e_\Theta, e_Z)$ basis) is then

$$ X = \begin{pmatrix} R \\ 0 \\ Z \end{pmatrix}, \qquad (5) $$

is transported by the map to the following deformed position, which can be expressed in both the deformed and reference basis (omitting the time dependency) as

$$ \begin{aligned} \mathbf{x}[\zeta](\mathbf{X}) &= \begin{pmatrix} r \\ 0 \\ z \end{pmatrix}_{(e_r, e_\Theta, e_z)} \\ &= \begin{pmatrix} \cos(\theta - \Theta) & -\sin(\theta - \Theta) & 0 \\ \sin(\theta - \Theta) & \cos(\theta - \Theta) & 0 \\ 0 & 0 & 1 \end{pmatrix}_{(e_R, e_\Theta, e_Z)} \otimes \begin{pmatrix} r \\ 0 \\ z \end{pmatrix}_{(e_r, e_\Theta, e_z)} \cdot \begin{pmatrix} r \\ 0 \\ z \end{pmatrix}_{(e_r, e_\Theta, e_z)} = \begin{pmatrix} \cos(\beta Z + \varphi(R))(R + \rho(R)) \\ \sin(\beta Z + \varphi(R))(R + \rho(R)) \\ (1 + \varepsilon) Z + \eta(R) \end{pmatrix}_{(e_R, e_\Theta, e_Z)}. \end{aligned} \qquad (6) $$

Consequently, the expression of the displacement field in the reference basis is

$$ \mathbf{u}[\zeta](\mathbf{X}) := \mathbf{x}[\zeta](\mathbf{X}) - \mathbf{X} = \begin{pmatrix} \cos(\beta Z + \varphi(R))(R + \rho(R)) - R \\ \sin(\beta Z + \varphi(R))(R + \rho(R)) \\ \varepsilon Z + \eta(R) \end{pmatrix}. \qquad (7) $$

We introduce similar notations (omitting both spatial and temporal dependencies) for the velocity (which is the push forward of $\dot{\zeta}$ by the kinematics) and acceleration fields such that

$$ \left\{ \begin{aligned} \mathbf{v}[\zeta] &:= \dot{\mathbf{u}}[\zeta] = \mathrm{D}\mathbf{u}[\zeta](\dot{\zeta}) \\ \mathbf{a}[\zeta] &:= \dot{\mathbf{v}}[\zeta] = \mathrm{D}^2\mathbf{u}[\zeta](\dot{\zeta}, \dot{\zeta}) + \mathrm{D}\mathbf{u}[\zeta](\ddot{\zeta}) \end{aligned} \right. , \qquad (8) $$

involving the first

$$ \mathrm{D}\mathbf{u}[\zeta](\hat{\zeta}) = \begin{pmatrix} \cos(\beta Z + \varphi) \hat{\rho} - \sin(\beta Z + \varphi) (\hat{\rho}Z + \hat{\varphi}) (R + \rho) \\ \sin(\beta Z + \varphi) \hat{\rho} + \cos(\beta Z + \varphi) (\hat{\rho}Z + \hat{\varphi}) (R + \rho) \\ Z \hat{\varepsilon} + \hat{\eta} \end{pmatrix} \quad \forall \hat{\zeta} \qquad (9) $$

and second derivative

$$ \mathrm{D}^2\mathbf{u}[\zeta](\bar{\zeta}, \zeta) = \begin{pmatrix} -\sin(\beta Z + \varphi) [\hat{\rho}(\bar{\beta}Z + \bar{\varphi}) + \bar{\rho}(\hat{\beta}Z + \hat{\varphi})] - \cos(\beta Z + \varphi)(R + \rho)(\hat{\beta}Z + \hat{\varphi})(\bar{\beta}Z + \bar{\varphi}) \\ \cos(\beta Z + \varphi) [\hat{\rho}(\bar{\beta}Z + \bar{\varphi}) + \bar{\rho}(\hat{\beta}Z + \hat{\varphi})] - \sin(\beta Z + \varphi)(R + \rho)(\hat{\beta}Z + \hat{\varphi})(\bar{\beta}Z + \bar{\varphi}) \\ 0 \end{pmatrix} \quad \forall \bar{\zeta}, \zeta \qquad (10) $$

of the displacement $\mathbf{u}$ constrained by the mapping $\psi$.
1 It is convenient to express the deformation gradient associated with the chosen kinematics in mixed basis, i.e., as the matrix
2 $(F_{iJ})_{i \in \{r, \theta, z\}}$ such as for any vector field $\mathbf{u} \in \mathbb{R}^3$ we have $F \cdot \mathbf{u} = \sum_{i \in \{r, \theta, z\}} \sum_{J \in \{R, \Theta, Z\}} F_{iJ} (e_J \cdot \mathbf{u}) e_i$. The differential of the
3 reference position is

$$
\delta \mathbf{X}(R, \Theta, Z; \delta R, \delta \Theta, \delta Z) := \lim_{h \to 0} \left( \frac{\partial}{\partial h} \left( \mathbf{X}(R + h \delta R, \Theta + h \delta \Theta, Z + h \delta Z) \right) \right) = \begin{pmatrix} \delta R \\ R \delta \Theta \\ \delta Z \end{pmatrix}, \quad (11)
$$

4 in the reference basis. For the deformed position x, its differential expressed in the deformed basis reads

$$
\delta \mathbf{x}(R, \Theta, Z; \delta R, \delta \Theta, \delta Z) = \begin{pmatrix} \delta r \\ r \delta \theta \\ \delta z \end{pmatrix} = \begin{pmatrix} r_{,R} \delta R + r_{,\Theta} \delta \Theta + r_{,Z} \delta Z \\ r (\theta_{,R} \delta R + \theta_{,\Theta} \delta \Theta + \theta_{,Z} \delta Z) \\ z_{,R} \delta R + z_{,\Theta} \delta \Theta + z_{,Z} \delta Z \end{pmatrix}. \quad (12)
$$

5 Recalling that by definition $F \cdot \delta X = \delta x$, we thus infer

$$
F[\zeta] = \begin{pmatrix} r_{,R} & \frac{r_{,\Theta}}{R} & r_{,Z} \\ r \theta_{,R} & \frac{r \theta_{,\Theta}}{R} & r \theta_{,Z} \\ z_{,R} & \frac{z_{,\Theta}}{R} & z_{,Z} \end{pmatrix} = \begin{pmatrix} 1 + \rho' & 0 & 0 \\ R \left(1 + \frac{\rho}{R}\right) \varphi' & 1 + \frac{\rho}{R} & R \left(1 + \frac{\rho}{R}\right) \beta \\ \eta' & 0 & 1 + \varepsilon \end{pmatrix}, \quad (13)
$$

6 in the above-specified mixed basis. Consequently, the volume ratio within the wall reads

$$
J[\zeta] := \det(F[\zeta]) = (1 + \rho')(1 + \frac{\rho}{R})(1 + \varepsilon). \quad (14)
$$

7 All the remaining tensor quantities will be expressed in the reference basis ($e_R$, $e_\Theta$, $e_Z$), including the right Cauchy-Green
defomation tensor
8

$$
C[\zeta] := F[\zeta]^T \cdot F[\zeta] = \begin{pmatrix} (1+\rho')^2 + (R(1+\frac{\rho}{R})\varphi')^2 + \eta'^2 & R(1+\frac{\rho}{R})^2 \varphi' (R(1+\frac{\rho}{R})^2) \varphi'\beta + \eta'(1+\varepsilon) \\ (1+\frac{\rho}{R})^2 & R(1+\frac{\rho}{R})^2 \beta \\ \text{sym} & (R(1+\frac{\rho}{R})\beta)^2 + (1+\varepsilon)^2 \end{pmatrix} \quad (15)
$$

9 and the Green-Lagrange strain tensor

$$
E[\zeta] := \frac{1}{2} (C[\zeta] - \mathbf{1}) = \frac{1}{2} \begin{pmatrix} (1+\rho')^2 + (R(1+\frac{\rho}{R})\varphi')^2 + \eta'^2 - 1 & R(1+\frac{\rho}{R})^2 \varphi' (R(1+\frac{\rho}{R})^2) \varphi'\beta + \eta'(1+\varepsilon) \\ (1+\frac{\rho}{R})^2 - 1 & R(1+\frac{\rho}{R})^2 \beta \\ \text{sym} & (R(1+\frac{\rho}{R})\beta)^2 + (1+\varepsilon)^2 - 1 \end{pmatrix}, \quad (16)
$$

10 illustrating the role of each parameter of the chosen kinematics in the cylinder deformation.

**2.2 | Dynamics equation**

11 As the deformation of the cylindrical model is parametrized in V, we need to write the dynamics equation on the manifold
12 of admissible deformations. For the sake of generality and conciseness, we propose to express the dynamical equilibrium in a
13 global, variational way, by invoking the principle of virtual power (a.k.a. principle of virtual work). We refer to Appendix A
14 for a complete presentation of our notation and of how we formally write the principle of virtual power onto the admissible
deformation manifold, extending the work of [Mardare 2015] to dynamics.
15
16 Thus, as detailed in Appendix A, in the context of the considered kinematically constrained reduced model the principle of
17 virtual work reads

18 Thus, as detailed in Appendix A, in the context of the considered kinematically constrained reduced model the principle of virtual work reads
19 The virtual power of acceleration forces is given by

$$
\text{Find } \zeta \text{ such that } \mathcal{P}_a[\zeta](\hat{\zeta}) + \mathcal{P}_i[\zeta](\hat{\zeta}) = \mathcal{P}_e[\zeta](\hat{\zeta}) \quad \forall \hat{\zeta}. \qquad (17)
$$

19 The virtual power of acceleration forces is given by

$$
\mathcal{P}_a[\zeta](\hat{\zeta}) := \int_{\Omega} \varrho_0 a[\zeta] \cdot \hat{v}[\zeta; \hat{\zeta}] \, d\Omega, \quad (18)
$$

20 where the acceleration has been given in (8), and the virtual velocity $\hat{v}$ is the push forward of the virtual parameter set $\hat{\xi}$ by the
21 kinematics

$$
\hat{v}[\zeta; \hat{\zeta}] := D\mathbf{u}[\zeta](\hat{\zeta}). \qquad (19)
$$
1 The virtual power of internal forces is given by

$$P_i[\zeta](\hat{\zeta}) := \int_{\Omega} \Sigma[\zeta] : D E[\zeta](\hat{\zeta}) \, d\Omega, \quad (20)$$

2 where the expression of the second Piola-Kirchhoff stress tensor $\Sigma[\zeta]$ as a function of strain and strain rate will be detailed in

3 Section 2.3, and the derivative of the Green-Lagrange strain tensor is given by

$$ \left\{ \begin{aligned} DE[\zeta](\hat{\zeta})_{RR} &= (1 + \rho') \hat{\rho}' + R \left(1 + \frac{\rho}{R}\right) \varphi'^2 \hat{\rho} + R^2 \left(1 + \frac{\rho}{R}\right)^2 \varphi' \hat{\varphi}' + \eta' \hat{\eta}' \\ DE[\zeta](\hat{\zeta})_{\Theta\Theta} &= \frac{1}{R} \left(1 + \frac{\rho}{R}\right) \hat{\rho} \\ DE[\zeta](\hat{\zeta})_{ZZ} &= R \left(1 + \frac{\rho}{R}\right) \beta^2 \hat{\rho} + R^2 \left(1 + \frac{\rho}{R}\right)^2 \beta \hat{\beta} + (1 + \varepsilon) \hat{\varepsilon} \\ DE[\zeta](\hat{\zeta})_{R\Theta} &= DE[\zeta](\hat{\zeta})_{\Theta R} = \left(1 + \frac{\rho}{R}\right) \varphi' \hat{\rho} + \frac{R}{2} \left(1 + \frac{\rho}{R}\right)^2 \hat{\varphi}' \\ DE[\zeta](\hat{\zeta})_{\Theta Z} &= DE[\zeta](\hat{\zeta})_{Z\Theta} = \left(1 + \frac{\rho}{R}\right) \beta \hat{\rho} + \frac{R}{2} \left(1 + \frac{\rho}{R}\right)^2 \hat{\beta} \\ DE[\zeta](\hat{\zeta})_{ZR} &= DE[\zeta](\hat{\zeta})_{RZ} = R \left(1 + \frac{\rho}{R}\right) \beta \varphi' \hat{\rho} + \frac{R^2}{2} \left(1 + \frac{\rho}{R}\right)^2 \varphi' \hat{\beta} + \frac{R^2}{2} \left(1 + \frac{\rho}{R}\right)^2 \beta \hat{\varphi}' + \frac{\eta'}{2} \hat{\varepsilon} + \frac{1+\varepsilon}{2} \hat{\eta}' \end{aligned} \right. . \quad (21) $$

4 The virtual power of external forces will be detailed in Section 2.4.

## 2.3 | Constitutive law

6 The myocardium is known to behave as a nonlinear, quasi-incompressible, viscous, active, anisotropic and heterogeneous solid [Humphrey 2002]. Many laws have been proposed in the literature [Holzapfel & Ogden 2009], focusing on various aspects of this complex behavior, and could be used in the context of the reduced model presented here. To illustrate the generality of our approach to reduced modeling, we propose to use the constitutive law proposed in [Chapelle et al. 2012], based on previous work from [Bestel, Clément, & Sorine 2001; Sainte-Marie, Chapelle, Cimrman, & Sorine 2006] and recently extended in [Kimmig, Chapelle, & Moireau 2019], which contains all the aforementioned behavioral specificities. We now briefly summarize its various ingredients.

13 We start by the symmetry class and the distribution of material orientation. The myocardium is known to be orthotropic [Dokos, Smaill, Young, & LeGrice 2002; Humphrey 2002], but is often approximated as transversely isotropic [Holzapfel & Ogden 2009; Tueni, Allain, & Genet 2022]. In this case, its local orientation, which varies throughout the ventricle, is defined by a single direction, i.e., the local myofiber orientation, which we denote by $e_F$, a unit vector. Following classical microstructural models [LeGrice, Hunter, & Smaill 1997], we assume that $e_F$ lies in the circumferential-longitudinal plane ($e_\Theta, e_Z$), and that its orientation depends only on the radial position $R$ in the form

$$e_F := \begin{pmatrix} 0 \\ \cos(\alpha(R)) \\ \sin(\alpha(R)) \end{pmatrix}, \quad (22)$$

19 where $\alpha(R) := \alpha_i \frac{R_e - R}{R_e - R_i} + \alpha_e \frac{R - R_i}{R_e - R_i}$ is a linear function of $R$.

20 The considered constitutive law is based on the rheology represented in Figure 3, which assumes that the passive and active components are in parallel, leading to the following decomposition of the stress (here, the second Piola-Kirchhoff stress tensor)

$$\Sigma = \Sigma^p + \Sigma^a. \quad (23)$$

22 The passive stress is itself assumed to be decomposed into three parts, namely the deviatoric, bulk and viscous parts

$$\Sigma^p = \Sigma^d + \Sigma^b + \Sigma^v. \quad (24)$$

23 Based on the previously mentioned symmetry considerations, the deviatoric stress is assumed to derive from the following transversely isotropic strain energy potential

$$\Sigma^d = \frac{\partial W^e}{\partial E} \quad \text{with} \quad W^e(\bar{I}_1, \bar{I}_2, \bar{I}_4) := C_1(\bar{I}_1 - 3) + C_2(\bar{I}_2 - 3) + C_3 e^{C_4(\bar{I}_1 - 3)^2} + C_5 e^{C_6(\bar{I}_4 - 1)^2}, \quad (25)$$

25 where $C_{1-6}$ are material parameters, and $\bar{I}_1, \bar{I}_2$ and $\bar{I}_4$ are reduced invariants of the transformation – namely $\bar{I}_1 := \bar{C} := \mathbf{1}$, $\bar{I}_2 := \frac{1}{2}(\bar{I}_1^2 - \bar{C}^2 := \mathbf{1})$ and $\bar{I}_4 := \bar{C} := (e_F \otimes e_F)$ with $\bar{C} := \bar{F}^\top \cdot \bar{F}, \bar{F} := J^{-1/3} \bar{F}, J := \det(\bar{F})$. Note that compared to
previous presentations of the model in [Caruel et al. 2014; Chapelle et al. 2012], simple neo-Hookean and Mooney-Rivlin terms were added to prevent ill-conditioning issues in the reference configuration, since the other terms have zero initial stiffness.

The myocardium is known to be slightly compressible [Humphrey 2002], but can be approximated as fully incompressible [Holzapfel & Ogden 2009]. In this case, the stress is not fully determined by the deformation; specifically, the hydrostatic pressure remains undetermined. We classically define

$$ \Sigma^b = -pJ\Sigma^{-1}, \qquad (26) $$

where $p$ is the Lagrange multiplier associated with the incompressibility constraint ($J = 1$), which is resolved through the following additional weak equation

$$ \int_{\Omega} (J - 1) \hat{p} \, d\Omega = 0 \quad \forall \hat{p}. \qquad (27) $$

Finally, the viscous behavior of the myocardium is less documented in the literature [Humphrey 2002]. As a consequence, we assume that the viscous part of the stress derives from the following simple isotropic dissipation pseudo-potential [Chapelle et al. 2012]

$$ \Sigma^v = \frac{\partial W^v}{\partial \dot{E}} \quad \text{with} \quad W^v := \frac{\gamma}{2} |\dot{E}|^2, \qquad (28) $$

where $\gamma$ is a material damping parameter.

Many models of the myocardial contraction have been proposed, from purely phenomenological [Guccione & McCulloch 1993] to multiscale [Kimmig et al. 2019], resulting in very different computational complexity, from simple explicit expressions [Guccione & McCulloch 1993] to very large systems of differential equations [Kimmig et al. 2019]. Here we consider the model developed in [Chapelle et al. 2012], which has a strong physical basis and good predictive power while requiring the resolution of two local ordinary differential equations only. Following the hypothesis that active forces are generated only along the local direction of the myofibers (although other hypotheses have been considered [Genet, Lee, et al. 2015; Guccione & McCulloch 1993; Guccione et al. 1993]), we have

$$ \Sigma^a = \sigma_{1D} e_F \otimes e_F, \qquad (29) $$

where, according to the rheology presented in Figure 3, we have

$$ \sigma_{1D} = \frac{k_s (e_{fib} - e_c)}{1 + e_{fib}}, \qquad (30) $$

where $k_s$ is the stiffness of the active component of the rheology, $e_c$ the local sarcomere deformation due to the actin-myosin binding, and $e_{fib}$ the fiber strain satisfying

$$ 1 + e_{fib} = (1 + 2 e_F \cdot E \cdot e_F)^{\frac{1}{2}} = \left( (1 + \varepsilon)^2 \sin^2 \alpha + \left( \beta(R + \rho) \sin \alpha + \frac{R + \rho}{R} \cos \alpha \right)^2 \right)^{\frac{1}{2}}. \qquad (31) $$

Moreover for the active contraction part,

$$ T_{fib} := k_s (e_{fib} - e_c) = \tau_c + \mu \dot{\epsilon}_c, \qquad (32) $$

where $\tau_c$ is the active stress developed by the actin-myosin bonds, and $\mu$ the dissipation of the active component of the rheology. Then, a distribution moment approximation applied to Huxley's formulation [Bestel et al. 2001; Chapelle et al. 2012; Kimmig et al. 2019; Zahalak 1981] allows to specify the dynamics of $\tau_c$ in terms of an ordinary differential equation linking $\tau_c$ to the active stiffness variable $k_c$ by

$$ \begin{cases} \dot{\tau}_c = -(|v| + \alpha|\dot{e}_c|) \tau_c + n_0(e_c) \sigma_0|v|_+ + k_c \dot{e}_c, \\ \dot{k}_c = -(|v| + \alpha|\dot{e}_c|) k_c + n_0(e_c) k_0|v|_+ \end{cases}, \qquad (33) $$

where $| \cdot |_+$ denotes the positive part. In (33), the parameter $k_0$ denotes the maximum active stiffness parameter and $\sigma_0$ the corresponding maximum active stress, $\alpha$ is a time constant, $n_0(e_c)$ a function accounting for the Frank-Starling mechanism, and $v([Ca^{2+}])$ a function triggering the contraction, typically when $[Ca^{2+}] > c_{th}$ with $c_{th}$ a given threshold.

This concludes the description of the material law chosen to illustrate the ventricular reduction principle using a cylindrical model. More details can be found in [Chapelle et al. 2012; Kimmig et al. 2019]. Even if our choice of constitutive law is motivated by the thermodynamical properties sustaining this macroscopic behavior model, any choice of passive, viscous and active law could be made in the context the proposed cylindrical ventricular model.
## 2.4 | Boundary conditions and loadings

The ventricular mechanics is tightly coupled to the blood flow and the vasculature mechanics. Since the focus of the current paper is on ventricular reduced modeling, we propose to couple it to a simple circulation model; note, however, that there is no theoretical limitation in that regard, and that reduced ventricular models can be coupled to highly complex circulation models.

Thus, here blood is only taken into account through a pressure applied onto the ventricle internal surface, denoted by $P_v$. This corresponds to the following virtual power term

$$P_e^i[\zeta](\hat{\xi}) := \int_{\Gamma_i} P_v e_R \cdot F^{-1} \cdot Du[\zeta](\hat{\xi}) J d\Gamma. \quad (34)$$

The ventricular pressure is regulated by reduced models of the valves and the general circulation. The valves are modeled through the following relationship between ventricular outgoing blood flow $Q := -\dot{V}$ and ventricular ($P_v$), atrial ($P_{at}$) and aortic ($P_{ar}$) pressures

$$ \left\{ \begin{array}{ll} \text{Filling:} & P_v \le P_{at} \Rightarrow Q = K_{at} (P_v - P_{at}) \\ \text{Isovol. contraction/relaxation:} & P_{at} \le P_v \le P_{ar} \Rightarrow Q = K_{iso} (P_v - P_{at}) \\ \text{Ejection:} & P_{ar} \le P_v \Rightarrow Q = K_{iso} (P_v - P_{at}) + K_{ar} (P_v - P_{ar}) \end{array} \right. , \quad (35) $$

where $K_{at}$, $K_{iso}$ and $K_{ar}$ are material parameters. Only the arterial circulation is taken into account, through a two-stage Windkessel model [Hunter, Pullan, & Smaill 2003]

$$ \left\{ \begin{array}{l} C_{\text{valve}} \dot{P}_v + \dot{V} + K_{\text{ar}} |P_v - P_{\text{ar}}|_+ + K_{\text{iso}} |P_v - P_{\text{at}}|_+ - K_{\text{at}} |P_{\text{at}} - P_v|_+ = 0 \\ C_{\text{ar}} \dot{P}_{\text{ar}} - K_{\text{ar}} |P_v - P_{\text{ar}}|_+ + \frac{P_{\text{ar}} - P_{\text{d}}}{R_{\text{p}}} = 0, \\ C_{\text{d}} \dot{P}_{\text{d}} + \frac{P_{\text{d}} - P_{\text{ar}}}{R_{\text{p}}} - \frac{P_{\text{vs}} - P_{\text{d}}}{R_{\text{d}}} = 0 \end{array} \right. , \quad (36) $$

where $C_{valve}$, $C_{ar}$, $R_p$, $C_d$ and $R_d$ are material parameters. The venous circulation is not modeled explicitly, but considered through the prescribed venal ($P_{vs}$) and atrial ($P_{at}$) pressures. More details can be found in [Sainte-Marie et al. 2006].

In order to be consistent with physiology, and represent phenomena such as longitudinal lengthening during diastole, ventricular pressure increase during systole, etc., the cylinder must be closed at both ends. Here we propose a formulation that is significantly different from that of [Guccione et al. 1991]. We consider that the bottom and top lids are infinitely thin, infinitely stiff in the longitudinal direction, and infinitely compliant in the radial direction, see Figure 2. Thus, the lids do not have any contribution to the virtual power of inertial or internal forces, and only contribute to the virtual power of external forces as follows

$$P_e^{\pm} := \int_{\gamma_{\pm}} -P_v \mathbf{n}^{\pm} \cdot \hat{v}^{\pm} d\gamma, \quad (37)$$

with $\mathbf{n}^{\pm} = \mp e_Z$. Because of symmetry conditions in the model, the lids remain flat and horizontal. Furthermore, because of the nature of the applied load, i.e., a normal pressure, we see that only the vertical component of the virtual displacement is present in the virtual power. Indeed, we have

$$P_e^{\pm} = \pm P_v \pi r_i^2 \hat{v}_Z^{\pm} \quad (38)$$

A choice must be made on how to connect the lids and the ventricle, i.e., how to express $\hat{v}_Z^{\pm}$ in terms of $\hat{v}[\zeta; \hat{\xi}]$. One possibility is to connect the lids to the internal line of the faces, i.e.,

$$\hat{v}_Z^{\pm} = \hat{v}_Z[\zeta; \hat{\xi}] (R = R_i) = \hat{\epsilon}Z + \hat{\eta}(R_i), \quad (39)$$

such that

$$\hat{v}_Z^{\pm} - \hat{v}_Z = \hat{\epsilon}L. \quad (40)$$

Another possibility would be to connect the lids to the average position of the faces, i.e.,

$$\hat{v}_Z^{\pm} = \frac{1}{|\gamma_{\pm}|} \int_{\gamma_{\pm}} \hat{v}_Z[\zeta; \hat{\xi}] d\gamma = \frac{\int_{\Gamma_{\pm}} (\hat{\epsilon}Z + \hat{\eta}) \|F^{-1} \cdot \mathbf{N}^{\pm}\| J d\Gamma}{\int_{\Gamma_{\pm}} \|F^{-1} \cdot \mathbf{N}^{\pm}\| J d\Gamma}. \quad (41)$$
1 However, because of the constrained kinematics, F and thus J do not depend on Z, so we have again

$$ \hat{v}_Z^+ - \hat{v}_Z^- = \hat{\varepsilon}L. \qquad (42) $$

2 Thus, the two considered options are actually equivalent. Finally, we obtain

$$ P_e^+ + P_e^- = P_v \pi L (R_i + \rho(R_i))^2 \hat{\varepsilon} \qquad (43) $$

3 and the complete virtual power of external forces reads

$$ P_e := P_e^i + P_e^+ + P_e^- \qquad (44) $$

4 Other types of boundary conditions could be considered, such as visco-elastic conditions on the external surface to model the presence of the pericardial fluid [Chabiniok et al. 2011; Pfaller et al. 2019], but will not be considered in this paper.

<figure><img src="image_3.png" /><figcaption>**FIGURE 2** Schematic of the cylindrical model closure. The lid is shown in purple (in reference and deformed configuration), its normal in blue; the displacement field in green.</figcaption></figure>

6 During the static simulations, the rigid body motions allowed by the kinematics, i.e., translation along $e_Z$ and rotation around

7 $e_Z$, must be blocked. To do so, we impose $\eta(R = R_i) = 0$ and $\varphi(R = R_i) = 0$. These constraints are not necessary in the

8 dynamic setting.
## 2.5 | Mathematical problem

The complete cylindrical cardiac model finally reads

Find $(\zeta, p, e_c, \tau_c, k_c, P_v, P_{ar}, P_d)$ such that

$$
\left\{
\begin{array}{ll}
\text{Kinematics: Section 2.1} & \\
\text{Equilibrium: } P_a[\zeta](\hat{\zeta}) + P_i[\zeta, p](\hat{\zeta}) = P_e[\zeta](\hat{\zeta}) \quad \forall \hat{\zeta} & \\
\text{with} & \\
\quad P_a[\zeta](\hat{\zeta}) = \int_{\Omega} \rho_0 \ddot{u}[\zeta] \cdot \mathbf{Du}[\zeta](\hat{\zeta}) \, d\Omega & \\
\quad P_i[\zeta, p](\hat{\zeta}) = \int_{\Omega} \Sigma[\zeta, p] : \mathbf{DE}[\zeta](\hat{\zeta}) \, d\Omega & \\
\text{Behavior: Section 2.3} & \\
\quad P_e[\zeta](\hat{\zeta}) = \int_{\Gamma_i} P_v e_R \cdot \mathbf{F}^{-1} \cdot \mathbf{Du}[\zeta](\hat{\zeta}) J \, d\Gamma + P_v \pi L (R_i + \rho(R_i))^2 \hat{\varepsilon} & \\
\text{Incompressibility: } \int_{\Omega} (J[\zeta] - 1) \hat{p} \, d\Omega = 0 \quad \forall \hat{p} & \\
\text{Valves: (35)} & \\
\text{Circulation: (36)} &
\end{array}
\right.
\tag{45}
$$

Note that we have kept the domain and boundary integrals only for the sake of conciseness. Indeed, thanks to the chosen parametrization the integrals can be simplified, and the problem (45) can be reduced to a problem in only one spatial dimension,

i.e. the radial dimension, see Appendix B for details on this integration. To simulate a full cardiac cycle (or multiple ones), we consider an unloaded, stress-free consideration (thus neglecting any pre-stress [Genet, Rausch, et al. 2015; Guccione et al. 1991]), load it statically to the internal pressure $P_{at}$ (corresponding to a end-diastole state), and then start the dynamic simulation.

The active contraction is driven by the time function $v$, which is prescribed as detailed in [Chapelle et al. 2012; Kimmig et al. 2019]. The different phases of the cardiac cycle are automatically handled by the model, notably through the switch function (35). Multiple cycles must be run in order to reach periodicity, although faster methods could be used [Khristenko & Le Tallec 2018].

## 2.6 | Energy balance

One of the fundamental properties of our formulation based on an Hellinger-Reissner Lagrangian is to be directly compatible with an energy balance similar to what was presented in the seminal work of Chapelle et al. [2012] but here with a constrained kinematics. For instance, by considering the actual velocity as a virtual velocity in the principle of virtual work (17), we obtain

$$ \frac{d}{dt} \left( \mathcal{K}[\zeta] + \int_{\Omega_0} \left[ \frac{E_s}{2}(e_{\text{fib}} - e_c)^2 + U_c \right] d\Omega + \mathcal{E}[\zeta] \right) = - \int_{\Omega_0} (|v| + \alpha |\dot{e}_c|) U_c \, d\Omega - \int_{\Omega_0} \mu(\dot{e}_c)^2 \, d\Omega - \int_{\Omega_0} \frac{\partial W_v}{\partial \dot{e}} : \dot{e} \, d\Omega + \int_{\Omega_0} n_0 U_0 |v|_+ \, d\Omega + \mathcal{P}_e u[\zeta], \quad (46) $$

with the energy contributions given by the kinetic energy and hyperelastic energy

$$ \mathcal{K} = \frac{1}{2} \int_{\Omega_0} \rho |v|^2 \, d\Omega \quad \text{and} \quad \mathcal{E} = \int_{\Omega_0} W_e \, d\Omega, $$

while along the fiber direction, we have an additional passive elastic energy $\frac{1}{2} \int_{\Omega_0} E_s (e_{\text{fib}} - e_c)^2 \, d\Omega$, but more fundamentally $U_c$ the microscopic elastic energy created by the actin-myosin bridges, solution of

$$ \dot{U}_c = -(|u| + \alpha |\dot{e}_c|)U_c + \dot{e}_c \tau_c + n_0 U_0 |v|_+ $$
and controlling the active stiffness and stress with $U_c \ge \tau_c^2/k_c$ as proven by Chapelle et al. [2012]. Moreover, the dissipative terms are

$$ \int_{\Omega_0} \frac{\partial W_v}{\partial \dot{\epsilon}} : \dot{\epsilon} \, d\Omega $$

while $\int_{\Omega_0} \mu(\dot{\epsilon}_c)^2 \, d\Omega$ is an additional viscous dissipation along the fiber, and $\int_{\Omega_0} (|u| + \alpha |\dot{\epsilon}_c|) U_c \, d\Omega$ the energy dissipated by actin-myosin bridges creation and destruction. Finally, the source terms are given by $\int_{\Omega_0} n_0 U_0 |u|_+ \, d\Omega$, i.e. the positive power associated with the *actin-myosin engine*, while $P_e$ is the power of external forces associated with the coupling with the external circulation. The energy exchanges between the ventricle and the circulation are conservative as well, since we have

$$ P_e[u|\zeta, v|\zeta] + P_v Q[\zeta] = 0 \qquad (47) $$

with $Q[\zeta] := -\dot{V}[\zeta]$ the flux of blood leaving the ventricle, i.e., entering the circulation. Following Manganotti et al. [2021], we have the *Valve + Windkessel energy balance*

$$ \frac{d}{dt} \left[ \frac{C_{\text{valve}}}{2} P_v^2 + \frac{C_{\text{ar}}}{2} P_{\text{ar}}^2 \right] = -K_{\text{ar}} |P_v - P_{\text{ar}}|_+^2 - K_{\text{iso}} |P_v - P_{\text{at}}|_+^2 - K_{\text{at}} |P_{\text{at}} - P_v|_+^2 + (K_{\text{at}} |P_{\text{at}} - P_v|_+ - K_{\text{iso}} |P_v - P_{\text{at}}|_+) P_{\text{at}} + P_v Q[\zeta] \quad (48) $$

where $P_v Q[\zeta]$ is, indeed, the exchange term with the cardiac cavity, $\frac{C_{\text{valve}}}{2} P_v^2$ and $\frac{C_{\text{ar}}}{2} P_{\text{ar}}^2$ are the energies stored in the capacitances, $K_{\text{at}} |P_v - P_{\text{at}}|_+ P_{\text{at}}$ and $K_{\text{iso}} |P_v - P_{\text{at}}|_+ P_{\text{at}}$ are source terms associated with the preload (valve open or small regurgitation), and all the other terms are dissipation due to the resistances.

These energy balances are critical, and special care will be taken to derive discretizations that preserve them at the discrete level, as detailed in Section 2.7.1.

## 2.7 | Numerical resolution

We now present our discretization choices, in terms of both temporal and spatial discretizations. Simulations presented in the paper were performed using the MATLAB code CardiacLab (InterDeposit Digital Number: IDDN.FR.001.470012.000.S.P.2016.000.31235), which is available upon request for academic collaborations.

### 2.7.1 | Time discretization

We first introduce time discretization. In terms of notation, we consider a sequence of time instants {$t^n$}$_{n \in \mathbb{N}}$ and define the associated time steps $\Delta t^n := t^{n+1} - t^n$. For known quantities, $\cdot^n$, $\cdot^{n+1}$ and $\cdot^{n+\frac{1}{2}}$ denote their values at time $t^n$, $t^{n+1}$ and $t^{n+\frac{1}{2}} := \frac{t^n+t^{n+1}}{2}$. For the unknowns of the problem, $\cdot^n$ denotes their value at time $t^n$ (which is considered known at a given time step), $\cdot^{n+1}$ their value at time $t^{n+1}$ (which must be determined), and we often use the notation $\cdot^{n+\frac{1}{2}} := \frac{\cdot^n+\cdot^{n+1}}{2}$.

A major objective of the proposed discretization scheme is to satisfy, at the discrete level, beside standard consistency and convergence properties, the energy balance holding at the continuous level, in the spirit of previous works in model reduction [Caruel et al. 2014; Chapelle et al. 2012; Manganotti et al. 2021]. For the internal variables associated with the active stress, as well as for the circulation variables, we rely on the discretization of [Chapelle et al. 2012], which is not recalled here, as we focus on the difficulty associated with the constrained kinematics. To start with, we propose to introduce an auxiliary velocity field $v$ over $\Omega$, and turn the second-order system (45) into the following first order system

$$ \text{Find } (\zeta, v, p) \text{ such that } \left\{ \begin{array}{l} P_a[\zeta, v](\hat{\zeta}) + P_j[\zeta, p](\hat{\zeta}) = P_e[\zeta](\hat{\zeta}) \quad \forall \hat{\zeta} \\ v = v[\zeta] \quad \forall X \\ \int_{\Omega} (J[\zeta] - 1) \hat{p} \, d\Omega = 0 \quad \forall \hat{p} \end{array} \right. \qquad (49) $$

where

$$ P_a[\zeta, v](\hat{\zeta}) := \int_{\Omega} \rho_0 \dot{v} \cdot \hat{v}[\zeta; \hat{\zeta}] \, d\Omega \qquad (50) $$
<figure><img src="image_4.png" /><figcaption>FIGURE 3 Schematic of the complete cylindrical model, including the active model rheology, the windkessel afterload model, and the valve model.</figcaption></figure>

1 We discretize this system using a mid-point scheme for the equilibrium equations and a backward scheme for the incompressibility equation, leading to the following discrete system

$$ \text{Find } (\xi^{n+1}, \boldsymbol{v}^{n+1}, p^{n+1}) \text{ such that } \left\{ \begin{array}{l} P_a^{n+\frac{1}{2}\#}[\zeta^{n+1}, \boldsymbol{v}^{n+1}](\hat{\zeta}) + P_i^{n+\frac{1}{2}\#}[\zeta^{n+1}, p^{n+1}](\hat{\zeta}) = P_e^{n+\frac{1}{2}\#}[\zeta^{n+1}](\hat{\zeta}) \quad \forall \hat{\zeta} \\ \boldsymbol{v}^{n+\frac{1}{2}} = \mathrm{D}\boldsymbol{u}[\zeta^{n+\frac{1}{2}}]\left(\frac{\zeta^{n+1}-\zeta^n}{\Delta t^n}\right) \quad \forall \boldsymbol{X} \\ \int_{\Omega} (J[\zeta^{n+1}] - 1) \hat{p} \, d\Omega = 0 \quad \forall \hat{p} \end{array} \right. \quad (51) $$

3 where $P_a^{n+\frac{1}{2}\#}[\zeta^{n+1}, \boldsymbol{v}^{n+1}](\hat{\zeta})$, $P_i^{n+\frac{1}{2}\#}[\zeta^{n+1}, p^{n+1}](\hat{\zeta})$ and $P_e^{n+\frac{1}{2}\#}[\zeta^{n+1}](\hat{\zeta})$ are now detailed.

4 Regarding the virtual power of acceleration forces, we propose

$$ P_a^{n+\frac{1}{2}\#}[\zeta^{n+1}, \boldsymbol{v}^{n+1}](\hat{\zeta}) := \int_{\Omega} \rho_0 \frac{\boldsymbol{v}^{n+1} - \boldsymbol{v}^n}{\Delta t^n} \cdot \mathrm{D}\boldsymbol{u}[\zeta^{n+\frac{1}{2}}](\hat{\zeta}) \, d\Omega. \quad (52) $$

5 This implies

$$ \frac{\Delta \mathcal{K}}{\Delta t^n} = P_a^{n+\frac{1}{2}\#}[\zeta^{n+1}, \boldsymbol{v}^{n+1}](\frac{\zeta^{n+1} - \zeta^n}{\Delta t^n}), \quad (53) $$

6 where $\Delta \mathcal{K} := \mathcal{K}^{n+1} - \mathcal{K}^n$ with $\mathcal{K}^n := \int_{\Omega} \rho_0 (\boldsymbol{v}^n)^2 \, d\Omega$. Thus, we have the first part of the discrete energy balance, i.e., the discrete counterpart to the continuous energy balance (46).
1 Regarding the virtual power of internal forces, we propose the following discretization

$$P_i^{n+\frac{1}{2}\#}[\zeta^{n+1}, p^{n+1}](\hat{\xi}) := \int_{\Omega} \Sigma^{n+\frac{1}{2}\#}[\zeta^{n+1}, p^{n+1}] : DE[\zeta^{n+\frac{1}{2}}](\hat{\xi}) \, d\Omega, \quad (54)$$

2 where

$$\Sigma^{n+\frac{1}{2}\#} := \Sigma_d^{n+\frac{1}{2}\#} + \Sigma_b^{n+\frac{1}{2}\#} + \gamma \dot{E}^{n+\frac{1}{2}\#} + \sigma_{1D}^{n+\frac{1}{2}\#} e_F \otimes e_F \quad (55)$$

3 with $\Sigma_d^{n+\frac{1}{2}\#}$, $\Sigma_b^{n+\frac{1}{2}\#}$ and $\sigma_{1D}^{n+\frac{1}{2}\#}$ adequate discretization of the passive and active stress tensors given below, and

$$\dot{E}^{n+\frac{1}{2}\#} := DE[\zeta^{n+\frac{1}{2}}](\frac{\zeta^{n+1} - \zeta^n}{\Delta t^n}). \quad (56)$$

4 Concerning the passive stress, one could follow the conservative approach proposed by [Gonzalez 2000; Hauret & Le Tallec 2006], leading to an exact discrete analog to the second part of the continuous energy balance (46). However, as the system is strongly dissipative, the conservative correction introduced in [Gonzalez 2000] can be omitted—see the discussion in [Hauret & Le Tallec 2006]. Thus, for the sake of simplicity, here we propose the following discretizations

8 and

$$\Sigma_d^{n+\frac{1}{2}\#} := \frac{\partial W_e}{\partial E}(E^{n+\frac{1}{2}\#}) \quad (57)$$

$$\Sigma_b^{n+\frac{1}{2}\#} := -p^{n+1} J^{n+\frac{1}{2}\#} C^{n+\frac{1}{2}\#}^{-1}. \quad (58)$$

9 Concerning the active stress, the discretization reproduces exactly what was proposed in [Kimmig et al. 2019], namely,

$$\sigma_{1D}^{n+\frac{1}{2}\#} := \frac{T_{\text{fib}}^{n+\frac{1}{2}\#}}{1 + 2 e_F \cdot E^n \cdot e_F} \frac{e_F \cdot E^{n+1} - E^n \cdot e_F}{e_F \cdot \Delta t^n \dot{E}^{n+\frac{1}{2}\#} \cdot e_F} \quad (59)$$

10 with

$$T_{\text{fib}}^{n+\frac{1}{2}\#} := \sqrt{k_c^{n+1}} \left( \frac{\tau_c}{\sqrt{k_c}} \right)^{n+1} + \mu \frac{e_c^{n+1} - e_c^n}{\Delta t^n} \quad (60)$$

11 and $\sqrt{k_c^n}$, $(\tau_c/\sqrt{k_c})^n$ and $e_c^n$ following the time-scheme proposed in [Chapelle et al. 2012]. As detailed in [Chapelle et al. 2012; Kimmig et al. 2019], this discretization leads to the following relation

$$\frac{\Delta \mathcal{E}}{\Delta t^n} = P_i^{n+\frac{1}{2}\#}[\zeta^{n+1}, \boldsymbol{v}^{n+1}](\frac{\zeta^{n+1} - \zeta^n}{\Delta t^n}), \quad (61)$$

13 where $\Delta \mathcal{E} := \mathcal{E}^{n+1} - \mathcal{E}^n$ with $\mathcal{E}^n$ the discrete elastic energy. Thus, combining (53) and (61), as well as the first equation in (51), we see that the proposed discretization conserves the energy balance (46) at the discrete level.

15 Finally, concerning the external loading, we propose

$$P_e^{n+\frac{1}{2}\#}[\zeta^{n+1}](\hat{\xi}) := \int_{\Gamma_i} P_v^{n+\frac{1}{2}} e_R \cdot F^{n+\frac{1}{2}\#}^{-1} \cdot \mathrm{D}u[\zeta^{n+\frac{1}{2}}](\hat{\xi}) \, d\Gamma + P_v^{n+\frac{1}{2}} \pi L \left( R_i + \rho^{n+\frac{1}{2}} (R_i) \right)^2 \hat{\varepsilon}, \quad (62)$$

16 which, considering the proposed kinematics, can be integrated as

$$P_e^{n+\frac{1}{2}\#}[\zeta^{n+1}](\hat{\xi}) = P_v^{n+\frac{1}{2}} 2\pi L \left(R_i + \rho^{n+\frac{1}{2}}(R_i)\right) \left(1 + \varepsilon^{n+\frac{1}{2}}\right) \hat{\rho} + P_v^{n+\frac{1}{2}} \pi L \left(R_i + \rho^{n+\frac{1}{2}}(R_i)\right)^2 \hat{\varepsilon}. \quad (63)$$

17 This leads to the following relation

$$P_e^{n+\frac{1}{2}\#}[\zeta^{n+1}](\frac{\zeta^{n+1} - \zeta^n}{\Delta t^n}) + P_v^{n+\frac{1}{2}} Q^{n+\frac{1}{2}\#}[\zeta^{n+1}] = 0 \quad (64)$$

18 with

$$Q^{n+\frac{1}{2}\#}[\zeta^{n+1}] := 2\pi L \left(R_i + \rho^{n+\frac{1}{2}}(R_i)\right) \left(1 + \varepsilon^{n+\frac{1}{2}}\right) \frac{\rho^{n+1}(R_i) - \rho^n(R_i)}{\Delta t^n} + \pi L \left(R_i + \rho^{n+\frac{1}{2}}(R_i)\right)^2 \frac{\varepsilon^{n+1} - \varepsilon^n}{\Delta t^n}. \quad (65)$$

19 Thus, the proposed scheme verifies the discrete counterpart to the energy balance between the ventricle and the circulation (47).
## 2.7.2 | Nonlinear iterations and spatial discretization

At each time step, we propose to solve system (51) in two stages. The first stage is obtained by introducing the second equation into the first, such that the system becomes

$$ \text{Find } (\xi^{n+1}, p^{n+1}) \text{ such that} \quad \left\{ \int_{\Omega} (J[\xi^{n+1}] - 1) \hat{p} \, d\Omega = 0 \quad \forall \hat{p} \right. \tag{66} $$

with

$$ \begin{aligned} \bar{P}_a^{n+\frac{1}{2}\#}[\xi^{n+1}](\hat{\xi}) &:= P_a^{n+\frac{1}{2}\#}[\xi^{n+1}, 2Du[\xi^{n+\frac{1}{2}}](\frac{\xi^{n+1}-\xi^n}{\Delta t^n}) - v^n](\hat{\xi}) \\ &= \int_{\Omega} 2\rho_0 \left( Du[\xi^{n+\frac{1}{2}}](\frac{\xi^{n+1}-\xi^n}{\Delta t^{n^2}}) - \frac{v^n}{\Delta t^n} \right) \cdot Du[\xi^{n+\frac{1}{2}}](\hat{\xi}) \, d\Omega, \end{aligned} \tag{67} $$

while the second stage is a simple local update of the velocity field, as follows

$$ v^{n+1} = 2Du[\xi^{n+\frac{1}{2}}](\frac{\xi^{n+1}-\xi^n}{\Delta t^n}) - v^n \quad \forall X. \tag{68} $$

The system (66) being highly nonlinear, we solve it through Newton-Raphson iterations. Integrals are simplified by analytically integrating over the circumference and length of the cylinder. For a complete description of the integrated system with the chosen time scheme, see Appendix B. Thus, the only remaining spatial coordinate is the radial position. We then discretize this 1D problem with finite elements. $\rho$, $\varphi$ and $\eta$ are discretized with standard Lagrange $P_k$ elements, and the pressure $p$ with $P_{k-1}$ elements so as to prevent any numerical locking potentially induced by the incompressibility constraint. $\beta$ and $\epsilon$ are simple scalars and do not need to be discretized. The auxiliary field $v$ is discretized at the quadrature point level, similarly to internal variables. In practice, we consider a single finite element across the cylinder thickness, albeit of high order. Moreover, a high quadrature degree is used, so as to adequately integrate the fields, including the fiber distribution, across the thickness of the cylinder.

# 3 | RESULTS AND DISCUSSION: MODEL RESPONSE AND SENSITIVITY

We now illustrate the response of the proposed reduced cylindrical model, and its sensitivity to various parameters. We start with the discretization parameters, in order to study the numerical convergence of the model response. We then focus on the impact of the various physical (i.e., geometrical, structural and mechanical) parameters on the model response.

## 3.1 | Sensitivity to numerical parameters

The three main parameters of the numerical scheme are (i) the interpolation order, controlling the spatial discretization, (ii) the quadrature order, determined by the number of integration points through the thickness of the cylindrical ventricle, and (iii) the time step, controlling the time discretization. Moreover, multiple cycles are required to obtain a periodic solution, so the number of simulated cycles is another input parameter. In order to analyze their impact on the solution, we ran the model for various values of these parameters and, taking the finest value (i.e., largest interpolation/quadrature order, smallest time step, final cardiac cycle) as reference, we computed multiple convergence metrics. Since the main quantities of interest of the model are ventricular pressure and volume, we defined normalized errors of such quantities over a cardiac cycle (duration denoted by $T$, discretized with $N$ time steps) as

$$ \left\{ \begin{aligned} \text{error}_p &:= \frac{\sqrt{\frac{1}{T} \int_0^T (P_v - P_{v,\text{ref}})^2} \, dt}{\sqrt{\frac{1}{T} \int_0^T (P_{v,\text{ref}})^2} \, dt} &\approx \sqrt{\frac{\frac{1}{N} \sum_{i=1}^N (P_v^n - P_{v,\text{ref}}^n)^2}{\frac{1}{N} \sum_{i=1}^N (P_{v,\text{ref}}^n)^2}} \\ \text{error}_v &:= \frac{\sqrt{\frac{1}{T} \int_0^T (V - V_{\text{ref}})^2} \, dt}{\sqrt{\frac{1}{T} \int_0^T (V_{\text{ref}})^2} \, dt} &\approx \sqrt{\frac{\frac{1}{N} \sum_{i=1}^N (V^n - V_{\text{ref}}^n)^2}{\frac{1}{N} \sum_{i=1}^N (V_{\text{ref}}^n)^2}} \end{aligned} \right. \tag{69} $$
<figure><img src="image_5.png" /><figcaption>FIGURE 4 Numerical convergence plots for the sensitivity of the model response (PV: work performed during a cardiac cycle, P: ventricular pressure evolution throughout a cardiac cycle, V: ventricular volume evolution throughout a cardiac cycle) to numerical parameters (upper left: interpolation order (using a quadrature order of 11 and a time step of 1 ms), upper right: quadrature order (using an interpolation order of 3 and a time step of 1 ms), lower left: time step (using an interpolation order of 3 and a quadrature order of 5), lower right: cardiac cycle index (using an interpolation order of 3, a quadrature order of 5 and a time step of 1 ms)). For each parameter, the model response with the finest value (i.e., upper left: interpolation order of 5, upper right: quadrature order of 11, lower left: time step of 0.125 ms, lower right: cardiac cycle index of 6) is taken as reference. All remaining simulations will be performed with an interpolation order of 3, a quadrature order of 5, a time step of 1 ms, and considering the solution of the fifth cardiac cycle, in order to maintain the discretization and periodicity errors on the main quantities of interest below 1 %.</figcaption></figure>

1 We also defined an error metric in the work produced by the ventricle onto the blood stream during a cardiac cycle, combining both blood pressure and blood flow data, as follows

$$ \text{error}_{PV} := \left| \frac{W - W_{\text{ref}}}{W_{\text{ref}}} \right| \quad \text{with} \quad \begin{cases} W := \int_{0}^{T} P_{v} \dot{V} dt \approx \sum_{i=1}^{N} P_{v}^{n+\frac{1}{2}} (V^{n+1} - V^{n}) \\ W_{\text{ref}} := \int_{0}^{T} P_{v,\text{ref}} \dot{V}_{\text{ref}} dt \approx \sum_{i=1}^{N} P_{v,\text{ref}}^{n+\frac{1}{2}} (V_{\text{ref}}^{n+1} - V_{\text{ref}}^{n}) \end{cases} \quad (70) $$

3 Note that because of the different normalizations used for the pressure, volume and pressure-volume errors, there is no strict ordering of the curves. Convergence results are shown in Figure 4. In order to have limited (< 1 %) discretization and periodicity errors while maintaining a fast computation time, all remaining simulations will be performed with an interpolation order of 3, a quadrature order of 5, a time step of 1 ms, and we will consider the solution of the fifth simulated cardiac cycle.
## 3.2 | Baseline response

Considering normal values of material parameters, as identified in previous works [Caruel et al. 2014; Chapelle et al. 2012; Kimmig et al. 2019] and reported in Table 1, we computed the baseline response of the model, shown in Figure 5. Pressure and volume temporal evolutions, as well as the pressure-volume curve, are highly physiological, and resemble the one produced by an existing reduced model based on a spherical geometry [Caruel et al. 2014]. The ventricular twist temporal evolution is a key output of the proposed cylindrical model that does not exist in the spherical counterpart. Many twist curves have been reported in the literature, most often measured through magnetic resonance [Kowallick et al. 2015; Rademakers et al. 1992] or ultrasound [Burns, La Gerche, Prior, & MacIsaac 2009; Sengupta et al. 2008] imaging and feature tracking. Globally, it was found that the ventricle starts twisting at the beginning of systole, during isovolumic contraction phase, then reaches a peak torsion towards the end of systole, and untwisting occurs during the isovolumic relaxation and early diastolic filling phases, which is consistent with our model prediction. Moreover, peak twist predictions are consistent with measured values [Kowallick et al. 2015; Sengupta et al. 2008]. Note that the negative initial torsion is simply due to the fact that the initial state of the cardiac cycle simulation corresponds to the end of filling, a state that is computed from the reference configuration (where there is no torsion at all) through an internal pressure loading.

However, most experimental curves are more peaked than our predicted twist, which shows a twist plateau, although the experimental peak is more or less pronounced from one paper to another, see for instance [Sengupta et al. 2008] (very pronounced) and [Kowallick et al. 2015] (less pronounced). This can be explained by the cylindrical assumption and associated boundary conditions simplifications, as well as the chosen kinematics in which the global ventricular rotation is arbitrary and only the rotation differential plays a role. Indeed, apical and basal rotations are actually quite different, which cannot be described by the current model. It is nevertheless interesting to see that apical rotation curves reported for instance in [Burns et al. 2009] show a plateau very similar to our model predictions. A subsequent analysis will be performed to characterize the influence of various model parameters on the twist response.

<table><caption>TABLE 1 Summary of physical model parameters, and values used for the baseline simulation.</caption><thead><tr><th></th><th></th><th></th><th></th></tr></thead><tbody><tr><td colspan="4"><strong>Geometry</strong></td></tr><tr><td>Reference ventricular chamber volume (V)</td><td>65 mL</td><td>Reference ventricular external radius (R<sub>e</sub>)</td><td>33.0 mm</td></tr><tr><td>Reference ventricular wall volume</td><td>130 mL</td><td>Reference ventricular thickness</td><td>13.9 mm</td></tr><tr><td>Reference ventricular aspect ratio</td><td>1.5</td><td>Reference ventricular length (L)</td><td>57.1 mm</td></tr><tr><td>Reference ventricular internal radius (R<sub>i</sub>)</td><td>19.0 mm</td><td></td><td></td></tr><tr><td colspan="4"><strong>Microstructure</strong></td></tr><tr><td>Myofiber helix angle at endocardium (α<sub>i</sub>)</td><td>+60°</td><td>Myofiber helix angle at epicardium (α<sub>e</sub>)</td><td>-60°</td></tr><tr><td colspan="4"><strong>Passive behavior</strong></td></tr><tr><td>Reference mass density (ρ<sub>0</sub>)</td><td>1 kg/L</td><td>C<sub>4</sub></td><td>2</td></tr><tr><td>C<sub>1</sub></td><td>7 Pa</td><td>C<sub>5</sub></td><td>50 Pa</td></tr><tr><td>C<sub>2</sub></td><td>0 Pa</td><td>C<sub>6</sub></td><td>4</td></tr><tr><td>C<sub>3</sub></td><td>700 Pa</td><td>Material viscosity (γ)</td><td>70 Pa</td></tr><tr><td colspan="4"><strong>Active behavior</strong></td></tr><tr><td>Spring stiffness (k<sub>s</sub>)</td><td>10<sup>8</sup> Pa</td><td>Maximum active stress (σ<sub>0</sub>)</td><td>65 10<sup>3</sup> Pa</td></tr><tr><td>μ</td><td>70 Pa.s</td><td>Time constant (α)</td><td>12.0</td></tr><tr><td>Maximum active stiffness (k<sub>0</sub>)</td><td>260 10<sup>3</sup> Pa</td><td>n<sub>0</sub> (e<sub>c</sub>)</td><td>[Chapelle et al. 2012]</td></tr><tr><td colspan="4"><strong>Valve model</strong></td></tr><tr><td>Atrial pressure (P<sub>at</sub>)</td><td>[Chapelle et al. 2012]</td><td>K<sup>-1</sup><sub>ar</sub></td><td>7.6923 10<sup>4</sup> Pa.s.m<sup>-3</sup></td></tr><tr><td>K<sup>-1</sup><sub>at</sub></td><td>1.1111 10<sup>5</sup> Pa.s.m<sup>-3</sup></td><td>K<sup>-1</sup><sub>iso</sub></td><td>2.0000 10<sup>9</sup> Pa.s.m<sup>-3</sup></td></tr><tr><td>C<sub>valve</sub></td><td>9.0000 10<sup>-9</sup> m<sup>3</sup>.Pa<sup>-1</sup></td><td></td><td></td></tr><tr><td colspan="4"><strong>Circulation model</strong></td></tr><tr><td>R<sub>p</sub></td><td>1.3500 10<sup>7</sup> Pa.s.m<sup>-3</sup></td><td>C<sub>ar</sub></td><td>2.1150 10<sup>-10</sup> m<sup>3</sup>.Pa<sup>-1</sup></td></tr><tr><td>R<sub>d</sub></td><td>1.0949 10<sup>8</sup> Pa.s.m<sup>-3</sup></td><td>C<sub>d</sub></td><td>2.0158 10<sup>-8</sup> m<sup>3</sup>.Pa<sup>-1</sup></td></tr></tbody></table>
<figure><img src="image_6.png" /></figure>

**FIGURE 5** Baseline response of the model, in terms of pressure (upper left) and volume (upper right) temporal evolutions, pressure-volume curve (lower left), as well as ventricular twist temporal evolution (lower right).

**3.3 | Sensitivity to physical parameters**

We now analyze the influence of various physical parameters on the model response, including geometrical, microstructural and material parameters.

**3.3.1 | Geometry**

Starting with geometrical parameters, Figures 6, 7 and 8 show the model response, in terms of ventricular pressure, volume and twist evolutions as well as pressure-volume curves, for various values of the ventricular aspect ratio (keeping internal and wall volumes at their baseline values), ventricular internal volume (keeping aspect ratio and wall volume at their baseline values) and ventricular wall volume (keeping aspect ratio and internal volume at their baseline values), respectively.

Regarding ventricular aspect ratio (i.e., length divided by internal diameter), it can be seen in Figure 6 that it has almost no influence on the model pressure and volume response. Ventricular torsion is proportional to ventricular aspect ratio, i.e., to the length of the ventricle, showing that the torsion per unit length is mostly driven by the material law and not affected by the ventricular shape. Such a direct impact of the ventricular length on ventricular torsion has been observed for instance in [van Dalen et al. 2010], although in vivo studies do not allow a strict uncoupling between ventricular aspect ratio and other structural parameters such as length or diameter. Modeling studies found that ventricular output was mostly influenced by wall curvature [Taber et al. 1996], which is not present in our straight cylindrical model.

Regarding ventricular volume, Figure 7 shows that it induces an expected shift in the volume and pressure-volume curves. The model also shows a slight impact on the peak systolic pressure and torsion, due to the fact that ventricles with smaller internal volume have smaller internal radius and length (to preserve aspect ratio) but larger wall thickness (to preserve wall volume).
<figure><img src="image_7.png" /><figcaption>FIGURE 6 Influence of the ventricular aspect ratio on the model response.</figcaption></figure>

Finally, regarding ventricular wall volume, Figure 8 shows an expected increase in peak systolic pressure as well as a decrease in end-systolic volume (hence ejection fraction) when increasing wall volume (i.e., ventricular mass). Ventricular twist, however, is not significantly affected, since only the wall thickness changes, not the internal radius or the ventricular length.

As a summary, these results show that the volume curve is mostly affected by internal radius, the pressure curve by ventricular thickness, and the twist curve by ventricular length.

### 3.3.2 | Fiber orientation

One key strength of the proposed cylindrical model compared to previous models based on the spherical symmetry assumption, is its ability to describe the myofiber orientation distribution throughout the ventricular wall. Figure 9 shows the drastic impact of the fiber orientation on the model response. Interestingly, the physiological distribution (+60° at endocardium, -60° at epi-cardium) allows to optimize peak systolic pressure as well as work produced by the ventricle, and leads to the second largest torsion.

## 4 | CONCLUSIONS AND PERSPECTIVES

In this paper we introduced a new level in our cardiac models hierarchy, with a reduced cylindrical model based on the geometry and kinematics proposed by [Guccione et al. 1991]. The series of models already contains a full-order 3D model [Chapelle et al. 2012] and a reduced “0D” spherical model [Caruel et al. 2014], and all models are based on a unique myocardial tissue constitutive framework [Chapelle et al. 2012]. For a still much reduced computational cost compared to the 3D model, the cylindrical model offers some interesting additional physiological features compared to the OD model, namely the myofiber orientation distribution and the ventricular torsion. Thus, the key strengths of the reduced-order model compared to the full-order model are its simplicity (all model parameters and variables have a clear, interpretable, meaning) and its efficiency (the model boils down to a 1D only mathematical problem). We proposed a fully consistent formulation of the model, including boundary conditions and the cylinder closure lid, in a variational framework based on the principle of virtual work written in the
<figure><img src="image_8.png" /><figcaption>FIGURE 7 Influence of the ventricular internal volume on the model response.</figcaption></figure>

1

2

3

4

manifold generated by the chosen kinematics. We also proposed an energy-preserving (except for the passive stress part, which is a common approximation for dissipative laws [Hauret & Le Tallec 2006]) discretization approach based on finite elements and a mid-point integration scheme. The model was shown to have a physiological response, including its sensitivity to various model parameters.

5

6

7

8

9

10

11

12

13

There are multiple perspectives to this work. In terms of modeling, we plan to investigate the cylindrical kinematics itself. Making $\beta$ and $\epsilon$ functions of $R$ (like $\rho$, $\varphi$ and $\eta$) should enrich the cylinder motion, potentially making its response more physiological. Similarly, we plan to encode better longitudinal dependency, for instance to be able to distinguish apical from basal rotations, and to make the model ventricular torsion response closer to measured ones. Additionally, more complex pericardial boundary conditions will be considered, and there impact characterized. We also plan to perform a systematic comparison between the responses of the 3D full-order model, the 1D cylindrical model and the OD spherical model. This will allow to characterize the intrinsic scaling between the various models, induced by the underlying structural hypothesis. In fine, this will pave the way for a hierarchical model parameter identification pipeline, where cheap reduced models are used in a first, exploratory step, and full-order models are only used in a finalization step.

14

<figure><img src="image_9.png" /></figure>

15

# APPENDIX

16

## A FROM VARIATIONAL PRINCIPLE TO THE PRINCIPLE OF VIRTUAL WORK

17

18

19

20

Inspired by the work of [Mardare 2015] about static elasticity in a Riemannian manifold, we want to propose a principle of virtual work adapted to our constrained kinematics formulation derived from a variational principle, albeit valid for elastodynamics. We restrict our formal justification to a passive hyperelastic material without any damping subjected to a uniform external pressure on its internal cavity. Then, a formal extension will be proposed for active and damped materials such as the cardiac tissue.
<figure><img src="image_10.png" /><figcaption>FIGURE 8 Influence of the ventricular wall volume on the model response.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIGURE 9 Influence of the fiber orientation distribution (fiber helix angle varies linearly from endocardium to epicardium, symmetrically with respect to mid-wall) on the model response.</figcaption></figure>
We consider in this Appendix that the map $\psi$ defines a transformation from the kinematics parametrization space to the Euclidian space $\mathbb{R}^3$, namely obliterating for now the cylindrical coordinates. To avoid any confusion, we denote this mapping $\psi$. We introduce the Hellinger-Reissner Lagrangian [Gravouil, Combescure, & Brun 2014; Washizu 1975] as

$$ \mathcal{L}(\zeta, \boldsymbol{v}) = \int_0^t \left[ \int_{\Omega} \rho_0 \left( \boldsymbol{v} \cdot (\mathrm{D}\boldsymbol{u} \cdot \zeta) - \frac{1}{2} |\boldsymbol{v}|^2 \right) \mathrm{d}\Omega - \mathcal{W}(E) + P_v \mathcal{V}_v \right] \mathrm{d}s. \quad (A1) $$

where $\mathcal{W}(E) = \int_{\Omega} W(E) \mathrm{d}\Omega$ is the hyperelastic potential, $P_v$ is the homogeneous pressure in the cavity and $\mathcal{V}_v$ the volume of the cavity delimited by $\Gamma_i$. Considering the stationarity of $\mathcal{L}$ with respect to $\zeta \in \mathcal{V}$ and $\boldsymbol{v} \in \mathcal{T}\mathcal{M}$, we find that

$$ \frac{d\boldsymbol{v}}{dt} = \mathrm{D}\boldsymbol{u} \cdot \dot{\zeta} \quad (A2) $$

and for all $\dot{\zeta} \in \mathcal{V}$,

$$ \int_{\Omega} \left[ \rho_0 \frac{d\boldsymbol{v}}{dt} \cdot (\mathrm{D}\boldsymbol{u} \cdot \dot{\zeta}) + \frac{\partial W(E)}{\partial E} : \mathrm{D}E \cdot \dot{\zeta} \right] \mathrm{d}\Omega = - \int_{\Gamma_i(t)} P_v \boldsymbol{n} \cdot (\mathrm{D}\boldsymbol{u} \cdot \dot{\zeta}) \mathrm{d}\Gamma \quad (A3) $$

where $\mathrm{D}\boldsymbol{u} \cdot \dot{\zeta}$ is the derivative of the displacement field $\boldsymbol{u}$ with respect to $\zeta$ and applied to a parametrization $\dot{\zeta} \in \mathcal{V}$, constrained by the mapping $\psi$.

In the same manner, $\mathrm{D}E \cdot \dot{\zeta}$ is the derivative of the Green-Lagrange strain tensor with respect to the parametrization applied to $\dot{\zeta}$.

When, the constitutive law contains an active part and a damping part, we therefore extend (A3) into a general principle of virtual work, namely,

$$ \forall \dot{\zeta} \in \mathcal{V}, \quad \int_{\Omega} \left[ \rho_0 \frac{d\boldsymbol{v}}{dt} \cdot (\mathrm{D}\boldsymbol{u} \cdot \dot{\zeta}) + \boldsymbol{\Sigma} : \mathrm{D}E \cdot \dot{\zeta} \right] \mathrm{d}\Omega = - \int_{\Gamma_i(t)} P_v \boldsymbol{n} \cdot (\mathrm{D}\boldsymbol{u} \cdot \dot{\zeta}) \mathrm{d}\Gamma \quad (A4) $$

where $\boldsymbol{\Sigma}$ is the second Piola-Kirchhoff stress tensor, namely the energy conjugate of the Green-Lagrange deformation tensor $E$.

## B ONE DIMENSIONAL REDUCTION OF THE MATHEMATICAL PROBLEM

The analytical integration along the $\Theta$ and $Z$ directions of the system described in (66) will be detailed hereafter.

Noting that the derivative of the mapping of the current position can be factorised as

$$ \mathrm{D}\boldsymbol{u}[\dot{\zeta}](\dot{\zeta}) = \begin{pmatrix} +\cos(\beta Z + \varphi) & -\sin(\beta Z + \varphi)Z(R+\rho) & -\sin(\beta Z + \varphi)(R+\rho) & 0 & 0 \\ +\sin(\beta Z + \varphi) & +\cos(\beta Z + \varphi)Z(R+\rho) & +\cos(\beta Z + \varphi)(R+\rho) & 0 & 0 \\ 0 & 0 & 0 & Z & 1 \end{pmatrix} \cdot \dot{\zeta}, \quad (B5) $$

the inertia forces can be rewritten as

$$ \begin{aligned} \bar{P}_a^{n+\frac{1}{2}\#}[\zeta^{n+1}](\dot{\zeta}) &= 2 \int_{\Omega} \rho_0 \dot{\zeta}^{\intercal} \cdot (\mathrm{D}\boldsymbol{u}[\zeta^{n+\frac{1}{2}}])^{\intercal} \cdot (\mathrm{D}\boldsymbol{u}[\zeta^{n+\frac{1}{2}}]) \cdot \left( \frac{\zeta^{n+1} - \zeta^n}{\Delta t^{n^2}} - \frac{\boldsymbol{v}^n}{\Delta t^n} \right) \mathrm{d}\Omega \\ &= 2 \int_{\Omega} \rho_0 \dot{\zeta}^{\intercal} \cdot (\mathrm{D}\boldsymbol{u}[\zeta^{n+\frac{1}{2}}])^{\intercal} \cdot (\mathrm{D}\boldsymbol{u}[\zeta^{n+\frac{1}{2}}]) \cdot \left( \frac{\zeta^{n+1} - \zeta^n}{\Delta t^{n^2}} - \rho_0 \dot{\zeta}^{\intercal} \cdot (\mathrm{D}\boldsymbol{u}[\zeta^{n+\frac{1}{2}}])^{\intercal} \frac{\boldsymbol{v}^n}{\Delta t^n} \right) \mathrm{d}\Omega. \end{aligned} \quad (B6) $$

Integrating over $\Theta$ and $Z$ gives

$$ \bar{P}_a^{n+\frac{1}{2}\#}[\zeta^{n+1}](\dot{\zeta}) = 4\pi \int_{R_i}^{R_e} \rho_0 \dot{\zeta}^{\intercal} \cdot \left( \mathcal{M} \cdot \left( \frac{\zeta^{n+1} - \zeta^n}{\Delta t^{n^2}} \right) - \mathcal{N}^{\intercal} \cdot \frac{\boldsymbol{v}^n}{\Delta t^n} \right) \mathrm{d}R, \quad (B7) $$

where

$$ \mathcal{M} = \begin{pmatrix} L & 0 & 0 & 0 & 0 \\ 0 & \frac{L^3}{3}(R+\rho^{n+\frac{1}{2}})^2 & \frac{L^2}{2}(R+\rho^{n+\frac{1}{2}})^2 & 0 & 0 \\ 0 & \frac{L^2}{2}(R+\rho^{n+\frac{1}{2}})^2 & L(R+\rho^{n+\frac{1}{2}})^2 & 0 & 0 \\ 0 & 0 & 0 & \frac{L^3}{3} & \frac{L^2}{2} \\ 0 & 0 & 0 & \frac{L^2}{2} & L \end{pmatrix}, \quad (B8) $$
and

$$ \mathcal{N} = \begin{pmatrix} \frac{1}{\beta^{n+\frac{1}{2}}} \left( \sin(\beta^{n+\frac{1}{2}} L + \varphi^{n+\frac{1}{2}}) - \sin \varphi^{n+\frac{1}{2}} \right) & \frac{R_\rho^{n+\frac{1}{2}}}{\beta^{n+\frac{1}{2}} \beta^2} \left( \sin \varphi^{n+\frac{1}{2}} - \sin(\beta^{n+\frac{1}{2}} L + \varphi^{n+\frac{1}{2}}) + \beta^{n+\frac{1}{2}} L \cos(\beta^{n+\frac{1}{2}} L + \varphi^{n+\frac{1}{2}}) \right) & \frac{R_\rho^{n+\frac{1}{2}}}{\beta^{n+\frac{1}{2}}} \left( \cos(\beta^{n+\frac{1}{2}} L + \varphi^{n+\frac{1}{2}}) - \cos \varphi^{n+\frac{1}{2}} \right) & 0 & 0 \\ \frac{1}{\beta^{n+\frac{1}{2}}} \left( \cos \varphi^{n+\frac{1}{2}} - \cos(\beta^{n+\frac{1}{2}} L + \varphi^{n+\frac{1}{2}}) \right) & \frac{R_\rho^{n+\frac{1}{2}}}{\beta^{n+\frac{1}{2}} \beta^2} \left( \cos(\beta^{n+\frac{1}{2}} L + \varphi^{n+\frac{1}{2}}) - \cos \varphi^{n+\frac{1}{2}} + \beta^{n+\frac{1}{2}} L \sin(\beta^{n+\frac{1}{2}} L + \varphi^{n+\frac{1}{2}}) \right) & \frac{R_\rho^{n+\frac{1}{2}}}{\beta^{n+\frac{1}{2}}} \left( \sin(\beta^{n+\frac{1}{2}} L + \varphi^{n+\frac{1}{2}}) - \sin \varphi^{n+\frac{1}{2}} \right) & 0 & 0 \\ 0 & 0 & 0 & \frac{L^2}{2} & L \end{pmatrix} \quad (B9) $$

For the internal forces, the spatial dependency is only along the radial direction as both the derivative of the Green-Lagrange strain tensor given in (21) and the second Piola-Kirchhoff stress tensor, defined through the chain rule involving the invariants of the Cauchy-Green dilation tensor (15), do not have a $\Theta$ or $Z$ dependency. Thus the integration simply reduces to

$$ \mathcal{P}_i^{n+\frac{1}{2}\#}[\zeta^{n+1}, p^{n+1}](\hat{\xi}) := 2\pi L \int_{R_i}^{R_e} \Sigma^{n+\frac{1}{2}\#}[\zeta^{n+1}, p^{n+1}] : \mathrm{DE}[\zeta^{n+\frac{1}{2}}](\hat{\xi}) \, dR. \quad (B10) $$

For the external forces, its integrated expression has already been given in (63).

The incompressibility constraint does not depend on $\Theta$ or $Z$ and simply reads

$$ \int_{\Omega} (J[\zeta^{n+1}] - 1) \hat{p} \, d\Omega = 2\pi L \int_{R_i}^{R_e} \left( (1 + \rho'^{n+1}) \left(1 + \frac{\rho^{n+1}}{R}\right) (1 + \varepsilon^{n+1}) - 1 \right) \hat{p} \, dR. \quad (B11) $$

Finally, the integrated system (66) reduces to

Find $(\zeta^{n+1}, p^{n+1})$ such that

$$ \left\{ \begin{array}{l} 4\pi \int_{R_i}^{R_e} \rho_0 \hat{\xi}^T \cdot \left( \mathcal{M} \cdot \left( \frac{\zeta^{n+1} - \zeta^n}{\Delta t^n} \right) - \mathcal{N}^T \cdot \frac{\mathbf{v}^n}{\Delta t^n} \right) \, dR + 2\pi L \int_{R_i}^{R_e} \Sigma^{n+\frac{1}{2}\#}[\zeta^{n+1}, p^{n+1}] : \mathrm{DE}[\zeta^{n+\frac{1}{2}}](\hat{\xi}) \, dR \\ \qquad = P_v^{n+\frac{1}{2}} 2\pi L \left( R_i + \rho^{n+\frac{1}{2}}(R_i) \right) \left(1 + \varepsilon^{n+\frac{1}{2}}\right) \hat{p} + P_v^{n+\frac{1}{2}} \pi L \left( R_i + \rho^{n+\frac{1}{2}}(R_i) \right)^2 \hat{\varepsilon}, \\ \\ 2\pi L \int_{R_i}^{R_e} \left( (1 + \rho'^{n+1}) \left(1 + \frac{\rho^{n+1}}{R}\right) (1 + \varepsilon^{n+1}) - 1 \right) \hat{p} \, dR = 0. \end{array} \right. \quad (B12) $$

## References

9. Arts, T., Reneman, R. S., & Veenstra, P. C. (1979, May 01). A model of the mechanics of the left ventricle. *Annals of Biomedical Engineering*, 7(3), 299–318. doi: 10.1007/BF02364118

10. Berberoğlu, E., Stoeck, C. T., Kozerke, S., & Genet, M. (2022). Quantification of Left Ventricular Strain by Joint Analysis of 3D Tagging and Cine MR Images. *In revision at Medical Image Analysis*.

11. Bestel, J., Clément, F., & Sorine, M. (2001). A Biomechanical Model of Muscle Contraction. In G. Goos, J. Hartmanis, J. van Leeuwen, W. J. Niessen, & M. A. Viergever (Eds.), *Medical Image Computing and Computer-Assisted Intervention - MICCAI 2001* (Vol. 2208, pp. 1159–1161). Berlin, Heidelberg: Springer Berlin Heidelberg. Series Title: Lecture Notes in Computer Science. doi: 10.1007/3-540-45468-3_143

12. Brenner, P., Grivet-Talocia, S., Quarteroni, A., Rozza, G., Schilders, W., & Silveira, L. M. (Eds.). (2021). *Model Order Reduction Volume 2: Snapshot-Based Methods and Algorithms*. De Gruyter. doi: 10.1515/9783110671490

13. Burns, A. T., La Gerche, A., Prior, D. L., & Maclsaac, A. I. (2009). Left Ventricular Untwisting Is an Important Determinant of Early Diastolic Function. *JACC: Cardiovascular Imaging*, 2(6), 709–716. doi: 10.1016/j.jcmg.2009.01.015

14. Caruel, M., Chabiniok, R., Moireau, P., Lecarpentier, Y., & Chapelle, D. (2014, Aug). Dimensional reductions of a cardiac model for effective validation and calibration. *Biomech Model Mechanobiol*, 13(4), 897–914.

15. Castellanos, D. A., Škardová, K., Bhattaru, A., Berberoğlu, E., Greil, G., Tandon, A., ... Chabiniok, R. (2021). Left Ventricular Torsion Obtained Using Equilibrated Warping in Patients with Repaired Tetralogy of Fallot. *Pediatric Cardiology*. doi: 10.1007/s00246-021-02608-y
1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

Chabiniok, R., Moireau, P., Lesault, P.-F., Rahmouni, A., Deux, J.-F., & Chapelle, D. (2011, jul). Estimation of tissue contrac- tility from cardiac cine-MRI using a biomechanical heart model. *Biomechanics and Modeling in Mechanobiology, 11*(5), 609–630. doi: 10.1007/s10237-011-0337-8

Chabiniok, R., Wang, V. Y., Hadjicharalambous, M., Asner, L., Lee, J., Sermesant, M., ... Nordsletten, D. A. (2016). Multi- physics and multiscale modelling, data–model fusion and integration of organ physiology in the clinic: Ventricular cardiac mechanics. *Interface Focus*, 6(2), 20150083. doi: 10.1098/rsfs.2015.0083

Chapelle, D., Le Tallec, P., Moireau, P., & Sorine, S. (2012). Energy-preserving muscle tissue model: formulation and compatible discretization. *International Journal for Multiscale Computational Engineering*, 10(2), 189–211. doi: 10.1615/intjmultcompeng.2011002360

Dokos, S., Smaill, B. H., Young, A. A., & LeGrice, I. J. (2002). Shear properties of passive ventricular myocardium. *American Journal of Physiology: Heart and Circulatory Physiology*, 283(6), H2650–9. doi: 10.1152/ajpheart.00111.2002

Genet, M., Lee, L. C., Ge, L., Acevedo-Bolton, G., Jeung, N., Martin, A. J., ... Guccione, J. M. (2015). A Novel Method for Quantifying Smooth Regional Variations in Myocardial Contractility Within an Infarcted Human Left Ventricle Based on Delay-Enhanced Magnetic Resonance Imaging. *Journal of Biomechanical Engineering*, 137(8). doi: 10.1115/1.4030667

Genet, M., Lee, L. C., Nguyen, R., Haraldsson, H., Acevedo-Bolton, G., Zhang, Z., ... Guccione, J. M. (2014). Distribution of normal human left ventricular myofiber stress at end diastole and end systole: A target for in silico design of heart failure treatments. *Journal of Applied Physiology*, 117, 142–152. doi: 10.1152/japplphysiol.00255.2014

Genet, M., Rausch, M. K., Lee, L. C., Choy, S., Zhao, X., Kassab, G. S., ... Kuhl, E. (2015). Heterogeneous growth-induced prestrain in the heart. *Journal of Biomechanics*, 48(10), 2080–2089. doi: 10.1016/j.jbiomech.2015.03.012

Genet, M., Stoeck, C. T., von Deuster, C., Lee, L. C., & Kozerke, S. (2018). Equilibrated Warping: Finite Element Image Registration with Finite Strain Equilibrium Gap Regularization. *Medical Image Analysis*, 50, 1–22. doi: 10.1016/j.media.2018.07.007

Gonzalez, O. (2000, January). Exact energy and momentum conserving algorithms for general models in nonlinear elasticity. *Computer Methods in Applied Mechanics and Engineering*, 190(13-14), 1763–1783.

Gravouil, A., Combescure, A., & Brun, M. (2014, December). Heterogeneous asynchronous time integrators for computational structural dynamics. *International Journal for Numerical Methods in Engineering*, 102(3-4), 202–232.

Guccione, J. M., & McCulloch, A. D. (1993). Mechanics of Active Contraction in Cardiac Muscle: Part I—Constitutive Relations for Fiber Stress that describe Deactivation. *Journal of Biomechanical Engineering*, 115(1), 72–81. doi: 10.1115/1.2895473

Guccione, J. M., McCulloch, A. D., & Waldman, L. K. (1991, February). Passive material properties of intact ventricular myocardium determined from a cylindrical model. *Journal of biomechanical engineering*, 113(1), 42–55.

Guccione, J. M., Waldman, L. K., & McCulloch, A. D. (1993). Mechanics of Active Contraction in Cardiac Muscle: Part II—Cylindrical Models of the Systolic Left Ventricle. *Journal of biomechanical engineering*, 115(1), 82–90. doi: 10.1115/1.2895474

Hanna, W. (1973). A Simulation of Human Heart Function. *Biophysical Journal*, 13(7), 603–621. doi: 10.1016/S0006-3495(73)86011-4

Hauret, P., & Le Tallec, P. (2006, January). Energy-controlling time integration methods for nonlinear elastodynamics and low-velocity impact. *Computer Methods in Applied Mechanics and Engineering*, 195(37-40), 4890–4916.

Holzapfel, G. A., & Ogden, R. W. (2009). Constitutive modelling of passive myocardium: A structurally based framework for material characterization. *Philosophical transactions of the Royal Society A: Mathematical, Physical & Engineering Sciences*, 367(1902), 3445–75. doi: 10.1098/rsta.2009.0091

Humphrey, J. D. (2002). *Cardiovascular Solid Mechanics: Cells, Tissues, and Organs*. New York: Springer-Verlag.

Hunter, P. J., Pullan, A. J., & Smaill, B. H. (2003). Modeling total heart function. *Annual Review of Materials Research*, 5, 147–77. doi: 10.1146/annurev.bioeng.5.040202.121537

Keener, J. P., & Sneyd, J. (2009). *Mathematical Physiology, Second Edition* (Springer ed.).

Khristenko, U., & Le Tallec, P. (2018). Newton-Krylov method for computing the cyclic steady states of evolution problems in nonlinear mechanics: Newton-Krylov method for computing the cyclic steady states of evolution problems in nonlinear mechanics. *International Journal for Numerical Methods in Engineering*, 116(3), 178–201. doi: 10.1002/nme.5920

Kimmig, F., Chapelle, D., & Moireau, P. (2019). Thermodynamic properties of muscle contraction models and associated discrete-time principles. *Advanced Modeling and Simulation in Engineering Sciences*. doi: 10.1186/s40323-019-0128-9

Kwalallick, J. T., Morton, G., Lamata, P., Jogiya, R., Kutty, S., Lotz, J., ... Schuster, A. (2015). Inter-study reproducibility of left
ventricular torsion and torsion rate quantification using MR myocardial feature tracking. *Journal of magnetic resonance imaging: JMRI*. doi: 10.1002/jmri.24979

Lee, L. C., Genet, M., Dang, A. B., Ge, L., Guccione, J. M., & Ratcliffe, M. B. (2014). Applications of computational modeling in cardiac surgery. *Journal of Cardiac Surgery*, 29(3), 293-302. doi: 10.1111/jocs.12332

Le Gall, A., Vallée, F., Pushparajah, K., Hussain, T., Mebazaa, A., Chapelle, D., ... Chabiniok, R. (2020, May). Monitoring of cardiovascular physiology augmented by a patient-specific biomechanical model during general anesthesia. A proof of concept study. *PLoS ONE*. Retrieved from https://hal.inria.fr/hal-02561128 doi: 10.1371/journal.pone.0232830

LeGrice, I. J., Hunter, P. J., & Smaill, B. H. (1997). Laminar structure of the heart: A mathematical model. *The American journal of physiology*, 272(5 Pt 2), H2466-76.

Manganotti, J., Caforio, F., Kimmig, F., Moireau, P., & Imperiale, S. (2021). Coupling reduced-order blood flow and cardiac models through energy-consistent strategies: Modeling and discretization. *Advanced Modeling and Simulation in Engineering Sciences*, 39.

Mardare, C. (2015). Static Elasticity in a Riemannian Manifold. In *Differential geometry and continuum mechanics* (pp. 307-342). Cham: Springer, Cham.

Morton, S., Dickson, J., Chase, J., Docherty, P., Desaive, T., Howe, S., ... Tawhai, M. (2018). A virtual patient model for mechanical ventilation. *Computer Methods and Programs in Biomedicine*, 165, 77-87. doi: 10.1016/j.cmpb.2018.08.004

Patte, C., Brillet, P.-Y., Fetita, C., Gille, T., Bernaudin, J.-F., Nunes, H., ... Genet, M. (2022). Estimation of regional pulmonary compliance in idiopathic pulmonary fibrosis based on personalized lung poromechanical modeling. *Journal of Biomechanical Engineering*. doi: 10.1115/1.4054106

Pfaller, M. R., Hörmann, J. M., Weigl, M., Nagler, A., Chabiniok, R., Bertoglio, C., & Wall, W. A. (2019, April). The importance of the pericardium for cardiac biomechanics: from physiology to computational modeling. *Biomechanics and Modeling in Mechanobiology*, 18(2), 503-529. doi: 10.1007/s10237-018-1098-4

Quarteroni, A., Lassila, T., Rossi, S., & Ruiz-Baier, R. (2017). Integrated Heart—Coupling multiscale and multiphysics models for the simulation of the cardiac function. *Computer Methods in Applied Mechanics and Engineering*, 314, 345-407. doi: 10.1016/j.cma.2016.05.031

Rademakers, F. E., Buchalter, M. B., Rogers, W. J., Zerhouni, E. A., Weisfeldt, M. L., Weiss, J. L., & Shapiro, E. P. (1992). Dissociation between left ventricular untwisting and filling. Accentuation by catecholamines. *Circulation*, 85(4), 1572-1581. doi: 10.1161/01.CIR.85.4.1572

Sainte-Marie, J., Chapelle, D., Cimrman, R., & Sorine, M. (2006). Modeling and estimation of the cardiac electromechanical activity. *Computers & Structures*, 84(28), 1743-1759. doi: https://doi.org/10.1016/j.compstruc.2006.05.003

Sengupta, P. P., Tajik, A. J., Chandrasekaran, K., & Khandheria, B. K. (2008). Twist Mechanics of the Left Ventricle. *JACC: Cardiovascular Imaging*, 1(3), 366-376. doi: 10.1016/j.jcmg.2008.02.006

Taber, L. A., Yang, M., & Podszus, W. (1996). Mechanics of ventricular torsion. *Journal of Biomechanics*, 29(6), 745-752. doi: 10.1016/0021-9290(95)00129-8

Tueni, N., Allain, J.-M., & Genet, M. (2022). On the structural origin of the anisotropy in the myocardium: Multiscale modeling and analysis. *Submitted to the Journal of the Mechanical Behavior of Biomedical Materials*.

van Dalen, B. M., Kauer, F., Vletter, W. B., Soliman, O. I. I., van der Zwaan, H. B., ten Cate, F. J., & Geleijnse, M. L. (2010). Influence of cardiac shape on left ventricular twist. *Journal of Applied Physiology*, 108(1), 146-151. doi: 10.1152/japplphysiol.00419.2009

Vinje, V., Ringstad, G., Lindstrøm, E. K., Valnes, L. M., Rognes, M. E., Eide, P. K., & Mardal, K.-A. (2019). Respiratory influence on cerebrospinal fluid flow – a computational study based on long-term intracranial pressure measurements. *Scientific Reports*, 9(1), 9732. doi: 10.1038/s41598-019-46055-5

Wang, J., Khoury, D. S., Yue, Y., Torre-Amione, G., & Nagueh, S. F. (2008). Preserved left ventricular twist and circumferential deformation, but depressed longitudinal and radial deformation in patients with diastolic heart failure. *European heart journal*, 29(10), 1283-9. doi: 10.1093/eurheartj/ehn141

Washizu, K. (1975). *Variational methods in elasticity and plasticity* (Vol. 3). Pergamon press Oxford.

Young, A. A., & Cowan, B. R. (2012). Evaluation of left ventricular torsion by cardiovascular magnetic resonance. *Journal of Cardiovascular Magnetic Resonance*, 14(1), 49. doi: 10.1186/1532-429X-14-49

Zahalak, G. I. (1981). A distribution-moment approximation for kinetic theories of muscular contraction. *Mathematical Biosciences*, 55(1-2), 89-114. doi: 10.1016/0025-5564(81)90014-6