Simulation of action potential propagation based on the ghost structure method

In this paper, a ghost structure (GS) method is proposed to simulate the monodomain model in irregular computational domains using finite difference without regenerating body-fitted grids. In order to verify the validity of the GS method, it is first used to solve the Fitzhugh-Nagumo monodomain model in rectangular and circular regions at different states (the stationary and moving states). Then, the GS method is used to simulate the propagation of the action potential (AP) in transverse and longitudinal sections of a healthy human heart, and with left bundle branch block (LBBB). Finally, we analyze the AP and calcium concentration under healthy and LBBB conditions. Our numerical results show that the GS method can accurately simulate AP propagation with different computational domains either stationary or moving, and we also find that LBBB will cause the left ventricle to contract later than the right ventricle, which in turn affects synchronized contraction of the two ventricles.

The heart is a rhythmic pump that maintains blood circulation throughout the body 1 . The rhythmic beating of the heart is caused by the regular spread of action potential (AP) within the heart. Abnormal conduction of AP in the heart can cause arrhythmias. Symptoms of arrhythmia include extrasystole, tachycardia, ventricular fibrillation, etc., of which ventricular fibrillation is the leading cause of cardiac sudden death 2,3 . Sudden cardiac death accounts for 15% of global deaths, and about 80% of sudden cardiac death is the result of ventricular arrhythmias 4 . Furthermore, about a quarter of patients with heart failure are diagnosed with LBBB 5 , which causes asynchronous AP propagation and contraction of the left ventricle, and then potentially leads to the global left ventricle dysfunction 6 . Therefore, it is of great significance to study the mechanism of arrhythmia, such as through numerical modelling, which can explore extreme situations that is difficult to perform in experiments. For several decades, the electrical activity of the heart has been modeled by a system of singularly perturbed reaction-diffusion partial differential equations that couples a set of ordinary differential equations used to describe the cell membrane dynamics 7,8 . The effects of different types of the electrical stimulation on arrhythmia can then be studied by solving these differential equations numerically. At present, numerical simulation of electrical activity has become a powerful tool for studying and understanding cardiac electrophysiology and arrhythmia [9][10][11] .
To mathematically model cardiac action potential, a cardiomyocyte model is required. With the abundance of experimental data, myocyte models have been continuously improved. A large number of mammalian cardiomyocyte models already exist in the literature 12 , such as Beeler-Reuter model 13 , Luo-Rudy model 14 , Fenton-Karma model 15 , etc. In order to accurately study the human heart, a large number of human cardiomyocyte models have also been proposed, such as ten Tusscher model 16 , Grandi-Pasqualini-Bers (GPB) model 17 , etc. For example, the GPB model can be used to describe Ca 2+ handling and ionic currents in human ventricular myocytes, and its effectiveness has been validated against available experimental data 7,17,18 . In 1969, Schmitt et al. 8 proposed a bidomain model for AP propagation in tissue level, then was further developed in the late 1970s [19][20][21] . The bidomain model describes active cardiomyocytes on a macroscopic scale by membrane ion current, membrane potential and extracellular potential 22 . Based on a given membrane potential, the bidomain model can modelling both the extracellular potential and the body-surface potential 23 . Recently, Bendahmane et al. 24 introduced a "stochastically forced" version of the bidomain model that accounted for various random effects, and further a two-dimensional FitzHugh-Nagumo model in a rectangular and circular regions 40,53 . The FitzHugh-Nagumo model is where, u is a normalized transmembrane potential, v is a recovery variable. K x and K x are the components of diffusion coefficient K. The model parameters a = 0.1, ε = 0.01, β = 0.5, γ = 1, σ = 0. In this Fitzhugh-Nagumo model, when K x = K y = 10 −4 , Fig. 1 gives the spiral wave of the stable rotation solution at t = 1000. As can be seen from Fig. 1, for the rectangular region, the spiral wave of the model generates a clockwise rotation curve. Figures 1(a,b) give the spiral waves obtained by the GS method and Liu 40,53 , respectively. It can be found that the spiral wave structure obtained by the GS method is consistent with the results obtained by Liu 40 . Figure 2(a,b) show the numerical results obtained by the GS method and Liu et al. 40,53 with K x = 10 −4 , K y /K x = 0.25, and Fig. 2(c,d) are the results with K y = 10 −4 , K x /K y = 0.25. It can be found that for K x = 10 −4 , K y /K x = 0.25 and K y = 10 −4 , K x /K y = 0.25, the spiral wave structures obtained by the GS method are nearly identical as those obtained by Liu et al. 40,53 .
Secondly, using the same model parameters and boundary condition, the computational domain is changed to be a circle Ω = {(x, y)|(x − 1.25) 2 + (y − 1.25) 2 ≤ 1.25 2 } with the following initial conditions: The computational domain of the ghost structure remains the same as the rectangular case. Figure 3(a,b) show the numerical results obtained by the GS method and Liu 41 with K x = K y = 10 −4 at t = 1000. It can be seen from Fig. 3 that the results obtained by the GS method in the circular region are also identical to those obtained by Liu 41 .
Finally, we simulate the transmembrane potential propagation in moving regions by using the same model parameters and boundary condition. To do this, we assume the Lagrangian point in the rectangular and circular region mentioned above expands in the normal direction → =  AP propagation on ventricular section. The electrical conduction system of the heart triggers myocardial contraction by electrical impulses transmitted through the sinus node. As shown in Fig. 6, electrical pulses pass through the atrium to the atrioventricular node and enter the ventricle along the left bundle branch, right bundle branch and Purkinje fiber. When an electrical pulse is transmitted to the cardiomyocytes and triggers the AP production, an excitation-contraction coupling occurs in cardiomyocytes. Myocardial contraction is highly dependent on the dynamics of calcium in a single myocardial cell 54 , which is becasue myofilament contraction is regulated by an increase of the intracellular calcium transient (CaT). Therefore, the excitation-contraction coupling of myocytes essentially depends on the calcium-induced calcium release 55 .
In this section, we employ the GPB model to model the myocyte electrophysiology, which describes various ionic current I ion and Ca 2+ dynamics. The reasons for choosing the GPB model are that (1) it matches experimental data well 17 ; (2) it is adequate to analyse AP with detailed Ca 2+ dynamics, which plays a crucial role in excitation-contraction in myocardium 54 . In order to understand the AP propagation in ventricles, we select one transverse and one longitudinal sections obtained from a real human heart 51,56,57 . The computational domain for the ghost structure of the transverse and longitudinal sections of the ventricle are 130 mm × 110 mm and 140 mm × 140 mm, respectively, and the spatial step size in each direction is 0.14 mm. The discrete time step is 0.08 ms, the membrane capacitance C m = 1 μF/cm 2 and the surface-to-volume ratio 58 Am = 0.24 μm −1 . The human ventricular conductivity is listed in Table 1. σ i L and σ e L are the longitudinal intra-and extracellular conductivities, respectively. σ i t and σ e t are the transversal intra-and extracellular conductivities, respectively. As shown in Table 1, the ventricular conductivity depends on the region and direction of cardiomyocyte. In this study, we employ the data obtained by Potse et al. 28 in Table 1.
The right and left bundle branches and Purkinje are the main conduction system in the ventricles. The right and left boundle branches divide into a few major branches and subsequently into Purkinje fibers 59 as shown in Fig. 6. Purkinje fibers penetrate into the ventricular muscle, entangled in endocardium, and form a network. The network of Purkinje fibers do not contribute to the activation of ventricular muscles until it reaches the middle and lower third of the septum and ventricle 45 . Purkinje fibers allow for rapid, coordinated, and synchronous physiologic depolarization of the ventricles. Therefore, we consider the initial electrical stimulation to be located in the middle and lower third segments of the septum and endocardium 60 , as indicated by the blue region in Fig. 6(a). For a healthy heart in Fig. 6(a), the blue region of endocardial surface will have the electrical stimulus transmitted from the Purkinje network. The LBBB refers to the blockage of the left bundle branch conduction. In this study, we consider complete left bundle branch block. Thus, no electrical stimulus is transmitted in the left branch, but through the interventricular septum from the right ventricular endocardium to the left ventricle endocardium 45,61 . As shown in the Fig. 6(b), only the blue region in the right ventricle receives the electrical stimulus from the Purkinje fibers.
To understand the effects of LBBB, we now compare AP propagations in a healthy heart with and without LBBB. Figure 7 shows AP propagation in the transverse section of heart at different times. Figure 7(a-d) illustrate the state of AP propagation under healthy condition. Figure 7(e-h) shows AP propagation across the transverse section of heart with LBBB. Figure 8 shows AP propagation on the longitudinal section of heart. In the transverse section of heart, the healthy heart completes the AP propagation within approximately 347.2 ms, which is faster www.nature.com/scientificreports www.nature.com/scientificreports/ than the propagation with LBBB (540 ms). As shown in Fig. 7, under healthy conditions, when t = 40 ms, the AP spreads to most of the right ventricle and half of the left ventricular region. However, in the LBBB case, only half of the right ventricular region is excited, and the AP has not yet arrived at the left ventricular. When t = 120 ms, the whole healthy heart is repolarized, but still half of the left ventricle is not excited in the LBBB case. In the longitudinal section, the heart with LBBB needs 1011.2 ms to complete the AP propagation, which is 2.4 times the time (416 ms) that the healthy heart completes the propagation. As shown in Fig. 8, at t = 120 ms, almost the whole healthy heart is stimulated, but for the heart with LBBB, only the right ventricle is stimulated. When t = 360 ms, the AP only spreads to the apex of the heart with LBBB, while in the healthy heart, almost all regions return to the resting potential. Therefore, the LBBB causes significant delay in the activation of the left ventricle.
To further explain how LBBB affect ventricular contraction, we select four points on the transverse and longitudinal sections, as shown in Fig. 6. Transmembrane potential and calcium ion concentration at each point are then analyzed. As shown in Fig. 9, the points in the healthy case are stimulated almost simultaneously, but not in www.nature.com/scientificreports www.nature.com/scientificreports/ the heart with LBBB, much delayed at points 4, 7 and 8. Comparing Fig. 9(a) with Fig. 9(b), the activation time of the point in the interventricular septum (such as point 2) in the LBBB case is similar to that in the healthy case. However, the activation time is significantly affected by LBBB at points away from the right ventricle. The farther away from the right ventricle, the later the activation. Figure 10 shows the calcium ion concentration at selected points. Since myofilament contraction is regulated by intracellular calcium transient, it can be seen from Fig. 10(a,c) that all points in the healthy cases can contract at the same time. While in the heart with LBBB ( Fig. 10(a-d)), the points in the interventricular septum (e.g. points 2, 3 and 6) are essentially unaffected and will contract at nearly same time, but points 4 and 7 are about 250 ms late when they start to contract. For the point farthest from the right ventricle, namely point 8 at the lateral wall, the contraction time is about 600 ms later. Therefore, the LBBB will cause significant contraction delay in the left ventricular lateral wall, could potentially lead to heart failure in the long term. www.nature.com/scientificreports www.nature.com/scientificreports/

