A METHOD FOR ANALYZING THE STABILITY OF THE RESTING STATE FOR A MODEL OF PACEMAKER CELLS SURROUNDED BY STABLE CELLS

The purpose of this paper is to derive and analyze methods for examining the stability of solutions of partial differential equations modeling collections of excitable cells. In particular, we derive methods for estimating the principal eigenvalue of a linearized version of the Luo-Rudy I model close to an equilibrium solution. It has been suggested that the stability of a collection of unstable cells surrounded by a large collection of stable cells can be studied by considering only a collection of unstable cells equipped with a Dirichlet type boundary condition. This method has earlier been applied to analytically assess the stability of a reduced version the Luo-Rudy I model. In this paper we analyze the accuracy of this technique and apply it to the full Luo-Rudy I model. Furthermore, we extend the method to provide analytical results for the FitzHugh-Nagumo model in the case where a collection of unstable cells is surrounded by a collection of stable cells. All our analytical findings are complemented by numerical computations computing the principal eigenvalue of a discrete version of linearized models.


(Communicated by Mette Olufsen)
Abstract.The purpose of this paper is to derive and analyze methods for examining the stability of solutions of partial differential equations modeling collections of excitable cells.In particular, we derive methods for estimating the principal eigenvalue of a linearized version of the Luo-Rudy I model close to an equilibrium solution.It has been suggested that the stability of a collection of unstable cells surrounded by a large collection of stable cells can be studied by considering only a collection of unstable cells equipped with a Dirichlet type boundary condition.This method has earlier been applied to analytically assess the stability of a reduced version the Luo-Rudy I model.In this paper we analyze the accuracy of this technique and apply it to the full Luo-Rudy I model.Furthermore, we extend the method to provide analytical results for the FitzHugh-Nagumo model in the case where a collection of unstable cells is surrounded by a collection of stable cells.All our analytical findings are complemented by numerical computations computing the principal eigenvalue of a discrete version of linearized models.

