An intergrid transfer operator using radial basis functions with application to cardiac electromechanics

In the framework of efﬁcient partitioned numerical schemes for simulating multiphysics PDE problems, we propose using intergrid transfer operators based on radial basis functions to accurately exchange information among different PDEs deﬁned in the same computational domain. Different (potentially non-nested) meshes can be used for the space discretization of the PDEs. The projection of the (primary) variables that are shared by the different PDEs (through the coupling terms) is carried out with Rescaled Localized Radial Basis Functions. We validate our approach by a numerical test for which we also show the scalability of the intergrid transfer operator in the framework of high performance computing. Then, we apply it to the electromechanical model for the human heart function, and simulate a heartbeat of an idealized left ventricle. We show that our approach enables the solution of large-scale multiphysics problems, especially when the individual models exhibit very different spatial scales.


Introduction
Intergrid transfer operators, acting among different grids defined over the same computational domain or accross different geometries, are employed in a variety of applications [15], whereas radial basis functions (RBF) are used in several fields in both computer science and applied mathematics, such as neural networks [31] and mesh free methods for solving PDEs [30,88,89]. The application of RBF as intergrid transfer operators is documented in [29] with the aim of exchanging information on a fluid-structure interface between non-conforming meshes defined on two different computational domains. Here we aim at constructing a new appropriate intergrid transfer operator based on RBF that has the capability to interpolate in a fast and accurate way both scalar and vectorial fields between different meshes defined on the same computational domain. We deal with systems of B Matteo Salvador matteo1.salvador@polimi.it 1 PDEs whose solution components represent different physical variables, and we want to use non-nested grids on the same computational domain to represent these different numerical variables. For the sake of illustration we will first introduce and verify our method on a simple elliptic system with two variables u 1 and u 2 that will be approximated on two different grids. Then we will address our target application in the framework of cardiac electromechanics. To the best of our knowledge, in this context, only techniques involving nested meshes have been used. Nested meshes are less flexible to accommodate geometrical heterogeneity and to be tuned to different local accuracy requirements. Besides, they cannot be generated in a completely independent fashion [14,19]. Moreover, proofs of scalability are not provided for this kind of techniques; also, comparisons in terms of computational costs are lacking.
Our focus in cardiac electromechanics moves from the need to better understand cardiac function, especially in pathological cases. Cardiovascular diseases are indeed one of the most common causes of death globally [48], and several pathologies are still not completely understood. Mathematical modeling of the heart and numerical simulations allow a better understanding of the phenomena occurring both in physiological and pathological conditions [36,50,59,85,86]. Several types of processes occur in the human heart, such as propagation of an action potential in the myocardium, which contracts, together with the blood flowing in the four chambers (atria and ventricles) and through the valves [21,27,56,78,80,81]. This problem is challenging from the numerical standpoint [13,17,28,42] as it involves different temporal and spatial scales. Cardiac electrophysiology, blood fluid dynamics and myocardial mechanics require a different numerical resolution in space and time, going from a really detailed one for what concerns the electric part, moving towards a bigger one for the mechanical behaviour [60,61]. For this reason the application of multi-mesh and staggered methods is strongly justified from a physical viewpoint, and permits to solve in a really efficient way this multiphysics problem without loosing the details of the phenomena that we want to model.
In this work we focus on the electromechanical modeling of the left ventricle, which has been extensively studied over the past years [14,20,32], but still needs further investigations both from theoretical and numerical perspectives. For the electric part we consider the monodomain and Bueno-Orovio models [9,58,73]. For the mechanical part we use the Holzapfel-Ogden model [40] together with the active strain formulation [2,3], combined with a model for the transmural heterogeneous thickening of the myocardium [6]. Myocardial fibers contraction defines the bridge between electrophysiology and mechanics [71]. From the numerical viewpoint, the Finite Element Method (FEM) is used in order to discretize in space the continuous single core models (i.e. electrophysiology, activation and mechanics) by means of piecewise linear elements, whereas time discretization is carried on using Backward Differentiation Formulas (BDFs) of order 1 [62]. As already done in literature [5,14,45], the integration of the discrete core models leads to the formulation of segregated (i.e. which are all solved separately and sequentially) and staggered (i.e. using with different timesteps) algebraic equations. The timestep changes according to the time scale of the single problem. Semi-implicit or implicit time discretizations are considered depending on the stiffness of the models. Moreover, according to the required space resolution, we use different independent (i.e. non-nested) meshes and we perform intergrid fields transfer by means of RL-RBF. We use different simplified Cauchy models to describe the blood flow inside the left ventricle. In this way we can compute how pressure and volume evolve along the heartbeat [26,69,86]. We also consider a prestress technique, which is applied in the pre-processing phase, to estimate the internal stresses of the myocardium at the initial time of the simulation [41,82].
As observed before, electrophysiology and mechanics are solved separately, in a segregated (partitioned) fashion. We use different meshes and timesteps (staggered approach) for the two fields. After space and time discretization we solve each block linear system by means of generalized mini-mal residual (GMRES) method. We use our partitioned and staggered electromechanics solver to carry out a numerical simulation in the High Performance Computing framework of the whole cardiac cycle for an idealized left ventricle geometry, and we analyze the numerical results in terms of clinically relevant indicators: specifically, we produce the so-called pressure-volume loops in order to assess the left ventricle function and its properties.
The paper is organized as follows: in Sect. 2 we present the interpolant based on RBF that will act as intergrid transfer operator. In Sect. 3 we propose a test case with known exact solution to show the accuracy and the reliability of the operator introduced in Sect. 2. In Sect. 4 we apply the methodology that we have developed and tested to cardiac electromechanics, a complex framework in which different physics (electrophysiology, mechanical activation and mechanical deformation of the myocardium) and different space scales are present. We also propose a novel partitioned scheme for the time discretization of the single core models related to cardiac electromechanics, for the numerical simulation of one heartbeat of an idealized left ventricle. We finally draw our conclusions in Sect. 5.

