Reduced left ventricular dynamics modeling based on a cylindrical assumption

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, for example, for pre‐calibration of full‐order models, fast predictions, real‐time applications, and so forth. 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 Guccione, 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 aswell, with consistent spatial (finite element) and time discretizations. Finally, we analyze the sensitivity of the model response to various numerical and physical parameters, and study its physiological response.

are combined with clinical data through model personalizationor data assimilationprocedures in order to generate objective and quantitative diagnosis and/or prognosis information. 4,5owever, 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, 6,7 and even more for parameter identification. 8,9Reduced models are much valuable to be used, for example, for precalibration of full-order models, fast predictions, real-time applications such as critical care patient monitoring, 10 and so forth.
Many approaches to derive reduced models have been developed.Organ-scale "lumped parameters" models have been proposed, 11 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 12 ; 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, 13,14 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, 15,16 and is strongly influenced by the ventricular shape. 17It 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. 18,19It may also serve as an early indicator for various cardiac diseases and conditions. 20,21In practice, twist can be extracted from MRI or US imaging, 19,22 through motion tracking/image registration. 23,24educed models based on cylindrical geometry and kinematics have thus been proposed in the literature. 15,25,26Compared 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 Refs.25,26, while maintaining the constitutive ingredients of the full-order model of Ref. 6 and corresponding spherical reduction, 13,27 thus completing a tightly connected hierarchy of ventricular models.There are multiple major differences with the original cylindrical model of Refs.25,26: 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 Ref. 25,26 (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 Ref. 6 (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.

| Cylindrical geometry and kinematics
Following Refs.25,26, 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 and Γ À its internal, external, top and bottom boundaries, respectively.
The admissible deformation is constrained by the following map 25 which is parametrized by three scalar fields ρ, φ and η belonging to some functional space ℋ Ω R ð Þ assumed to be regular enough and two additional scalars β and ε.Basically, ρ drives the circumferential deformation (e.g., the circumferential shortening during systole), while its derivative with respect to Rwhich we denote by ρ 0drives the radial deformation (e.g., the radial thickening during systole), and ε the longitudinal deformation (e.g., the longitudinal shortening during systole).β describes the global twist of the ventricle, that is, the global ventricular torsion.φ and η describe more subtle deformation features, that is, the in-plane and out-of-plane shear of cylindrical slices.Let us denote by the set of quantities characterizing the kinematics, belonging to the tensor product space denoted by and henceforth denote the deformation map by ψ ζ ½ R, Θ, Z,t ð Þto explicitly indicate the dependency in ζ, which we will do for other quantities as well.For instance, for the internal volume of the cavity, we write In the sequel, every spatial entity will be expressed within the cylindrical coordinate system.We denote by e R ,e Θ , e Z ð Þ the cylindrical basis of reference, and e r , e θ , e z ð Þthe cylindrical basis associated with the map ψ (see Figure 1).Thus, a material point of cylindrical coordinates R, Θ, Z ð Þin the reference configuration, whose reference position in the reference basis (i.e., vector components in the e R , e Θ , e Z ð Þbasis) is then 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 Consequently, the expression of the displacement field in the reference basis is We introduce similar notations (omitting both spatial and temporal dependencies) for the velocity (which is the push forward of _ ζ by the kinematics) and acceleration fields such that second derivatives of the displacement u constrained by the mapping ψ.
It is convenient to express the deformation gradient associated with the chosen kinematics in mixed basis, that is, as the matrix such as for any vector field u ℝ 3 we have in the reference basis.For the deformed position x, its differential expressed in the deformed basis reads Recalling that by definition F Á δX ¼ δx, we thus infer in the above-specified mixed basis.Consequently, the volume ratio within the wall reads All the remaining tensor quantities will be expressed in the reference basis e R , e Θ , e Z ð Þ , including the right Cauchy-Green deformation tensor and the Green-Lagrange strain tensor illustrating the role of each parameter of the chosen kinematics in the cylinder deformation.

| Dynamics equation
As the deformation of the cylindrical model is parametrized in V, we need to write the dynamics equation on the manifold of admissible deformations.For the sake of generality and conciseness, we propose to express the dynamical equilibrium in a global, variational way, by invoking the principle of virtual power (a.k.a.principle of virtual work).We refer to Appendix A.1 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 Ref. 28 to dynamics.
Thus, as detailed in Appendix A.1, in the context of the considered kinematically constrained reduced model the principle of virtual work reads The virtual power of acceleration forces is given by where the acceleration has been given in (8), and the virtual velocity b v is the push forward of the virtual parameter set b ζ by the kinematics The virtual power of internal forces is given by where the expression of the second Piola-Kirchhoff stress tensor Σ ζ ½ as a function of strain and strain rate will be detailed in Section 2.3, and the derivative of the Green-Lagrange strain tensor is given by The virtual power of external forces will be detailed in Section 2.4.

| Constitutive law
The myocardium is known to behave as a nonlinear, quasi-incompressible, viscous, active, anisotropic and heterogeneous solid. 29Many laws have been proposed in the literature, 30 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 Ref. 6 based on previous work from Refs.31,32, and recently extended in Ref. 33 which contains all the aforementioned behavioral specificities.We now briefly summarize its various ingredients.We start by the symmetry class and the distribution of material orientation.The myocardium is known to be orthotropic, 29,34 but is often approximated as transversely isotropic. 30,35In this case, its local orientation, which varies throughout the ventricle, is defined by a single direction, that is, the local myofiber orientation, which we denote by e F , a unit vector.Following classical microstructural models, 36 we assume that e F lies in the circumferential-longitudinal plane e Θ ,e Z ð Þ, and that its orientation depends only on the radial position R in the form R e ÀR i is a linear function of R. 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) The passive stress is itself assumed to be decomposed into three parts, namely the deviatoric, bulk and viscous parts Based on the previously mentioned symmetry considerations, the deviatoric stress is assumed to derive from the following transversely isotropic strain energy potential where C 1À6 are material parameters, and I 1 , I 2 and I 4 are reduced invariants of the transformationnamely I 1 ≔ C : 1, Note that compared to previous presentations of the model in Refs.6,13, simple neo-Hookean and Mooney-Rivlin terms were added to prevent illconditioning issues in the reference configuration, since the other terms have zero initial stiffness.
The myocardium is known to be slightly compressible, 29 but can be approximated as fully incompressible. 30In this case, the stress is not fully determined by the deformation; specifically, the hydrostatic pressure remains undetermined.We classically define where p is the Lagrange multiplier associated with the incompressibility constraint (J ¼ 1), which is resolved through the following additional weak equation Finally, the viscous behavior of the myocardium is less documented in the literature. 29As a consequence, we assume that the viscous part of the stress derives from the following simple isotropic dissipation pseudo-potential 6 where γ is a material damping parameter.Many models of the myocardial contraction have been proposed, from purely phenomenological 37 to multiscale, 33 resulting in very different computational complexity, from simple explicit expressions 37 to very large systems of differential equations. 33Here we consider the model developed in Ref. 6 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 26,37,38 ), we have where, according to the rheology presented in Figure 3, we have where k s is the stiffness of the active component of the rheology, e c the local sarcomere deformation due to the actinmyosin binding, and e fib the fiber strain satisfying Moreover for the active contraction part, where τ c is the active stress developed by the actin-myosin bonds, and μ the dissipation of the active component of the rheology.Then, a distribution moment approximation applied to Huxley's formulation 6,31,33,39 allows to specify the dynamics of τ c in terms of an ordinary differential equation linking τ c to the active stiffness variable k c by where j Á j þ denotes the positive part.In (33), the parameter k 0 denotes the maximum active stiffness parameter and σ 0 the corresponding maximum active stress, α is a time constant, n 0 e c ð Þ a function accounting for the Frank-Starling mechanism, and ν 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. 6,33Even 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.

