{'fraction_non_alphanumeric': 0.03532863227471149, 'fraction_numerical': 0.012459572663167382, 'mean_word_length': 4.441815733794961, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Predicting critical ignition in slow-fast excitable models

Linearization around unstable travelling waves in excitable systems can be used to approximate strength-extent curves in the problem of initiation of excitation waves for a family of spatially confined perturbations to the rest state. This theory relies on the knowledge of the unstable travelling wave solution as well as the leading left and right eigenfunctions of its linearization. We investigate the asymptotics of these ingredients, and utility of the resulting approximations of the strength-extent curves, in the slow-fast limit in two-component excitable systems of FitzHugh-Nagumo type, and test those on four illustrative models. Of these, two are with degenerate dependence of the fast kinetic on the slow variable, a feature which is motivated by a particular model found in the literature. In both cases, the unstable travelling wave solution converges to a stationary ``critical nucleus'' of the corresponding one-component fast subsystem. We observe that in the full system, the asymptotics of the left and right eigenspaces are distinct. In particular, the slow component of the left eigenfunction corresponding to the translational symmetry does not become negligible in the asymptotic limit. This has a significant detrimental effect on the critical curve predictions. The theory as formulated previously uses an heuristic to address a difficulty related to the translational invariance. We describe two alternatives to that heuristic, which do not use the misbehaving eigenfunction component. These new heuristics show much better predictive properties, including in the asymptotic limit, in all four examples.


I. INTRODUCTION
Excitable systems represent a distinct modeling legacy in the context of biological systems, especially regarding cardiac and neuronal cells. These models utilize the existence of a typically state-dependent threshold to distinguish between relaxation and transiently amplified dynamics. In classic FitzHugh-Nagumo type models this threshold is explicitly the unstable branch of the fast variable nullcline, volumes have been dedicated to investigations around the response of excitable models to driving.
In the nomenclature of previous efforts [1], we focus the formalism of a "stimulation by voltage", which corresponds to an initial value problem for reaction diffusion models of the form, where u = [u 1 , u 2 ], D = diag (1, 0), Γ = diag (1, γ). When ∂ xū = 0 and f (ū) = 0 uniquely, thenū is the asymptotically stable rest-state. The time-scale separation of the dynamics of u 1 and u 2 is controlled by γ, the ratio of relevant kinetic timescales. The limit γ → 0 designates the fast subsystem and the transition from moving solutions (γ > 0, c = 0) to stationary solutions (γ = 0, c = 0) of (1). In this work we concern ourselves primarily with the transient dynamics in the vicinity of the traveling wave solutions of these slow-fast systems. Traveling wave solutions of (1) satisfy a nonlinear eigenvalue problem posed on the real line x ∈ R, for the wave solutionũ and the associated wave speedc = 0. Additionally, we require that the wave approach the rest state asymptotically,ũ(x → ±∞) →ū so that the solution is localized in space, and a homoclinic connection to and from the rest state in the co-moving frame. The simplest solutions to (2) are two single-pulse waves: a faster, stable waveǔ, and the slower, unstable waveû, with |ĉ| < |č|. When speaking about general solutions of (2), we shall refer to generic wavesũ and associated generic speedsc. The linear stability of these waves is determined by the eigenspectrum (σ i , v i ) of the operator L, the linearization about (2), which is guaranteed to have a marginal eigenfunction, Lv = 0v, where v = ∂ xũ . Further, for the unstable asymptotic wave solutionû, the operator L must have one unstable mode, Re(σ) > 0, whose shape describes the fastest-growing mode in the co-moving frame. Similarly, the inner-product over the domain, defines the set of adjoint eigenfunctions w j with eigenvalues σ j satisfying the biorthogonality condition (σ j − σ i ) w j | v i = (σ j − σ i )δ ij . Equivalently, we may consider the right eigenfunctions of the adjoint linear operator, L † , where (·) stands for transposed and (·) † for Hermitian conjugate.
The transient dynamics in the vicinity of the unstable solutionû is important phenomenologically for understanding the initiation of e.g., electrical waves in heart tissue. Ref. [1] used the properties of the unstable wave to predict the minimal perturbations to the rest state which lead to the generation of new excitation waves, using a particular notion of 'smallness' for the constructed perturbations which relies on assumptions about the behavior of the solutionû and its linearization in the limit γ → 0. In this work, we investigate the limitations of those assumptions, and test alternative methods for constructing those minimal perturbations. We find that in several cases, the applicability of the existing theory to traveling wave solutions of simple slow-fast models predicts at the very least sub-optimal or occasionally unrealistically large minimal perturbations. We have in mind three specific examples of the models: the classical FitzHugh-Nagumo model in the original formulation [2] (FHN), the two-variable Karma 1994 model [3] (Karma) and two-variable reduction of the Fenton-Karma model due to Mitchell and Schaeffer [4] (MS). To our surprise we found that Karma model is different from the other two in its asymptotics in γ, which has proved to be due to a specific form of dependence of the fast kinetics f 1 on the slow variable u 2 , namely via u 2 M , with M > 1, whereas in the resting stateū 2 = 0. To illustrate further the specifics of such dependence, we have also considered a variation of the FitzHugh-Nagumo model which also has this feature, that is, f 1 depends on u 2 3 ; we shall call it "FHN with cubic recovery", or FHNCR for short.
The outline of this paper is as follows. First we review the basics of sub-and super-threshold response for classic excitable models of FitzHugh-Nagumo type, the complication from embedding them in spatially extended media, and how to distinguish sub-and super-threshold excitations in space. Second, we review the essential ingredients of the linear theory of critical excitation, in the context of the asymptotics in the slow-fast time scale separation between the activator and inhibitor subsystems. We show that the approach applied to slow-fast systems relies on misleading assumptions about the asymptotic structure of the leading left and right eigenspaces and particularly the form of the adjoint eigenfunctions. We detail the computation of the slow wave solutions, their eigenspectra, and the solution of transient trajectories by direct numerical simulation. Third, we propose heuristics based on minimization principles which do not rely on the slow mode corresponding to translational symmetry, and which outperform the previously suggested skew-product motivated heuristic. Finally we reinforce this conclusion with numerical examples using several slow-fast systems, and use these results to infer properties of the leading eigenspace and the effect of degenerate nonlinearities on observing other nearby saddle solutions.

II. THEORY
Given an initial condition u(0, x) =ū +h(x; x s , U s ) parameterized by the properties (x s , U s ) of the perturbation to the rest state for which configurations will the initial conditions eventually recruit the entire domainexcite the medium -and for which configurations will it return directly to the rest state. In the language of coherent structures, this corresponds to identifying the boundary of the basin of attraction isolating the stable wave solution from the uniform rest state, and projecting this infinite-dimensional manifold down to the shapemodifying parameter space of (x s , U s ). Throughout this work we shall use a parameterized perturbation to the rest state,h where e 1 = 1, 0 and H(x) is the Heaviside distribution, so that

A. Linear theory of critical excitations
Here we present a brief motivation for the linear theory of critical excitations, and recount the assumptions of the method. Given the initial state u(0, x) =û(x) +ĥ(x), whereĥ(x) is understood to be perturbatively small, then the dynamics of the state subject to (1) can be understood through the linearization aboutû(ξ), with ξ = x −ĉt − s the co-moving frame coordinate. Expressing the linearized dynamics in terms of the spectral expansion, and recalling that Re(σ 1 ) > 0, then we can consider the requirement that the sole unstable mode is not excited due to the perturbation, where we identify the modal amplitudes of the linearization, a i (0) = w i ĥ (ξ) . At time t = 0, we express the initial perturbationĥ to the critical solutionû in terms of a perturbationh to the rest stateū so that ξ + s = x, where the invariance of the rest state with respect to translational symmetry manifests as a freedom in the origin of the rest-state perturbation, the shift parameter s. The freedom in choosing this origin must be dealt with and the simplest functional method relies on computing the root of a scalar function, whose form is heuristically determined based on assumptions about the asymptotic structure of the eigenfunctions. The combined system is, to be solved for U s and s as a self-consistent system, for each chosen x s , with functional K constructed heuristically for = 1, 2, 3, in the next section.

B. Shift selection
Here we describe three heuristic arguments which lead to different forms of K , resulting in three different values of the shift s, and ultimately three different predictions for the critical amplitude, U s , per chosen extent of the perturbation x s .
The first heuristic seeks to minimize the amplitude U s chosen across all the possible choices of the shift, s. DefiningX(ξ; x s ) according to (4), so thatX(ξ; x s ) is normalized in the L ∞ -norm, rearranging (6) for the amplitude of the perturbation, leads to the maximization of the denominator (as the numerator is independent of s), and the resulting condition, with ∂ 2 s w 1 (ξ) X (ξ + s; x s ) < 0. Computing the derivative of this function with respect to s can be simplified using the definition of the inner product, which yields w 1 (ξ) X (ξ + s; x s ) = 0, so that we must find only the roots of a simple scalar equation, which when combined with (6) creates an appropriate solvability condition for U s .
The second heuristic assumes that the limiting feature of the linear theory is the magnitude of the perturbation to the critical wave,ĥ(ξ + s; x s , U s ), and a minimization of the perturbation norm in the L 2 sense, with respect to the shift, ensures the dynamics are appropriately linear and chooses the appropriate origin s. Expressing the perturbation to the wave in terms of the perturbation to the rest state and minimizing with respect to s, This reaffirms the importance of the Goldstone mode, and as with heuristic 1, the condition is combined with (6) to create an appropriate solvability condition for U s .
The third heuristic requires that a 2 (0) = 0, following the same formulation as (5), and is equivalent to the form used in Ref. [1]. As with heuristic 2, the condition is then combined with (6) to create an appropriate solvability condition for U s . The values of K l (ξ) and N l are summarized, Returning to (6), the solvability condition for U s is given by which determines the value of the shift s for a given choice of . We shall refer to the functions Φ as "shift selectors". The determination of s defines the solution of (6) for U s . Note that according to (8), scaling Φ by a nonzero constant factor, or even a factor that is a function which is finite and nonzero everywhere, does not change the answer. In the subsequent we shall silently use this property to simplify the expressions where convenient.
Here we list the explicit forms of the shift selectors dropping dependence on ξ, In the following we will use each heuristic-based shift selectors to predict the critical excitation curve for several slow-fast models.
Note that since v 2 = ∂ xû , the components required for computing these shift selectors, in addition to the critical pulse solutionû, are w 1 and w 2 . Also, the stucture of (10) means that v 2 and w 1,2 are required up to a nonzero constant scaling factor.
We look for the nonlinear wave solution as an expansion in , with the wave speedĉ = 0 + ĉ (1) + O 2 , and ξ = x −ĉt. Substituting into the traveling wave equation (2) and expanding in , we have in O 0 , so we have u (0) 1 = φ and the naive assumption is true. In where can be obtained from here in quadratures, provided that f 12 ≡ 0. We will not need the explicit expressions here, and return to these details when considering the degenerate case, characterized by f 12 ≡ 0. The linearization in the comoving frame (3) is similarly expanded, L = L 0 + L 1 + O 2 , with where f ijk = f ijk (ξ) = ∂ 2 f i /∂u j ∂u k evaluated at u = u(ξ, 0) = u (0) 1 (ξ),ū 2 , with the adjoint L † defined by the inner product w | Lv = L † w v . The eigenfunctions of the operators L, L † are also expanded for small , In the leading order, this gives by components 1,2 = 0, hence the naive assumption holds for v 1 but not w 1 . Note that we can similarly argue that v

So we have v
whereas the leading order for the left eigenfunction gives where the last (trivial) equation for w (0) 2,2 is understood to mean that any function satisfies the asymptotic eigenproblem at this stage, so long as f 12 (ξ) = 0. So again the naive assumption holds for v 2 but not w 2 .
For the sake of comparing the asymptotics with the numerics, we would like to know the asymptotic order of v 1,2 . The first-order correction v (1) 1 , using standard perturbation theory, is obtained as a linear combination of v (0) j for all j = 1. We have seen that inasmuch as i,2 = 0, and besides, v (0) 2,2 = 0 from symmetry considerations, hence we conclude that v (1) 1,2 = 0. Finally, the O 1 order for w 2 gives Assuming w (0) 2,2 (−∞) = 0, we find, up to a normalization constant, To summarize, the expected scaling of the key ingredients of the theory in the limit of → 0 is: The behaviour of these ingredients for the FitzHugh-Nagumo system obtained numerically is illustrated below in fig. 2, where we have used the empirically established scaling v 1,2 = O 2 . Taking into the account the structure of the initial perturbation given by (4), of practical importance are the "voltage" components of the shift selectors, Φ ,1 Using the definitions (10), we find where A and B are some constants; for reference, Observe that since φ(ξ) and χ(ξ) are even functions, we have that Φ 1,1 (ξ) and Φ 2,1 (ξ) are odd in the limit → 0, which guarantees the availability of the choice s = 0 for these selectors, as would be expected. At the same time, since A and B are typically both nonzero, lim →0 Φ 3,1 (ξ) is not odd, and the choice s = 0 is not available in this case. Though the limit 0 = lim →0 s exists for Φ 3,1 (ξ).

D. Asymptotic structure for the degenerate case
In the case of the Karma model and also for the cubic recovery variant of the FitzHugh-Nagumo model, the standard asymptotics described above do not work. More precisely, it fails for any model in which f 12 (u 1 , u 2 ) ≡ 0. To see why, let us consider in more detail the O 1 cor- (14) and the asymptotic boundary condition u 2 (ξ → +∞) → 0 it follows that This gives the leading order O 1 of the slow component of the critical pulse. The value ofĉ (1) can be obtained if we multiply (13) by u (0) 1 and integrate, by exploiting (11) and the boundary conditions u 1 (±∞)) = 0. In the degenerate case, f 12 (u 1 , u 2 ) ≡ 0, and according to (16) we haveĉ (1) = 0, and consequently no answer for u (1) 2 which hasĉ (1) in the denominator in (15). Therefore the asymptotics are to be determined separately, taking into account the specific dependence of f 1 on u 2 . We consider the dependence of the form where ω(u 2 ) = u 2 q and q > 0, so the the problem for the critical pulse is and postulate u 1 = u 2 , andĉ = γ νcĉ(1) to leading order. Substitution into the traveling wave equations yields the expected equation (11) for the critical nucleus solution, while the second equation relates two terms which must match to leading order in γ, from which we conclude that ν c + ν 2 = 1.
Considering the next-to-leading order in γ from (17) we have for which the balance to leading order in γ is achieved for ν 1 = ν c = qν 2 . Combined with the previous results for ν 2 and ν c , this gives Introducing = γ 1/(q+1) , we summarise that the nonlinear solution scales as While we recover ν 1 = ν 2 = ν c = 1/2 in the classical case (q = 1), we have more exotic scaling for different values of q.
We now determine the scaling of the solutions of the linearised problems, focussing on q > 1. The leading terms of the linearization operator in the comoving frame are and of its adjoint where the derivatives of the kinetic terms, including ϕ ω ≡ ∂ϕ/∂ω, are understood to be evaluated at the critical solution.
For the leading order for the first eigenpair we have , the latter being down to our arbitrary choice of normalization. In the leading order, the equation for the first component decouples from the second equation, . The second component is then to be obtained from the second equation, with the first component and the eigenvalue considered as given: Note that the v 1,2 -dependent terms on the left-hand side are asymptotically smaller than the right-hand side, hence the balance is achieved via so that the leading order contribution of v 1,2 is making self-constent our earlier assumption that v 1,2 does not exceed O( ). Curiously, we observe that the leading term in v 1,2 = O q+1 is smaller than the first- The second eigenpair is different in that we know σ 2 = 0 exactly and v 2 = u due to translational symmetry. Otherwise we proceed as before,