Intergrid transfer operator
Our aim is to transfer the values of a certain function from one mesh to another one. However we present in this section our intergrid transfer operator in the most general framework. We exploit several properties of RBF to perform the interpolation task considering a general function f , following the idea developed in [22].
Let f : R 3 → R d be a scalar (d = 1) or possibly vector field (d > 1) in the 3D setting. Given a set of M nodes of the general field f by means of RBF in the following way: with {γ f m } M m=1 set of the interpolation weights. RBF are denoted by π(·, r ), which can be either globally or locally supported according to the choice of the radius r .
Other options, involving also globally supported basis functions, are available as well [10,25]. We introduce an interpolation matrix Φ int ∈ R M×M such that Φ int i, j = π(||ξ i − ξ j ||, r ) with i, j = 1, . . . , M. We call f Ξ the evaluation of the field f in all the M interpolation nodes that belong to the set Ξ .
The interpolation constraint is algebraically expressed as follows: with γ f = {γ f m } M m=1 solution of linear system (2). Both fields f and Π f (x) have the same value at the interpolation nodes, i.e.
The choice of local RBF leads to a sparse pattern of the matrix Φ int . At this point, once Π f (x) is determined, we can evaluate the interpolant on a set Λ = {λ n } N n=1 of N different points with respect to the interpolation nodes contained in Ξ : In our application Ξ and Λ will be two different sets of nodes of two independent triangulations of the computational This sparse matrix is used to determine f Λ , i.e. the evaluation of the RBF interpolant Π f on Λ: In order to obtain a smoother interpolant that is able to interpolate exactly any constant field and that is accurate for small values of the radius r [22], we rescale Π f (x) by the interpolant Π g (x) of the constant function g(x) = 1, which assumes a value equal to one at each interpolation point: We formulate in this way Rescaled Localized Radial Basis Functions (RL-RBF). From the algebraic perspective, the interpolation problem associated with (5) can be written in the following form: where γ g = {γ g m } M m=1 and 1 Ξ vector of ones on the interpolation nodes defined in Ξ . Linear systems (6) and (7) are solved separately. The evaluation of interpolantΠ f in a specific point x is: where

