On complex dynamics in a Purkinje and a ventricular cardiac cell model

Cardiac muscle cells can exhibit complex patterns including irregular behaviour such as chaos or (chaotic) early afterdepolarisations (EADs), which can lead to sudden cardiac death. Suitable mathematical models and their analysis help to predict the occurrence of such phenomena and to decode their mechanisms. The focus of this paper is the investigation of dynamics of cardiac muscle cells described by systems of ordinary differential equations. This is generically performed by studying a Purkinje cell model and a modified ventricular cell model. We find chaotic dynamics with respect to the leak current in the Purkinje cell model, and EADs and chaos with respect to a reduced fast potassium current and an enhanced calcium current in the ventricular cell model -- features that have been experimentally observed and are known to exist in some models, but are new to the models under present consideration. We also investigate the related monodomain models of both systems to study synchronisation and the behaviour of the cells on macro-scale in connection with the discovered features. The models show qualitatively the same behaviour to what has been experimentally observed. However, for certain parameter settings the dynamics occur within a non-physiological range.


Introduction
Nowadays, mathematical modelling and numerical simulations are essential and standard approaches to study and analyse real world problems and phenomena in life science. One major aim is the understanding of complex dynamics and behaviour of these systems. For this purpose, bifurcation theory has proven to be a very helpful and powerful tool in order to investigate dynamical systems and their (complex) dynamics, see [1,2,3,4,5] for an overview. Furthermore, numerical bifurcation analysis has become a profitable tool in the study of (for instance) climate, neuronal and cardiac models.
Cardiac muscle cells can exhibit complex patterns of oscillations like spiking and bursting, which is related to ion current interactions of the considered cell. Aside from normal action potentials of a cardiac muscle cell, certain kinds of cardiac arrhythmia can occur. This includes specific Email addresses: andreerh@math.uio.no (André H. Erhardt), susanne.solem@ntnu.no (Susanne Solem) 1 The authors have contributed equally to the content of this manuscript. types of abnormal heart rhythms, which can lead to sudden cardiac death. In addition, irregular behaviour, such as (deterministic) chaos or chaotic early afterdepolarisations, has been observed in experimental as well as in computational studies, see [6,7] and the references therein. It is therefore highly interesting and important to understand the complex behaviour and mechanism of such biological phenomena. Moreover, cardiac dynamics or heart rhythms can be quite sensitive to the influence of certain drugs, which has been investigated experimentally but also in computational studies, see e.g. [8,9,10]. In later years, the focus has been to continuously move towards interdisciplinary research, including biology, computer science, and mathematics, to tackle these issues. As a consequence, the number of existing mathematical models based on experimental data is also continuously increasing.
The development of a good and precise mathematical model is essential to design numerical experiments for the study of cardiac dynamics, but also for the investigation of the influence of certain external effects such as drugs or oxidative stress. To this end, mathematical analysis is key to decode occurring phenomena and to validate a derived model in all details. The newly gained information of the considered model, can then be either used to improve the model or to proceed with the original aspiration, e.g. the investigation of optimal properties of drugs [9].
To this end, bifurcation theory has been utilised to investigate the dynamics of cardiac muscle cells in recent years, see e.g. [11,12,13,14]. Continuing on this line of research, this paper highlights how useful bifurcation theory can be for the understanding of complex cardiac dynamics and how it can be applied to find hidden features and dynamics of cardiac single cell models described by ordinary differential equations (ODEs) through numerical investigations.
In the end, it is the synchronisation of a large group of cells that decides whether a cardiac arrhythmia spreads or dies out. For this reason, a brief study of how the micro-scale single cell features of these models affect the behaviour of an ensemble of cells at the macro-scale level (cm) is provided. This is done by an up-scaling of the ODE system to a PDE-ODE monodomain model. All of the above will be done by an in-depth mathematical and numerical investigation of the two cardiac cell models introduced in [15,16], where one is a model of a Purkinje cell, and the other a model of a human ventricular cell. In particular, we find new features of the considered models, such as chaos and early afterdepolarisations, and investigate how this affects groups of cells at the tissue level.
We find chaotic dynamics in both models considered. Similar chaotic dynamics to what we discover can be observed in experiments, cf. [6]. However, the dynamics seems to appear in a non-physiological range. This either requires the improvement of the models or the experimental validation.
In addition, we discover EADs in the human ventricular cell model. This behaviour does seem to be within the physiological range [17], which is a validation of the model in this case.
Finally, we show that the (in)validity of the models in terms of being within the physiological range carries over to the synchronisation effects in the corresponding monodomain models.
These findings clearly highlight the advantages of bifurcation theory in the analysis of cardiac muscle cell dynamics by detecting unexpected or maybe non-physiological behaviour of the model.
Outline of the paper. In Section 2 we start with a mathematical and biological description of the models and problems under investigation. A brief background on the modelling of cardiac muscle cells and on up-scaling to a monodomain model at the tissue level (cf. [18,19,20]) is provided. Furthermore, we perform a stability analysis of the ODE system from [15], and show how to extend this analysis to the macro-scale monodomain model. The approach is explained in detail for the four dimensional model [15]. However, the ansatz can also be used for more complex models, see [21,22].
The structure of the systems in [15,16] is similar. Nevertheless, the behaviour and dynamics that they display can be quite disparate due to different complexity and parameter settings. In Section 3, we apply a numerical bifurcation analysis in order to derive a complete understanding of the dynamics of the model from [15] with respect to certain parameters. The analysis is then extended to the corresponding macro-scale monodomain model. Based on the results in Section 3, we then continue our analysis by studying a ten dimensional version of the model from [16] in Section 4. Finally, we close our paper with a discussion in Section 5 and then a conclusion.