Discussion
The validity of the GS method is verified by some standard examples of the FitzHugh-Nagumo monodomain models. Based on the GS method, we capture the patterns of heterogeneity and complex connectivity of electrophysiological dynamics in biological tissues by solving the Fitzhugh-Nagumo monodomain model in rectangular and circular regions. The spiral wave obtained by the GS method is the same as that obtained by others, and it is verified that the GS method can effectively solve the monodomain model in the rectangular and circular regions. By comparing with the results in the stationary region, propagation velocity and shape of spiral waves in the moving region change. The propagation velocity in the moving region is higher than the velocity in the stationary region. The width of the spiral wave in the moving region is also slightly larger than the width in the stationary region. Since the GS method permits nonconforming discretization of the transmembrane potential and membrane dynamics, the monodomain model can be directly solved using finite different method in a regular ghost structure. Another advantage of the GS method is in dealing with the moving regions. Compared to Liu's method 40 , the GS method needs to use the delta function to transform the Lagrangian and Eulerian variables, and the membrane dynamics need to be solved in a finer Lagrangian grid. In this sense, the GS method will require higher computational resource than Liu's implementation 40 . For the stationary rectangular region, the computational time of the GS method is about 3 hours, since we have not fully optimized the implementation for high performance computing but only employed OpenMP functionality for dealing with "for". It is expected that once we employ MPI and GPU computing, the computational time can be reduced significantly.
The AP propagation on the transverse and longitudinal sections of human heart is successfully simulated by the GS method. In a real heart, the AP transmits involves varieties of conduction cells, such as myocyte, sinoatrial node cells, atrioventricular nodal cells, Purkinje fibers, and fibroblasts. The electrical activity of the heart begins   www.nature.com/scientificreports www.nature.com/scientificreports/ with the sinoatrial node at the right atrium. Pulses from the sinoatrial node travel through the left and right atrium and meet at the atrioventricular node. From the atrioventricular node, electrical impulses travel along the bundle and are transmitted to the right and left ventricles through the right and left bundle branches. Finally, the bundle at the end of the bundle branch is divided into millions of Purkinje fibers. Nowadays, many researchers have begun to study electrophysiology by including various conduction cells. The sinoatrial node is the normal pacemaker of the mammalian heart. There are a few mathematical models of sinoatrial nodes. For example, based on the Severi-DiFrancesco model of a rabbit sinoatrial node cell and the electrophysiological data from human sinoatrial node cells, Fabbri et al. 62 proposed a comprehensive model of the electrical activity of a human sinoatrial node cell. The AP and CaT obtained by Fabbri were close to experimentally recorded values. In order to illustrate the functional role of various genetic isoforms of ion channels in generating cardiac pace-making AP, Kharche et al. 63 developed a mathematical model for spontaneous AP of mouse sinoatrial node cells. In that model, biophysical properties of membrane ionic currents and intracellular mechanisms were considered. Results showed that their model could reproduce the physiological exceptionally short AP and high pacing rates of mouse sinoatrial node cells effectively.
Because of the importance of Purkinje system in both normal ventricular excitation and ventricular arrhythmias, modelling of the Purkinje system is essential for a realistic ventricle model of the heart 61 . Recently, inclusion of Purkinje network in AP modelling has attracted much attention 64 . For example, Oleg et al. 65 developed a detailed model of the canine Purkinje-ventricular junction and varied its heterogeneity parameters to determine www.nature.com/scientificreports www.nature.com/scientificreports/ the relationship between wave conduction velocity, tissue structure, and safety of discontinuous conduction at nonuniform junctions. Oleg found that fast or slow conduction was unsafe, and there existed an optimal velocity that provided the maximum safety factor for conduction through the junction. Vergara et al. 66 developed a model for the electrophysiology in the heart to handle the electrical propagation in the Purkinje system and myocardium. Their results illustrated the importance of using physiologically realistic Purkinje-trees for simulating cardiac activation. However, the majority of current anatomical models have not included models of the Purkinje network 61 . Instead, a simplified approach is adopted by applying electrical stimulus in the middle and lower third segments of the septum and endocardium 60 . The same approach is used in this study. This is partially due to the fact that extensive branching of the Purkinje fibers makes modelling Purkinje network extremely challenging as suggested by Tawara 59 , who carried out a formidable study lasting over 2 years to reconstruct the conduction system from experimental data. Since the focus of this study is to develop a numerical method for AP within myocardium, like many other studies 67,68 , we only consider myocytes, and other conduction cells are not simulated.
When the bundle branch is injured after myocardial infarction, or cardiac surgery, it may stop transmitting electrical impulses completely. This will result in a change in the path of ventricular depolarization. According to the anatomical location of the defect that leads to the bundle branch block, the block is further divided into the right bundle branch block and the left bundle branch block. Since the electrical pulse can no longer use the preferred path through the bundle branch, it can only spread through muscle fibers, which slows down the electrical propagation and changes the directional propagation of the electrical pulse. Lange et al. 44   www.nature.com/scientificreports www.nature.com/scientificreports/ tendon could compensate for the propagation delay caused by the LBBB. As demonstrated in this study, LBBB leads to delayed triggering of electrical excitation of the left ventricle, resulting in the loss of ventricular electrical synchrony, and potentially causing mechanical discoordination.

