Strength-Duration Relationship in an Excitable Medium

We consider the strength-duration relationship in one-dimensional spatially extended excitable media. In a previous study [Idris and Biktashev 2008] set out to separate initial (or boundary) conditions leading to propagation wave solutions from those leading to decay solutions, an analytical criterion based on an approximation of the (center-)stable manifold of a certain critical solution was presented. The theoretical prediction in the case of strength-extent curve was later on extended to cover a wider class of excitable systems including multicomponent reaction-diffusion systems, systems with non-self-adjoint linearized operators and in particular, systems with moving critical solutions (critical fronts and critical pulses) [Bezekci et al. 2015]. In the present work, we consider extension of the theory to the case of strength-duration curve.


Motivation
The threshold phenomenon "deals with the minimal, an event, or stimulus just strong enough to be perceived or to produce a response" [56] and the presence of it "imposes the restriction on the types of mathematical model suitable to describe" biological/chemical systems [25]. Its extreme importance can be highlighted through examples. For instance, propagation of excitation in the heart involves action potential and threshold value controls if an applied stimulus is sufficient enough to generate an action potential. Understanding the mechanisms of initiation of propagating is extremely crucial as successful propagation enables continuous electrical and chemical communication between cells and failure may lead to serious medical conditions [64]. Threshold phenomenon also plays a key role in understanding many age related diseases such as Alzheimer and Parkinson. Studies on neuronal changes in brain suggest that the threshold hypothesis helps to explain "some of the associations between clinical and pathological findings" [53,1].
Originally, the term excitability has come to be used to refer to the "property of living organisms to respond strongly to the action of a relatively weak external stimulus" [65]. A well-known example of excitability is the ability of nerve cells to generate and propagate electrical activity. By definition, an excitable medium is a spatially distributed system, each element of which possesses the property of excitability and it is usually defined as nonlinear reaction-diffusion system, where the reaction term defines how the constituents of the system are transformed into each other, and the diffusion part provides propagation of information [3,29,65]. There are a wide variety of areas where the term "excitable medium" has been used repeatably for decades in many fields including physical, chemical and biological systems and so on [18,40,37,65,22,59]. * Corresponding author burhanbezekci@kilis.edu.tr (B. Bezekci); v.n.biktashev@exeter.ac.uk (V.N. Biktashev) ORCID(s):

Problem statement
We consider the problems of initiation of propagating waves in one-dimensional reaction-diffusion systems, where ∶ ℝ × ℝ → ℝ is a -component reagents field, ≥ 1, defined for ∈ ℝ and ∈ ℝ + , vector-function ∶ ℝ → ℝ describes the reaction rates and ∈ ℝ × is the matrix of diffusivity. Equation (1) is assumed to describe an excitable medium as a system "composed of elementary segments or cells, each of which possesses the following properties: 1. a well-defined rest state, 2. a threshold for excitation, and 3. a diffusive-type coupling to its nearest neighbors. . . . Stimuli below the threshold are damped out and produce no persistent change in the system, . . . stimuli above the threshold induce the cell to change from its rest state to an excited state." [24] A closely related class are bistable systems: whereas an excitable system proper returns to the resting state after spending some time in the excitable state, a bistable system remains in the excitable state for ever.
A definitive feature of an excitable or bistable medium is existence of traveling wave solutions of (1). These can be described by transforming the system of partial differential equations (PDEs) to a moving frame of reference, We are interested in the solutions that are stationary in this frame of reference, for a fixed , i.e.
If the velocity = 0, then the traveling wave is called the standing wave. The traveling wave is a front if (−∞) and (∞) exist and different from each other (this is typical for bistable systems), and it is a pulse if (∞) = (−∞) = (this happens in excitable systems).
Travelling wave solutions of (1) have been a topic of intense study. For applications, for instance modelling of biological media and chemical processes, the question of particular importance is emergence of such solutions as a result of a perturbation of the resting state, localized in space where and describe the shapes of the initial and boundary profiles, and s and s are the strengths of those profiles. The cases of non-homogeneous initial condition and non-homogeneous boundary condition are usually handled separately. In electrophysiological terms, these can be described as follows: 1. Stimulation by current: s = 0, s ≠ 0. This is the case when the current is injected at the boundary point = 0 during some time interval. For a fixed boundary profile ( ), there exist a corresponding threshold strength value * s such that the solution tends to propagating wave ("ignition") as → ∞ whenever s > * s , and the solution tends to resting state ("failure") otherwise. For a one-parametric family of profiles, parametrized by the stimulus duration s , the corresponding curve * s ( s ) is called a strength-duration curve (see fig. 1). 2. Stimulation by voltage: s ≠ 0, s = 0. Here the perturbation is instantaneous at = 0, but is spread in space. For a fixed initial profile ( ), there exist a corresponding threshold strength value * s such that the solution tends to propagating wave as → ∞ whenever s > * s , and to resting state otherwise. For a one-parametric family of profiles, parametrized by the stimulus extent s , we shall have the corresponding critical curve * s ( s ), called a strength-extent curve. In our previous paper [7] we have analysed some analytical and semi-analytical approaches to description of the strength-extent curves. In this paper, we focus on the strength-duration curves. In all specific examples below we shall consider a rectangular profile of duration s , where the fixed vector determines which reagents are being injected, and H(⋅) is the Heaviside step function.