Biological and mathematical background
The history of mathematical modelling of action potentials (APs) of excitable biological cells like neurons and cardiac cells starts with the famous and pioneering Hodgkin-Huxley (HH) model from 1952 [23]. In [23], the authors established a mathematical approach that can be used to model APs of excitable biological cells by a system of ODEs. The first model of a cardiac cell is the Noble model from 1962 [15] of a generic Purkinje cell. In 1991, Luo and Rudy published an ionic model for cardiac action potential in guinea pig ventricular cells [24]. In the last decades, there has been an immense development in the modelling of cardiac muscle cells, see e.g. [25,26,27,28]. These conductance-based models represent a minimal biophysical interpretation of an excitable biological cell in which current flow across the membrane is due to charging of the membrane capacitance and movement of ions across ion channels, cf. Figure 1. Ion channels are selective for particular ionic species. In general, an AP is a temporary, characteristic variance in the membrane potential of an excitable biological cell from its resting potential. The molecular mechanism of an AP is based on the interaction of voltage-sensitive ion channels. The reason for the formation and the special properties of the AP is established in the properties of different groups of ion channels in the plasma membrane. An initial stimulus activates the ion channels as soon as a certain threshold potential is reached. Then, these ion channels break open and/or up allowing an ion current flow, which changes the membrane potential. A normal AP is always uniform and the cardiac muscle cell AP is typically divided into four phases: the resting phase, the upstroke phase, the (long) plateau phase and the repolarisation phase. This mechanism is based on several different currents. One example is the potassium current I K which is usually divided into a fast (I Kr ) and a slow current (I Ks ), cf. scheme in Figure 1(a).
This electrophysiological behaviour can be described by the ordinary differential equation: where V denotes the voltage (in mV ) and t the time (in ms), while I ion is the sum of all transmembrane ionic currents. I stimulus represents the externally applied stimulus and C m denotes membrane capacitance. The model in [15] contains a single potassium current I K , while the authors of [16] merged a fast (I Kr ) and a slow current (I Ks ) to derive their model, which is based on the system in [29]. In this paper, we will slightly modify the model from [16], i.e. we replace I K by the currents I Kr and I Ks from [29]. Furthermore, the different ion currents may depend on different gating variables, individual ionic conductances G i and Nernst potentials E i , i = Na, K, Ca, etc., cf. Section 2.1.
cytoplasm SR (a) Scheme of a cardiac muscle cell.
Physical system of a cardiac muscle cell. Figure 1: (a) Scheme of a cardiac muscle cell [26], where SR denotes the sarcoplasmic reticulum, INaCa = Na + /Ca 2+ exchanger current and INaK = Na + /K + pump current. (b) Physical system of a cardiac muscle cell. The dashed lines denote the ion currents, which are not included due to the lack of space.
We want to highlight that cardiac cell models usually have different time scales and may exhibit so-called mixed-mode oscillations [30], cf. [31], and/or chaotic behaviour, cf. [12,32,33,34,35], which can be linked to certain cardiac arrhythmia.
For instance, if there are depolarising variations of the membrane voltage, we are speaking about afterdepolarisations (ADs). These ADs are divided into early (EADs) and delayed afterdepolarisations (DADs). This division depends on the timing obtaining the AP. EADs occur either in the plateau or in the repolarisation phase of the AP and are benefited by an elongation of the AP, while DADs occur after the repolarisation phase is completed. EADs are resulting, for example, from a reduction of the repolarising K + currents or an enhancement in Ca 2+ currents, see e.g. [33]. Triggers for this are congenital disorders of ion channels or the ingestion of medicaments. In general, EADs are additional small amplitude spikes (mathematically speaking mixed-mode oscillations), i.e. pathological voltage oscillations, during the plateau or repolarisation phase. They are caused by ion channel diseases, oxidative stress or drugs. Furthermore, the presence of EADs strongly correlates with the onset of dangerous cardiac arrhythmias, including torsades de pointes (TdP), which is a specific type of abnormal heart rhythm that can lead to sudden cardiac death, see [36,37,17]. Thus, it is highly important to understand the complex behaviour of such biological phenomena [38].
Finally, we want to point out that not all of the existing cardiac cell models may include complex dynamics such as EADs or chaos.