Conclusion
In this study, we have developed a GS method by immerse the actual irregular electrophysiology computational domain into a larger rectangular region. Action potential propagation using the monodomain model is solved successfully with the GS method. In a rectangular and a circular regions, by using the GS method to solve the FitzHugh-Nagumo monodomain model, we capture the patterns of heterogeneity and complex connectivity of electrophysiological dynamics in biological tissues, and demonstrate the validity of the method. Numerical results show that the GS method can effectively simulate the AP propagation in irregular region. Furthermore, we employ the GS method to simulate the transmembrane propagation in moving regions and analyze the influence of moving region on transmembrane propagation. Our results show that the moving regions affect not only the propagation velocity but also the shape of spiral waves. Subsequently, we simulate the AP propagation on the transverse and longitudinal sections of a healthy heart and a heart with LBBB by using the GS method. The numerical results demonstrate how LBBB affects action prorogation in ventricles.

Model Introduction
Monodomain model. Simulating myocardial electrical activity needs to describe the anisotropic excitation conduction based on the ion channel model of myocardial cells. In general, a bidomain model based on the Poisson equation is used to describe the electrical coupling between myocytes and the electrical conduction cells in the tissue. At the microscopic level, myocardium can be seen as consisting of two separate regions separated by the cell membrane: the intracellular space (Ω i ) and the extracellular space (Ω e ). The bidomain model consists of the equations for the intracellular potentials (φ i ) and the extracellular potentials (φ e ), thus the transmembrane potential is V m = φ i − φ e . The governing equations of the bidomain model are www.nature.com/scientificreports www.nature.com/scientificreports/ e l e t are the conductivity tensor. A m is the surface-to-volume ratio, i.e., the amount of membrane found in a given volume of tissue. I m is the transmembrane current density. C m is the membrane capacitance, I ion is the ionic current through a number of different types of ion channels. I s is an imposed stimulation current.
Assuming the anisotropy ratios in intracellular and extracellular spaces are the same, let σ e = λσ i , then Eq. (8) will be reduced to Substituting Eq. (10) into Eq. (9), we will obtain the governing equation of the monodomain model, that is systems. They are singularly perturbed systems for model parameters and reasonable initial data 69 . When the diffusion phenomenon is included in the system, the threshold phenomenon ensures the stability of the traveling wave solution. The so-called threshold phenomenon, that is, there is a threshold value of the transmembrane potential V m in the uniform space, for the electrophysiology model Eq. (9) or Eq. (11), when the potential is lower than the threshold value, it quickly returns to the stable state; when the potential is higher than the threshold value, it will produce a large excursion before it returns to the stable state. In cardiac tissue, an initial perturbation with a sufficiently large transmembrane potential V m triggers AP propagation. www.nature.com/scientificreports www.nature.com/scientificreports/ Figure 11 shows a typical AP curve for a human myocyte, which consists of four phases, the depolarization phase (phase 0), the early repolarization phase (phase 1), the plateau phase (phase 2), the repolarization (phase 3), and the resting phase (phase 4). In a cardiac cycle, once a myocyte is excited, it can not be excited again for a period of time, the so-called effective refractory period (ERP). During this period, the depolarization of adjacent cardiomyocytes does not trigger already-excited myocytes. When entering the resting phase (phase 4), myocytes are ready for next excitation. ERP is usually characterized by the interval between the depolarization (phase 0) and repolarization (phase 3) phases. As a protective mechanism, ERP can control the heart rate, prevent arrhythmias and coordinate muscle contractions.
The ghost structure method. In this ghost structure method, the transmembrane potential V m is described by an Eulerian form and discretized with a regular Cartesian grid, while the membrane dynamics is described by a Lagrangian form and calculated by the cell membrane model. As shown in Fig. 12, the entire computational domain (the ghost structure region) is represented by Ω, where X = (X 1 , X 2 ) ∈ Ω c is Eulerian (physical) coordinates. The region where myocytes sit is expressed as Ω c , X = (X 1 , X 2 ) ∈ Ω c are fixed Lagrangian (material) coordinates. The mapping χ(X, t) ∈ Ω gives the physical position of each Lagrangian point at time t. Therefore, the physical region occupied by myocardium at time t is denoted as Ω c (t) = χ(X, t), and the region of non-myocardial cells at time t is denoted as Ω non (t) = Ω − Ω c . The Lagrangian and Eulerian variables are transformed by the integral transformation of the delta function. Finally, the full governing equations of the monodomain model in the GS method are  www.nature.com/scientificreports www.nature.com/scientificreports/ s s c where V m (x, t) is the Eulerian transmembrane potential, I ion (x, t) is the Eulerian ionic current, and I s (x, t) is the Eulerian stimulation current density.
s are the transmembrane potential, ionic current and stimulation current density in Lagrangian form, respectively. y is a Lagrangian vector of ionic fluxes and their corresponding channel gating variables are described by the suitable ordinary differential equations Eq. (14). f represents the right hand side of the ordinary differential equations used to describe ion channels. g represents a nonlinear function that relates the ionic flux to the total ionic current. Eq. (14) and Eq. (15) are cardiac membrane models used to solve the current-voltage relationship. δ(x) is a two-dimensional delta function that transforms the transmembrane potential and ion current between the Eulerian and Lagrangian coordinates. In this study, we calculate ion currents by Eq. (13) and Eq. (14) through the GPB model, and then convert  I t X ( , ) ion into the Eulerian ion current in the ghost structure region. Finally, the Eulerian transmembrane potential V m (x, t) is obtained by solving Eq. (12). Spatial discretization. The regular ghost structure region Ω is discrete by a N 1 × N 2 Cartesian grid with spatial steps ∆ = . The gradient of V m is approximated at the center point of the grid again using the following difference scheme, www.nature.com/scientificreports www.nature.com/scientificreports/ The divergence of σ i ∇V m is defined also at the centre of the grid as The initial value of V m at each Eulerian point in the entire computational domain needs to be approximated based on the transmembrane potential at the Lagrangian point and the integral transformation of the delta function. In order to solve more accurately, the transmembrane potential at the Lagrangian point in Ω c and the delta function are required to update the transmembrane potential at the Eulerian point outside the structure Ω c at regular intervals. The mutual transformation between the Eulerian variable and the Lagrangian variable will be explained in the subsection "Lagrangian-Eulerian interaction".
Time discretization. For the ordinary differential equations Eq. (14), the third-order TVD Runge-Kutta method 70 is used to solve the system: n m n n m n m n 1 ( 2) 1 1 (1) 1 (2) For the nonlinear partial differential equation Eq. (12), it can be rewritten as . In this paper, the third-order TVD Runge-Kutta method is also used to solve the nonlinear reaction-diffusion equation Eq. (21) In this paper, nodes of the mesh are regarded as Lagrangian points. The Lagrangian points must be finer than the Cartesian points to avoid leaks 71 . In the calculation process, the integral transformation form of delta function is used to realize the conversion between the Eulerian variable and www.nature.com/scientificreports www.nature.com/scientificreports/ Based on the above delta function, the approximate values of physical quantities at each Eulerian point can be obtained directly by using the values of Lagrangian points around the Eulerian point. In order to obtain the approximate value of physical quantity in Eulerian coordinate system more accurately, we use the Gaussian quadrature rule with quadrature points X Q e , where ∈ Ω X Q e e , and accociated weights ω = ... Q  N  (  1, , ) Q e e . In the GS method, the value of each Gaussian integral point is obtained through the basis function of the finite element mesh, and then the approximate value of the physical quantity at the Eulerian point is obtained by using the integral transformation form of the delta function. Taking I ion (x, t) as an example, the current  I ion of the Gaussian integration point in the element is obtained by the value of the Lagrangian point through the element basis function of the finite element method, and the approximate value I ion of the  I ion at the Eulerian point is Similarly, we can use the integral transformation form of the corresponding delta function to obtain the approximate value of V m at each Gaussian integral point by using the transmembrane potential V m at Eulerian points, that is In order to update the ion current at the next time step through the cell membrane model, we need to obtain the transmembrane potential ∼ V m at the Lagrangian point or the grid nodes. The transmembrane potential ∼ V m at grid nodes is obtained by using the approximate value of Gaussian integral point from Eq. 27 and a L 2 projection method 48 .
All simulations are performed on a windows workstation with Intel(R) Xeon(R) Gold 5115 (20 cores, 2.40 GHz, 64 GB memory), implemented in C++.