Bifurcations, chaos, and sensitivity to parameter variations in the Sato cardiac cell model

The dynamics of a detailed ionic cardiac cell model proposed by Sato et al. (2009) is investigated in terms of periodic and chaotic action potentials, bifurcation scenarios, and coexistence of attractors. Starting from the model’s standard parameter values bifurcation diagrams are computed to evaluate the model’s robustness with respect to (small) parameter changes. While for some parameters the dynamics turns out to be practically independent from their values, even minor changes of other parameters have a very strong impact and cause qualitative changes due to bifurcations or transitions to coexisting attractors. Implications of this lack of robustness are discussed. © 2016 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Mathematical modelling has become an important tool in the life sciences to address problems that are not approachable experimentally. For example, to investigate cardiac arrhythmias or sudden cardiac death, models have been developed to describe cardiac dynamics from gene to organ level [2,3]. Therefore, many models exist to describe the action potentials (AP) of single ventricular cells. For better comparison with experimental results, these models are often developed to represent dynamics of specific mammals: For example, the well-known Luo-Rudy model [4] can be used to model guinea pig ventricular cells, while the model used by Wang and Sobie [5] describes mouse ventricular action potentials, and the one used by Sato et al. [1] (first AP model, in the following referred to as the Sato model) models ventricular rabbit myocytes.
However, many of these models share the feature of consisting of a great number of equations and parameters. For example, the Sato model uses 27 variables whose calculation requires 118 equations and 177 parameters. Therefore, the exact implementation of these models as given in the original publications is quite error-prone, and since changes in parameters can dramatically alter the dynamics of the model, reproducing results from former studies can be a challenging task.
Studies analysing the parameter sensitivity of electrophysiological models are still rare [6,7]. Therefore, in this paper we use the Sato model to investigate the sensitivity of the dynamics to parameter variations. This cardiac cell model was used to provide an explanation for ventricular tachycardia and ventricular fibrillation originating from early afterdepolarisations at the cellular level by a chaos synchronization mechanism [1,8]. As will be shown in the following sections this model exhibits chaotic action potentials as well as coexistence of different types of periodic and chaotic attractors. To investigate the robustness of the dynamics given by the standard parameter set of this model [1], bifurcation diagrams are computed showing dynamical changes when going below or above the standard parameter values and it turns out that for some parameters already minor deviations from the standard value lead to qualitatively different dynamics.