A Purkinje cardiac cell model
First, we focus on the model from [15], which reads as follows: with the membrane capacitance C m = 12 µF cm 2 and the ion currents I Na (sodium), I K (potasssium), and I L (leak current), described by I Na = (G Na m 3 h + 0.14)(V − E Na ), and I L = G L (V − E L ), respectively. The individual ionic conductances are given by G K 1 = G K 2 = 1.2 mS cm 2 , G Na = 400 mS cm 2 and G L = 0.075 mS cm 2 , while the Nernst potentials are given by E K = −100 mV , E Na = 40 mV and E L = −60 mV . Furthermore, the different gating variables m, h and n satisfy the differential equation where y represents the gating variables m, h and n, while y ∞ := y ∞ (V ) = ay ay+by denotes the equilibrium of the gating variable y and τ y := τ y (V ) = 1 ay+by its relaxation time constant with Notice that the gating variables are important for the activation and inactivation of the different ion currents, which is needed for the ion current flows and the resulting action potential. The Noble model (1) describes the long lasting action and pace-maker potentials of the Purkinje fibres of the heart based on the Hodgkin-Huxley formulation [23]. While the sodium current is very similar to the one from [23], the potassium current differs in its formulation and the calcium current is missing. Nevertheless, the solutions to system (1) closely resembles the Purkinje fibre action and pace-maker potentials. It is shown that its behaviour in response to applied currents and to changes in ionic permeability corresponds fairly well with that observed experimentally. Furthermore, the Noble model (1) is one of the earliest mathematical models of cardiac APs which is able to produce a certain type of cardiac arrhythmia, so-called alternans in AP duration (APD). These alternans in APD is a result of a period doubling bifurcation with respect to the cycle length, see [39,40]. A period doubling bifurcation is a creation or destruction of a periodic orbit with double the period of the original orbit.
In [15] the author numerically studied the influence of the conductance G L of the leak current on the trajectory of the Noble model (1). We will extend this numerical study by using bifurcation analysis to derive a more detailed insight into the behaviour of the solutions to (1) by varying G L , see Section 3.1.

Stability analysis
The steady state or equilibrium of system (1) is determined by h ≡ h ∞ (V ), m ≡ m ∞ (V ), n ≡ n ∞ (V ) and solving the algebraic equation This yields where X ∞ is depending on several system parameters, cf. system (1). Furthermore, the Jacobian of the right hand side of system (1) evaluated at X ∞ is given by where we used the fact that Note that the location and stability of X ∞ is depending on the different system parameters, e.g. varying G L changes the location and the stability of X ∞ , while varying C m changes only the stability. To determine the stability of the equilibrium one has to calculate the solution(s) of the characteristic polynomial i.e. the eigenvalues of J , where In addition, the Routh-Hurwitz criterion implies that all characteristic exponents λ i , i = 1, . . . , 4 have negative real parts if and only if the conditions [3,4,41]. Furthermore, if all Hurwitz minors satisfy ∆ i > 0 for i = 1, · · · , n − 1 and ∆ n = 0, where n denotes the dimension of the system, we know that the system exhibits an Andronov-Hopf bifurcation. Using ∆ 4 = 0, we get i.e. the equilibrium has a pair of purely imaginary eigenvalues λ 1,2 = iω 0 with ω 0 = a 3 a 1 > 0 and two further eigenvalues In general, an Andronov-Hopf bifurcation corresponds to the birth of a limit cycle, when the equilibrium changes stability via a pair of purely imaginary eigenvalues. Usually, an Andronov-Hopf bifurcation is considered as a trigger to oscillatory behaviour in dynamical systems and may cause normal AP and cardiac arrhythmia in a cardiac cell model. In case that a n = 0, then the system exhibits a fold or saddle-node or limit point bifurcation, i.e. the equilibrium has at least one eigenvalue which is equal to zero. A limit point bifurcation is a collision and disappearance of two equilibria in dynamical systems, which is a turning point.

A monodomain model of the Purkinje cardiac cell model
Beside the study of cardiac cell models an important focus is the behaviour and dynamics of several cells, i.e. the dynamics on a macro-scale (cm), where many cells are connected together and interacting with each other. To this end, we consider the following monodomain model, i.e. extension of the ODE model (1) to a PDE-ODE model including an additional diffusion term: in Ω, where Ω is a bounded domain, M i denotes the intracellular conductivity tensor, λ the extra-to intracellular conductivity ratio, χ is the membrane surface area per unit volume and ν is the unit normal, cf. [42,43,44]. System (6) is a monodomain model, meaning that equal anisotropy rates, i.e. M e = λM i in mS cm , are assumed in the more complex bidomain model. Here, λ ∈ R is constant and M e denotes the intracellular conductivity tensor. Furthermore, we use λ 1+λ M i χ = 1 360 mS, which is the diffusion constant originally used in [16], unless otherwise stated.
For a better understanding of the behaviour of cells on a macro-scale level, we follow the approach from [45,46,47] to derive a linearised system of model (6). First of all, we know that on rectangle-like domains Ω := [0, 1 ] × · · · × [0, n ] ⊂ R n the eigenvalues and eigenfunctions of the Neumann problem see [48]. In general, the eigenvalue of the spectral problem (7) can be derived by multiplying the first equation of system (7) by u k , integrating over Ω, using Green's formula and applying the boundary condition. Then, one gets for the Neumann problem (7) the following eigenvalues: . Now, we consider the linearised system of the monodomain equation (6) around an equilibrium X ∞ of the ODE system (3). It has the form where D denotes the 4 × 4 diffusion matrix with almost everywhere zero entries except the first one, which is while J is the Jacobian as stated in (4). Define the linear operator L as Then, consider the following characteristic equation where (ψ 1 , · · · , ψ 4 ) T is the eigenfunction of L corresponding to the eigenvalue µ. Thus, let where V k , h k , m k and n k are time-dependent coefficients. We can then conclude that where Keep in mind that system (7) has eigenvalues Hence, the linearised system (8) has infinitely many Jacobians J k . Furthermore, we have linearised our PDE-ODE system (6) around an equilibria of the ODE system (1), such that there is an underlying assumption of the solutions of system (6) being close to the equilibria of system (1). We continue by deriving an ODE model to represent the behaviour of the linearised system (8) for each mode k = 0, 1, 2, 3, . . . in close proximity to the equilibrium V ∞ . This is done by ensuring that the resulting system has the same equilibrium as the Noble model (1) and the Jacobian J k .
Then, we are in a situation where we can analyse the behaviour of a single cell on a macro-scale including the influence of the diffusion term (dependent on k) of the monodomain model (6). This allows us to gain a better understanding of the cell interactions. The resulting ODE system reads as follows: where V ∞ is an equilibria for equation (3). We have designed the location of the equilibrium of system (10) to be the same as for the Noble model (1). However, considering the stability analysis from Section 2.1, the stability of the equilibrium may be different. The first entry of the Jacobian J in (4) is replaced by J k 11 . Thus, also the coefficients a j , j = 1, . . . , 4 of the characteristic polynomial (5) are changed. Obviously, this will affect the stability of the system (10) dependent on the parameters λ, M i , χ, and k. This indicates that the cellular behaviour of model (1) is not (necessarily) one-to-one transferred to the behaviour and dynamics of the monodomain equation (6). Nevertheless, it is a good starting point to study the dynamics of a single cell before extending the analysis to the macro-scale. Do however note that the mode k = 0 does give us the same dynamics as that of the ODE system (1).
In the discrete setting, if we choose k = 1 h 2 , where h denotes the grid size, and k = 1 h 2 is the biggest eigenvalue for the discrete laplacian, the term in (10) tends to −∞ and blows up as h → 0, which should stabilise (10) for large modes k (or refined grid size h). Looking at the coefficients a j , j = 1, . . . , 4 of the characteristic polynomial (5), we can see that a 1 , a 2 and a 3 will tend to ∞, while a 4 will tend to −∞ as h tends to zero and k tends to ∞. This implies that 1 a 4 > 0 and thus, the numerics of the PDE will stabilise for decreasing grid size, which is to be expected.
The same analysis can be used on 2D domains Ω = [0, ] × [0, * ], , * > 0, by considering the slightly different term A similar modification applies to a 3D cube or cuboid. For a more general geometry one can expect that the additional term deriving from the linearisation of the monodomain model (6) is more complicated. However, the general discussion above is expected to remain valid.