Introduction.
A normal heart beat is initiated in the Sinoatrial (SA) node, which sets up an electrical wave propagating through the entire cardiac muscle and thereby sets off the contraction of the muscle.This contraction is the pumping mechanism of the heart and thus vital for any human being.The process is repeated about three billions times during an average life-time so the machinery is extremely robust, but not without difficulties; every year about 300,000 US citizens die of ventricular fibrillation, which is a complete break-down of the normal propagation of the electrical wave starting in the SA-node.Understanding atrial and ventricular fibrillation is a major scientific challenge and has been a subject of intense research for centuries.Over the past fifty years, mathematical models have been extensively used to understand what happens when the electrical flow turns from a smooth, predictable and well synchronized flow enabling a forceful contraction of the cardiac muscle to a seemingly chaotic situation which disables the contraction and ultimately, in the ventricular case, is lethal.For a comprehensive introduction to mathematical models in this field we refer to [3,11,9].
The main triggering point of electrical activity in the heart is the SA-node.The SA-node consists of cells that are automatic in the sense they keep on firing signals without being stimulated.Also other parts of the heart contain automatic cells, but most of the heart consists of cells that are stable in the sense that if they are in equilibrium, they will remain quiescent until they are stimulated.To some extent extra automatic cells represent a back-up system that will take over the pacemaker role if the SA-node is in trouble, but there are also locations where automatic cells may generate difficulties.For instance, around the pulmonary veins automatic cells have been identified and this is a possible cause of atrial fibrillation; see the review by Khan [5].The ability of a set of automatic cells to initiate waves depend on a number of parameters such as the strength of the automaticity, the number of automatic cells and the electrical coupling between the cells; see [5,10].If such an ectopic wave is generated, it may trigger a re-entrant electrical wave that again may lead to fibrillation; see e.g.Weiss et al. [13,14].
It is generally accepted that the Monodomain equations (see [4]) combined with a system of ordinary differential equations modeling the cell dynamics can model the electrical properties of the cardiac tissue adequately.In such a model, an ectopic wave can be realized as an instability of the resting state.Jacquemet [2] studied the stability of the resting state for single cell models, whereas the spatially coupled case was studied by e.g.Pumir et al. [10] and in the paper [12].In [12] we derived a condition for setting off ectopic waves in models of excitable cells where automatic cells were surrounded by stable cells.This was done by studying when the automatic cells were able to break the resting state, depending on different model parameters.In that paper, we provided numerical and analytical evidence indicating that the stability of the resting state of the problem could be analyzed by studying a simplified problem exclusively consisting of automatic cells.The reasoning was as follows: If the automatic cells are surrounded by strongly stable cells, the stable cells will act like a Dirichlet boundary condition for the automatic cells and thus the problem can be approximated by a reduced model.That reduced model can in turn be linearized and we can derive analytical bound for the stability of the steady solution to the problem.
The aim of the present paper is to assess the validity of the approach presented in [12].More precisely, we consider the Monodomain system written on the form Here v denotes the transmembrane potential, r denotes a vector containing concentrations and gating variables, δ is the spatial diffusion coefficient, and I and R are given non-linear functions explicitly depending on position x (in order to model different cell collections).The system is here equipped with homogeneous Dirichlet boundary conditions on the domain Ω, and the domain consists of two parts; Ω N and Ω A where the first is assumed to surround the latter.In Ω N we assume that all cells are normal in the sense that the equilibrium state of individual cells are stable with respect to small perturbations.More precisely we assume that in the case of δ = 0, there is a state (v 0 N , r 0 N ) which is stable in the sense that the principal eigenvalue of the associated Jacobian has real part less than zero.Similarly, we assume that the individual cells in Ω A are automatic in the sense that a single cell in that region would act as a pacemaker.Furthermore, we assume that there is a steady state solution (v 0 A , r 0 A ) for the single cell case, and that this state is unstable.We let a denote the diameter of Ω A , and recall that the stability of the stationary solution to the coupled system (1) depends on the factor δ/a 2 ; see [4,10,12].
If the automatic cells are indeed capable of stimulating the normal excitable heart cells they should be able to break the resting state.As mentioned previously this ability is dependent on a number of parameters.In this paper we will with analytical techniques investigate how the stability of the steady state depends on different model parameters.Note that stability of the stationary solution indicates that the automatic collection is not able to generate an ectopic wave, since in this case all cells remain quiescent.In the paper [12], we demonstrated that a good estimate for the stability region graphed in the (a, δ)−plane could be achieved by considering the problem (1) defined on Ω A using (v 0 A , r 0 A ) as a Dirichlet-type boundary condition.The purpose of the present paper is to assess the accuracy of this technique.Ideally, one would like to validate the reduced approach on realistic three-dimensional domains.It is certainly possible to obtain numerical evidence for the validity of the method by computer simulations, but this will not give a rigorous justification of its accuracy.Here we aim to assess the accuracy of the method by mathematical analysis of the PDEs.To this end we will need to solve the governing eigenvalue problems analytically to find where the principal eigenvalue crosses the imaginary axis depending on model parameters.In one spatial dimension this is a feasible approach and we will therefore consider one-dimensional problems in the following.
We start by considering a caricature model and strengthen the computational evidence provided in [12] with an analytical assessment of the accuracy.Similar arguments are given for a version of the FitzHugh-Nagumo equations adopted to the setting of the present paper.Finally, we analyze the stability of the Luo-Rudy I model; see [8].We derive a linearized version of the model and analytically compute an estimate of the stability region, and show that this estimate fits very well to the results of numerical computations.
2. Caricature model.In a first stage we consider the Caricature model from [12]; it has the purpose of mimicking the situation of automatic cells surrounded by normal cells.By utilizing techniques from elementary quantum mechanics for solving potential well problems we obtain useful insight into the stability properties of the Caricature model.In the next section we will extend and apply the method to the FitzHugh-Nagumo model.
Consider the Caricature model where δ > 0 is the diffusion coefficient and Here 2a is the size of the automaticity region, 2M a with M ≥ 1 is the domain size, α > 0 corresponds to the degree of automaticity and β > 0 controls the stability of surrounding cells.The PDE (2) has the equilibrium equation and the stability of the null solution to (2) is governed by the corresponding eigenvalue problem δu ′′ + (p(x) − λ)u = 0 , u(±M a) = 0 .
For convenience we introduce the scaling x = x/a and obtain the equivalent eigenvalue problem where 2.1.Infinite domain.Our goal is to gain insight into how the stability properties of (3) depend on the parameters a, α, β and δ; later on, we will consider the M dependence as well.To begin with however, we assume that the number of normal cells is much larger than the number of automatic cells; i.e.M ≫ 1.Thus, we will start by considering (3) on R and look for square integrable solutions to The problem (4) resembles the finite well problem in quantum mechanics and can be analyzed by adapting standard techniques from that field; see e.g.[6].
The eigenvalue problem (4) on the infinite domain can be written (omitting the bars) where the function g is piecewise constant in Here are functions defined on the normal and automatic domains, respectively.Now, let k = a g A (λ)/δ and κ = a g N (λ)/δ and make the following ansatz for the first eigenfunction First we must ensure that the eigenfunction is continuous at x = ±1.This is readily done by putting c = e −κ cos k .
Next, the ansatz function should be differentiable at x = ±1.Due to symmetry it is sufficient to impose the following condition at x = 1 A straight forward calculation gives that κ = k tan k must hold for differentiability.From this condition and the ansatz for k and κ we obtain an equation for the eigenvalue λ which in the present case takes the explicit form In the Caricature model, the stability of the surrounding cells is given by the size of β; the cells become more stable as this value is increased.From (11) we note that as β goes to infinity, we get since arctan(x) → π/2 as x → ∞.This coincides with the principal eigenvalue of the purely automatic problem given by Thus, as stability of the surrounding cells increases we approach the purely automatic case equipped with homogeneous Dirichlet boundary conditions.
2.1.1.Stability: infinite domain.We begin by considering the Caricature model on the infinite domain.Equation ( 4) is self-adjoint and has therefore real eigenvalues.Hence, stability is lost when the largest eigenvalue becomes positive.With λ = 0 in (11) we can solve for η ∞ = δ/a 2 (the subscript indicates the infinite domain) to obtain η ∞ = α The eigenvalue λ in ( 11) changes sign at η ∞ and we can solve the eigenvalue equation (11) for an η = η ∞ to find that λ < 0 for η > η ∞ .This motivates the definition of a stability function δ Stability functions separate the stable and unstable cases.For the purely automatic case ( 13) we obtain from ( 12) that and the corresponding stability function is The stability of the two cases can be compared by considering the distance between the stability functions which is positive since 4 arctan 2 (ξ) monotonically approaches π 2 from zero as ξ grows from zero.So, the case of automatic cells surrounded by an infinite collection of normal cells is less stable than the corresponding purely automatic case, in the sense that it requires a larger diffusion coefficient δ to stabilize.Also note that, since ξ = β/α, the distance between the stability functions decreases with decreasing automaticity α and increasing stability β of the normal cells.
2.2.Finite domain.We turn to the eigenvalue problem (3) that can be written in the form ( 5)-( 6) with boundary conditions u(±M ) = 0.As before we have g N (λ) = β + λ and g A (λ) = α − λ in the function ( 6), and we let k = a g A (λ)/δ and κ = a g N (λ)/δ.The ansatz for the eigenfunction reads Continuity at x = ±1 gives and a differentiability condition on w at x = ±1 leads to .
From this condition and the ansatz for k and κ we obtain the following equation for the eigenvalue λ which in the present case takes the form We remark that homogeneous Neumann conditions can be treated by replacing the sinh-functions in (18) by cosh-functions.
This equation defines implicitly the stability function δ M = η M a 2 through η M , which can be found numerically by using e.g.Newton's method.However, we can show that δ M is bounded below by δ a and bounded above by δ ∞ .To obtain the bounds on δ M we first calculate the limit Comparing this to (14) we conclude that lim M→∞ η M = η ∞ .Similarly by ( 16) and 0 = lim Note from (21 which is positive for the considered parameters.So, (22) holds and therefore δ M is bounded below by δ a and bounded above by δ ∞ .