and the equation for the first component gives
in accordance with the leading asymptotic expansion of u 1 . Then the second equation gives Comparison of the first and second term here shows that v 2,2 is of a higher asymptotic order than v 2,1 . We therefore can neglect the third term in comparison with the second, which leads to For the first left eigenfunction, we have Hence we can take w , and then Finally, for the second left eigenfunction, we have and therefore To summarise, the eigenfunction components scale as Note that the generic case asymptotics are recovered by setting q = 1 and correspondingly = γ 1/2 , including the scaling of v 1,2 which in the generic case was not established conclusively.

III. METHODS
Throughout the remainder of this paper we shall deal with two-variable systems of partial differential equations of the form given in (1), distinguished by the details of the functional form of f . The essential ingredients of the linear theory for predicting the critical excitation strengthextent relationship remain the same across these models, however, and our methods for computing these are likewise similar. We begin by writing the system of partial differential equations in the frame moving with speed c yields a system of three ordinary differential equations, whose unique equilibrium is given by [ū, 0]. Continuing the rest state for increasing current forcing J connects to a Hopf bifurcation, from which a family of periodic orbits emanate. Continuing this family of periodic orbits to large periods with c = 0 followed by decreasing current forcing yields an unforced periodic orbit of the autonomous system, (20), equivalently a traveling wave solution of (1) with periodic boundary conditions on a domain x ∈ [0, L). To compute asymptotic traveling wave solutions of (1) we continue the unforced periodic orbit in the (L, c)-plane. The asymptotic solutions correspond to large domain sizes, for which c becomes constant: lim L→+∞ c (L) = 0. In this asymptotic regime, lim L→+∞ c(L) is multi-valued, and the lowest of the speeds, c =ĉ, designates the critical solution, (û,ĉ). In the large-L limit the periodic solution approximates the homoclinic originating from the rest state, and as a practical matter this limit is numerically inaccessible, particularly for small γ. In this context, we will approximate the homoclinic solution with projection boundary conditions [5], which permits aperiodic solutions by enforcing orthogonality to the stable/unstable eigenspaces of the rest state at the boundaries of the domain.
The periodic critical solution is computed on an adaptive collocation grid using Auto [6], at large γ and interpolated onto a Chebyshev grid of size M × DA representing M Chebyshev modes and a dealiasing factor DA ≥ 1. A nonlinear boundary value problem is constructed which corresponds to equations (20) and projection boundary conditions. The projection boundary conditions require the eigenvectors of the Jacobian evaluated at the rest state. The Jacobian of the ODE system is given by with unstable and stable subspaces spanned by the right eigenvectors satisfying E u , E s diag (Λ u , Λ s ) = J E u , E s . We require that the unstable pulse travels to the right (ĉ > 0) and thus that the perturbation from the rest state on the right hand side be orthogonal to the stable subspace (guaranteeing excitation dynamics), while the perturbation from the rest state on the left hand side be orthogonal to the unstable subspace (guaranteeing a relaxation to the rest state). The projectors of E u and E s are the corresponding left eigenvectors of J , so that P u , P s = E u , E s −1 , and the boundary conditions are where P u gives one condition at the left boundary (x = 0) and P s gives two conditions at the right boundary (x = L). The boundary value problem is solved using the Newton solver in the open-source Dedalus framework [7] until the update is smaller than the tolerance of The linearization utilizes the boundary value problem solution within the forcing terms, and is likewise discretized using M × DA grid points M Chebyshev modes and a dealiasing factor of DA ≥ 1. The projection boundary conditions for the linearization follow the same logic as presented for the nonlinear boundary value problem. The eigenproblem is solved by calling the sparse eigensolver package Arpack through scipy.linalg.eigs, with up to 4096 iterations retaining 64 Lanczos basis vectors to resolve the leading eigenmodes (significantly fewer than M ). Note that the adjoint linearization equations and adjoint boundary conditions are distinct from the forward linearization equations and boundary conditions; in particular, while the forward linearization maintains one boundary condition at the left boundary and two conditions at the right boundary, the adjoint linearization problem applies two conditions at the left boundary and one condition at the right boundary. The Jacobian for the forward eigenproblem is the same as the nonlinear problem Jacobian, (21), and thus uses the same boundary conditions applied to the first-order form of v i . The Jacobian for the adjoint eigenproblem is given by , and the corresponding boundary conditions are Since the forward and adjoint linearized eigenproblems are formulated independently, the eigenvalues (σ j ) resulting from the calculation of the left eigenfunctions (w j ) and the eigenvalues (σ i ) resulting from the calculation of the right eigenfunctions (v i ) are compared and matched pairwise, and the left and right sets of eigenfunctions are used to verify the biorthogonality conditions, The prediction of the critical curve is conveniently done using the Fourier transform, so the solutions are sampled on a uniform grid of sufficiently large size N , for the nonlinear as well as the linear problems. It begins by forming the shift selector Φ , explicitly computing the inner products and appropriate sums using the normalized eigenfunctions, such that w j | v i = δ ij . To determine the shift value for a prescribed perturbation shape,X(x), two cross-correlation integrals are computed. First, the perturbation is cross-correlated with Φ (x) using the product of the Fourier transforms, ],X(q) = F X (ξ) , and the Fourier transform and its inverse are defined aŝ The roots of µ are computed by checking for sequential differences in the sign of the elements of µ , and refined using a Newton method applied to a locally adapted spline interpolant of µ . The linear prediction for the critical strength as a function of the shift (7) is also computed using cross-correlation via Fourier transform, and selecting the global minimum of |U s (s)| > 0 for refinement by applying a scalar Newton method to a function which computes µ (s) = 0 to sub-grid accuracy by a cubic interpolant. The root of µ (s) determines the position of the perturbation at time t = 0, The linear theory prediction is compared to direct numerical simulations (DNS), using initial conditions corresponding to the perturbation shape, u(0, x) =ū + U sX (x). The problem is solved on the domain x ∈ [0, l] where the domain size is no larger than half the asymptotic domain size l ≤ L/2, and over a time interval t ∈ [0, T ] where T is comparable toč/l whereč is the speed of the fastest isolated asymptotic wave. The problem is discretized in space using a Chebyshev spectral method using M modes, and time-stepped using a third-order implicit-explicit method (RK443 [8]), as implemented in the numerical package Dedalus [7]. The boundary conditions are chosen as ∂ x u 1 (t, 0) = ∂ x u 1 (t, l) = 0, and enforced using the Chebyshev-tau method [9].
For an initial condition parameterized by U s , a wave may ignite and recruit the entirety of the tissue (U s > U sc ), or it may immediately decay and approach the rest state (U s < U sc ). To numerically distinguish between these events it is necessary to track the state as it evolves over time, and in particular for the tracking to unambiguously characterize these outcomes. To this end we define a distance function, which compares the amplitude of the initial condition in the voltage channel to the amplitude of the state in the voltage channel at all later times. Considering only the final value ψ(T ) gives an effective scalar function, f : R → R, and for sufficiently smooth flows f is continuous. For initial conditions set to the critical wave solution, u(0, x) =û(x), it is clear that ψ(t) = 0, as it is indeed for all initial conditions which are equilibria of the underlying partial differential equations (1). For generic initial conditions parameterized by U s , and considering an interval U s ∈ [U s , U s ] where f (U s ) < 0 and f (U s ) > 0, we can compute the roots of f using a iterative bisection procedure. Additionally parameterizing the perturbation to the rest state by the width, x s , then for each sampled width x s the range of U s is determined, and the bisecting procedure is applied to determine the pair (x s , U sc ) which defines the critical perturbation.