A brief history of the mathematical approaches
Mathematically, the problem of determining the conditions of initiation of propagating waves in excitable or bistable media is spatially-distributed, nonstationary, nonlinear and has generally no helpful symmetries, so the accurate treatment is feasible only numerically. However, the practical value of these conditions is so high that analytical answers, even if very approximate, are in high demand. Historically, there have been numerous attempts to obtain such answers, based on various phenomenological and heuristic approaches. Here we review some of these attempts, in chronological order.
Phenomenological models describing experimental relationship between the minimum stimulus amplitude required to excite an axon and the duration for which the stimulus is applied first appeared well before the physical mechanisms of biological excitability have been discovered. The study of the charge-duration relation was first carried out by Weiss [62] who experimentally derived the following linear equation where is the threshold charge and and are fitted parameters. In his original formula, Weiss did not interpret the constants and physically and hence they were later on replaced by rheobasic current and chronaxie ℎ , so that = ℎ and = ℎ [14] so (6) becomes which is known as the Weiss excitation law for the charge. An empirical equation developed by Lapicque [38,15]) reiterated Weiss's equation in a different form, for the relation between the stimulus strength and duration, i.e. the strength-duration curve. Lapicque observed that the strength of the current s required to stimulate an action potential increased as the duration s was decreased. Lapicque proposed the following current law for excitation s = ℎ 1 + s (8) which is equivalent to (6) as = s s .
Note that the rheobase current, ℎ may be defined as the minimal current amplitude of infinite duration for which threshold can be reached and that the chronaxie time of cell, refers to the value of the stimulus duration s at twice the rheobase current.
An alternative expression for threshold stimulating current was based on the idea that the nerve cell membrane could be represented by a parallel resistance and capacitance . The same Lapicque paper, and also later Blair [11,12] discussed a speculative model relying on an network to formulate the strength-duration curve. This resulted in the strength-duration relationship of the form Lapicque-Blair's exponential strength-duration curve looks similar to that given by the hyperbolic Weiss-Lapicque law (7), (8), and they nearly fit the same data [12]. Lapicque-Blair's model combines fairly accurately fit experimental outcomes with mathematical simplicity. Thus, a number of researchers began to focus on it, among which three important ones are Rashevsky [51], Monnier [44] and Hill [30]. Their results are equivalent up to the interpretation, but Hill's article is the most cited. Hill examined the relationship between the stimulus, the excitability of the tissue, and its accommodation, where the term "accommodation" [46] is used to describe the membrane potential response to a sufficiently slow increase in the stimulating current without exciting. A plausible mathematical description for a speculative dynamic variable describing the accommodation resulted in Hill's two time-constant model where , are the time constant of excitation and the time constant of accommodation, respectively. When → ∞ and = , Hill's equation (10) reduces to Lapicque-Blair's equation (9). All above approaches were phenomenological and the parameters in the strength-duration relationships were to be fitted to experimental data rather than derived from "first principles".
Study of spatial aspects of the initiation problem dates back at least to 1937, when Rushton [55] introduced the concept of the "liminal length", to represent the idea that in order to be successful, the stimulating current should excite a sufficiently large portion of the excitable cable. He supported this idea by a mathematical model of the nerve axon, which was of course linear (passive), and the active character of the membrane and the existence of the excitation threshold were taken into account in a speculative, axiomatic manner.
The situation of course changed radically after Hodgkin and Huxley have succeeded in producing a mathematical model describing the work of a nerve membrane based on experimentally established physical mechanisms. There were no more need for speculative modelling, but the real equations were strongly nonlinear, apparenly making analytical studies unfeasible and necessitating use of numerical methods. We note one such early study done by Noble and Stein [48], who used a simplification of the Hodgkin-Huxley membrane model to explore numerically the influence of the membrane activation time and accommodation on the strength-duration curve. They also deduced that in the spatially-extended context, the strength-duration curve is highly dependent on the geometry of the stimulus.
The analytical, or at least partly analytical, approaches started looking less unfeasible after the papers by McKean and Moll [42] and Flores [26]. They worked with one simple class, scalar bistable models. Considering the corresponding PDE on half-line with homogeneous boundary conditions as a dynamical system in a functional space, they have identified a critical role of one special solution, dubbed "standing wave" or, later, "critical nucleus". This is a stationary, spatially non-uniform, unstable solution, with exactly one unstable eigenvalue. Its special role is that its stable manifold forms the boundary between basins of attraction of "successful" and "unsuccessful" outcomes of a stimulation attempt.
An example of an approach seeking to take advantage of this understanding is the work by Neu et al. [47]. They considered a Galerkin projection of the infinite-dimensional dynamical system described by the PDE onto a twodimensional manifold of spatial profiles, "resembling" the shape of a developing excitation wave on a half-line (specifically, they used Gaussians). This resulted in a second-order system of ordinary differential equations (ODEs), in which the critical nucleus is represented by a saddle point, and the boundary between the basins is the stable separatrix of this saddle point. This however still left an open problem of describing analytically this separatrix. One possibility was explored by Idris [34, 2.5.2], who approximated this separatrix by the stable space of the saddle, which yielded an analytical expression for the strength-extent curve. This expression, however, produced a result that was only qualitatively correct.
The next step was using the linear approximation of the stable manifold right in the functional space. This idea was implemented by Idris and Biktashev [36] and rendered 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 excludes, for instance, cardiac excitable models. Hence a question arises, what if anything is the equivalent of the critical nucleus in such systems. An answer to this question was proposed in another work by Idris and Biktashev [35]. It happens that in excitable systems, the role of the critical nucleus is played not by stationary, but by moving solutions with one unstable eigenvalue, "critical pulses" and "critical fronts". These are unstable "counterparts" of the propagating wave solutions, existence of which was realised long before.
Finally in this brief review, the approach of [36] was extended to the moving critical solutions by Bezekci et al. [7], who also explored the possibility of using quadratic rather than linear approximation of the stable manifold. However, that paper only considered the strength-extent curve, which in the dynamical systems parlance is easier as it is about an autonomous systems. The question of strength-duration curve involves studying non-autonomous systems, and it is the subject of the present contribution.