Methods
The Sato model, which is given in detail in Appendix A, describes the action potentials in the membrane voltage V(t) of ventricular rabbit myocytes. It uses 27 variables that can be described by 118 equations involving 177 parameters. Besides ionic concentrations, channel conductances or physical constants, many of these parameters are coefficients, obtained from fitting mathematical models to experimental data. Since the sensitivity of the model dynamics to changes of these fitted parameters is to be examined as well, all fit coefficients were labelled (see Tables A.3-A. 6) and included in the total list of model parameters.
The cardiac cell model is given as a set of ordinary differential equations (ODE's) where C m is the membrane capacity and I stim the external stimulation current with pulses of 1 ms duration and amplitudes of −40 μA/cm 2 . I ion is the sum of all considered transmembrane or intracellular ionic currents. The state vector h includes all additionally needed, time dependent variables like gating variables or ionic concentrations. The ODE system (1) was solved in C++, using the Nordsieck BDF method with adaptive time steps implemented in the GNU Scientific Library [9]. The maximum absolute and relative error tolerances were set to 10 -8 with a maximum allowed time step of 0.5 s. Since this ODE solver requires the Jacobian matrix of the ODE system (1), the symbolic Jacobian was calculated using the GiNaC framework [10]. For the Sato model, no initial conditions were given in the original publications [1,11]. We obtained steady state initial conditions by setting all 27 variables to 0.1 prior to pacing the system with a constant PCL of 0.8 s until a steady state was reached, which took about 300 s of simulated cell activity. The initial conditions obtained by this procedure are given in Table A.2.
The model was originally designed to reproduce cardiac dynamics at fast pacing [11] and includes modifications to be capable of generating early afterdepolarisations (EADs) [1], which are known to be potential triggers of lethal cardiac arrhythmias [12,13]. A detailed description of EADs and the underlying ionic mechanisms can, for example, be found in Ref. [14]. Briefly, an EAD is an abnormal depolarisation during the plateau phase of an AP, which essentially prolongs the action potential duration (APD). The APD for each AP is defined here as the time duration where V > −85 mV. To analyse the action potentials and the appearance of EADs, we also introduced a Poincaré section to the 27-dimensional state space such that every action potential is represented by the values of the 27 variables at a fixed phase after the initial depolarisation of the AP. These chosen state space points need to reflect the dynamics during the plateau phase of the action potential. Since the shape and duration of the action potential are PCL-dependent, we defined the phase space point v p n of the nth action potential of variable v as v p n := v(t 0 + (n + 0.3)PCL) (2) where t 0 is the time of the depolarisation of the first AP considered. Thus, the phase increases linearly with PCL. Furthermore, for calculating Lyapunov exponents, a discrete QR decomposition based method as described in [15] was implemented.   Fig. 2 shows a three dimensional projection of the chaotic attractor underlying the time series shown in the middle panel of Fig. 1. Colours indicate the average concentration c j of free Ca 2+ in the sarcoplasmic reticulum (SR).

Periodic and chaotic action potentials
To investigate the occurrence of EADs for a larger PCL range, we calculated the APD for 200 APs for 1.1s < PCLs < 1.4 s with a step size of PCL = 0.5 ms. Prior to the APD calculation, we paced the system for 500 s of cell activity and for the nth PCL n , we used the final state vector of the previous PCL n−1 = PCL n − PCL as initial condition. The resulting bifurcation diagram in Fig. 3A shows that for a certain PCL range, the APD takes many different values, which shows the irregular behaviour in the appearance of EADs. Note that an APD value larger than about 0.5 s represents an EAD. For smaller and larger PCL values the APD does not vary, which indicates periodic behaviour. Furthermore, the irregular behaviour in the intermediate PCL range is interrupted by periodic windows. Fig. 3B shows a bifurcation diagram where instead of APDs   Fig. 1). Plotted is a three dimensional projection into a subspace of the state space spanned by the membrane voltage V, the intracellular Na + concentration [Na + ] i , and the average concentration c i of free Ca 2+ in the cytosol. The colour encodes the average concentration c j of free Ca 2+ in the sarcoplasmic reticulum (SR).
values of the membrane voltage V p in the Poincaré section are plotted which are calculated according to Eq. (2). Note that a V p value larger than about −30 mV represents an EAD. While the shape of the two bifurcation diagrams shown in Fig. 3A and B differ due to the different quantities used, the same irregular behaviour in the intermediate PCL range is visible. Furthermore, in the bifurcation diagram using the membrane voltages V p in the Poincaré section the periodic windows are more clearly visible: a period-2, 3, 4 and 5 window can easily be identified, which is not the case in the APD bifurcation diagram. Since the latter contains no further information that the V p bifurcation diagram lacks, we will use the values V p in the Poincaré section for further analysis. We found that for 0.3 s < PCL < 4.0 s, there is no other PCL range exhibiting irregular V p behaviour except for the one described above.
To test whether the irregular appearance of EADs is due to dynamical chaos, the three largest Lyapunov exponents λ were calculated and are shown for 1.1 s < PCLs < 1.4 s in Fig. 3C: The PCL values where APD and V p show irregular behaviour correspond to a positive largest Lyapunov exponent λ max , which shows that the irregularity is dynamical chaos.
All of the results above are similar to the results by Sato et al. [1] but the chaotic PCL ranges differ quantitatively. To the best of knowledge, the model was implemented exactly as described in [1,11]. Due to the number of equations and respectively, the APD is constant. The irregular behaviour is discontinued by periodic windows, a typical characteristic for chaotic models. B: Same as A, but with the dynamics described by phase space points V p . C: The three largest Lyapunov exponents λ over PCL: When APD and V p fluctuate irregularly, the largest Lyapunov exponent becomes positive, which shows that the irregularity is dynamical chaos. parameters, the possibility of errors while transferring the equations from paper to computer or vice versa cannot fully be excluded. In the following, we therefore analyse the sensitivity of the model dynamics to parameter changes.

