
Citation: Diogo Gomes, Julian Gutierrez, Ricardo Ribeiro. A mean field game price model with noise[J]. Mathematics in Engineering, 2021, 3(4): 1-14. doi: 10.3934/mine.2021028
[1] | Christian Cortés García, Jasmidt Vera Cuenca . Impact of alternative food on predator diet in a Leslie-Gower model with prey refuge and Holling Ⅱ functional response. Mathematical Biosciences and Engineering, 2023, 20(8): 13681-13703. doi: 10.3934/mbe.2023610 |
[2] | Claudio Arancibia–Ibarra, José Flores . Modelling and analysis of a modified May-Holling-Tanner predator-prey model with Allee effect in the prey and an alternative food source for the predator. Mathematical Biosciences and Engineering, 2020, 17(6): 8052-8073. doi: 10.3934/mbe.2020408 |
[3] | Kunlun Huang, Xintian Jia, Cuiping Li . Analysis of modified Holling-Tanner model with strong Allee effect. Mathematical Biosciences and Engineering, 2023, 20(8): 15524-15543. doi: 10.3934/mbe.2023693 |
[4] | Peter A. Braza . A dominant predator, a predator, and a prey. Mathematical Biosciences and Engineering, 2008, 5(1): 61-73. doi: 10.3934/mbe.2008.5.61 |
[5] | Gunog Seo, Mark Kot . The dynamics of a simple Laissez-Faire model with two predators. Mathematical Biosciences and Engineering, 2009, 6(1): 145-172. doi: 10.3934/mbe.2009.6.145 |
[6] | Shuangte Wang, Hengguo Yu . Stability and bifurcation analysis of the Bazykin's predator-prey ecosystem with Holling type Ⅱ functional response. Mathematical Biosciences and Engineering, 2021, 18(6): 7877-7918. doi: 10.3934/mbe.2021391 |
[7] | Eduardo González-Olivares, Betsabé González-Yañez, Jaime Mena-Lorca, José D. Flores . Uniqueness of limit cycles and multiple attractors in a Gause-typepredator-prey model with nonmonotonic functional response and Allee effecton prey. Mathematical Biosciences and Engineering, 2013, 10(2): 345-367. doi: 10.3934/mbe.2013.10.345 |
[8] | Christian Cortés García . Bifurcations in a discontinuous Leslie-Gower model with harvesting and alternative food for predators and constant prey refuge at low density. Mathematical Biosciences and Engineering, 2022, 19(12): 14029-14055. doi: 10.3934/mbe.2022653 |
[9] | Saheb Pal, Nikhil Pal, Sudip Samanta, Joydev Chattopadhyay . Fear effect in prey and hunting cooperation among predators in a Leslie-Gower model. Mathematical Biosciences and Engineering, 2019, 16(5): 5146-5179. doi: 10.3934/mbe.2019258 |
[10] | Jinxing Zhao, Yuanfu Shao . Bifurcations of a prey-predator system with fear, refuge and additional food. Mathematical Biosciences and Engineering, 2023, 20(2): 3700-3720. doi: 10.3934/mbe.2023173 |
The heart plays the crucial role of pumping blood into the circulatory system by providing all sort of vital substances to the cells. Moreover, cardiovascular related diseases represent the leading causes of death in the whole world [60]. While advancements in medical practice are continuously improving patients conditions and diseases outcomes, recent progresses in mathematical modeling and computational mechanics allow to perform realistic cardiovascular numerical simulations [38,65,76,77,84,97,105,107], thus potentially providing medical doctors and clinicians with valuable diagnostic and predictive tools. Moreover, clinical data and the application of image segmentation techniques to Magnetic Resonance Imaging (MRI) and Computed Tomography (CT) scans feed, as inputs, the mathematical models, thus allowing numerical simulations in subject- and patient-specific frameworks [23,55]; in addition, uncertainty quantification techniques allow to estimate the models parameters and data, and to cope with their variability [30,59].
The mathematical modeling of the heart involves several challenges intrinsically related to the complexity of its function [17,20,33]. A satisfactory model must be able to describe a wide range of different processes, such as the evolution of the transmembrane potential in the myocardium, the deformation caused by the muscles contraction, and the dynamics of the blood in the heart chambers and through the valves [31,71,76,77,96,98,99]. All these processes feature different temporal and spatial scales, yielding the so-called tyranny of scales [8]. A satisfactory knowledge of the models for isolated processes - the electrophysiology, the active and the passive mechanics - which we refer to as "single core models", is nowadays quite established; however, further theoretical and numerical studies are still necessary to better understand their mutual interactions [18,22,30,82]. In this work, we use state-of-the-art models in passive myocardial tissue modeling – the Holzapfel-Ogden model [46] – together with the active strain formulation [2,3] in combination with a recently proposed model for the transmurally heterogeneous thickening of the myocardium [7]; the latter is for the first time used in this paper for the integrated electromechanics modeling of the left ventricle (LV). The link between electrophysiology and mechanics is finally established through a model describing the shortening of the myocardial fibers [88], which is in turn triggered by a change in the ionic concentrations in the cardiac cells.
From the numerical viewpoint, we discretize in space the continuous single core models by means of the Finite Element Method (FEM) with both linear and quadratic elements, while the time discretization is carried on using Backward Differentiation Formulas (BDFs) of orders 1 and 2 [78]. BDFs are used both within implicit and semi-implicit schemes, the latter consisting in equal-order extrapolation of the unknowns in the nonlinear terms[16,37]; however, the semi-implicit scheme poses strong limitations on the choice of the timestep size with respect to the implicit scheme, even if the latter yields a nonlinear problem, which we approximate by means of the Newton method. Integrating the discrete single core models leads to the formulation of a monolithic algebraic coupled problem. Numerical coupling of cardiac electromechanics is however typically performed by means of segregated or staggered approaches, for which the electrophysiology and mechanics solvers are not applied simultaneously [4,5,18,38,54,84,107]. It is however known that partitioned schemes (segregated or staggered) do not generally guarantee unconditional stability [15], an issue not concerning instead the monolithic scheme.
Monolithic schemes for cardiac electromechanics of the heart were firstly proposed and analyzed in [24,25], where only the contraction phase was simulated and the blood pressure at the endocardium (the LV inner wall) was neglected. More recently, a monolithic scheme including a 0D model for the vascular system, heart valve, and atrial chamber model was proposed in [44]. However, other than being less explored, monolithic schemes have been used so far only for cardiac electromechanics within the framework of the active stress formulation and only for idealized or simplified LV geometries. Instead, in this paper we propose and successfully use a monolithic algorithm for the numerical simulation of cardiac electromechanics in the framework of the active strain formulation; in addition, we simulate full cardiac cycles (one heartbeat) for both idealized and patient-specific geometries. Essential components of our framework are: modeling FSI at the LV endocardium wall by means of
We solve the large, sparse block linear system stemming from the monolithic strategy by means of the GMRES method and, in order to speed-up the convergence of the linear solver, we employ a novel (right) preconditioner [89]. Our preconditioning strategy extends the one proposed in [27] for FSI problems and is based on the factorization of a block Gauss-Seidel preconditioner which exposes the multiphysics nature of the integrated problem, thus allowing to adapt its action on each single factor of the preconditioner. This strategy can be easily extended to a wide range of integrated multiphysics problems: the factors are indeed independently preconditioned by using algebraic multigrid [11] or domain decomposition [14,104] preconditioners tailored by exploiting physics-specific information (such as the number of equations) for each block; this local information would be instead lost by using a global monolithic preconditioner.
We use our monolithic electromechanics solver to perform numerical simulations in the High Performance Computing framework of the whole cardiac cycle, both for idealized and patient-specific LV geometries, 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 LV function and its properties.
The paper is organized as follows: in Sec. 2 we introduce the mathematical models for the electrophysiology, the mechanics and the activation of the myocardium; we then integrate them thus obtaining a continuous integrated model. In Sec. 3 we carry on the space and time numerical approximations of the single core models; in Sec. 4 we introduce the preconditioner used to solve the monolithic linear system arising after the time and space discretizations. In Sec. 5 we report and discuss the numerical results obtained with the proposed methods, and we finally draw our conclusions in Sec. 6.
The myocardium is a complex tissue composed of cardiomyocytes, organized in fibers and laminar collagen sheets [94], characterising the LV orthotropic internal structure. Therefore, the material properties of the LV are strongly dependent on the direction of fibers and sheets. First, the electrical conductivity of the cardiomyocytes is much larger in the longitudinal direction than transversally [52,75]; hence, at the macroscopic level, the transmembrane electric potential travels faster along the fibers. Secondly, it has been experimentally observed in [46] that the response of internal stresses to an external load significantly varies when measured among different directions since the myocardium is stiffer along the fibers. Finally, the contraction of the LV is made possible thanks to an active force acting on the cardiomyocytes lined up along the fibers direction [86]. In order to mathematically define fibers and sheets, we identify a local frame of reference by defining the mutually orthonormal vector fields
Electrophysiology models take into account the electrochemical reactions occuring in the myocardium [21] triggered by an electric impulse originated at the sino-atrial node and then conveyed through the Purkinje fibers [20]. Such signal causes a quick depolarization of the LV cardiomyocytes, meaning that the transmembrane potential (i.e. the difference in electric potential between the interior and the exterior of a cell) changes sign in few microseconds. The potential drives a change in the concentration of different ionic species flowing through the so-called ionic gates located on the cellular membrane; the concentration of these ionic species, in turn, influences the potential thus slowly repolarizing the cells. The continuous interaction between the ions concentration and the potential causes a cascade effect for which a fast traveling wave known as action potential propagates in the whole myocardium [57].
In this work, we consider the monodomain equation for the description of the evolution of the cellular transmembrane potential
{χ(Cm∂V∂t+Iion(V,w))=∇⋅(JF−1DmF−T∇V)+Iapp(t)in Ω0×(0,T),(JF−1DmF−T∇V)⋅N=0on ∂Ω0×(0,T),V=V0in Ω0×{0}. | (2.1) |
Here,
{dwidt=αi(V)(w∞i(V)−wi)+βi(V)wiin (0,T),wi(0)=wi,0,at t=0,for i=1,…,NI, | (2.2) |
where the unknowns
{∂V∂t+∑q∈{fi,so,si}Iq(V,w)=∇⋅(JF−1DmF−T∇V)+Iapp(t)in Ω0×(0,T),(JF−1DmF−T∇V)⋅N=0on ∂Ω0×(0,T),∂w∂t=α(V)(w∞(V)−w)+β(V)win Ω0×(0,T),V=V0,w=(1,1,0)Tin Ω0×{0}, | (2.3) |
since
α(V)=diag(1−HV1(V)τ−1(V),1−HV2(V)τ−2(V),1τ3(V)),β(V)=diag(−HV1(V)τ+1,−HV2(V)τ+2,0),w∞(V)=(1−HV−1(V),HVo(V)(w∞∗−1+Vτ∞2)+1−Vτ∞2,Hk3V3(V))T,Ifi(V,w1)=−HV1(V)(V−V1)(ˆV−V)τfiw1,Iso(V)=(1−HV2(V))(V−Vo)τo(V)+HV2(V)τso(V), |
and
An adequate mechanical model for the description of the myocardium's displacement must account for the tissue's complex behavior. Firstly, the internal stresses induced by a prescribed deformation are highly anisotropic [41] and in our formulation depend on the directions
We recall the momentum conservation equation in the reference configuration
{ρ∂2d∂t2−∇⋅P(d,γf)=0in Ω0×(0,T),q(d,∂d∂t)+P(d,γf)N=0on Γη0×(0,T),P(d,γf)N=pendo(t)Non Γendo0×(0,T),d=d0,∂d∂t=˙d0in Ω0×{0}. | (2.4) |
Here
P(d)=∂W(C)∂F. |
A very popular model for the myocardial tissue is the Holzapfel-Ogden [46] one. This is obtained by considering different contributions and by taking into account the anisotropic nature of the muscle:
˜W(C)=W1(I1)+W4f(I4f)+W4s(I4s)+W8fs(I8fs)=a2beb(I1−3)+∑i∈{f,s}ai2bi[eb⟨I4i−1⟩2−1]+afs2bfs[ebI28fs−1], |
where
The mechanical model should also account for the volumetric change to which the myocardium undergoes during the cardiac cycle. It has been observed in [19,113] that albeit this change is moderate, still it significantly ranges from 2% to 15%. For these reasons we use a nearly-incompressible formulation, which allows for small volume variations [28]. We multiplicatively decompose the deformation gradient
F=Fv¯F,Fv=J13I, | (2.5) |
where
Wvol(J)=B2(J−1)log(J), | (2.6) |
such that the larger is the bulk modulus
W=˜W+Wvol=W1(J−23I1)+W4f(I4f)+W4s(I4s)+W8fs(I8fs)+Wvol(J). |
We now proceed by modeling the active behavior of the myocardium driven by the stretching of the fibers by means of the active strain approach [2,3,63,86]. A virtual intermediate state
P=det(FA)PEF−TA,with PE=∂W(CE,J)∂FE. |
In Sec. 2.3, we will provide the explicit form of the tensor
P(d,FA)=aeb(J−23IE1−3)J−23(FEF−1A−IE13(FEFA)−T)+2afebf⟨IE4f−1⟩2⟨IE4f−1⟩(fE⊗fA)+2asebs⟨IE4s−1⟩2⟨IE4s−1⟩(sE⊗sA)+afsebfs(IE8fs)2IE8fs(fE⊗sA+sE⊗fA)+B2(J+Jlog(J)−1)(FEFA)−T, | (2.7) |
where
The activation model of the myocardium represents the link between the electrophysiology and the tissue (passive) mechanics. Instead of cellular models describing the complex dynamic taking place inside the sarcomeres [81], we exploit a phenomenological model for the local shortening of the fibers
{g(c)∂γf∂t−εΔγf=Φ(c,γf,d)in Ω0×(0,T),∇γf⋅N=0on ∂Ω0×(0,T),γf=0in Ω0×{0}. | (2.8) |
Here,
The contraction of the tissue is triggered by the calcium concentration exceeding a threshold value
Being
FA=I+γff0⊗f0+γss0⊗s0+γnn0⊗n0, |
which is symmetric and we set
γn=k′(λ)(1√1+γf−1). |
Indeed, the latter depends on an independent variable
k′(λ)=¯k′(¯kendoλ−λepiλendo−λepi+¯kepiλ−λendoλepi−λendo). | (2.9) |
Finally, in order to have
Writing together Eqs.(2.3), (2.4), and (2.8), the coupled electromechanics problem finally reads:
{∂V∂t+∑q∈{fi,so,si}Iq(V,w)=∇⋅(JF−1DmF−T∇V)−Iapp(t)in Ω0×(0,T),∂w∂t=α(V)(w∞(V)−w)+β(V)win Ω0×(0,T),ρ∂2d∂t2−∇0⋅P(d,γf)=0in Ω0×(0,T),g(w3)∂γf∂t−εΔγf=Φ(w3,γf,d)in Ω0×(0,T),(JF−1DmF−T∇V)⋅N=0on ∂Ω0×(0,T),q(d,∂d∂t)+P(d)N=0on Γη0×(0,T),P(d)N=pendo(t)Non Γendo0×(0,T),∇γf⋅N=0on ∂Ω0×(0,T),V=V0,w=(1,1,0)T,γf=0in Ω0×{0},d=d0,∂d∂t=˙d0in Ω0×{0}. | (2.10) |
We remark that at this stage the pressure
A common problem which has to be considered in biological modeling involving a fluid-structure interface is that the reference geometry
Pressure preload ([30,84,101,103]): the reference geometry
Pressure prestress ([47,100]): we compute internal stresses distribution such that the reference geometry is in equilibrium with the blood pressure
We adopt the pressure prestress approach since in the preload one the procedure returns a configuration which might be significantly different with respect to the initial geometry as shown in [47]. Indeed, we already know the loaded configuration from medical images. In order to compute
{∇0⋅P(d)=−∇0⋅P0in Ω0,(N⊗N)Kη⊥d+(I−N⊗N)Kη∥d+P(d)N=0on Γη0,P(d)N=pendo(0)Non Γendo0. | (2.11) |
Eq. (2.11) is the steady, passive version of problem (2.4) with the additive decomposition of the strain tensor. Moreover, we set:
pk=kS¯pendo,k=1,…,S, |
where
Geometry | Mesh | # Vertices | # Elements | |
Idealized(Fig. 6) | (a) | 4.6 mm | 1'827 | 3'010 |
(b) | 2.3 mm | 11'658 | 12'040 | |
(c) | 1.2 mm | 81'335 | 416'000 | |
Patient-specific(Fig. 7) | (d) | 0.8 mm | 126'031 | 637'379 |
In this section we describe the numerical approximation of the electromechanics problem (2.10). The FEM [78] is used for the space discretization of the PDEs, while BDF [78] are used for the time discretization.
We consider a mesh composed of tetrahedrons
*The MRI images are provided by Prof. J. Schwitter (Chief physician at the Centre Hospitalier Universitaire Vaudois CHUV, Lausanne) and Dr. P. Masci in the framework of the collaboration CMCS@EPFL-CHUV hospital.
We use the FEM for the spatial approximation of the monodomain equation (2.1). We first introduce the finite dimensional space
∫Ω0˙VhψidΩ0+∫Ω0(JF−1hDmF−Th∇Vh)⋅∇ψidΩ0+∑q∈{fi,so,si}∫Ω0Iq(Vh,wh)ψidΩ0=∫Ω0Iapp(t)ψidΩ0,∀i=1,…,NdofV, | (3.1) |
with
{M˙Vh(t)+K(dh(t))Vh(t)+Iion(Vh(t),wh(t))=Iapp(t)t∈(0,T],Vh(0)=V0,h, | (3.2) |
where
(Iion(Vh,wh))i=∑q∈{fi,so,si}∫Ω0Iq(Vh,wh)ψidΩ0,(Iapp(t))i=∫Ω0Iapp(t)ψidΩ0. | (3.3) |
We will discuss in Sec. 3.1.3 two different strategies for the approximation of the nonlinear term
The ionic model (2.2) is a system of ODEs that indirectly depends on the space variable through the transmembrane potential
wh(t)={wlh(t)}NIl=1andwlh(t)={wlj(t)}Ndofwj=1. | (3.4) |
The problem thus obtained reads: given
˙wlj+(αl(Vj)−βl(Vj))wlj=αl(Vj)w∞l(Vj), | (3.5) |
with
{˙wh(t)+U(Vh(t))wh(t)=Q(Vh(t)),t∈(0,T],wh(0)=w0,h, |
where the block matrix
As anticipated in Sec. 3.1.1, different choices can be made for the approximation of the term
State Variable Interpolation (SVI): the variables
∫Ω0Iq(Vh,wh)ψidΩ0=∑K∈Th(Nq∑s=1Iq(NdofV∑j=1Vj(t)ψj(xKs),NdofV∑j=1wj(t)ψj(xKs))ψi(xKs)ωKs). |
This approach corresponds to the standard FEM approximation.
Ionic Currents Interpolation (ICI): for each element,
∫Ω0Iq(Vh,wh)ψidΩ0≈∑K∈Th(Nq∑s=1NdofV∑j=1Iq(Vj(t),wj(t))ψj(xKs)ψi(xKs)ωKs). |
The two approaches are equivalent if
As previously done for the monodomain equation, we use the FEM to approximate the momentum equation (2.4). We denote by
∫Ω0ρs∂2dh∂t2⋅ψidΩ0+∫Ω0P(dh,γf,h):∇0ψidΩ0+∑η∈{epi,base}∫Γη0q(dh,∂dh∂t)⋅ψidΓ0=∫Γendo0pendo(t)N⋅ψidΓ0,∀i=1,…,Ndofd, |
with
We rewrite it in algebraic form as:
{ρsM¨dh(t)+F˙dh(t)+Gdh(t)+S(dh(t),γf,h(t))=pendo(t)t∈(0,T],dh(0)=d0,h,˙dh(0)=˙d0,h, |
where
Finally, we once again use FEM to discretize in space the equation for the unknown
∫Ω0∂γf,h∂tψidΩ0+ε∫Ω01g(ch)∇γf,h⋅∇ψidΩ0−∫Ω01g(ch)Φ(ch,γf,h,dh)ψidΩ0=0,∀i=1,…,Ndofγf, |
with
{M˙γf,h(t)+εK(c(t))γf,h(t)+Φ(c(t),γf,h(t),dh(t))=0t∈(0,T],γf,h(0)=0. |
After having applied the space discretization to the single core models, we obtain the semi-discretized formulation of problem (2.10) in the form of a nonlinear system of ODEs. We denote by
{Mz(t)+T(z(t))=h(t)t∈(0,T],z(0)=z0,˙zd(0)=˙d0,h, | (3.6) |
where
˙zi(tn+1)≈1Δt(ϑ(Ⅰ)0zn+1i−zRHSi),zRHSi=σ∑k=1ϑ(Ⅰ)kzn−k+1i,¨zi(tn+1)≈1(Δt)2(ϑ(Ⅱ)0zn+1i−zRHSi),zRHSi=σ+1∑k=1ϑ(Ⅱ)kzn−k+1i, |
where
In the implicit case, we obtain the following nonlinear system:
A(zn+1)=bn+1n=σ,…,NT−1, | (3.7) |
with
In the semi-implicit case, on the other hand, we extrapolate the variables in the nonlinear term
zi(tn+1)≈z∗i=σ∑k=1βkzn−k+1i. |
We then avoid the nonlinear terms by partially evaluating
A(zn+1)≈A(z∗)zn+1+gn+1for n=σ,…,NT−1. |
By recalling Eq.(3.7), we hence obtain a system in the form:
A(z∗)zn+1=˜bn+1n=σ,…,NT−1, | (3.8) |
with
By applying the Newton method [78] to (3.7), we iteratively solve for
{Jn+1kδzn+1k+1=−rn+1k,zn+1k+1=zn+1k+δzn+1k+1, | (3.9) |
for
Jn+1k=[Jw(Vn+1k)JwV(wn+1k,Vn+1k)00JVw(wn+1k,Vn+1k)JV(dn+1k)0JVd(Vn+1k,dn+1k)Jγfw(wn+1k,γfn+1k,dn+1k)0Jγf(wn+1k,γfn+1k,dn+1k)Jγfd(wn+1k,γfn+1k,dn+1k)00Jdγf(γfn+1k,dn+1k)Jd(γfn+1k,dn+1k)], |
while the residual is defined as
In this case, the block structure of the matrix associated to linear system (3.8) to be solved is sparser with respect to the implicit one, as
\begin{equation} \label{eq:semi_implicit_pattern} \mathbb{A}({\bf z}^*) = \left[ \begin{array}{c|c|c|c} \mathbb{A}_{{\bf w}}({\bf V}^*)& 0&0&0 \\ \hline \mathbb{A}_{{\bf V}{\bf w}}({\bf w}^*, {\bf V}^*) &\mathbb{A}_{{\bf V}}({\bf d}^*§)& 0&0 \\ \hline 0&0 & \mathbb{A}_{\mathit{\boldsymbol{\gamma}}_f}({\bf w}^*, {\mathit{\boldsymbol{\gamma}}_f}^*, {\bf d}^*)&0 \\ \hline 0&0&0& \mathbb{A}_{{\bf d}}({\mathit{\boldsymbol{\gamma}}_f}^*, {\bf d}^*) \end{array} \right]. \end{equation} | (3.10) |
The matrix (3.10) is block lower triangular thanks to the ordering of the variables in the electrophysiology model (3.6): indeed, by swapping the ionic and the electric blocks we would have obtained a block upper triangular matrix. The semi-implicit scheme has the clear advantage of avoiding, at each timestep, the solution of a nonlinear problem. Moreover, both the assembly of the system matrix
On the other hand, semi-implicit schemes may impose strong limitations on the timestep size
We will further investigate this issue in the future; however, for the next numerical simulations of this paper we will consider only the implicit scheme.
We rewrite, for the sake of simplicity, the linear system associated to a single Newton iteration for the implicit scheme (3.9) in the following form:
\begin{equation} \label{linear_problem} \mathbb{J} \, \delta {\bf z} = -{\bf r}, \qquad \text{with } \mathbb{J} = \left[ \begin{array}{c|c|c|c} J_{11}& J_{12}&0&0 \\ \hline J_{21} & J_{22}&0 & J_{24} \\ \hline J_{31}&0 & J_{33} & J_{34} \\ \hline 0&0 & J_{43}& J_{44} \end{array} \right]. \end{equation} | (4.1) |
We use the GMRES method [89] for the solution of the linear problem (4.1) as
The choice of the preconditioner is critical in order to ensure the convergence of the linear solver; while this is true in general, it is even more relevant in the case of monolithic multiphysics problems [51]. Indeed, using a black-box algebraic preconditioner for problem (4.1) the information related to each differential problem associated to a single core model would be neglected; we instead consider a strategy exploiting such information at the block level, that is we create a preconditioner that exploits the "physics" of the coupled problem at the level of the block structure. Following [26,27,35] for FSI problems, we propose a block Gauss-Seidel preconditioner
\begin{equation*} \mathcal{P} = \left[ \begin{array}{c|c|c|c} J_{11}&0&0&0 \\ \hline J_{21} & J_{22}&0&0 \\ \hline J_{31}&0 & J_{33}&0 \\ \hline J&0 & J_{43}&J_{44} \end{array} \right]. \end{equation*} |
\begin{equation*} \mathcal{P} = \underbrace{ \left[ \begin{array}{c|c|c|c} J_{11}&0&0&0 \\ \hline 0&I&0&0 \\ \hline 0&0&I&0 \\ \hline 0&0&0&I \end{array} \right]}_{\mathcal{P}_1 = \mathcal{P}_{ion}} \underbrace{ \left[ \begin{array}{c|c|c|c} I&0&0&0 \\ \hline J_{21} & J_{22}&0&0 \\ \hline 0&0&I&0 \\ \hline 0&0&0&I \end{array} \right]}_{\mathcal{P}_2 = \mathcal{P}_{pot}} \underbrace{ \left[ \begin{array}{c|c|c|c} I&0&0&0 \\ \hline 0&I&0&0 \\ \hline J_{31}&0 & J_{33}&0 \\ \hline 0&0&0&I \end{array} \right]}_{\mathcal{P}_3 = \mathcal{P}_{act}} \underbrace{ \left[ \begin{array}{c|c|c|c} I&0&0&0 \\ \hline 0&I&0&0 \\ \hline 0&0&I&0 \\ \hline 0&0&J_{43}& J_{44} \end{array} \right].}_{\mathcal{P}_4 = \mathcal{P}_{mec}} \end{equation*} |
This decomposition can also be interpreted as an inexact factorization of matrix
\begin{equation} \label{eq:precond_linear_problem} \left\{ \begin{aligned} &\mathbb{J} \mathcal{P}^{-1}{\bf y} = -{\bf r} \\ &\mathcal{P} \delta {\bf z} = {\bf y}. \end{aligned} \right. \end{equation} | (4.2) |
Since each diagonal block
\begin{equation} \label{precond_application} \left\{ \begin{aligned} &\delta {\bf z}_{{\bf w}} && = J_{11}^{-1} {\bf y}_{{\bf w}}, \\ &\delta {\bf z}_V && = J_{22}^{-1} ({\bf y}_V - J_{21} \delta {\bf z}_{{\bf w}} ), \\ &\delta {\bf z}_{{\bf d}} && = J_{33}^{-1} ({\bf y}_{{\bf d}} - J_{31} \delta {\bf z}_{{\bf w}} ), \\ &\delta {\bf z}_{\mathit{\boldsymbol{\gamma}}_f} && = J_{44}^{-1} ({\bf y}_{\mathit{\boldsymbol{\gamma}}_f} - J_{43} \delta {\bf z}_{{\bf d}} ); \end{aligned} \right. \end{equation} | (4.3) |
the solution of each linear system in (4.3) is carried out by using again the GMRES method and by exploiting Algebraic Multigrid (AMG) preconditioners [11]. This new preconditioner
For our simulations we will consider a full heartbeat, by taking the conventional duration of
(1) an isovolumic contraction phase in which the endocardial pressure
(2) an ejection phase characterized by a decrement in the ventricular volume
(3) an isovolumic relaxation phase during which
(4) a filling phase starting with the opening of the mitral valve causing an increment of
When necessary, we compute the volume
We then model the endocardial pressure
1) Isovolumic phase I: At each timestep, we iteratively solve problem (2.10) by updating at each subiteration the pressure in the following manner:
\begin{equation} \label{eq:isochoric} p^{n+1}_{k+1} = p^{n+1}_k+\frac{V^{n+1}_k-V^n}{C_p} \qquad k = 0, \dots \end{equation} | (4.4) |
with
2) Ejection phase: A two-element Windkessel model [110] of the form:
\begin{equation} \label{eq:windkessel} \left\{ \begin{aligned} &C\frac{dp}{dt} = -\frac{p}{R}-\frac{dV}{dt} \quad t \in (T', T"] \\[0.1cm] &p(T') = \overline{p} \end{aligned} \right. \end{equation} | (4.5) |
is solved in the pressure variable with a BDF scheme of order
3) Isovolumic phase II: Modeled as the first isovolumic phase. This phase ends when the pressure reaches the minimal value of
4) Filling phase: The pressure is linearly increased at each timestep so that it reaches the EDP value at the final time
Even if (4.5) is strongly coupled with (2.10) during the ejection phase, for simplicity we solve the two problems in a segregated fashion.
We now describe the geometries used for the numerical simulation of the integrated electromechanics problem. As anticipated in Sec. 3, we use both an idealized prolate geometry and a patient-specific geometry segmented from MRI images as detailed in Sec. 5.1. Moreover, in Sec. 5.2, we outline the procedure that is exploited to define the fibers and the sheets fields on the different meshes. We conclude this section by showing and analyzing the results of the simulations.
Image segmentation is the process of locating regions of interest (ROI) in the form of a subset of pixels [42]. In biomedical applications, this accounts to assign different flags to regions containing different types of tissues and/or fluids. This result, depending on the properties of the biological material which is to be segmented, can often be achieved through semi-automatic or automatic procedures (see e.g. [32,49,106] for arteries and blood vessels and [10,109] for the Purkinje network), exploiting a large set of different techniques. However, since the development of algorithms for the segmentation of the chambers of the heart [72,114,115] is beyond the scope of this work, we used a manual procedure based on the brightness of the pixels of the aforementioned 3D MRI images. We first apply a threshold filter to select the smallest set of pixels containing the whole myocardium; then, we manually remove artifacts and features that we decide to neglect such as the papillary muscles and the upper part of the LV. Finally, a Gaussian smoothing procedure [111] was performed to improve the quality of the mesh. See Figure 5.
The image was taken at the end of the diastolic phase, that is when the LV reaches its maximum enlargement. The internal volume of the patient-specific myocardium at this stage of the heartbeat measures approximately 95 ml.
Unlike the geometry, the fibers and sheets field distribution in the myocardium may not be extracted from medical images such as MRIs unless special procedures are applied [80]. For this reason, several mathematically rule-based definition of the fields have been used in literature [38,56,61,86], which try to approximate their orientation. At the epicardium and at the endocardium, the fiber direction is tangential to the boundary, while the sheet direction belongs to the plane identified by the normal and the ventricle centerline. Then, in the most general case, angles
We setup a total of five numerical simulations, each for a single heartbeat, using the four meshes shown in Figures 6-7. The monodomain, the activation, and the mechanics equations are approximated with
Geometry | Mesh | # Vertices | # Elements | |
Idealized(Fig. 6) | (a) | 4.6 mm | 1'827 | 3'010 |
(b) | 2.3 mm | 11'658 | 12'040 | |
(c) | 1.2 mm | 81'335 | 416'000 | |
Patient-specific(Fig. 7) | (d) | 0.8 mm | 126'031 | 637'379 |
Geometry | Mesh | # CPUs | ||
Idealized | (a) | 1 | 14'616 | 1 |
(b) | 1 | 93'264 | 6 | |
(c) | 1 | 650'680 | 45 | |
(a) | 2 | 93'264 | 6 | |
Patient-specific | (d) | 1 | 1'008'248 | 72 |
We use LifeV†, an open-source finite element library for the solution of problems described by PDEs in a High Performance Computing framework, to implement the numerical methods. All the computations were carried out using Piz Daint, a Cray XC50/XC40 supercomputer installed at the Swiss National Supercomputing Center (CSCS)‡.
As detailed in Sec. 4 we exploit the Newton method for the solution of the nonlinear system (3.7) by setting a maximum number of 5 iterations per timestep. Moreover we use the relative residual as stopping criterion, with a tolerance equal to
For simplicity, we initiate the electric signal propagation by applying at the same time 2-milliseconds long stimuli at three distinct points on the endocardium; we remark that, even if this pacing strategy is not fully realistic since the electrical signal delivered by the Purkinje network activates asynchronously hundreds of cells [83,108], it provides physically meaningful results. The other parameters used for the simulations are reported in Table 3.
| | | | | | | | | | | | | |
| | | | | | | | | | | 12 | 500 | |
| | | | | | | | | | | | | |
50 | 10 | 1500 | 0 | | 5 | 1 | 0 | 0 | | | | |
We compare in Figure 8 the transmembrane potential obtained using
We observe from Figure 9 a difference of several milliseconds in the complete activation of the myocardium, which is significant in the electrophysiology characteristic time scales. However, because of the multiscale nature of the integrated problem, the mechanics is not affected by such delay as can be observed in Figure 11. Here, the only remarkable difference between the three results is the underestimation of the displacement magnitude in the simulation with the coarsest mesh (a). In this case, the number of degrees of freedom is not sufficient to either fullfill the nearly-incompressible constraint nor to represent the large deformations observable in the other two cases. The deformation of the patient-specific mesh (d) at the same times is shown in Figure 12.
In both the idealized and the patient-specific cases, a significant thickening of the myocardium walls takes place, which is in accordance with experimental observations [79]. The model, with the set of parameters in Table 3, produces however a significant rotation of the LV: a recent work [74] suggests that this behavior is related to the choice of the incompressibility constraint, the bulk modulus
In order to better evaluate the ability to describe the LV activity of the integrated numerical model for a full cardiac cycle, pressure-volume (pV) loops are also represented. A comparison with in-vivo measurements would be pointless in the case of the ideal geometry, however we stress the fact that the pV-loops reported in Figure 14 are in both qualitative and quantitative agreement with those observed experimentally as e.g. in [87]; moreover, the maximal endocardial pressure, which is reached during the ejection phase, is approximately 110 mmHg and hence lays well within the physiological range of average individuals [91]. The difference between the four cases in Figure 14 accounts to a maximum variation of about 2 % in the minimal ventricle volume during the second isochoric phase.
The pV-loop for the patient-specific case is reported in Figure 15. Experimental measurements of the pressure and of the volume for the patient under study were not available, therefore we did not carry on a quantitative comparison of our results against patient-specific data. However, we once again observe that the values of the pV-loop are in line with real data [87].
We then report in Table 4, the timestep used for each simulation, the average number of Newton iterations, the average number of fixed point iterations for each of the two isochoric phases, and the average time spent for different operations during a single timestep. The overall time required for the complete simulations is also reported. Table 4 clearly shows that finer meshes require larger computational resources for a full heartbeat simulation (
Geometry | Mesh | |||||||||
Idealized | (a) | 1 | 3.02 | 2.31 | 2.55 | 2.15 | 0.023 | 0.006 | 5h 41m | |
(b) | 1 | 3.37 | 2.14 | 3.44 | 6.67 | 0.163 | 0.019 | 18h 21m | ||
(c) | 1 | 4.23 | 2.49 | 3.26 | 9.19 | 0.632 | 0.040 | 46h 42m | ||
(a) | 2 | 3.29 | 1.68 | 2.11 | 6.47 | 0.178 | 0.025 | 30h 18m | ||
Patient-specific | (d) | 1 | 4.44 | 1.96 | 4.51 | 22.33 | 1.692 | 0.076 | 68h 37m |
Finally, we test the strong scalability of the proposed monolithic solver for cardiac electromechanics. With this aim, we further refine the ideal mesh of Figure 6, thus obtaining a finer one featuring 4'629'817 vertices and 26'624'000 tetrahedra, which is more adequate for the strong scalability test. We then set
§Performing the strong scalability test with a number of CPUs smaller than 800 was not feasible with this mesh; indeed, memory limitations occurred due the excessive number of unknowns per single CPU.
In this work we proposed a novel monolithic solver for the simulation of the integrated electromechanics problem for the LV of the human heart. We coupled the monodomain equation, the ionic minimal model, the activation model for the fibers contraction, and the myocardial mechanics in the active strain framework with a transmurally variable activation. The interaction of the myocardium with the blood inside the ventricle was taken into account by prestressing and by solving additional 0D problems for the fluid.
We considered both implicit and semi-implicit BDF schemes for our monolithic solver together with a preconditioning strategy based on block Gauss-Seidel preconditioner which exploits the multiphysics structure of the coupled problem. The strong limitations on the timestep size in the semi-implicit case are such that the implicit scheme is by far the most efficient for numerical simulations lasting for one or multiple heartbeats. We solve the coupled problem in the High Performance Computing framework and we performed a scalability analysis for a benchmark test involving more than 37'000'000 unknowns.
The results obtained by simulating a complete heartbeat with both ideal and patient-specific geometries are qualitatively in agreement with physiological values measured in healthy patients. Moreover, the comparison of the simulations on the idealized geometry with different mesh sizes shows that the proposed numerical model is able to capture the physical solution with satisfactory accuracy even in the case of relatively coarse meshes. From our simulations we observe that choosing a mesh with maximum size
This research was partially supported by the Swiss Platform for Advanced Scientific Computing (PASC, project "Integrative HPC Framework for Coupled Cardiac Simulations"). We also gratefully acknowledge the Swiss National Supercomputing Center (CSCS) for providing the CPU resources for the numerical simulations under project ID s635. Finally, we acknowledge the support from the ERC Advanced Grant iHEART, “An Integrated Heart Model for the simulation of the cardiac function”, 2017-2022, P.I.A. Quarteroni (ERC-2016-ADG, project ID: 740132).
We acknowledge Prof. J. Schwitter and Dr. P. Masci (CHUV) for providing the MRI images used in this work and for the constructive discussions on the topic from a physician's perspective. We also thank Prof. P. Tozzi (CHUV) for the invaluable insights in the functioning of the human heart. Finally, we acknowledge both Dr. D. Forti and MER S. Deparis (EPFL) for the useful discussions and their precious help with the LifeV library.
The authors declare no conflict of interest.
[1] | Burger M, Caffarelli LA, Markowich PA, et al. (2013) On a Boltzmann-type price formation model. P R Soc Lond Ser A Math Phys Eng Sci 469: 1-20. |
[2] |
Caffarelli LA, Markowich PA, Pietschmann JF (2011) On a price formation free boundary model by Lasry and Lions. C R Math Acad Sci Paris 349: 621-624. doi: 10.1016/j.crma.2011.05.011
![]() |
[3] | Carmona R, Delarue F (2014) The master equation for large population equilibriums. arXiv:1404.4694. |
[4] |
Carmona R, Delarue F, Lacker D (2016) Mean field games with common noise. Ann Probab 44: 3740-3803. doi: 10.1214/15-AOP1060
![]() |
[5] | Alasseur C, Taher IB, Matoussi A (2017) An extended mean field game for storage in smart grids. J Optim Theory Appl 184: 644-670. |
[6] |
Couillet R, Perlaza SM, Tembine H, et al. (2012) Electrical vehicles in the smart grid: A mean field game analysis. IEEE J Sel Area Commun 30: 1086-1096. doi: 10.1109/JSAC.2012.120707
![]() |
[7] | De Paola A, Angeli D, Strbac G (2016) Distributed control of micro-storage devices with mean field games. IEEE T Smart Grid 7: 1119-1127. |
[8] |
De Paola A, Trovato V, Angeli D (2019) A mean field game approach for distributed control of thermostatic loads acting in simultaneous energy-frequency response markets. IEEE T Smart Grid, 10: 5987-5999. doi: 10.1109/TSG.2019.2895247
![]() |
[9] | Gomes D, Lafleche L, Nurbekyan L (2016) A mean-field game economic growth model, In: Proceedings of the American Control Conference, 4693-4698. |
[10] | Gomes DA, Saúde J (2020) A mean-field game approach to price formation. Dyn Games Appl, To appear. |
[11] | Graber J, Mouzouni C (2017) Variational mean field games for market competition. arXiv:1707.07853. |
[12] | Guéant O, Lasry JM, Lions PL (2011) Mean field games and applications. In: Paris-Princeton Lectures on Mathematical Finance 2010, Berlin: Springer, 205-266. |
[13] | Kizilkale AC, Malhame RP (2014) A class of collective target tracking problems in energy systems: Cooperative versus non-cooperative mean field control solutions. In: Proceedings of the IEEE Conference on Decision and Control, 3493-3498. |
[14] | Kizilkale AC, Malhame RP (2014) Collective target tracking mean field control for electric space heaters. In: 2014 22nd Mediterranean Conference on Control and Automation, MED 2014, 829-834. |
[15] |
Kizilkale AC, Malhame RP (2014) Collective target tracking mean field control for markovian jump-driven models of electric water heating loads. IFAC Proceedings Volumes 47: 1867-1872. doi: 10.3182/20140824-6-ZA-1003.00630
![]() |
[16] |
Malhamé R, Chong CY (1988) On the statistical properties of a cyclic diffusion process arising in the modeling of thermostat-controlled electric power system loads. SIAM J Appl Math 48: 465-480. doi: 10.1137/0148026
![]() |
[17] | Malhamé R, Kamoun S, Dochain D (1989) On-line identification of electric load models for load management. In: Advances in Computing and Control, Berlin: Springer, 290-304. |
[18] |
Markowich PA, Matevosyan N, Pietschmann JF, et al. (2009) On a parabolic free boundary equation modeling price formation. Math Mod Meth Appl Sci 19: 1929-1957. doi: 10.1142/S0218202509003978
![]() |
1. | Francesco Regazzoni, Luca Dedè, Alfio Quarteroni, Active Force Generation in Cardiac Muscle Cells: Mathematical Modeling and Numerical Simulation of the Actin-Myosin Interaction, 2020, 2305-221X, 10.1007/s10013-020-00433-z | |
2. | Matteo Salvador, Luca Dede’, Alfio Quarteroni, An intergrid transfer operator using radial basis functions with application to cardiac electromechanics, 2020, 66, 0178-7675, 491, 10.1007/s00466-020-01861-x | |
3. | F. Regazzoni, A. Quarteroni, An oscillation-free fully staggered algorithm for velocity-dependent active models of cardiac mechanics, 2021, 373, 00457825, 113506, 10.1016/j.cma.2020.113506 | |
4. | F. Regazzoni, L. Dedè, A. Quarteroni, Machine learning of multiscale active force generation models for the efficient simulation of cardiac electromechanics, 2020, 370, 00457825, 113268, 10.1016/j.cma.2020.113268 | |
5. | Stefano Pagani, Luca Dede’, Andrea Manzoni, Alfio Quarteroni, Data integration for the numerical simulation of cardiac electrophysiology, 2021, 0147-8389, 10.1111/pace.14198 | |
6. | Francesco Regazzoni, Luca Dedè, Alfio Quarteroni, Daniel A Beard, Biophysically detailed mathematical models of multiscale cardiac active mechanics, 2020, 16, 1553-7358, e1008294, 10.1371/journal.pcbi.1008294 | |
7. | Francesco Regazzoni, Luca Dedè, Alfio Quarteroni, Active contraction of cardiac cells: a reduced model for sarcomere dynamics with cooperative interactions, 2018, 17, 1617-7959, 1663, 10.1007/s10237-018-1049-0 | |
8. | Luca Barbarotta, Simone Rossi, Luca Dedè, Alfio Quarteroni, A transmurally heterogeneous orthotropic activation model for ventricular contraction and its numerical validation, 2018, 34, 2040-7939, 10.1002/cnm.3137 | |
9. | Luca Pegolotti, Luca Dedè, Alfio Quarteroni, Isogeometric Analysis of the electrophysiology in the human heart: Numerical simulation of the bidomain equations on the atria, 2019, 343, 00457825, 52, 10.1016/j.cma.2018.08.032 | |
10. | L. Dede’, A. Gerbi, A. Quarteroni, 2020, Chapter 3, 978-3-030-45196-7, 81, 10.1007/978-3-030-45197-4_3 | |
11. | Mahmoud Moustafa, Mohd Hafiz Mohd, Ahmad Izani Ismail, Farah Aini Abdullah, Dynamical analysis of a fractional order eco-epidemiological model with nonlinear incidence rate and prey refuge, 2021, 65, 1598-5865, 623, 10.1007/s12190-020-01408-6 | |
12. | Hongyan Zhao, Fangxin Zhou, Huaping Liu, Deep learning for diplomatic video analysis, 2020, 79, 1380-7501, 4811, 10.1007/s11042-018-6650-9 | |
13. | Zhaolun Wang, Numerical Analysis of Deformation Control of Deep Foundation Pit in Ulanqab City, 2021, 0960-3182, 10.1007/s10706-021-01836-6 | |
14. | Alessandro Nitti, Josef Kiendl, Alessio Gizzi, Alessandro Reali, Marco D. de Tullio, A curvilinear isogeometric framework for the electromechanical activation of thin muscular tissues, 2021, 382, 00457825, 113877, 10.1016/j.cma.2021.113877 | |
15. | Balaram Murthy Chintakindi, Mohammad Farukh Hashmi, SSAD: Single-Shot Multi-Scale Attentive Detector for Autonomous Driving, 2023, 0256-4602, 1, 10.1080/02564602.2023.2176932 | |
16. | Chengsheng Luo, Linjun Wang, Youxiang Xie, Baojia Chen, A New Conjugate Gradient Method for Moving Force Identification of Vehicle–Bridge System, 2022, 2523-3920, 10.1007/s42417-022-00824-1 | |
17. | Junhua Liu, Xiaoli Liu, Fu Wang, Xuejiao Ma, Xiaohang Yue, Equilibrium decisions and coordination of a logistics service supply chain with bilateral corporate social responsibility, 2022, 56, 0399-0559, 2595, 10.1051/ro/2022115 | |
18. | F. Regazzoni, A. Quarteroni, Accelerating the convergence to a limit cycle in 3D cardiac electromechanical simulations through a data-driven 0D emulator, 2021, 135, 00104825, 104641, 10.1016/j.compbiomed.2021.104641 | |
19. | F. Regazzoni, M. Salvador, P.C. Africa, M. Fedele, L. Dedè, A. Quarteroni, A cardiac electromechanical model coupled with a lumped-parameter model for closed-loop blood circulation, 2022, 457, 00219991, 111083, 10.1016/j.jcp.2022.111083 | |
20. | Roberto Piersanti, Francesco Regazzoni, Matteo Salvador, Antonio F. Corno, Luca Dede’, Christian Vergara, Alfio Quarteroni, 3D–0D closed-loop model for the simulation of cardiac biventricular electromechanics, 2022, 391, 00457825, 114607, 10.1016/j.cma.2022.114607 | |
21. | Michele Bucelli, Alberto Zingaro, Pasquale Claudio Africa, Ivan Fumagalli, Luca Dede', Alfio Quarteroni, A mathematical model that integrates cardiac electrophysiology, mechanics, and fluid dynamics: Application to the human left heart, 2023, 2040-7939, 10.1002/cnm.3678 | |
22. | Barış Cansız, Michael Kaliske, A comparative study of fully implicit staggered and monolithic solution methods. Part I: Coupled bidomain equations of cardiac electrophysiology, 2022, 407, 03770427, 114021, 10.1016/j.cam.2021.114021 | |
23. | Alfio Quarteroni, Luca Dedè, Francesco Regazzoni, Modeling the cardiac electromechanical function: A mathematical journey, 2022, 59, 0273-0979, 371, 10.1090/bull/1738 | |
24. | F. Regazzoni, M. Salvador, L. Dede’, A. Quarteroni, A machine learning method for real-time numerical simulations of cardiac electromechanics, 2022, 393, 00457825, 114825, 10.1016/j.cma.2022.114825 | |
25. | Feng Wang, Guanghui Song, Jingqi Xuan, Han Wu, 2021, Chapter 201, 978-3-030-70664-7, 1860, 10.1007/978-3-030-70665-4_201 | |
26. | Tobias Gerach, Steffen Schuler, Jonathan Fröhlich, Laura Lindner, Ekaterina Kovacheva, Robin Moss, Eike Moritz Wülfers, Gunnar Seemann, Christian Wieners, Axel Loewe, Electro-Mechanical Whole-Heart Digital Twins: A Fully Coupled Multi-Physics Approach, 2021, 9, 2227-7390, 1247, 10.3390/math9111247 | |
27. | Matteo Salvador, Marco Fedele, Pasquale Claudio Africa, Eric Sung, Luca Dede', Adityo Prakosa, Jonathan Chrispin, Natalia Trayanova, Alfio Quarteroni, Electromechanical modeling of human ventricles with ischemic cardiomyopathy: numerical simulations in sinus rhythm and under arrhythmia, 2021, 136, 00104825, 104674, 10.1016/j.compbiomed.2021.104674 | |
28. | Matteo Salvador, Francesco Regazzoni, Stefano Pagani, Luca Dede', Natalia Trayanova, Alfio Quarteroni, The role of mechano-electric feedbacks and hemodynamic coupling in scar-related ventricular tachycardia, 2022, 142, 00104825, 105203, 10.1016/j.compbiomed.2021.105203 | |
29. | Marco Fedele, Roberto Piersanti, Francesco Regazzoni, Matteo Salvador, Pasquale Claudio Africa, Michele Bucelli, Alberto Zingaro, Luca Dede’, Alfio Quarteroni, A comprehensive and biophysically detailed computational model of the whole human heart electromechanics, 2023, 410, 00457825, 115983, 10.1016/j.cma.2023.115983 | |
30. | Ludovica Cicci, Stefania Fresca, Elena Zappon, Stefano Pagani, Francesco Regazzoni, Luca Dede', Andrea Manzoni, Alfio Quarteroni, 2023, 9780323899673, 403, 10.1016/B978-0-32-389967-3.00028-7 | |
31. | Michele Bucelli, Martin Geraint Gabriel, Alfio Quarteroni, Giacomo Gigante, Christian Vergara, A stable loosely-coupled scheme for cardiac electro-fluid-structure interaction, 2023, 490, 00219991, 112326, 10.1016/j.jcp.2023.112326 | |
32. | Michele Bucelli, Francesco Regazzoni, Luca Dede’, Alfio Quarteroni, Preserving the positivity of the deformation gradient determinant in intergrid interpolation by combining RBFs and SVD: Application to cardiac electromechanics, 2023, 00457825, 116292, 10.1016/j.cma.2023.116292 | |
33. | Ludovica Cicci, Stefania Fresca, Andrea Manzoni, Alfio Quarteroni, Efficient approximation of cardiac mechanics through reduced‐order modeling with deep learning‐based operator approximation, 2023, 2040-7939, 10.1002/cnm.3783 | |
34. | Benjamin Maier, Dominik Göddeke, Felix Huber, Thomas Klotz, Oliver Röhrle, Miriam Schulte, OpenDiHu: An efficient and scalable framework for biophysical simulations of the neuromuscular system, 2024, 18777503, 102291, 10.1016/j.jocs.2024.102291 | |
35. | Barış Cansız, Michael Kaliske, A comparative study of fully implicit staggered and monolithic solution methods. Part II: Coupled excitation–contraction equations of cardiac electromechanics, 2024, 442, 03770427, 115704, 10.1016/j.cam.2023.115704 | |
36. | Elena Zappon, Matteo Salvador, Roberto Piersanti, Francesco Regazzoni, Luca Dede’, Alfio Quarteroni, An integrated heart–torso electromechanical model for the simulation of electrophysiological outputs accounting for myocardial deformation, 2024, 427, 00457825, 117077, 10.1016/j.cma.2024.117077 | |
37. | Giorgos Troulliotis, Alison Duncan, Xiao Yun Xu, Alessandro Gandaglia, Fillipo Naso, Hendrik Versteeg, Saeed Mirsadraee, Sotiris Korossis, Effect of Excitation Sequence of Myocardial Contraction on the Mechanical Response of the Left Ventricle, 2024, 13504533, 104255, 10.1016/j.medengphy.2024.104255 | |
38. | Francesco Regazzoni, Stabilization of loosely coupled schemes for 0D–3D fluid–structure interaction problems with application to cardiovascular modelling, 2025, 157, 0029-599X, 249, 10.1007/s00211-025-01452-z | |
39. | Jiří Vaverka, Jiří Burša, A modification of Holzapfel–Ogden hyperelastic model of myocardium better describing its passive mechanical behavior, 2025, 111, 09977538, 105586, 10.1016/j.euromechsol.2025.105586 |
Geometry | Mesh | # Vertices | # Elements | |
Idealized(Fig. 6) | (a) | 4.6 mm | 1'827 | 3'010 |
(b) | 2.3 mm | 11'658 | 12'040 | |
(c) | 1.2 mm | 81'335 | 416'000 | |
Patient-specific(Fig. 7) | (d) | 0.8 mm | 126'031 | 637'379 |
Geometry | Mesh | # Vertices | # Elements | |
Idealized(Fig. 6) | (a) | 4.6 mm | 1'827 | 3'010 |
(b) | 2.3 mm | 11'658 | 12'040 | |
(c) | 1.2 mm | 81'335 | 416'000 | |
Patient-specific(Fig. 7) | (d) | 0.8 mm | 126'031 | 637'379 |
Geometry | Mesh | # CPUs | ||
Idealized | (a) | 1 | 14'616 | 1 |
(b) | 1 | 93'264 | 6 | |
(c) | 1 | 650'680 | 45 | |
(a) | 2 | 93'264 | 6 | |
Patient-specific | (d) | 1 | 1'008'248 | 72 |
| | | | | | | | | | | | | |
| | | | | | | | | | | 12 | 500 | |
| | | | | | | | | | | | | |
50 | 10 | 1500 | 0 | | 5 | 1 | 0 | 0 | | | | |
Geometry | Mesh | |||||||||
Idealized | (a) | 1 | 3.02 | 2.31 | 2.55 | 2.15 | 0.023 | 0.006 | 5h 41m | |
(b) | 1 | 3.37 | 2.14 | 3.44 | 6.67 | 0.163 | 0.019 | 18h 21m | ||
(c) | 1 | 4.23 | 2.49 | 3.26 | 9.19 | 0.632 | 0.040 | 46h 42m | ||
(a) | 2 | 3.29 | 1.68 | 2.11 | 6.47 | 0.178 | 0.025 | 30h 18m | ||
Patient-specific | (d) | 1 | 4.44 | 1.96 | 4.51 | 22.33 | 1.692 | 0.076 | 68h 37m |
Geometry | Mesh | # Vertices | # Elements | |
Idealized(Fig. 6) | (a) | 4.6 mm | 1'827 | 3'010 |
(b) | 2.3 mm | 11'658 | 12'040 | |
(c) | 1.2 mm | 81'335 | 416'000 | |
Patient-specific(Fig. 7) | (d) | 0.8 mm | 126'031 | 637'379 |
Geometry | Mesh | # Vertices | # Elements | |
Idealized(Fig. 6) | (a) | 4.6 mm | 1'827 | 3'010 |
(b) | 2.3 mm | 11'658 | 12'040 | |
(c) | 1.2 mm | 81'335 | 416'000 | |
Patient-specific(Fig. 7) | (d) | 0.8 mm | 126'031 | 637'379 |
Geometry | Mesh | # CPUs | ||
Idealized | (a) | 1 | 14'616 | 1 |
(b) | 1 | 93'264 | 6 | |
(c) | 1 | 650'680 | 45 | |
(a) | 2 | 93'264 | 6 | |
Patient-specific | (d) | 1 | 1'008'248 | 72 |
| | | | | | | | | | | | | |
| | | | | | | | | | | 12 | 500 | |
| | | | | | | | | | | | | |
50 | 10 | 1500 | 0 | | 5 | 1 | 0 | 0 | | | | |
Geometry | Mesh | |||||||||
Idealized | (a) | 1 | 3.02 | 2.31 | 2.55 | 2.15 | 0.023 | 0.006 | 5h 41m | |
(b) | 1 | 3.37 | 2.14 | 3.44 | 6.67 | 0.163 | 0.019 | 18h 21m | ||
(c) | 1 | 4.23 | 2.49 | 3.26 | 9.19 | 0.632 | 0.040 | 46h 42m | ||
(a) | 2 | 3.29 | 1.68 | 2.11 | 6.47 | 0.178 | 0.025 | 30h 18m | ||
Patient-specific | (d) | 1 | 4.44 | 1.96 | 4.51 | 22.33 | 1.692 | 0.076 | 68h 37m |