Aims
The purpose of this paper is to quantify the strengthduration curves, as an extension of the study [36]. We investigate how the quality of approximation produced by our method depends on the parameters that define various test systems. Moreover, we investigate the feasibility of improving the accuracy by using a quadratic rather than a linear approximation of the critical manifold, and related problems. Finally, we extend the method to the case where there are no critical nucleus solutions. This is observed in multicomponent reaction-diffusion systems, where it has been previously demonstrated that instead of a critical nucleus, one has unstable propagating waves, such as critical pulses [27] or critical fronts [35].
The structure of the paper is as follows. After this introductory section, Section 2 describes the proposed analytical approach to the problem of the ignition of propagation waves in one-dimensional bistable or excitable media, from one-component with a critical nucleus to multicomponent systems having moving critical fronts and pulses, including linear and quadratic approximation of the critical manifold. The strongly non-linear nature of the equations makes it unlikely that the ingredients of these approximations can be found analytically; thus, a short outline of the numerical techniques used in the "hybrid" approach is given in the section 3. The applicability of the approach will be illistrated for five different models from onecomponent examples Zeldovich-Frank-Kamenetsky (ZFK) and McKean detailed in section 4, to multicomponent examples Na -caricature, FitzHugh-Nagumo (FHN) and Beeler-Reuter (BR) detailed in section 5. Finally, section 6 concludes the paper with a short review of the results and some possible further research directions.

Analytical theory
We aim at classification of the solutions of the system (1) set on ∈ [0, ∞), ∈ [0, ∞), supplied with the following initial and boundary conditions, in terms of their behaviour as → ∞: whether it will approach the propagating wave solution ("ignition") or the resting state ("failure"). We find it convenient to formalize the initiation problem as one posed on the whole real line ∈ ℝ, where the boundary condition at = 0 in (11) is formally represented by where (⋅) is the Dirac delta function. The principal assumption of our approach is existence of a critical solution, which is defined as a self-similar solution, which is unstable with one unstable eigenvalue. Here − is the asymptotic state behind the critical solution: − = for a critical nucleus or critical pulse, but − ≠ for a critical front. Similar to the stable wave solution, there is then a whole one-parametric family of critical solutions, Due to this translation invariance, the critical solution always has one zero eigenvalue. Hence its stable manifold has codimension two, whereas its center-stable manifold has codimension one and as such it can partition the phase space, i.e. it can serve as a boundary between basins of different attractors. Our strategy is to approximate this center-stable manifold. In the first instance, we consider a linear approximation, and in selected cases, we also explore the feasibility of the quadratic approximation.

Linear approximation
Let us rewrite the system (12) in a frame of reference moving with a constant speed , so that We linearize this equation on the critical solution, which is stationary in the moving framẽ The linearization gives where is the Jacobian matrix of the kinetic term, evaluated at the critical solution and is the forcing term as measured in the moving frame of reference. Equation (17) is a linear non-homogeneous equation, with time-independent linear operator, For simplicity of the argument, we assume that the eigenfunctions of , are simple and form a basis in an appropriate functional space, and the same is true for the adjoint  + 1 . Another assumption, which simplifies formulas and is true for all examples considered, is that all eigenvalues important for the theory are real. We shall enumerate the eigenpairs in the decreasing order of , so by assumption we always have The solution of (23) is By assumption, 1 > 0, and due to translational symmetry, 2 = 0, and the rest of the spectrum is assumed within the left half-plane. Since the stimulation is supposed to be finite in time, ℎ ( ) ≡ 0 for > s . Therefore, the condition of criticality is from which we seek to obtain the critical curve based on this linear approximation.

General Setting
Using the definitions of 1 (0) and ℎ 1 ( ′ ), we have, in terms of the original model, The forcing term is defined as hence (27) gives This is a finite equation relating s and s so, in principle, gives the answer. However, it contains the parameter which still has to be decided upon. This question occurs already for the strength-extent curve, and we refer the reader to [7] for a detailed discussion. The new issue here is that for ∈ [0, s ] we are now dealing with a time-dependent problem.

The case of critical nucleus
This is the case when = 0, i.e. the critical solution is stationary, and moreover it is even in . Then there is a natural choice of = 0 prescribed by symmetry. Hence, (29) gives the classical Lapicque-Blair formula [38,11] where the rheobase is

The case of moving critical solution
In this case there is no → − symmetry and the choice of is no longer trivial. We follow the ideas discussed in [7]. According to those, we assume that the linear approximation works best if the initial value for is the smallest in some sense. Our heuristic is that it will be the smallest in the 2 norm if not only 1 = 0, but also 2 = 0. Moreover, small shifts of the critical solutions are equivalent to adding a small amount of 2 . Hence an appropriate shift can achieve that 2 = 0. So we adopt 2 = 0 as the criterion of selecting . The only modification of this idea for the present case is that we apply this condition not at = 0, but at the moment from which the system is autonomous, i.e. at = s . Employing both 1 s = 0 and 2 s = 0 results in following system of equations where the right hand sides  1 and  2 are constants, defined entirely by the properties of the model, System (32) is a nonlinear system of two equations for two unknown parameters, s and . By eliminating the parameter s , we find the compatibility condition as follows: This can be further simplified by using the following change of variable, that leads to The finite equation (35) imposes a connection between the shift and the stimulus duration s . This connection defines as an implicit function of s . After finding the value of , one only needs to employ this value in one of the compatibility conditions in (32) in order to find the amplitude s since both produce the same result. This completes the construction of the linear approximation of the strength-duration curve s ( s ).