Bifurcations
A simple, but insightful approach to qualitatively analyse the parameter sensitivity of the model is to vary every model parameter independently and to examine and compare the resulting bifurcation diagrams. For our purposes, a variation for each parameter P by up to 100 % around its standard valueP as given in Tables A.3-A.6 is considered to be sufficient. More precisely, for every P, we calculated the membrane voltages V p for 200 APs for (P/100) < P < 2 ·P with a step size of P = (P/100). Prior to the V p calculation, we paced the system for at least 300 APs and for the nth parameter value P n , we used the final state vector of the previous P n−1 = P n − P as initial condition. For P 0 we used the initial conditions as given in Table A.2. This "upward" calculation for every P (i.e. from small to large parameter values) was followed by an analogous "downward" calculation (i.e. from large to small parameter values) where again we used the final state vector of the nth parameter value calculation P n as initial condition for the following parameter value P n+1 . This calculation was carried out for PCL = 1.282 s such that the model should yield chaotic dynamics if P =P according to Figs. 2 and 3C. In Fig. 4, a selection of the bifurcation diagrams is shown. The standard valueP for each parameter is represented by a vertical line, and the differently coloured points representing the membrane voltage V p in the Poincaré section reflect the previously explained "upward" (light blue) and "downward" (dark red) calculation. Note that the "downward" calculation is plotted on top of the "upward" calculation, such that the dark red points can essentially hide the light blue points, meaning that the dynamics in the affected parts of the bifurcation diagram are equal for both calculations. The bifurcation diagrams show that a wide range of different behaviour can be observed under individual variation of different parameters: • Fig. 4A shows that the chaotic dynamics, which is expected for PCL = 1.282 s according to Figs. 2 and 3, is stable under parameter variation of parameter eq48P1, independent of whether the variation occurs from small to large values or vice versa.  period-5 orbit (with two nearby points at V p ≈ −3 mV) for a critical value of eq92P1 ≈ 33. If now decreased again, the period-5 orbit does not become unstable at eq92P1 ≈ 33 but remains stable up to eq92P1 ≈ 12.4: In this parameter range, at least two attractors coexist with periodic and chaotic behaviour, respectively. Note that this range also includes the standard value of eq92P1 = 20. At eq92P1 ≈ 12.4, the period-5 orbit becomes unstable and the dynamics of "upward" and "downward" calculation equal.
• Fig. 4C: The variation of eq73P2 shows that the "upward" variation of this parameter yields chaotic dynamics for almost all parameter values (with two exceptions for eq73P2 ≈ 3 and eq73P2 ≈ 13), while the "downward" variation yields periodic dynamics with period 5 (with two nearby points at V p ≈ −3 mV) except for very large values of eq73P2: Therefore, once the trajectory is bound to one of these attractors, both are relatively stable to parameter variations of eq73P2.
• Fig. 4D: Except for very small values of the dissociation constant K mem , the dynamics is bound for both "upward" and "downward" calculation to the already mentioned periodic attractor with period 5: Once bound to this attractor, the dynamics is very stable to variations in K mem and chaotic motion is never observed in this case. For decreasing values of K mem , the period-5 orbit becomes unstable and changes over to a period-3 orbit, while the "upward" calculation exhibited different behaviour in this regime.
• Fig. 4E shows that for increasing values of eq89P1 the dynamics follows periodic orbits of different periods and subsequently shows chaos for a fairly narrow parameter range around the standard value of eq89P1 = 3. This chaotic behaviour becomes unstable if eq89P1 is further increased and eventually leads to a period-2 orbit. Note that the chaotic behaviour does not occur for parameter values smaller than the standard value, meaning that the motion is highly sensitive to slightly decreased values of eq89P1. In this bifurcation diagram both "upward" and "downward" calculation yield the same dynamics for all parameter values larger than the standard value. However, for smaller parameter values it is clearly visible how the dynamics depends on the "direction" of parameter variation: For decreasing parameter values, the motion is bound to the respective periodic attractor longer than for the "upward" calculation before it becomes unstable and goes over to a different periodic attractor. • Fig. 4F: The bifurcation diagram for K mNao shows that if this parameter is only slightly varied to either side around its standard value of K mNao = 87.5, the dynamics becomes periodic. Several more chaotic regions exist for larger and smaller parameter values, respectively. Also, in the "downward" calculation, the motion is periodic even for the standard value. Therefore, the dynamics is highly sensitive to changes in K mNao . • Fig. 4G: For small values of the threshold for leak onset k j , the motion exhibits chaotic behaviour which is interrupted by periodic windows for increasing values of k j . Period doubling bifurcations are clearly visible within these periodic windows. Subsequently, chaotic motion also appears around the standard value of k j = 50, but quickly becomes periodic on either side of this value. For even larger values of k j , the motion goes over to the already mentioned period-5 attractor which is stable over a fairly broad parameter range. If decreased again, the dynamical behaviour is nearly identical with the exception that the period-5 attractor remains stable for values of k j even smaller than the standard value. These examples illustrate that the model parameters can strongly influence the dynamics of the system and even small deviations from the given standard parameter values can lead to very different dynamical behaviour. In addition, for the standard parameter values a periodic attractor with period-5 coexists to the chaotic attractor which occurred for PCL = 1.282 s in Fig. 3. The possibility of more coexisting attractors cannot be excluded. As can be seen from the bifurcation diagrams, some parameters influence the dynamics of the system while others do not. Also, the sensitivity of the dynamics to changes in the parameters varies from parameter to parameter. We therefore heuristically defined five parameter classes in order to group parameters with similar influence on the dynamics: The first class contains parameters which do not alter the dynamics of the system if varied by up to 100%, as eq48P1 in Fig. 4A. The second and third class contain parameters which qualitatively influence the dynamics if varied by more (as eq92P1 in Fig. 4B) or less (as eq89P1 in Fig. 4E) than 10 %, respectively. All parameters of Table A.7 are excluded from this grouping as they consist of physical constants or other physical parameters like the cell volume or the temperature. These parameters form the fourth class. The fifth class contains all those parameters where numerical difficulties (e.g., divergence of solutions) occurred for parameter values P n far fromP and the "upward" or "downward" calculation of the bifurcation diagrams failed at some point. 1 The parameter classes, together with a short description and the number of parameters per class are shown in Table 1. Here, three PCLs (PCL = 1.1 s, PCL = 1.282 s, and PCL = 1.37 s) are considered and the corresponding parameter bifurcation diagrams (similar to Fig. 4) are used to classify the model parameters. Furthermore, we grouped the parameters into different categories: Ca 2+ , Na + , K + , and the rest. This categorisation was motivated by the hypothesis that the impact of each parameter may depend on the ion dynamics it is involved in. Table 1 Parameter classes and number of parameters per class. For each PCL, the number of parameters in class 0 is the number of parameters needed to give a grand total of 177 parameters (including the constants).