Numerical tests.
In the present experiment the Caricature model ( 2) is discretized as follows where These discrete equations can be written U ′ = AU , where the matrix A is symmetric and therefore has real eigenvalues.The largest eigenvalue of A determines stability of the discrete system.We perform experiments with N = 400.To assure that enough points are used in the computations in this paper we have run some numerical experiments using larger values of N , e.g.N = 800, and we have found that the output agrees with the results reported here.The stability of the discrete system is recorded for different values of a and δ using certain fixed parameters α, β and M .The results of the simulations are visualized in plots as follows: if the largest eigenvalue is negative the corresponding location is marked with an 'o' to indicate stability, and if not it is marked with an '•' to denote instability.In addition, the stability function δ a = η a a 2 in (17) for the purely automatic case and δ M = η M a 2 for bounded domains are plotted.Recall that η a and η M are the scaled diffusion coefficients for which the principal eigenvalue of the equations ( 13) and ( 3) is zero, respectively.Note also that δ M = η M a 2 is obtained by solving (21) for η M .We will investigate the influence of automaticity α and stability β of normal cells on the stability of (2).
First, we will numerically verify the effect of β on the stability.The upper row in Fig. 1 shows simulations on the domain [−M a, M a] with M = 8 and automatic cells on [−a, a].We observe that as β increases, the purely automatic case approximates the mixed cell problem more accurately.In the lower row, the points η = δ/a 2 where stability changes are graphed as a function of β.The results confirm that the mixed cell problem on [−M a, M a], automatic on [−a, a], has points η M that satisfy η a ≤ η M ≤ η ∞ .That is, the purely automatic case is more stable, and the mixed cell problem on the infinite domain is less stable than the mixed problem on the finite domain.Notice that η a , η M and η ∞ coincide for large β.
Next, the effect of α is considered.From the results visualized in the upper row in Fig. 2 we conclude that the purely automatic case approximates the mixed cell problem more accurately as α decreases.The lower row of the figure depicts η a , η M and η ∞ for the three cases.We observe that η a ≤ η M ≤ η ∞ and that they are closer for smaller α.
We remark that scaling p in ( 3) and ( 4) by a factor q (keeping the ratio β/α fixed) results in the same scaling of the corresponding stability functions δ ∞ and δ M ; put λ = 0 in the equations to see this.A comparison of the top rows in Fig. 1 and Fig. 2 confirms that changing α and β but keeping the ratio β/α fixed, results in a corresponding scaling of δ.  showed that δ = ηa 2 , where a is the size of the automatic region.It is therefore sufficient to study the scaled problem and its scaled diffusion coefficient η.We have determined scaled diffusion coefficients η such that the equilibrium equation of the Caricature model (3) has the principal eigenvalue λ * = 0; thus stability of the null solution changes for these values of η.
It is found that the scaled diffusion coefficient η M for the problem with two cell types on a finite domain satisfies (22), η a ≤ η M ≤ η ∞ , while keeping the cell characteristics fixed.Here, η a and η ∞ are the diffusion coefficients for the purely automatic problem and the mixed cell problem on the infinite domain, respectively.From (10) and ( 12) these coefficients can be written where λ * = 0 and g N and g A are given in ( 7)-( 8).Because of ( 22), the purely automatic case is a better stability indicator for the mixed cell problem on a finite domain when η ∞ ≈ η a .η ∞ and η a will be close when η a is small: Let 0 < g A (λ * ) =: z and fix the normal cells, 0 < g N (λ * ) =: C. Then a Taylor expansion reveals that When g A (λ * ) is small, the scaled diffusion coefficient η a in (24) is small and the automatic cells are not very unstable.In this case the purely automatic case is a good stability indicator for the mixed cell problem on a finite domain.
3. FitzHugh-Nagumo model.In this section we consider the FitzHugh-Nagumo model [1,4].Our objective is to assess the stability of the equilibrium solution when two types of cells are present.In a first step the equilibrium equation with only automatic parameters is studied; this is followed by the case of automatic cells surrounded by normal cells.
3.1.Purely automatic case.The FitzHugh-Nagumo model with diffusion reads and is augmented with the boundary conditions v(±a, t) = 0.Here δ > 0, 1 ǫ ≫ 1, γ > 0 and for the automatic case α > 0. This PDE has the equilibrium equation with v satisfying v(±a) = 0.The stability properties of the equilibrium solution v 0 = w 0 = 0 are uncovered by linearizing the equilibrium equation and solving the corresponding eigenvalue problem with v(±a) = 0.Here the partial derivatives of f are given by so the eigenvalue problem simplifies to By noting that w = v/(λ + γ) and substituting this into the first equation we obtain the scalar equation With the boundary conditions v(±a) = 0 the corresponding eigenvalue equation reads where as before η = δ/a 2 .Stability of the null solution is lost when Re λ > 0.
To find the possible values of η for which the eigenvalues cross the real axis we put λ = ib, with b real, in (28); one complex equation, or two real equations are obtained.The equation for the imaginary part is Solving this equation for b we get and from (28) the corresponding values of η are Note that at η 1 the eigenvalues are purely imaginary.We have found the η-values that have Re λ = 0, so the real part of the eigenvalues can only change sign at these values of η.We can easily find out if both eigenvalues have negative real part (indicating stability) for large η by solving (28) for an η > η a where In the following we will fix the parameters ǫ = 1/100 and γ = 1, and consider the automatic case α ∈ (1/100, 1).These parameters give η 0 = 400(α − 1)/π 2 < 0 for λ = 0, and η 1 = 4(100α − 1)/π 2 > 0 for λ = ib * = ±i √ 99.So, we have and the corresponding stability function is (30)