Numerical test for elliptic PDEs
We want to illustrate the properties of our intergrid transfer operator in terms of scalability and convergence for both the L 2 and H 1 norms in a given domain Ω ⊂ R d . With this aim, we propose the following test case: with Ω = (−1, 1) 3 and v = 1. The forcing terms are: sin(π y)sin(π z) + π cos(π x)sin(π y)sin(π z) + π sin(π x)cos(π y)sin(π z) + π sin(π x)sin(π y)cos(π z).
This one-way 2-field coupled system involving Laplace and diffusion-advection-reaction PDEs is well posed and is endowed with an exact solution, u 1 =sin(π x)sin(π y)sin(π z) and u 2 = −sin(π x)sin(π y)sin(π z). We use the Finite Element method to solve numerically this test case. We provide two different meshes T h 1 and T h 2 of the computational domain made by tetrahedrons, with h 1 and h 2 representing the maximum size of the element K in terms of maximum mean diameter of the circumscribed circumference, with ∪ K ∈T h 1 K = ∪ K ∈T h 2 K = Ω. Both structured and unstructured meshes can be potentially employed, either in a nested or a non-nested fashion. We denote N u 1 and N u 2 the number of degrees of freedom for u 1 and u 2 respectively (these are the internal finite element nodes). We also introduce two finite dimensional spaces X r respectively. The discretized formulation of (9) reads: find where u 1,h 2 = N u 2 j=1 u 1, j ψ j is the interpolation of u 1,h 1 on the second mesh, and again v = 1. For this purpose, we use the RL-RBF introduced in Sect. 2, so that, at a continuous level, u 1,h 2 (x) =Π u 1 (x). In order to take into account the distribution of mesh points of T h 1 and T h 2 , we define an adaptive strategy to select the radius of the support of Beckert & Wendland basis functions by means of the number of links that a certain vertex of the mesh has with the surrounding neighborhood. In this way we exploit the structure of the mesh to build a variable and local support of the RBF that keeps into account the level of refinement of the mesh in each region of the computational domain. For more details about this technique we refer to [22]. A number of links equal to 1, as shown in Fig. 1, is sufficient to obtain a good interpolated solution. The derivatives of u 1,h 2 with respect to x, y and z are calculated using the Zienkiewicz-Zhu gradient recovery method, which is known to be efficient and superconvergent [92,93].
We have performed strong scalability tests and convergence analysis of both L 2 and H 1 errors for u 1,h 1 and u 2,h 2 , considering structured, nested unstructured and non-nested unstructured grids. P 1 finite elements are used in all sim- ulations. In Fig. 2 we show a comparison over the entire computational domain and over a cross-section of u 1,h 1 versus u 1 and u 2,h 2 versus u 2 , considering, as depicted in Fig. 3, two different unstructured grids that are non-nested and providing in both cases a good match of the numerical and analytical solutions. We see that our approach based on a RBF interpolant whose support is chosen by the number of links is able to pass properly the information of u 1,h 1 inside the equation of u 2,h 2 . This is confirmed by the fact that the exact solutions u 1 and u 2 match really well the corresponding numerical solutions u 1,h 1 and u 2,h 2 . We also test the scalability of our RL-RBF operator, as shown in Fig. 4, where we see that the time that we need to interpolate u 1,h 1 inside the equation for u 2,h 2 decreases with the number of CPUs with all types of grids. We want to underline that this behaviour is observed in the case the coarse mesh is used for u 1,h 1 and the fine one for u 2,h 2 , and viceversa. Moreover, due to the generality that we want to obtain with our intergrid operator, it seems that, independently of the value of the mesh size, the interpolation process scales better with non-nested unstructured grids.
In Fig. 5 we show a convergence analysis of L 2 and H 1 errors where we consider a fixed grid for u 2,h 2 , with a certain mean diameter of the elements h u 2 = h 2,mean , and we gradually reduce the h u 1 = h 1,mean value, and viceversa. Table 1 reports the results in case we make a progressive increment of the number of mesh elements for both u 1,h 1 and u 2,h 2 in parallel: we start from a certain couple of nonnested unstructured/structured meshes with specific values of h 1,mean and h 2,mean , and then we halve these values and analyse the behaviour of both L 2 and H 1 errors. We highlight that the order of convergence that we expect theoretically (i.e 1 for the H 1 norm error and 2 for the L 2 norm error) is always recovered for every mesh used.

Cardiac electromechanics
Now we apply the approach developed and tested in Sects. 2 and 3 to solve the cardiac electromechanics problem.

Mathematical model
Heart tissue is made of cardiomyocytes, which determines the orthotropic structure of the left ventricle through the organization in fibers and laminar collagen sheets [76]. In this section, we review the model governing the electromechanical behaviour of the left ventricle by accounting the characteristics of myocardium [26].