A ventricular cardiac cell model
As mentioned, we first investigate the Noble model (1). Indeed, we will see that this model has some limitations. Therefore, we will also study a human ventricular cell model, which is more advanced due to the number of included ion currents. The description of the system in [16] for epicardial cells is similar to system (1), but it contains more ion currents and reads as follows: where I stimulus denotes an external stimulus and I ion = I Cab + I NaCa + I Nab + I Ca + I K + I NaK + I Na + I K1 + I to is depending on the fast sodium current I Na = G Na m 3 v 2 (V − E Na ), the slow calcium current In [16] the authors studied system (11) with a delayed rectifier current I K = G K x 2 (V − E K ), while we are considering the delayed rectifier current I K = I Kr + I Ks from [29]. Here, the current is the slowly activating current. Including the fast and slow potassium current will makes the dynamics more realistic.
Furthermore, system (11) contains the background currents I Cab and I Nab , the Na + /Ca 2+ exchanger current I NaCa and the Na + /K + pump current I NaK , cf. Figure 1(a). Notice that the system is depending on 9 gating variables, i.e. m, v, d, f , r, to, x r , x s and K1, satisfying the differential equation (2), where d, r and K1 are assumed to be equal to their steady states. We will consider all gating variables as state variables, therefore the dimension of the system is 10.
Moreover, we use C m = 1.534 µF cm 2 , cf. [29], while the individual ionic conductances are given by G Na = 16 mS cm 2 , G Ca = 0.064 mS cm 2 , G to = 0.3 mS cm 2 , G Kr = 0.015 mS cm 2 , G Ks = 0.02 mS cm 2 and G K1 = 2.5 mS cm 2 . The equilibria and the Jacobian are similarly determined as for the Noble model (1). The only difference is that we have 6 ODEs more to consider. This increases the 4 × 4 matrix J to a 10 × 10 matrix. Following the same approach as in Section 2.2, we can derive from the monodomain model related to system (11) with I stimulus = 40 µA cm 2 and a duration of 2 seconds, i.e.
in Ω with Neumann boundary condition ν · (M i ∇V ) = 0 on ∂Ω, the following ODE system d dt where V ∞ the equilibrium of the voltage V of system (11). Note that the stability analysis for system (10) and the previous discussion also holds true for system (13).