IV. RESULTS
Throughout this section we compare DNS results to the predictions computed using the linear theory, and pay particular attention to the variation in the critical amplitude predictions for different choices of shift-selecting heuristics, as well as the adequacy of the predictions as we approach the singular limit of these slow-fast systems parametrically. We also compare the critical amplitude to the unstable root of the kinetics, f 1 (u 1 ,ū 2 ) = 0, which is the asymptotic threshold for the double limit γ → 0, x s → ∞ for the class of perturbation shapes used in this work. We perform these comparisons for a small but representative set of nonlinear excitation models with two variables, with temporal and structural differences but falling under the slow-fast paradigm.

A. FitzHugh-Nagumo
FitzHugh-Nagumo is a prototypical model of excitation following the asymptotic reduction of the Hodgkin-Huxley model equations for the giant squid axon [2]. The FitzHugh-Nagumo kinetics for f are given below, with fixed α = 0.37, β = 0.131655. The speed ratio varied in the interval 10 −10 ≤ γ ≤ 10 −2 . The precise choice of β is such that it equals to its critical value β c = γ c + 2 √ αγ c , which corresponds to the transition between stable node and stable focus in the local kinetics, for γ = γ c = 10 −2 . The first step is to compute the unstable wave solution, u(x), and then its linear eigenspectrum. We repeat these calculations for varied γ, from critical γ = γ c down to very small value to test the singular limit. Note that the singular limit of traveling waves in the classical FitzHugh-Nagumo system is well-studied [10][11][12][13]. The scaling results of the previous section apply directly to this model, and thus serve as a method of veri-fying the numerical results under variations in the modal expansion length M , domain length L, and the application of projection boundary conditions, across several decades in γ. In particular, we find the expected asymptotic scaling of the wave components and pulse speed with γ for L ≥ 50 and M ≥ 32, with exponentially small corrections to the leading eigenvalues for M ≥ 256 when L ≥ 100. For sufficiently small γ, the solution of the marginal eigenvalue σ 2 = 0 is eventually corrupted by the proximity of the essential spectrum of the wave to this eigenvalue, and the eigenfunctions v 2 and w 2 may not be computed using these methods; while M and L dependent, we have found this to occur for γ < 10 −10 , placing a limit on the reliability of the numerical approach taken in this work. Figure 1 demonstrates the correct scaling of the pulse speedĉ with γ, with good agreement with the γ 1/2 behavior for small γ. Figure 2 summarizes the scaling of the components of the linear theory across six decades of γ, additionally showing the fast component of the constructed Φ . Note that since the perturbation is along the fast component, only the fast components of the shift selector, Φ ,1 , are relevant. Despite the correct observed scaling of the critical pulseû and leading eigenfunctions v 1 , v 2 , w 1 , w 2 , the construction of Φ indicates a subtlety in the selection of an optimal frame for the prediction of critical perturbations. While Φ 1 and Φ 2 converge to antisymmetric functions centered about the maximum ofû 1 , Φ 3 instead converges to an asymmetric function whose central root remains offset from ξ = ξ 0 in the limit γ → 0. Considering the components of Φ 3 as defined in (10), we recognize that the contribution proportional to w 1 does not monotonically vanish as γ → 0, as originally expected in the formulation of the linear theory [14] which assumed identity of the γ → 0 limit with the γ = 0 "critical pulse" case.
The deviation in the central root of Φ 3 has implications for the selection of optimal reference frames and the critical perturbation amplitudes associated with them. Figure 3 shows the critical strengths over three decades of the parameterized width x s , both predicted using Φ and computed using direct numerical simulation. The figure also shows the selection of the optimal frame shift s for each of the predictions, in particular the systematic offset of s 3 as x s → 0, which coincides with the central root of Φ 3 . This systematic offset in the selection of the optimal frame becomes more severe as γ → 0, indicating that the prediction due to Φ 3 does not converge (this is not evident from the figure as presented but confirmed by careful observation of raw data). We conclude that as far as Φ 3 is concerned, the slow-fast system with γ → 0 and the single-variable fast system dynamics of u 1 in isolation are qualitatively different, i.e., this is a singular limit. Further, it indicates that predictions made at large γ are in fact more accurate than those made for smaller γ, which is reflected in the figure. Notably, for all observed values of γ, Φ 3 fails to recover the large (x s → ∞) perturbation limit. A superficial interpretation of this small paradox is that in the γ → 0 limit, all the events that decide the fate of a particular perturbation happen at the time scale O(1), so the slow variable remains almost at its resting value, whereas the heuristic behind Φ 3 is based on the matching of the initial condition against the critical pulse in full, i.e. taking into account both components. In the critical pulse, the slow variable is different from the resting value; although this difference is small, the sensitivity of the pulse position to the perturbation in the slow component, measured by w 2,2 , is on the contrary large in this limit, hence the overall contribution of the slow component does not vanish in the limit γ → 0. A small, but nonetheless important note is that the predictive power of the linear theory is sensitive to the monotonicity and the root structure of Φ for some perturbation widths, to the extent that the dependence of predicted s and U s on x s is discontinuous. Such discontinuity occurs, for example, when a branch of µ (s) = 0, while continuing in x s , terminates in a fold. In the vicinity of the fold value of x s , the solution (s, U s ) on one side of the fold may differ significantly from the solution determined on the other side of the fold.
This appears in the linear theory prediction utilizing Φ 3 , i.e. fig. 3, for γ = 10 −2 . For x s 25 the linear theory makes no reasonable predictions for U s , and for 5 x s 25 the prediction for U s diverges from the asymptotic value correctly predicted by Φ 1 and Φ 2 . The mechanical reason for this deviation is that the correlation integral defining µ corresponds to a box filter or smoothing operation, and for sufficiently large x s this smoothing is destructive. This smoothing destroys the delicate root structure seen in Φ 3 for γ = 10 −2 ; so that while for x s → 0 the correlation integral forms a translation operation (preserving the number of roots), when the box filter width x s is comparable to the distance between successive roots of a function, the result of the correlation integral may have fewer roots, and specifically, it may only have roots which produce anomalously large predictions for U s .