Automatic cells surrounded by normal cells. The case of automatic cells surrounded by normal cells is modelled by
with homogeneous Dirichlet boundary conditions for the finite domain [−M a, M a], and a condition on square integrability in the case of infinite domain.Here the automaticity parameter p depends on x, Stability is governed by the eigenvalue problem for the linearized equilibrium equation for (31); linearizing around the null solution we find augmented with the appropriate boundary conditions.For convenience we introduce the scaling x = x/a and obtain the equivalent eigenvalue problem x > 1 .By substituting w = v/(λ + γ) in the v-equation we find, on omitting the bars, a scalar equation on the form (5) where g is a function on the form ( 6) Thus, we are now in the same setting as for the Caricature model.Note however, that the functions g A and g N are different for the Caricature and the FitzHugh-Nagumo models.In the present case g A is given in (27) and From the results on the Caricature model in Sect. 2 we have that k = a g A (λ)/δ and κ = a g N (λ)/δ in the eigenfunctions ( 9) and ( 18) on the infinite and finite domains, respectively.Since we are looking for real valued (and square integrable) solutions we must have κ real.Thus, it is required that δκ 2 = a 2 g N (λ) > 0. At the points where the real part of λ changes sign we have with b ∈ R .
To have g N > 0 the imaginary part of the eigenvalue must consequently be b 0 = 0 or b * 2 = 1 ǫ − γ 2 > 0. In order to satisfy the differentiability condition on the eigenfunction, k must be real as well.Otherwise we would have a cosh-function in the ansatz for the eigenfunction and that would make the ansatz non-differentiable at the boundaries between the automatic region and the normal regions.Therefore it is required that also g A is positive.With c ∈ R we calculate , As in the previous subsection we will fix the parameters ǫ = 1/100 and γ = 1; the present mixed cell problem is investigated for α ∈ (1/100, 1) and β > 0. Let us first consider the infinite domain case.