Quadratic approximation of the stable manifold
In this subsection, we restrict consideration to the case of a critical nucleus.
For simplicity, rather than using the matrix notation as in the linear approximation, we shall now proceed with an explicit notation for the components of the reaction-diffusion systems. We use Greek letters for superscripts to enumerate them, and adopt Einstein's summation convention for those indices. In this way we start from the generic reactiondiffusion system where overdots denote differentiation with respect to time, subscripts (⋅) denote differentiation with respect to space and Greek subscripts after a comma designate a partial differentiation by the corresponding reactive components. The right and left eigenfunctions are defined respectively by where ∈ {1, 2, 3, ⋯}, and the biorthogonality condition is We consider only even solutions, so in subsequent sums only those that correspond to even eigenfunctions are assumed.
We seek solutions in the form of generalized Fourier series in the right eigenfunctions, where the Fourier coefficients are defined by Time-differentiation of this giveṡ and We assume that eigenvalues are real and ordered from larger to smaller, 1 > 0, 2 = 0 is of course the eigenvalue corresponding to the translational symmetry and an odd eigenfunction 2 =̂ ′ , and < 0 for ≥ 3. Our task is to determine the conditions on the initial values of the Fourier coefficients that would ensure that 1 (∞) = 0, which means that the trajectory approaches the critical nucleus, so the initial condition is precisely at the threshold.
Let us rewrite the system (36) as an equivalent system of integral equations, . Successive approximations to the solution can be obtained by direct iterations of this system, Taking (0) = 0 for all , we have The requirement (1) 1 (∞) = 0 recovers the linear approximation. The next iteration produces (2) Note that e ( + − 1 ) → 0 as → ∞ because , ≤ 3 < 0 for , ≥ 3, so upon exchanging the order of intergration and summation, we have converging improper integrals. The requirement (2) 1 (∞) = 0 leads to a quadratic equation for s ,

Hybrid approach
With a few exceptions, the ingredients for the expressions used in the linear and quadratic approximations of the strength-duration curve, starting from the critical solution itself, are not available analytically and have to be found numerically. In this section we describe numerical methods we used to find these ingredients. We divided this section into two subsections. for models with self-adjoint and non-selfadjoint linearization operator.

Ingredients of the one-component systems
This corresponds to the case of the critical nucleus and for the linear approximation, we need to have the knowledge of̂ , 1 , 1 while for the quadratic approximation, ideally the whole spectrum of , , = 1, 3, 5, ⋯ is needed. Here we shortly describe the methods to obtain the mentioned ingredients in algorithmic forms. For a more detailed explanation, see [7].
In order to find the critical nucleus, we take advantage of the fact that its center-stable manifold has codimension one, and divides the phase space into two open sets, one leading to successful initiation and the other to decay [26,27,41,43,47,36]. This means that the critical trajectories, corresponding to the stimulus strength exactly equal to the threshold, tend towards the critical nucleus as → ∞, whereas the stimulus strength slightly above or below the threshold produces the solution that gets close to the critical nucleus and stays in its vicinity for a long time, before deviating from it to propagate or to collapse. This can be used to calculate an approximation of the critical nucleus, see Algorithm 1. The calculations are done for (1,11) for ∈ [0, ], where is chosen large enough for the results to be not significantly different from those for ∈ [0, ∞).
Input: Pre-found value of * s for a chosen s Output: Critical nucleuŝ • Find ( , ) by solving initial value problem (1,11).
We calculate the eigenpairs of the linearized operator  defined by (20) (note that in the present case = 0) using a variant of the power iteration method. We use a random number generator to assign linear independent initial guesses for 1 , 2 , . . . , choose a time domain ∈ (0, ), and then follow Algorithm 2 until the desired convergence criterion is fulfilled. The appropriate values of of course depend on the spectrum of  which may not be known a priori; our choice was entirely empirical. Note that the algorithm is described without prejudice to the choice of the norm used for then obvious simplifications are possible. Also, we use the convergence criterion based on the change in each eigenvalue; it can also be done in terms of the change, say in 2 -norm, in each eigenfunction.

Ingredients of the multi-component systems
The non-stationary critical solutions, observed in multicomponent systems, can be found using an appropriate modification of Algorithm 1, exploiting computations in a comoving frame of reference, as described in [7]. However, more accurate results can be obtained by continuation of the boundary-value problem (14), an autonomous system for vector-function̂ ( ) and scalar . For the critical pulses − = , and our strategy is to use periodic solutions with very long periods as approximation to pulses. We aim to calculate conduction velocity restitution curve [60], that is, a one-parametric family of solutions of the following periodic boundary-value problem: where is the spatial period of the waves. When the problem is well posed, (43) defines a curve in the , , ( ) space. In the limit → ∞, this curve splits into two branches, the upper branch with a stable propagating pulse solution, ( w , w ( )) and the lower branch with an unstable critical pulse solution, ( ,̂ ( )), which is of interest to us. We performed the continuation using AUTO [19]. To obtain the periodic solutions, we consider an extension of (43) by an extra parameter corresponding to "stimulation current" added to the transmembrane voltage equation. Starting from an initial guess of w , the continuation is done in accordance with Algorithm 3.
• Equilibrium, resting state: As a final step, we calculate the left and right eigenfunctions by means of the Gram-Schmidt orthogonalization process, modified with account of the fact that  is now not self-adjoint, Algorithm 4.

Zeldovich-Frank-Kamenetsky equation
Our first application example is the one-component reaction-diffusion equation, first introduced by Zeldovich and Frank-Kamenetsky (ZFK) [63] to describe propagation of flames; it is also known as "Nagumo equation" [41] and "Schlögl model" [57]: where we assume that ∈ (0, 1∕2). The critical nucleus solution̂ = ̂ for this equation can be found analytically [26,36] The other two components required for the definition of critical curves in the linear approximation are 1 and 1 = 1 = 1 which are solutions of 2 Actually, expressions given in both of these works contain typos.