| 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 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 ≔ À _ V and ventricular (P v ), atrial (P at ) and aortic (P ar ) pressures Filling : 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 40 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 Ref. 32.
In order to be consistent with physiology, and represent phenomena such as longitudinal lengthening during diastole, ventricular pressure increase during systole, and so forth, the cylinder must be closed at both ends.Here we propose a formulation that is significantly different from that of Ref. 25.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 with n AE ¼ ∓ e Z .Because of symmetry conditions in the model, the lids remain flat and horizontal.Furthermore, because of the nature of the applied load, that is, a normal pressure, we see that only the vertical component of the virtual displacement is present in the virtual power.Indeed, we have A choice must be made on how to connect the lids and the ventricle, that is, how to express b v Another possibility would be to connect the lids to the average position of the faces, that is, However, because of the constrained kinematics, F and thus J do not depend on Z, so we have again Thus, the two considered options are actually equivalent.Finally, we obtain and the complete virtual power of external forces reads 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, 8,41 but will not be considered in this paper.
During the static simulations, the rigid body motions allowed by the kinematics, that is, translation along e Z and rotation around e Z , must be blocked.To do so, we impose η R¼ R i ð Þ¼0 and φ R¼ R i ð Þ¼0.These constraints are not necessary in the dynamic setting.