B. Mitchell-Schaeffer
The Mitchell-Schaeffer [4] model is a popular semiconceptual model of cardiac cells, combining the simplicity of only two components with the relatively realistic description of action potential shape and restitution properties. Historically, it has been derived via an asymptotic reduction (adiabatic elimination of the fastest processes) of the more detailed Fenton-Karma model of atrial excitation, but still incorporates multiple decay timescales for the system. We re-write the kinetics of the model in the form where θ(u 1 ) = H k (u 1 − u g ) = {1 + tanh [k(u 1 − u g )]} /2 is a smoothed Heaviside distribution centered at u 1 = u g with width k −1 , and the timescales ratio in terms of the original parameters is γ = (τ i /τ c ) in this rescaling. The standard parameter values are τ i = 0.3 ms, τ o = 120 ms, τ u = 6 ms, τ c = 150 ms, u g = 0.03 and k = 100. In the following examples, the parameter ratios τ i /τ u = 0.05 and τ c /τ o = 1.25 are kept at standard values, while γ is treated as a free parameter. That is, for the purpose of the asymptotic theory, all of τ i /τ u , u g and k −1 are treated as finite even though they are "small" in layman's terms. The time t is dimensionless as presented, likewise we absorb the original diffusion coefficient in the nondimensional spatial scale, x. The scaling of the critical pulse and associated leading eigenfunctions in the Mitchell-Shaeffer model follow the expected scaling, see fig. 4 and fig. 5, similar to the FitzHugh-Nagumo model, while reproducing more realistic action potentials and gate switching dynamics. However, while the localization of the critical pulse in the fast variable (u 1 ) is not dissimilar to the critical pulse in FitzHugh-Nagumo, the relative scale of the components of the leading eigenfunctions are reversed. Namely, we note that in relative terms, the second components of v 1 and v 2 are two orders of magnitude larger than would be expected based on their asymptotic order in γ alone. This is related to the above mentioned non-asymptotic small parameters in the model, specifically, sharp switching of the f 2 across u 1 = u g .
The shape of Φ 1 and Φ 2 are standard -both are approximately antisymmetric about the peak of the wave, see the bottom row in fig. 5. However, the contribution in the term proportional to w 1 in Φ 3 is again significant, so we should expect a similar offset root and likewise inaccurate predictions for this shift selector.
The predicted critical excitation amplitudes U sc track the DNS results over three orders of magnitude for the extent of the perturbation, for the value computed using shift selectors Φ 1 and Φ 2 . As expected, the predictions using the shifts determined by Φ 3 are much less accurate: they are systematically larger than the DNS results by nearly an order of magnitude. Note that this error in U sc is caused by a deviation of s by only a fraction of the critical pulse width off the predictions of the other two shift selectors. Comparison of the predicted results for large γ indicates that it is the proximity of the central root of Φ to the position of the peak of the unstable pulse which determines the effectiveness of the prediction, c.f., fig. 5. However, careful examination of fig. 6 indicates that between γ = 10 −4 and γ = 10 −2 the primacy of the different shift selectors changes; that is, at γ = 10 −4 , we see that Φ 2 is the most accurate shift selector, while at γ = 10 −2 it is Φ 1 though no shift selector generates satisfactorily accurate predictions.
Taking note of fig. 6 for γ = 10 −2 , we observe that Φ 3 predicts critical excitation strengths which fail for x s > 11, i.e., for some sufficiently wide perturbations the branch of µ 3 (s) = 0 ends in a fold, just as with the FHN model results. This mechanism thus appears to be a generic feature of qualitatively similar inputs to the infrastructure of the linear theory. It remains to determine the conditions under which Φ 3,1 changes shape sufficiently that the convolution with the perturbation in the formation of µ 3 destroys the optimal root. In principle, we may consider a toy model of the dependence Φ 3,1 (s; γ), in the form ϕ(s; a) = aw 1,1 (s) + √ 1 − a 2 w 2,1 (s), whereby ϕ(s; 1) has no roots, while ϕ(s; 0) has a single root which persists under convolution with the perturbationX with a fixed, large, width x s . At some intermediate value of a = a c , the number of roots in the convolution changes, which may be distinct from the value of a = a c at which the number of roots within a central region of the domain of ϕ( · ; a) changes. The disparity between a c and a c suggests an analogous liminal region in which the time-scale separation γ yields Φ 3,1 with an appropriate root, but no corresponding root in µ 3 .
Indeed, fig. 6 suggests that for a given x s , we may select a frame which precisely reflects the DNS results, that is, the actual position of the critical nucleus observed as a long transient when the perturbation magnitude is at its closest to the critical value. We say that this makes a postdictive optimal frame through the determination of s DNS , effectively. The relevant value of γ = 10 −2 presents a situation in which all three predictive curves are distinct and each badly represents the DNS results. Iterating through perturbation extents we find that s DNS > s 2 > s 1 > s 3 over three decades of x s . In the limit of x s → 0, lim s DNS ≈ +3.47, while lim s 1 ≈ −5.94, lim s 2 ≈ 0.00, and lim s 3 ≈ −26.4. This exercise informs about the neighborhood of the predictive measures, however extending this observation to a general principle, i.e. generating a Φ DNS , is far from obvious.