Input: A linearly independent set
| ≤ tolerance ∀ Algorithm 4: Numerical computation of principal eigenpairs of non-self-adjoint operator  by "marching".
We have been unable to find solution of this eigenvalue problem analytically. We note, however, that̂ given by (45) is unimodal, thereforê ′ , which is the eigenfunction of  corresponding to = 0, has one root, hence by Sturm's oscillation theorem,̂ ′ = 2 and 2 = 0, and there is indeed exactly one simple eigenvalue 1 > 0 and the corresponding 1 solving (46) has no roots.

The small-threshold limit and the "fully analytical" result
In this subsection we extend the results of [36] in the parameter space and correct some typos found in the paper. For ≪ 1, the critical nucleus (45) is  ( ) uniformly in , and is approximatelŷ In the same limit, the nonlinearity can be approximated by ( ) ≈ ( − ). With these approximations, problem (46) has the solution and (30) then gives an explicit expression for the strengthduration curve in the form with the following rheobase and chronaxie form [36] rh = 45 64 3∕2 , Remark that the performance of the resulting approximation based on the analytical expression for the strengthduration curve (49) and (50) can be further improved by obtaining the essential ingredients numerically. This is done considering (49), in which the rheobase is, instead, defined according to (31), as The plot of the hybrid numeric-asymptotic prediction is compared with the direct numerical simulations as shown in fig. 3(b). As depicted in the figure, reasonable agreement between the two data sets is observed when the threshold parameter is small.
It should be noted that the strength-duration curve approximation remains above the a priori lower bound (42) , for all s .

Hybrid approach
Numerical computation of the essential ingredients needed for linear approximation of the critical curves is carried out using Algorithms 1 and 2 described in 3.1. Fig. 2 illustrates the processes of numerical computation of the critical nucleus and the first eigenmodes in the ZFK equation, for the threshold parameter varying from 0.05 to 0.45 with the increment 0.1. The stimulation is done by fixing the duration time at the value s = 1.5. To obtain the minimum of ( ) and # = argmin( ( )), the bisection loop is terminated as soon as the absolute difference between upper and lower estimate for the threshold is sufficiently small, i.e. | s − s | < 10 −5 . For each case, the solution ( , # ) of the nonlinear problem = + ( ) provides an estimate of the critical nucleus.

Quadratic theory
To estimate a few principal eigenmodes of the ZFK equation, we have considered a finite interval ∈ [0, ) as an approximation of ∈ [0, ∞). We find only a few approximate eigenvalues in the discrete spectrum, while the remaining ones are in the continuous spectrum. The eigenfunctions corresponding to the discrete eigenvalues are well localized whereas those corresponding to the continuous eigenvalues are evidently non-localized, and thus they cannot be taken into account in quadratic approximation. We observe that at increasing values of , the eigenfunctions 1 and 3 corresponding to the discrete eigenvalues 1 and 3 are well localized towards the left end of the interval ∈ [0, ], whereas those corresponding to the continuous eigenvalues

McKean equation 4.2.1. Model formulation
Our second example is a piece-wise linear version of the ZFK equation, considered by McKean in [41] and then also in [52]: where we assume that ∈ (0, 1∕2). The critical nucleus solution in this equation is found in a closed form, where * = obtained by the fact that̂ ( ) and its derivative are continuous at this point. The eigenvalue problem can be expressed as where the linearization operator contains the Dirac delta function: The principal eigenvalue and the corresponding eigenfunction can be written in the form where and W 0 (⋅) is the principal branch of the Lambert W-function as defined e.g. in [17].

Hybrid approach
In this model, since the exact analytical solution for the critical nucleus and the ignition eigenpair are known for an arbitrary ∈ (0, 1∕2), the "hybrid approach" is not necessary. However, for technical purposes, we address it here as well, to show the numerical computation of the essential ingredients based on Algorithms 1 and 2 works satisfactorily even for the models with discontinuous right hand sides.
Due to the discontinuous terms, the numerical computation of the ingredients of the McKean equation requires the finite-element treatment which was outlined in our previous paper [7]. Hence, we skip the details here. Fig. 5 illustrates the processes involved in obtaining the critical nucleus and ignition mode. The results of these ingredients are compared with their analytical counterparts and we see a good agreement between the two.

Linear theory
Linear approximation of the strength-duration curve can be found using the analytically derived expression given by (30). However, it must be noted that in this case, the rheobase is found as where  = sinh * + + 1 cosh * This linear prediction formalism compared with the direct numerical simulations is depicted in fig. 6(a). The a priori bound for these chosen threshold parameters is outside of the duration domain. As shown in the figure, the linear approximation for parameter values close to 1∕2 better fits to the numerical simulation than that for small values. This may be related to the fact that the leading eigenvalue is inversely proportional to the threshold parameter . Even for larger , there are still some deviations between the linear theory and numerical simulations, which can be reduced by considering second order approximation that will be outlined in the following subsection.

Quadratic theory
Linear approximation of the strength-duration curve can be improved by considering the second order approximation. For the quadratic approximation, the knowledge of the whole spectrum is ideally required. We know that the linearization spectrum of the critical nucleus has only one unstable eigenvalue 1 > 0, and due to translational symmetry, 2 = 0, and the rest of the spectrum lies entirely in the left half-plane. The first aim of this part of the subsection is to find these remaining stable eigenvalues and corresponding eigenfunctions. To obtain these eigenpairs, we replace the infinite interval [0, ∞) with a finite interval [0, ] with a homogeneous Dirichlet boundary condition at = , aiming to consider the limit → ∞. The solution of eigenvalue problem (56) in this case is which is dependent of the domain size as opposed to the eigenfunction expression. After finding analytical expressions for the eigenpairs, the next step is to calculate (37), (38) and (39), and then substitute them back in the coefficients of the quadratic equation for s (40) so that the second-order approximation of the strength-duration curve can be generated. Fig. 6(b) shows graphs of the linear and quadratic approximation of the strength-duration curve along with its numerical result for = 0.4. The quadratic approximation was obtained for = 10 and 287 eigenvalues. The accuracy of the second-order approximation is much closer to the direct numerical simulation compared to the first-order approximation.

