DAD characterization in electromechanical cardiac models

We investigate the possibility of modeling the delayed afterdepolarization (DAD) occurrence in the framework of the classical FitzHugh-Nagumo (FN) dynamical system, as well as in more recent electromechanically-coupled cardiac models. Within the FN model, we identify the domain in the constitutive parameters' space for which orbits exist which exhibit a sufficiently strong secondary impulse. We then address the question whether a locally-induced secondary pulse succeeds or not in originating a self-propagating traveling impulse. Our results evidence that, in the range where secondary impulses exceed the physiological threshold for DAD onset, a local impulse almost certainly causes a traveling impulse (mechanism known as all-or-none). We then consider a recently proposed electromechanically-coupled generalization of the FN model, and show that the mechanical coupling stabilizes the system, in the sense that the more strong the coupling, the less likely is DAD to occur.


Introduction.
A delayed afterdepolarization (DAD) is a depolarization of the cardiac myocytes, which occurs after the completion of the final rapid repolarization phase of the cardiac action potential, but before the membrane potential returns to its equilibrium value. Depending on their intensity, DAD events are labeled as sub-or supra-threshold. The former occur when the depolarization does not exceed the threshold fit to generate a widespread new action potential, while the latter arise when the membrane depolarization overcomes such threshold. Consequently, a suprathreshold DAD triggers a spontaneous action potential, and eventually induces an aftercontraction in the heart. From the biochemical point of view, DAD occurrence is tightly linked to alterations in Ca 2+ load in the cardiac cells. In particular, the ionic mechanism which induces DAD is the fluctuation in release of calcium from the sarcoplasmic reticulum under calcium overload conditions -a phenomenon called calcium-induced calcium release [1,2,3].
It is the main purpose of the present study to analyze -in the framework of the classical FitzHugh-Nagumo dynamical system and recently proposed generalizations of the same -the values of the constitutive parameters for which orbits exist that exhibit a secondary impulse (or spike) sufficiently strong to model the generation of a DAD. It is to be remarked that, if and when a new contraction is triggered, the evolution of the action potential undergoes such an abrupt change that it becomes meaningless to follow the orbit details after that event. Therefore, we will pay attention to the choice of the threshold for extra-pulses to be induced, while counting the exact number of pulses exhibited by a particular orbit is less important in order to address the possible onset of DAD-triggered action potentials. In what follows, we will therefore search for solutions which possess at least a secondary spike.
The FitzHugh-Nagumo (FN) dynamical system [4,5] models the evolution of the fast action potential v, and the slow gating variable w The constitutive parameters α, γ, ε are positive constants: α ∈ (0, 1) represents a threshold for the amplitude of the excitations, the small parameter ε represents the ratio between the characteristic time scales of the fast and the slow evolution variables, while γ is typically adjusted as a fitting parameter [6]. Finally, I appl models the total membrane current density. In [7] we put forward the bases for the study of the spike solutions which emerge from the dynamical system (1). By definition, a spike in an orbit is counted each time the action potential v overcomes a fixed threshold v thr . The value of the threshold is to be fixed in order to guarantee that the corresponding pulse has the capability to generate a propagating signal in a finite region. In order to allow for a direct comparison with the results previously derived for a piecewise-linearized FitzHugh-Nagumo model [8], the threshold value in [7] was chosen to be coincident with the minimizing value of the cubic f α . Such a choice, however, provided a v thr explicitly depending on α, and in particular vanishing in the α → 0 + limit. It emerged that, if the threshold is allowed to become arbitrarily small when the constitutive parameters are varied, more and more trajectories may trivially be found with many spikes. It was then suggested, and it will be here implemented, that the threshold for spike counting should be indeed fixed, independently of the value of the constitutive parameters. Indeed, the physiologically significant threshold value for the action potential must be set equal to the potential value fit to induce a mechanical contraction in the muscle fibers.
By making resort to the database of experimental data it is possible to ascertain that the potential value necessary to trigger an action potential, once the Ca 2+ current is suitably regulated, ranges over the interval (−65, −60) mV. To provide specific examples, in [9] the authors set this value to −65 mV, in [10] to −62.9 mV, while in [11] the threshold depolarization required to trigger a spontaneous action potential is computed in terms of [Ca 2+ ] in the sarcoplasmic reticulum. From these results we can conclude that the potential value the DAD has to exceed in order to generate an aftercontraction is about 1 8 of the total amplitude of an action potential. More precisely, as the characteristic amplitude of the action potential associated with the dynamical system (1) with I appl = 0 is of order unity [12,13], we fix the threshold value for the appearance of spontaneous action potential to be equal to v thr = 0.125. Clearly, most of the following analysis could be easily adapted to different threshold or applied current values, with no dramatic consequences on the qualitative results, as long as a single stable equilibrium configuration exists [14].
In what follows, we will search for the critical values of the constitutive parameters associated with orbits which originate extra-spikes. More precisely, in Section 2 we study the original FN model, and in Section 3 we move towards a more precise simulation of the cardiac dynamics. A concluding section summarizes and discusses our results.
2. Extra-spikes in FitzHugh-Nagumo orbits. Let us consider the dynamical system (1). We fix a threshold in the phase plane, v = v thr . By varying the constitutive parameters ε, α, γ, and by applying similar techniques as the ones described in [7] (see in particular Section 3 therein), we look for the region in the parameter space where orbits exist which exhibit at least an extra-spike. To be precise, for any given set of values of the constitutive parameters, we integrate backwards in time the trajectory passing through the point (v, w) = (v thr , f α (v thr )). The proof in [7] evidences that orbits with more than one spike exist if and only if such a trajectory crosses at least once the half-line {(v thr , w), w < f α (v thr )}.
In Figure 1(a) several curves v thr,cr (α) have been plotted, each associated with a different value of γ. Their meaning is the following. For each γ, the reported curve separates the parameter space (α, v thr ) in two regions such that when v thr < v thr,cr (α) it exist choices of ε and the initial condition such that the resulting orbit exhibits at least an extra-spike, while only single-spike trajectories can be found for any ε if v thr > v thr,cr (α). In words, if we focus on a particular value of v thr (e.g., the physiologically relevant value v thr = 1 8 ), it is important to ascertain whether the displayed curves lie above or below the chosen value. Whenever the curves stay below the relevant threshold, we can be certain that no secondary spikes are possible for any value of and the initial point. To this aim, it is interesting to report that, in correspondence to the physiological values (α, γ) phys = (0.1, 0.5), the critical value v thr,cr (0.1) is close to the value 1 8 evidenced above. More precisely v thr,cr (0.1) is just smaller than 1 8 , thus signaling on the one hand that no orbits should exhibit suprathreshold spikes if (α, γ) attain their typical values, but extra spikes may arise as soon as either α or γ slightly decrease their values.   curve, while no extra spikes can be found for the values placed at the right-side of the curve. Therefore, to make an example, if v thr is greater than right-most value attained by a specific curve, no extra-spikes can be found for any . The shape of the curve evidences that must exceed a minimum value for extra-spikes to arise, which is in agreement with our physical intuition that very small values of induce a well-known dynamical process in which the state variables evolve very close to the nullcline w = f α (v), and no room is left for finite-amplitude oscillations around the equilibrium configuration. It is, in our opinion, less significative the existence of an upper bound in for extraspikes to exist, as such upper bound is certainly much greater than the (small) physiological values of .
In summary, the reported analysis suggests that the FN dynamical system predicts no multi-spike solutions, as long as the constitutive parameters maintain their physiological values, but that only small perturbations of the parameters themselves are required for extra-spikes to appear.