C. Modified "cubic recovery" FitzHugh-Nagumo
Here we consider a modification of the FitzHugh-Nagumo model (22) which is motivated by the nonlinear dependence of f 1 on u 2 in the Karma model considered in the next subsection. Our modified FHN model is Obviously, (22) corresponds to the case q = 1. In this subsection, we shall look at q = 3 instead. We refer to this model as "FHN with cubic recovery" or FHNCR for short. We keep β = 0.131655 as in (22), but take α = 1.2 in order to keep the kinetics excitable, i.e. have only one equilibrium, see fig. 7(a).
Due to the degenerate dependence of f 1 on u 2 , this model has different asymptotic properties, discussed in Subsection II D. As expected from the results of that subsection, the pulse speed scales as γ 3/4 , giving relatively fast convergence of the pulse to a static nucleus state, see fig. 7(b). Meanwhile, u 2 ∼ γ 1/4 , see fig. 8 (top row) -a much slower convergence than standard FHN (q = 1). Recall the pulse solutionû asymptotically converges to the critical nucleus solution u = [φ, 0] as γ → 0. Similarly, slower convergence is observed for components v 2,1 , w 2,1 , and the component w 1,1 asymptotically vanishes unlike its counterpart in FHN model. This is all in agreement with the predictions from Subsection II D. We are of course mostly interested in how this difference affects the accuracy of the lineary theory predictions for the critical curves.
The FHNCR model results should be read in dialogue with the FHN (q = 1) model results. As a point of concrete comparison, consider γ = 10 −8 and x s 10 in the FHN model, fig. 3, where the frame selector Φ 3 predicts an amplitude which is approximately twice as large as the DNS result. Compare to the same configuration for the FHNCR model results for γ = 10 −7 and γ = 10 −9 , fig. 9, which predicts a critical amplitude which is an order of magnitude larger than the DNS result, and five times larger than the linear model prediction. While the frame shifts selected by the Φ 3 condition prove significantly worse for FHNCR than for linear FHN, the frames selected by Φ 1 and Φ 2 are perfectly adequate predictors for the critical amplitude for both models. That is, the predictive power of the linear theory is not affected by the relative scaling of the components so long as a frame selecting heuristic is chosen carefully. Further, numerical experiments with the FHNCR (not shown) suggest that, as q increases further, we should expect the predictive power of Φ 1 and Φ 2 to outstrip that of Φ 3 more generally. Some slow-fast models of cardiac excitation, in particular, form a highly nonlinear (and indeed, parametric) dependence of f 1 on u 2 . We consider one such model in the next subsection.