Model formulation
Our next example is the caricature model of an Nadriven cardiac excitation front suggested in [8]. It is a twocomponent reaction-diffusion system (1) with = ( , ℎ) ⊤ , and H(⋅) is the Heaviside step function. The component of the solution corresponds to the nondimensionalized transmembrane voltage, and the component ℎ describes the inactivation gate of the fast sodium current, which is known in electrophysiology as Na and which is mainly responsible for the propagation of excitation in cardiac muscle in the norm. A special feature of this model is that there is a continuum of potential resting/pre-front states, and a continuum of potential post-front states, so any front solution connects a point from one continuum to a point from another continuum. The critical solution̂ = (̂ ,ĥ) ⊤ is described bŷ where the post-front voltage and front thickness Δ are given by and the front speed is defined by an implicit equation or equivalently where = 2 , = ∕( + 1).
For the analytical expression of the first two left and right eigenfunctions, please see [7,34].

Hybrid approach
Even though we know the ingredients of the linear theory for this model analytically, we still found them numerically as well. The hybrid approach is needed not only because it helps to validate the analytical result but also because, in some cases, it is the only option as the analytical derivation is not always possible. Due to the discontinuous right-hand sides, it is essential to use the standard finite element method, at least when dealing with these discontinuous terms. The complete discretization formula for the critical front of the model is presented in Appendix B.
For two initial conditions, fig. 7 shows the evolution of component in the comoving frame of reference. For each case, the solution approaches the critical front, i.e. the solution at = 120 in the figure and then gives rise to the stable propagating wave if the initial condition is above the threshold or decays back to the resting state otherwise. Fig. 8 gives the comparison of the numerical critical front obtained using operator splitting method and its analytical closed-form solution given by (64). We can see that the shooting procedure provides a good approximation of the critical front for the selected parameters.
Numerical experiments suggest that there are two values of the speed , and , satisfying 0 < < < ∞, such that the faster fronts are higher and stable, and the slower fronts are lower and unstable, hence a slower front either dissipates or increases in the speed and magnitude to the fast branch solution depending on the initial condition being below-and above-threshold, respectively [8,32]. This can be seen in the right panel of fig. 7 where the blue circle and green square symbols represent fast and slow front speeds for the selected pair of = 1, = 8.2, and the red line indicates how the speed of the front changes in time. For initial conditions slightly above threshold, the front speed gets closer to the slow speed and stays in the vicinity of it for a long hnum. time before developing into the fast speed while the initial condition slightly below threshold results in the front speed to drop to zero eventually. The next step after finding the critical front of the Nacaricature model is the determination of the right and left eigenfunctions along with the corresponding eigenvalues employing Algorithm 4 detailed in 3.1. In fig. 9, the eigenfunctions obtained using the hybrid method fairly resemble exact analytical eigenfunctions. The largest difference between the numerical and analytical eigenfunctions is observed in the vicinity of the discontinuous values, = 0 and = −Δ. This is totally expected as the numerical scheme used to calculate the eigenfunctions is the second-order accurate, except near discontinuities, where it reduces to first order accuracy and introduces spurious oscillations due to the Beam-Warming method [21].

Linear theory
The result of the calculation of the strength-duration curve for Na -caricature model is visualized in fig. 10, where the linear approximation is based on the formulas (32) and (35). For every chosen stimulus duration, s , we calculate the zeros of (35) in order to find the value of the shift . We then substitute this value of into one of the equations in (32) (both produce the same result) to get the corresponding value of s . In the simulations, we choose two different set of the model parameters, = 7.8, = 2∕3 and = 8, = 9∕11 from which the resulting curves are respectively shown in fig. 10(a,c) and fig. 10(b,d). The shape of ( ) is rather similar for both cases and the main difference between the two is the closeness of the values for two different duration of stimulus values, s = 3 and s = 10. We observe that the theoretical critical curve for the first set of parameter values is well adapted to the direct numerical simulation threshold curve for smaller values of s , and then bends down dramatically as the value of s increases. The theoretical prediction for the second set of parameters values, however, gets better with s .

Model formulation
The FitzHugh-Nagumo (FHN) system is a twocomponent reaction-diffusion system, which could be considered as a ZFK equation extended by adding a second, slow variable, describing inhibition of excitation. It is probably the single historically most important model describing excitable media. We consider it in the form

Hybrid approach
System (69) has an unstable propagating pulse solution as opposed to its reduced form, the ZFK equation with nontrivial stationary solution. It is known (see e.g. [27] and references therein) that in the limit ↘ 0, the cricial pulse solution whose -component is small and -component is close to the critical pulse of the corresponding ZFK equation. This makes it feasible to obtain explicit analytical solutions by using perturbation techniques in the double limit ↘ 0, ↘ 0. These asymptotics will be described in a separate publication, and here we describe only the hybrid approach. The critical pulse is obtained by applying Algorithm 3 by means of AUTO. The corresponding CV restitution curves are illustrated in fig. 11. For the critical pulses, we take the solutions at lower branches at > 7.5 ⋅ 10 3 (see top row of fig. 12). Other essential ingredients of the theory are the first two left eigenfunctions and the first leading eigenvalue, which are computed using Algorithm 4. The eigenfunctions for the two selected cases look rather similar, as shown in fig. 12. Fig. 13 illustrates the calculation of the strength-duration curve for FHN model for = 0.05 and = 0.13 according to the formulas (32) and (35). The equation (35) has two roots, one of them is negative close to zero and the other is positive. We find that in both cases the smaller root denoted by blue circle and red square points in fig. 13(a,b) gives the corresponding value of . The critical curves compared with those obtained from direct numerical simulation are sketched in fig. 13(c,d). From this plot, it can be seen that the theoretical prediction for both values of parameter is almost equally close to the numerical prediction.