Dynamics of the Noble model
In [15], the author mentioned a change in the dynamics of system (1) by varying the leak conductance G L . We will show that changing G L has influence on the period of the AP as well as if the system converges into a stable equilibrium or not. In Figure 2, the trajectory of system (1) is given for two different values of the leak conductance G L , i.e. G L = 0.075 mS cm 2 and G L = 0.18 mS cm 2 . In Figure 2(a), system (1) reveals a periodic trajectory representing two APs of a cardiac single cell, while in Figure 2(b) the trajectory needs certain amount of time to reach a stable periodic pattern also representing APs. This shows that system (1) is sensitive with respect to G L , which may have an influence on the amplitude and the period T of the trajectory given by V (t + T ) − V (t) = 0. In addition, the initial value also has an influence on the appearing pattern (at least locally).  3.1. Bifurcation analysis of system (1) with respect to G L The behaviour described above, we now investigate systematically utilising numerical bifurcation analysis as done in [31,52]. Our first step is to determine a bifurcation diagram with respect to the leak current, i.e. we choose the leak conductance G L as the bifurcation parameter. Notice that we will only consider positive values of G L , since they are the physiologically relevant ones. Mathematically and numerically we can easily extend the bifurcation diagram also to negative values of G L . We start by determining the equilibrium curve, i.e. we calculate for different values of G L the corresponding equilibria, which gives us the desired equilibrium curve. Moreover, we determine the stability of each equilibria. The equilibrium curve we can easily calculate with a self-written routine. However, to derive a detailed enough bifurcation diagram efficiently, it is advisable to use a robust continuation algorithm as the one mentioned in Section 2.4.
The bifurcation diagram of system (1) related to G L exhibits an unstable (black dashed line) and a stable (black solid line) equilibrium branch, see negative first Lyapunov coefficient. From the supercritical Andronov-Hopf bifurcation a stable limit cycle branch (solid blue line) bifurcates, which becomes unstable (dashed blue line) via a period doubling bifurcation (solid red square, G L ≈ 0.187785 mS cm 2 ). Then, the limit cycle branch gains stability again via a limit point of cycle bifurcation (solid green square, G L ≈ 0.193546 mS cm 2 ). Furthermore, this limit cycle branch contains a second period doubling bifurcation, which is also connected to the first one via a second limit cycle branch, cf. Figure 3 and Figure 4.
The bifurcation diagram in Figure 3 only includes the first two limit cycle branches, as including further details would not be particularly visible. Indeed, further branches exist as we will see below.
Together with Figure 4 and Figure 5, a nice graphical explanation of the observations in [15] appear. We can identify values of G L for which system (1) oscillates or has a stable equilibrium. Even more, we can determine the maximal amplitude of an oscillation corresponding to G L . In addition, we also get the corresponding period of each limit cycle. We have that the period T for G L = 0.075 mS cm 2 is approximately 564.1345ms, while for G L = 0 mS cm 2 we have T ≈ 839.5015ms and for G L = 0.18 mS cm 2 we have T ≈ 324.2749ms. Thus, we know that the pattern of the trajectory repeats faster if G L increases before it converges into a stable equilibrium for G L too large, which explains the observations from [15]. Furthermore, from Figure 3 we observe that the maximal amplitude is decreasing while G L increases.
In Figure 4 a third limit cycle branch is also included and we focus on a smaller range of G L values, where interesting dynamics may appear. The next additional limit cycle branches behave very similarly to the third one, i.e. they are bifurcating from a period doubling bifurcation, they lose stability via a further period doubling bifurcation, and stay unstable until they converge into a period doubling bifurcation of the previous limit cycle branch. Only the second limit cycle branch behaves slightly different -it bifurcates from the first period doubling bifurcation (G L ≈ 0.187785 mS cm 2 ), becomes unstable via a limit point of cycle bifurcation, gains stability again via a second limit point of cycle bifurcation and finally, loses stability via a further period doubling bifurcation (G L ≈ 0.187308 mS cm 2 ). Then, it stays unstable until it converges into a period doubling bifurcation of the first limit cycle branch, cf.  Starting from the second limit cycle branch, system (1) exhibits a stable period doubling cascade, which is usually a route to chaos, cf. e.g. [35]. However, in the bifurcation diagram in Figure 5, where several trajectories (red lines) are highlighted for different values of G L , no irregular or chaotic behaviour nor additional oscillations can be seen.
The reason is quite simple: The system exhibits a stable limit cycle branch from the first limit cycle branch, which "surrounds" the interesting dynamics of the Noble model (1), cf. Figure 3 and Figure 4. Hence, for the standard initial condition we do not get any sudden change in the dynamics. The trajectory will jump on to the "outer" stable limit cycle branch and stays there. However, if we choose an initial value from inside our limit cycle branches, e.g. one of the unstable equilibria, we will have a sudden change in the dynamics of the system close to the period doubling bifurcations and the period doubling cascade. For instance, Figure 6 presents a simulation of system (1) over 150s, which shows an irregular and chaotic behaviour. Hence, not only the choice of bifurcation parameter is essential, but also the choice of initial values, cf. Figure 7. 3.2. Bifurcation analysis with respect to the potassium I K current Next, we analyse the dynamics of system (1) with respect to the potassium I K current. To this end, we choose one of the two potassium conductances G K 1 and G K 2 . Notice that usually the I K current depends on one gating variable, cf. [14,16,24] or one has a splitting into a fast I Kr and a slow I Ks current, cf. [26,27,29]. Here the situation is different, since I K is modelled as the sum of two currents, which is not split into a fast I Kr and slow I Ks part, nor into a potassium current and a background current, cf. [16,26,29]. Nevertheless, we will study the influence of a deficit in the potassium current in (1).
It is well known that a deficit in the potassium current may induce EADs. This was among others verified for simplistic models in [12,14]. Using bifurcation analysis we can conclude from Figure 8 that no EADs appear via a deficit in G K 2 . In the case where we combine G K 1 and G K 2 , i.e. consider only one potassium conductance G K = G K 1 = G K 2 , system (1) behaves similarly. Figure 9 shows that the behaviour of system (1) does not change dramatically under variations of G K 1 compared to variations of G K 2 . We see that the model (1) does not exhibit EADs via a reduced potassium current. This was also checked for G L between 0 mS cm 2 and 0.2 mS cm 2 . One explanation is that the calcium current I Ca has an important influence on the behaviour of a cardiac cell, cf. [52,31]. This current is missing in the Noble model (1). For each additional current one gets additional ion current interactions, and one gains additional system parameters influencing the dynamics of the system. From [52] it is known that the ion current interaction between the potassium and the calcium current is important for the occurrence of EADs. This might be a reason for the appearance of different behaviour from what one would expect in a cardiac muscle cell.  In conclusion, we see that we need a more detailed model to study more diverse behaviour (including EADs), see e.g. [53]. The above investigation shows that Hodgkin-Huxley (type) models are sensitive with respect to their system parameters, and also with respect to their initial values. Thus, to capture all dynamics it is essential to consider both the system parameters and the initial values. The challenge is to derive a model which exhibits all dynamics of interest, and which at the same time is simple enough to allow its behaviour to be studied efficiently. This study indicates that model (1) exhibits physiologically relevant APs or very fast oscillations with small amplitudes with respect to the potassium current I K . The system also exhibits chaotic behaviour with respect to the leak current I L . However, the voltage is always negative, see Figure 6, and this seems to be non-physiological [6].
In Section 4 we will see that by slightly increasing the complexity of the model considered, more of the expected dynamics appear.