| Mathematical problem
The complete cylindrical cardiac model finally reads Find ζ, p, e c ,τ c ,k c , P v , P ar , P d ð Þ such that Kinematics : Section 2:1 Behavior : Section 2:12:3 Incompressibility : 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, that is, the radial dimension, see Appendix A.2 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 25,42 ), load it statically to the internal pressure P at (corresponding to an end-diastole state), and then start the dynamic simulation.The active contraction is driven by the time function ν, which is prescribed as detailed in. 6,33The 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. 43

| 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 6 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 with the energy contributions given by the kinetic energy and hyperelastic energy while along the fiber direction, we have an additional passive elastic energy 1 2 dΩ, but more fundamentally U c the microscopic elastic energy created by the actin-myosin bridges, solution of and controlling the active stiffness and stress with U c ≥ τ 2 c =k c as proven by. 6Moreover, the dissipative terms are 2 dΩ is an additional viscous dissipation along the fiber, and Þ U c dΩ the energy dissipated by actin-myosin bridges creation and destruction.Finally, the source terms are given by R Ω 0 n 0 U 0 juj þ dΩ, that is, the positive power associated with by 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 where P v Q ζ ½ is, indeed, the exchange term with the cardiac cavity, C valve 2 P 2 v and C ar 2 P 2 ar are the energies stored in the capacitances, K at jP at À P v j þ P at and K iso jP v À P at j þ P 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.

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

| Time discretization
We first introduce time discretization.In terms of notation, we consider a sequence of time instants t n f g n ℕ and define the associated time steps Δt n ≔ t nþ1 À t n .For known quantities, : n , : nþ1 and : nþ 1 2 denote their values at time t n , t nþ1 and t nþ 1 2 ≔ t n þt nþ1 2 .For the unknowns of the problem, : n denotes their value at time t n (which is considered known at a given time step), : nþ1 their value at time t nþ1 (which must be determined), and we often use the notation : nþ 1 2 ≔ : n þ: 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. 6,13,27For the internal variables associated with the active stress, as well as for the circulation variables, we rely on the discretization of, 6 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 Ω, and turn the second-order system (45) into the following first order system Find ζ, v, p ð Þsuch that where 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 Find ζ nþ1 , v nþ1 , p nþ1 À Á such that where

Regarding the virtual power of acceleration forces, we propose
This implies where . Thus, we have the first part of the discrete energy balance, that is, the discrete counterpart to the continuous energy balance (46).
Regarding the virtual power of internal forces, we propose the following discretization where with Σ 1D adequate discretization of the passive and active stress tensors given below, and Concerning the passive stress, one could follow the conservative approach proposed by, 44,45 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 44 can be omitted-see the discussion in. 45Thus, for the sake of simplicity, here we propose the following discretizations and Concerning the active stress, the discretization reproduces exactly what was proposed in, 33 namely, with and e n c following the time-scheme proposed in Refs.6.As detailed in Refs.6,33, this discretization leads to the following relation where ΔE ≔ E nþ1 À E n with 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.Finally, concerning the external loading, we propose which, considering the proposed kinematics, can be integrated as This leads to the following relation with Thus, the proposed scheme verifies the discrete counterpart to the energy balance between the ventricle and the circulation (47).