Electrophysiology
The contraction of a single cardiac cell is initiated by an electric activation due to an action potential, a depolarizing phase that raises the so called transmembrane potential of an excitable cell from its resting value ranging between -90 and -80 mV to slightly positive values, followed by a peak, a plateau, and a repolarization, that returns the transmembrane potential to its resting value [17,18]. The electric activity of the heart starts at the so called sinoatrial (SA) node and propagates throughout the right atrium [18]. Thanks to the Bachmann's bundle and some other preferential lines of trasmission, the electric signal reaches the left atrium [17]. Then, the activation front arrives at the atrioventricular (AV) node, which conducts the action potential through the nonexcitable atrioventricular septum, activating the special-ized fibers of the bundle of His and the Purkinje network, that spreads as a tree-like left and right bundle branches ending on the endocardial surface of the ventricles. These Purkinje terminations transmit the action potentials to the ventricular walls and cardiac excitation then propagates throughout the ventricles [46,75].
The monodomain model is a diffusion-reaction PDE system able to describe the electric properties of cardiac muscle cells, assuming the same anisotropy ratios between the intracellular and extracellular spaces [58,73]. It is a continuum model, which means that it is used to capture average properties of many cardiomyocytes, and not the behaviour of single cells. The monodomain model reads: Ω 0 ⊂ R 3 is the domain in the reference configuration, represented in our case by an idealized left ventricle at the end of the diastolic phase. T > 0 is the final time. C m is the total membrane capacitance and χ is the area of cell membrane per tissue volume. u is the dimensionless transmembrane potential, whereas the vector z = {z 1 , z 2 , . . . , z k } expresses k recovery (or gating) variables, which play the role of probability density functions and model the fraction of open ionic channels across the membrane of a single cell. D M = σ t I + (σ l − σ t ) f 0 ⊗ f 0 refers to the diffusion tensor, being f 0 the vector field expressing the fibers direction and σ l , σ t ∈ R + the longitudinal and transversal conductivities respectively [72]. By defining X and x as the reference and deformed coordinates, we introduce the , keeping into account the effect of the mechanical displacement d s on the electrophysiology. I app (t) is an external applied current, which should simulate in our case the behaviour of the Purkinje network. Indeed we use it to trigger the action potential at specific points of the myocardium. I ion (u, z) is the feedback from the cellular scale into the tissue one, and strictly depends on the chosen ionic model. A Neumann boundary condition is applied all over the boundary and defines the condition of electrically isolated domain.
In literature there are several possible choices of ionic models [1,47,84]. We present here the monodomain equation supplemented by the Bueno-Orovio minimal ionic model, which is specifically designed for the human left ventricle to describe from a phenomenological perspective the microscopical details of each single cardiomyocyte [9]: where is an Heaviside function centered at x 0 , z = {v, w, s} is the vector of gating variables and evolves over the time thanks to a system of 3 ODEs, which are solved at each point of the reference domain. The ionic current I ion (u, z) is given by: The function I ion (u, z) and all the parameters related to the Bueno-Orovio ionic model are taken from [9], by considering values that permit reproducing the behaviour of the ten Tusscher-Noble-Noble-Panfilov (TNNP) model [83].

Mechanical activation
The mechanical activation bridges electrophysiology and passive mechanics. We consider a phenomenological model that keeps into account the local shortening of the fibers γ f at the macroscopical level [33,69,71]. In this way we are able to describe the behaviour of the fibers in a faster but more approximated way, with respect to more complex and accurate models that are able to describe the dynamic of sarcomeres [65]. Myocardial displacement feedback d s and concentration of calcium ions Ca 2+ play an important role in the time evolution of γ f . As already done in literature [33,68,69], the gating variable s substitutes the concentration of Ca 2+ , due to the fact that they have a pointwise similar time pattern, even if the order of magnitude is different. The formulation is the following: where ] is a truncated Fourier series expressing the sarcomere force-length relationship [37]. Both α and μ A should be calibrated according to the specific case under investigation. The active deformation is calculated exploiting the following orthotropic form [3]: where s 0 and n 0 represent the sheets and their normal direction respectively, with γ s and γ n corresponding local shortening or elongation [6,52]: Here λ represents a transmural coordinate, varying from λ endo at the endocardium to λ epi at the epicardium, which permits to have a transversely non-homogeneous thickening of the left ventricular wall. γ s set like (16) yield to det(F A ) = 1 (Fig. 6).