Infinite domain. Recall the eigenvalue equation (10)
where now the functions g A and g N are given by ( 27) and (33).Putting λ = 0 in this equation and solving for η = δ/a 2 we get which is not a real number for the parameter ranges, α ∈ (1/100, 1) and β > 0. For λ = ib * = ±i √ 99 we calculate which is a positive number for the considered parameter ranges.Therefore we have that As in Sect.2.1.1 on the Caricature model, we will compare the stability of the null solution in the case of automatic cells surrounded by an infinite collection of normal cells to the purely automatic case.The distance between the stability function and the stability function (30) for the purely automatic case is which is positive, since 4 arctan 2 (ξ) monotonically approaches π 2 from zero as ξ grows from zero.Thus, as for the Caricature model, the purely automatic case is more stable.The argument of the arctan-function is ξ = √ 100β + 1/ √ 100α − 1, and we conclude that the stability functions are closer for increasing β of the normal cells and decreasing α of the automatic cells, α → ǫ = 1/100 from above.Notice from (29) that η a → 0 as α → ǫ and the cells become less automatic.These results are analogous to the ones for the Caricature model in Sect. 2.
Next, we turn to the case of the finite domain, [−M, M ] after rescaling, with homogeneous Dirichlet boundary conditions at ±M .

Finite domain. Recall the eigenvalue equation (19)
where g A and g N in the present case are given by ( 27) and ( 33).With λ = ib * = ±i √ 99 and η = δ/a 2 we obtain which implicitly determines δ M = η M a 2 through η M ; this η M can be found numerically by using e.g.Newton's method.Following the calculations in Sect.2.2.1 we find the bounds on δ M by first calculating the limit Table 1 quantifies our results.Shown is a comparison of the critical diffusion coefficients obtained for the mixed cell problem (δ M from (36)) and the purely automatic case (δ a from ( 29)) when the model parameters α and β are varied; the parameters ǫ = 1/100, γ = 1 and M = 8 are fixed.Examining Table 1 closely, we may draw some formal conclusions: From the table we find that δ M /δ a → 1 + as 1/β → 0 with the rate ∼ 1/2 for fixed α.Similarly, for fixed β we have that δ M /δ a → 1 + as µ → 0, where µ = α − ǫ, with a rate around 1/2.Furthermore, we observe that (η M − η a ) → 0 + as µ → 0 with a rate ∼ 3/2 for fixed β, while (η M − η a ) → 0 + as 1/β → 0 with a rate ∼ 1/2 for fixed α.  (34) the results (23) and (24) hold with λ * = ib * as well: where g A is given in (27) and g N in (33).The situation is completely analogous to the Caricature model: equation (25) holds with z = g A (λ * ) and We observe from (24) that η a → 0 as z → 0 and conclude from (25) that the automatic cells cannot be very unstable if the purely automatic problem should accurately determine the stability of the null solution to the mixed cell problem.
4. The Luo-Rudy I model.Although interesting from a mathematical point of view, the FitzHugh-Nagumo model does not provide an accurate description of the action potential of a cardiomyocyte.Over the past decades, a large number of increasingly accurate models have been developed; see e.g.[15] where an overview is given.In the present paper we consider the Luo-Rudy I model [8] which is a fairly accurate model, but still simple enough to be analyzed.We are able to analyze the purely automatic case with respect to stability.However, for the problem involving two types of cells we will resort to numerical computations.Our purpose is to find out if the purely automatic case serves as a good stability indicator for the more complicated mixed problem, as has been the case for the previously analyzed simpler models.
Equipped with diffusion the Luo-Rudy I kinetics read where v is the transmembrane potential, c a scaled calcium concentration and g = (x, m, h, j, f, d) T is a vector of gating variables.We are interested in the stability of the equilibrium solution (v 0 , c 0 , g 0 ) T to (37).To this end we augment the equilibrium solution with a small perturbation that vanishes on the boundary, put the result in the Luo-Rudy 1 heart tissue model (37) and linearize around the equilibrium to obtain the relevant eigenvalue problem.This eigenvalue problem governs the stability of the stationary solution and can be written as follows, with homogeneous Dirichlet boundary conditions at x = ±a.Here superscript 0 indicates evaluation at the component-wise constant equilibrium solution (v 0 , c 0 , g 0 ) T , and φ i are the functions given by the right hand side in the ODEs for the gate variables in (37).Introducing the notation u = (v, c, g 1 , . . ., g 6 ) T for the unknowns and letting s ij be the factors multiplying u j in the left hand side of the i th equation of (38), the system can be written with u i (±a) = 0 for all components.Note that the s ij s are nothing but the entries of the Jacobian of the single cell LRI-system, evaluated at the equilibrium solution.Solution components 3 to 8 of u in terms of the first one are given by To find the η-points where the real parts of λ cross zero we can proceed as in Sect.3.1.We start by calculating F (ib, η) = 0 for b, η ∈ R to obtain equations for the real and imaginary parts of (40) F r (b, η) = 0 (41) Here, F r and F i are linear in η and polynomials of degree 8 and 7 in b with only even and odd powers, respectively.We will have 8 + 7 = 15 roots, (b * , η * ), which can be found by using computer algebra software; e.g.Maple.Only real roots are interesting; in particular the one with the largest η * , say η a .Now we can solve (40) for an η bigger than η a to verify that all λ have negative real part for all η > η a .Before we investigate the stability of the PDE (37) we will parameterize the single cell Luo-Rudy I model to describe automatic cells of different strengths.In [12] the extracellular potassium concentration, [K] o = 5.4, and the m-gate rate functions were modified in order to model an automatic cell.The Luo-Rudy I rate functions for a normal cell read [8], where µ = 47.13 and ν = 11.0.To model an automatic cell we follow the suggestions in [12,7] and modify the Luo-Rudy I parameters [K] o , µ and ν in the following way, where α is the parameter controlling the automaticity; note that α = 0 corresponds to the unperturbed Luo-Rudy I model.Let us proceed to find numerical values of η a = δ/a 2 for the PDE (37) on [−a, a] for the purely automatic case parameterized with α.We determine the non-zero Jacobian matrix entries s ij at the ODE equilibrium solutions tabulated in Table 2 and obtain from equation (41) the results listed in Table 3 below.Since the equilibrium solution to the PDE modelling automatic cells surrounded by normal cells is not piecewise constant, we are not able to analytically assess the stability for this case.Points η a = δ/a 2 where stability changes for different automaticities α; the corresponding eigenvalue with the largest real part of the Jacobian is the purely imaginary λ * .In the next section numerical experiments will be performed to see if δ a = η a a 2 is a good approximation of the unknown stability function δ M for (37) on [−M a, M a] with automatic cells on [−a, a].Different degrees of automaticity (controlled by α) of the automatic cells will be investigated.Recall that the previous sections have shown, both for the Caricature model and the FitzHugh-Nagumo model, that δ a is a good approximation of δ M when the automatic cells are not very unstable.37) is linearized around its equilibrium solution and discretized on [−M a, M a], using Dirichlet boundary conditions given by the equilibrium states for the normal cell.On [−a, a] we have a collection of automatic cells.The discretization is similar to the one for the FitzHugh-Nagumo model in Sect.3.3.We now need to compute the non-zero and smooth equilibrium solution to the present PDE for the given parameters.By linearizing around this solution we are able to formulate a semi-discrete system U ′ (t) = CU (t).The eigenvalue with the largest real part of the system matrix C governs stability of the equilibrium solution.We present the results by plotting the stability functions δ = η a a 2 from Table 3, and indicating stability of the discretized system with the symbol 'o'.N = 200 points are used in all simulations, and the domain has M = 8.
In a first experiment three simulations are run on [−M a, M a] with different strength of the automatic collection on [−a, a].A smaller α means smaller η a and weaker automatic cells; the surrounding normal cells have α = α N = 0.The results are depicted in Fig. 5.We observe that the purely automatic case is a better approximation to the mixed cell case when η a (and α) is smaller; just as for the previously studied models.However, in contrast to the other models the purely automatic case does not act as a lower stability bound.We believe the reason for this is the non-constant smooth equilibrium solution to the present PDE.The solution differs from the piecewise constant equilibrium states suggested by the single cell ODE for the normal and automatic cases.This discrepancy is bigger for problems with larger α which will need a larger diffusion coefficient δ to stabilize.We recall from Table 2 that α = 0.803 corresponds to a stable state where the single cell equilibrium state is closer to the ones for the automatic cases.In a second experiment we therefore conduct simulations where the surrounding "normal" cells have α N = 0.803.This will make the PDE equilibrium solution being close to component-wise constant.The results shown in Fig. 6 agree qualitatively with the output for the previous models.The main point is, however, that the purely automatic case gives us a good indication of stability for the mixed cell case, especially so for collections of less automatic cells.5. Summary and conclusions.We have assessed the stability of equilibrium solutions to different equations modelling collections of excitable and automatic cells.Depending on model parameters we have analyzed when the automatic cells are able to break the resting state in the presence of normal excitable cells.Stability of the stationary solution indicates that the automatic collection is not able to generate an ectopic wave.
It is found that for all these models, including Luo-Rudy I, the stability for the purely automatic case is a good indicator of the stability for the case of unstable cells surrounded by normal cells; in particular if the automatic cells are not very unstable and the normal cells are stable.For the Monodomain model with FitzHugh-Nagumo cell dynamics we established quantitatively the approximation error in terms of cell characteristics.An advantage is that this indicator can be extended to multidimensional problems.For a simple Caricature model and the FitzHugh-Nagumo model we found that the mixed cell case is becoming less stable with a larger collection of surrounding normal cells, and that the corresponding purely automatic case gives us a lower stability bound.

