Strength-Duration Relationship in an Excitable Medium
SStrength-Duration Relationship in an Excitable Medium
B. Bezekci and V. N. Biktashev ∗ College of Engineering, Mathematics and Physical Sciences, University of Exeter, Exeter EX4 4QF, UK
We consider the strength-duration relationship in one-dimensional spatially extended excitablemedia. In a previous study [1] set out to separate initial (or boundary) conditions leading topropagation wave solutions from those leading to decay solutions, an analytical criterion based onan approximation of the (center-)stable manifold of a certain critical solution was presented. Thetheoretical prediction in the case of strength-extent curve was later on extended to cover a widerclass of excitable systems including multicomponent reaction-diffusion systems, systems with non-self-adjoint linearized operators and in particular, systems with moving critical solutions (criticalfronts and critical pulses) [2]. In the present work, we consider extension of the theory to the caseof strength-duration curve.
CONTENTS
I. Introduction 1A. Motivation 1B. Problem statement 2C. A brief history of the mathematicalapproaches 3D. Aims 4II. Analytical theory 5A. Linear approximation 5General Setting 61. The case of critical nucleus 62. The case of moving critical solution 6B. Quadratic approximation of the stablemanifold 7C. A priori bound for critical nucleus case 8III. Hybrid approach 9A. Ingredients of the one-component systems 9B. Ingredients of the multi-component systems 9IV. One component systems 10A. Zeldovich-Frank-Kamenetsky equation 101. The small-threshold limit and the “fullyanalytical” result 102. Hybrid approach 113. Quadratic theory 11B. McKean equation 121. Model formulation 122. Hybrid approach 123. Linear theory 134. Quadratic theory 13V. Multi component systems 14A. I Na -caricature model 141. Model formulation 142. Hybrid approach 143. Linear theory 15B. FitzHugh-Nagumo system 16 ∗ Corresponding author:[email protected]
1. Model formulation 162. Hybrid approach 163. Linear theory 16C. The modified Beeler-Reuter model 161. Model formulation 162. Hybrid approach 173. Linear theory 17VI. Discussion 18Appendices 19A. Discretization formula for strength-durationcurve 191. Finite Difference Discretization Formula 192. Finite Element Discretization Formula 19Discretization Formula for McKeanequation 19Discretization Formula for the I Na -caricature model 193. Threshold curve 21B. Numerical methods for simplified cardiacexcitation model 211. Discretization formula for the critical front ofthe model 222. Discretization Formula for the LinearizedProblem 233. Discretization Formula for the AdjointLinearized Problem 25References 25 I. INTRODUCTIONA. Motivation
The threshold phenomenon “deals with the minimal,an event, or stimulus just strong enough to be perceivedor to produce a response” [3] and the presence of it “im-poses the restriction on the types of mathematical modelsuitable to describe” biological/chemical systems [4]. Its a r X i v : . [ n li n . PS ] D ec extreme importance can be highlighted through exam-ples. For instance, propagation of excitation in the heartinvolves action potential and threshold value controls ifan applied stimulus is sufficient enough to generate anaction potential. Understanding the mechanisms of ini-tiation of propagating is extremely crucial as successfulpropagation enables continuous electrical and chemicalcommunication between cells and failure may lead to se-rious medical conditions [5]. Threshold phenomenon alsoplays a key role in understanding many age related dis-eases such as Alzheimer and Parkinson. Studies on neu-ronal changes in brain suggest that the threshold hypoth-esis helps to explain “some of the associations betweenclinical and pathological findings” [6, 7].Originally, the term excitability has come to be usedto refer to the “property of living organisms to respondstrongly to the action of a relatively weak external stim-ulus” [8]. A well-known example of excitability is theability of nerve cells to generate and propagate electricalactivity. By definition, an excitable medium is a spatiallydistributed system, each element of which possesses theproperty of excitability and it is usually defined as nonlin-ear reaction-diffusion system, where the reaction term de-fines how the constituents of the system are transformedinto each other, and the diffusion part provides propaga-tion of information [8–10]. There are a wide variety ofareas where the term “excitable medium” has been usedrepeatably for decades in many fields including physical,chemical and biological systems and so on [8, 11–15]. B. Problem statement
We consider the problems of initiation of propagatingwaves in one-dimensional reaction-diffusion systems, ∂ u ∂t = D ∂ u ∂x + f ( u ) , (1)where u : R × R → R d is a d -component reagents field, d ≥
1, defined for x ∈ R and t ∈ R + , vector-function f : R d → R d describes the reaction rates and D ∈ R d × d is the matrix of diffusivity. Equation (1) is assumed todescribe an excitable medium as a system “composed ofelementary segments or cells, each of which possesses thefollowing properties: 1. a well-defined rest state, 2. athreshold for excitation, and 3. a diffusive-type couplingto its nearest neighbors. . . . Stimuli below the thresholdare damped out and produce no persistent change in thesystem, . . . stimuli above the threshold induce the cellto change from its rest state to an excited state.” [16] Aclosely related class are bistable systems: whereas an ex-citable system proper returns to the resting state afterspending some time in the excitable state, a bistable sys-tem remains in the excitable state for ever.A definitive feature of an excitable or bistable mediumis existence of traveling wave solutions of (1). These canbe described by transforming the system of partial differ- ential equations (PDEs) to a moving frame of reference, ξ = x − ct, u ( x, t ) = U ( ξ, t ) . (2)We are interested in the solutions that are stationary inthis frame of reference, for a fixed c , i.e. D ∂ ξξ U + c∂ ξ U + f ( U ) = 0 . (3)If the velocity c = 0, then the traveling wave is called thestanding wave. The traveling wave is a front if U ( −∞ )and U ( ∞ ) exist and different from each other (this istypical for bistable systems), and it is a pulse if U ( ∞ ) = U ( −∞ ) = u r (this happens in excitable systems).Travelling wave solutions of (1) have been a topic ofintense study. For applications, for instance modelling ofbiological media and chemical processes, the question ofparticular importance is emergence of such solutions asa result of a perturbation of the resting state, localizedin space and time. For a problem set on the half-infiniteinterval x ∈ [0 , ∞ ), this can be formalized by initial andboundary conditions u ( x,
0) = u r + U s X ( x ) , Du x (0 , t ) = − I s T ( t ) , (4)where X and T describe the shapes of the initial andboundary profiles, and U s and I s are the strengths ofthose profiles. The cases of non-homogeneous initialcondition and non-homogeneous boundary condition areusually handled separately. In electrophysiological terms,these can be described as follows:1. Stimulation by current: U s = 0 , I s (cid:54) = 0. This isthe case when the current is injected at the bound-ary point x = 0 during some time interval. Fora fixed boundary profile T ( t ), there exist a corre-sponding threshold strength value I ∗ s such that thesolution tends to propagating wave (“ignition”) as t → ∞ whenever I s > I ∗ s , and the solution tendsto resting state (“failure”) otherwise. For a one-parametric family of profiles, parametrized by thestimulus duration t s , the corresponding curve I ∗ s ( t s )is called a strength-duration curve (see fig. 1).2. Stimulation by voltage: U s (cid:54) = 0 , I s = 0. Herethe perturbation is instantaneous at t = 0, but isspread in space. For a fixed initial profile X ( x ),there exist a corresponding threshold strengthvalue U ∗ s such that the solution tends to propa-gating wave as t → ∞ whenever U s > U ∗ s , and toresting state otherwise. For a one-parametric fam-ily of profiles, parametrized by the stimulus extent x s , we shall have the corresponding critical curve U ∗ s ( x s ), called a strength-extent curve.In our previous paper [2] we have analysed some analyt-ical and semi-analytical approaches to description of thestrength-extent curves. In this paper, we focus on thestrength-duration curves. In all specific examples belowwe shall consider a rectangular profile of duration t s , T ( t ) = H( t s − t ) e , (5)where the fixed vector e determines which reagents arebeing injected, and H( · ) is the Heaviside step function. C. A brief history of the mathematical approaches
Mathematically, the problem of determining the con-ditions of initiation of propagating waves in excitableor bistable media is spatially-distributed, nonstationary,nonlinear and has generally no helpful symmetries, so theaccurate treatment is feasible only numerically. However,the practical value of these conditions is so high that ana-lytical answers, even if very approximate, are in high de-mand. Historically, there have been numerous attemptsto obtain such answers, based on various phenomenolog-ical and heuristic approaches. Here we review some ofthese attempts, in chronological order.Phenomenological models describing experimental re-lationship between the minimum stimulus amplitude re-quired to excite an axon and the duration for which thestimulus is applied first appeared well before the physicalmechanisms of biological excitability have been discov-ered. The study of the charge-duration relation was firstcarried out by Weiss [17] who experimentally derived thefollowing linear equation Q = a + bt s , (6)where Q is the threshold charge and a and b are fittedparameters. In his original formula, Weiss did not inter-pret the constants a and b physically and hence they werelater on replaced by rheobasic current τ and chronaxie I rh , so that a = τ I rh and b = I rh [18] so (6) becomes Q = I rh ( τ + t s ) , (7)which is known as the Weiss excitation law for the charge.An empirical equation developed by Lapicque [19, 20])reiterated Weiss’s equation in a different form, for therelation between the pulse strength and duration, i.e. the strength-duration curve. Lapicque observed that thestrength of the current I s required to stimulate an ac-tion potential increased as the duration t s was decreased.Lapicque proposed the following current law for excita-tion I s = I rh (cid:18) τt s (cid:19) (8)which is equivalent to (6) as Q = I s t s .Note that the rheobase current, I rh may be definedas the minimal current amplitude of infinite duration forwhich threshold can be reached and that the chronaxietime of cell, τ refers to the value of the pulse duration t s at twice the rheobase current.An alternative expression for threshold stimulatingcurrent was based on the idea that the nerve cell mem-brane could be represented by a parallel resistance R andcapacitance C . The same Lapicque paper, and also later Blair [21, 22] discussed a speculative model relying on an RC network to formulate the strength-duration curve.This resulted in the strength-duration relationship of theform I s = I rh − exp ( − t s /τ ) . (9)Lapicque-Blair’s exponential strength-duration curvelooks similar to that given by the hyperbolic Weiss-Lapicque law (7), (8), and they nearly fit the samedata [22].Lapicque-Blair’s model combines fairly accurately fitexperimental outcomes with mathematical simplicity.Thus, a number of researchers began to focus on it,among which three important ones are Rashevsky [23],Monnier [24] and Hill [25]. Their results are equivalent upto the interpretation, but Hill’s article is the most cited.Hill examined the relationship among the stimulus, theexcitability of the tissue, and its accommodation, wherethe term “accommodation” [26] is used to describe themembrane potential response to a sufficiently slow in-crease in zthe stimulating current without exciting. Aplausible mathematical description for a speculative dy-namic variable describing the accommodation resulted inHill’s two time-constant model I s = I rh (1 − κ/λ )exp ( − t s /λ ) − exp ( − t s /κ ) , (10)where κ , λ are the time constant of excitation and thetime constant of accommodation, respectively. When λ → ∞ and κ = τ , Hill’s equation (10) reduces toLapicque-Blair’s equation (9).All above approaches were phenomenological and theparameters in the strength-duration relationships were tobe fitted to experimental data rather than derived from“first principles”.Study of spatial aspects of the initiation problem datesback at least to 1937, when Rushton [27] introduced theconcept of the “liminal length”, to represent the idea thatin order to be successful, the stimulating current shouldexcite a sufficiently large portion of the excitable cable.He supported this idea by a mathematical model of thenerve axon, which was of course linear (passive), andthe active character of the membrane and the existenceof the excitation threshold were taken into account in aspeculative, axiomatic manner.The situation of course changed radically afterHodgkin and Huxley have succeeded in producing amathematical model describing the work of a nervemembrane based on experimentally established physicalmechanisms. There were no more need for speculativemodelling, but the real equations were strongly nonlin-ear, apparenly making analytical studies unfeasible andnecessitating use of numerical methods. We note onesuch early study done by Noble and Stein [28], who used asimplification of the Hodgkin-Huxley membrane model toexplore numerically the influence of the membrane acti-vation time and accommodation on the strength-duration u x failure t = 1 t = 10 t = 100 t = 200 t = 300ˆ u u x ignition failure ignition I s t s I ∗ s ( t s ) (a) (b) (c) FIG. 1. Response to a below- and above-threshold initial perturbation in ZFK equation. Parameter values: θ = 0 . U s = 0, t s = 0 . I s = 0 . I s = 0 . curve. They also deduced that in the spatially-extendedcontext, the strength-duration curve is highly dependenton the geometry of the stimulus.The analytical, or at least partly analytical, ap-proaches started looking less unfeasible after the papersby McKean and Moll [29] and Flores [30]. They workedwith one simple class, scalar bistable models. Consid-ering the corresponding PDE on half-line with homoge-neous boundary conditions as a dynamical system in afunctional space, they have identified a critical role ofone special solution, dubbed “standing wave” or, later,“critical nucleus”. This is a stationary, spatially non-uniform, unstable solution, with exactly one unstableeigenvalue. Its special role is that its stable manifoldforms the boundary between basins of attraction of “suc-cessful” and “unsuccessful” outcomes of a stimulation at-tempt.An example of an approach seeking to take advantageof this understanding is the work by Neu et al. [31].They considered a Galerkin projection of the infinite-dimensional dynamical system described by the PDEonto a two-dimensional manifold of spatial profiles, “re-sembling” the shape of a developing excitation wave ona half-line (specifically, they used Gaussians). This re-sulted in a second-order system of ordinary differentialequations (ODEs), in which the critical nucleus is rep-resented by a saddle point, and the boundary betweenthe basins is the stable separatrix of this saddle point.This however still left an open problem of describing an-alytically this separatrix. One possibility was explored byIdris [32, 2.5.2], who approximated this separatrix by thestable space of the saddle, which yielded an analytical ex-pression for the strength-extent curve. This expression,however, produced a result that was only qualitativelycorrect.The next step was using the linear approximation ofthe stable manifold right in the functional space. Thisidea was implemented by Idris and Biktashev [1] and ren-dered surprisingly good approximations of both strength-extent and strength-duration curves.Naturally, this approach is applicable only to systems where there exists a critical nucleus. It obviously ex-cludes, for instance, all excitable models. Hence a ques-tion arises, what if anything is the equivalent of the crit-ical nucleus in such systems. An answer to this ques-tion was proposed in another work by Idris and Bikta-shev [33]. It happens that in excitable systems, the roleof the critical nucleus is played not by stationary, but bymoving solutions with one unstable eigenvalue, “criticalpulses” and “critical fronts”. These are unstable “coun-terparts” of the propagating wave solutions, existence ofwhich was realised long before.Finally in this brief review, the approach of [1] wasextended to the moving critical solutions by Bezekci et al. [2], who also explored the possibility of usingquadratic rather than linear approximation of the sta-ble manifold. However, that paper only considered thestrength-extent curve, which in the dynamical systemsparlance is easier as it is about an autonomous systems.The question of strength-duration curve involves study-ing non-autonomous systems, and it is the subject of thepresent contribution. D. Aims
The purpose of this paper is to quantify the strength-duration curves, as an extension of the study [1]. We in-vestigate how the quality of approximation produced byour method depends on the parameters that define vari-ous test systems. Moreover, we investigate the feasibilityof improving the accuracy by using a quadratic ratherthan a linear approximation of the critical manifold, andrelated problems. Finally, we extend the method to thecase where there are no critical nucleus solutions. This isobserved in multicomponent reaction-diffusion systems,where it has been previously demonstrated that insteadof a critical nucleus, one has unstable propagating waves,such as critical pulses [34] or critical fronts [33].The structure of the paper is as follows. After this in-troductory section, Section II describes the proposed an-alytical approach to the problem of the ignition of prop-agation waves in one-dimensional bistable or excitablemedia, from one-component with a critical nucleus tomulti-component systems having moving critical frontsand pulses, including linear and quadratic approximationof the critical manifold. The strongly non-linear natureof the equations makes it unlikely that the ingredientsof these approximations can be found analytically; thus,a short outline of the numerical techniques used in the“hybrid” approach is given in the section III. The appli-cability of the approach will be illistrated for five differentmodels from one-component examples Zeldovich-Frank-Kamenetsky (ZFK) and McKean detailed in section IV,to multicomponent examples I Na -caricature, FitzHugh-Nagumo (FHN) and Beeler-Reuter (BR) detailed in sec-tion V. Finally, section VI concludes the paper with ashort review of the results and some possible further re-search directions. II. ANALYTICAL THEORY
We aim at classification of the solutions of the system(1) set on x ∈ [0 , ∞ ), t ∈ [0 , ∞ ), supplied with the fol-lowing initial and boundary conditions, u ( x,
0) = u r , Du x (0 , t ) = − I s ( t ) , x, t > , (11)in terms of their behaviour as t → ∞ : whether it willapproach the propagating wave solution (“ignition”) orthe resting state (“failure”). We find it convenient to for-malize the initiation problem as one posed on the wholereal line x ∈ R , ∂ u ∂t = D ∂ u ∂x + f ( u ) + h ( x, t ) , ( x, t ) ∈ R × R + , u ( x,
0) = u r , h ( x, t ) ≡ for t > t s , (12)where the boundary condition at x = 0 in (11) is formallyrepresented by h ( x, t ) = 2 I s T ( t ) δ ( x ) , (13)where δ ( · ) is the Dirac delta function.The principal assumption of our approach is existenceof a critical solution , which is defined as a self-similarsolution, u ( x, t ) = ˆ u ( x − ct ) , = D d ˆ u d ξ + c dˆ u d ξ + f (ˆ u ) , ˆ u ( ∞ ) = u r , ˆ u ( −∞ ) = ˆ u − (14)which is unstable with one unstable eigenvalue. Here ˆ u − may be different from u − but in our examples ˆ u − = u r when u − = u r .Similar to the stable wave solution, there is then awhole one-parametric family of critical solutions, u ( x, t ) = ˆ u ( x − ct − s ) , s ∈ R . (15) Due to this translation invariance, the critical solutionalways has one zero eigenvalue. Hence its stable mani-fold has codimension two, whereas its center-stable man-ifold has codimension one and as such it can partitionthe phase space, i.e. it can serve as a boundary betweenbasins of different attractors. Our strategy is to approxi-mate this center-stable manifol. In the first instance, weconsider a linear approximation, and in selected cases,we also explore the feasibility of the quadratic approxi-mation. A. Linear approximation
Let us rewrite the system (12) in a frame of referencemoving with a constant speed c , so that u ( x, t ) = ˜ u ( ξ, τ ), ξ = x − ct − s , τ = t , ∂ ˜ u ∂τ = D ∂ ˜ u ∂ξ + c ∂ ˜ u ∂ξ + f (˜ u ) + h ( ξ + cτ + s, τ ) , ˜ u ( ξ,
0) = u r . We linearize this equation on the critical solution,which is stationary in the moving frame˜ u ( ξ, τ ) = ˆ u ( ξ ) + v ( ξ, τ ) . (16)The linearization gives ∂ v ∂τ = D ∂ v ∂ξ + c ∂ v ∂ξ + F ( ξ ) v + ˜ h ( ξ, τ ) , v ( ξ,
0) = u r − ˆ u ( ξ ) , (17)where F ( ξ ) = ∂ f ∂ u (cid:12)(cid:12)(cid:12)(cid:12) u =ˆ u ( ξ ) (18)is the Jacobian matrix of the kinetic term, evaluated atthe critical solution and˜ h ( ξ, τ ) = h ( ξ + cτ + s, τ ) (19)is the forcing term as measured in the moving frame ofreference. Equation (17) is a linear non-homogeneousequation, with time-independent linear operator, ∂ τ v = L v + ˜ h , L (cid:44) D ∂ ∂ξ + c ∂∂ξ + F ( ξ ) . (20)For simplicity of the argument, we assume that the eigen-functions of L , L V j ( ξ ) = λ j V j ( ξ ) (21)are simple and form a basis in an appropriate functionalspace, and the same is true for the adjoint L + [35]. An-other assumption, which simplifies formulas and is truefor all examples considered, is that all eigenvalues im-portant for the theory are real. We shall enumerate theeigenpairs in the decreasing order of λ j , so by assump-tion we always have λ > λ = 0 > λ > . . . . Thenthe general solution of problem (17) in that space can bewritten as a generalized Fourier series v ( ξ, τ ) = (cid:88) j a j ( τ ) V j ( ξ ) . (22)The coefficients a j will then satisfy decoupled ODEs,d a j d τ = λ j a j + h j ( τ ) ,a j (0) = (cid:68) W j ( ξ ) (cid:12)(cid:12)(cid:12) v ( ξ, (cid:69) , (23)where h j ( τ ) = (cid:68) W j ( ξ ) (cid:12)(cid:12)(cid:12) ˜ h ( ξ, τ ) (cid:69) , (24)the scalar product (cid:68) · (cid:12)(cid:12)(cid:12) · (cid:69) is defined as (cid:68) a (cid:12)(cid:12)(cid:12) b (cid:69) = ∞ (cid:90) −∞ a (cid:62) b d ξ, and W j are eigenfunctions of the adjoint operator, L + W j = λ j W j , L + = D (cid:62) ∂ ∂ξ − c ∂∂ξ + F (cid:62) ( ξ ) , (25)normalized so that (cid:68) W j (cid:12)(cid:12)(cid:12) V k (cid:69) = δ j,k . (26)The solution of (23) is a j ( τ ) = e λ j τ a j (0) + τ (cid:90) h j ( τ (cid:48) )e − λ j τ (cid:48) d τ (cid:48) . By assumption, λ >
0, and due to translational sym-metry, λ = 0, and the rest of the spectrum is assumedwithin the left half-plane. Since the stimulation is sup-posed to be finite in time, h j ( τ ) ≡ τ > t s . There-fore, the condition of criticality is a ( t s ) = 0which implies a (0) + t s (cid:90) h ( τ (cid:48) ) e − λ τ (cid:48) d τ (cid:48) = 0 , from which we seek to obtain the critical curve based onthis linear approximation. General Setting
Using the definitions of a (0) and h ( τ (cid:48) ), we have, interms of the original model, t s (cid:90) e − λ τ (cid:48) (cid:68) W ( ξ ) (cid:12)(cid:12)(cid:12) h ( ξ + cτ (cid:48) + s, τ (cid:48) ) (cid:69) d τ (cid:48) = (cid:68) W ( ξ ) (cid:12)(cid:12)(cid:12) ˆ u ( ξ ) − u r (cid:69) . (27)The forcing term is defined as h ( x, t ) = 2 I s e H( t s − t ) δ ( x ) , (28)hence (27) gives2 I s t s (cid:90) e − λ τ (cid:48) W ( − cτ (cid:48) − s ) (cid:62) e d τ (cid:48) = ∞ (cid:90) −∞ W ( ξ ) (cid:62) (ˆ u ( ξ ) − u r ) d ξ. (29)This is a finite equation relating t s and I s so, in principle,gives the answer. However, it contains the parameter s which still has to be decided upon. This question occursalready for the strength-extent curve, and we refer thereader to [2] for a detailed discussion. The new issuehere is that for t ∈ [0 , t s ] we are now dealing with a time-dependent problem.
1. The case of critical nucleus
This is the case when c = 0, i.e. the critical solution isstationary, and moreover it is even in x . Then there is anatural choice of s = 0 prescribed by symmetry. Hence,(29) gives the classical Lapicque-Blair formula [19, 21] I s = I rh − e − λ t s , (30)where the rheobase is I rh = λ ∞ (cid:82) W ( ξ ) (cid:62) (ˆ u ( ξ ) − u r ) d ξ W (0) (cid:62) e . (31)
2. The case of moving critical solution
In this case there is no x → − x symmetry and thechoice of s is no longer trivial. Following the ideas dis-cussed in [2], we assume that the linear approximationworks best if the initial value for v is the smallest, and itis the smallest in L if not only a = 0 but also a = 0.This helps to fix s as small shifts of the critical solutionsare equivalent to adding a small amount of V . The onlymodification of this idea for the present case is that weapply this condition not at t = 0, but at the moment fromwhich the system is autonomous, i.e. at t = t s . Employ-ing both a ( t s ) = 0 and a ( t s ) = 0 results in followingsystem of equations I s t s (cid:82) e − λ τ (cid:48) W ( − cτ (cid:48) − s ) (cid:62) e d τ (cid:48) = N , I s t s (cid:82) e − λ τ (cid:48) W ( − cτ (cid:48) − s ) (cid:62) e d τ (cid:48) = N , (32)where the right hand sides N and N are constants,defined entirely by the properties of the model, N l = (cid:68) W l ( ξ ) (cid:12)(cid:12)(cid:12) ˆ u ( ξ ) − u r (cid:69) , l = 1 , . (33)System (32) is a nonlinear system of two equations fortwo unknown parameters, I s and s . By eliminating theparameter I s , we find the compatibility condition as fol-lows: t s (cid:90) (cid:104) N e − λ τ (cid:48) W ( − cτ (cid:48) − s ) (cid:62) e −N e − λ τ (cid:48) W ( − cτ (cid:48) − s ) (cid:62) e (cid:105) d τ (cid:48) = 0 . (34)This can be further simplified by using the followingchange of variable, τ (cid:48) = − ( ζ + s ) c that leads to µ ( s ) (cid:44) N e λ s/c c − ct s − s (cid:90) − s e λ ζ/c W ( ζ ) (cid:62) e d ζ − N e λ s/c c − ct s − s (cid:90) − s e λ ζ/c W ( ζ ) (cid:62) e d ζ = 0 . (35)Finite equation (35) defines the shift s for a given t s , orvice versa. After finding the value of s , one only needs toemploy this value in one of the compatibility conditionsin (32) in order to find the amplitude I s since both pro-duce the same result. This completes the construction ofthe linear approximation of the strength-duration curve I s ( t s ). B. Quadratic approximation of the stable manifold
In this subsection, we restrict consideration to the caseof a critical nucleus.For simplicity, rather than using the matrix notationas in the linear approximation, we shall now proceed withan explicit notation for the components of the reaction-diffusion systems. We use Greek letters for superscripts to enumerate them, and adopt Einsteins summation con-vention for those indices. In this way we start from thegeneric reaction-diffusion system ∂u α ∂t = D αβ ∂ u β ∂x + f α ( u β ) + 2 I s e α H( t s − t ) δ ( x )then consider the deviation v α of the solution u α fromthe critical nucleus ˆ u α , u α ( x, t ) = ˆ u α ( x ) + v α ( x, t ) , the equation defining the critical nucleus, D αβ ∂ ˆ u β ∂x + f α (ˆ u ) = 0and the Taylor expansion of the equation for the devia-tion, ˙ v α = D αβ v βxx + f α,β (ˆ u ) v β + f α,βγ (ˆ u ) v β v γ + 2 I s e α H( t s − t ) δ ( x ) + . . . , where overdots denote differentiation with respect totime, subscripts ( · ) x denote differentiation with respectto space and Greek subscripts after a comma designate apartial differentiation by the corresponding reactive com-ponents. The right and left eigenfunctions are definedrespectively by D αβ ∂ xx V βj ( x ) + f α,β ( x ) V βj ( x ) = λ j V αj ( x )and D βα ∂ xx W βj ( x ) + f β,α ( x ) W βj ( x ) = λ j W αj ( x )where j ∈ { , , , · · · } , and the biorthogonality conditionis (cid:68) W j (cid:12)(cid:12)(cid:12) V k (cid:69) (cid:44) ∞ (cid:90) −∞ W αj ( x ) V αk ( x ) d x = δ j,k . We consider only even solutions, so in subsequent sumsonly those j that correspond to even eigenfunctions areassumed. We seek solutions in the form of generalizedFourier series in the right eigenfunctions, v α ( x, t ) = (cid:88) j a j ( t ) V αj ( x )where the Fourier coefficients are defined by a j ( t ) = (cid:68) W j ( x ) (cid:12)(cid:12)(cid:12) v ( x, t ) (cid:69) (cid:44) ∞ (cid:90) −∞ W αj ( x ) v α ( x, t ) d x. Time-differentiation of this gives˙ a j ( t ) = (cid:68) W αj ( x ) (cid:12)(cid:12)(cid:12) ˙ v α ( x, t ) (cid:69) = λ j a j + (cid:88) m,n Q jm,n a m a n + 2 I s E j H( t s − t ) (36)where Q jm,n = Q jn,m (cid:44) ∞ (cid:90) −∞ W αj ( x ) f α,βγ (ˆ u ( x )) V βm ( x ) V γn ( x ) d x, (37)and E j = W αj (0) e α . (38)We assume that eigenvalues are real and ordered fromlarger to smaller, λ > λ = 0 is of course the eigen-value corresponding to the translational symmetry andan odd eigenfunction V = ˆ u (cid:48) , and λ j < j ≥ A j (cid:44) a j (0) = ∞ (cid:90) −∞ W αj ( x ) v α ( x,
0) d x (39)that would ensure that a ( ∞ ) = 0 , which means that the trajectory approaches the criticalnucleus, so the initial condition is precisely at the thresh-old.Let us rewrite the system (36) as an equivalent systemof integral equations, a j ( t ) =e λ j t (cid:20) A j + 2 I s E j λ j (cid:16) − e − λ j ( t ∧ t s ) (cid:17) + t (cid:90) e − λ j t (cid:48) (cid:88) m,n Q jm,n a m ( t (cid:48) ) a n ( t (cid:48) ) d t (cid:48) where we use the notation ( a ∧ b ) (cid:44) min( a, b ). Succes-sive approximations to the solution can be obtained bydirect iterations of this system, a ( i +1) j ( t ) =e λ j t (cid:20) A j + 2 I s E j λ j (cid:16) − e − λ j ( t ∧ t s ) (cid:17) + t (cid:90) e − λ j t (cid:48) (cid:88) m,n Q jm,n a ( i ) m ( t (cid:48) ) a ( i ) n ( t (cid:48) ) d t (cid:48) . Taking a (0) j = 0 for all j , we have a (1) j ( t ) = e λ j t (cid:20) A j + I s E j λ j (cid:16) − e − λ j ( t ∧ t s ) (cid:17)(cid:21) . The requirement a (1)1 ( ∞ ) = 0 recovers the linear approx- imation. The next iteration produces a (2) j ( t ) = e λ j t (cid:26) A j + I s E j λ j (cid:16) − e − λ j ( t ∧ t s ) (cid:17) + t (cid:90) d t (cid:48) e − λ j t (cid:48) (cid:88) m,n Q jm,n × e λ m t (cid:48) (cid:20) A m + I s E m λ m (cid:16) − e − λ m ( t (cid:48) ∧ t s ) (cid:17)(cid:21) × e λ n t (cid:48) (cid:20) A n + I s E n λ n (cid:16) − e − λ n ( t (cid:48) ∧ t s ) (cid:17)(cid:21)(cid:27) . Note that e ( λ m + λ n − λ ) t → t → ∞ because λ m,n ≤ λ < m, n ≥
3, so upon exchanging the orderof intergration and summation, we have converging im-proper integrals. The requirement a (2)1 ( ∞ ) = 0 leads toa quadratic equation for I s , ζ I s2 + ζ I s + ζ = 0 , (40)where ζ = 4 (cid:88) m,n Q m,n (cid:20) E m E n λ m λ n (cid:26) − e − λ t s λ − e ( λ n − λ ) t s − λ n − λ − e − λ t s − e ( λ m − λ ) t s − e ( λ n − λ ) t s + 1 λ m + λ n − λ − e ( λ m − λ ) t s − λ m − λ (cid:27)(cid:21) ,ζ = − E (cid:0) e − λ t s − (cid:1) λ + 2 (cid:88) m,n Q m,n (cid:20) A m E n λ n (cid:26) e ( λ m − λ ) t s − λ m + λ n − λ − e ( λ m − λ ) t s − λ m − λ (cid:27) + A n E m λ m (cid:26) e ( λ n − λ ) t s − λ m + λ n − λ − e ( λ n − λ ) t s − λ n − λ (cid:27)(cid:21) ,ζ = A − (cid:88) m,n Q m,n A m A n λ m + λ n − λ . C. A priori bound for critical nucleus case
We conclude this section with an a priori bound forthe strength-duration curve, obtained by Mornev [36]. Itis applicable to scalar equations, d = 1, u = (cid:0) u (cid:1) , f = (cid:0) f (cid:1) ,such that f ( u j ) = 0, j = 1 , , u r = u < u < u = u − , f (cid:48) ( u , ) < f (cid:48) ( u ) >
0. Then I ∗ s ( t s ) (cid:38) I ∗ s , t s → ∞ , (41)where I ∗ s = max | ˆ u (cid:48) ( x ) | = | ˆ u (cid:48) ( x ∗ ) | = (cid:18) − (cid:90) u u f ( u ) d u (cid:19) / , (42)and x ∗ is the coordinate of the inflexion point of thegraph of ˆ u ( x ), i.e. ˆ u ( x ∗ ) = u . III. HYBRID APPROACH
With a few exceptions, the ingredients for the expres-sions used in the linear and quadratic approximations ofthe strength-duration curve, starting from the critical so-lution itself, are not available analytically and have to befound numerically. In this section we describe numericalmethods we used to find these ingredients. We dividedthis section into two subsections. for models with self-adjoint and non-self-adjoint linearization operator.
A. Ingredients of the one-component systems
This corresponds to the case of the critical nucleus andfor the linear approximation, we need to have the knowl-edge of ˆ u , λ , V while for the quadratic approximation,ideally the whole spectrum of λ (cid:96) , V (cid:96) , (cid:96) = 1 , , , · · · isneeded. Here we shortly describe the methods to obtainthe mentioned ingredients in algorithmic forms. For amore detailed explanation, see [2].In order to find the critical nucleus, we take advantageof the fact that its center-stable manifold has codimen-sion one, and divides the phase space into two open sets,one leading to successful initiation and the other to de-cay [1, 30, 31, 34, 37, 38]. This means that the criticaltrajectories, corresponding to the stimulus strength ex-actly equal to the threshold, tend towards the criticalnucleus as t → ∞ , whereas the stimulus strength slightlyabove or below the threshold produces the solution thatgets close to the critical nucleus and stays in its vicinityfor a long time, before deviating from it to propagate orto collapse. This can be used to calculate an approx-imation of the critical nucleus, see Algorithm 1. Thecalculations are done for (1,11) for x ∈ [0 , L ], where L ischosen large enough for the results to be not significantlydifferent from thos for x ∈ [0 , ∞ ). Input:
Pre-found value of I ∗ s for a chosen t s Output:
Critical nucleus ˆ u • Find u ( x, t ) by solving initial value problem (1,11). • S ( t ) ← || ˙ u || L = (cid:82) L u t ( x, t ) d x • t ← argmin( S ( t )) . • ˆ u ( x ) ← u ( x, t ) Algorithm 1:
Numerical critical nucleus by“shooting”.We calculate the eigenpairs of the linearized operator L defined by (20) (note that in the present case c = 0) usinga variant of the power iteration method. We use randoma number generator to assign linear independent initialguesses for V , V , . . . , choose a time domain t ∈ (0 , T ),and then follow Algorithm 2 until the desired convergencecriterion is fulfilled. Note that the algorithm is describedwithout prejudice to the choice of the norm used for thenormalization; if (cid:107) V (cid:107) = (cid:107) V (cid:107) L = (cid:68) V (cid:12)(cid:12)(cid:12) V (cid:69) / then obvi-ous simplifications are possible. Also, we use the conver- gence criterion based on the change in each eigenvalue; itcan also be done in terms of the change, say in L -norm,in the each eigenfunction. Input:
A linearly independent set (cid:0) V , V , · · · , V n (cid:1) Output: ( λ , V ) , ( λ , V ) , · · · , ( λ n , V n ) • λ , λ , · · · , λ n ← • i ← repeat • i ← i + 1 for k = 1 , , . . . , n do • Solve IVP: V ik ← exp ( L T ) V i − k • Orthogonalize: V ik ← V ik − k − (cid:80) m =1 (cid:42) V ik (cid:12)(cid:12)(cid:12) V im (cid:43)(cid:42) V im (cid:12)(cid:12)(cid:12) V im (cid:43) V im • Normalize: V ik ← V ik (cid:107) V ik (cid:107)• Eigenvalue: λ ik ← T ln (cid:16)(cid:68) V ik (cid:12)(cid:12)(cid:12) V i − k (cid:69)(cid:17) end foruntil | λ ik − λ i − k | ≤ tolerance ∀ k Algorithm 2:
Numerical computation of n principaleigenpairs of self-adjoint operator L by “marching”. B. Ingredients of the multi-component systems
The non-stationary critical solutions, observed inmulti-component systems, can be using an appropriatemodification of Algorithm 1, exploiting computations ina co-moving frame of reference, as described in [2]. How-ever, more accurate results can be obtained by continua-tion of the boundary-value problem (14), an autonomoussystem for vector-function ˆ u ( ξ ) and scalar c . We aim tocalculate conduction velocity restitution curve [39], thatis, a one-parametric family of solutions of the followingperiodic boundary-value problem: = D d u P d ξ + c P d u P d ξ + f ( u P ) , u P ( ξ + P ) ≡ u P ( ξ ) , (43)where P is the spatial period of the waves. Whenthe problem is well posed, (43) defines a curve in the( P, c P , u P ( ξ )) space. In the limit P → ∞ , this curvesplits into two branches, the upper branch with a sta-ble propagating pulse solution, ( c w , u w ( ξ )) and the lowerbranch with an unstable critical pulse solution, ( c, ˆ u ( ξ )),which is of interest to us. We performed the continua-tion using AUTO [40]. To obtain the periodic solutions,we consider an extension of (43) by an extra parame-ter corresponding to “stimulation current” added to thetransmembrane voltage equation. Starting from an ini-tial guess of c w , the continuation is done in accordancewith Algorithm 3.As a final step, we calculate the left and right eigen-functions by means of the Gram-Schmidt orthogonaliza-0 Input:
An initial guess of c w . Output: c , ˆ u ( ξ ) . D d u P d ξ + c P d u P d ξ + f ( u P ) + I ext e = 0 , u P ( ξ + P ) ≡ u P ( ξ ) . • Equilibrium, resting state: I ext ← , c ← c w , ˆ u ( ξ ) ← u r . • Continue the equilibrium by increasing I ext , until aHopf bifurcation is reached. • Continue the periodic orbit from the Hopf bifurcationin the ( c P , P ) plane, down by c P until the fold isreached. • Continue the periodic orbit in the ( I ext , c P ) plane,down by I ext until I ext = 0 is reached. • Continue the periodic orbit in the ( c P , P ) plane bothways. • For the branch with smaller c , select a sufficientlylarge P , and take c ← c P , ˆ u ( ξ ) ← u P ( ξ ) , in a suitablychosen interval of ξ . Algorithm 3:
Numerical critical pulse by AUTO.tion process, modified with account of the fact ath L isnow not self-adjoint, Algorithm 4. Input:
A linearly independent set (cid:0) V , W (cid:1) , (cid:0) V , W (cid:1) , · · · , (cid:0) V n , W n (cid:1) Output: ( λ , V , W ) , ( λ , V , W ) , . . . , ( λ n , V n , W n ) • λ , λ , · · · , λ n ← • i ← repeat • i ← i + 1 for k = 1 , , . . . , n do • Solve IVP: V ik ← exp ( L T ) V i − k • Solve IVP: W ik ← exp (cid:0) L + T (cid:1) W i − k • Biorthogonality: V ik ← V ik − k − (cid:80) m =1 (cid:42) W im (cid:12)(cid:12)(cid:12) V ik (cid:43)(cid:42) W im (cid:12)(cid:12)(cid:12) V im (cid:43) V im • Biorthogonality: W ik ← W ik − k − (cid:80) m =1 (cid:42) W ik (cid:12)(cid:12)(cid:12) V im (cid:43)(cid:42) W im (cid:12)(cid:12)(cid:12) V im (cid:43) W im • Normalization: V ik ← V ik (cid:107) V ik (cid:107)• Normalization: W ik ← W ik (cid:107) W ik (cid:107)• Eigenvalue: λ ik ← T ln (cid:16)(cid:68) V ik (cid:12)(cid:12)(cid:12) V i − k (cid:69)(cid:17) end foruntil | λ ik − λ i − k | ≤ tolerance ∀ k Algorithm 4:
Numerical computation of n principaleigenpairs of non-self-adjoint operator L by“marching”. IV. ONE COMPONENT SYSTEMSA. Zeldovich-Frank-Kamenetsky equation
Our first application example is the one-componentreaction-diffusion equation, first introduced by Zeldovichand Frank-Kamenetsky (ZFK) [41] to describe propaga-tion of flames; it is also known as “Nagumo equation” [37]and “Schl¨ogl model” [42]: d = 1 , D = (cid:0) (cid:1) , u = (cid:0) u (cid:1) , f ( u ) = (cid:0) f ( u ) (cid:1) , f ( u ) = u ( u − θ )(1 − u ) , (44)where we assume that θ ∈ (0 , / u = (cid:0) ˆ u (cid:1) for this equation can be found analyti-cally [1, 30] [43]ˆ u ( x ) = 3 θ √ θ ) √ x √ θ ) √ − θ + 2 θ . (45)The other two components required for the definition ofcritical curves in the linear approximation are λ and W = V = (cid:0) V (cid:1) which are solutions ofd V d x + (cid:0) − u + 2( θ + 1)ˆ u − θ (cid:1) V = λ V ,λ > , V ( ±∞ ) = 0 . (46)We have been unable to find solution of this eigenvalueproblem analytically. We note, however, that ˆ u given by(45) is unimodal, therefore ˆ u (cid:48) , which is the eigenfunctionof L corresponding to λ = 0, has one root, hence bySturm’s oscillation theorem, ˆ u (cid:48) = V and λ = 0, andthere is indeed exactly one simple eigenvalue λ > V solving (46) has no roots.
1. The small-threshold limit and the “fully analytical”result
In this subsection we extend the results of [1] in the pa-rameter space and correct some typos found in the paper.For θ (cid:28)
1, the critical nucleus (45) is O ( θ ) uniformly in x , and is approximatelyˆ u ( x ) ≈ θ x √ θ ) = 32 θ sech ( x √ θ/ . (47)In the same limit, the nonlinearity can be approximatedby f ( u ) ≈ u ( u − θ ). With these approximations, problem(46) has the solution λ ≈ θ, V ≈ sech ( x √ θ/ , (48)and (30) then gives an explicit expression for thestrength-duration curve in the form I s = I rh − e − t s /τ (49)1 -16 -12 -8 -4 tSθ = 0 . . . . . θ t I ∗ s .
05 145 .
76 0 . .
15 62 .
34 0 . .
25 48 .
14 0 . .
35 48 .
17 1 . .
45 67 .
70 15 . ˆ u θ ≪ V θ = 0 . . . . . -1.5-1-0.5 0 0.5 1 0 5 10 15 xV xV FIG. 2. ZFK ingredients for different threshold parame-ters. Top row: The illustration of the typical function S ( t )along with at near-threshold boundary conditions. Middlerow: Comparison between analytical and numerical criticalnuclei and first eigenfunction. Bottom row: Second and thirdeigenfunctions obtained using Gram-Schmidt orthogonaliza-tion method. Parameters: θ = 0 . , . , . , , . x = 0 .
03, ∆ t = 4∆ x / L = 30, t s = 1 .
5, tolerance = 10 − . with the following rheobase and chronaxie form [1] I rh = 4564 πθ / , τ = ( λ ) − = 45 θ . (50)Fig. 3(a) illustrates this approximate strength-durationcurve, compared to the direct numerical simulations. Forchosen parameter values, the comparison is significantlybetter with smaller values of θ , as expected.Remark that the performance of the resulting ap-proximation based on the analytical expression for thestrength-duration curve (49) and (50) can be further im-proved by obtaining the essential ingredients numerically.This is done considering (49), in which the rheobase is,instead, defined according to (31), as I rh = λ ∞ (cid:82) V ( x )ˆ u ( x ) d xV (0) . (51)The plot of the hybrid numeric-asymptotic predictionis compared with the direct numerical simulations asshown in fig. 3(b). As depicted in the figure, reasonableagreement between the two data sets is observed whenthe threshold parameter is small. I ∗ s
0 3 6 9hybridnumerical t s (a) (b) FIG. 3. Strength-duration curves for the ZFK model, for θ = 0 .
05, 0 .
15, 0 .
25, 0 .
35, 0 .
45 (bottom to top), compari-son of direct numerical simulations (lines with symbols) withtheoretical predictions (dashed lines), (a) for the explicit an-alytical answers in the θ (cid:28) x = 0 . t = 4∆ x / L = 100, tolerance = 10 − . It should be noted that the strength-duration curveapproximation remains above the a priori lower bound(42) I ∗ s = − θ (cid:90) u ( u − θ )(1 − u ) d u / = θ / (cid:112) (2 − θ ) √ , for all t s .
2. Hybrid approach
Numerical computation of the essential ingredientsneeded for linear approximation of the critical curves iscarried out using Algorithms 1 and 2 described in III A.Fig. 2 illustrates the processes of numerical computationof the critical nucleus and the first eigenmodes in theZFK equation, for the threshold parameter varying from0 .
05 to 0 .
45 with the increment 0 .
1. The stimulation isdone by fixing the duration time at the value t s = 1 . S ( t ) and t = argmin( S ( t )),the bisection loop is terminated as soon as the abso-lute difference between upper and lower estimate for thethreshold is sufficiently small, i.e. | I s − I s | < − . Foreach case, the solution u ( x, t ) of the nonlinear prob-lem u t = u xx + f ( u ) provides an estimate of the criticalnucleus.
3. Quadratic theory
To estimate a few principal eigenmodes of the ZFKequation, we have considered a finite interval x ∈ [0 , L )as an approximation of x ∈ [0 , ∞ ). We find only a fewapproximate eigenvalues in the discrete spectrum, whilethe remaining ones are in the continuous spectrum. Theeigenfunctions corresponding to the discrete eigenvaluesare well localized whereas those corresponding to the2 I ∗ s t s FIG. 4. Quadratic approximation of the strength-extent curvefor ZFK model for θ = 0 .
05 (green line) compared with di-rect numerical simulations (lines with symbols) and linear ap-proximation (blue line). Parameters: L = 100, ∆ x = 0 . t = 4∆ x / continuous eigenvalues are evidently non-localized, andthus they cannot be taken into account in quadratic ap-proximation. We observe that at increasing values of L ,the eigenfunctions V and V corresponding to the dis-crete eigenvalues λ and λ are well localized towardsthe left end of the interval x ∈ [0 , L ], whereas those corre-sponding to the continuous eigenvalues are evidently non-localized, i.e. vary significantly throughout x ∈ [0 , L ].Thus, for the quadratic approximation, we retain in(40) only the leading term. Setting n = m = 3 gives aclosed expression for the critical curve in the strength-duration plane, ζ I ∗ s 2 + ζ I ∗ s + ζ = 0 , (52)where ζ = 4 Q , E λ (cid:26) − e − λ t s λ − ( λ − λ ) t s − λ − λ − e − λ t s − ( λ − λ ) t s + 12 λ − λ (cid:27) ,ζ = 4 Q , A E (cid:0) − e ( λ − λ ) t s (cid:1) (2 λ − λ ) ( λ − λ ) − E (cid:0) e − λ t s − (cid:1) λ ,ζ = − Q , A λ − λ + A , the coefficients in which are defined by (37), (38) and(39). Fig. 4 shows the comparison between quadraticapproximation of the critical curves and the numericalcurves. Compared to the linear approximation, one cansee some significant improvement for θ = 0 .
05, 0 .
15, 0 . .
35, while the discrepancy between analytical and nu-merical results continues for θ = 0 . B. McKean equation
1. Model formulation
Our second example is a piece-wise linear version of theZFK equation, considered by McKean in [37] and then also in [44]: d = 1 , D = (cid:0) (cid:1) , u = (cid:0) u (cid:1) , f ( u ) = (cid:0) f ( u ) (cid:1) , f ( u ) = − u + H( u − a ) , (53)where we assume that a ∈ (0 , / u ( x ) = − (1 − a ) cosh( x )cosh( x ∗ ) , x ≤ x ∗ ,a exp( x ∗ − x ) , x ≥ x ∗ , (54)where x ∗ = 12 ln (cid:18) − a (cid:19) (55)obtained by the fact that ˆ u ( x ) and its derivative are con-tinuous at this point. The eigenvalue problem can beexpressed as L V = λV (56)where the linearization operator contains the Dirac deltafunction: L (cid:44) ∂ ∂ξ − − a δ ( x − x ∗ ) . (57)The principal eigenvalue and the corresponding eigen-function can be written in the form λ = − κ ,V = cosh ( κx )cosh ( κx ∗ ) , x ≤ x ∗ , exp ( κ ( x ∗ − x )) , x ≥ x ∗ , (58)where κ = 12 a + 12 x ∗ W (cid:16) x ∗ a e − x ∗ /a (cid:17) (59)and W ( · ) is the principal branch of the Lambert W-function as defined e.g. in [45].
2. Hybrid approach
In this model, since the exact analytical solution forthe critical nucleus and the ignition eigenpair are knownfor an arbitrary a ∈ (0 , / -10 -8 -6 -4 -2 ttSS a = 0 . . . . . a t I ∗ s .
05 1 .
19 0 . .
15 1 .
48 0 . .
25 2 .
34 0 . .
35 3 .
94 1 . .
45 7 .
88 1 . ˆ u xa = 0 . . . . .
0 2 4 0.20.40.60.81 xV a = 0 . . . . . FIG. 5. Illustration of the numerical computation of the crit-ical nucleus and ignition mode by “shooting” and “march-ing” in McKean. Top panel: Typical function S ( t ) at nearthreshold boundary conditions. Bottom panel: Critical nu-cleus solutions and ignition modes of the McKean model(53) for various values of the parameter a . Parameters: a = 0 . , . , . , . , .
45, ∆ x = 0 .
03, ∆ t = 4∆ x / L = 10, t s = 1 .
2, tolerance = 10 − .
3. Linear theory
Linear approximation of the strength-duration curvecan be found using the analytically derived expressiongiven by (30). However, it must be noted that in thiscase, the rheobase is found as I rh = λ N , where N = sinh ( κx ∗ ) κ + aκ + 1 cosh ( κx ∗ ) − − a x ∗ ) (cid:18) sinh (( κ + 1) x ∗ ) κ + 1 + sinh (( κ − x ∗ ) κ − (cid:19) . This linear prediction formalism compared with the di-rect numerical simulations is depicted in fig. 6(a). The a priori bound for these chosen threshold parameters isoutside of the duration domain. As shown in the figure,the linear approximation for parameter a values close to1 / a . Even for larger a , there are still somedeviations between the linear theory and numerical sim-ulations, which can be reduced by considering second or-der approximation that will be outlined in the followingsubsection. I ∗ s I ∗ s t s (a) (b) FIG. 6. Strength-duration curves in McKean model: directnumerical simulations (red circles) vs (a) linear theory, for a = 0 .
35 at the bottom , a = 0 . a = 0 .
45 to a = 0 .
48 atthe top, and (b) linear and quadratic theories, for a = 0 . ?? ).Green short-dashed lines: the predictions given by quadratictheory. Discretization: ∆ x = 0 .
03, ∆ t = 4∆ x / L = 10,tolerance = 10 − .
4. Quadratic theory
Linear approximation of the strength-duration curvecan be improved by considering the second order approx-imation. For the quadratic approximation, the knowl-edge of the whole spectrum is ideally required. We knowthat the linearization spectrum of the critical nucleus hasonly one unstable eigenvalue λ >
0, and due to trans-lational symmetry, λ = 0, and the rest of the spectrumlies entirely in the left half-plane. The first aim of thispart of the subsection is to find these remaining stableeigenvalues and corresponding eigenfunctions. To obtainthese eigenpairs, we replace the infinite interval [0 , ∞ )with a finite interval [0 , L ] with a homogeneous Dirich-let boundary condition at x = L , aiming to consider thelimit L → ∞ .The solution of eigenvalue problem (56) in this case is V ( x ; λ ) = ρ cos ( ρx ) − cos ( ρx ∗ ) sin ( ρ ( x − x ∗ )) H( x − x ∗ ) where ρ = √− − λ and the eigenvalues are expressedin terms of the following transcendental equation, h ( ρ ) = tan ( ρL ) − tan ( ρx ∗ ) − aρ cos ( ρx ∗ ) = 0 , (60)which is dependent of the domain size L as opposed tothe eigenfunction expression. After finding analytical ex-pressions for the eigenpairs, the next step is to calculate(37), (38) and (39), and then substitute them back in thecoefficients of the quadratic equation for I s (40) so thatthe second-order approximation of the strength-durationcurve can be generated.Fig. 6(b) shows graphs of the linear and quadratic ap-proximation of the strength-duration curve along withits numerical result for a = 0 .
4. The quadratic approx-imation was obtained for L = 10 and 287 eigenvalues.The accuracy of the second-order approximation is muchcloser to the direct numerical simulation compared to thefirst-order approximation.4 V. MULTI COMPONENT SYSTEMSA. I Na -caricature model
1. Model formulation
Our next example is the caricature model of an I Na -driven cardiac excitation front suggested in [46]. It is atwo-component reaction-diffusion system (1) with u =( E, h ) (cid:62) , D = (cid:18) (cid:19) and f = ( f E , f h ) (cid:62) , where f E ( E, h ) = H( E − h,f h ( E, h ) = 1 τ (H( − E ) − h ) , (61)and H( · ) is the Heaviside step function. The component E of the solution corresponds to the nondimensionalizedtransmembrane voltage, and the component h describesthe inactivation gate of the fast sodium current, which isknown in electrophysiology as I Na and which is mainlyresponsible for the propagation of excitation in cardiacmuscle in the norm.A special feature of this model is that there is a con-tinuum of potential resting/pre-front states, u r = lim ξ →∞ ˆ u = (cid:0) − α, (cid:1) (cid:62) , α > , and a continuum of potential post-front states, u − = lim ξ →−∞ ˆ u = (cid:0) ω, (cid:1) (cid:62) , ω > , so any front solution connects a point from one contin-uum to a point from another continuum. The criticalsolution ˆ u = ( ˆ E, ˆ h ) (cid:62) is described byˆ E ( ξ ) = ω − τ c τ c e ξ/ ( τc ) , ξ ≤ − ∆ , − α + α e − cξ , ξ ≥ − ∆ , ˆ h ( ξ ) = e ξ/ ( τc ) , ξ ≤ , , ξ ≥ , (62)where the post-front voltage ω and front thickness ∆ aregiven by ω = 1 + τ c (1 + α ) , ∆ = 1 c ln (cid:18) αα (cid:19) , (63)and the front speed c is defined by an implicit equation τ c ln (cid:18) (1 + α )(1 + τ c ) τ (cid:19) + ln (cid:18) α + 1 α (cid:19) = 0 , (64)or equivalently τ = g ( β, σ ) (cid:44) σ − β β − /σ , (65) -1 0 1 2 3 4 5 fastslow fastslow E t = 020 c -1 0 1 2 3 4 5-20 -10 0 10fastslow fastslow ξ t = 0201203001000
250 500 750 0 0.1 0.2 0.3 0.4 0.5numfastslow t FIG. 7. Evolution of E component of the I Na -caricaturemodel with chosen sub- and super-threshold initial conditionin the comoving frame of reference. Parameters used: ∆ x =0 .
05, ∆ t = 4∆ x / α = 1, τ = 8 . x s = 1 . U s = 2 . U s = 2 . where σ = τ c , β = α/ ( α + 1) . (66)For the analytical expression of the first two left and righteigenfunctions, please see [2, 32].
2. Hybrid approach
Even though we know the ingredients of the linear the-ory for this model analytically, we still found them nu-merically as well. The hybrid approach is needed notonly because it helps to validate the analytical resultbut also because, in some cases, it is the only optionas the analytical derivation is not always possible. Dueto the discontinuous right-hand sides, it is essential touse the standard finite element method, at least whendealing with these discontinuous terms. The completediscretization formula for the critical front of the modelis presented in Appendix B.For two initial conditions, fig. 7 shows the evolutionof E component in the comoving frame of reference.For each case, the solution approaches the critical front, i.e. the solution at t = 120 in the figure and then givesrise to the stable propagating wave if the initial condi-tion is above the threshold or decays back to the restingstate otherwise. Fig. 8 gives the comparison of the nu-merical critical front obtained using operator splittingmethod and its analytical closed-form solution given by(62). We can see that the shooting procedure provides agood approximation of the critical front for the selectedparameters.Numerical experiments suggest that there are two val-ues of the speed c , c slow and c fast , satisfying 0 < c slow < -1 0 1 2 3 -30 -15 0 15 30analyticalnumerical E -30 -15 0 15 30 0 0.2 0.4 0.6 0.8 1analyticalnumerical ξh (a) (b) FIG. 8. (a) Comparison between analytical and numericalcritical front of the I Na -caricature model. For the numericalfront, we used following discretization parameters: α = 1, τ =8 .
2, ∆ x = 0 .
05, ∆ t = 4∆ x / L = 20. c fast < ∞ , such that the faster fronts are higher and sta-ble, and the slower fronts are lower and unstable, hencea slower front either dissipates or increases in the speedand magnitude to the fast branch solution depending onthe initial condition being below- and above-threshold,respectively [46, 47]. This can be seen in the right panelof fig. 7 where the blue circle and green square symbolsrepresent fast and slow front speeds for the selected pairof α = 1, τ = 8 .
2, and the red line indicates how thespeed of the front changes in time. For initial conditionsslightly above threshold, the front speed gets closer tothe slow speed and stays in the vicinity of it for a longtime before developing into the fast speed while the ini-tial condition slightly below threshold results in the frontspeed to drop to zero eventually.The next step after finding the critical front of the I Na -caricature model is the determination of the right and lefteigenfunctions along with the corresponding eigenvaluesemploying Algorithm 4 detailed in III A. In fig. 9, theeigenfunctions obtained using hybrid method fairly re-semble exact analytical eigenfunctions. The largest dif-ference between the numerical and analytical eigenfunc-tions is observed in the vicinity of the discontinuous val-ues, ξ = 0 and ξ = − ∆. This is totally expected asthe numerical scheme used to calculate the eigenfunc-tions is the second-order accurate, except near disconti-nuities, where it reduces to first order accuracy and in-troduces spurious oscillations due to the Beam-Warmingmethod [48].
3. Linear theory
The result of the calculation of the strength-durationcurve for I Na -caricature model is visualized in fig. 10,where the linear approximation is based on the formu-las (32) and (35). For every chosen stimulus duration, t s , we calculate the zeros of (35) in order to find thevalue of the shift s . We then substitute this value of s into one of the equations in (32) (both produce the sameresult) to get the corresponding value of I s . In the simu-lations, we choose two different set of the model param- -0.75-0.5-0.25 0 0.25 V V V V V V V V V E ana. h ana. E num. h num. -0.25 0 0.25 0.5 V V V ξ W -50 -25 0 25 -1-0.8-0.6-0.4-0.2 0 ξ W FIG. 9. Comparison between first two right and left eigen-functions of the I Na -caricature model. Parameters same asprevious figure. -3-2-1 0-30 -15 0 15 30 t s = 10 t s = 3 µ ( s ) -30 -15 0 15 30 t s = 10 t s = 3 s (a) (b) I ∗ s
0 1 2 3analyticalnumerical t s (c) (d) FIG. 10. Comparison of analytical and numerical strengthduration curve for I Na -caricature model for the choice of pa-rameters τ = 7 . α = 2 / τ = 8, α = 9 /
11 (b,d).Panels (a) and (b) show functions µ ( s ) defined by (35) for twoselected values of t s and their roots. Panels (c) and (d) arestrength-duration curves. eters, τ = 7 . α = 2 / τ = 8, α = 9 /
11 from whichthe resulting curves are respectively shown in fig. 10(a,c)and fig. 10(b,d). The shape of µ ( s ) is rather similar forboth cases and the main difference between the two isthe closeness of the s values for two different durationof stimulus values, t s = 3 and t s = 10. We observe thatthe theoretical critical curve for the first set of parametervalues is well adapted to the direct numerical simulationthreshold curve for smaller values of t s , and then bendsdown dramatically as the value of t s increases. The theo-retical prediction for the second set of parameters values,however, gets better with t s .6 c P . . . . β .
05 0 .
13c 0 . . λ . . λ ± · − ± · − (a) (b) FIG. 11. CV restitution curves for the FHN model for twoselected values of the model parameter.
B. FitzHugh-Nagumo system
1. Model formulation
The FitzHugh-Nagumo (FHN) system is a two-component reaction-diffusion system, which could beconsidered as a ZFK equation extended by adding a sec-ond, slow variable, describing inhibition of excitation. Itis probably the single historically most important modeldescribing excitable media. We consider it in the form d = 2 , D = diag(1 , , u = (cid:0) u, v (cid:1) (cid:62) , f ( u ) = (cid:0) f u ( u, v ) , f v ( u, v ) (cid:1) (cid:62) ,f u ( u, v ) = u ( u − β )(1 − u ) − v,f v ( u, v ) = γ ( αu − v ) . (67)for fixed values of the slow dynamics parameters, γ =0 .
01 and α = 0 .
37, and two values of the excitationthreshold for the fast dynamics, β = 0 .
05 and β = 0 .
2. Hybrid approach
System (67) has an unstable propagating pulse solutionas opposed to its reduced form, the ZFK equation withnontrivial stationary solution. It is known (see e.g. [34]and references therein) that in the limit γ (cid:38)
0, thecricial pulse solution whose v -component is small and u -component is close to the critical pulse of the corre-sponding ZFK equation. This makes it feasible to ob-tain explicit analytical solutions by using perturbationtechniques in the double limit γ (cid:38) β (cid:38)
0. Theseasymptotics will be described in a separate publication,and here we describe only the hybrid approach. The crit-ical pulse is obtained by applying Algorithm 3 by meansof AUTO. The corresponding CV restitution curves areillustrated in fig. 11. For the critical pulses, we take thesolutions at lower branches at
P > . · (see toprow of fig. 12). Other essential ingredients of the theoryare the first two left eigenfunctions and the first lead-ing eigenvalue, which are computed using Algorithm 4.The eigenfunctions for the two selected cases look rathersimilar, as shown in fig. 12. u ˆ u v ˆ u -0.1 0 0.1 0.2 W W -0.1 0 0.1 -40 -20 0 20 40 W -40 -20 0 20 40 ξ W FIG. 12. FHN theory ingredients for (a) α = 0 .
05 and (b) α = 0 .
13. Shown are components of scaled vector functions,indicated in top right corner of each panel, where ˆ u = S ˆ u , W j = S − W j , and S = diag(1 , ξ = 0 at the maximum of ˆ u . Correspondenceof lines with components is according to the legends at thetop.
3. Linear theory
Fig. 13 illustrates the calculation of the strength-duration curve for FHN model for α = 0 .
05 and α = 0 . s . Thecritical curves compared with those obtained from directnumerical simulation are sketched in fig. 13(c,d). Fromthis plot, it can be seen that the theoretical predictionfor both values of parameter is almost equally close tothe numerical prediction. C. The modified Beeler-Reuter model
1. Model formulation
Here, we looked at a variant of the classical Beeler-Reuter (BR) model of mammalian ventricular cardiacmyocytes [49], modified to describe phenomenologicallythe dynamics of neonatal rat cells [50–53]: d = 7 , D = diag(1 , , , , , , , (68) u = (cid:0) V, h, j, x , d, f, Ca (cid:1) (cid:62) , (69)7 -2-1 0 1 2 -50 -25 0 25 50 t s = 10 t s = 310 µ ( s ) -50 -25 0 25 50 t s = 10 t s = 3 s (a) (b) I ∗ s
0 1 2 3 4analyticalnumerical t s (c) (d) FIG. 13. Equation (35) that defines the shift s in termsof t s and comparison of analytical and numerical strengthduration-curve for the FHN model for β = 0 .
05 and β = 0 . γ = 0 . α = 0 .
37, ∆ x = 0 .
03, ∆ t =4∆ x / L = 60. f ( u ) = − ( I K + I x + I Na + I s ) α h (1 − h ) − β h hα j (1 − j ) − β j jα x (1 − x ) − β x x α d (1 − d ) − β d dα f (1 − f ) − β f f − − I s + 0 . − − Ca) . (70)For the detailed description of the components of f ( u ),please see the appendix of [2].
2. Hybrid approach
As in the FitzHugh-Nagumo system, the critical so-lution is a moving pulse, and thus, we obtain the CVrestitution curves and the critical pulse in a similar way.The CV restitution curves for the modified Beeler-Reutermodel is sketched the corresponding propagation speedsare shown in fig. 14. Apart from the critical pulse, thesolution at lower branches, the knowledge of the first twoleft eigenfunctions and the first leading eigenvalue arealso required. These ingredients have been found by themarching method given by Algorithm 4. The essentialingredients of the theory for two different data sets aresketched in fig. 15.
3. Linear theory
Fig. 16 exhibits the strength-duration threshold curveanalysis for BR model for two different excitability pa- c P . . . . α .
105 0 . . . λ . . λ ± · − ± · − (a) (b) FIG. 14. CV restitution curves for the BR model for twoselected values of the model parameter. Stable (upper) andunstable (lower) branches are shown by different line types. -0.5 0 0.5 1 ˆ uV h j x ˆ ud f Ca -1 0 1 2 3 W W -8-6-4-2 0 2 4 6 -50 -25 0 25 50 W -50 -25 0 25 50 ξ W FIG. 15. BR theory ingredients for (a) α = 0 .
105 and (b) α = 0 . u = S ˆ u , W j = 10 S − W j , and S = diag(10 − , , , , , , ). Thespace coordinate is chosen so that ξ = 0 at the maximum ofˆ V . Correspondence of lines with components is according tothe legends at the top. rameters, α = 0 .
105 and α = 0 . s are deter-mined by the transcendental equation (35) and comparedto FHN system, it is easier to detect the zeros of thisequation, two of which are shown in the top panel ofthe figure for t s = 3 and t s = 10. Then, the remainingpart is to insert the found value of s back into theoreticalthreshold curve generated by (32). The bottom panel ofthe figure shows these threshold curves being comparedwith numerical critical curves. As can be seen from thefigure, the analytical estimate for α = 0 .
115 provides asomewhat better approximation than that for α = 0 . -1 0 1 2 3 -50 -25 0 25 50 t s = 10 t s = 310 µ ( s ) -50 -25 0 25 50 t s = 10 t s = 3 s (a) (b)
10 100 0 3 6 9analyticalnumerical I ∗ s
0 3 6 9analyticalnumerical t s (c) (d) FIG. 16. Comparison of analytical and numerical strengthduration curve for BR for α = 0 .
105 (a) and α = 0 . x = 0 .
03, ∆ t = 4∆ x / L = 30,tolerance = 10 − . VI. DISCUSSION
The main aim of this thesis was to extend the methodproposed in [1] for analytical description of the thresholdcurves that separate the basins of attraction of propagat-ing wave solutions and of decaying solutions of certainreaction-diffusion models of spatially-extended excitablemedia. Specific aims are: • Extending the proposed theory to analysis of awider class of excitable systems, including mul-ticomponent reaction-diffusion systems, systemswith non-self-adjoint linearized operators and inparticular, systems with moving critical solutions(critical fronts and critical pulses). • Building an extension of this method from a linearto a quadratic approximation of the (center-)stablemanifold of the critical solution to demonstrate thediscrepancy between the analytical based on linearapproximation and numerical threshold curves en-countered when considering this quadratic approx-imation.The essential ingredients of the theory are the criti-cal solution itself, and the eigenfunctions of the corre-sponding linearized operator. For the linear approxima-tion in the critical nucleus case, we need the leading left(adjoint) eigenfunction; in the moving critical solutioncase, we need two leading left eigenfunctions; and for thequadratic approximations we require as many eigenvaluesand left and right eigenfunctions as possible to achievebetter accuracy. Of course, closed analytical formulasfor these ingredients can only be obtained in exceptional cases, and in a more typical situation a “hybrid” ap-proach is required, where these ingredients are obtainednumerically. We thus have provided insight into how thenumerical computation of these essential ingredients.The theory have been demonstrated on five differ-ent test problems ranging from one-component reaction-diffusion systems where the critical solution is the criticalnucleus to the multicomponent test problems with eithercritical front or critical pulse solution. In all models, theanalytical threshold curves are compared with the numer-ical simulations obtained using the bisection algorithm.We have applied both linear and quadratic approxima-tions for one-component test problems, ZFK and McK-ean models. The quadratic approximation agreed muchbetter with numerical threshold curves compared to thelinear approximation’s results, as would be expected.The accuracy and efficiency of the hybrid computationof the essential ingredients is dependent on the numeri-cal scheme and mesh resolution. It is obvious that someof the numerical schemes discussed in this paper do notoutperform other long-running and mathematically morecomplicated numerical schemes. In particular, the nu-merical study of I Na -caricature model has introduced thespurious oscillations and first-order accurate result nearthe discontinuities due to Beam-Warming method, whichcan be tackled using some advanced shape-preserving ad-vection schemes (see, for example, [54]). Hence, such ap-proaches that demand high computational cost can becarried out as an interesting direction for future researchif higher accuracy is required.As the results of the theory pointed out, our methodprovides more accurate results for some parameter val-ues than the others, especially in linear analysis. Eventhough the quadratic approximation offered for one-component test problems with the critical nucleus so-lutions provides more accurate estimates, still it is notfully understood why the choice of parameter values sig-nificantly matters. Hence, this remains quite importantresearch line of research. On the other hand, the pro-posed theory based on the moving critical solutions in-volves only the linear approximation. As an additionalconsideration, quadratic approximation for the cases ofmoving critical solutions can also be carried out.The theory established in this thesis is limited to onespatial dimension. Therefore, it could be of interest forfurther researches to adopt the theory to two and threedimensions.Throughout the theory, we have made the assumptionthat the spectrum is real. This is, however, not neces-sarily the case for the non-self-adjoint problems, whichremains an interesting direction for future research.Another extension of the work worth considering wouldbe to investigate the theory on some up-to-date realisticcardiac excitation models [55], simplified cardiac modelswith unusual properties [56, 57], and other excitable me-dia such as combustible media [58–60], pipe flow [61, 62] etc. APPENDICESAppendix A: Discretization formula forstrength-duration curve
With the aim of comparing the explicit approximationsfor the strength-duration curve, we describe how thethreshold curve is obtained using direct numerical simu-lations in this section. The numerical strength-durationthreshold curve is computed by solving the nonlinear sys-tem (1) for the initial and boundary conditions given by(4)-(5) using standard finite difference or finite elementdiscretization. More specifically, for ZFK, FHN and BRmodels we use finite differences, and for McKean and I Na -caricature models we implement finite element methodinstead.
1. Finite Difference Discretization Formula
The discretization formula for the generic form ofinitial-boundary value problem u t = Du xx + f ( u ) , u ( x,
0) = u r , Du x (0 , t ) = − I s H( t s − t ) e , x, t > , (A1)We employ explicit first order Euler forward differenceapproximation in time and explicit second order cen-tered finite difference approximation in space and plug-ging these into (A1), we obtainˆ u j +1 i = ˆ u ji + D ∆ t ∆ x (cid:16) ˆ u ji − − u ji + ˆ u ji +1 (cid:17) + ∆ t f (cid:16) ˆ u ji (cid:17) , (A2)in conjunction with its initial and boundary conditionsˆ u i = u r , ˆ u = ˆ u + 2∆ x I s H( t s − t ) e , ˆ u N +1 = ˆ u N − , (A3)ˆ u j +10 = ˆ u j +12 , ˆ u j +1 N +1 = ˆ u j +1 N − .
2. Finite Element Discretization Formula
Two of our test problems namely, McKean and I Na -caricature models, contain discontinuous kinetic terms. Discretization Formula for McKean equation
We use even extension of the problem and start withone component McKean model ∂u∂t = ∂ u∂x − u + H( u − a ) + 2 I s H( t s − t ) δ ( x ) ,u ( x,
0) = 0 , ( x, t ) ∈ R × R + . (A4) In the Galerkin finite element method, this can be writtenin the vector form, for ˇ u j , A dˇ u d t + ( A + B ) ˇ u = F + 2 I s H( t s − t ) D , (A5)where A is the the mass matrix, B is the stiffness matrixand F is the load vector (see [2] for a crude introduc-tion to finite element method and the derivation of thematrices). The vector D has only one nonzero entry bydefinition of the Dirac delta function.Finally, we employ the generalized trapezoidal rule(also known as θ scheme) [63], in which the residual isevaluated at j + θ , with this notation implyingˇ u n + θ = θ ˇ u n +1 + (1 − θ ) ˇ u n where 0 ≤ θ ≤ A + ∆ t θ ( A + B )] ˇ u n +1 = [ A − ∆ t (1 − θ ) ( A + B )] ˇ u n + F + 2 I s H( t s − t ) D . (A6)For θ = 0 , the linear system (A6) is the explicit Eulermethod that has a stability condition to be satisfied andits truncation error is O (∆ t ) + O (cid:0) ∆ x (cid:1) , θ = 1 / O (cid:0) ∆ t (cid:1) + O (cid:0) ∆ x (cid:1) , and θ = 1 gives the first-order accurate implicit Euler rulethat is also unconditionally stable and its truncation er-ror is O (∆ t ) + O (cid:0) ∆ x (cid:1) [63]. Discretization Formula for the I Na -caricature model Even extended version of the model is in the followingform ∂E∂t = ∂ E∂x + H( E − h + 2 I s H( t s − t ) δ ( x ) ,∂h∂t = 1 τ (H( − E ) − h ) , (A7) E ( x,
0) = − α, h ( x,
0) = 1 , ( x, t ) ∈ R × R + . The fully discretized finite element discretization formulafor the model is[ A + ∆ t θ B ] ˇ E n +1 = [ A − ∆ t (1 − θ ) B ] ˇ E n + ∆ t S ˇ h n + 2∆ t I s H( t s − t ) D [ τ + ∆ t θ ] A ˇ h n +1 = [ τ − ∆ t (1 − θ )] A ˇ h n + ∆ t G, (A8)where S = [ s i,j ] = (cid:90) L − L H (cid:0) ˇ E − (cid:1) Φ i ( ξ )Φ j ( ξ ) d ξ,G = [ g i ] = (cid:90) L − L Φ i ( ξ )H (cid:0) − ˇ E (cid:1) d ξ. S is a tridiagonal matrix with following di-agonal elements: s i,i = (cid:90) L − L H (cid:0) ˇ E − (cid:1) Φ i ( ξ ) d ξ = 1∆ ξ (cid:32)(cid:90) ξ i ξ i − H (cid:0) ˇ E − (cid:1) ( ξ − ξ i − ) d ξ + (cid:90) ξ i +1 ξ i H (cid:0) ˇ E − (cid:1) ( ξ i +1 − ξ ) d ξ (cid:33) = I i + I i , for i = 2 , , · · · , N − I i = 13∆ ξ ∆ ξ , E t i − , (cid:0) ˇ ξ i − − ξ i − (cid:1) , E m i − , ∆ ξ − (cid:0) ˇ ξ i − − ξ i − (cid:1) , E b i − , , otherwise, I i = 13∆ ξ ∆ ξ , E t i , ∆ ξ − (cid:0) ξ i +1 − ˇ ξ i +1 (cid:1) , E m i , (cid:0) ξ i +1 − ˇ ξ i +1 (cid:1) , E b i , , otherwise. s i,i +1 = (cid:90) L − L H (cid:0) ˇ E − (cid:1) Φ i ( ξ )Φ i +1 ( ξ )d ξ = 1∆ ξ (cid:90) ξ i +1 ξ i H (cid:0) ˇ E − (cid:1) ( ξ i +1 − ξ ) ( ξ − ξ i ) d ξ = 16∆ ξ ∆ ξ , E t i , ξ i +1 ˇ ξ i +1 + 3 ξ i +1 ξ i − ξ i +1 − ξ i , E m i , − ξ i +1 ξ i ˇ ξ i +1 + 3 ξ i ˇ ξ i +1 ξ i +1 − ξ i +1 ˇ ξ i +1 + 2 ˇ ξ i +1 − ξ i ˇ ξ i +1 , E b i , − ξ i ˇ ξ i +1 + 6 ξ i ξ i +1 ˇ ξ i +1 , otherwise. s i − ,i = (cid:90) L − L H (cid:0) ˇ E − (cid:1) Φ i ( ξ )Φ i − ( ξ )d ξ = 1∆ ξ (cid:90) ξ i ξ − H (cid:0) ˇ E − (cid:1) ( ξ − ξ i − ) ( ξ i − ξ ) d ξ = 16∆ ξ ∆ ξ , E t i − , ξ i ˇ ξ i − + 3 ξ i − ˇ ξ i − + 3 ξ i ξ i − − ξ i − , E m i − , − ξ i − ξ i ˇ ξ i − − ξ i − ξ i + 6 ξ i ξ i − ˇ ξ i − + 2 ˇ ξ i − − ξ i − ξ i , E b i − , − ξ i ˇ ξ i − − ξ i − ˇ ξ i − , otherwise.We need to implement no flux boundary condition s , = 13∆ ξ ξ , E t , ξ − (cid:0) ξ − ˇ ξ (cid:1) − (cid:0) ˇ ξ − ξ + 2∆ ξ (cid:1) , E m , (cid:0) ˇ ξ − ξ + 2∆ ξ (cid:1) + (cid:0) ξ − ˇ ξ (cid:1) , E b , , otherwise. s N,N = 13∆ ξ ξ , E t N − , (cid:0) ξ N − + 2∆ ξ − ˇ ξ N − (cid:1) + (cid:0) ˇ ξ N − − ξ N − (cid:1) , E m N − , − (cid:0) ξ N − + 2∆ ξ − ˇ ξ N − (cid:1) ξ − (cid:0) ˇ ξ N − − ξ N − (cid:1) , E b N − , , otherwise.In these formulas, we use a shorthand notation, E t o E m o E b o = ˇ E o > , ˇ E o +1 > , ˇ E o > , ˇ E o +1 < , ˇ E o < , ˇ E o +1 > . Having in mind that the tent functions we have chosenare piecewise linear, the points ˇ ξ i − , ˇ ξ i +1 , ˇ ξ and ˇ ξ N − are then obtained from the linear interpolation methodˇ ξ p +1 = (cid:2) ˇ E ( ξ p +1 ) − (cid:3) ξ p + (cid:2) − ˇ E ( ξ p ) (cid:3) ξ p +1 ˇ E ( ξ p +1 ) − ˇ E ( ξ p ) , for p = 1 , i − , i, N −
2. The vector G has entries for i = 2 , , · · · , N − g i = 1∆ ξ (cid:32)(cid:90) ξ i ξ i − H (cid:0) − ˇ E (cid:1) ( ξ − ξ i − ) d ξ + (cid:90) ξ i +1 ξ i H (cid:0) − ˇ E (cid:1) ( ξ i +1 − ξ ) d ξ = I i + I i (cid:33) , I i = 12∆ ξ ∆ ξ , E t i − , (cid:0) ξ i − − ξ i − (cid:1) , E m i − , ∆ ξ − (cid:0) ξ i − − ξ i − (cid:1) , E b i − , , otherwise.and I i = 12∆ ξ ∆ ξ , E t i , ∆ ξ − (cid:0) ξ i +1 − ξ i +1 (cid:1) , E m i , (cid:0) ξ i +1 − ξ i +1 (cid:1) , E b i , , otherwise.and on the boundaries g = 12∆ ξ x , E t , ξ − (cid:0) ξ − ξ + 2∆ ξ (cid:1) − (cid:0) ξ − ξ (cid:1) , E m , (cid:0) ξ − ξ + 2∆ ξ (cid:1) + (cid:0) ξ − ξ (cid:1) , E b , , otherwise. g N = 12∆ ξ x , E t N − , (cid:0) ξ N − + 2∆ ξ − ξ N − (cid:1) + (cid:0) ξ N − − ξ N − (cid:1) , E m N − , ξ − (cid:0) ξ N − − ξ N − (cid:1) − (cid:0) ξ N − + 2∆ ξ − ξ N − (cid:1) , E b N − , , otherwise.The shorthand notations used above are, E t o E m o E b o = ˇ E o < , ˇ E o +1 < , ˇ E o < , ˇ E o +1 > , ˇ E o > , ˇ E o +1 < , and ξ are found via interpolation method ξ p +1 = ˇ E ( ξ p +1 ) ξ p − ˇ E ( ξ p ) ξ p +1 ˇ E ( ξ p +1 ) − ˇ E ( ξ p ) , for p = 1 , i − , i, N −
3. Threshold curve
To obtain the threshold curve in the stimulus strength-duration plane, we solve a sequence of the “stimulation bycurrent” initial-value problem (12) and (13). The choiceof the numerical scheme changes according to the spe-cific model as defined above. The computation is doneby fixing stimulation time duration t s and varying thestrength of the current I s . For any initial upper estimate I s (superthreshold), known to be sufficient for ignition,and lower estimate I s , known to fail to ignite, the follow-ing bisection algorithm gives the threshold value I ∗ s : Input: t s , I s , I s Output: I ∗ s while (cid:0) | I s − I s | ≥ tolerance (cid:1) do • I s ← (cid:0) I s + I s (cid:1) / • Solve (12)-(13) with I s = I s if ignition then • I s ← I s else • I s ← I s end ifend while • I ∗ s ← I s Algorithm 5:
Bisection loop for finding the strengthof the current I ∗ s for a fixed parameter t s .This procedure is repeated as many times for different t s as necessary to obtain the strength-duration curve. Appendix B: Numerical methods for simplifiedcardiac excitation model
This appendix contains the discretization formula forthe I Na -caricature model. Specifically, we aim to providethe numerical procedure for finding the critical front andfirst two leading eigenvalues and the corresponding leftand right eigenfunctions accordingly. To begin with, weintroduce co-moving frame of reference by setting ξ = x − x f ( t ) and T = t such that E ( ξ, T ) is the voltage profilein the standard position and x f ( t ) is the movement ofthis profile. Problem (A7) then becomes ∂E∂T = ∂ E∂ξ + x f (cid:48) ( t ) ∂E∂ξ + H ( E − h,∂h∂T = x f (cid:48) ( t ) ∂h∂ξ + 1 τ (H ( − E ) − h ) . (B1)Here, we do the numerical computation through the ini-tial condition E ( ξ,
0) = E s H ( − ξ ) H ( ξ + 2 x s ) − α. (B2)We also need to impose a pinning condition in order tofind the value of x f (cid:48) which varies at each time step. Thiscan be achieved by considering the shape of E component2of the critical front solution. A common way to definesuch condition is to choose a constant E ∗ representedonce in the front profile for every time step, E ( x f ( t ) , T ) = E ∗ . (B3)For simplicity, we take E ∗ = 0.
1. Discretization formula for the critical front ofthe model
Numerical analysis of the critical front for I Na -caricature model is based on a combination of finite ele-ment and finite difference methods by using the operatorsplitting method (see e.g. [64]). Finite element method isused to handle the right hand sides of the equations withdiscontinuity terms involving the Heaviside step function.For the standard finite difference discretization, we useBeam-Warming scheme for the first spatial derivatives ofboth E and h . We set the domain of ξ and T coordinatesto be − L ≤ ξ ≤ L , 0 ≤ T ≤ T f so that the grid points ξ i , T j are ξ i = − L + i ∆ ξ , T j = j ∆ T , where ∆ ξ > T > i = 0 , , · · · , N , j = 0 , , · · · , M for N, M > E ji , h ji as the numerical approximation of E ( ξ i , T j )and h ( ξ i , T j ), E ji ∗ = 0 as our pinning condition whichwill be further explained later, and finally c j = x f (cid:48) ( T j )as the speed at j − th time step. Using these representa-tions, we solve (B1) numerically using operator splittingmethod approach in the following seven steps: Step 1(Finite element method):
As a first step, wesolve E equation without the advection term as it is mul-tiplied by the speed which is not determined yet ∂E∂T = ∂ E∂ξ + H ( E − h. We have to employ the finite element method due to thediscontinuous Heaviside step function which gives( A + ∆ T B θ ) E j + i = ( A − ∆ T B (1 − θ )) E ji + ∆ T S h ji . (B4) Step 2(Thomas algorithm):
Inserting the elements ofthe matrices A , B , D and the discretized solution E ji and h ji (both known) into (B4) yields a tridiagonal system of N equations in a following form β γ . . . α β γ . . . α β . . . ...... ... . . . . . . γ N − . . . α N β N E j + ......... E j + N = δ ......... δ N . This can be solved using a standard Gaussian eliminationmethod such as Thomas algorithm and using such algo-rithm is sometimes crucial as it leads to a reduced com-putational cost. The back substitution procedure (see e.g. [65] for more detailed explanation) generates the so-lution as γ (cid:48) i = (cid:40) γ i β i , i = 1 , γ i β i − α i γ (cid:48) i − , i = 2 , , · · · , N − ,δ (cid:48) i = (cid:40) δ i β i , i = 1 , δ i − α i δ (cid:48) i − β i − α i − γ (cid:48) i − , i = 2 , , · · · , N,E j + N = δ (cid:48) N , (B5) E j + i = δ (cid:48) i − γ (cid:48) i E j + i +1 , i = N − , N − , · · · , . Step 3(Finding the value of the speed):
We havedivided E equation in (B1) into two parts and it remainsto find the solution of the advection step in the splittedscheme. Before we update the solution, it is necessaryto find the value of the speed according to the pinningcondition E j + i ∗ = 0, where the index i ∗ corresponds toan integer value, indicating the spatial position at whichthe solution is equal to zero initially, i.e. ξ i ∗ = 0. Asthe Beam-Warming method is second order accurate, wefind the speed value using the Beam-Warming scheme bymeans of the discretized solution found in the previousstep as, c j = − (cid:114)(cid:104) ∆ T ξ (cid:16) E j + i ∗ − E j + i ∗ − + E j + i ∗ − (cid:17)(cid:105) − T ∆ ξ (cid:16) E j + i ∗ − E j + i ∗ − + E j + i ∗ − (cid:17) (cid:16) E j + i ∗ − E ∗ (cid:17) ∆ T ∆ ξ (cid:16) E j + i ∗ − E j + i ∗ − + E j + i ∗ − (cid:17) − ∆ ξ (cid:16) E j + i ∗ − E j + i ∗ − + E j + i ∗ − (cid:17) T (cid:16) E j + i ∗ − E j + i ∗ − + E j + i ∗ − (cid:17) . (B6)3This formula is actually derived from the following stan-dard Beam-Warming discretization, which is quadraticin c . Step 4(Beam-Warming scheme):
The next step is touse Beam-Warming scheme to update the solution of ad-vection term of E component at step j + 1, E j +1 i = E j + i + c j ∆ T ξ (cid:16) E j + i − E j + i − + E j + i − (cid:17) + (cid:32) c j ∆ T √ ξ (cid:33) (cid:16) E j + i − E j + i − + E j + i − (cid:17) . (B7) Step 5(Finite element method):
In a similar manner,the equation of h component can be divided into two partand once again, we use finite element method to solve h equation with the advection term removed ∂h∂T = 1 τ (H ( − E ) − h ) , which gives A (cid:20) T θτ (cid:21) h j + i = A (cid:20) − ∆ T (1 − θ ) τ (cid:21) h ji + ∆ T Gτ . (B8)
Step 6 (Thomas algorithm) : As (B8) is also a tridi-agonal matrix, we can employ the Thomas algorithm hereas well, similarly to the case of E equation. Step 7( Beam-Warming scheme):
Once again, weuse second order accurate Beam-Warming scheme withtruncation error O (cid:0) ∆ T , ∆ ξ (cid:1) for the advection term ofthe h component h j +1 i = h j + i + c j ∆ T ξ (cid:16) h j + i − h j + i − + h j + i − (cid:17) + (cid:32) c j ∆ T √ ξ (cid:33) (cid:16) h j + i − h j + i − + h j + i − (cid:17) . (B9)The numerical computation of the critical front isachieved by solving (B1) according to above seven-step procedure. As the I Na -caricature model is a two-component system, to find the numerical estimation ofthe critical front, we calculate S ( T ) as S ( T ) = L (cid:90) − L (cid:0) E T ( ξ, T ) + h T ( ξ, T ) (cid:1) d ξ. (B10)
2. Discretization Formula for the LinearizedProblem
We linearize the system (B1) about the critical front (cid:16) ˆ E, ˆ h (cid:17) using E ( ξ, T ) = ˆ E ( ξ ) + (cid:15)E ( ξ, T ) ,h ( ξ, T ) = ˆ h ( ξ ) + (cid:15)h ( ξ, T ) , (B11) where (cid:15) (cid:28) , | E ( ξ, T ) | (cid:28) , | h ( ξ, T ) | (cid:28)
1. Hence, wehave the following system of equations: ∂E∂T = ∂ E∂ξ + c ∂E∂ξ − E (cid:48) ( − ∆) δ ( ξ + ∆) ˆ hE + H ( − ∆ − ξ ) h,∂h∂T = c ∂h∂ξ + (cid:32) E (cid:48) (0) δ ( ξ ) E − h (cid:33) /τ. (B12)We solve this linearized equation with the operator split-ting technique, by splitting the system into four equa-tions. We use either the finite element or finite differencemethods to obtain the solution of each of these four equa-tions as follows: Step 1 (Finite element method):
First of all, wesolve ∂E∂T = ∂ E∂ξ − E (cid:48) ( − ∆) δ ( ξ + ∆) ˆ hE + H ( − ∆ − ξ ) h, using the finite element method this yields (cid:34) A + ∆ T θ (cid:32) B + L ˆ E (cid:48) ( − ∆) (cid:33)(cid:35) E j + i (B13)= (cid:34) A − ∆ T (1 − θ ) (cid:32) B + L ˆ E (cid:48) ( − ∆) (cid:33)(cid:35) E ji + ∆ T K h ji , where the matrices K and L are K = [ k i,j ] = (cid:90) L − L H ( − ∆ − ξ ) Φ i ( ξ )Φ j ( ξ ) d ξ, L = [ l i,j ] = (cid:90) L − L δ ( ξ + ∆) ˆ h ( ξ )Φ i ( x )Φ j ( ξ ) d ξ = ˆ h ( − ∆) Φ i ( − ∆) Φ j ( − ∆) . The matrix L has exactly 4 non-zero elements and theseare l i,j = − ˆ h ( − ∆)∆ ξ − ( ξ l +1 + ∆) , i = j = l, ( ξ l +1 + ∆) (∆ + ξ l ) , i = l, j = l + 1 , (∆ + ξ l ) ( ξ l +1 + ∆) , i = l + 1 , j = l, − (∆ + ξ l ) , i = j = l + 1 , , otherwise , where ξ l ≤ ∆ ≤ ξ l +1 . Actually, this can be further sim-plified by discretizating the spatial domain in the waythat ∆ is situated exactly on the grid that makes L withonly one non-zero element. On the other hand, the diag-onal entries of the matrix K are found as k i,i = (cid:90) L − L H ( − ∆ − ξ ) Φ i ( ξ ) d ξ = I i + I i , where I i = 13∆ ξ ∆ ξ , E t − ∆ i − , − (∆ + ξ i − ) , E m − ∆ i − , ∆ ξ + (∆ + ξ i − ) , E b − ∆ i − , , otherwise,4and I i = 13∆ ξ ∆ ξ , E t − ∆ i , ∆ ξ − ( ξ i +1 + ∆) , E m − ∆ i , ( ξ i +1 + ∆) , E b − ∆ i , , otherwise.The supradiagonal elements of K are evaluated as k i,i +1 = (cid:90) L − L H ( − ∆ − ξ ) Φ i ( ξ )Φ i +1 ( ξ )d ξ = 1∆ ξ (cid:90) ξ i +1 ξ i H ( − ∆ − ξ ) ( ξ i +1 − ξ ) ( ξ − ξ i ) d ξ = 16∆ ξ ∆ ξ , E t − ∆ i , ξ i +1 ∆ + 6 ξ i +1 ξ i ∆+ 2∆ + 3 ξ i ∆ , E m − ∆ i , + 3 ξ i +1 ξ i − ∆ ,ξ i +1 − ξ i ξ i +1 − − ξ i ∆ , E b − ∆ i , − ξ i +1 ∆ − ξ i ξ i +1 ∆0 , otherwise.The subdiagonal elements of K are also found as k i − ,i = (cid:90) L − L H ( − ∆ − ξ ) Φ i ( ξ )Φ i − ( ξ )d ξ = 1∆ ξ (cid:90) ξ i ξ i − H ( − ∆ − ξ ) ( ξ − ξ i − ) ( ξ i − ξ ) d ξ = 16∆ ξ ∆ ξ , E t − ∆ i − , ξ i ∆ + 2∆ + 3 ξ i ξ i − − ξ i − , E m − ∆ i − , + 6 ξ i − ξ i ∆ + 3 ξ i − ∆ ξ i − ξ i − ξ i − ξ i ξ i − ∆ − , E b − ∆ i − , − ξ i ∆ − ξ i − ∆ , otherwiseand on the boundaries, we have k , = 13∆ ξ ξ , E t − ∆ , ξ − (2∆ ξ − ∆ − ξ ) − ( ξ + ∆) , E m − ∆ , (2∆ ξ − ∆ − ξ ) + ∆ ξ , E b − ∆ , , otherwise, k N,N = 13∆ ξ ξ , E t − ∆ N − , ( ξ N − + 2∆ ξ + ∆) − (∆ + ξ N − ) , E m − ∆ N − , ξ + (∆ − ξ N − ) − ( ξ N − + 2∆ ξ + ∆) , E b − ∆ N − , , otherwise.The shorthand notations used above are, E t − ∆ o E m − ∆ o E b − ∆ o = ˇ E o < − ∆ , ˇ E o +1 < − ∆ , ˇ E o < − ∆ , ˇ E o +1 > − ∆ , ˇ E o > − ∆ , ˇ E o +1 < − ∆ . Step 2 (Beam Warming scheme):
The secondstep is to solve the pure advection equation using theBeam-Warming scheme giving E j +1 i = E j + i + c ∆ T ξ (cid:16) E j + i +1 − E j + i − E j + i +2 (cid:17) + (cid:32) c ∆ T √ ξ (cid:33) (cid:16) E j + i − E j + i +1 + E j + i +2 (cid:17) , (B14)that finishes the numerical scheme of E component of thelinearized problem. Step 3 (Finite element method):
Again the finite el-ement method is conveniently employed to numericallysolve the first equation of h component, ∂h∂T = (cid:32) E (cid:48) (0) δ ( ξ ) E − h (cid:33) /τ, that results in (cid:20) A + ∆ T θ A τ (cid:21) h j + i = (cid:20) A − ∆ T (1 − θ ) A τ (cid:21) h ji + ∆ T M τ ˆ E (cid:48) (0) E ji , (B15)where M = [ m i,j ] = (cid:90) L − L δ ( ξ )Φ i ( ξ )Φ j ( ξ ) d ξ = Φ i (0)Φ j (0)= 1∆ ξ ξ r +1 , i = j = r, − ξ r +1 ξ r , i = r, j = r + 1 , − ξ r ξ r +1 , i = r + 1 , j = r,ξ r , i = j = r + 1 , , otherwise,such that ξ r ≤ ≤ ξ r +1 . Step 4 (Beam Warming scheme):
Similar to E h component, h j +1 i = h j + i + c ∆ T ξ (cid:16) h j + i +1 − h j + i − h j + i +2 (cid:17) + (cid:32) c ∆ T √ ξ (cid:33) (cid:16) h j + i − h j + i +1 + h j + i +2 (cid:17) , that completes the numerical solution of the linearizedequation.
3. Discretization Formula for the AdjointLinearized Problem
The adjoint linearized problem for the I Na -caricaturemodel is ∂ ¯ E∂T = ∂ ¯ E∂ξ − c ∂ ¯ E∂ξ − E (cid:48) ( − ∆) δ ( ξ + ∆) ˆ h ¯ E + δ ( ξ ) τ ˆ E (cid:48) (0) ¯ h,∂ ¯ h∂T = − c ∂ ¯ h∂ξ + H ( − ∆ − ξ ) ¯ E − ¯ hτ . (B16)This problem can also be solved in 4-steps as follows: Step 1 (Finite element method):
We begin with theequation of ¯ E component and solve first the following, ∂ ¯ E∂T = ∂ ¯ E∂ξ − E (cid:48) ( − ∆) δ ( ξ + ∆) ˆ h ¯ E + δ ( ξ ) τ ˆ E (cid:48) (0) ¯ h with the solution based on the finite element method, (cid:34) A + ∆ T θ (cid:32) B + L ˆ E (cid:48) ( − ∆) (cid:33)(cid:35) ¯ E j + i = ∆ T M τ ˆ E (cid:48) (0) ¯ h ji + (cid:34) A − ∆ T (1 − θ ) (cid:32) B + L ˆ E (cid:48) ( − ∆) (cid:33)(cid:35) ¯ E ji . (B17) Step 2 (Beam Warming scheme):
As the advectionterm has negative sign in the front, the Beam-Warmingnumerical scheme is in this case,¯ E j +1 i = ¯ E j + i − c ∆ T ξ (cid:16) E j + i − E j + i − + ¯ E j + i − (cid:17) + (cid:32) c ∆ T √ ξ (cid:33) (cid:16) ¯ E j + i − E j + i − + ¯ E j + i − (cid:17) . Step 3 (Finite element method):
Using the finite el-ement method to solve, ∂ ¯ h∂T = H ( − ξ − ∆) ¯ E − ¯ hτ , let us obtain, (cid:20) A + A ∆ T θτ (cid:21) ¯ h j + i = (cid:20) A − ∆ T A (1 − θ ) τ (cid:21) ¯ h ji +∆ T M ¯ E ji . (B18) Step 4 (Beam Warming scheme):
We conclude thenumerical solution of the adjoint problem by solving thepure advection equation for ¯ h component,¯ h j +1 i =¯ h j + i − c ∆ T ξ (cid:16) h j + i − h j + i − + ¯ h j + i − (cid:17) + (cid:32) c ∆ T √ ξ (cid:33) (cid:16) ¯ h j + i − h j + i − + ¯ h j + i − (cid:17) . We note that the Thomas algorithm can also be appliedto the diagonal systems (B13), (B15), (B17) and (B18)in a similar manner to the case of critical front steps.Hence, we skip the similar derivations here.For this model, the linear approximation of the criticalcurves requires the knowledge of the critical front as wellas first two leading eigenvalues and corresponding leftand right eigenfunctions. The numerical calculating ofthe eigenpairs are determined by the method discussed inSection III, and these two last subsections are dedicatedto how to derive the first two eigenmodes as a numericalsolution of the linearized and adjoint linearized problemsby means of the operator splitting method. Alternatively,the implicitly restarted Arnoldi method [66] can be usedto estimate these essential ingredients, in which case weuse the matrix representations of the discretized versionsof the equations (B12) and (B16). [1] I. Idris and V. N. Biktashev, Physical Review Letters , 244101 (2008).[2] B. Bezekci, I. Idris, R. D. Simitev, and V. N. Biktashev, Physical Review E , 042917 (2015).[3] M. B. Schiffer, Advances in archaeological method andtheory (New York: Academic Press, 1985). [4] R. FitzHugh, The Bulletin of Mathematical Biophysics , 257 (1955).[5] D. Zipes and J. Jalife, Cardiac electrophysiology: fromcell to bedside (WB Saunders CO, 2000).[6] M. Roth, British Medical Bulletin , 42 (1986).[7] T. Arendt and V. Bigl, Neurobiology of Aging , 552(1987).[8] V. S. Zykov, Scholarpedia , 1834 (2008).[9] D. Barkley, Physica D: Nonlinear Phenomena , 61(1991).[10] Y. Guo, Y. Zhao, S. A. Billings, D. Coca, R. I. Ristic,and L. DeMatos, International Journal of Bifurcation andChaos , 2137 (2010).[11] M. C. Cross and P. C. Hohenberg, Reviews of ModernPhysics , 851 (1993).[12] W. D. McCormick, Z. Noszticzius, and H. L. Swinney,The Journal of Chemical Physics , 2159 (1991).[13] D. T. Kaplan, J. R. Clay, T. Manning, L. Glass, M. R.Guevara, and A. Shrier, Physical Review Letters ,4074 (1996).[14] I. Farkas, D. Helbing, and T. Vicsek, Nature , 131(2002).[15] G. Seiden and S. Curland, New Journal of Physics ,033049 (2015).[16] F. H. Fenton, E. M. Cherry, H. M. Hastings, and S. J.Evans, Biosystems , 73 (2002).[17] G. Weiss, Archives Italiennes de Biologie , 413 (1990).[18] H. Bostock, The Journal of Physiology , 59 (1983).[19] L. Lapicque, J. Physiol. Pathol. Gen , 620 (1907).[20] N. Brunel and M. C. W. van Rossum, Biological Cyber-netics , 341 (2007).[21] H. A. Blair, The Journal of General Physiology , 709(1932).[22] H. A. Blair, The Journal of General Physiology , 731(1932).[23] N. Rashevsky, Protoplasma , 42 (1933).[24] A. M. Monnier, Hermann (1934).[25] A. V. Hill, Proceedings of the Royal Society of London.Series B, Biological Sciences , 305 (1936).[26] W. Nernst, Pfl¨ugers Archiv European Journal of Physi-ology , 275 (1908).[27] W. A. H. Rushton, Proceedings of the Royal Society ofLondon. Series B, Biological Sciences , 210 (1937).[28] D. Noble and R. B. Stein, J Physiol , 129 (1966).[29] H. P. McKean and V. Moll, Bulletin of the AmericanMathematical Society , 255 (1985).[30] G. Flores, Journal of Differential Equations , 306(1989).[31] J. C. Neu, R. S. Preissig, and W. Krassowska, PhysicaD: Nonlinear Phenomena , 285 (1997).[32] I. Idris, Initiation Of Excitation Waves , Ph.D. thesis,University of Liverpool (2008).[33] I. Idris and V. N. Biktashev, Physical Review E ,021906 (2007).[34] G. Flores, SIAM Journal on Mathematical Analysis ,392 (1991).[35] This assumption will, of course, have to be verified ineach particular case.[36] O. A. Mornev, “On the conditions of excitation of one-dimensional autowave media,” in Autowave processes insystems with diffusion (Institute of Applied Physics ofthe USSR Academy of Sciences, Gorky, 1981) pp. 92–98.[37] H. P. McKean, Advances in Mathematics , 209 (1970). [38] V. Moll and S. I. Rosencrans, SIAM Journal on AppliedMathematics , 1419 (1990).[39] R. D. Simitev and V. N. Biktashev, Bulletin of Mathe-matical Biology , 72 (2011).[40] E. Doedel and J. P. Kernevez, AUTO, software for con-tinuation and bifurcation problems in ordinary differen-tial equations (California Institute of Technology, 1986).[41] Y. B. Zel’dovich and D. A. Frank-Kamenetsky, in