Effects on the macro-scale
In Section 2.2 the analysis of the linearised system (10) showed that one cannot expect that the cellular behaviour of a single cell model is one-to-one transferred to the behaviour and dynamics of the corresponding monodomain equation. Therefore, we briefly visualise how the interaction of an ensemble of cells belonging to two different regimes might play out in this section. Based on the discussion in the last two paragraphs of Section 2.2, we focus on a 1 cm one-dimensional cable, i.e. x ∈ [0, 1], for simplicity. To see the additional effects of the cell interactions at the macro-scale, we set the initial condition to partly belong to the chaotic regime of the Noble ODE model (1), and partly to the stable regime: where D ∈ [0, 1]. Note that this is not an equilibrium of the monodomain model (6). We fix G L = 0.1845 mS cm 2 to allow for both chaotic and stable behaviour. Figure 10 shows that the stable behaviour suppresses the chaotic one and a travelling wave dynamic appear as time progresses. Increasing the initial chaotic domain slightly to D = [0, √ 2 − 0.99), one can observe from Figure 11 that the chaotic behaviour prevails. Moving the chaotic regime to the middle of the interval, the system can tolerate a much larger area of initially chaotic cells, see Figure 12(a). Here the initial data is set to (14) with D = (0.1, 0.5) ∪ (0.51, 0.9). However, reducing the diffusion parameter, the behaviour of the system turns chaotic, see Figure 12(b).
One can see from above that there is a critical mass of initially chaotic cells, depending on the diffusion parameter and the placement of the initially chaotic region, for inducing chaotic behaviour along the whole cable. These results are in concert with earlier observations noting that enough cells have to be triggered for chaotic behaviour to prevail [33,34]. Furthermore, although induced differently, similar behaviour to what we see in Figure 12 has been observed in experiments [6]. However, as noted in the previous section, the voltages seem to be within a non-physiological range.