| 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 Find ζ nþ1 , p nþ1 À Á such that with while the second stage is a simple local update of the velocity field, as follows 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 2. Thus, the only remaining spatial coordinate is the radial position.We then discretize this 1D problem with finite elements.ρ, φ and η 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.β and ε 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.

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

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

| Baseline response
Considering normal values of material parameters, as identified in previous works Refs.6,13,33 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. 13The 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 18,46 or ultrasound 19,47 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 toward 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. 19,46Note 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 19 (very pronounced) and 46 (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 Ref. 47 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.
F I G U R E 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%.

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

| Geometry
Starting with geometrical parameters, Figures 6, 7 (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, that is, 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 Ref. 17, 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, 16 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).
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 curve by ventricular thickness, and the twist curve by ventricular length.

| 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 epicardium) allows to optimize peak systolic pressure as well as work produced by the ventricle, and leads to the second largest torsion.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 Ref. 25.The series of models already contains a full-order 3D model 6 and a reduced "0D" spherical model, 13 and all models are based on a unique myocardial tissue constitutive framework. 6For 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 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 45 ) 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.
There are multiple perspectives to this work.In terms of modeling, we plan to investigate the cylindrical kinematics itself.Making β and ε functions of R (like ρ, φ and η) 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 their 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.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.

ORCID
Martin Genet https://orcid.org/0000-0003-2204-201XF I G U R E 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.

ℒ ζ, v
ð Þ¼ where W E ð Þ¼ R Ω W E ð ÞdΩ is the hyperelastic potential, P v is the homogeneous pressure in the cavity and V v the volume of the cavity delimited by Γ i .Considering the stationarity of ℒ with respect to ζ V and v T ℳ, we find that and for all b ζ V, where Du Á b ζ is the derivative of the displacement field u with respect to ζ and applied to a parametrization b ζ V, constrained by the mapping ψ.
In the same manner, DE Á b ζ is the derivative of the Green-Lagrange strain tensor with respect to the parametrization applied to b ζ.
When, the constitutive law contains an active part and a damping part, we therefore extend (A3) into a general principle of virtual work, namely, where Σ is the second Piola-Kirchhoff stress tensor, namely the energy conjugate of the Green-Lagrange deformation tensor E.

A.2 | ONE DIMENSIONAL REDUCTION OF THE MATHEMATICAL PROBLEM
The analytical integration along the Θ 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 factorized as the inertia forces can be rewritten as Δt n dΩ: Integrating over Θ and Z gives where and 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 Θ or Z dependency.Thus the integration simply reduces to For the external forces, its integrated expression has already been given in (63).The incompressibility constraint does not depend on Θ or Z and simply reads Finally, the integrated system (66) reduces to.Find ζ nþ1 ,p nþ1 À Á such that 4π

AEZ
in terms of b v ζ; b ζ h i .One possibility is to connect the lids to the internal line of the faces, that is, F I G U R E 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. b Schematic of the complete cylindrical model, including the active model rheology, the windkessel afterload model, and the valve model.

F U R E 5
Baseline response of the model, in terms of pressure (upper left) and volume (upper right) temporal evolutions, pressurevolume curve (lower left), as well as ventricular twist temporal evolution (lower right).

T A B L E 1
Summary of physical model parameters, and values used for the baseline simulation.at endocardium (α i ) +60Myofiber helix angle at epicardium (α e ) 10 À8 m 3 :Pa À1 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 F G U R E 7 Influence of the ventricular internal volume on the model response.FI G U R E 6 Influence of the ventricular aspect ratio on the model response.

F I U R E 8
Influence of the ventricular wall volume on the model response.