Abstract
This study aims to investigate the computationally complex cardiac dynamics of the well-known human ventricular model of Ten Tusscher and Panfilov from 2006. The corresponding physical system of the cellular model is described by a set of nonlinear differential equations containing various system parameters. If one or a set of specific system parameters crosses a certain threshold, the system is forced to change dynamics, which might result in dangerous cardiac dynamics and can be precursors to cardiac death. To perform an efficient numerical analysis, the original model is revised and simplified such that the modified models perfectly match the trajectory (time-dependent cardiac potential) of the original model. Moreover, it is demonstrated that the reduced models have the same dynamics. Furthermore, using the lowest dimensional model, it is shown by means of bifurcation analysis that combinations of reduced slow and rapid potassium currents and an enhanced calcium current may lead to early afterdepolarizations, which are pathological voltage oscillations during the repolarization or plateau phase of cardiac action potentials and are considered potential precursors to cardiac arrhythmia. Finally, to outline synchronization effects and pattern formation on a larger scale (macro scale), a two-dimensional epicardial monodomain equation is studied.
1 Introduction
Computational physiology and medicine, computational physics, biology and biophysics, mathematics for healthcare, and modeling of biomedical applications have gained importance in numerous interdisciplinary and multidisciplinary research projects [1–4]. Mathematical modeling has become an integral part and contributes significantly to a better understanding of real-world phenomena, such as cardiac or neuronal dynamics [5–14]. In addition to mathematical modeling and simulations, the analysis of these complex systems has increasingly become the focus of research. In the field of mathematical and computational cardiology, a strong focus is on the investigation of certain cardiac arrhythmias, such as early afterdepolarizations (EADs). EADs are pathological voltage oscillations during the repolarization or plateau phase of cardiac action potentials (APs) (cf. Figure 1) and are considered potential precursors to cardiac arrhythmia, often associated with potassium deficiency or elevated calcium or sodium currents caused, for example, by ion channel diseases, drugs, or oxidative stress.
FIGURE 1
The presence of EADs strongly correlates with the onset of dangerous cardiac arrhythmia, including torsades de pointes, polymorphic ventricular tachycardia, and fibrillation, which are specific types of abnormal heart rhythms that can lead to sudden cardiac death [16, 17]. This potential spreading of EADs from cells to the whole heart, resulting in dangerous cardiac arrhythmias, provided the number of affected cells is large enough, is called EAD synchronization [16]. It is therefore highly important to fully understand the mechanism of such a phenomenon on the cellular as well as the tissue/organ level.
The main emphasis in this research field has been the modeling of cardiac muscle cells and the heart and the computational studies of these models [18, 19]. However, nowadays the focus moves towards systematic and more mathematical investigations of phenomena on the cellular level, such as EADs [20–23]. Bifurcation theory has proven to be a good tool for analyzing these systems by investigating changes in the qualitative or topological structure of a particular family of dynamic systems. It seeks to explain how the behavior of a system (e.g., its stable states or periodic cycles) changes with variation in a system parameter. In the context of cardiac dynamics, this often means studying how cardiac muscle cell dynamics such as action potentials or the heart rhythm can change—from normal sinus rhythm to arrhythmias such as tachycardia or fibrillation—when physiological parameters (such as ionic conductance, and pacing rate) are altered. Bifurcation theory is useful because it helps identify thresholds or critical points where the dynamics of the cardiac cell or the heart shift from one behavior to another, such as a small change in a physiological parameter leading to sudden cardiac arrhythmia. Bifurcation theory is particularly suited to systems with nonlinear dynamics, where other approaches fail; thus, bifurcation analysis is utilized in this manuscript.
This approach has mainly been used for studies of low- (up to four) dimensional ODE models [21, 24, 25]. Potential reasons for this are that bifurcation analysis of high-dimensional nonlinear systems is more challenging and might fail more easily due to imprecise settings. Further reasons are issues occurring due to the modeling of cell models themselves since they might not be smooth and regular enough. However, one main issue considering such a reduced problem is the loss of information [26, 27]. On the other hand, one can argue (for example, in [24, 25, 28]) that some sort of time-scale separation argument [29] is used to identify the gating variable of the current as the slowest variable to reduce the dimension of the system to , by using this gating variable as a bifurcation parameter. The advantage of this approach is that the analysis of the rapid subsystem is easier, but again, one may lose information cf. [27]. Moreover, this approach is not applicable to studying the effect of the current or the current since their gating variables are faster than that of the current.
The aim of the paper is to provide a mathematically feasible and physiologically realistic model for analyzing the electrophysiological dynamics of a cardiac muscle cell. Based on the reduced model, the bifurcation analysis enables the systematic investigation of the ionic current interactions, particularly the interaction of the calcium current with the slow and fast potassium currents. This analysis provides a deeper insight into how an enhanced calcium current in combination with reduced slow and fast potassium currents leads to EADs. In addition, the result of the bifurcation analysis allows quantification of these dangerous interactions. This study is divided into three parts. 1) To decode the dynamical structure that leads to EADs, we perform model reductions for the highly nonlinear and high-dimensional cell model. The reduced models are compared, and their advantages are emphasized. Furthermore, we show that the simplified models reproduce the dynamics of the full model. The aim is to obtain a model that is as simple and numerically efficient as possible to simulate and analyze it without losing information about the behavior of the voltage of the original model. These reduced models retain the essential nonlinear features of the full system but are better suited for mathematical analysis, such as bifurcation theory. 2) We investigate the reduced model by means of numerical bifurcation analysis and determine the complex behavior of the cardiac cell model. Furthermore, we will show that the occurrence of EADs in both the full and reduced model is determined by the same underlying bifurcation mechanism—a period doubling cascade. This shows that EADs arise on a generic route to chaos, which can also be captured and explained in simplified models. In the case of a stable period doubling cascade, chaos can also be simulated. Using bifurcation analysis, we identify critical transitions in the behavior of the system when key parameters are varied. In particular, we focus on the occurrence of EADs through a period doubling cascade. This approach allows (i) confirmation of the validity of the reduction; (ii) accurate determination of the conditions for the occurrence of EADs; (iii) predictions that can serve as a basis for further numerical or experimental studies. 3) Finally, based on these results, we will perform numerical simulations for the two-dimensional monodomain model to determine macroscale synchronization effects.
The motivation for this is as follows. Cardiac myocytes can exhibit complex oscillatory patterns, such as spiking and bursting, that are related to ion–current interactions. In addition to normal APs of a cardiac myocyte, certain types of cardiac arrhythmias can occur. These include certain types of cardiac arrhythmias that can lead to sudden cardiac death. Furthermore, irregular behavior such as (deterministic) chaos or chaotic EADs has been observed in both experimental and computational studies (cf. [30–33] and the references contained therein). Moreover, heart dynamics or the heart rhythm can react very sensitively to the influence of certain medications, and computational studies can provide new insights and help make better predictions cf. [34–39]. To this end, we modified and simplified the cardiac muscle cell model of Ten Tusscher and Panfilov from 2006 [15, 40] to perform numerical bifurcation analysis. This allows the investigation of the dynamics of the TP06 model and prediction of the occurrence of normal APs or cardiac arrhythmias, such as EADs. The focus of this study is on the model reduction and enhancement of the efficiency of the numerical bifurcation analysis without loss of information to the original TP06 model.
2 Methods
This section focuses on the methodology used in this manuscript, including the modeling of cardiac single models and the heart, as well as model reduction. Moreover, the bifurcation analysis and the numerical methods are discussed.
2.1 Cardiac single-cell modeling
Here, we briefly describe the cardiac model we will investigate in this paper. The mathematical modeling of action potentials (APs) of excitable biological cells, such as cardiac myocytes, has its origin in the Hodgkin–Huxley model [
41]. Here, an approach was developed to model the APs of excitable biological cells through a system of ordinary differential equations (ODEs). These conductance-based models represent a minimal biophysical interpretation of excitable biological cells in which current flow across the membrane is due to the charging of membrane capacitance and movement of ions through ion channels that are selective for certain ion species. An initial stimulus activates the ion channels once a certain threshold potential is reached. Then, these ion channels open and/or close, allowing an ion current to flow that changes the membrane potential. This electrophysiological behavior of a cardiac myocyte is represented by the following ODE:
where
denotes the voltage (in
) and
the time (in
), while
is the sum of all transmembrane ionic currents. The (epi-, mid-myo,- and endocardial) human ventricular TP06 model contains several different ion currents, ion pump, ion exchanger, and background currents:
These currents depend on individual ionic conductances
and Nernst potentials
. Moreover, they may depend on gating variables which are important for the activation and inactivation of ion currents. The gating variables satisfy the differential equation
where
represents the gating variables, while
denotes the equilibrium of the gating variable
and
its time scale. Furthermore, the ionic concentrations of the TP06 model from [
15,
40] read as:
Note that we do not specify the constants, but they are available in [
15,
40] or in the code provided in [
42]. Furthermore, the main ion currents except the sodium current are listed in the following.
L-type calcium current with gating variables , , and :
transient outward current with gating variables and :
slow delayed rectifier current with gating variable :
rapid delayed rectifier current with gating variables and :
inward rectifier current:
where the gating variables satisfy the differential
Equation 1. Finally, note the difference between the epi-, mid-myo-, and endocardial human ventricular TP06 model. There are two parts of the model which differ for these three cell types.
1. The modeling of the gating variable of the transient outward current is different:
and
2. The values of and differ:
and
For a full description of the variables involved, please see [
15,
40].
This also implies that we only need one proper bifurcation analysis with respect to for epicardial and M cells since the models are different only in . Furthermore, we need a separate analysis for endocardial cells. This also explains why the trajectories in Figures 2, 3 are different. For the simulations in Figures 2, 3, we used the same setting as in [15], except for the value , which we set to of the original value and the value , which we have chosen as five times bigger than that in [15]. Therefore, the M cell enters a region where EADs occur.
FIGURE 2
FIGURE 3
2.2 Monodomain model
To model the heart, different types of reaction–diffusion type systems are available, where the focus is on different research questions, and their complexities vary (cf. e.g., [43–48]). This paper focuses on the two-dimensional monodomain model since it is used to model the original TP06 model (cf. [15, 40]). In addition to this, monodomain models have been shown to be a good approximation for wave propagation in cardiac tissue (e.g., [49]) and are well studied (e.g., [50, 51]). The monodomain model reads as follows:with equal anisotropy rates , where is a scalar constant, is the spatial gradient, and denotes the intra- and extracellular conductivity tensors, and is the membrane surface area per unit volume. Moreover, the system is equipped with Neumann boundary conditions, , where denotes the normal vector. In the fashion of [39], we write the monodomain Equation 2 as a reaction–diffusion system, which it is:where denotes the diffusion coefficient . Moreover, we use the same mesh, time, and space discretization and stimulus, as in [39] (Table 1).
TABLE 1
| Spatial grid size | Time step | Stimulation duration | Stimulus strength | Second stimulus | Diffusion coefficient | Integration domain size |
|---|---|---|---|---|---|---|
Mesh, stimulation, and diffusion parameters.
Notice that for normal AP simulations, a time step of is used, while for the simulations of EAD settings, it is .
2.3 Model reduction
The aim of model reduction is to simplify a complex mathematical model while retaining its essential behaviors and properties so that the reduced model can be calculated more quickly and analyzed more easily. Therefore, we modify and simplify the TP06 cell model and then compare the dynamics of the resulting three modified TP06 models—one 18-dimensional, one 16-dimensional, and one 14-dimensional system of ordinary differential equations (ODEs). The 18-dimensional model is in the fashion of [52, 53] and perfectly represents the original model, while the 16-dimensional model gives a perfect approximation of the original model in case we modify the initial stimulus. The same applies to the 14-dimensional model. Notice that without the modification of the initial stimulus, the models show very similar dynamics compared to the original one; however, the trajectories do not perfectly coincide with that from the TP06 model.
2.4 Bifurcation theory
Very powerful tools to systematically analyze the dynamics of cardiac myocytes are bifurcation theory [54] and geometric singular perturbation theory [55] (see for instance, recent studies in [21, 26, 56–65] and the references contained therein). These computational studies can be used to develop new therapies and help in improving clinical decisions [66–70].
In general, the state of a physical system can be observed when it is stable, and one expects that a small change in a system parameter should not change its dynamics. Rather, stable solutions can be expected to continuously change in unique ways. No dramatic change is observed when varying any parameter, as long as a continuous solution branch retains its stability. However, when a certain physical parameter exceeds a threshold, the physical system may be forced to change its dynamics, and complex behavior may result. Therefore, bifurcation theory is used to explain certain phenomena and dynamics of the TP06 model of Ten Tusscher and Panfilov from 2006 (cf. [15]).
One requirement to be able to apply bifurcation theory is that the investigated system is sufficiently smooth. However, this is not the case for the original TP06 model (reported, for instance, in [53, 71]). In [53], the authors modified the sodium current similar to the approach used in [52], while in [71] the sodium current was modified based on [72]. Both approaches enable (numerical) bifurcation analysis with slight differences. The advantage of the approach used in [53] is that the remodeled system almost always perfectly fits the original TP06 model without changing the stimulus, while the ansatz in [71] leads to an almost perfect remodeling of the original TP06 model. On the other hand, the ansatz in [71] results in a lower-dimensional model, and therefore the bifurcation analysis is more efficient. These advantages and disadvantages must be considered.
The important bifurcations in this study are 1) Andronov–Hopf bifurcations; 2) torus bifurcations; 3) period-doubling bifurcations [54]. An Andronov–Hopf bifurcation occurs when a fixed point of a dynamic system changes its stability and a limit cycle (periodic orbit) occurs. In a supercritical Andronov–Hopf bifurcation, a stable limit cycle occurs. In the subcritical case, an unstable limit cycle appears. The Lyapunov coefficient determines the nature of the Andronov–Hopf bifurcation. Its sign indicates whether the bifurcation is supercritical (negative coefficient, stable limit cycle) or subcritical (positive coefficient, unstable limit cycle). It quantifies the nonlinear effects near the bifurcation point. A torus bifurcation (or Neimark–Sacker bifurcation) occurs when a limit cycle loses stability and a quasiperiodic orbit emerges on a torus. The system’s dynamics then evolve with two incommensurate frequencies, leading to motion on a 2D torus instead of a simple periodic loop. At period-doubling bifurcations, the periodic behavior of the system changes, giving rise to a new orbit with twice the period. Moreover, a period-doubling cascade is clearly linked to the occurrence of EADs and, in the case of a supercritical period-doubling cascade, is clearly linked to chaos (e.g., [30]).
2.5 Numerical methods
For our simulations, we utilized MATLAB R2023b and Python 3.9 with FEniCS [73, 74]. To solve the ODE system, we mainly use the MATLAB ode solver ode15s with a relative tolerance of and an absolute tolerance of . For the monodomain model, the two-dimensional simulations are done in FEniCS, where the coupled ODE–PDE system is solved using cbcbeat [75] or fenics-beat [76], which are Python/FEniCS-based libraries for running cardiac electrophysiology simulations. Note that cbcbeat and fenics-beat use a second-order splitting scheme with the Crank–Nicolson method for time stepping. Moreover, the Rush and Larsen scheme [77] is used to integrate the gating variables in time. For the bifurcation analysis, we used CL_MATCONT, a continuation toolbox for MATLAB [78, 79]. Codes are available via [42].
3 Results
The modifications we present here reduce the complexity of the 19-dimensional TP06 model in several steps to an 18-, 16-, and 14-dimensional system with astonishing consistency with the original model (Figures 2, 3). Here, we changed the conductance of the slow delayed rectifier current , rapid delayed rectifier current , and the L-type calcium current such that .
As well as the modifications and model reduction which we will present in the next section, we have to adjust the initial stimulus for the 16-dimensional model in Figure 2, using instead of . Moreover, the effect of the different stimuli can be clearly observed in Figure 3. In the first row of Figure 3, an initial stimulus of is applied to the 16-dimensional model. Although the trajectories do not match perfectly, the dynamics are similar. In the second row of Figure 3, perfectly matching trajectories of the TP06 model and its different versions are shown, highlighting whether a different stimulus is applied. Notice that the MATLAB code to produce Figures 2, 3 is provided in [42]. Using this code, the readers can also convince themselves that the modified systems effectively approximate the original system.
As well as the model reduction, the focus is on the complex dynamics these models can develop. Numerical bifurcation analysis is applied, highlighting that a reduction in the slow and rapid delayed rectifier current, and , in combination with an enhanced L-type calcium current , force the system to develop EADs. This can be linked to the existence of a subcritical period-doubling cascade. Moreover, it is shown that no chaos may appear.
Finally, we investigate the 2D synchronization effects. In the 2D case, we highlight that the monodomain equation’s dependence on the system parameters develops stable spiral waves, but wave break-up may also appear (cf. Figure 10) due to the EAD developed by the cellular model.
3.1 Model modification and reduction
First, we remove the ODE representing the intracellular ion concentration and set constant equal to the initial concentration of the TP06 model, . In [80], it was postulated that models for cardiac cells that account for changes in intracellular ion concentrations violate a conservation principle. Moreover, it turned out that the variation in the concentration is rather small (cf. [71]). Consequently, these systems never reach a steady state—a resting potential—which is required to be able to apply bifurcation theory. All modifications we present below are tabularly collected and compared in the first section of the supplementary material.
3.1.1 18 dimensional model
As indicated in [71], the first issue is in the modeling of the sodium current in (40, 15)where denotes the ionic conductance and the Nernst potential, while the different gating variables , , and satisfy differential Equation 1. To be more precise, the issue appears in the modeling of the gating variables and since the voltage-dependent functions , , , and are not continuous:
(cf. [40]). In the fashion of [52, 53], we introduce a new function
(cf. Supplementary Figure S8) and remodel the rate constants and as follows:Now, the modified model is sufficiently smooth without any discontinuity. Therefore, bifurcation analysis is applicable, and no further modifications are needed to fit the original model. Note that the function is modeled in such a way that for to represent the switching modeled in [40] of and at , while the factor “” might be improvable (cf. Supplementary Figure S8–10).
3.1.2 16 dimensional model
A further way to avoid the issue with the discontinuity of the rate constants and is to notice that the equilibrium of and is equal. In [71], the gating variables and are reformulated to one new gating variable , the sodium current is modified to
and a new time scale is introduced such that satisfies (Equation 2) with
The time relaxation constant is given bywhich could probably be improved. However, we realized that instead of considering a 17-dimensional model, it is more practical to set and to adjust the initial stimulus. This removes the discontinuity, and in addition, we reduce the model by a further dimension.
3.1.3 14 dimensional model
Finally, in this fashion and motivated by [81], we also fixed two further gating variables equal to their steady-state solutions—that is, we fixed the gating variable of the transient outward current and the gating variable of the rapid delayed rectifier current such that
This again requires modification of the initial stimulus, which must be higher than the 16-dimensional model. However, this simplified model shows a remarkably accurate approximation of the original model, and additionally, it is much less challenging to analyze, which we will highlight later in more detail.
3.2 Bifurcation analysis of the cardiac single-cell model
The aim of the study is to make statements about the behavior of the trajectories and the dynamics of the TP06 system. To this end, it investigates the stability and bifurcations of the modified systems (cf. [30, 53, 54, 82]). Moreover, we provide the corresponding codes in [42].
A bifurcation of a dynamic system is a qualitative change in its dynamics caused by the change of parameters. Therefore, we study the modified TP06 model by means of (numerical) bifurcation analysis using CL_MATCONT, a continuation toolbox for MATLAB [78, 79]. We consider an autonomous system of ordinary differential equations, where the right-hand side of this system depends on several state variables and parameters. Therefore, we consider the modified TP06 models, containing 18, 16, or 14 state variables, without initial stimulus for the bifurcation analysis, and we will focus our study on the effect of a deficit in the fast and slow potassium currents and , respectively, and an enhancement in the L-type calcium current .
The starting point of the bifurcation analysis is to determine an equilibrium of the modified autonomous system. To this end, the algebraic equations are solved thus:with , where represents the gating variables. The next step is to derive the stability of the equilibrium to determine the eigenvalues of the Jacobian or using the Routh–Hurwitz criterion [30, 82]. If we do so, we are able to determine the stability of the modified systems dependent on the ionic conductances , , and (Figure 4). Here, we have stable (black surfaces) and unstable regions (gray surfaces).
FIGURE 4
The unstable region allows the system to oscillate; after an initial stimulus the trajectory—the time-dependent voltage variation of the ODE model—can either develop normal APs or other oscillatory behavior such as EADs or some sort of ventricular tachycardia, which cannot be specified at this stage. However, in case the trajectory enters the stable region, it will reach a stable equilibrium—in this case not the resting potential, the potential becomes . This corresponds to sudden death, because if multiple cardiac cells lose their membrane potential, the heart can no longer generate electrical impulses, leading to asystole (cardiac arrest). Figure 4 also indicates that for increasing and decreasing and , the stable area becomes larger and the risk of sudden death increases, where the possibility that the voltage becomes continuously 0 is increased. This result is also compatible with the current state of knowledge. Moreover, Figure 4 also shows two red lines—the Andronov–Hopf bifurcation curves. At an Andronov-Hopf bifurcation, the system changes stability via a pair of purely imaginary eigenvalues— and —and a limit cycle appears.
Note that the standard value is . Therefore, in case does not increase too much and remains its standard value, the risk that the cardiac myocyte TP06 model develops sudden death is small to almost 0. However, if increases, such as by a factor 5, and is small, the risk of a sudden death increases and the M cells already develop EADs for the standard value of (Figure 3). Again, this is reasonable and compatible with the current state of knowledge since it is known that EADs can be related to the limit cycle bifurcating from an Andronov–Hopf bifurcation (e.g., [30, 71, 82]), and from a medical point of view, EADs can be precursors to sudden cardiac death.
Our goal is to examine the parameter set and to better understand the behavior of the (endocardial) cell model and compare the stationary dynamics of the 18-, 16-, and 14-dimensional model. A more practical way of analyzing a dynamic system by means of bifurcation theory is to use the continuation algorithm from [78, 79]. Choosing the parameter set and as a bifurcation parameter, we obtain two supercritical Andronov–Hopf bifurcations (Tables 2, 3) which both generate a stable limit cycle, and thus the bifurcating limit cycles are both stable. Between these Andronov–Hopf bifurcations, the system exhibits the stable equilibrium branch. Tables 2, 3 show that the three different models have almost identical Andronov–Hopf bifurcations.
TABLE 2
| value | Lyapunov exponent | |
|---|---|---|
| 18-dimensional | 0.027907858929580 | −2.683800872009436 |
| 16-dimensional | 0.027907858929578 | −2.683800872002161 |
| 14-dimensional | 0.027907284533354 | −2.667544014101886 |
Comparison: First supercritical Andronov–Hopf bifurcation.
TABLE 3
| value | Lyapunov exponent | |
|---|---|---|
| 18-dimensional | 0.071030406847997 | |
| 16-dimensional | 0.071030406847997 | |
| 14-dimensional | 0.071026913636226 |
Comparison: second supercritical Andronov–Hopf bifurcation.
We give all eigenvalues corresponding to the Andronov–Hopf bifurcations from Tables 2, 3 in Supplementary Tables S1, 2.
Starting a limit cycle continuation from the second supercritical Andronov–Hopf bifurcation, we derive a limit cycle branch containing a torus bifurcation and a (first) period-doubling bifurcation (cf. Table 4). Note that this limit cycle branch is the reason for the occurrence of both AP and EADs, which we will highlight later in more detail.
TABLE 4
| value of the torus bifurcation | value of the PD bifurcation | |
|---|---|---|
| 18-dimensional | 0.074450869484859 | |
| 16-dimensional | 0.074450869481487 | |
| 14-dimensional | 0.074412915016840 |
Comparison: torus bifurcation and period-doubling (PD) bifurcation of the first limit cycle branch.
Again, we give all characteristic multipliers corresponding to the period-doubling bifurcations from Table 4 in Supplementary Table S3. It is remarkable that, for 18, 16, and 14 dimensions, not only are the equilibrium curves identical, which is obvious, but the stability and the bifurcation points are also identical to a certain degree, which is visible in Tables 2–4. This is also reflected in Figure 5, where no are apparent differences. The first limit cycle branch (red surface) of the three simplified models bifurcates from a supercritical Andronov–Hopf bifurcation (red dot). At this bifurcation, the equilibrium curve (black line: the dashed part represents the unstable branch, the solid part the stable branch) also changes stability.
FIGURE 5
Note that fixing more gating variables equal to its steady-state solution will not change the equilibrium curve, but it will change its stability and potential bifurcations points and, therefore, the dynamics of the whole system. Furthermore, Figure 6 shows the limitation of the 18-dimensional model since for the M cells, the trajectory does not fit perfectly. Moreover, the stimuli for the 16- and 14-dimensional models differ from those in Figure 3. In any case, normal AP is well represented by the simplified models. In the case of more complex patterns such as EADs, the general dynamics of all models are identical, but finding the correct basin of attraction modifying the stimulus is more difficult.
FIGURE 6
Our next aim is to study the 14-dimensional epicardial TP06 model in more detail by means of bifurcation analysis. Our focus is on a situation like in Figure 6, where we use as the bifurcation parameter and consider value shifts and . We thus have a reduced rapid delayed rectifier current and an enhanced L-type calcium current. Utilizing the bifurcation analysis, we can extract the behavior and dynamics of the system. Starting a numerical continuation from a steady state solution, we find two supercritical Andronov–Hopf bifurcations (cf. Table 5).
TABLE 5
| Bifurcation points | value | Description | Period |
|---|---|---|---|
| Andronov–Hopf | Supercritical | ||
| Supercritical | |||
| Torus | |||
| Period-doubling | Subcritical | ||
| Subcritical | |||
| Subcritical |
Bifurcation for the 14-dimensional epicardial cell model with .
From these bifurcation points, a stable limit cycle branch bifurcates each. Following the limit cycle branch from the second Andronov–Hopf bifurcation , we find a torus bifurcation and a first period-doubling bifurcation. Furthermore, the limit cycle branch changes stability; finally, it will collide with the equilibrium curve and terminate there (cf. Figure 7 and Supplementary Figure S7).
FIGURE 7
This behavior is similar for all investigated combinations of . However, the position of the points varies (cf. Tables 5, 6). Starting a continuation from the first period-doubling bifurcation, we find a second limit cycle branch containing a second period-doubling bifurcation (cf. again Figure 7). Taking this approach further results in an unstable period-doubling cascade; in Figure 8, we zoom in on the region with the first four limit cycle branches.
TABLE 6
| bifurcation points | value | Description | Period |
|---|---|---|---|
| Andronov–Hopf | Supercritical | ||
| Supercritical | |||
| Torus | |||
| Period-doubling | Subcritical | ||
| Subcritical |
Bifurcation for the 14-dimensional epicardial cell model with .
FIGURE 8
To highlight the transition at each limit cycle bifurcation, Figure 9 shows four trajectories at the torus bifurcation and the first three period-doubling bifurcations (from left to right).
FIGURE 9
We start a simulation with an initial value on the limit cycle and with the corresponding value. Further details are in Supplementary Figures S1–S6.
Based on the bifurcation diagram, we can identify a region where EADs may appear. This region is clearly linked to the period-doubling cascade; however, whether EADs occur or not is also dependent on the initial values and the initial stimulus. Thus, if the initial values and stimulus are chosen in such a way that the trajectory is not able to enter the basin of attraction of the period-doubling cascade, then no EAD will appear, even though the system is in the dangerous region. This may happen, for instance, if the initial stimulus is high, which one may link to the usage of a pacemaker.
Finally, we provide for comparison the first five important bifurcations in the situation of the parameter shift (cf. Table 6).
3.3 Simulation of the monodomain model
The final focus of this study is to outline the synchronization and pattern formation of cardiac cells. Studying a monodomain model after analyzing its cellular model helps us understand how individual cellular electrophysiological behaviors scale and interact across tissue-level structures to simulate realistic electrical propagation in the heart. The simulations of the monodomain Equation 3 are done with FEniCS [73, 74], and the coupled ODE–PDE system is solved using cbcbeat (described in [75]) or fenics-beat [76]. Again, codes are provided in [42].
As in [39], we apply the S1–S2 protocol to generate spiral waves, meaning that we first apply to a thin strip on the left side of the square region. This induces a plane wave that propagates to the right edge of the square. Then, after , we apply a second stimulus to the lower half plane ( and , where denotes domain size).
We restrict ourselves to the epicardial human ventricular tissue model of [15] and to two different parameter settings: the normal AP setting and the parameter set . The last setting either means that we simulate the mid-myocardial model with reduced and enhanced or the epicardial model with reduced and and enhanced .
In Figure 10, we represent ten different time points, where each column corresponds to one time point, with spiral-wave pattern formations for a normal AP setting. The first row presents the normal AP for time points and , and the second for time points and . A stable spiral wave is thus initiated. However, in the case of EADs, we see a wave break-up, where the spiral wave disappears, the AP reaches its resting potential, and no further activity is recognized (cf. Figure 11).
FIGURE 10
FIGURE 11
4 Discussion
This paper aims to investigate the complex dynamics of the TP06 model [15, 40], such as the occurrence of action potentials and certain cardiac arrhythmia. While the focus is on early afterdepolarizations (EADs), the approach presented is also applicable to other diseases. A very efficient and beneficial method of studying the behavior of dynamical systems is delivered by bifurcation theory. Thus, one main goal of the paper is the numerical bifurcation analysis of the cellular TP06 model. To this end, we reported the difficulties in modeling the original model, which is also highlighted in [53, 71]. In addition, we discussed several possible modifications and model reductions of the TP06 model which allow numerical bifurcation analysis and apply the numerical continuation algorithm provided in [78, 79]. The aim of model reduction is to simplify a complex mathematical model while retaining its essential behaviors and properties so that the reduced model can be calculated more quickly and analyzed more easily. From an analytical point of view, it is desirable to derive a very simple model; indeed, certain details can be lost, such as on specific ion currents or concentrations, due to the focus on the voltage potential. One specific example is fixing the intracellular potassium concentration . A time-dependent variation of has only a small effect on the voltage in this model and causes problems for the bifurcation analysis; however, it might be of importance in certain applications. Thus, we present several steps/levels of model reduction and are able to reduce the 19-dimensional TP06 model to a 14-dimensional version, which has (almost) identical dynamics and trajectories as the original model (cf. Figure 5 and Tables 2–4). This model reduction allows not only performance of numerical bifurcation but also decreases the numerical effort and time. Apart from the steady-state dynamics, which are (almost) identical, it is necessary to adjust the initial stimulus such that the trajectory is able to enter the same basin of attraction. Moreover, Figure 6 illustrates that finding the correct basin of attraction is more difficult for complex dynamics and patterns. On the other hand, we can predict the dynamics of the original TP06 model extract very accurately by analyzing the reduced model. Thus, using the results from a mathematical model which is more feasible and low dimensional for the original model avoids the problem of possible loss of information as long as the approximation reflects the dynamics of the model sufficiently well—which is the case in our study.
After the remodeling, a systematical bifurcation analysis is provided using the set of parameters , which are known to induce EADs for certain combinations of values (e.g., [7, 8]). This is one aspect for focusing on EADs; another is they are complex voltage oscillations, and therefore this focus provides a better test for the accuracy of our model reduction. Thus, the bifurcation analysis can be seen not only as a tool for analyzing the dynamics of the system but can also be used to test the accuracy of our model reduction. The (numerical) bifurcation analysis shows that the oscillatory behavior of modified TP06 models is induced by an Andronov–Hopf bifurcation, and the occurrence of EADs is related to the existence of an unstable period-doubling bifurcation cascade. In addition, we can conclude oscillatory behavior for the TP06 model from the behavior of the modified models. Moreover, deterministic chaos is not observed, which is in line with [30, 71]. We show that even though the system is in a dangerous region and EADs potentially occur, the initial values and stimulus decide whether the trajectory enters the basin of transient attraction of the period-doubling bifurcation cascade. Therefore, we have three driving forces that induce EADs: the steady-state dynamics of the system related to a combination of a reduced slow and a reduced fast potassium current— and —with an enhanced calcium current , the initial values, and the initial stimulus. This knowledge can now be used to prevent the occurrence of EADs.
Furthermore, to outline the synchronization effects of cardiac cells, we performed numerical experiments highlight how networks of cells 2D-synchronize. Here, we can extract from our investigation how stable patterns appear (Figure 10). In addition, we showed that EADs may induce a wave break-up (Figure 11). While these results show a correlation between EADs and spiral wave breakup, the underlying mechanisms remain complex. A possible explanation is that the prolongation of AP, which often occurs in EADs, increases the heterogeneity of tissue cohesion and refractoriness. This may increase the spatial dispersion of repolarization and thus create the conditions for wavefront fragmentation. In addition, prolonged APs may shorten the effective wavelength of reentrant waves, increasing the chance of wavefront curvature and instability. It is also possible that EADs themselves act as transient depolarization triggers that disrupt the spiral core or interact with the refractory tail of the wavefront to promote breakup. Future research using targeted perturbations and parameter sensitivity analysis will be required to reveal the contributions of these different mechanisms. Moreover, it is highly important to analyze how and when EADs induce wave break-ups resulting in cardiac death. The heart (model) itself contains properties that allow it to stabilize and recover phenomena such as EADs. Here, it is not only the cellular behavior that plays a key role in the geometry of the heart (model) but also its discretization, its diffusivity, and the amount of affect cells that develop EADs.
The main focus of this study has been the model reduction of the cardiac cell model and providing an approach and codes to efficiently investigate cardiac dynamics numerically. These codes can be utilized to investigate different types of arrhythmia and cardiac dynamics on single-cell level. Moreover, it can be used for drug design or to make virtual experiments on the response of cardiac cells or the heart on external influences, such as a pacemaker. We showed on a tissue level that the local steady-state dynamics of the system include certain pattern formations. Furthermore, we discovered that the diffusivity of the system, the initial configuration, and the initial stimulus is highly important for pattern formation and network dynamics. This should be studied in more detail and is part of a future study.
Statements
Data availability statement
The raw data supporting the conclusions of this article will be made available by the author without undue reservation.
Author contributions
AE: writing – original draft and writing – review and editing.
Funding
The author(s) declare that financial support was received for the research and/or publication of this article. The author acknowledges the funding by the German Research Foundation (DFG) within the DFG Priority Program SPP 2171 Dynamic Wetting of Flexible, Adaptive, and Switchable Substrates through the project 422786086.
Acknowledgments
The author wishes to thank the referee for the careful reading of the original manuscript and the comments that eventually led to an improved presentation.
Conflict of interest
The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Generative AI statement
The author(s) declare that no Generative AI was used in the creation of this manuscript.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphy.2025.1569121/full#supplementary-material
Footnotes
1.^Mathematics Subject Classification: 92B25, 92C05, 37G15, 37M05, 37N25
References
1.
Tsaneva-AtanasovaKDiaz-ZuccariniV. Editorial: mathematics for healthcare as part of computational medicine. Front Physiol (2018) 9:985. 10.3389/fphys.2018.00985
2.
XiaLZhangHZhengD. Editorial: multi-scale computational cardiology. Front Physiol (2022) 13:847118. 10.3389/.fphys.2022.847118
3.
ErhardtAHTsaneva-AtanasovaKLinesGTMartensEA. Editorial: dynamical systems, pdes and networks for biomedical applications: mathematical modeling, analysis and simulations. Front Phys (2023) 10. 10.338910.3389/fphy.2022.1101756fphy.2022.1101756
4.
ErhardtAHPeschkaDDazziCSchmellerLPetersenAChecaSet alModeling cellular self-organization in strain-stiffening hydrogels. Comput Mech (2025) 75:875–96. 10.1007/s00466-024-02536-7
5.
CherryEMFentonFHKrogh-MadsenTLutherSParlitzU. Introduction to focus issue: complex cardiac dynamics. Chaos: Interdiscip J Nonlinear Sci (2017) 27:093701. 10.1063/1.5003940
6.
ErhardtAHTsaneva-AtanasovaKLinesGTMartensEA. In: Dynamical systems, PDEs and networks for biomedical applications: mathematical modeling, analysis and simulations. Lausanne: Frontiers Media SA (2023). 10.3389/978-2-83251-458-0
7.
Van NieuwenhuyseESeemannGPanfilovAVVandersickelN. Effects of early afterdepolarizations on excitation patterns in an accurate model of the human ventricles. PLOS ONE (2017) 12:e0188867–19. 10.1371/journal.pone.0188867
8.
ZimikSVandersickelNNayakARPanfilovAVPanditR. A comparative study of early afterdepolarization-mediated fibrillation in two mathematical models for human ventricular cells. PLOS ONE (2015) 10:e0130632–20. 10.1371/journal.pone.0130632
9.
TveitoAMardalKARognesME (2021). Modeling excitable tissue. Springer. 10.1007/978-3-030-61157-6
10.
MardalKARognesMEThompsonTBValnesLM (2022). Mathematical modeling of the human brain. Springer. 10.1007/978-3-030-95136-8
11.
NiedererSALumensJTrayanovaNA. Computational models in cardiology. Nat Rev Cardiol (2019) 16:100–11. 10.1038/s41569-018-0104-y
12.
SachseFBComputational cardiology: modeling of anatomy, electrophysiology, and mechanics, 2966. Springer Science and Business Media (2004).
13.
Bueno-OrovioACherryEMFentonFH. Minimal model for human ventricular action potentials in tissue. J Theor Biol (2008) 253:544–60. 10.1016/j.jtbi.2008.03.029
14.
ClaytonRBernusOCherryEDierckxHFentonFMirabellaLet alModels of cardiac tissue electrophysiology: progress, challenges and open questions. Prog Biophys Mol Biol (2011) 104:22–48.
15.
ten TusscherKHWJPanfilovAV. Alternans and spiral breakup in a human ventricular tissue model. Am J Physiology-Heart Circulatory Physiol (2006) 291:H1088–100. PMID: 16565318. 10.1152/ajpheart.00109.2006
16.
WeissJGarfinkelAKaragueuzianHChenPQuZ. Early afterdepolarizations and cardiac arrhythmias. Heart Rhythm (2010) 7:1891ff–9. 10.1016/j.hrthm.2010.09.017
17.
RodenDViswanathanP. Genetics of acquired long QT syndromeg. J Clin Invest (2005) 115:2025–32. 10.1172/.JCI25539
18.
FinkMNiedererSCherryEFentonFKoivumäkiJSeemannGet alCardiac cell modelling: observations from the heart of the cardiac physiome project. Prog Biophys Mol Biol (2011) 104:2–21. 10.1016/j.pbiomolbio.2010.03.002
19.
ClaytonRBernusOCherryEDierckxHFentonFMirabellaLet alModels of cardiac tissue electrophysiology: progress, challenges and open questions. Prog Biophys Mol Biol (2011) 104:22–48. 10.1016/j.pbiomolbio.2010.05.008
20.
KimSSatoD. Mathematical analysis of the role of heterogeneous distribution of excitable and non-excitable cells on early afterdepolarizations. Front Phys (2018) 6:117. 10.3389/fphy.2018.00117
21.
BarrioRMartínezMASerranoSPueyoE. Dynamical mechanism for generation of arrhythmogenic early afterdepolarizations in cardiac myocytes: insights from in silico electrophysiological models. Phys Rev E (2022) 106:024402. 10.1103/PhysRevE.106.024402
22.
KimreyJVoTBertramR. Canards underlie both electrical and ca2+-induced early afterdepolarizations in a model for cardiac myocytes. SIAM J Appl Dyn Syst (2022) 21:1059–91. 10.1137/22M147757X
23.
BarrioRJover-GaltierJAMartínezMPérezLSerranoS. Mathematical birth of early afterdepolarizations in a cardiomyocyte model. Math Biosciences (2023) 366:109088.
24.
TranDSatoDYochelisAWeissJGarfinkelAQuZ. Bifurcation and chaos in a model of cardiac early afterdepolarizations. Phys Rev Lett (2009) 102:258103. 10.1103/PhysRevLett.102.258103
25.
XieYIzuLBersDSatoD. Arrhythmogenic transient dynamics in cardiac myocytes. Biophys J (2014) 106:1391–7. 10.1016/j.bpj.2013.12.050
26.
ErhardtAH. Early afterdepolarisations induced by an enhancement in the calcium current. Processes (2019) 7:20. 10.3390/pr7010020
27.
ErhardtAHMardalKASchreinerJE. Dynamics of a neuron–glia system: the occurrence of seizures and the influence of electroconvulsive stimuli. J Comput Neurosci (2020) 48:229–51. 10.1007/s10827-020-00746-5
28.
KitajimaHYazawaTBarrioR. Fast–slow analysis and bifurcations in the generation of the early afterdepolarization phenomenon in a realistic mathematical human ventricular myocyte model. Chaos: Interdiscip J Nonlinear Sci (2024) 34:123108. 10.1063/5.0230834
29.
KuehnC. Multiple time scale dynamics. Appl Math Sci (2015) 191. Available online at: https://link.springer.com/book/10.1007/978-3-319-12316-5
30.
ErhardtAHSolemS. On complex dynamics in a purkinje and a ventricular cardiac cell model. Commun Nonlinear Sci Numer Simulation (2021) 93:105511. 10.1016/j.cnsns.2020.105511
31.
QuZ. Chaos in the genesis and maintenance of cardiac arrhythmias. Prog Biophys Mol Biol (2011) 105:247–57. 10.1016/j.pbiomolbio.2010.11.001
32.
DixMLilienkampTLutherSParlitzU. Influence of conduction heterogeneities on transient spatiotemporal chaos in cardiac excitable media. Phys Rev E (2024) 110:044207. 10.1103/PhysRevE.110.044207
33.
SuthDLutherSLilienkampT. Chaos control in cardiac dynamics: terminating chaotic states with local minima pacing. Front Netw Physiol (2024) 4:1401661. 10.3389/fnetp.2024.1401661
34.
JægerKHCharwatVCharrezBFinsbergHMaleckarMMWallSet alImproved computational identification of drug response using optical measurements of human stem cell derived cardiomyocytes in microphysiological systems. Front Pharmacol (2020) 10:1648. 10.3389/fphar.2019.01648
35.
TveitoAMaleckarMMLinesGT. Computing optimal properties of drugs using mathematical models of single channel dynamics. Comput Math Biophys6 (2018) 41–64. 10.1515/cmb-2018-0004
36.
ZhouXQuYPassiniEBueno-OrovioALiuYVargasHMet alBlinded in silico drug trial reveals the minimum set of ion channels for torsades de pointes risk assessment. Front Pharmacol (2020) 10:1643. 10.3389/fphar.2019.01643
37.
MargaraFWangZJLevrero-FlorencioFSantiagoAVázquezMBueno-OrovioAet alIn-silico human electro-mechanical ventricular modelling and simulation for drug-induced pro-arrhythmia and inotropic risk assessment. Prog Biophys Mol Biol (2021) 159:58–74. 10.1016/j.pbiomolbio.2020.06.007
38.
PassiniEZhouXTrovatoCBrittonOJBueno-OrovioARodriguezB. The virtual assay software for human in silico drug trials to augment drug cardiac testing. J Comput Sci (2021) 52:101202. 10.1016/j.jocs.2020.101202
39.
PravdinSFEpanchintsevTIPanfilovAV. Overdrive pacing of spiral waves in a model of human ventricular tissue. Scientific Rep (2020) 10:20632. 10.1038/s41598-020-77314-5
40.
ten TusscherKHWJNobleDNoblePJPanfilovAV. A model for human ventricular tissue. Am J Physiology-Heart Circulatory Physiol (2004) 286:H1573–89. PMID: 14656705. 10.1152/ajpheart.00794.2003
41.
HodgkinALHuxleyAF. A quantitative description of membrane current and its application to conduction and excitation in nerve. J Physiol (1952) 117:500–44. 10.1113/jphysiol.1952.sp004764
42.
ErhardtAH. Cardiac dynamics of a human ventricular tissue model with focus on early afterdepolarizations (2024). Available online at: https://github.com/andreerhardt/cardiac-dynamics-of-a-human-ventricular-tissue-model-with-focus-on-early-afterdepolarizations.
43.
SundnesJLinesGNielsenBMardalKATveitoA (2006). Computing the electrical activity in the heart. Springer. 10.1007/3-540-33437-8
44.
EllingsrudAJSolbråAEinevollGTHalnesGRognesME. Finite element simulation of ionic electrodiffusion in cellular geometries. Front Neuroinformatics (2020) 14:11. 10.3389/fninf.2020.00011
45.
BenedusiPEllingsrudAJHerlyngHRognesME. Scalable approximation and solvers for ionic electrodiffusion in cellular geometries. SIAM J Scientific Comput (2024) 46:B725–51. 10.1137/24M1644717
46.
TveitoAJægerKHKuchtaMMardalKARognesME. A cell-based framework for numerical modeling of electrical conduction in cardiac tissue. Front Phys (2017) 5. 10.3389/fphy.2017.00048
47.
JægerKHHustadKGCaiXTveitoA. Efficient numerical solution of the emi model representing the extracellular space (e), cell membrane (m) and intracellular space (i) of a collection of cardiac cells. Front Phys (2021) 8. 10.3389/fphy.2020.579461
48.
JægerKHTveitoA. Deriving the bidomain model of cardiac electrophysiology from a cell-based model; properties and comparisons. Front Physiol (2022) 12:811029. 10.3389/fphys.2021.811029
49.
MulimaniMKZimikSPanditR. An in silico study of electrophysiological parameters that affect the spiral-wave frequency in mathematical models for cardiac tissue. Front Phys (2022) 9. 10.3389/fphy.2021.819873
50.
LawsonBAdos SantosRWTurnerIWBueno-OrovioABurragePBurrageK. Homogenisation for the monodomain model in the presence of microscopic fibrotic structures. Commun Nonlinear Sci Numer Simulation (2023) 116:106794. 10.1016/j.cnsns.2022.106794
51.
AlonsoSBärMEchebarriaB. Nonlinear physics of electrical wave propagation in the heart: a review. Rep Prog Phys (2016) 79:096601. 10.1088/0034-4885/79/9/096601
52.
LuoCHRudyY. A model of the ventricular cardiac action potential. depolarization, repolarization, and their interaction. Circ Res (1991) 68:1501–26. 10.1161/01.RES.68.6.1501
53.
RosePSchleimerJHSchreiberS. Modifications of sodium channel voltage dependence induce arrhythmia-favouring dynamics of cardiac action potentials. PLOS ONE (2020) 15:e0236949–14. 10.1371/journal.pone.0236949
54.
KuznetsovYA. Elements of applied bifurcation theory. New York: Springer-Verlag (2004). 10.1007/.978-1-4757-3978-7
55.
FenichelN. Geometric singular perturbation theory for ordinary differential equations. J Differ Equ (1979) 31:53–98. 10.1016/0022-0396(79)90152-9
56.
TsumotoKShimamotoTAojiYHimenoYKudaYTanidaMet alTheoretical prediction of early afterdepolarization-evoked triggered activity formation initiating ventricular reentrant arrhythmias. Computer Methods Programs Biomed (2023) 240:107722. 10.1016/j.cmpb.2023.107722
57.
BarrioRJover-GaltierJMartínezMPérezLSerranoS. Mathematical birth of early afterdepolarizations in a cardiomyocyte model. Math Biosciences (2023) 366:109088. 10.1016/j.mbs.2023.109088
58.
TsumotoKKurataY. Bifurcations and proarrhythmic behaviors in cardiac electrical excitations. Biomolecules (2022) 12:459. 10.3390/biom12030459
59.
KurataYTsumotoKHayashiKHisatomeIKudaYTanidaM. Multiple dynamical mechanisms of phase-2 early afterdepolarizations in a human ventricular myocyte model: involvement of spontaneous sr ca2+release. Front Physiol (2020) 10:1545. 10.3389/fphys.2019.01545
60.
KitajimaHYazawaT. Bifurcation analysis on a generation of early afterdepolarization in a mathematical cardiac model. Int J Bifurcation Chaos (2021) 31:2150179. 10.1142/S0218127421501790
61.
TsumotoKKurataYFurutaniKKurachiY. Hysteretic dynamics of multi-stable early afterdepolarisations with repolarisation reserve attenuation: a potential dynamical mechanism for cardiac arrhythmias. Scientific Rep (2017) 7:10771. 10.1038/s41598-017-11355-1
62.
KimreyJVoTBertramR. Big ducks in the heart: canard analysis can explain large early afterdepolarizations in cardiomyocytes. SIAM J Appl Dynamical Syst (2020) 19:1701–35. 10.1137/19M1300777
63.
BarrioRMartínezMAPueyoESerranoS. Dynamical analysis of early afterdepolarization patterns in a biophysically detailed cardiac model. Chaos: Interdiscip J Nonlinear Sci (2021) 31:073137. 10.1063/5.0055965
64.
DiekmanCOWeiN. Circadian rhythms of early afterdepolarizations and ventricular arrhythmias in a cardiomyocyte model. Biophysical J (2021) 120:319–33. 10.1016/j.bpj.2020.11.2264
65.
OtteSBergSLutherSParlitzU. Bifurcations, chaos, and sensitivity to parameter variations in the sato cardiac cell model. Commun Nonlinear Sci Numer Simulation (2016) 37:265–81. 10.1016/.j.cnsns.2016.01.014
66.
PeirlinckMCostabalFYaoJGuccioneJMTripathySWangYet alPrecision medicine in human heart modeling. Biomech Model Mechanobiol (2021) 20:803–31. 10.1007/s10237-021-01421-z
67.
TrayanovaNABoylePM. Advances in modeling ventricular arrhythmias: from mechanisms to the clinic. WIREs Syst Biol Med (2014) 6:209–24. 10.1002/wsbm.1256
68.
TrayanovaNADoshiANPrakosaA. How personalized heart modeling can help treatment of lethal arrhythmias: a focus on ventricular tachycardia ablation strategies in post-infarction patients. WIREs Syst Biol Med (2020) 12:e1477. 10.1002/wsbm.1477
69.
DasíARoyASachettoRCampsJBueno-OrovioARodriguezB. In-silico drug trials for precision medicine in atrial fibrillation: from ionic mechanisms to electrocardiogram-based predictions in structurally-healthy human atria. Front Physiol (2022) 13:966046. 10.3389/fphys.2022.966046
70.
TrayanovaNALyonAShadeJHeijmanJ. Computational modeling of cardiac electrophysiology and arrhythmogenesis: toward clinical translation. Physiol Rev (2023) 104:1265–333. 10.1152/physrev.00017.2023
71.
ErhardtAHSolemS. Bifurcation analysis of a modified cardiac cell model. SIAM J Appl Dynamical Syst (2022) 21:231–47. 10.1137/21M1425359
72.
BernusOWildersRZemlinCWVerscheldeHPanfilovAV. A computationally efficient electrophysiological model of human ventricular cells. Am J Physiology-Heart Circulatory Physiol (2002) 282:H2296–308. 10.1152/ajpheart.00731.2001
73.
AlnaesMSBlechtaJHakeJJohanssonAKehletBLoggAet alThe FEniCS project version 1.5, 3. Archive of Numerical Software (2015). 10.11588/ans.2015.100.20553
74.
LoggAMardalKWellsGN. Automated solution of differential equations by the finite element method. Springer (2012). 10.1007/978-3-642-23099-8
75.
E RognesME FarrellPW FunkeSE HakeJMc MaleckarM. cbcbeat: an adjoint-enabled framework for computational cardiac electrophysiology. The J Open Source Softw (2017) 2:224. 10.21105/joss.00224
76.
FinsbergH. fenics-beat. zenodo (2024). 10.5281/zenodo.13323643
77.
RushSLarsenH. A practical algorithm for solving dynamic membrane equations. IEEE Trans Biomed Eng, 25 (1978) 389–92. 10.1109/TBME.1978.326270
78.
DhoogeAGovaertsWKuznetsovYA. Matcont: a matlab package for numerical bifurcation analysis of odes. ACM Trans Math Softw (2003) 29:141–64. 10.1145/779359.779362
79.
DhoogeAGovaertsWKuznetsovYAMeijerHGESautoisB. New features of the software matcont for bifurcation analysis of dynamical systems. Math Comput Model Dyn Syst (2008) 14:147–75. 10.1080/13873950701742754
80.
HundTJKuceraJPOtaniNFRudyY. Ionic charge conservation and long-term steady state in the Luo–rudy dynamic cell model. Biophysical J (2001) 81:3324–31. 10.1016/s0006-3495(01)75965-6
81.
TusscherKHWJTPanfilovAV. Cell model for efficient simulation of wave propagation in human ventricular tissue under normal and pathological conditions. Phys Med and Biol (2006) 51:6141–56. 10.1088/0031-9155/51/23/014
82.
ErhardtAH. Bifurcation analysis of a certain Hodgkin–Huxley model depending on multiple bifurcation parameters. Mathematics (2018) 6:103–15. 10.3390/math6060103
Summary
Keywords
cardiac dynamics, early afterdepolarization, model reduction, bifurcation analysis, monodomain equation, pattern formation
Citation
Erhardt AH (2025) Cardiac dynamics of a human ventricular tissue model with a focus on early afterdepolarizations. Front. Phys. 13:1569121. doi: 10.3389/fphy.2025.1569121
Received
31 January 2025
Accepted
13 May 2025
Published
30 July 2025
Volume
13 - 2025
Edited by
Aslak Tveito, Simula Research Laboratory, Norway
Reviewed by
Alexandre Lewalle, Imperial College London, United Kingdom
Kunichika Tsumoto, Kanazawa Medical University, Japan
Updates
Copyright
© 2025 Erhardt.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: André H. Erhardt, andre.erhardt@wias-berlin.de
ORCID: André H. Erhardt, orcid.org/0000-0003-4389-8554
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.