Model formulation
Here, we looked at a variant of the classical Beeler-Reuter (BR) model of mammalian ventricular cardiac myocytes [6], modified to describe phenomenologically the dy- For the detailed description of the components of ( ), please see the appendix of [7].

Hybrid approach
As in the FitzHugh-Nagumo system, the critical solution is a moving pulse, and thus, we obtain the CV restitution curves and the critical pulse in a similar way. The CV restitution curves for the modified Beeler-Reuter model are shown in fig. 14. Apart from the critical pulse, the solution at lower branches, the knowledge of the first two left eigenfunctions and the first leading eigenvalue are also required. These ingredients have been found by the marching method given by Algorithm 4. The essential ingredients of the theory for two different data sets are sketched in fig. 15.   (32) and (35). Firstly, the values of are determined by the transcendental equation (35) and compared to FHN system, it is easier to detect the zeros of this equation, two of which are shown in the top panel of the figure for s = 3 and s = 10. Then, the remaining part is to insert the found value of back into theoretical threshold curve generated by (32). The bottom panel of the figure shows these threshold curves being compared with numerical critical curves. As can be seen from the figure, the analytical estimate for = 0.115 provides a somewhat better approximation than that for = 0.105.  1, 1, 1, 1, 1, 10 5 ). The space coordinate is chosen so that = 0 at the maximum of̂ . Correspondence of lines with components is according to the legends at the top.

Discussion
The main aim of this paper was to extend the method proposed in [36] for analytical description of the threshold curves that separate the basins of attraction of propagating wave solutions and of decaying solutions of certain reactiondiffusion models of spatially-extended excitable media. Specific aims are: • Extending the proposed theory to analysis of a wider class of excitable systems, including multicomponent reaction-diffusion systems, systems with nonself-adjoint linearized operators and in particular, systems with moving critical solutions (critical fronts and critical pulses).
• Building an extension of this method from a linear to a quadratic approximation of the (center-)stable manifold of the critical solution to demonstrate the discrepancy between the analytical based on linear approximation and numerical threshold curves encountered when considering this quadratic approximation.
The essential ingredients of the theory are the critical solution itself, and the eigenfunctions of the corresponding linearized operator. For the linear approximation in the critical nucleus case, we need the leading left (adjoint) eigenfunction; in the moving critical solution case, we need two leading left eigenfunctions; and for the quadratic approximations we require as many eigenvalues and left and right eigenfunctions as possible to achieve better accuracy. Of course, closed analytical formulas for these ingredients can only be obtained in exceptional cases, and in a more typical situation a "hybrid" approach is required, where these ingredients are obtained numerically. We thus have provided insight into how the numerical computation of these essential ingredients can be done. The theory have been demonstrated on five different test problems ranging from one-component reaction-diffusion systems where the critical solution is the critical nucleus to the multicomponent test problems with either critical front or critical pulse solution. In all models, the analytical threshold curves are compared with the numerical simulations obtained using the bisection algorithm. We have applied both linear and quadratic approximations for one-component test problems, ZFK and McKean models. The quadratic approximation agreed much better with numerical threshold curves compared to the linear approximation's results, as would be expected.
The accuracy and efficiency of the hybrid computation of the essential ingredients is dependent on the numerical scheme and mesh resolution. It is obvious that some of the numerical schemes discussed in this paper do not outperform other long-running and mathematically more complicated numerical schemes. In particular, the numerical study of Na -caricature model has introduced the spurious oscillations and first-order accurate result near the discontinuities due to Beam-Warming method, which can be tackled using some advanced shape-preserving advection schemes (see, for example, [54]). Hence, such approaches that demand high computational cost can be carried out as an interesting direction for future research if higher accuracy is required.
As the results of the theory pointed out, our method provides more accurate results for some parameter values than the others, especially in linear analysis. Even though the quadratic approximation offered for one-component test problems with the critical nucleus solutions provides more accurate estimates, still it is not fully understood why the choice of parameter values significantly matters. Hence, this remains quite important line of research. On the other hand, the proposed theory based on the moving critical solutions involves only the linear approximation. As an additional consideration, quadratic approximation for the cases of moving critical solutions can also be carried out.
The theory established in this paper is limited to one spatial dimension. Therefore, it could be of interest for further research to adapt the theory to two and three dimensions.
Throughout the theory, we have made the assumption that the spectrum is real. This is, however, not necessarily the case for the non-self-adjoint problems, which remains an interesting direction for future research.
Another extension of the work worth considering would be to investigate the theory on some up-to-date realistic cardiac excitation models [16], simplified cardiac models with unusual properties [20,31], and other excitable media such as combustible media [33,58,39], pipe flow [4,5] etc.In the context of this Special Issue, of particular interest is the link to the problem of threshold of front propagation in neural fields, discussed in the recent paper [23], where the Lapicque-Blair strength-duration relationships also naturally occurs under certain simplifying assumptions. Data statement: The research data supporting this publication are provided within this paper.

A. Discretization formula for strength-duration curve
With the aim of comparing the explicit approximations for the strength-duration curve, we describe how the threshold curve is obtained using direct numerical simulations in this section. The numerical strength-duration threshold curve is computed by solving the nonlinear system (1) for the initial and boundary conditions given by (4)-(5) using standard finite difference or finite element discretization. More specifically, for ZFK, FHN and BR models we use finite differences, and for McKean and Na -caricature models we implement finite element method instead.