Active and passive mechanics
In order to describe the displacement of the myocardium along the entire heartbeat, we model the properties of the tissue by means of fibers f 0 , sheets s 0 and their normals n 0 , which permit to obtain highly anisotropic internal stresses associated with a prescribed deformation [39]. We use a nearly-incompressible formulation by weakly penalizing large volumetric variations [77]. Moreover, the dimensionless variable γ f provides a link between electrophysiology and mechanics. The momentum conservation equation with boundary and initial conditions reads [51]: A Robin boundary condition is prescribed at the epicardium to account for the action of the pericardial sac, so that the presence of the pericardium is also addressed and modelled [57]. K epi ⊥ , K epi , C epi ⊥ , C epi ∈ R + are respectively local values of stiffness and viscosity constants of the epicardial tissue in the normal or tangential directions. A Neumann boundary condition is defined at the base and acts as a free contraction condition. p endo (t) is the internal pressure of the ventricular chamber and is modelled through a 0D model, as will be detailed in Sect. 4.1.5. The boundary condition at the endocardium should depend on the deformation field [64]: we however neglect this feedback in this work as we deem this to be an acceptable approximation aimed at reducing the overall complexity of our model. The Piola-Kirchhoff tensor P = P(d s , γ f ) incorporates both the passive and active mechanical properties of the tissue. After defining the symmetric positive definite right Cauchy-Green tensor C = F T F, with F = I + ∇d s deformation tensor, the strain energy function W : R 3×3 − → R provides a link between the strain and the energy of the material. Under the hyperelasticity assumption, the strain energy function can be differentiated with respect to the deformation tensor F in order to obtain P: There are several models in literature which are able to describe the anisotropic nature of the tissue, such as the Guccione [38] or the Holzapfel-Ogden laws [40]. We use the second one, which formulates the following additive decomposition of the strain energy function: W such that J = 1 is its global minimum. In this way, we are able to set the nearly-incompressible constraint [16,23,91]. B ∈ R + is the bulk modulus, which weakly enforced the incompressibility constraint [33]. Following [74], we model the evaluation of the isotropic term W 1 in J − 2 3 I 1 and not directly in the first invariant (Fig. 7).
The active part, which is due to the electrical impulse, is formulated in the active strain framework [2,3,49,69]. In addition to the reference configuration Ω 0 and the deformed one Ω, we introduce an intermediate stateΩ, which represents the active part of the deformation. By means of the 2 nd order tensor F A , we map Ω 0 intoΩ, whereas the F E tensor takes the role to finally transformΩ into Ω. We finally reach, the multiplicative decomposition of F = F E F A . The first Piola-Kirchhoff strain tensor P reads: For additional details on the final form of tensor P, we refer the reader to [33]. he inertial term in (17) allows us to properly account for the evolution in time of the mechanical problem. Moreover, the division by a small timestep of the mass matrix yields a strong zero-th order diagonal contribution in this non-linear system, easing the numerical resolution of this non-linear problem. Nevertheless, the term depending on the Piola-Kirchhoff tensor is dimensionally predominant.

The multifield coupled problem
The fully coupled multifield electromechanical problem is presented here below: where p endo (t) is the unknown of a 0D fluid problem that changes over the time with the different phases that involve a heartbeat.

Cardiac cycle
We model the evolution of both endocardial pressure and volume along the heartbeat, considering a total duration T = 0.8 s. With this aim, we consider ordinary differential equations to model the fluid in the left ventricle chamber, where we assume that the endocardial pressure is uniform over the domain. The heartbeat is conventionally split into four phases, in which we solve four different ODEs [33]: -Isovolumetric contraction: increase of p endo from the End Diastolic Pressure (EDP, about 10 mmHg) to the aortic pressure (about 85 mmHg) while the volume V endo remains almost constant, according to the following equation: where V endo is set to the initial volume, and T 1 = T 1 ( p endo ) is the earliest occurrence in time at which p endo ≥ p ao , that forces the aortic valve opening. -Ejection: the left ventricle contracts and pushes blood to flow through the aortic valve, which closes itself at the end of this phase (when p endo ≤ p ao ). V endo decreases following a two 0D element Windkessel model [90]: with p endo (T 1 ) = p ao , C, R > 0 are two parameters representing the capacitance and the resistance of the electric circuit that mimics the blood flowing in the aorta. T 2 = T 2 (V endo ) is the first instant of time in which dV endo dt ≥ 0.
-Isovolumetric relaxation: decrease of p endo with constant value of V endo , treated from the mathematical perspective in the same way as the isovolumetric contraction. T 3 = T 3 ( p endo ) is the first occurrence in time at which p endo ≤ p endo min = 5 mmHg. -Filling: the mitral valve opens due to the pressure drop, blood flows inside the left ventricle causing an increment of V endo until both p endo and V endo attain the EDP values. We increase linearly p endo until it reaches the initial value p endo E D P atT 3 = 0.7 s and we keep it constant from T 3 to the final time T = 0.8 s by setting: where The endocardial volume V endo (t) is computed exploiting the following formula: which is derived in [67], with ξ vector directed in the normal direction of the left ventricle base.