Figure 1 .
Figure 1.Caricature model.Varying β.Top: The dashed red line is the stability function for the purely automatic case; the solid blue line is for the domain [−M a, M a], automatic on [−a, a].The marking 'o' indicates that the largest eigenvalue of the Jacobian of the discretized system is negative.Parameter values; α = 0.1, M = 8 and from left to right: β = 1, β = 10 and β = 100.Bottom: The points η = δ/a 2 where stability changes, plotted as a function of β; the three cases are (i) purely automatic (dashed red line), (ii) domain [−M a, M a], automatic on [−a, a] (solid blue line), (iii) infinite domain, automatic on [−a, a] (broken green line).Parameter values; α = 0.1, a = 1 and from left to right: M = 2, M = 4 and M = 8.

Figure 2 .
Figure 2. Caricature model.Varying α.Top: The dashed red line is the stability function for the purely automatic case; the solid blue line is for the domain [−M a, M a], automatic on [−a, a].The marking 'o' indicates that the largest eigenvalue of the Jacobian of the discretized system is negative.Parameter values; β = 1, M = 8 and from left to right: α = 0.1, α = 0.01 and α = 0.001.Bottom: The points η = δ/a 2 where stability changes, plotted as a function of α; the three cases are (i) purely automatic (dashed red line), (ii) domain [−M a, M a], automatic on [−a, a] (solid blue line), (iii) infinite domain, automatic on [−a, a] (broken green line).Parameter values; β = 1, a = 1 and from left to right: M = 2, M = 4 and M = 8.