A.1. Finite Difference Discretization Formula
The discretization formula for the generic form of initialboundary value problem We employ explicit first order Euler forward difference approximation in time and explicit second order centered finite difference approximation in space and plugging these into (73), we obtain in conjunction with its initial and boundary conditionŝ

A.2. Finite Element Discretization Formula
Two of our test problems namely, McKean and Nacaricature models, contain discontinuous kinetic terms.

Discretization Formula for McKean equation
We use even extension of the problem and start with one component McKean model In the Galerkin finite element method, this can be written in the vector form, foř , where is the the mass matrix, is the stiffness matrix and is the load vector (see [7] for a crude introduction to finite element method and the derivation of the matrices).
The vector has only one nonzero entry by definition of the Dirac delta function. Finally, we employ the generalized trapezoidal rule (also known as scheme) [13], in which the residual is evaluated at + , with this notation implyinǧ where 0 ≤ ≤ 1 is a real parameter. Based on the above, the fully discrete problem (77) becomes For = 0 , the linear system (78) is the explicit Euler method that has a stability condition to be satisfied and its truncation error is  Δ +  Δ 2 , = 1∕2 gives the second-order unconditionally stable Crank-Nicolson method with truncation error  Δ 2 + Δ 2 , and = 1 gives the first-order accurate implicit Euler rule that is also unconditionally stable and its truncation error is  Δ +  Δ 2 [13].

Discretization Formula for the Na -caricature model
Even extended version of the model is in the following form The fully discretized finite element discretization formula for the model is The matrix is a tridiagonal matrix with following diagonal elements: We need to implement no flux boundary condition In these formulas, we use a shorthand notation, Having in mind that the tent functions we have chosen are piecewise linear, the pointš −1 ,̌ +1 ,̌ 2 anď −1 are then obtained from the linear interpolation methoď The vector has entries for = 2, 3, ⋯ , − 1 and on the boundaries The shorthand notations used above are,

A.3. Threshold curve
To obtain the threshold curve in the stimulus strengthduration plane, we solve a sequence of the "stimulation by current" initial-value problem (12) and (13). The choice of the numerical scheme changes according to the specific model as defined above. The computation is done by fixing stimulation time duration s and varying the strength of the current s . For any initial upper estimate s (superthreshold), known to be sufficient for ignition, and lower estimate right hand sides of the equations with discontinuity terms involving the Heaviside step function. For the standard finite difference discretization, we use Beam-Warming scheme for the first spatial derivatives of both and ℎ. We set the domain of and coordinates to be − ≤ ≤ , 0 ≤ ≤ so that the grid points , are where Δ > 0 and Δ > 0 are fixed space and time steps, and = 0, 1, ⋯ , , = 0, 1, ⋯ , for , > 0. For convenience, we use the following shorthand notations: , ℎ as the numerical approximation of , and ℎ , , * = 0 as our pinning condition which will be further explained later, and finally = ′ as the speed at -th time step. Using these representations, we solve (81) numerically using operator splitting method approach in the following seven steps: Step 1(Finite element method): As a first step, we solve equation without the advection term as it is multiplied by the speed which is not determined yet We have to employ the finite element method due to the discontinuous Heaviside step function which gives Step 2(Thomas algorithm): Inserting the elements of the matrices , , and the discretized solution and ℎ (both known) into (84) yields a tridiagonal system of equations in a following form This can be solved using a standard Gaussian elimination method such as Thomas algorithm and using such algorithm is sometimes crucial as it leads to a reduced computational cost. The back substitution procedure (see e.g. [61] for more detailed explanation) generates the solution as Step 3(Finding the value of the speed) : We have divided equation in (81) into two parts and it remains to find the solution of the advection step in the splitted scheme. Before we update the solution, it is necessary to find the value of the speed according to the pinning condition where the index * corresponds to an integer value, indicating the spatial position at which the solution is equal to zero initially, i.e. * = 0. As the Beam-Warming method is second order accurate, we find the speed value using the Beam-Warming scheme by means of the discretized solution found in the previous step as, Step 6 (Thomas algorithm): As (88) is also a tridiagonal matrix, we can employ the Thomas algorithm here as well, similarly to the case of equation. Step The numerical computation of the critical front is achieved by solving (81) according to above seven-step procedure. As the Na -caricature model is a two-component system, to find the numerical estimation of the critical front, we calculate ( ) as
The subdiagonal elements of are also found as and on the boundaries, we have The shorthand notations used above are, Step 2 (Beam Warming scheme): The second step is to solve the pure advection equation using the Beam-Warming , (94) that finishes the numerical scheme of component of the linearized problem.
Step 3 (Finite element method): Again the finite element method is conveniently employed to numerically solve the first equation of ℎ component, that results in  , that completes the numerical solution of the linearized equation.

B.3. Discretization Formula for the Adjoint Linearized Problem
The adjoint linearized problem for the Na -caricature model is ̄ = This problem can also be solved in 4-steps as follows: Step 1 (Finite element method) : We begin with the equation of̄ component and solve first the following, with the solution based on the finite element method, Step 2 (Beam Warming scheme): As the advection term has negative sign in the front, the Beam-Warming numerical scheme is in this case, .
We note that the Thomas algorithm can also be applied to the diagonal systems (93), (95), (97) and (98) 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 critical curves requires the knowledge of the critical front as well as first two leading eigenvalues and corresponding left and right eigenfunctions. The numerical calculating of the eigenpairs are determined by the method discussed in Section 3, and these two last subsections are dedicated to how to derive the first two eigenmodes as a numerical solution of the linearized and adjoint linearized problems by means of the operator splitting method. Alternatively, the implicitly restarted Arnoldi method [50] can be used to estimate these essential ingredients, in which case we use the matrix representations of the discretized versions of the equations (92) and (96).