D. Karma
The Karma-1994 [3] model is a qualitative model of cardiac excitation, designed to reproduce chosen restitution curves. The Karma model kinetics are given by the following functions, where we will consider the parameter set u * = 1.5415, M = 4, u d = 3.0, β = 1.389, u a = 0.5, τ 1 = 2.5 ms, τ 2 = 250 ms, and in terms of the original notations of [3] we define γ ≡ τ 1 /τ 2 as the ratio of timescales, so that t is dimensionless. The original diffusion coefficient of the model is dimensional, we absorb this quantity in the non-dimensionalization of the spatial scale, x. Figure 10(a) sketches a phase portrait of the kinetics and fig. 10(b) shows the scaling of the speed of the unstable pulse solution for the Karma model. The pulse speed scales as γ 4/5 , which matches the predictions of Subsection II D; we of course have q = M here. Figure 11 summarizes the unstable pulse solution and leading eigenfunctions. As for the other models discussed above, the amplitude of u 2 decreases as γ → 0; however, for Karma the convergence is significantly slower than the FHN or MS convergence rates. The convergence of the slow components of the leading left eigenfunctions is likewise complicated by the quartic dependence on u 2 . In particular, the convergence rate for w 1,2 and for w 2,2 are markedly different, suggesting that for some intermedi- ate value of γ the dominant contribution to Φ 3,1 switches from the term proportional to w 2,1 to the term proportional to w 2,2 , and that predictions made on one side of the scale of γ will contradict predictions made on the other, or asymptotically. A simple calculation suggests that this occurs for very large γ, when the asymptotic argument no longer holds. While it is generically true that the structure of the leading eigenfunctions of the unstable pulse changes as γ → 0, for the generic models considered previously the localization of the eigenfunction components are convergent for sufficiently small γ. For the FHNCR and Karma models, the leading right eigenfunctions (v 1 , v 2 ) and the leading left eigenfunction (w 1 ) are well-behaved in the limit of small γ. Both the leading right eigenfunctions (v 1 , v 2 ) are localized to the same region as the pulse, i.e., the leading right eigenfunctions inherit their localization from the nonlinear solution. Likewise, w 1 is localized like v 1 and the first components of these eigenfunctions coincide in the limit of small γ. However, as γ → 0, the slow component of w 2 , w 2,2 , behaves qualitatively differently. While for the generic models, w 2,2 converges to a localized function which decays quickly outside of the central region of the unstable pulse, for the degenerate models it does not. As γ → 0, w 2,2 delocalizes asymmetrically, such that |w 2,2 (ξ)| > 0 ahead of the excited region of u 1 , while |w 2,2 (ξ)| → 0 behind the peak.
The Karma model results represent a stress-test for the asymptotic scaling, both in terms of the stiffness of the kinetics and the parametric nonlinear dependence of f 1 on u 2 . The former specifically marks the asymptotic structure of the pulse as an approximately piecewise linear curve, and the small-γ structure of the leading right eigenfunctions as wildly deviant from the other, smoother, models described in this work. The latter presents an opportunity to test the efficacy of the asymptotic analysis for severe nonlinearity.
In addition to the quartic nonlinearity (M = 4) results presented here, we also computed the asymptotic scaling of the linear theory ingredients for M = 2, M = 8, M = 16, and M = 32. In each instance of the parameterized model, the scaling of each component follows the predicted asymptotics described in Subsection II D. As with the comparison of the cubic and linear FHN model critical curve predictions, as the nonlinearity M increases, the less accurate Φ 3 becomes, while Φ 1 and Φ 2 maintain their predictive power across several decades of x s and γ.
V. DISCUSSION a. The mathematical problem addressed here is that of conditions required for initiation of a propagating excitation wave by an instant perturbation from the resting state by a stimulus of a certain spatial extent. The gist of our approach, previously exposed in [1,[14][15][16][17][18], is in linearization of the PDE system around a critical solution. A delicate issue is translational invariance of the problem which generates a one-parametric family of critical solutions, and poses a problem of identification of the member of that family that corresponds to a given initial perturbation. In the previous works, this issue was addressed by an heuristic suggesting that the initial condition of the linearized problem should not contain the shift mode. In this framework, the essential ingredients of the linearized theory, apart from the critical solution itself, are the right and left eigenfunctions of the linearization, corresponding to the first two eigenvalues, the first positive eigenvalue responsible for the instability of the critical solution, and the second being the zero eigenvalue corresponding to the translational symmetry. In particular, the left eigenfunction corresponding to the zero eigenvalue serves as the projector onto the shift mode. In this paper, the focus is on four selected twocomponent systems with slow-fast time scale separation, like in FitzHugh-Nagumo model (FHN), including FHN model itself. In this class of models, the critical solution is a slow unstable propagating pulse. The results presented here were originally thought of as no more than further tests of applicability of the above mentioned approach in a particular class of models. This aim has been broadly achieved, however certain unexpected aspects have been revealed, which may serve as valuable lessons both for the specific problem but also at large for the theory of slow-fast systems.
b. The expectation was to verify an asymptotic theory [14] based on the would-be obvious assumption that in the asymptotic limit, the events in the fast subsystem dominate, and therefore the predictions of the linearized theory should converge to those for the one-component system, in which the slow variable is frozen at the resting state value. Specifically, it was expected that the critical pulse solution converge to the stationary unstable nonuniform "critical nucleus" solution of the fast subsystem, and the eigenfunctions correspondingly converge to those of the critical nucleus as far as the fast components are concerned, and the slow components be negligible in the asymptotic limit.
c. Lesson one: this expectation has proved wrong. Although the critical pulse solution does indeed converge as expected, and the fast components of the eigenfunctions indeed converge as expected, but the slow component paradoxically does not become negligible. Specifically, the slow component of the second left eigenfunctions, which is the projector to the translational mode and is therefore material for the heuristic used for critical pulse selection, grows large in the slow-fast asymptotic limit. This eigenfunction corresponds to sensitivity of the speed of the propagating pulse to perturbations of the slow component ahead of it, and increases and spreads out in sync with slow-down of that component. As a result, the overall contribution of the slow component in the overlap integral, coming from the left and right eigenfunction, does not vanish in the asymptotic limit, and although the predictions based on ignoring the slow components work reasonably well, as it appeared in the analysis done in [14], the non-vanishing contribution from the slow component creates a systematic error in identifying the critical pulse, which considerably spoils prediction and in some cases renders the theory inappicable in principle. d. Lesson two: there are other heuristics which do withstand the asymptotic limit. Heuristics are required in our approach because of its leading idea to use linearization in the situation where the solution in question is in fact not small. The old heuristic, which has been in use since [16], was to make the linearization "more applicable", by making sure that its initial condition is "as small as possible" in the sense that in the critical situation, i.e. at the margin between successful and unsuccessful initiation, not only the first generalized Fourier component vanishes, which is a condition of criticality, but the second component vanishes, too, which is always achievable by an appropriate translation of the critial solution with respect to the initiation stimulus. There are of course may other senses in which the initial condition can be made "as small as possible", offering alternative heuristics. In here we have explored two of them, one that minimizes the initial condition of the linearized problem in L 2 norm, and the other that minimizes the predicted threshold given by the criticality condition. Both new heuristics do not depend on the second left eigenfunction, and both have shown expectable convergence in the asymptotic limit, i.e. to the critical nucleus results, and good predictive ability, i.e. correspondence with the direct numerical simulations.
e. Lesson three: asymptotics may not give good predictions even when the small parameter is large. This of course can happen if the small parameter in question, i.e. the ratio of time scales of the fast and slow component, is not the only small parameter in the problem, but there are others, such as other time scale ratios, or sharpness of transitions in the reaction kinetics, as we have seen in the examples of Mitchell-Schaeffer and Karma-1994 models. In such cases, the range of applicability of the asymptotics depends on those other parameters, and may be well away from realistic parameter ranges.
f. Lesson four: not all asymptotics are the same. The natural genericity assumptions about the dependence of the kinetics terms on the dynamic variables may fail, leading to completely different asymptotics in the slow-fast limits. A priori this possibility might seem remote, but the fact that we have stumbled on such a failure "by accident" in a popular, even if simplified, cardiac excitation model [3] suggests that this possibility should be kept in mind. The slow rate of convergence in the fast/slow time scale separation parameter, particularly together with other small parameters present in the model, may render the fast/slow asymptotics irrelevant, in the sense that the behaviour of the solution at the original parameter value may be rather far from the asymptotic one, even though the original parameter value appear rather small.
g. Further directions. Straightforward extension of this study would be to slow-fast systems with more than one fast and/or more than one slow component with similar, Tikhonov type occurrence of the small parameter. A more intriguing possibility is about asymptotics in case of non-Tikhonov slow-fast models, of the kind discussed in [19][20][21]. The difference from systems with the Tikhonov structure is that the asymptotic limit of the critical solution is not a stationary "nucleus", but a moving front. Applicability of the linearized theory has been tested on a conceptual model of such critical front in [1,16]. However, convergence of the ingredients of the linearized theory and of corresponding predictions for the critical curves has not been explored so far to our knowledge. Given the lessons from the FitzHugh-Nagumo type systems discussed above, one should not take such convergence for granted. This direction is particularly important since non-Tikhonov asymptotics have been argued to better represent the properties of realistic ionic models of cardiac excitation, than FitzHugh-Nagumo type systems, particularly at the margins of propagation.
Implications of the observations presented here are also relevant for exotic solutions of one dimensional excitable models. We have noted that the convergence of the slow component of w 2 is different for the generic and degenerate models, not only in the scaling of the solution near the asymptotic, but the asymptotic shape itself for small γ. To reiterate, w 2,2 is asymmetrically extended for x ahead of the pulse peak, while decays to zero quickly behind the pulse peak. This feature of w 2 in conjunction with the localization of w 1 suggests that the slow dynamics (the physically relevant dynamics) of the unstable pulse are nearly insensitive to perturbations positioned post-peak, but very sensitive to being slowed by perturbations almost arbitrarily ahead of the pulse peak. This may play an important role in the develop-ment of "back-initiation", or the observability of a "onedimensional spiral" solution generally [22]. The extension of w 2,2 should lead to acceleration of newly created pulse formations ignited by back-initation and increase the potential for local collapse to the resting state. As we know such an unstable solution exists in FHN-type models, one would expect that the increased sensitivity of an extended w 2,2 may suppress the formation of these dynamics, suggesting that degenerate models may have more complex saddle structures.
One natural extension of the existing program is to the "critical quenching problem", that is, of cessation of stable wave propagation in an excitable medium by addition of minimally invasive perturbations. The application to quenching is an inversion of the application to ignition, though the central ingredients can be the same and rely on the same linearization about the unstable pulse state. Crucially, as quenching considers an equivariant state in the form of the stable wave, as compared to the invariant quiescent state, the problem involves the consideration of an additional parameter which fixes the additional translational symmetry. This problem will be addressed in a forthcoming paper.