Class
Definition 0: not grouped Parameters that caused numerical difficulties and were not grouped with this procedure 1: non sensitive Parameters whose variation up to 100 % does not effect the qualitative dynamics (periodic, chaotic, …) of the system 2: mildly sensitive Parameters whose variation of more than 10 % and less than 100 % does effect the dynamics of the system qualitatively 3: highly sensitive Parameters whose variation of less or equal than 10 % does effect the dynamics of the system qualitatively 4: constants Physical constants and ionic concentrations (11 parameters  For PCL = 1.1 s 72 parameters belong to classes 2 or 3 and therefore qualitatively alter the dynamics if varied by less than 100 %. From these, 14 are "highly sensitive" parameters, i.e. a variation of less than 10 % of their values changes the dynamics qualitatively. For PCLs 1.282 s, resulting in chaotic dynamics, and 1.37 s the number of sensitive parameters is even larger (137 and 122, respectively). To address the question whether particular parameters exhibit for all three PCLs (high) sensitivity and a major impact on the dynamics of the model Fig. 5 shows the class membership of all parameters and all three PCLs (differently coloured symbols). There are only 23 parameters in classes 1, 2, or 3 which belong for all three PCLs to the same group (see parameters marked with an asterisk in Tables A.3-A.6). In most cases, the class membership changes with the PCL and physiological factors determining the class membership (for arbitrary PCLs) are very difficult to identify. Even when grouping the parameters according to their role in modelling the dynamics of different ions (see Fig. 5) no direct relation between ion type and sensitivity class is visible. Therefore, we conclude that the observed bifurcation scenarios and the impact of parameter variations are not governed by biophysical or physiological effects but mainly due to the very large number of parameters (and model terms). From Statistical Learning Theory [16] it is known that overly detailed models ("overfitting") may suffer from adverse features like low generalisation ability and poor prediction properties. This is a general phenomenon that does not rely on the particular (physical) context of the model and appears to be relevant for cardiac cell modelling, too.

Conclusion
Many ionic cardiac cell models aim at detailed "realistic" modelling taking into account any known biophysical detail. As a consequence these models consist of a large number of variables and an even larger number of parameters. A large number of variables implies a large dimension of the state space. This is a priori not a problem but perhaps not necessary, because even the chaotic attractors of such cell models are usually quite low dimensional and therefore, a four or five dimensional state space, for example, would be sufficient to "host" the relevant attractors. Much more problematic is the very large number of parameters. The values of some parameters can be precisely measured but others are difficult to determine or are just imported from previous research (under partly different conditions or with different species). In any case, there may be limited information about the true and proper values of parameters and this raises the question of what happens if these parameters vary within a small interval (enclosing the conjectured but unknown "true" value)? The bifurcation diagrams obtained for the parameters of the Sato model clearly demonstrate that this is a severe issue. While some parameters have practically no influence on the dynamics others have to be known very precisely to achieve the expected or desired dynamics. This (extreme) sensitivity reduces the robustness of the model and its ability to generalise. These findings are in good agreement with results from statistical learning theory and machine learning where overfitting leads to strong degradation of regression and classification results [16]. Our attempt to interpret the class membership of each parameter in a physiologically meaningful sense failed, because most parameters occur in different classes depending on the PCL applied. There is also no evidence that parameters determining a specific ion dynamics are more (or less) sensitive than others. We conjecture that this feature of low robustness is a general "overfitting effect" and not due to some physiological or biophysical processes. In this sense we expect that it is shared by many (if not all) high dimensional ionic cell models and has to be taken into account when using these models for simulating or even predicting dynamical events like the onset of arrhythmias.

Appendix A. The Sato model
The following section contains the equations, parameters and initial conditions of the cardiac myocyte model investigated in this article. It is based on a model developed and described by Mahajan et al. [11], which is also available online as a cellML implementation [17]. As these versions differ slightly, the equations from [17] were adopted. Furthermore, the model is modified as described by Sato et al. [1] (first AP model), and further ambiguities where eliminated after correspondance with the authors. Thus, some of the following model equations contain terms that are not given or differ from the ones in the original publications. In the following, these terms are boxed . Furthermore, for parameter eq6P2 a value of 0.5 was used instead of 0.05 [1], to avoid very large derivatives dP r /dV when computing and applying the Jacobian matrix of the model.

Nonlinear buffering
The fast sodium current (I Na ) For V ≥ −40 mV: Averaged Ca 2+ dynamics in the dyadic space Inward rectifier K + current (I K1 ) The rapid component of the delayed rectifier K + current (I Kr ) The slow component of the delayed rectifier K + current (I Ks ) The Na + /K + pump current (I NaK ) The fast component of the rapid inward K + current (I to, f ) The slow component of the rapid outward K + current (I to, s )