Figure 3 . 3 . 4 .Figure 4 .
Figure 3. FitzHugh-Nagumo model.Varying β.Top: The dashed red line is the stability function for the purely automatic case; the solid blue line is for the domain [−M a, M a], automatic on [−a, a].The marking 'o' indicates that the largest real part of the eigenvalues of the Jacobian of the discretized system is negative.Parameter values; α = 0.02, M = 8 and from left to right: β = 0.2, β = 2 and β = 20.Bottom: The points η = δ/a 2 where stability changes, plotted as a function of β; the three cases are (i) purely automatic (dashed red line), (ii) domain [−M a, M a], automatic on [−a, a] (solid blue line), (iii) infinite domain, automatic on [−a, a] (broken green line).Parameter values; α = 0.02, a = 1 and from left to right: M = 2, M = 4 and M = 8.

δ a =η a a 2 Figure 5 .
Figure 5. Luo-Rudy I model.Varying α.Normal cells have α N = 0.The dashed red line is the stability function for the purely automatic case.The marking 'o' indicates that the largest real part of the eigenvalues of the Jacobian of the discretized system is negative.Parameter values; M = 8 and from left to right: α = 0.805, α = 0.8045 and α = 0.80445.

2 Figure 6 .
Figure 6.Luo-Rudy I model.Varying α.Surrounding "normal" cells have α N = 0.803.The dashed red line is the stability function for the purely automatic case.The marking 'o' indicates that the largest real part of the eigenvalues of the Jacobian of the discretized system is negative.Parameter values; M = 8 and from left to right: α = 0.805, α = 0.8045 and α = 0.80445.

Table 1 .
FitzHugh-Nagumo model.Comparison of the critical diffusion coefficients δ M (automatic and normal excitable cells) and δ a (only automatic cells) where the stationary solution loses stability for different values of the model parameters.

Table 2 .
[12]e2lists, for different values of α, the single cell equilibrium states and the eigenvalue with largest real part of the corresponding Jacobian.The computation of the equilibrium state is explained in[12].Single cell Luo-Rudy I (ODE).Equilibrium states for different automaticities α; note that α = 0 (unperturbed Luo-Rudy I) and α = 0.803 correspond to stable states since λ + < 0.Here λ + = max{Re λ} of the Luo-Rudy I Jacobian evaluated at the steady state.