Numerical discretization
In this section we numerically discretize the problem (21) using the Finite Element Method in space and Backward Differentiation Formula in time. Our numerical method permits to separate and manage properly the space and time scales related to cardiac electromechanics. We try to reduce at maximum the number of interpolations of the fields defined on different grids and the number of time advancements of the different physics while preserving accuracy and stability.

Space discretization
We provide two meshes T h 1 and T h 2 of the computational domain made of tetrahedrons. h 1 and h 2 (with h 1 < h 2 ) represent the maximum size of the element K , say the maximum of the mean diameter of the circumscribed circumference.
T h 1 , which is the fine mesh, is used for electrophysiology. The coarser T h 2 mesh is employed for both activation and mechanics. This is motivated by the fact that we need a higher resolution for the electric part, going down to the cellular level, whereas both cardiac mechanics and activation evolve on larger space scales. The geometry for the left ventricle is represented by a prolate ellipsoid (as often done in literature [26,33,38,67]). We denote by N u , N z , N γ f and N d s the number of degrees of freedom for dimensionless transmembrane potential, gating variables, mechanical activation and displacement respectively. We define the set of polynomials with degree smaller than or equal to r over an element K of the mesh with P r (K ), and we introduce the finite dimensional spaces X r The semidiscretized formulation of the Monodomain equation reads: find and F h 1 are the semidiscretized versions of the gating variables and the interpolated deformation tensor respectively, whereas u h 1 (t) = N u j=1 u j,h 1 (t)φ j is the finite element solution that approximates u(t).
F h 1 is obtained with the following procedure: -RL-RBF are employed for the interpolation of d s,h 2 , which is obtained from (17). The interpolantΠ d s (x) is built on T h 2 and it is used to obtain d s,h 1 on T h 1 following the procedure explained in Sect. 2 and exploiting formula (5). -We get F h 1 = I + ∇d s,h 1 through the adaptation of the Zienkiewicz-Zhu gradient recovery method [92,93] to the tensor case.
We can rewrite equation (27) as a system of non linear ODEs by setting u h 1 (t) = {u j,h 1 (t)} N u j=1 : where · ∇φ i dΩ 0 and: In order to avoid numerical instabilities, we use a lumped M L approximation of the mass matrix M [11]. There are several strategies in literature for the treatment of the non linear term I ion (u h 1 (t), z h 1 (t)) [33,43,[53][54][55]. We use the so-called state variable interpolation (SVI) approach. By denoting with {x K s } N q s=1 and {ω K s } N q s=1 the quadrature nodes and weights of a generic element of the mesh K ∈ T h 1 , both u h 1 and z h 1 are evaluated at the quadrature nodes as follows: Bueno-Orovio ionic model The ionic model is a system of ODEs which indirectly depends on the space variable through the transmembrane potential u. The vector of gating variables, which is a state vector located at the nodes of the mesh, is rearranged as follows: The semi-discrete formulation can be written as follows: Activation model The Galerkin formulation related to the equation for γ f reads: given d s,h 2 (t) and s h 2 (t) =Π s h 1 (x) interpolation of s h 1 (t) by means of RL-RBF, i.e. formula (5), By introducing the proper matrices, the following system of ODEs is obtained: its basis. The semi-discretized version of (17) reads: given with d s,h 2 The algebraic formulation reads: with: Equations (27), (31), (32) and (34) provide the semidiscretization of (21) in a splitted fashion ready for a partitioned and staggered time discretization.

Time discretization
In order to fully discretize the electromechanical problem, we introduce a block vector y = {y z , y u , y γ f , y d s } containing all the unknowns of the problem and we reformulate the semidiscrete problem as follows: where T and H terms represent the core models. The time discretization performed with the BDF scheme of general order σ [62] reads as follows: where Δt = T N T is the timestep, N T is the total number of timesteps, and θ I k , θ I I k , k = 0, . . . , σ/σ + 1 depend on the order of the BDF scheme. We use σ = 1.

Segregated scheme
Implicit monolithic strategies are stable and lead to accurate results, as shown in [33]. However there are two types of drawbacks connected with this approach. One is we are forced to use the same timestep for both electrical and mechanical parts even if the timescale of electrophysiological phenomena is much smaller than the one of the myocardial activation/mechanics. Moreover, the calculation of the Jacobian matrix J E M is quite demanding and this matrix requires a significant amount of RAM. In order to overcome these issues, we propose a segregated strategy based on the Godunov splitting scheme [35]. This approach permits to advance in time faster and consuming less memory, at the expense of introducing first order error on the solution [32,34]. This splitting error, along with the first order accuracy of the BDF scheme with σ = 1, degrades the overall convergence of the numerical solution in time by another term of order one. Due to the fact that the time integration scheme does not guarantee unconditional stability in general [62], we have to impose a limitation on the timestep size Δt. For both electrophysiology and activation we employ, as done in [32], a semi-implicit scheme. Mechanics is instead numerically discretized in time implicitly, due to the fact that the highly non-linear (exponential) terms of the strain energy function W would need a restrictive Δt in both the semi-implicit and explicit contexts. We refer from now on to (I SI ) − (E SI ) − (A SI ) for the fully-segregatedsemi-implicit scheme applied to the ionic, monodomain and activation models, and to M I for the implicit scheme applied to the mechanical core. We employ a staggered approach in which two different timesteps are used for (I SI ) − (E SI ) and (A SI ) − (M I ). We indicate Δt as the timestep for both activation and mechanics. τ = Δt N sub is the timestep for electrophysiology (ionic and monodomain models), with N sub ∈ N number of intermediate substeps that must be done by (I SI ) − (E SI ) before a timestep Δt of (A SI ) − (M I ) is performed. The time advancement that has been just described is sketched in Fig. 8.
We set t for m = 1, . . . , N sub , whereū I h 1 ,z I h 1 are evaluated by using the variables on the fine mesh T h 1 at times t n , t n −τ, . . . , t n − (σ − 1)τ , with σ = 1 in our case. After solving (38) and (39) for N sub steps, treat (A SI ) − (M I ) at t n+1 in the following way: - -d n+1 s,h 2 is obtained by solving: by means of the Newton method [62].  [12].
Finally we solve the ordinary differential equations representing the fluid in the left ventricle chamber, again in a segregated fashion. The endocardial volume V endo,n+1 h is

Numerical results
In this subsection we present a numerical simulation of a full heartbeat lasting 0.8 s in the electromechanical framework, by considering an idealized left ventricle. We use LifeV [44], an open-source finite element library, for the resolution of the electromechanical problem in a High Performance Computing framework and for the implementation of the numerical method. All the computations were carried out using a full node (32 Intel ® Xeon ® E5-4610 v2 2.3GHz cores) of the HPC centre available at MOX.
In Figs. 9 and 10 we show the meshes that we have used. They have been generated using VMTK -the Vascular Modelling Toolkit [4,79], and their information is reported in Table 2. With VMTK we can choose a certain diameter h mean and we can generate fine enough meshes featuring elements with similar aspect ratio. This toolkit is also used for the tagging procedure to determine which elements belong to the base, the endocardium and the epicardium. Different tags are assigned to these three different regions of the computational domain to apply in a proper way boundary conditions. With our methodology we are also able to deal with different grids that are completely independent. We use P 1 finite elements in order to approximate monodomain, activation and mechanics equations, so that the unknowns are actually calculated on the vertices of each tetrahedron. The ionic model evolves directly on the vertices of each element. We employ a first order BDF scheme in time [62]. The timestep for electrophysiology is Δt = 50 µs, whereas the one for activation and mechanics is Δt = 250 µs. We also have to solve several linear systems coming from:  We perform these tasks using the GMRES method with a stopping criterion based on the relative residual and a tolerance given by 10 −8 . We apply a current I app (term of the monodomain model) for 2 ms in three different points of the myocardium to trigger the electrical signal in the left ventricle. Even if we do not keep into account precisely the propagation of stimuli inside the Purkinje network [66,87], the electrical activation of the tissue that we propose is known to provide physically acceptable results [32] (Table 3).
In literature there are several rule-based techniques to generate fibers and sheets distribution [7,24,69] for both idealized and patient-specific cases. We use the strategy proposed in [33], where fibers rotate transmurally from the epicardium, with α epi = −60 • to the endocardium, with α endo = +60 • . In order to properly take into account the internal stresses of the myocardium at the beginning of the electromechanical simulation, we have to apply the so called prestress to our computational domain. We compute a distribution of stresses such that the reference geometry is in equilibrium with the blood pressure p endo at the end of the diastolic phase. An additive decomposition of tensor P = P(d s ) + P 0 is operated, where the prestress tensor P 0 is determine to ensure a null displacement d s,0 in correspondance of the initial pressure at the endocardium. For more information about this technique we refer to [32][33][34].
In Figs. 11 and 12 we observe the evolution in space and time of the transmembrane potential V , activation variable γ f and total displacement magnitude |d s | over one entire heartbeat. The conduction velocity is overestimated due to the fact that, even if the electrophysiogical mesh is fine, we should use a smaller value of h mean to describe properly all the space scales and to have a convergent velocity of the wavefront [8,33,67]. The activation is slightly delayed with respect to the propagation of the action potential because it is driven by the calcium concentration (here approximated with the gating variable s), that evolves in time more slowly than the transmembrane potential. The myocardial tissue undergoes a significant thickening, which is in accordance with experimental observations [63]. A high value of the bulk modulus B permits to obtain a significant torsion of the left ventricle and to impose the quasi-incompressibility constraint [33]. With the choice of parameters for the Robin boundary condition at the epicardium in the mechanics problem (17), we are able to properly keep into account the effect of the pericardium [57]: in this way we can reduce the movement of the apex and increase the contraction of the base.
The pV-loop is reported in Fig. 13. Even if a comparison with in-vivo measurements would be meaningless, due to the fact that we are dealing with an idealized framework, we can say that the pV-loop developed over the simulation is in accordance with those observed experimentally [70], at least qualitatively.
Finally, with reference to [33], we have also performed this simulation in a monolithic fashion, considering only the first isochoric phase, for a total time T = 50 ms. In this case we are forced to use the timestep of electrophysiology (here Δt = 50 µs) also for both activation and mechanics, due to the fact that action potential and calcium dynamics need a higher resolution in time. For what concerns space discretization, we are again forced to use only the mesh of electrophysiology (here the one of the first row in Table  2), even if we do not need such a high number of ele-  Fig. 15 Comparison between monolithic and segregated-staggeredintergrid approaches for γ f over the computational domain at t = 0.03 s ments for the mechanical problem. In Table 4, we report the time needed to perform one single timestep with both monolithic and segregated-staggered-intergrid approaches. We place ourselves in the worst scenario where we consider a timestep for the segregated-staggered-intergrid scheme in which we have to solve both the mechanical problem and interpolate the calcium concentration. This does not occur in all the times of the simulation, as explained in Sect. 4.2.3. We see that the use of two different timesteps and two different meshes (allowed by our intergrid transfer operator), yields a significant gain with respect to the monolithic strategy.
Monolithic schemes are indeed poorer in terms of performance with respect to segregated ones. With the approach presented in this paper we observe a 10x speed-up in the final computational time with respect to [33], while guaranteeing a proper and accurate capture of both the time and space scales of this multiphysics problem. Indeed, in Figs. 14 and 15, we observe small differences in the velocity of propagation of the wavefronts. This little discrepancy is due to the splitting error of the segregated scheme, to the small error in space introduced by the interpolation process, and finally to the larger timestep used for activation and mechanics. Same conclusions hold also for the electrophysiology problem, due to the feedback that involves the deformation tensor F.

Conclusions
In this work, we proposed a novel segregated solver, which makes use of an accurate and efficient intergrid transfer operator based on RBF to simulate the electromechanical activity of an idealized left ventricle. We consider a coupling between the monodomain equation and the Bueno-Orovio minimal ionic model for the electric part. Cardiac mechanics is modelled in the active strain formulation using the Holzapfel-Ogden constitutive law and a transmurally variable activation. Electrophysiology and mechanics are linked by means of a phenomenological model that reproduces fibers contraction. Prestress and additional 0D problems for the fluid permit to keep into account the interaction of the myocardium with the blood inside the left ventricle. We solve this multifield coupled problem in the High Performance Computing framework. The use of two different meshes is motivated by the fact that the level of resolution required by the two physics is not the same, going from an element size with the order of 0.1 mm for electrophysiology, to 1 mm for mechanics. This allows to obtain a significant speed-up in the numerical simulation with respect to the case in which we use one single mesh that could be forced to adopt the same resolution of the electrophysiological model. The same considerations lead to the choice of different timesteps and to solve the electromechanical problem in a staggered fashion, which is again faster and less memory-demanding than the monolithic approach.