An integrative smoothed particle hydrodynamics framework for modeling cardiac function
Chi Zhang, Jianhang Wang, Massoud Rezavand, Dong Wu, Xiangyu Hu
AAn integrative smoothed particle hydrodynamicsframework for modeling cardiac function
Chi Zhang, Jianhang Wang, Massoud Rezavand, Dong Wu, Xiangyu Hu ∗ Department of Mechanical Engineering, Technical University of Munich, 85748 Garching,Germany
Abstract
Mathematical modeling of cardiac function can provide augmented simulation-based diagnosis tool for complementing and extending human understanding ofcardiac diseases which represent the most common cause of worldwide death. Asthe realistic starting-point for developing an unified meshless approach for totalheart modeling, herein we propose an integrative smoothed particle hydrody-namics (SPH) framework for addressing the simulation of the principle aspectsof cardiac function, including cardiac electrophysiology, passive mechanical re-sponse and electromechanical coupling. To that end, several algorithms, e.g.,splitting reaction-by-reaction method combined with quasi-steady-state (QSS)solver , anisotropic SPH-diffusion discretization and total Lagrangian SPH for-mulation, are introduced and exploited for dealing with the fundamental chal-lenges of developing integrative SPH framework for simulating cardiac function,namely, (i) the correct capturing of the stiff dynamics of the transmembranepotential and the gating variables , (ii) the stable predicting of the large defor-mations and the strongly anisotropic behavior of the myocardium, and (iii) theproper coupling of electrophysiology and tissue mechanics for electromechanicalfeedback. A set of numerical examples demonstrate the effectiveness and ro-bustness of the present SPH framework, and render it a potential and powerful ∗ Corresponding author.
Email addresses: [email protected] (Chi Zhang), [email protected] (JianhangWang), [email protected] (Massoud Rezavand), [email protected] (Dong Wu), [email protected] (Xiangyu Hu )
Preprint submitted to Comput. Methods Appl. Mech. Engrg. September 9, 2020 a r X i v : . [ c s . C E ] S e p lternative that can augment current lines of total cardiac modeling and clinicalapplications. Keywords:
Cardiac function, Electrophysiology, Electromechanics, Smoothedparticle hydrodynamics
1. Introduction
The heart is one of our not only most vital, but also most complex organs.The four chambers and four valves act precisely in concert to regulate the heart’sfilling, ejecting and overall pump function, by the interplay of electrical and me-chanical (including solid and fluid) dynamics. Cardiac diseases which effect theheart through complex mechanisms represent one of the most important cate-gory of problems in public health, effecting millions of people each year accordingto the reports of World Health Organization (WHO) [1]. Mathematical model-ing of the heart and its function can complement and expand our understandingof cardiac diseases and the clinical practice of cardiology [2]. In cardiac research,computational modeling and simulation have received tremendous efforts andare recognized as the community’s next microscope, only better [3]. However,the integrative model capable of simulating the fully coupled cardiac functionis still in its infancy. While the state-of-art computational models are able tosimulate the coupled electromechanics or the coupled fluid-solid dynamics, theyface serious difficulties on integrating all the three dynamics due to the conflictsbetween their limited modeling flexibility with respect to the complex physicalprocesses involved.To date, there are mainly two computational approaches [4] that developedfro integrative cardiac modeling, viz. the finite-element method (FEM) [5] andthe immersed-boundary method (IBM) [6]. However, the ability of FEM ishindered by the coupling between the solid and fluid mechanics. Typical dif-ficulties include the treatment of the convective terms, the incompressibilityconstraint and the updating of the mesh; especially when the opening and clos-ing of valves are taken into account. The IBM has been developed to compute2SI problems, which are the main difficulty for the FEM modeling. However,the fairly weak coupling formulation of distributing the forces computed on thedeformable Lagrangian mesh to the Eulerian mesh using a kernel function inIBM can lead to the Lagrangian-Eulerian mismatches on the kinematics. Suchdifficulty becomes very serious when the material properties and active stress arecomplex as the ones in cardiac myocardium. An alternative approach, the mesh-less methods, e. g. smoothed particle hydrodynamics (SPH) [7, 8, 9, 10], hasshown peculiar advantages in handling multi-physics and multi-scale problems[11, 12, 13, 14, 15, 16]. These advantages render the SPH method an interestingand potent alternative in the integrative simulation of cardiac function.As a realistic starting-point for developing an integrative meshless approachfor cardiac modeling, the main objective of this work is to present a SPH frame-work for simulating the fundamental and indispensable components, e. g. thecardiac electrophysiology, passive mechanical response and electromechanicalcoupling (active mechanical response), of cardiac function. The cardiac electro-physiology describes the myocardial electrical activation sequence, based on theion currents and tissue conductivity [2]. The cardiac fibers contract due to thepropagation of electrical stimuli initiated in the sinoatrial node. This electricalstimuli produces a sharp rise (depolarization) followed by a sudden fall (repo-larization) of the transmembrane potential. This phenomenon can be mathe-matically modeled by means of a reaction-diffusion equation where the sourceterm encapsulates the cellular ion exchange. In this paper, we consider thesimple monodomain approach derived with the assumption of equal anisotropicconductivities in the intra- and extra-cellular compartments. The monodomainapproach has been widely used for three-dimensional simulations consideringionic models ranging from simple FitzHugh-Nagumo variants [17, 18] to themore complex Luo-Rudy model [19]. Regarding the cardiac passive mechani-cal response, finite elasticity models are needed to describe cardiac contractionand relaxation due to the fact that the cardiac cells change in length by up to20 −
30% during a physiological contraction. Furthermore, the suitable elas-ticity model should replicate the anisotropic passive behavior determined via a3et of collagen fibers and sheets duo to the extremely complex and heteroge-neous material property of cardiac myocardium. In this work, we consider theclassic invariant-based presentation of the strain energy proposed by Holzapfeland Ogden [20] for the characterization of the passive mechanical response ofthe cardiac myocardium. The Holzapfel-Ogden model in which the cardiac my-ocardium is treated as a non-homogeneous, thick-walled and nonlinearly elasticmaterial is a structural based model that accounts for the muscle fiber direc-tion and the myocyte sheet structure. The electromechanical coupling can bephenomenologically described by means of activation models. In general, twoapproaches, namely, active stress [21] and active strain [22], can be followed forthe definition of activation models. Here, we consider the active stress [21] ap-proach with an evolution equation for active cardiomyocite concentration stress[23].As the first attempt towards a integrative meshless model for cardiac mod-eling, the proposed SPH framework should accurately characterize its interest-ing critical aspects, including electrophysiology, passive and active mechanicalresponses. At first, we adopt an operator splitting scheme for the reaction-diffusion equation to split reaction and diffusion to ensure numerical stabilityand accuracy. This consideration also leads to a much larger time step size com-pared to the simple forward Euler method [2]. Then, we introduce a splittingreaction-by-reaction method [24] combined with quasi-steady-state (QSS) solverto capture the stiff dynamics of the transmembrane potential and the gatingvariables of the ionic model governed by nonlinear ordinary differential equations(ODEs). Furthermore, an anisotropic SPH discretization for diffusion equationderived by Tran-Duc et al. [25] is modified by introducing a linear operator toimprove the computational efficiency and using a correction kernel matrix toimprove the numerical accuracy. The total Lagrangian SPH formulation is em-ployed to predict the large deformations and the strongly anisotropic behavior ofthe myocardium. Ultimately, the proposed SPH discretization for monodomainequation is integrated to predict the active response of myocardium by imple-menting the active stress approach [21]. A comprehensive set of numerical ex-4mples, viz. benchmarks on iso- and aniso-tropic diffusion, the transmembranepotential propagation in the manner of free-pulse and spiral wave, passive andactive responses of myocardium, and electrophysiology and electromechanics ina generic biventricular heart, are computed to demonstrate the accuracy, ro-bustness and feasibility of the proposed SPH framework. Base on the presentdevelopments and the previous achievements of the SPH method [14, 10, 12],the proposed framework will shed light on the multi-physics and multi-scaletotal cardiac modeling, in particular with regards to the fluid-electro-structureinteractions. For a better comparison and future openings for in-depth studies,all the computational codes and data-sets accompanying this work are availableonline at https://github.com/Xiangyu-Hu/SPHinXsys .This manuscript is organized as follows. Section 2 introduces the basic prin-ciples of the kinetics and the governing equations describing the evolution ofthe transmembrane potential and the mechanical response. Section 3 presentsthe constitutive laws with respect to the monodomain approach, the passiveand the active electromechanical responses. In Section 4, the proposed SPHframework is fully described. A comprehensive set of examples are included inSection 5 and the concluding remarks and a summary of the key contributionsof this paper are given in Section 6.
2. Kinematics and governing equations
To characterize the kinematics of the finite deformation, the deformationmap ϕ which maps a material point r from the initial reference configurationΩ ⊂ R d to the point r = ϕ (cid:0) r , t (cid:1) in the deformed configuration Ω = ϕ (cid:0) Ω (cid:1) isintroduced. In this work, we will use superscript ( • ) to denote the quantitiesin the initial reference configuration. The deformation tensor F can then bedefined by its derivative with respect to the initial reference configuration as F = ∇ ϕ = ∂ϕ∂ r = ∂ r ∂ r . (1)5ote that the deformation tensor F can also be calculated from the displacement u = r − r through F = ∇ u − I , (2)where I represents the unit matrix. For an incompressible material, we have theconstraint J = det ( F ) ≡ . (3)Associated with F are the right and left Cauchy-Green deformation tensorsdefined by C = F T · F and B = F · F T , (4)respectively. Then, four typical invariants of C (and also of B ) can be definedas I I = tr ( C ) , I ff = f · (cid:0) C f (cid:1) , I ss = s · (cid:0) C s (cid:1) , I fs = f · (cid:0) C s (cid:1) , (5)where f and s are the undeformed myocardial fiber and sheet unit direction,respectively. Here, I I is the first principal invariant, structure-based invariants I ff and I ss are the isochoric fiber and sheet stretch squared as the squaredlengths of the deformed fiber and sheet vectors, i. e. f = F f and s = F s , while I fs indicates the fiber-sheet shear [20].We consider a coupled system of partial differential equations (PDEs) govern-ing the motion of the material point r and the evolution of the transmembranepotential V m . The time dependent evolution of the transmembrane potential inLagrangian framework is characterized by the normalized monodomain equation[2] C m d V m d t = ∇ · ( D ∇ V m ) + I ion in Ω × [0 , T ] , (6)where C m is the capacitance of the cell membrane and I ion the ionic current.Note that the conductivity tensor is defined by D = d iso I + d ani f ⊗ f with d iso denoting the isotropic contribution and d ani the anisotropic contributionto account for faster conductivity along fiber direction f .6n a total Lagrangian framework, the conservation of the mass and linearmomentum corresponding to the cardiac mechanics can be expressed as ρ = ρ J ρ v d t = ∇ · P T Ω × [0 , T ] , (7)where ρ is the density and P the first Piola-Kirchhoff stress tensor and P = FS with S denoting the second Piola-Kirchhoff stress tensor. Note that the bodyforce is neglected in Eq. (7).
3. Constitutive equations
To close the systems of Eqs.(6) and (7), we specify herein the constitutivelaws for the ionic current I ion and the first Piola-Kirchhoff stress P . To close the monodomain equation Eq. (6), a model for the ionic current isrequired. Following Refs [17, 19], we consider the so-called reduced-ionic modelin which I ion ( V m , w ) is a function of the trasnmembrane potential V m and thegating variable w which represents the percentage of the open channels perunit area of the membrane. The most widely used reduced-ionic model is theFitzhugh-Nagumo model [17] and the variant Aliev-Panfilow model [18] whichonly have two currents, viz. inward and outward, one gating variable and noexplicit ionic concentration variables.The Fitzhugh-Nagumo model reads [17] I ion ( V m , w ) = − V m ( V m − a )( V m − − w ˙ w = g ( V m , w ) = (cid:15) ( βV m − γw − σ ) , (8)where (cid:15) , β , γ and σ are suitable constant parameters will be given specifically.As a variant of Fitzhugh-Nagumo model, the Aliev-Panfilow model [18] hasbeen successfully implemented in previous simulations of ventricular fibrillationin real geometries [26] and it is particularly suitable for applications where7lectrical activity of the heart is the main interest. The Aliev-Panfilow modelreads I ion ( V m , w ) = − kV m ( V m − a )( V m − − wV m ˙ w = g ( V m , w ) = (cid:15) ( V m , w )( − kV m ( V m − b − − w ) , (9)where (cid:15) ( V m , w ) = (cid:15) + µ w/ ( µ + V m ) and k , a , b , (cid:15) , µ and µ are suitableconstant parameters to be fixed later.Note that both Fitzhugh-Nagumo and Aliev-Panfilow models involve di-mensionless variable V m , w and t . The actual transmembrane potential E indimension mV and time T in dimension ms can be obtained through [18] E = 100 V m − T = 12 . t . (10) Following the work of Nash and Panfilov [21], we couple the stress tensorwith the transmembrane potential V m through the active stress approach whichdecomposes the first Piola-Kirchhoff stress P into passive and active parts P = P passive + P active . (11)Here, the passive component P passive describes the stress required to obtain agiven deformation of the passive myocardium, and an active component P active denotes the tension generated by the depolarization of the propagating trans-membrane potential.For the passive mechanical response, we consider the Holzapfel-Odgen modelwhich proposed the following strain energy function, considering different contri-butions and taking the anisotropic nature of the myocardium into account. Toensure that the stress vanishes in the reference configuration and encompassesthe finite extensibility, we modify the strain-energy function as W = a b exp [ b ( I − − a ln J + λ J ) + (cid:88) i = f,s a i b i { exp (cid:104) b i ( I ii − (cid:105) − } + a fs b fs { exp (cid:2) b fs I fs (cid:3) − } , (12)8here a , b , a f , b f , a s , b s , a fs and b fs are eight positive material constants, withthe a parameters having dimension of stress and b parameters being dimension-less. Here, the second Piola-Kirchhoff stress S being defined by S = 2 ∂ W ∂ C − p C − = 2 (cid:88) j ∂ W ∂ I j ∂ I j ∂ C − p C − j = I, f f, ss, f s, (13)where ∂ I ∂ C = I , ∂ I ff ∂ C = f ⊗ f , ∂ I ss ∂ C = f ⊗ f , ∂ I fs ∂ C = f ⊗ s + s ⊗ f , (14)and p is the Lagrange multiplier arising from the imposition of incompressibility.Substituting Eqs. (13) and (14) into Eq.(12) the second Piola-Kirchhoff stressis given as S = a exp [ b ( I I − { λ ln J − a } C − + 2 a f ( I f −
1) exp (cid:104) b f ( I f − (cid:105) f ⊗ f + 2 a s ( I s −
1) exp (cid:104) b s ( I s − (cid:105) s ⊗ s + a f s I fs exp (cid:104) b f s ( I fs ) (cid:105) fs ⊗ fs . (15)The cardiac electrical activation stem from two processes: the generation ofionic currents which produces the transmembrane potential at the microscopicscales and the traveling of the transmembrane potential from cell to cell at themacroscopic scales. The propagation of the transmembrane potential can bedescribed by means of PDEs, suitably coupled with ODEs governing the ioniccurrents. In particular, a monodomain equation can be defined with the contin-uum assumption of the coexistence of extra- and intra-cellular information atevery point. Following the active stress approach proposed by Nash and Pan-filov [21], the active component provides the internal active contraction stressby P active = T a F f ⊗ f , (16)where T a represents the active magnitude of the stress and its evolution is givenby an ODE as ˙ T a = (cid:15) ( V m ) [ k a ( V m − V r ) − T a ] , (17)9here parameters k a and V r control the maximum active force, the restingtransmembrane potential and the activation function [23] (cid:15) ( V m ) = (cid:15) + ( (cid:15) ∞ − (cid:15) −∞ ) exp {− exp (cid:2) − ξ (cid:0) V m − V m (cid:1)(cid:3) } . (18)Here, the limiting values (cid:15) −∞ at V m → −∞ and (cid:15) ∞ at V m → ∞ , the phase shift V m and the transition slope ξ will ensure a smooth activation of the muscletraction.
4. SPH method for cardiac eletrophysiology and electromechanics
In this section, the proposed SPH method for cardiac eletrophysiology, pas-sive mechanical response and the electromemchanical coupling is presented.
Before moving on to the SPH discretization, we first briefly summarize thetheory and fundamentals of the SPH method. For more details the readers arereferred to the comprehensive review in Ref. [27].By introducing a Dirac delta function δ ( r − ´ r ) around r , a continuousfunction f ( r ) and its approximation, i.e., the smoothing kernel function W ( r − ´ r , h ) with smoothing length h defining the support domain, has the relation f ( r ) = (cid:90) Ω f (´ r ) δ ( r − ´ r ) d ´ r ≈ (cid:90) Ω f (´ r ) W ( r − ´ r , h ) d ´ r , (19)where Ω denotes the volume of the integral domain. Here, the introduction ofsmoothing kernel function [9, 28] establishes a discrete model due to the finitesize of the smoothed length h . From Eq. (19) the gradient of function f can beapproximated by ∇ f ( r ) ≈ (cid:90) Ω ∇ f (´ r ) W ( r − ´ r , h ) dV (´ r ) . (20)Integrating by parts of Eq. ((20)) and applying Gauss theorem yields ∇ f ( r ) ≈ (cid:90) ∂ Ω f (´ r ) W ( r − ´ r , h ) n dS (´ r ) − (cid:90) Ω f (´ r ) ∇ W ( r − ´ r , h ) dV (´ r ) . (21)10f the computational domain is discretized by a set of particles, the gradient of f can be approximated as in SPH form as the first term vanishes due to compactsupport of the kernel function in the right hand side of Eq. (21) ∇ f ( r ) − N (cid:88) j =1 m j ρ j f ( r j ) ∇ i W ( r i − r j , h ) . (22)Note that m i /ρ i is defined to express the differential volume element dV i . As mentioned in Section 3.1, the monodomain equation consists of a coupledsystem of PDE and ODE. The former governs the diffusion of the transmem-brane potential and the latter the reactive kinetics of the gating variable. Inthis paper, we employ the operator splitting method [2] which results in a PDEof anisotropic diffusion C m d V m d t = ∇ · ( D ∇ V m ) , (23)and two ODEs C m d V m d t = I ion ( V m , w ) d w d t = g ( V m , w ) , (24)where I ion ( V m , w ) and g ( V m , w ) are defined by FitzHugh-Nagumo Eq. (8) orAliev-Panfilow model Eq. (9). Different from the previous strategies for the discretization of diffusion equa-tion [29, 30], we employ and modify the anisotropic SPH dicretization proposedby Tran-Duc et al. [25]. Following Ref. [25], the diffusion tensor D is consideredto be a symmetric positive-definite matrix and can be decomposed by Choleskydecomposition as D = LL T (25)where L is a lower triangular matrix with real and positive diagonal entriesand L T denotes the transpose of L . The diffusion operator in Eq. (23) can be11ewritten to isotropic form by ∇ · ( D ∇ ) = ∇ · ( LL T ∇ ) = ( L T ∇ ) T · ( L T ∇ ) = (cid:101) ∇ , (26)where (cid:101) ∇ = L T ∇ . Then, the new isotropic diffusion operator is approximatedby the following kernel integral with neglecting the high-order term (cid:101) ∇ · ( (cid:101) ∇ )Φ = 2 (cid:90) Ω Φ( (cid:101) r ) − Φ( (cid:101) r (cid:48) ) | (cid:101) r − (cid:101) r (cid:48) | ∂W (cid:16)(cid:101) r − (cid:101) r (cid:48) , (cid:101) h (cid:17) ∂ | (cid:101) r − (cid:101) r (cid:48) | d (cid:101) r , (27)where (cid:101) r = L − r and (cid:101) h = L − h . Upon the coordinate transformation, thekernel gradient can be rewritten as ∂W (cid:16)(cid:101) r − (cid:101) r (cid:48) , (cid:101) h (cid:17) ∂ ( (cid:101) r − (cid:101) r (cid:48) ) = 1 | L − || L − e (cid:101) rr | ∂W (cid:16) r − r (cid:48) (cid:17) ∂ | r − r (cid:48) | (28)with e r (cid:48) r = r (cid:48) − r | r (cid:48) − r | . At this stage, Eq. (26) can be discretized in SPH form as (cid:101) ∇ Φ ≈ N (cid:88) j m j ρ j (cid:18) Φ( r i ) − Φ( r j ) (cid:19) L − ij e ij ) r ij ∂W ij ∂r ij , (29)where e ij = r ij r ij , D ij = L ij L Tij and D ij = D i D j D i + D j , which ensure the antisymmetricproperty of the physical diffusion phenomenon. Note that Eq. (29) is excessivecomputational expensive due to the fact that one time of Cholesky decomposi-tion and the corresponding matrix inverse is required for each pair of particleinteraction. To optimize the computational efficiency, we modify Eq. (29) byreplacing the term L − ij with its linear approximation (cid:101) L ij given by (cid:101) L ij = (cid:101) L i (cid:101) L j (cid:101) L i + (cid:101) L j (30)where (cid:101) L i is defined as (cid:101) D i = (cid:16)(cid:101) L − i (cid:17) (cid:16)(cid:101) L − i (cid:17) T . (31)In this case, the Cholesky decomposition and the corresponding matrix inverseare computed once for each particle before the simulation. Also note that Eq.(29) can be further improved by introducing a kernel correction matrix to im-prove the numerical accuracy which will be detailed in the following section.12 .2.2. Reaction-by-reaction splitting The system of ODEs defined by Eq. (24) are generally stiff, therefore numer-ical instability occurs where the integration time step is not sufficiently small.In this work, we employ a reaction-by-reaction splitting method proposed byWang et al. [24]. The multi-reaction system can be decoupled, e.g. second-order accurate Strange splitting, as R (∆ t ) = R ( ∆ t ) V ◦ R ( ∆ t ) w ◦ R ( ∆ t ) w ◦ R ( ∆ t ) V , (32)where the ◦ symbol separates each reaction and indicates that the operator R (∆ t ) V is applied after R (∆ t ) w . Note that the reaction-by-reaction splitting methodologycan be extended to more complex ionic models, e.g. the Tusscher-Panfilov model[31].Following Ref. [24], we rewrite an ODE in Eq. (24) in the following formd y d t = q ( y, t ) − p ( y, t ) y, (33)where q ( y, t ) is the production rate and p ( y, t ) y is the loss rate [24]. The generalform of Eq. (33), where the analytical solution is not explicitly known or difficultto derive, can be solved by using the quasi-steady-state (QSS) method for anapproximate solution as y n +1 = y n e − p ( y n ,t )∆ t + q ( y n , t ) p ( y n , t ) (cid:16) − e − p ( y n ,t )∆ t (cid:17) . (34)Note that QSS-based method is unconditionally stable due to the analytic form,and thus a larger time step is allowed for the splitting method, leading to a highercomputational efficiency. The elastic response of the soft myocardium is highly nonlinear and theirdeformation under working load are intrinsically large, therefore a robust nu-merical method is required. In this work, we adopt the total Lagrangian SPHformulation, i.e., the initial reference configuration is used for finding the neigh-boring particles and the set of neighboring particles is not altered, to ensure thefirst-order consistency and eliminates the tensile-instability,13ollowing the work of Vignjevic et al. [32], a correction matrix B is firstintroduced as B i = (cid:88) j V j (cid:0) r j − r i (cid:1) ⊗ ∇ i W ij − , (35)where ∇ i W ij = ∂W (cid:0) | r ij | , h (cid:1) ∂ | r ij | e ij (36)stands for the gradient of the kernel function evaluated at the initial referenceconfiguration. Again, the correction matrix is computed in the initial config-uration and therefore, it is calculated only once before the simulation. UsingEqs. (22) and (35), the linear momentum conservation equation, Eq.(7), can bediscretized in the following form¨ r i = 2 m i (cid:88) j V i V j ˜ P ij ∇ i W ij , (37)where the inter-particle averaged first Piola-Kirchhoff stress ˜ P is given as˜ P ij = 12 (cid:0) P i B i + P j B j (cid:1) . (38)Here, the first Piola-Kirchhoff stress tensor is computed with the constitutivelaw where the deformation tensor F is computed by F i = (cid:88) j V j ( u j − u i ) ⊗ ∇ i W ij B i − I . (39)It worth noting that, we apply the renormalization kernel correction to im-prove the accuracy and consistency of Eq. 29 as C m ˙ V m = 2 (cid:88) j m j ρ j (cid:18) B i V m ( r i ) − B j V m ( r j ) (cid:19) · e ij · e ij ( (cid:101) L ij e ij ) r ij ∂W ij ∂r ij . (40) Here we describe the details of the implementation of the proposed SPHframework for integrating the monodomain equation with mechanical responseof myocardium. To maintain the numerical stability, the time step size forsolving the monodomain equation is restricted by the diffusion coefficient∆ t p = 0 . (cid:18) h d | D | (cid:19) , (41)14here d is the dimension number and | D | the trace of the diffusion tensor. Forthe passive elastic response, he Courant-Friedichs-Levy (CFL) condition is givenas ∆ t m = 0 . (cid:32) hc + | v | max , (cid:115) h | d v d t | max (cid:33) . (42)The final time step size is chosen from∆ t = min (∆ t p , ∆ t m ) . (43)We denote the values at the beginning of a time step by the superscript n ,at the mid-point by n + and eventually at the end of the time-step by n + 1.Following the splitting method, the transmembrane potential V m and the gatevariable w are first updated in sequence for a half time step as V n + m = V nm + ∆ t p (cid:0) d V m d t (cid:1) reaction w n + = w n + ∆ t p (cid:0) d w d t (cid:1) reaction . (44)Here, the diffusive operator is applied and the transmembrane potential V m isupdated for a time step V ∗ m = V n + m + ∆ t p (cid:18) d V m d t (cid:19) diffusion . (45)Then the transmembrane potential V m and the gate variable w are updated ininverse sequence for another half time step as w n + = w n + ∆ t p (cid:0) d w d t (cid:1) reaction V n + m = V ∗ m + ∆ t p (cid:0) d V m d t (cid:1) reaction . (46)At this point, the active cardiomyocite contraction stress T a is updates forone time step if active response is taken into consideration. Following Ref.[14], a position-based Verlet scheme is applied for the time integration of themechanical response. At first, the deformation tensor, density and particleposition are updated to the midpoint as F n + = F n + ∆ t d F d t ρ n + = ρ J r n + = r n + ∆ t v n . (47)15hen the velocity is updated by v n +1 = v n + ∆ t d v dt . (48)Finally, the deformation tensor and position of solid particles are updated tothe new time step of the solid structure with F n +1 = F n + + ∆ t d F d t ρ n +1 = ρ J r n +1 = r n + + ∆ t v n +1 . (49)An overview of the complete solution strategy is presented in Algorithm 1 inAppendix A.
5. Numerical experiments
This section is devoted to present a comprehensive set of numerical exam-ples for validating the integrative SPH framework for the simulation of cardiacfunction with respect to electrophysioldoy, passive mechanical response and theelectromechanical coupling. We start with the benchmarks on both iso- andaniso-tropic diffusion process. We then validate the present method for cardiacelectrophysiology by solving the monodomain equation on regular and irregularcomputational domain with iso- and aniso-tropic diffusion coefficients. Then theaccuracy, robustness and applicability of the total Lagrangian SPH method formodeling the passive and active mechanical responses of myocardium are vali-dated. Having the validation studies presented, the excitation and excitation-contraction of three-dimensional generic biventricular heart are to show thepotential of the proposed SPH framework. In all the following simulations, the5 th -order Wendland smoothing kernel function [28] with a smoothing length of h = 1 . dp is employed, where dp denotes the initial particle spacing. Following Refs [33, 30], the first problem studied herein of the one-dimensionalisotropic diffusion is depicted in Figure 1 where a 0 . ×
1m rectangle is filled16 ! 𝑧 " 𝑧 C Figure 1: Schematic diagram for the one-dimensional diffusion problem with constant initialdistribution of the pollutant concentration. with water and a finite horizontal band of pollutant is located in the middleof the rectangle confined to the region of z ≤ z ≤ z . The initial pollutantconcentration is equal to C = 1kg · m − in the band and zero elsewhere, andthe diffusion coefficient is set as d iso = 1 . × − . According to Ref. [34], theanalytical solution is C ( z, t ) = C erfc (cid:16) z − z √ d iso t (cid:17) for z ≤ z C ( z, t ) = C erfc (cid:16) z − z √ d iso t (cid:17) for z > z . (50)where z = 0 . z = 0 .
45m and z = 0 . t = 1 . s is shown in Figure 2 (a). The bell shaped distribution of the concentration is inagreement with the analytical solution. It can also be observed that the presentresults converges with increasing spatial resolution.This problem is further considered by setting an exponential initial pollutantconcentration distribution as C ( z, t = 0) = exp (cid:18) − ( z − z ) d iso t (cid:19) , (51)where t = 1s and z = 0 . C ( z, t = 0) = C √ t + t exp (cid:18) − ( z − z ) d iso ( t + t ) (cid:19) . (52)Figure 2 (b) illustrates the comparison of the present predictions of the concen-tration distribution against the analytical solution at t = 1 . s . Again, a goodagreement is noted and the convergence of the concentration distributions withincreasing resolution is observed. 17 (m) C ( kg . m ) Analytical20 x 5040 x 10080 x 200 (a) constant initial concentration z(m) C ( kg . m ) Analytical20 x 5040 x 10080 x 200 (b) exponential intial concentrationFigure 2: Comparison of the concentration distribution between present results and the ana-lytical solutions for isotropic diffusion processes with constant and exponential initial concen-tration distributions.
In this section, we consider the anisotropic diffusion process from a con-taminant source in water. Following the work of Tran-Duc et al. [25], thecontaminant source is located in a two dimensional 200m × C ( z, t ) = 14 πt (cid:81) i =1 √ D ii (cid:89) i =1 exp (cid:20) − ( x i − x i, ) t D ii (cid:21) . (53)The initial condition for numerical solution is set at time t = 120s.In the first case, the anisotropic diffusion tensor is D = .
09 00 0 . (cid:0) m · s − (cid:1) . (54)Figure 3 shows the numerical and analytical distributions at time t = 1920 s . Itcan be observed that introducing the renormalized kernel correction can improvethe computational accuracy. In general, the concentration distributions are likeellipses with major axis in x -direction and minor axis in y -direction due to thefact that the diffusion rate in x -direction is larger than that in y -direction. Also,18he numerical solution is in agreement with the analytical one in both shape andvalue profiles. Figure 4 gives the present numerical concentration distributionsat horizontal cross section at x = 100m (Figure 4a) and vertical cross sectionat y = 100m (Figure 4b) and the corresponding comparison with analyticalsolutions. The present SPH approximated concentration profiles are in goodagreement with the analytical solution. Also, the present SPH results convergesto the analytical solution as the spatial resolution increases. Compared with theresults obtained by Tran-Duc et al. [25] with resolution 400 ×
400 (see Figure 3in their work), present results shows similar accuracy even as lower resolution200 ×
200 is used with the introduction of the renormalized kernel correction.Similar to Ref. [25], the present results reduce the anisotropy level of diffusionprocess and show a bit less anisotropic compared with the analytical one. Thisdiscrepancy induced by the isotropic property of kernel function in SPH whichaverages and smoothes the concentration function independent of direction.
Figure 3: Concentration distribution at 1920 s diffusion with the anisotropic diffusion tensor D . In the second test, a higher anisotropic ratio is considered by setting thediffusion tensor as D = . . (cid:0) m · s − (cid:1) . (55)Figure 5 shows the comparison between the simulated and analytical concentra-tion distributions at time t = 1920s. Again, the renormalized kernel correctionshows improved computational accuracy. As expected, the concentration distri-butions are also ellipses but with a higher ratio of major axis to minor comparedwith the results depicted in Figure 3. Again, the simulated distributions are in19 C Analytical50 x 50 (no kernel correction)50 x 50 (kernel correction)100 x 100 (kernel correction)200 x 200 (kernel correction) (a) at x = 100 m y C Analytical50 x 50 (no kernel correction)50 x 50 (kernel correction)100 x 100 (kernel correction)200 x 200 (kernel correction) (b) at y = 100 m Figure 4: Concentration distribution after 1800 s diffusion with the anisotropic diffusion tensor D : comparison with analytical solution by using three different spatial resolutions. consistent with the analytical one in both shape and value profiles. The presentnumerical concentration distributions at horizontal cross section at x = 100mand vertical cross section at y = 100m are given in Figs. 6a and 6b, respectively. Figure 5: Concentration distribution at 1920 s diffusion with the anisotropic diffusion tensor D . Following the work of Ratti and Verani [35], we consider a problem on thepropagation of transmembrane potential. It is assumed that the transmem-brane potential propagates in a two dimensional isotropic tissue in a squaredomain of (0 , and the nondimensional time interval is set as (0 , V m = exp (cid:104) − ( x − . + y . (cid:105) w = 0 . , (56)20 C Analytical50 x 50 (no kernel correction)50 x 50 (kernel correction)100 x 100 (kernel correction)200 x 200 (kernel correction) (a) at x = 100 m y C Analytical50 x 50 (no kernel correction)50 x 50 (kernel correction)100 x 100 (kernel correction)200 x 200 (kernel correction) (b) at y = 100 m Figure 6: Concentration distribution after 1800 s diffusion with the anisotropic diffusion tensor D : comparison with analytical solution by using three different spatial resolutions. and the nondimensional diffusion coefficients are d iso = 1 . d ani = 0 . Table 1: Parameters for the Aliev-Panfilow model. k a b (cid:15) µ µ . , . We now validate the SPH method in reproducing the spiral waves by solv-ing the monodomain equations with the FitzHugh-Nagumo model. The spiralwaves, which consists of complicated patterns of the transmembrane potentialalong with simple unidirectionally propagating pulses, are suitable choices for21 V m Figure 7: Time evolution of the transmembrane potential at point (0 . , .
7) in comparisonwith the results reported by Ratti and Verani [35]. validating the numerical solution of solving the reaction-diffusion equation. Inthis work, we consider both isotropic and anisotropic tissues in two dimensionalrectangular and circular geometries with the given parameters of the FitzHugh-Nagumo model in Table 2
Table 2: Parameters for the FitzHugh-Nagumo model. a (cid:15) β γ σ Following Wang et al. [36] and Liu et al. [37], the rectangular computationalregion is set as [0 , . × [0 , .
5] and the transmembrane potential and gatevariable are initialized by V m = . , < x ≤ .
25; 0 < y < . . , elsewhere , (57)and w = . , < x ≤ .
25; 1 . ≤ y < . . , . ≤ x < .
5; 0 < y < . , elsewhere . (58)22n the first test, we consider an isotropic tissue with the nondimensional dif-fusion D where d ios = 1 . × − and d ani = 0 .
0. Figure 8 (upper panel) showsa spiral wave of the stable rotation solution at five different time instants. Asobserved, the spiral wave generates a clockwise rotation curve in the rectangularregion. Note that the spiral wave profiles obtained by the present SPH methodis consistent with those obtained by Wang et al. [36] and Liu et al. [37] (seeFigure 1 (a) and (b) in Ref. [36]).Figure 8 (middle panel) shows the numerical results with an anisotropicdiffusion tensor D = . × −
00 2 . × − . (59)It can be observed that the spiral wave now propagates with elliptical patternsas also seen in Refs. [36, 37]. Again, a good agreement with those given in Refs.[36, 37] (see Figure 2 (a) and (b) in Ref. [36]) is noted.Another test with a larger diffusion ratio, given by D = . × −
00 1 . × − , (60)is reported in Figure 8 (lower panel). It can be found that the spiral wave nowhas a slightly smaller width compared with the one shown in Figure 8 (middlepanel). Also, the elliptical propagation shape of the spiral wave has a higherratio between the major and minor axes. In this part, the computational domain is changed to a nonuniform geometry,i.e. a circle of radius R = 1 .
25 and centered at r = ( R, R ). The transmembranepotential and gate variable are initialized by [36, 37] V m = . , R − (cid:113) R − ( R − y ) < x ≤ RR − (cid:113) R − ( R − x ) < y ≤ R . , elsewhere , (61)23 igure 8: Spiral waves of the FitzHugh-Nagumo model in rectangular computational domainwith isotropic diffusion tensor D (upper panel), D (middle panel) and D (lower panel) atfive time instants. and w = . , R − (cid:113) R − ( R − y ) < x < R + (cid:113) R − ( R − y ) R ≤ y < R + (cid:113) R − ( R − x ) . , elsewhere . (62)Figure 9 (upper panel) shows the contours of the stable rotating spiral wavewith isotropic diffusion tensor D given by previous section. As expected, thespiral wave generates a curve and rotates clockwise as reported in [36, 37].Again, a good agreement with those of Refs. [36, 37] is noted (see Figure 3 (a)and (b) in Ref. [36]).For the anisotropic diffusion tensor D given in Eq. (59), the transmember-ane potential propagation at four different time instants are shown in Figure9 (middle panel). Now, the spiral wave follows an elliptical pattern in the cir-24ular region. Figure 9 shows the propagation pattern of the spiral wave withanisotropic diffusion tensor D given in Eq. (60). With larger diffusion ratio,the propagation pattern of the spiral waves is effected and the width of the spiralwave is also slightly smaller than the one reported in Figure 9 (lower panel). Figure 9: Spiral waves of the FitzHugh-Nagumo model in circular computational domain withthe isotropic diffusion tensor D (upper panel), D (middle panel) and D (lower panel) atfive time instants. In this section, two benchmarks are investigated to validate the accuracy,robustness and applicability of the present SPH framework for modeling thepassive and active mechanical responses of the cardiac myocardium.
In this part, we consider the passive mechanical response of the myocardiumin the form of bending cantilever. Following Aguirre et al.[38], a three-dimensionalrubber-like cantilever whose bottom face is clamped to the ground and its25ody is allowed to bend freely by imposing an initial uniform velocity v =(5 √ , , T m · s − is considered (see Figure 10). For in-depth comparisons,both neo-Hookean and Holzapfel-Odgen material models are applied. For theneo-Hookean model, the strain-energy density function [39] is defined as W = µ tr ( E ) − µ ln J + λ J ) , (63)where λ and µ are Lam´ e parameters, K = λ + (2 µ/
3) is the bulk modulus and G = µ is the shear modulus. The relation between the two modulus is given by E = 2 G (1 + 2 ν ) = 3 K (1 − ν ) , (64)with E denoting the Young’s modulus and ν the Poisson ratio. Here, the Youngs’modulus is E = 1 . × Pa, Poisson ratio ν = 0 .
45 and density ρ = 1 . × kg · m − . For the Holzapfel-Odgen model, the material parameters are givenin Table 3 and the anisotropic terms are varying accordingly. Table 3: Parameters for the Holzapfel-Ogden constitution model (For the isotropic material,the anisotropic terms are set to zero). a = 5 .
860 MPa a f = ka a s = 0 a fs = 0 . b = 1 . b f = 0 b s = 0 b fs = 0Figure 11 shows the deformed configuration colored with von Mieses stresscontours. Compared with the results reported in Ref. [38] (see Figure 24 intheir work), good agreement in the deformation is observed. Also note thatboth material models predict almost the same deformed configurations. Quan-titative comparisons are given in Figure 12 which plots the time history of thevertical displacement of point S and a good agreement is noted. Figure 13shows the convergence study of this example with isotropic Holzapfel-Odgenmaterial model. The convergence of the solution is illustrated with increasedspatial resolution.We further demonstrate the applicability of the present method by studyingthis example considering the anisotropic Holzapfel-Odgen material model. For26 v xx zz yyss Figure 10: Passive response of a three-dimensional bending cantilever: initial configuration.Figure 11: Passive response of a three-dimensional bending cantilever: time evolution of thevon Mieses stress distribution in the deformed configuration for Neo-Hookean and Holzapfel-Ogden materials. The spatial particle discretization is h/dp = 12. ime(s) x Figure 12: Passive response of a three-dimensional bending cantilever: time history of thevertical position at node S . Results of isotropic Neo-Hookean and Holzapfel-Ogden materialsare compared with that of Aguirre et al. [38]. The spatial particle discretization is h/dp = 12. Time(s) x Figure 13: Passive response of a three-dimensional bending cantilever: time history of thevertical position at node S and the convergence study for the isotropic Holzapfel-Ogden ma-terial. the anisotropic material, we set he fibre and sheet directions f aligned with x and y directions, respectively. Three tests with different aniostropic ratios,viz. a f /a = 0 . a f /a = 0 . a f /a = 1 .
0, are studied. Figure 14 shows thedeformed configuration while Figure 15 plots the time history of the verticaldisplacement of point S . It can be observed that the deformation is reduced asthe aniostropic ratio increases. 28 igure 14: Passive response of a three-dimensional bending cantilever: time evolution of thevon Mieses stress distribution in the deformed configuration for anisotropic Holzapfel-Ogdenmaterial with an spatial particle discretization of h/dp = 12. The aniostropy ratio is set to(a) a f /a = 0 . a f /a = 0 . a f /a = 1 .
0. Note that the corresponding verticaldisplacement of point S is ploted in Figure 15. Following Garcia-Blanco et al. [40], we consider a unit cube of myocardiumwith an orthogonal material direction. The myocardium has the fiber and sheetdirections parallel to the global coordinates and the constitutive law describingthe passive response is the Holzapfel-Ogden model with the material parametersgiven in Table 4. To initiate the excitation-induced response [40], the transmem-brane potential is linearly distributed along the vertical direction with V m = 0and V m = 30 at bottom and top faces, respectively. For simplicity, the timevariation of transmembrane potential is neglected and an ad-hoc activation lawof active stress is given by T a = − . V m . (65)29 ime(s) x f / a = 0.1a f / a = 0.5a f / a = 1.0 Figure 15: Passive response of a three-dimensional bending cantilever: time history of thevertical position at node S . Study for anisotropic properties of the Holzapfel-Ogden materialmodel with three aniostropic ratios : a f /a = 0 . a f /a = 0 . a f /a = 1 .
0. The spatialparticle discretization is h/dp = 12.
Two different tests with iso- and aniso-tropic models are considered herein.Figure 16 shows the deformed configuration of the cubic myocardium. Com-pared with the results reported in Ref. [40] (see Figure 7 in their work), aqualitative good agreement is noted for the isotropic test. Furthermore, thepresent simulation shows that the displacement of the top face is 0 .
53, whichis in good agreement with that of 0 .
535 given in Ref. [40]. For the anisotropictest, the deformation is reduced due to the existence of the fiber and the sheet.
Table 4: Parameters for the Holzapfel-Ogden constitution model (For the isotropic material,the anisotropic terms are set to zero). a = 0 .
059 kPa a f = 18 .
472 kPa a s = 2 .
841 kPa a fs = 0 .
216 kPa b = 8 . b f = 16 . b s = 11 . b fs = 11 . To demonstrate the abilities of the present SPH framework in total cardiacsimulation, we consider the transmembrane potential propagation as free pulsestogether with scroll waves and the corresponding excitation-contraction in three-dimensional generic biventricular heart.30 igure 16: Active response of the unit cubic myocardium with both isotropic and anisotropicmaterial properties: contour plot of the transmembrane potential.
Following the work of [41], the inner surface of the left and right ventriclesof the generic biventricular heart are described by two ellipsoids x a lv + y b lv + z c lv = 1 , x a rv + y b rv + z c rv = 1 (66)where a lv = 45 mm , b lv = 54 mm , c lv = 24 mm and a rv = 18 mm , b rv = 58 mm , c rv = 18 mm . The ellipsoids are truncated from apex-to-centroid as shownin Figure 17 (a). We impose a wall thickness of 6 and 12 on the left andright ventricle, respectively. To discretize the generic biventricular, particlesare initialized through a relaxation-based algorithm, and the fiber and sheetdirections are computed approximately by a coupling level-set and rule-basedalgorithm. Before moving onto the simulation of biventricular heart, we introduce arelaxation-based technique to generate isotropic initial particle distribution. Acoupled level-set and rule-based algorithm is also introduced for fiber and sheetreconstruction.For solid dynamics, two approaches, viz, direct particle generation based ona lattice structure [42] and particle generation based on a volume element mesh[43], are commonly used in the SPH community. In the former approach, parti-31les are positioned directly on a cubic lattice and equispaced particle distributionis obtained. Accurate surface description, in particular complex geometries, re-quires a fine resolution in this approach thereby rendering it limited to rathersimple geometries [44]. The second approach convertes each volume element ofa tetra or hexahedral mesh into a particle. This approach shows advantages indescribing complex geometries, however, yields significantly non-uniform par-ticle distributions regarding the particles spacing and size. In this work, weintroduce an approach initialized from standard triangle language (STL) inputfiles, which uses the relaxation-based algorithm, proposed by Fu et al. [45] formesh generation, to generate the initial particle distribution of the biventric-ular heart. Following Ref. [45], a level-set field on a Cartesian backgroundmesh is required for particle relaxation. In the present work, the geometry isdescribed in the STL format as shown in Figure 17 and a passer is used forreading data from the STL files. Then the geometry surface is represented bythe zero level-set of a signed-distance function,Γ = { ( x, y, z ) | φ ( x, y, z ) = 0 } . (67)Here, the distance from a mesh point to the geometry surface is determined byfinding the nearest point on all vertices and a positive phase is defined if themesh point is located inside the object, otherwise a negative phase is marked.Then, the particle evolution is conducted for a number of steps following thestrategy proposed by Fu et al. [45] (see Section 5.2 in their work). Note thatin this work a constant particle smoothing length and constant backgroundpressure and density are used. Also note that the singularities are not takeninto consideration, i.e. the surface particles are only constrained on the geometrysurface. Figure 17 (b) shows the particle distribution for a biventricular heartafter 5000 steps of relaxation with a background pressure of p = 2 . ρ = 1 .
0. As expected, an isotropic particle configuration is obtainedand the geometry surface is reasonably well prescribed.Following the particle initialization, the fiber and sheet reconstructions areconducted. Assuming that the sheets are aligned with the transmural direction,32he sheet direction can be approximated directly from the level-set function s = sign( N , e y ) N , (68)where N is the normal direction obtained from N = ∇ φ |∇ φ | , (69)and e y is the normal vector parallel to the ventricular centerline, pointing fromapex to base. For each particle, the sheet direction is interpolated from the level-set field by using the trilinear interpolation. Following the work of Quarterioniet al. [46], the initial flat fiber direction of each particle can be defined by (cid:101) f = s × e y . (70)Then, the fiber direction f can be defined by rotating (cid:101) f with respect to the s axis according to the following rotation formula f = (cid:101) f cos ( θ ) + s × (cid:101) f sin ( θ ) + s (cid:16) s · (cid:101) f (cid:17) [1 − cos ( θ )] , (71)where the rotation angle θ is computed from θ = ( θ epi − θ endo ) ψ + θ endo . (72)Here θ epi = − o and θ endo = 80 o are the rotation angles at the epicardiumand endocardium, respectively. The pseudo-distance ψ is given by solving asimple Laplace equation of ψ by imposing the boundary condition ψ | epi = 1and ψ | endo = 0 [46]. Figure 17 shows the fiber direction of plane located at y = −
10 and y = −
30 in epicardium (c) and endocardium (d), respectively.
In this section, we consider the transmembrane potential propagation in themanner of a free pulse and a scroll wave for the generic biventricular heart withiso- and aniso-tropic material properties. The Aliev-Panfilow model is appliedfor monodomain equation with the constant parameters given in Table 5 and thediffusion coefficients are set as d iso = 1 . · ms − and d ani = 0 . · ms − .33 igure 17: Generic biventricluar heart: (a) visualization of an STL file, (b) particle distri-bution, (c) fibre orientations of epicardium and (d) endocardium at cross sections located at y = −
10 and y = − k a b (cid:15) µ µ S
1, is initiated by externally stim-ulating the particles located at the upper part of the septum (wall separatingthe ventricles) as indicated by the partially depolarized region at t = 0 . V m = 0 .
92. The stimulus generates the depolar-ization through the heart as shown in Figure 18. It can be observed that thetransmembrane potential propagates in the similar pattern for both iso- and34niso-tropic material model. For more comprehensive comparison, the trans-membrane potentials recorded at apex are plotted in Figure 19. As expected,the transmembrane potential propagates in anisotropic model faster than thatin isotropic model. For both iso- and aniso-tropic materials, the transmembranepotential profiles show a good agreement with the results reported in Ref. [18].
Figure 18: Generic biventricluar heart: transmembrane potential propagates in the heart inthe free-pulse pattern. The snapshots depict contours of the transmembrane potential V m . Time V m IsoAniso
Figure 19: Generic biventricluar heart: time evolution of the transmembrane potentialrecorded on the apex.
As mentioned in Section 5.4, the present method shows good accuracy inreproducing two-dimensional spiral waves. In this section, we will demonstratethe present method’s ability to reproduce the formation of scroll waves in a more35omplex biventricular heart. We consider two tests: one with single scroll waveand another with two scroll waves interacting with each other when propagating.To generate the scroll wave, the S1-S2 protocol, where a second broken stimulus(S2) is triggered during the repolarization phase of the S1 wave, is applied.For a single scroll wave, the S2 stimulus is initiated at a small region ofΩ = { ( x, y ) | ≤ x ≤ ∧ − ≤ y ≤ } (73)located at the anterior ventricular wall from time t = 105 to t = 105 . V m = 0 .
95. Figure 22 shows the formation and evolutionof the vortex wave re-entry in both iso- and aniso-tropic material models. It canbe observed that the combination of the complex biventricular geometry, thenon-symmetric perturbation and the inhomogeneous fiber and sheet orientationclearly triggers a chaotic non-stationary wave pattern with the center of thescroll moving in the septal basal region. Figure 21 gives the recorded profilesof the transmembrane potential at the apex. After the depolarization of the S1wave, self-oscillatory transmembrane potential is noted due to the propagationof scroll waves. Note that the self-oscillatory transmembrane potential showshigher frequency in anisotropic model compared to the isotropic model.
Figure 20: Generic biventricluar heart: transmembrane potential propagation on both iso-and aniso-tropic models. The snapshots depict contours of the transmembrane potential V m . For the two scroll waves, the initiated region of S2 is extended toΩ = { ( x, y ) | ≤ x ≤ ∧ − ≤ y ≤ } , (74)36 ime V m IsoAniso
Figure 21: Generic biventricluar heart: time evolution of the transmembrane potentialrecorded at the apex. and the initiated time is changed to the interval between t = 125 and t = 125 . Figure 22: Generic biventricluar heart: transmembrane potential propagation on both iso-and aniso-tropic models. The snapshots depict contours of the transmembrane potential V m . ime V m IsoAniso
Figure 23: Generic biventricluar heart: time evolution of transmembrane potential recordedat the apex.
In this section, we demonstrate that the basic feature of the cardiac func-tion can be captured by the present SPH framework in a reasonable mannerby modeling the excitation-contraction of the biventricular heart through elec-tromechanical coupling. Three different excitations, including not only the freepulse but also the scroll waves represented in Section 5.8.2, are considered. As amatter of fact, the scroll wave may correspond to pathological heart diseases, i.e.cardiac arrhythmias can be related to the presence of wavefront spirals whichlead to an irregular contraction of the cardiac muscle. Therefore, reproducingthe excitation-contraction under the scroll wave may extend human understand-ing of cardiac arrhythmias. For simplicity, the displacement degrees of freedomon the top base are constrained and the whole heart surface is assumed to beflux-free. To record the heart displacement, three nodes, namely A located at(0 , − , T mm, B at (0 , − , T mm and C at ( − , − , T mm, are used.Moreover, the constant parameters of Holzapfel-Ogden model are given in theTable 4 and the active contraction stress is T a = 0 .
15 kPa.In the first test, we consider excitation-contraction under the transmembranepotential propagation as a free-pulse. Figure 24 shows the resulting excitation-contraction of the heart with the transmembrane potential contours and thecorresponding cross sections. It can be observed that excitation-contraction38ives rise to the upward motion of the apex as the depolarization front travel-ing through the heart. Also, the apex’s upward motion is accompanied by thephysiologically observed wall thickening and the overall torsional motion of theheart as shown in the cross sections of Figure 24. This physiologically activeresponse through the non-uniform contraction of myofibers is due to the inho-mogeneous myocyte orientation distribution incorporated with the anisotropicmaterial model. Figure 25 shows the time evolution of the x , y and z compo-nents of the displacement at points A, B and C, respectively. At the end of thedepolarization process, the reference configuration is recovered. Figure 24: Generic biventricular heart: coupled excitation-contraction induced by the trans-membrane potential as a free pulse. Snapshots of the deformed body depict the transmem-brane potential contours at different stages of the depolarization and the corresponding crosssections.
In the second test, we consider the excitation-contraction corresponding tothe single scroll wave. The resulting excitation-contraction of the heart is shownin Figure 26. As the propagation of the scroll wave, the myocytes show oscilla-tory excitation-contraction and the heart is under contracted state during thesimulation. Figure 27 shows the time evolution of the x , y and z components ofthe displacement at points A, B and C, respectively. Different from the previousresults for a free pulse, the motion of the observed points are highly oscillatoryand non-recoverable.In the third test, the two scroll wave excitation-induced contraction is in-39 ime D i s p l ace m e n t x disp at Ay disp at Bz disp at C Figure 25: Generic biventricular heart: coupled excitation-contraction induced by the trans-membrane potential as a free pulse. The time histories of displacement at nodes A, B andC.Figure 26: Generic biventricular heart: coupled excitation-contraction induced by the trans-membrane potential as a single scroll wave. Snapshots of the deformed body depict the trans-membrane potential contours at different stages of the depolarization and the correspondingcross sections. vestigated and the results are shown in Figure 28 with the transmembranepotential field and the corresponding cross sections. Also here, the myocyteshows oscillatory excitation-contraction resulting the contracted state of theheart. As the propagation of the scroll wave, the myocytes show oscillatoryexcitation-contraction and the heart is under contracted state during the simu-lation. Figure 29 shows the time evolution of the x , y and z components of the40 ime D i s p l ace m e n t x disp at Ay disp at Bz disp at C Figure 27: Generic biventricluar heart: coupled excitation-contraction induced by the trans-membrane potential as a single scroll wave. The time histories of displacement at nodes A, Band C. displacement at points A, B and C, respectively. Different from the previousresults for a single scroll wave, the amplitude of the oscillation is decreased whilethe frequency is slightly increased.
Figure 28: Generic biventricluar heart: coupled excitation-contraction induced by the trans-membrane potentials double scroll waves. Snapshots of the deformed body depict the trans-membrane potential contours at different stages of the depolarization and the correspondingcross sections. ime D i s p l ace m e n t x disp at Ay disp at Bz disp at C Figure 29: Generic biventricluar heart: coupled excitation-contraction induced by the trans-membrane potential as two scroll waves. The time histories of displacement at nodes A, Band C.
6. Concluding remarks
As a realistic starting-point for developing a unified SPH approach for simu-lating total heart function, this paper address the numerical modeling of manychallenging aspects of heart function, including cardiac electrophysiology, pas-sive mechanical response and the electromechanical feedback. For electrophys-iology, we solve the monodomain equation by introducing a splitting reaction-by-reaction method combined with quasi-steady-state (QSS) solver to capturethe stiff waves. For stable prediction of the large deformations and the stronglyanisotropic behavior of the myocardium, we employ the total Lagrangian SPHformulation. Then, the coupling of electrophysiology and tissue mechanics forelectromechanical feedback is conducted in unified SPH framework. A compre-hensive and rigorous study of iso- and aniso-tropic diffusion process, transmem-brane potential propagation in the free-pulse and spiral wave pattern, passiveand active responses of myocardium, electrophysiology and electromechanics ina generic biventricular heart model has been conducted. The results demon-strate the robustness, accuracy and feasibility of the proposed SPH frameworkfor cardiac electrophysiology and electromechanics.The SPH methods developed in this work are the main components of anunified meshless approach for multi-physics modeling of total cardiac function.42n the future work, an open-heart simulator based on our opensource SPHinXsyslibrary will be developed for total human heart modeling. In particular, fullycoupled fluid-electro-structure interactions, which consist of four chambers andfour valves as electrically excitable, deformable and electroactive bodies inter-acting with the blood flows will be modeled. These studies will be beneficial forunderstanding fundamental mechanisms of the total cardiac function.
7. Acknowledgement
The authors would like to thank Dr. Luca Ratti for sharing the dataset usedin Section 5.3 to validate our results and express their gratitude to DeutscheForschungsgemeinschaft for their sponsorship of this research under grant num-ber DFG HU1572/10-1 and DFG HU1527/12-1.
ReferencesReferences [1] W. H. Organization, The top 10 causes of death, , [On-line; accessed 7-Jan-2020] (2018).[2] A. Quarteroni, A. Manzoni, C. Vergara, The cardiovascular system: math-ematical modelling, numerical algorithms and clinical applications, ActaNumerica 26 (2017) 365–590.[3] N. A. Trayanova, Whole-heart modeling: applications to cardiac electro-physiology and electromechanics, Circulation research 108 (1) (2011) 113–128.[4] P. J. Hunter, A. J. Pullan, B. H. Smaill, Modeling total heart function,Annual review of biomedical engineering 5 (1) (2003) 147–177.[5] O. C. Zienkiewicz, R. L. Taylor, P. Nithiarasu, J. Zhu, The finite elementmethod, Vol. 3, McGraw-hill London, 1977.436] C. S. Peskin, The immersed boundary method, Acta numerica 11 (2002)479–517.[7] L. B. Lucy, A numerical approach to the testing of the fission hypothesis,The Astronomical Journal 82 (1977) 1013–1024.[8] R. A. Gingold, J. J. Monaghan, Smoothed particle hydrodynamics: theoryand application to non-spherical stars, Mon. Not. R. Astron. Soc. 181 (3)(1977) 375–389.[9] X. Y. Hu, N. A. Adams, A multi-phase SPH method for macroscopic andmesoscopic flows, J. Comput. Phys. 213 (2) (2006) 844–861.[10] C. Zhang, X. Hu, N. A. Adams, A weakly compressible SPH method basedon a low-dissipation riemann solver, J. Comput. Phys. 335 (2017) 605–620.[11] C. Zhang, X. Y. Hu, N. A. Adams, A generalized transport-velocity formu-lation for smoothed particle hydrodynamics, J. Comput. Phys. 337 (2017)216–232.[12] M. Rezavand, C. Zhang, X. Hu, A weakly compressible sph method forviolent multi-phase flows with high density ratio, Journal of ComputationalPhysics 402 (2020) 109092.[13] C. Zhang, M. Rezavand, X. Hu, Dual-criteria time stepping for weaklycompressible smoothed particle hydrodynamics, Journal of ComputationalPhysics 404 (2020) 109135.[14] C. Zhang, M. Rezavand, X. Hu, A multi-resolution sph method for fluid-structure interactions, arXiv preprint arXiv:1911.13255.[15] M. Liu, Z. Zhang, Smoothed particle hydrodynamics (sph) for modelingfluid-structure interactions, SCIENCE CHINA Physics, Mechanics & As-tronomy 62 (8) (2019) 984701. 4416] X. Bian, Z. Li, G. E. Karniadakis, Multi-resolution flow simulations bysmoothed particle hydrodynamics via domain decomposition, Journal ofComputational Physics 297 (2015) 132–155.[17] R. FitzHugh, Impulses and physiological states in theoretical models ofnerve membrane, Biophysical journal 1 (6) (1961) 445–466.[18] R. R. Aliev, A. V. Panfilov, A simple two-variable model of cardiac excita-tion, Chaos, Solitons & Fractals 7 (3) (1996) 293–301.[19] P. C. Franzone, L. F. Pavarino, S. Scacchi, Mathematical cardiac electro-physiology, Vol. 13, Springer, 2014.[20] G. A. Holzapfel, R. W. Ogden, Constitutive modelling of passive my-ocardium: a structurally based framework for material characterization,Philosophical Transactions of the Royal Society A: Mathematical, Physicaland Engineering Sciences 367 (1902) (2009) 3445–3475.[21] M. P. Nash, A. V. Panfilov, Electromechanical model of excitable tissue tostudy reentrant cardiac arrhythmias, Progress in biophysics and molecularbiology 85 (2-3) (2004) 501–522.[22] L. A. Taber, R. Perucchio, Modeling heart development, Journal of elas-ticity and the physical science of solids 61 (1-3) (2000) 165–197.[23] J. Wong, S. G¨oktepe, E. Kuhl, Computational modeling of electrochemicalcoupling: a novel finite element approach towards ionic models for cardiacelectrophysiology, Computer methods in applied mechanics and engineering200 (45-46) (2011) 3139–3158.[24] J.-H. Wang, S. Pan, X. Y. Hu, N. A. Adams, A split random time-steppingmethod for stiff and nonstiff detonation capturing, Combustion and Flame204 (2019) 397–413.[25] T. Tran-Duc, E. Bertevas, N. Phan-Thien, B. C. Khoo, Simulation ofanisotropic diffusion processes in fluids with smoothed particle hydrody-45amics, International Journal for Numerical Methods in Fluids 82 (11)(2016) 730–747.[26] A. Panfilov, Three-dimensional organization of electrical turbulence in theheart, Physical Review E 59 (6) (1999) R6251.[27] J. J. Monaghan, Smoothed particle hydrodynamics, Annual review of as-tronomy and astrophysics 30 (1) (1992) 543–574.[28] H. Wendland, Piecewise polynomial, positive definite and compactly sup-ported radial functions of minimal degree, Advances in computationalMathematics 4 (1) (1995) 389–396.[29] S. Biriukov, D. J. Price, Stable anisotropic heat conduction in smoothedparticle hydrodynamics, Monthly Notices of the Royal Astronomical Soci-ety 483 (4) (2018) 4901–4909.[30] M. Rezavand, D. Winkler, J. Sappl, L. Seiler, M. Meister, W. Rauch, Afully Lagrangian computational model for the integration of mixing and bio-chemical reactions in anaerobic digestion, Computers & Fluids 181 (2019)224–235.[31] K. Ten Tusscher, D. Noble, P.-J. Noble, A. V. Panfilov, A model for humanventricular tissue, American Journal of Physiology-Heart and CirculatoryPhysiology 286 (4) (2004) H1573–H1589.[32] R. Vignjevic, J. R. Reveles, J. Campbell, Sph in a total lagrangian formal-ism, CMC-Tech Science Press- 4 (3) (2006) 181.[33] Y. Zhu, P. J. Fox, Smoothed particle hydrodynamics model for diffusionthrough porous media, Transport in Porous Media 43 (3) (2001) 441–471.[34] J. Crank, et al., The mathematics of diffusion, Oxford university press,1979. 4635] L. Ratti, M. Verani, A posteriori error estimates for the mon-odomain model in cardiac electrophysiology, Calcolo 56. doi:10.1007/s10092-019-0327-2 .[36] Y. Wang, L. Cai, X. Luo, W. Ying, H. Gao, Simulation of action potentialpropagation based on the ghost structure method, Scientific reports 9 (1)(2019) 10927.[37] F. Liu, I. Turner, V. Anh, Q. Yang, K. Burrage, A numerical method forthe fractional fitzhugh–nagumo monodomain model, Anziam Journal 54(2012) 608–629.[38] M. Aguirre, A. J. Gil, J. Bonet, A. A. Carre˜no, A vertex centred finite vol-ume jameson–schmidt–turkel (jst) algorithm for a mixed conservation for-mulation in solid dynamics, Journal of Computational Physics 259 (2014)672–699.[39] R. W. Ogden, Non-linear elastic deformations, Courier Corporation, 1997.[40] E. Garcia-Blanco, R. Ortigosa, A. J. Gil, C. H. Lee, J. Bonet, A new compu-tational framework for electro-activation in cardiac mechanics, ComputerMethods in Applied Mechanics and Engineering 348 (2019) 796–845.[41] M. Sermesant, K. Rhode, G. I. Sanchez-Ortiz, O. Camara, R. Andriantsimi-avona, S. Hegde, D. Rueckert, P. Lambiase, C. Bucknall, E. Rosenthal,et al., Simulation of cardiac pathologies using an electromechanical biven-tricular model and xmr interventional imaging, Medical image analysis 9 (5)(2005) 467–480.[42] J. L. Lacome, Smoothed particle hydrodynamics method in ls-dyna, in: 3rdGerman LS-DYNA forum, Bamberg, Germany, 2004.[43] A. F. Johnson, M. Holzapfel, Modelling soft body impact on compositestructures, Composite Structures 61 (1-2) (2003) 103–113.4744] R. Hedayati, S. Ziaei-Rad, A new bird model and the effect of bird geometryin impacts from various orientations, Aerospace Science and Technology28 (1) (2013) 9–20.[45] L. Fu, L. Han, X. Y. Hu, N. A. Adams, An isotropic unstructured meshgeneration method based on a fluid relaxation analogy, Computer Methodsin Applied Mechanics and Engineering 350 (2019) 396–431.[46] A. Quarteroni, T. Lassila, S. Rossi, R. Ruiz-Baier, Integrated heart—coupling multiscale and multiphysics models for the simulation of the car-diac function, Computer Methods in Applied Mechanics and Engineering314 (2017) 345–407.
Appendices
Appendix A : Algorithms of the SPH framework for the simulation of cardiacfunctions. lgorithm 1: Algorithms of the newly proposed SPH framework for thesimulation of cardiac functions. Setup parameters and initialize the simulation; Compute the correction matrix B for each particle; while simulation is not finished do Compute time-step size ∆ t p ; Integrate the ODE ˙ V m = C m I ion for half time step ∆ t p ; Integrate the ODE ˙ w = g for half time step ∆ t p ; Integrat the diffusive operator ˙ V m = C m ∇ · ( D ∇ V m ) for a time step∆ t p ; Integrate the ODE ˙ w = g for half time step ∆ t ; Integrate the ODE ˙ V m = C m I ion for half time step ∆ t p ; if active response is considered then Compute time-step size ∆ t m and choose the time step∆ t = min (∆ t p , ∆ t m ); Integrate the ODE ˙ T a = f with the QSS method for a time step ∆ t ; Compute the active first Piola-Kirchhoff stress; Integrate the elastic equations; end16