Dynamics of the modified Bernus model
In this section we investigate the slightly modified version of the human ventricular cardiac cell model from [16], i.e. system (11). As described in Section 2.3, this model contains more ion currents, pumps and exchangers compared to model (1), including the missing calcium current I Ca and the fast and slow potassium current, I Kr and I Ks . Therefore, one may expect that this model is more realistic and exhibits different and diverse dynamics. System (11) contains the important I Ca current as well as the fast potassium current I Kr , which are important for EADs to establish. Hence, system (11) may exhibits EADs dependent on the choice of system parameters.  Furthermore, we know that a reduced (fast) potassium current and/or an enhanced calcium current may leads to EADs, cf. [37]. In addition, in [52] it is shown how combinations of reduced and enhanced potassium and calcium currents increases the risk of the appearance of EADs, or conversely, may control the pattern of the APs. To illustrate the behaviour of system (11) Figure 13 shows a comparison of a normal AP and the occurrence of EAD patterns for different combinations of a reduced fast potassium and an enhanced calcium current. It is well known that EADs occur either in the plateau or in the repolarisation phase of the AP and they are benefited by an elongation of the AP. This may happen by an increase of one or more inward currents and/or a decrease in one or more outward currents [33]. In addition, it is well established that the calcium current I Ca plays an important rule during the plateau phase, while the potassium current I K during the repolarisation phase, cf. [52]. Figure 13(b) shows that EADs appear as a combination of a reduced fast potassium current I Kr and an enhanced calcium current I Ca . Therefore, we restrict our analysis to the case where we have a 80% block of the fast potassium current by introducing a new conductanceḠ Kr = 0.2 · G Kr and choosing G Ca as bifurcation parameter.
Similar to the analysis of the Noble model (1) we start by determining the equilibrium curve, cf. Figure 14(b) and (c). Here, we have again an unstable (black dashed line) and a stable (black solid line) equilibrium branch. The equilibrium curve changes stability via a subcritical Andronov-Hopf bifurcation (red dot, G Ca ≈ 0.096017 mS cm 2 ) with a positive first Lyapunov coefficient. From the subcritical Andronov-Hopf bifurcation an unstable limit cycle branch (dashed red line) bifurcates, cf. Figure 14 (11) including a 80% block of IK r , i.e.ḠK r = 0.2 · GK r , using GCa as bifurcation parameter.
The first limit cycle branch contains a limit point of cycle bifurcation (solid green square, G Ca ≈ 0.096259 mS cm 2 ) and a period doubling bifurcation (solid red square, G Ca ≈ 0.096255 mS cm 2 ) from which a stable period doubling cascade bifurcates. The first limit cycle branch changes stability via the limit point of cycle bifurcation and again via the (first) period doubling bifurcation. The bifurcation diagram in Figure 14 contains the first two limit cycle branches, where both limit cycles terminate at the unstable equilibrium branch. The limit cycle branches are mostly unstable and therefore, not attracting. Nevertheless, they influence the dynamics of the system, i.e they at least prolong the plateau phase and may cause EADs, provided the initial stimulus is strong enough to establish an AP (as it is in the standard setting).
In addition to Figure 14(b) and (c), we provide in Figure 14(d) the corresponding 3D bifurcation diagram including three different trajectories. From this it is obvious that the trajectories curl around the limit cycles resulting in a prolongation of the AP and/or EADs. Indeed, as soon as the trajectories enters the inside of the limit cycle branches, they will converge into the stable equilibrium branch. Notice that the stable equilibriums close to the Andronov-Hopf bifurcation are less attracting than others, since some of the negative eigenvalues are very small, but still negative. However, if the trajectory of system (11) is in the basin of attraction of the stable period doubling cascade (the stable attracting parts of the limit cycle branches), the system develops self-oscillating behaviour (Figure 15), or chaos ( Figure 16). A setting for this to happen is the combination of a G Ca value of the period doubling cascade with initial values in the basin of attraction and I stimulus = 0. Figure 15 contains four simulations of system (11) at the first four period doubling bifurcations of the period doubling cascade. The black fine line denotes the trajectory over 10000 ms, while the red line indicates the length of the trajectory with period T , i.e. V (T ) − V 0 = 0. Additionally, Figure 15 provides the corresponding phase space (x r , V ) to these simulations showing a closed curve starting from a red dot and terminating at a blue one. Note that the red dot is overlaid by the blue one due to V (T ) = V 0 and therefore, barely or not visible.
Finally, depending on the initial values the period doubling cascade is again a route to chaos similar to situation for the Noble model (1), see Figure 16. Again, it is clear that besides the system parameters also the initial values, and additionally the external stimulus, play a crucial role for the occurrence of certain dynamics and patterns such as (normal) AP, EADs or chaos. This indicates that a disorder in the external stimulus may also initiate a sudden death (at least on the cellular level).
Our analysis shows on the one hand that the dynamics of a single cell model are sensitive to its system parameters, and on the other hand they are sensitive to the choice of initial values. This is in accord with our previous analysis. However, the behaviours of the Noble model (1) and system (11) are quite different. Note that the strength of the initial stimulus influences the initial values and influences the dynamics of system (11).
Notably, the occurring EADs in Figure 14(a) appear in a physiological feasible range, cf. [17], acting as a validation of the model. However, the chaotic behaviour in Figure 16(a) is most likely non-physiological due to the small voltage range, and we would expect cardiac death as in Figure  16(b) to happen in the real cell.

Effects on the macro-scale
Based on the analysis in the previous section, we study the synchronisation behaviour of an ensemble of cells along a 1D cable of 1 cm to gain an intuition on whether EADs can spread or not in cardiac tissue. As for the monodomain model (6), we split the cable into two parts D and [0, 1] \ D. In the domain D, we set the calcium conductance to G Ca = 0.096229 mS cm 2 to ensure that an EAD with six additional oscillations (the pink line in Figure 14(a)) would occur in the ODE model (11), while in the remaining part it is set to G Ca = 0.09616 mS cm 2 (close to the occurence of EADs). We keep the diffusion small and set the diffusion constant to 0.00005 mS. Figure 17 shows the dynamics of an ensemble of 128 cells where 1%, 2%, and 50% of the cells are set to the six-oscillation setting for the ODE. We observe that no EADs occur when 1% of the cells are EAD prone (Figure 17(a)), while there is one additional small oscillation on parts of the cable for 2% EAD prone cells (Figure 17(b)). Increasing the percentage of EAD prone cells to 50%, we see that there are small additional oscillations along the whole cable ( Figure 17(c)). All three experiments show fewer small additional oscillations than in the ODE case.
In Figure 17(a)-(c) the cells surrounding the EAD prone cells (x ∈ [0, 1] \ D) are very close to establishing EAD behaviour. Hence, only a very small percentage of EAD prone cells are needed for EADs to occur along the cable. If the surrounding cells are further away from the EAD setting (G Ca = 0.064 mS cm 2 , i.e. the standard setting), we can observe that no EADs occur even if 50% of the cells are set to EAD inducing behaviour, see Figure 17(d).
Finally, we briefly study the effects of initially setting the cells in the region D = [0.2, 0.7) to be prone to chaotic behaviour. From Figure 18 we observe that chaotic behaviour does not spread in the three cases considered. In the first two simulations ( Figure 18)(a)-(b)) the surrounding cells are set to produce normal APs, which they indeed do for both the diffusion constant 1 360 mS and 0.00005 mS. However, if the surrounding cells are set to produce normal APs but are close to the EAD setting, we see that the dynamics die out (V → 0). In conclusion, chaotic behaviour does not spread even if 50% of the cells are initially chaotic. However, cardiac death can occur in the model (V → 0), see Figure 18(c). In Figure 18(a) the system produces normal APs (apart from the fact that the middle cells are initially set to zero). However, reducing the diffusion constant, the non-physiological behaviour that we saw in the previous section appears, see Figure 18(b). This indicates a critical lower bound on the diffusion constant in this model for the dynamics to be physiologically relevant.
Concerning EADs, the above results indicate that the spreading of EADs on the tissue level is depending on 1. the number of cells prone to establish EADs, 2. how close to the normal setting the surrounding cells are, and 3. the diffusivity of the monodomain model (12).