Traveling pulses.
If a DAD-triggered action potential generates a potential wave which spreads through the heart conduction system, the cardiac muscle may undergo a contraction (commonly called extrasystole) that interferes with the regular cardiac cycle and may result in a dangerous alteration of the heart beats. From this point of view a complete DAD-predictive model must provide information on the physiological conditions which allow for the propagation of the spike-induced electric pulse and eventually induce the afore mentioned arrhythmic phenomena. Such a request implies that we have to enrichen the model with a diffusion term accounting for the spreading of the potential waves through the tissue.
We will now show that -at least in a 1D framework -the FN dynamical system succeeds in modeling the occurrence of suprathreshold perturbations of the electric signal, and the possible generation of extrasystole as well. More precisely, we now focus on the question whether an extra-spike occurring at a precise position and time is able to originate a traveling pulse, which would model the suprathreshold DAD.
We develop this issue by considering the following one-dimensional extension of the dynamical system (1) Without loss of generality, the diffusion coefficient D in (2) may be set equal to 1 by suitably re-defining the unit length. In (2) we have adopted the single-diffusivity approximation, by taking into account only the diffusion of the potential into a cardiac fiber, under the assumption that the recovery process is not affected by the spatial distribution of the cardiac cells. A further, more restrictive assumption is to accept that the results derived in the previous section, i.e. the critical values of the constitutive parameters that specify the occurrence of secondary spikes, hold for the spatially-dependent model. This approximation neglects the fact that the propagation of the membrane potential in a spatial domain may induce dissipative mechanisms, due to the electrical interaction between adjacent myocytes and the transmission of the depolarization current [10]. Such phenomena may in turn affect the amplitude of the potential upstroke, and consequently the size of the oscillations occurring after the main pulse. As a result, the spatial coupling is expected to reduce the amplitude of the secondary pulse, and therefore to restrict the parameter domain associated with the emergence of multiple-spikes solutions. We must consequently interpret the following results as a series of necessary conditions for a suprathreshold DAD to emerge. We now adapt to our present setting the analysis put forward in [15], which can be retrieved from our results by simply setting ε = 1 (see [16] for a detailed derivation of the present calculations). To search for traveling solutions we introduce the moving coordinate ξ = x+ct, where c denotes the wave speed to be determined, and consequently transform equations (2) in a system of two ordinary differential equations by introducing v(x, t) = ϕ(x + ct), and w(x, t) = ψ(x + ct). According to the definition of a traveling pulse solution [15,17,6], we further require that and that ϕ, ψ are smooth enough. Appropriate manipulations of the transformed equations show then that traveling solutions for the dynamical system (2) exist only if the wave speed is bounded from above and below by two parameter-dependent limits. Such bounds represent further necessary requirements on the constitutive parameters that are to be satisfied for traveling pulse to exist. The lower bound is achieved by integrating from ξ = −∞ to ξ = +∞ the two ODEs which describe the traveling-solutions. Algebraic manipulations and integration by parts provide the two equations [15,16] If γ < γ cr = 4/(1 − α) 2 (this condition implying that system (2) has one stable equilibrium point [7]), the left hand side in (4) must be positive, since the integrand can be written as Then c 2 > γ/ε must hold, which, substituted in (5), implies that the left-hand side of this latter equation must be negative. Therefore Condition (7) is a second-grade polynomial inequality in ϕ(ξ), and admits solutions if and only if c 2 > 4/(ε(1 − α) 2 ) > γ/ε, which provides the lower bound for c The upper bound for c is deduced by deriving from the dynamical system the following differential inequality for the function H(ξ) = e −cε/(2ξ) ξ −∞ ϕ 2 (η) dη [15,16]: As H(ξ) ≥ 0 by definition, were ε 2 c 2 ≥ 2(1 − α) 2 , H would be convex. However, it follows from the definition of H that it can be convex if and only if it vanishes, thus implying that if c 2 ε 2 ≥ 2(1 − α) 2 the only existing traveling wave solution for (2) is a constant. We have thus proved by contradiction that the existence of a non-constant traveling pulse provides a strict upper bound for the pulse speed Moreover, by comparing relations (8) and (10) we come to the final expression from which we deduce that a critical value ε cr,w (α) exists such that pulses are prevented from propagating through the fibers if ε ≥ ε cr,w . This critical value needs to be compared with the ε-values which allow for the occurrence of extraspikes, presented in Figure 1(b). If, to make an example, we choose α = 0.15, from (11) we can compute ε cr,w (0.15) . = 0.261. The value of v thr (evidenced in Figure  2) corresponding to ε = 0.261 on the curve associated with α = 0.15 marks the distinction between the region where a DAD can trigger an extrasystole (grey area in Figure 2) from the region where the DAD does not propagate. More precisely, if the time scale parameter ε exceeds the critical value ε cr,w two circumstances may arise: either an after-depolarization appears but the triggered action potential is not transmitted to the adjacent cells, or the potential does not oscillate after the main impulse. The physical interpretation of this phenomenon is quite clear. If, in a portion of tissue, the dynamical time-scale of the transmembrane potential is too similar to that of the recovery process, the possible depolarization stimulus generated in this region is damped down by the nearest cells since their characteristic response time is too slow.   details, we stress that we are combining necessary conditions for the existence of suprathreshold DADs. Therefore, the displayed results are to be read in the sense that the failure of any of the above conditions suffices for preventing the emergence of aftercontractions. In other words, by analyzing Figure 3 we are able to identify a parameter region where suprathreshold DADs (and therefore aftercontractions) are certainly prevented, as such events are requested both to overcome a physiological threshold and to guarantee that the signal is able to propagate. Clearly, this is not to be read as a complete characterization, as even in the critical region (where both multi-spike orbits are allowed and pulses may travel) we can say that the pathological behavior may potentially emerge, not that it will certainly emerge. By taking into account the necessary condition for wave propagation (11), we define a combined threshold v thr,cr,w (α, γ) as the supremum of the threshold values v thr such that there exist orbits with the specified values of (α, γ), and a value of ε allowed by (11), which exhibit a secondary spike stronger than v thr . Figure 3 evidences that v thr,cr,w is a decreasing function of α, which vanishes for sufficiently large values of α. Its physiological meaning is that DAD events are prevented, in the framework of the spatially-dependent FitzHugh-Nagumo model, whenever the physiological threshold for contractions to arise happens to be greater than v thr,cr,w . Clearly v thr,cr,w ≤ v thr,cr , as the latter represents the critical threshold for secondary spikes to appear locally, while the former embeds also the possibility for a perturbation to propagate in space. Figure 3 shows further that it exists a range of values of α -close to, but not including the physiological value α = 0.1 -such that v thr,cr,w > v thr,phy , and therefore suprathreshold DADs may emerge. It is to be noticed that in this range the two curves v thr,cr,w and v thr,cr are very close to one another, thus indicating that if the depolarization reaches the physiological activation threshold, it is very likely that a traveling pulse will follow, and then the heart will undergo an extrasystole. This effect is in agreement with the so called all or none response typical for both the neuronic and the cardiac electrical activation [18]. Indeed the upstroke initiation is not related to the amplitude of the stimulus, but rather to the fact that the activation threshold is reached or not. Therefore, if the potential depolarization approaches v thr,phy , the myocyte is excited and such a stimulation easily activates the surrounding tissue. Results in Figure 3 have been obtained by setting the fitting parameter γ to its physiologically relevant value γ = 1 2 . Nevertheless, all the qualitative features stand for a wide range of (strictly positive) values of γ.
3. Electromechanical cardiac models. The action potential evolution occurs in a mechanical environment, the cardiac muscle, which strongly interacts with the electric signal itself. Much effort is therefore presently spent on studying generalized FN models which explicitly include the effects of the interaction between the electric potential and the elastic properties of the surrounding muscle. It is in particular essential to include in the model the mechanisms through which the electric potential induces the contraction of the heart tissue. To this aim we here follow the recent proposal put forward in [19,20] where, along with the electro-mechanical coupling, the authors take into account the spatial coupling of the excited cardiac cells as well as the diffusion of the action potential through the cellular interconnections (gap junctions). As a result, the propagation velocity of the action potential is influenced not only by the transmembrane currents but also by the electric flow spreading into the heart tissue. More precisely, we here follow [20] with the explicit aim of tracing the effects left by the electro-mechanical interaction at least in a simplified setting. The electro-mechanical generalization of the FN model [20] is represented by the following reaction-diffusion system where D represents the diffusion tensor, and f α (v) = v(1−v)(v −α). Quantitatively different, but qualitatively similar proposals can be found in [19,21,22,23], but we here refer to [20] for a more complete discussion about the modeling details. Both the divergence and the gradient in (12) are defined with respect to the spatial coordinates, which characterize the present configuration of the system. We now perform a pull-back operation of all scalar and vector fields, to rephrase equation (12) in a material frame of reference [20,24]. As a result we obtain where the Div (and Grad thereinafter) operator denotes now the divergence (gradient) with respect to the material coordinates, and J = det F, F being the deformation gradient. We now assume that D is isotropic, i.e. D = DI, with D a uniform diffusivity and I the identity tensor. If we further recall the Piola identity Div(JF − ) = 0 (which applies since the deformation acts as a vector potential for The mechanical deformations within the cardiac muscle are mainly due to the elongation/contraction of the cardiac fibers. Consequently, we assume that a fiber distribution n(X) is given in the reference frame, and that the deformation gradient is a simple uniaxial extension in the fiber's direction where g is a function of the action potential, and ⊗ denotes the tensor product. Moreover, as no deformation is expected when the action potential is at rest (v = 0), we perform a linear approximation of the function g by assuming that g(v) = βv.
As the maximum fiber contraction is of the order of 30% at the action potential peak (v 1), typical values of the parameter β are close to β = 0.3. A further simplification, already introduced in [20], arises if we assume that all fibers are parallel in the reference configuration. When this is the case, and defining the Xcoordinate along the fibers' direction, the system (14) reduces to an evolutionary problem in 1+1 dimensions, and the only relevant component of the deformation gradient turns out to be and system (14) reduces to Equation (17) shows that the electro-mechanical interaction (β = 0) bears significant consequences, even if the diffusion process is neglected. More precisely, and in order to quantify the comparison with the results obtained in the spatiallyindependent framework of equations (1), we now neglect the spreading of the action potential through the gap junctions (that is we let D → 0) and we obtain the final system of equations The phase plane analysis of the dynamical system (18) is certainly more involved than that of (1). First of all, it is to be noted that it exists a singular value of the potential which induces critical behavior in (18): v = 1 2 β −1 . This unphysical behaviour is to be related with the linear dependence assumed for the functional relation F ax (v) in (16). Indeed, equation (16) is intended to represent the leading order of an expansion expected to be valid close to the equilibrium value v = 0. A different (and, in particular, suitably bounded from below) functional dependence F ax (v) would avoid this problem. It is however to be noted that the singular value 1 2 β −1 is certainly greater than the typical maximum value of v along an action 2660 PAOLO BISCARI AND CHIARA LELLI potential so that restricting our analysis to the region {v ≤ 1 2 β −1 } does not involve any lack of physiologically significant information.
The v-nullclines in (18) are Equation (19) is quadratic in w, and therefore may identify up to two reals solutions w 1 = f 1 (v), and w 2 = f 2 (v). It is a long but straightforward calculation to show that the real or complex character of w 1 and w 2 depends on the sign of the following sixth order polynomial in v More precisely, negative values of P (v) are associated with complex nullclines f i (v). The polynomial P is certainly positive when the action potential is sufficiently large in absolute value. At finite potentials, P (v) may or may not exhibit changes in sign. If we, for example, consider the physiologically significant case γ = 0.5, β = 0.3, ε = 0.01 (see Figure 4), we notice that for any α ∈ (0, 1) it exists an interval of (negative) values of the action potentials for which P (v) becomes negative. This feature signals the presence of values of v for which the trajectories are never horizontal in the (v, w) phase plane. Figure 5 illustrates a phase-plane representation of the dynamical system (18), exhibiting the presence of a range of values of v in which the w-nullcines are not defined.  The dynamical system (18) may admit up to five different equilibrium points: P eq,1 = (0, 0), P eq,2-3 deriving from the intersection points between the cubic vnullcline and the w-nullclines, and P eq,4-5 resulting form the crossing between the w-nullclines and the critical line v = β −1 . The latter equilibria are not of interest in our study, as we have agreed above in restricting our area of analysis to the region {v < 1 2 β −1 }. Moreover, in the physiological range γ < γ cr = 4/(1 − α) 2 (see [7]), P eq,2-3 od not have real coordinates, while the origin P eq,1 is a spiral sink. These considerations altogether evidence that, if γ < γ cr , the trajectories associated with the dynamical system (18) exhibit qualitatively similar properties as the well-studied converging trajectories of the classical FN dynamical system. After an external impulse that moves the potential away from the equilibrium value, the solution proceeds on the right branch of the v-nullcline compatibly with the vector field, then it moves on the left region, follows the stable branch of the fast nullcline and finally achieves the resting state. Since the origin is a spiral sink, orbits experience many oscillations before reaching the equilibrium, so that we are again allowed to look for multiple-spikes solutions of system (18). Figure 6(a) reports the values of the threshold potential v thr which allow for the onset of spikes solutions in the trajectories associated with (18). The reported plots are all obtained by setting the typical value of the cardiac fiber contraction (β = 0.3) and for different values of γ. Their meaning is that multiple spikes arise (for suitable values of ε and of the initial condition) if the parameters (α, v thr ) lie below the curve corresponding to the correct value of γ. Figure 6(b) completes the analysis by showing the region in the plane (v thr , ε) where the existence of two-spikes trajectories is guaranteed, for γ = 0.5 and β = 0.3. It is interesting to remark that, by comparing the plots reported in Figures 1 and 6, we observe that the corresponding critical values of the threshold potential are smaller when computed on solutions of (18) than those obtained for (1). This result evidences that taking into account the electro-mechanical coupling in the cardiac evolution limits the possibility of onset of multiple-spikes solutions. As this emerges as one of the most interesting predictions of the present study, we next focus on how the possibility of finding multiple-spikes solutions depends on the choice of the fibershortening parameter β. Figure 7 summarizes the results of our numerical investigations. For several values of β ∈ [0.1, 0.7], and the reported values α = 0.05, 0.10, 0.15, and 0.25, we have identified the value of v thr,cr (β, α) such that if the actual v thr is smaller than the reported critical value, multiple-spikes solutions may exist, while if the actual v thr exceeds the reported v thr,cr no secondary spikes may arise. Clearly the above plots share some common information. To make an example the critical threshold v thr,cr (α) for α = 0.25, which can be read from Figure 6(a) for γ = 0.5, coincides with the abscissa of the right-most point (marked with +) in the α = 0.25 curve in Figure 6(b), and can be retrieved from Figure 7 as the ordinate of the α = 0.25 curve when β = 0.3 (also marked with +).  The physiological interpretation of the results evidenced in Figure 7 is that increasing the fiber-contraction ability helps in preventing the onset of spontaneous after-contractions. Indeed, the critical depolarization threshold decreases when β increases, which in turn increases the probability for the physiological threshold to fall in the no-secondary spikes domain. To better explain this idea we have reported in Figure 7 a horizontal dotted line to evidence the level v thr = 0.125. As we have already pointed out, this is roughly the physiological threshold for a suprathreshold DAD to arise. Figure 7 shows that for the physiological value α = 0.1 the dotted line falls entirely in the no-spikes region. Moreover, as the contraction parameter β increases, smaller and smaller values of α are required to allow for the onset of secondary spikes.
The prediction above is in agreement by the available experimental data. In [25], for instance, the authors evidenced the crucial role of the Na + -Ca 2+ exchanger in determining the concurrent appearance of contractile dysfunctions and arrhythmogenesis in heart failure events. Since an increase in the activity of Na + -Ca 2+ exchanger in the myocyte induces a reduction in supply of calcium to the sarcomere, the contraction ability of the myocyte decreases, corresponding to the decrease of β. Simultaneously, as in heart failures the threshold sarcoplasmic reticulum Ca 2+ load that causes the Ca 2+ -induced Ca 2+ release decreases, the storage of calcium in the sarcoplasmic reticulum is sufficient to induce the outward flux of Ca 2+ , thus contributing to the occurrence of arrhythmias initiated by DADs (in agreement with the growth of v thr,cr in Figure 7).

Conclusions.
We have illustrated how the delayed afterdepolarization event, that may occur in the cardiac action potential, can be modeled in the framework of generalized FitzHugh-Nagumo dynamical systems.
To predict DAD onset, we have defined the number of spikes associated with each trajectory, and studied how such number depends on the constitutive parameters of the model. In particular, we have focused attention on the possibility of finding orbits with more than one spike, becomes of their potential ability to trigger arrhythmia in the heart activity.
Within the classical FN model, we first identified the region in parameter space (see Figure 1) where orbits can be found that exceed the threshold potential. In this simplified setting we have also been able to address the question whether a local impulse is able or not to trigger an aftercontraction, modeled as a traveling impulse. The result we found (see Figure 3) is in perfect agreement with the experimental predictions that the aftercontraction event is a all-or-none mechanism, in the sense that once a spike occurs at a specific position, it is very likely that it will give rise to a traveling signal [18].
We have then moved to the study of more recent electromechanically-coupled generalizations of the FitzHugh-Nagumo model. The main result of this analysis, summarized in Figure 7, is that the mechanical contraction in the heart plays a stabilizing role in preventing the onset of secondary and unwanted impulses, a result in agreement with experimental data as well [25].
Extra-care must be clearly exercised before giving physiological significance to the quantitative predictions reported in our results (e.g., the exact values of the constitutive parameters associated with the onset of extra-spikes). It is indeed well established that there is no hope to obtain an accurate description of the very many processes (chemical, electrical, mechanical, geometrical, and fluido-dynamical, among others) simultaneously occurring in the heart evolution within the framework of a simple dynamical model [26]. Nevertheless, in our opinion some interesting qualitative results emerge from the present analysis. In particular we remark that no secondary spikes emerge when the constitutive parameters are set equal to the values typically chosen as physiological in the literature. Nevertheless, suitable small perturbations from the physiological range, and in particular small decreases in the parameter α, immediately bring extra spikes into play.