E ﬀ ect of ﬁbre orientation and bulk modulus on the electromechanical modelling of human ventricles

: This work concerns the mathematical and numerical modeling of the heart. The aim is to enhance the understanding of the cardiac function in both physiological and pathological conditions. Along this road, a challenge arises from the multi-scale and multi-physics nature of the mathematical problem at hand. In this paper, we propose an electromechanical model that, in bi-ventricle geometries, combines the monodomain equation, the Bueno-Orovio minimal ionic model, and the Holzapfel-Ogden strain energy function for the passive myocardial tissue modelling together with the active strain approach combined with a model for the transmurally heterogeneous thickening of the myocardium. Since the distribution of the electric signal is dependent on the ﬁbres orientation of the ventricles, we use a Laplace-Dirichlet


Introduction
The heart has the role of pumping deoxigenated blood to the lungs to get oxygen and, simultaneously, delivers blood rich of all sort of vital substances to tissues and organs through the arterial circulatory system [32]. Cardiovascular diseases are the leading cause of death in the global population, with millions of new patients every year [31]. Despite new experimental discoveries and improvements in the medical care, we are still distant from fully understanding these pathologies. Computational modelling has shown to be a potential mean of addressing this problem. Recent achievements in this field provide additional insights for the understanding of cardiac functions and disease treatments. However, many difficulties arise when approximating the modelling of cardiac function (including e.g., the mismatch of model parameters and spatio-temporal scales, associated to the difficult task of retrieving in-vivo measurements from the tissue) and when trying to simulate their joint electromechanical behaviour as a coupled multi-physics and multi-scale problem [22]. Advances in the fields of experimental and theoretical biology, physics, computer science and clinical data allow the improvement of the mathematical models, allowing numerical simulations in patient-specific framework [13,28,29]. Moreover, application of image segmentation techniques to Magnetic Resonance Imaging (MRI) and Computed Tomography (CT) scans provide data to proceed with personalized computational modelling and patient-specific diagnoses and therapies.
The simulations of human heart functioning involve several challenges intrinsically related to its complexity [11,12,16]. At the moment, several cardiac modelling studies at the whole heart level are still restricted to simulating particular components, such as, e.g., the electrophysiology or the electromechanics. Only few works presented a fully coupled model taking into account the integrated electro-mechano-fluid behaviour [17,44].
Heart physiology can be briefly summarized as follows: an electric potential propagates across the membrane of the cardiac muscle cells (cardiomyocytes) and induces complex biochemical reactions inside the cytosol that releases calcium from the sarcoplasmic reticulum, thus resulting in the generation of force within the sarcomeres (the basic contractile units within cardiac muscle cells), finally causing the individual cells to contract and the muscle to deform. The contraction of the muscle yields a rapid increase of pressure inside the ventricular cavities, which allows the heart valves to open and close in precise sequence and induces the periodic filling and ejection of blood from the ventricles and the atria [15,33,35,36,[45][46][47]. This physiological process underline an intrinsic multi-scale nature. Indeed, the solutions of the corresponding mathematical models feature a wide range of spatio-temporal scales, which need to be accurately captured by computationally efficient numerical methods. As a matter of fact, ion channels on the cell membrane and the excitation-contraction mechanism are typically modelled through systems of ODEs to be solved for each individual cell (length scales of µm) and are endowed with time scales of 10 −4 ms. On the other hand, the description of fluid dynamics and solid mechanics of the tissue at the organ level (length scales of cm) and of the propagation of the electrical signal (monodomain) are translated in non-linear partial differential equations of either parabolic or parabolic-hyperbolic type with time scales of 0.1 ms (for the fluid dynamics and solid mechanics) [36].
In this work, we analyse the effect of fibre direction and bulk modulus on the pressure-volume relationship through electromechanical modelling for both the left and right ventricles [4]. This is an extension of the model presented in [19,20] that was however applied to the left ventricle only. We use state-of-the-art models in passive myocardial tissue modelling (the Holzapfel-Ogden model [23]) together with the active strain approach [1,2] in combination with a recently proposed model for the transmurally heterogeneous thickening of the myocardium [5]; the latter is used in the integrated electromechanics context for both the left (LV) and the right (RV) human ventricles. Once established, the active mechanics is coupled with the electrophysiology through a model describing the shortening of the myocardial fibres [41], which is in turn triggered by a change in the ionic concentrations in the cardiac cells. We use a Laplace-Dirichlet Rule-Based algorithm [6] to reproduce the ventricles fibres and sheets configuration.
From the numerical point of view, we discretize in space the models by means of the Finite Element Method (FEM) with piecewise linear polynomials of degree one (P 1 ) for both the myocardium displacement and the transmembrane potential, while the time discretization is carried out by means of Backward Differentiation Formulas (BDFs) of order 2 [37]. We then exploit a fully segregated algebraic problem, where the discretized core models are solved sequentially. This segregated strategy was proposed in [20] for the LV through a Godunov splitting scheme.
In this work, we successfully use a segregated scheme for the simulation of the full heartbeat, by realizing pressure-volume loops patient-specific bi-ventricles.
We propose an innovative electro-mechanical model of human ventricles employing an original modelling of transversal and isotropic mechanical activation. The fluid-structure interaction (FSI) between the blood contained in each chamber and the endocardial wall is addressed through a simple 0D model (spatially independent) for the pressure variable tailored for the different phases of the heartbeat [14,39,49]; a prestress phase is also implemented for both the LV and the RV in order to estimate the internal stresses of the myocardium at the initial time [25,48], i.e., in the reference configuration of the muscle. It is relevant to notice that, at telediastole (the final phase of the ventricular diastole), the endocardial pressure is significantly different in the LV and RV chambers, thus this method has to be able to compute the response of the tissue to the fluid in the corresponding chamber. All the solvers are implemented in the open source finite element library LifeV (www.lifev.org) [18].
The simulations are carried out on a bi-ventricle geometry extracted from a full heart atlas derived from multiple patients [27]. A set of various fibre orientation and bulk modulus is tested to show the effect of these changes on the mechanical behaviour of the ventricles. We show the obtained simulated pressure-volume loops to assess the influence on the various phases of a full heartbeat. Finally, the numerical results are compared against physiological data and clinical measurements. In our work, we show the influence of fibre direction on the electro-mechanical computational modelling of human ventricles.
The article is organized as follows: in Section 2 we introduce the mathematical models for the electrophysiology, the mechanics and the activation of the myocardium; we then couple them to obtain a integrated model. We then proceed in Section 3 presenting the characteristics of the cardiac cycle and the various phases performed by the ventricles composing a single heartbeat. In Section 4 we carry on the space and time numerical approximations of the single core models and their numerical coupling. In Section 5 we describe the method used to generate the bi-ventricle geometry and the fibres and sheets directions. We report and discuss in Section 6 the numerical results obtained with the proposed methods.
Finally, we draw our conclusions in Section 7.

Mathematical modelling
In this section we introduce the mathematical models used to represent the complex behaviour of the bi-ventricle electro-mechanics.

Ionic model and monodomain equation
Assuming the same electrophysiological behaviour for both left and right ventricles, we consider the monodomain equation for the description of the evolution of the cellular transmembrane potential V [12,20,26,34,43]. The system of differential equations modelling the electrophysiology reads: Here, Ω 0 is the reference computational domain (represented e.g., by the configuration of the ventricles at the end of the diastolic phase) and T > 0 is the final time of our simulation. The components w j , with j ∈ {1, . . . , M} are M so-called gating variables. F = I + ∇ 0 d is the strain tensor, d the displacement of the tissue and N is the normal vector. In order to take into account the anisotropic electrical conductance [42], we define the diffusion tensor as D m = σ t I + (σ l − σ t ) f 0 ⊗ f 0 , where σ t , σ l ∈ R + are the electric conductivities in the directions transversal and longitudinal with respect to the fibres, respectively. f 0 is the local fiber orientation, while s 0 is the direction longitudinal to the sheet orientation in the reference configuration. The function I app (t) represents an externally applied current, which stands for the electric stimulus injected at the endocardium by the terminal fibres of the Purkinje network; for our purposes, we consider it as a source term which triggers the electrophysiological activity. The terms I ion (V, w), w ∞ , α(V), β(V) are peculiar of the ionic model which, in our case, following the approach of [20], is the Bueno-Orovio minimal model [10].

Passive and active mechanics
We assume, in this section, the same constitutive laws for both the ventricles to describe both the active and passive mechanics of the tissue. We will highlight later differences and settings for the RV and LV. We recall the momentum conservation equation in the reference configuration Ω 0 , endowed with boundary and initial conditions, in the unknown displacement variable d [30]: We denote by K η ⊥ , K η , C η ⊥ , C η ∈ R + the parameters of generalized Robin conditions on each of these boundary subsets: the symbols ⊥ and identify either a parameter relative to the normal or the tangential direction, respectively. These boundary conditions are set to mimic the effect of the pericardium and surrounding tissues on the biventricular geometry. ρ is the density of the myocardium, γ f denotes the local shortening of the fibres, d 0 andḋ 0 are initial conditions. We notice that pressures p endo,RV (t) and p endo,LV (t) are still prescribed at this stage; we will explain in the Section 2.3 the manner in which they are determined. Γ η 0 , with η = {epi, base}, represents the subsets of the boundary corresponding to the epicardium and the base of the myocardium, while Γ endo,LV 0 and Γ endo,RV 0 are, respectively, the left and right endocardium, as depicted in Figure 1. Following [20], we model the myocardium as an hyperelastic material, thus: ∂F . For more details on the strain-energy density function used, we refer the interested reader to [20]. We emphasise that our final W is obtained by adding to the strain energy function a convex term W vol (J), where J = det(F v ) = det(F), being det(F) = 1, used to weakly enforce the incompressibility constraint, such that J = 1 is its global minimum; penalizing large variations in volume. We choose: hence, in this quasi-incompressible formulation, we obtain a "stronger" enforcement of the incompressibility constraint with a larger bulk modulus B ∈ R + . We proceed with the active strain approach [40] to take into account the active behaviour of the myocardial muscle. This formulation consists in a particular decomposition of the strain tensor: where F E is the isochoric component of the elastic (passive) part of the deformation, F A corresponds instead to the active part. We recall the orthotropic form for the tensor F A Moreover, to take into account the fact that the thickening of the ventricles' walls is transversely nonhomogeneous, we use the formulation proposed in [5], that is: Since the cardiomyocites stretching is driven by the sarcomeres dynamics due to the concentration of the calcium ions (here denoted as w 3 ), we will use the fibres' shortening model [19]: where µ A represent a quantity to be properly tuned for the case under consideration. For the explicit form of the term Φ(w 3 , γ f , d) we refer to [19].

Prestress
The bi-ventricle geometries in the configuration at t = 0 will not be stress-free since we assume that Ω 0 is at telediastole. Therefore, we have to take into account the pressures p endo,LV (t) and p endo,RV (t) on each endocardium wall exerted by the blood. We recall that these quantities are always larger than zero during the cardiac cycle. In fact, from medical measurements (e.g., [9,38]), one obtains the approximate values in Table 1, corresponding to an healthy individual. Solving problem (2.2) with physiological endocardial pressures p endo,LV > 0 and p endo,RV > 0 would give rise to non-physiological displacements as the internal stresses are not in equilibrium with the intraventricular blood's pressure of the ventricles. We indicate with p endo,LV and p endo,RV the pressure at telediastole and the stressed ventricles' configuration is determined in these conditions. To address this issue, we adapt the method proposed in [19], an extension of the one presented in [25], called pressure prestress [48]. We compute an internal stresses distribution such that the reference geometry is in equilibrium with the blood pressures p endo,LV and p endo,RV both in the left and in the right ventricles. With this aim we proceed with an additive decomposition of the strain tensor P = P(d) + P 0 , where the prestress tensor P 0 is determined to ensure a null displacement d 0 in correspondence of the assigned pressures p endo,LV and p endo,RV . In order to compute P 0 according to this approach, we proceed by adapting the method proposed in [25] to our model by first defining the following mechanical problem: Moreover, we set where S ∈ N is the number of steps of the continuation method that we exploit to gradually increase the pressure value inside each chamber. We recall that, in our case, the endocardium is subdivided into two ventricles so we need to initialize our pressures taking physiological values for the right and left chamber at telediastole. We can do that since we have a particular label for each ventricle. We set p endo,LV = 10 mmHg and p endo,RV = 5 mmHg. Eq (2.4) is the steady, passive version of problem (2.2) with decomposition of the strain tensor. Then, the fixed point iteration described in Algorithm 1 is applied to compute P 0 . First, we compute the approximation P 0 = Prestress(100, 10 −2 , p endo,LV , p endo,RV , 0) and finally we set P 0 = Prestress(1, 10 −5 , p endo,LV , p endo,RV , P 0 ). The additional step is performed to require a smaller tolerance when the pressures have already reached a closer value of the tensor P 0 to the target. We observe that, while does not converge to 0 but to a vector which we denote by d and will be our initial state for the displacement. However, we observe that the quantity d ∞ is negligible with respect to the endocardial walls thickness, and that this initial displacement ensures that the prestress is in equilibrium with the pressure at the epicardium. Therefore, we set d 0 = d in (2.2). P 0 = P m 0,k ; 8: update P m+1 0,k = P m 0,k + P d m+1 k , I ; 9: set m = m + 1; 10: set P 0,k = P m 0,k ; 12: end for 13: return P 0,S 14: end function

Cardiac cycle
For our simulations we will consider a full heartbeat, by taking the conventional duration of T = 0.8 s. With this aim, we take into account for the fluid-structure interaction of the endocardium with the blood by modeling the pressures p endo,LV and p endo,RV as in Eq (2.4). Before introducing the models used to describe the behaviour of the blood, we first recall that the heartbeat of each ventricle is conventionally split into four phases (see Figure 2): (1). an isovolumic contraction phase (the red one in Figure 2) in which the endocardial pressures p endo,RV (t) and p endo,LV (t) increase from the End Diastolic Pressure (EDP) to the value measured at the pulmonary artery or at the aorta. This increment is driven by the early stages of the ventricles contraction; (2). an ejection phase (the one in yellow) characterized by a decrement in the ventricular volumes V endo,LV (t) and V endo,RV (t) due to the contraction of the ventricles. It is called ejection phase since the contraction forces the blood to flow out from the ventricular chambers; (3). an isovolumic relaxation (green one) phase during which p endo,LV (t) and p endo,RV (t) decrease as a consequence of the ventricles early relaxation; (4). a filling phase (purple one) starting with the opening of the tricuspid valve or of the mitral valve causing an increment of the endocardial volume until p endo,LV (t) and p endo,RV (t) reach the EDP value. Before proceeding it is essential to evaluate the behaviour and interaction of the ventricles. The LV and RV are asynchronous, moreover the timing of each phase is strictly dependent on the ventricle, as we can see in Figure 3. Finally, we can rewrite the four phases depicted before to describe the coupled dynamics of the ventricles: (1). we begin with the phase in which both ventricles are in the isovolumic contraction phase. The right ventricular free wall shortens and moves toward the septum, the left ventricular chamber compresses and shortens [7]; (2). since the pulmonary valve opens before the aortic one (see Figure 3), the RV enters the ejection phase meanwhile the LV is still in the isovolumic phase; (3). the opening of the aortic valve determines the beginning of the ejection phase also for the LV. At this moment, both the chambers are ejecting blood. The ventricles keep contracting until there is fluid within; (4). the first valve to close is the aortic one (see Figure 3) and thus the LV enters in the second isovolumic phase, meanwhile the RV continues the ejection; (5). when the pulmonary valve closes too both ventricles are in the isovolumic relaxation; (6). since the tricuspid valve opening occurs before than the mitral one (see Figure 3), the RV is the first to enter the filling phase; (7). during the filling phase both pressures p endo,LV (t) and p endo,RV (t) are increasing until reaching the EDP value.

Numerical discretization
We concisely present the numerical discretization of each single continuous core models introduced in Section 2. We start with the space semi-discrete formulation and we will end up with the full time discretization. We perform our numerical coupling through a segregated algorithm strategy, instead of a monolithic approach.

Space discretization
As already noticed, we will approximate each of the single core models in the computational domain with a space discretization based on the FEM. We consider a mesh composed of tetrahedrons T h , with h representing the maximum size of the elements K ∈ T h , such that ∪ K∈T h K = Ω 0 ; the mesh elements are pairwise disjoint and their union Ω 0 ⊂ R 3 is the region of the space identified by the myocardium at the telediastolic phase of the heartbeat. We denote by N dof V , N dof w , N dof d , and N dof γ f the number of Degrees of Freedom (DoF) (i.e., the size of the discretized single core model) for the potential, ionic variables, displacement, and fibre shortening, respectively. The underlying number of nodes determined by the mesh T h is indicated as N h . Firstly we introduce the finite dimensional space is the set of polynomials of degree equal to 1 in the element K.
We denote with x j N dof w j=1 the set of the degrees of freedom and the value of the l-th ionic variable in x j by w l j (t). Similarly, we write V j (t) = V h (x j , t), and finally rearrange the unknowns in the vector w h (t) as w h (t) = w l h (t) The projection of the solution V(t) on the finite element space X 1 h can hence be written as V h (t) = N dof V j=1 V j (t)ψ j and the weak semidiscrete formulation of the problem reads: given w h (t) and d h (t), find, for all t ∈ (0, T ), V h (t) such that w + j, for l = 1, . . . , N I , j = 1, . . . , N dof w . We will use in Eq (4.1) a lumped mass matrix M L in place of M in order to avoid numerical oscillations as presented in [19].
To proceed with the FEM approximation of the momentum equation (2.2) we denote by [X 1 h ] 3 the finite dimensional subspace of vector valued functions and by {ψ i } N do f d i=1 its basis. Recall that, in our case, the endocardium splits right and left ventricle chambers, as defined in Figure 1, so we have to take into account two different pressure values. Our semi-discrete formulation for the mechanics reads: given where, in particular, Finally, we once again use the FEM to discretize in space the equation for the unknown γ f : given c h (t) and d h (t), find, for all t ∈ (0, T ), γ f,h (t) ∈ X h such that

Time discrezation
We naturally proceed with the time discretization of each spatially discretized core model, in the form of a non-linear system of ODEs. Following the notation defined in [19] we denote by z = (z w , z V , z γ f , z d ) T the block vector containing the unknowns of each single core problem and we write: where we define the differential operator This operator permits us to apply a first order time derivative to the ionic variables, the transmembrane potential and the fibres shortening, while a second order time derivative to the displacement. We use the BDF for the time approximation of Eq (4.2), in order to obtain a fully discretized formulation. Hence, we write:ż where ∆t = T N T is the timestep, N T being the number of timesteps, while the parameters ϑ (I) k , ϑ (II) k , k = 0, . . . , σ depend on the order σ of the BDF scheme. We will in particular use BDF of order σ = 2. For the electrophysiology we use a semi-implicit scheme. However, for the mechanical part, we proceed with an implicit scheme.
Eventually, in the implicit case, we obtain the following non-linear system: with z k assigned for k = 0, . . . , σ and we set for simplicityḋ 0,h = 0. Problem (4.3) is solved by exploiting the Newton method [37] at each timestep.

Fully segregated strategy
We proceed by solving a fully segregated scheme proposed in [20], instead of the previous monolithic approach used in [4] for the same electromechanical model of both ventricles.
In Figure 4 we graphically represent the time advancement. During each timestep t n , the algorithm performs in order, for m = 1, ..., N sub , the following steps: (1). find w (4.6) Figure 4. Graphical representation of the time advancement. Image taken from [20].
After N sub steps, we solve at t n+1 , with the Newton method, the implicit mechanics problem (M): We decided to use a staggered strategy since this has mainly two advantages that result in a faster computational solution of our problem. First, in a monolithic scheme, the time advancement must be carried out using the same timestep size ∆t for all the core models involved. Therefore, the small timestep size required by the electrophysiology is going to be used for the mechanics too, thus solving the latter more often than "necessary". Second, we will need both a large amount of memory since we have to store a highly dense Jacobian matrix and processes to assemble it. In order to overcome these issues, we follow a segregated algorithm exploiting the Godunov splitting scheme [21]. The latter consists of properly uncoupling the core problems solving the final system in two consecutive steps at each timestep.

Discretization of the cardiac cycle
We can now proceed with the description of the models used to reproduce the pressure and volume behaviour inside the ventricles. We will adapt the procedure presented only for the left ventricle in [19] to the case of both the ventricles. As we have seen rewriting the four phases for the bi-ventricle, we have to define a model which has to be able to take into account the fact that the ventricles are not aligned during the cardiac cycle, thus they can be in different phases modelled in a completely distinct way. We compute the volumes V endo,LV (t) and V endo,RV (t) of each chamber at time n by exploiting the formula which is rigorously derived in [39] and where ξ LV and ξ RV are vectors spanning the endocardium surface directed as the centerline of the LV and RV, respectively. We point out that the formulas 4.8 and 4.9 are valid only in the case in which the base is flat and orthogonal to the centerline. We then model the endocardial pressures p endo,LV (t) and p endo,RV (t) with different 0D models, proposed in [14,19,39,49], for each one of the phases presented above (in the following, we drop the "endo" superscript for simplicity and we leave "LV" or "RV" depending on the ventricle): (1). Both ventricles in isovolumic phase: We assume that both ventricles begin this phase at t = 0.
At the beginning the endocardial pressures are p LV = 10 mmHg and p RV = 5 mmHg. At each timestep, we update the pressure in the following manner: here we iterate on k until convergence, with p i,n+1 < ε is satisfied. The parameter ζ i < 0 has to be tuned, it has to be "sufficiently" large in order for the fixed point algorithm to converge but small "enough" to allow a quick convergence.The first ventricle to start the second phase is the right one, this happens when the pressure p RV,n+1 k+1 reaches the value p RV = 15 mmHg. (2). RV in ejection phase and LV in isovolumic: For the RV ejection phase we will use a two-element Windkessel model [51] of the form: p RV (T RV,1 ) = p RV (4.11) which is solved in the pressure variable, for simplicity, with a BDF scheme of order σ = 1 while the term dV RV dt is approximated at time n + 1 as V RV,n −V RV,n−1 ∆t , for simplicity. T RV,1 and T RV,2 are the initial and final time of this phase, for the right ventricle in this case, while the parameters C RV , R RV > 0 represent the capacitance and resistance of the equivalent electric circuit respectively. At the same time the left ventricular pressure keeps increasing till it reaches the value p LV = 85 mmHg.
(3). Both ventricles in ejection phase: Here both the pressures are updated through the Windkessel model: where i ∈ {LV, RV}. This phase ends when the (initially negative) term dV LV dt changes sign. (4). LV in isovolumic and RV in ejection phase: The second isovolumic phase for the left ventricle is modelled as the first isovolumic phase. Meanwhile, the RV ejection continues until the (initially negative) term dV RV dt changes sign. (5). Both ventricles in isovolumic phase: This phase is modelled as the first isovolumic phase and ends when the right ventricular pressure reaches the minimal value of p RV = 5 mmHg. (6). RV in filling and LV isovolumic phase: At each time step the right ventricular pressure is, for simplicity, linearly increased so that is reaches the EDP value at the final time T . The LV isovolumic phase finishes when p LV = 5 mmHg. (7). Both ventricles in filling phase: At each time step the pressures are linearly increased so that they reach the EDP value at the final time T .
We have here presented an adaptation of a model previously proposed for the left ventricle only [20], to the bi-ventricle case, in which the behaviour and the implementation of the various phases of the heartbeat are asynchronous.

Patient-specific mesh and fibres generation
Our geometry has been generated by image segmentation. We will take advantage of a complete human heart atlas presented in [24]. Since we are interested in modelling only the ventricles, we have extracted them from the whole heart, so obtaining the patient-specific mesh presented in Figure 5   To generate the fibres configuration we follow the Laplace-Dirichlet Rule-Based (LDRB) algorithm of [6], which is derived from histological and DTI data. The algorithm consists in extracting the transmural and apicobasal directions throughout the entire myocardium solving four Laplace equations with prescribed Dirichlet boundary conditions. The orthotropic fibers direction is obtained by a continuous interpolation throughout the myocardium passing through a straight-forward adaptation of quaternion spherical linear interpolation. Following the fibres orientation presented in [3], for a similar patient-specific bi-venticle geometry, we take different settings for α endo , α epi , β endo and β epi , these will result similar to the one found in human ventricles (see for examples the fibres obtained in [50]): By applying this method with the first set of angles to our geometry we obtain fibres and sheets displayed in Figure 6.

Numerical results
For our numerical simulations we employ LifeV * , an open-source finite element library. We apply the numerical methods implemented in a High Performance Computing framework and the computation are performed using Piz Daint, a Cray XC50/XC40 supercomputer installed at the Swiss National Supercomputing Center (CSCS) † . In Table 2 we present the common parameters used throughout the electromechanical simulations. Table 2. Parameters used in the electromechanical simulation: transversal and longitudinal conductivities ( mm 2 s ); transmurally heterogeneous wall thickening model parameters in Eq (2.2); activation model coefficients α (µM −2 ), c 0 , and µ A (µM 2 · s) in the four cardiac phases in Eq (2.2); density ρ ( g mm 3 ); Robin boundary condition coefficients ( kPa mm and kPa·s mm ) in Eqs (2.2)-(2.3); bulk modulus B (kPa), relaxation parameter for the two isochoric phases C I p and C II p ( kPa mm 3 ) in Eq (4.10); Windkessel model parameters C and R ( mm 3 kPa and kPa·s mm ) in Eq (4.12). We trigger the contraction of the ventricles simulating the propagation of the cardiac action potential along the myocardium. In physiological conditions the electrical wave follows the pathways prescribed by the fast conducting Purkinje fibres. To model this behaviour we decided to set delayed stimuli in five points in the bi-ventricle. The first stimulus is set at t = 0 ms at the centre of the septum, followed by an impulse next to the apex of each ventricle at t = 2.5 ms and one at the upper part of each free wall at t = 5.0 ms. In Figure 7 we can see the spreading of the transmembrane potential along the myocardium at different time steps. The electrical activation causes consequent changes in the Ca 2 + concentration resulting in mechanical stretching. We can notice how in Figure 8 the ion concentration and/or the fraction of open ionic channels on the cellular membrane follows clearly the transmembrane potential spreading along the cardiac tissue. The presence of calcium in the myocardium drives the mechanical activation and the consequent contraction and relaxation. In Figure 9 we present the time evolution of both endocardial pressure and ventricular internal volume together with the pressure-volume (pV) loops using different fibre setting. Furthermore, we decided to show the End Systolic Pressure-Volume Relationship (ESPVR), a feature representing the end-systolic elastance or inotropy, which provides a measure of the tissue contractility. This quantity is computed as the ratio between the ventricular pressure and volume at the end of the systolic phase. It can be easily evaluated as the slope of the tangent in the upper left corner of the pV loop. Even if ESPVR should be computed over different preloads and several heartbeats, we deem this biomarker to provide meaningful indications even in our simplified setting. In Figure 10 we find the pV loops together with the value of ESPVR of the left ventricle. It is worth to notice the effect of varying the fibre orientation on the slope of the ESPVR. We obtained the flattest slope with a fibre direction of 60 • , −60 • , this means that we have a decreased inotropy and so less contractility. All the other settings showed a steeper ESPVR, therefore varying the fibre orientation resulted in an increased contractility. Additionally, we wanted to see how an increasing bulk would affect the pressure and volume evolution along the heart beat, having the same fibre direction (in this case we chose 60 • , −60 • ). In Figure 11 we therefore report the pV loops and the bi-ventricle pressure/volume changes with three different bulk modulus (B = 1, B = 1.5 and B = 2). We highlight in Figure 12 the behaviour of the ESPVR with respect to the bulk modulus. As we expected, a lower bulk results in an increased contractility, since we are enforcing less our incompressibility constraint.
Other important quantities that are commonly exploited in clinical practice are stroke volume (SV) and ejection fraction (EF). The SV is computed by taking the difference between end-systolic volume (ESV) and end-diastolic volume (EDV) for a given ventricle. The EF is simply obtained dividing the SV with the EDV. In formula: In Tables 3 and 4 we report the values of SV, EF and ESPVR for the different fibre settings and the various bulk modulus. We notice how for LV the EF is increasing by varying the fibre direction from −55 • , 55 • to −70 • , 70 • . However, we obtained the highest value of EF for the right ventricle with a fibre orientation of −60 • , 60 • . From the simulations keeping fixed the fibres and changing the bulk modulus, we noticed an opposite behaviour happening in the two ventricles. Incrementing the bulk modulus resulted in a higher EF in LV and a smaller EF in RV.

Conclusions
In this work, we successfully implemented a segregated method for the electro-mechanical modelling of both human ventricles in which we exploited a transversal and isotropic model for the mechanical activation. We simulated one full heartbeat in a bi-ventricular geometry varying both active and passive mechanical parameters, such as fibre orientation and bulk, realizing pressure-volume loops.
We simulated the action potential propagation across the myocardium solving the monodomain equation coupled with the minimal ionic model. This leads to mechanical activation and contraction along the fibres following the action strain approach. We prestressed the geometry and solved a 0D circulatory models to take into account the blood-endocardium interaction.
We modelled a full heartbeat taking into consideration that the ventricles are asynchronous between each other, as a consequence each phase of the cardiac cycle is dependent on the ventricle. This complex joint behaviour is mutually influencing both ventricles through the septum. Future works should better address the filling stages to represent the full heartbeat in a detailed fashion including in the model both atria and the muscle contraction.
We performed our numerical simulations in a patient-specific geometry obtaining results qualitatively comparable with physiological data in healthy patients. Furthermore, we showed how changing fibres and bulk modulus influenced the pressure and volume evolutions over time of both ventricles in a different way.
We conclude that macroscopic mechanical indicators like pressure-volume loops, stroke volume, ejection fraction, etc. are hardly affected by fibers' orientation and values of bulk modulus (in a sufficiently wide range of values). Our findings somehow confirm the results of [8] regarding the pressure-volume loops of the left ventricle. We expect, however, that fibers direction and bulk modulus do have a large impact on local mechanical quantities in the tissue, such as stresses and strains.