Summary and discussion
In this paper, we investigated and analysed the behaviour of two mathematical models describing the action potentials of a Purkinje and a human ventricular cardiac muscle cell. To this end, we utilised bifurcation theory, numerical bifurcation analysis, and computational tools to establish an increased understanding of the dynamics of these models. This enabled us to find hidden features in both models considered. Furthermore, carrying out this analysis, we aimed at convincing the reader that 1. bifurcation analysis is very beneficial in the study of the dynamics of an ODE model and to detect hidden features of the considered system, 2. it is important to know how to interpret the corresponding bifurcation diagram, and 3. advancing cardiac cell research benefits from collaborations between mathematicians and physiologists/biologists.  First, we studied the dynamics of the Noble model (1) with respect to the leak current I L based on the discussion in [15]. In [15] the author already varied the leak current conductance G L , resulting in the observation that the conductance G L influences the period of the AP. Even more, if one chooses G L large enough, e.g. G L = 0.4 mS cm 2 , the system converges into a stable equilibrium and no AP can appear, cf. [15]. This behaviour was analysed in more detail using numerical bifurcation theory. It turns out that this system changes stability via a supercritical Andronov-Hopf bifurcation from which a stable limit cycle branch bifurcates. This limit cycle branch loses and wins stability via limit point of cycle bifurcations and a period doubling bifurcation, respectively. Moreover, from the first period doubling bifurcations of the second limit cycle branch a (stable) period doubling cascade bifurcates, which is also the route to chaos. Interestingly, every limit cycle branch contains two period doubling bifurcations, which are connected via two limit cycle branches. Dependent on the initial values and the choice of G L , the Noble model (1) exhibits complex patterns and chaos. However, although mathematically interesting, the chaos detected (see Figure 6) might be an artefact of the mathematical model and not within the physiologically relevant range.
We used the same approach to study the more complex model (11), more complex in the sense that it contains more ion currents, pumps and exchangers compared to model (1), including the missing calcium current I Ca and the fast and slow potassium current, I Kr and I Ks . On the single cell level it turned out that system (11) exhibits both chaotic behaviour and EADs via a combination of a reduction of the fast potassium current and an enhanced calcium current.
For system (11), we showed that for a 80% block of the fast potassium current and an en-hanced calcium current, system (11) exhibits a subcritical Andronov-Hopf bifurcation from which an unstable limit cycle branch bifurcates, which stays unstable and contains a period doubling bifurcation. From this period doubling bifurcation a stable period doubling cascade bifurcates, which causes both EADs and deterministic chaos. Note that also for other (fast) potassium block rates system (11) contains a subcritical Andronov-Hopf bifurcation and may exhibit EADs and/or chaos. As experiments have shown that this is how EADs can occur, we can interpret these findings as a further validation of the model. However, it is unclear whether the occurring chaotic and self-oscillatory patterns also appear in a biological cardiac muscle cell (or are modelling artefacts).
In particular, the dynamics in Figure 16(a) is likely not appearing within a physiologically relevant range. This has to be investigated further with accurate experiments. Besides comprehending the dynamics of the single cell models, it is crucial to understand how these dynamics affect the behaviour on the macro-scale (cm) due to the fact that multiple cardiac single cells may synchronise and cause arrhythmias. To this end, we introduced monodomain models for both systems, cf. (6) and (12). Based on the analysis of the single cell dynamics, we investigated cell synchronisation in both models. Both analyses showed that 1. the diffusivity of the model, 2. the number of cells, and 3. the placement of chaotic/EAD regions affect the global dynamics of the monodomain models. Furthermore, the non-physiological behaviour observed for single cells can transfer over to the macro-scale models, but this depends on the size of the diffusion constant. In particular, this warrants a bifurcation analysis with respect to the diffusion parameter for the two models considered.

Conclusion
Bifurcation theory is itself a powerful tool to study the behaviour of dynamical systems. It is used in many contexts and gets more and more attention also in (mathematical and computational) cardiac and neuroscience. In cardiac science, one can use this approach to establish a better understanding of cardiac arrhythmia [11,12,13,14,31,35,52]. In interdisciplinary research, bifurcation theory can also be an important component in successful treatment of human diseases.
We considered two specific cardiac cell models among a multitude. Although a large number of models exist, there is still a lot of future work to do to derive a complete understanding of all cardiac dynamics. One has to deal with several issues, e.g. complexity of realistic models, and numerical and computational problems. In the very end, the extension from the cellular level to the tissue level has to be understood [42,54,55,56,57]. In particular, the physiological relevance of complex dynamics at all levels of modelling, as the ones we consider in this paper, has to be proven by experimental data.
In conclusion, this research would benefit from close interdisciplinary collaborations, since 1. a good, robust and realistic mathematical model can be developed based on experimental data, 2. an in-depth mathematical analysis can validate the accuracy or display weaknesses of the model, 3. these new findings either would help to improve the modelling or reformulation of the system to derive a most realistic model, including all expected dynamics, 4. the analysis would also increase the understanding of the occurrence of certain phenomena, and 5. the new obtained knowledge would help to develop new potential treatments for human diseases.