Periodically kicked feedforward chains of simple excitable FitzHugh-Nagumo neurons

This article communicates results on regular depolarization cascades in periodically-kicked feedforward chains of excitable two-dimensional FitzHugh-Nagumo systems driven by sufficiently strong excitatory forcing at the front node. The study documents a parameter exploration by way of changes to the forcing period, upon which the dynamics undergoes a transition from simple depolarization to more complex behavior, including the emergence of mixed-mode oscillations. Both rigorous studies and careful numerical observations are presented. In particular, we provide rigorous proofs for existence and stability of periodic traveling waves of depolarization, as well as the existence and propagation of a simple mixed-mode oscillation that features depolarization and refraction in alternating fashion. Detailed numerical investigation reveals a mechanism for the emergence of complex mixed-mode oscillations featuring a potentially high number of large amplitude voltage spikes interspersed by an occasional small amplitude reset that fails to cross threshold. Further careful numerical investigation provides insights into the propagation of this complex phenomenology in the downstream, where we see an effective filtration property of the network; the latter amounts to a successive reduction in the complexity of mixed-mode oscillations down the chain.


Introduction
In mathematical neuroscience, the study of nerve depolarization goes back to the original Hodgkin-Huxley equations for electrical propagation in the squid giant axon [1]. That model features a nonlinear reaction-diffusion system that accounts for the ionic fluxes across the nerve cell membrane as well as the cell axon's properties as an electrical cable. A variety of models have since been derived from the Hodgkin-Huxley equations. Among them, the reduction proposed by FitzHugh and Nagumo [2,3] combines the effects of several auxiliary ion channels into a single recovery variable and thus leads to a simpler dynamical system that is more amenable to mathematical analysis; the resulting PDE has been analyzed extensively both theoretically and in applied context [4][5][6][7][8][9][10][11].
Some recent studies have focused on a (space-)non-homogeneous reactiondiffusion version of the FitzHugh-Nagumo model [12,13]; this work was done with particular emphasis on the interplay between the oscillatory and excitatory properties of the traditional two-dimensional ODE system. This interplay provides a setting wherein localized oscillations (realized by tuning a specific part of the medium to the oscillatory regime) diffuse into the surrounding excitable medium. The studies examined a class of Neumann boundary condition problems for the PDE (RD FHN) A crucial feature of this model was the incorporation of a spacial nonhomogeneity c (= c(x)) in the equation of the recovery variable v. This feature afforded a mechanism for the generation of a periodic signal in some small region of the spacial domain, which may in turn propagate or fail to do so depending on the excitability properties (which vary according to c(x)) of neighboring regions. From a mathematical perspective, the success or failure of the propagation of oscillations corresponds to an interesting bifurcation phenomenon in this infinite-dimensional setting. Similar properties have been exploited in a more applied context [14] where the authors have used a FitzHugh-Nagumo reaction-diffusion system with a non-homogeneous Laplacian term to model axon damage resulting from traumatic brain injury. The present article provides and analyzes a related simple model for wave propagation phenomena in spatially-discrete unidirectional chains of excitable nerve cells. Each cell in the chain has voltage dynamics modeled by a two-dimensional FitzHugh-Nagumo (FHN) ordinary differential equation. Continuous dependence in space is replaced by a simple coupling scheme through instantaneous kicks to the recovery variable. More precisely, with reference to fig. 1 we consider a feedforward chain of FitzHugh-Nagumo neurons {S j } j≥1 modeled by where t i = αi for some fixed period α, and ε is small but positive. The tsubscripts denote ordinary differentiation with respect to time. This chain is spatially homogeneous, in the sense that the function f as well as the parameters ε, c, A, and K are all j-independent. In the second set of equations, we employ the convention ∞ · 0 = 0. As in the standard FHN model, the function f is assumed to have a graph resembling the general shape of a negativeleading-coefficient cubic; for explicit calculations and numerical simulations, we assume f (u) = 3u − u 3 . The parameter A represents the size of Dirac stimuli presented in periodic fashion to the chain's front node S 1 , or else arriving at S j for j > 1 when the voltage variable u j−1 of the preceding cell S j−1 crosses (on upswing) a certain threshold value K (we also take K = 0 for numerical simulations); in either case, these stimuli are incorporated into the dynamics through the recovery variable v j . A detailed investigation of eq. (FHN FFC) is motivated in part by the necessity for a discrete-space alternative to eq. (RD FHN) that features comparable phenomenology by way of significantly simpler dynamical mechanisms. The most natural analog of eq. (RD FHN) studied in [12,13] involves a discretized Laplacian term in the right-hand side of the u equation as well as a dependence of c on the index / location of the cell within the network. This correspondence yields a coupled map lattice whose complexity is at least as high as that of the original partial differential equation (see [15], esp. relevant chapters [16][17][18]); accordingly, an analysis of the corresponding coupled map lattice would be just as much, if not more complicated than an analysis of the original reaction-diffusion model. While it is typical for neural coupling to enter the cell dynamics directly through the voltage variable [19,20], eq. (FHN FFC) features the well-studied dynamics of the isolated FitzHugh-Nagumo system as perturbed only occasionally by instantaneous kicks delivered to the recovery variable. The ensuing dynamical system is readily amenable to a mathematical treatment -we will discuss both rigorous results and detailed numerical work, all of which point to a variety of interesting and biophysically relevant bifurcation phenomena that can be established analytically. This simplification notwithstanding, the choice of coupling for our model is still appropriate since it presents a second-order cumulative effect to the voltage while serving to emphasize the threshold barrier that is a common feature of depolarization in the presence of excitability. (Specifically, for the example of the cubic nullcline featured in the FitzHugh-Nagumo system, the recovery variable's clearance of the local minimal point effectively induces a rapid jump in the voltage variable -this serves precisely the purpose of guaranteeing a mechanism for the generation of action potentials via excitation.) We note also that this simple coupling formalism is related to common modeling paradigms of collective network behavior in mathematical neuroscience. For instance, forcing the network's front node by an external input current, together with Dirac coupling between cells through the recovery variable is commensurate with, e.g., the usual integrate-and-fire coupling scheme where communication between cells is achieved by discrete kicks/jumps in the synaptic current variable prior to being incorporated into the voltage equation. Thus, the proposed model eq. (FHN FFC) features a simple coupling scheme analogous to those commonly employed in the formulation of spiking networks [21][22][23].

Further notations, definitions, and basic properties of FHN
The isolated FHN model S • with A ≡ 0 features a rest state at (c, f (c)) and is an example of a slow-fast system. The u nullcline v = f (u) is also the system's critical manifold, wherein trajectories are constrained for ε = 0. For 0 < ε 1, away from the critical points for the graph of f , singular perturbation theory guarantees the existence of a flow-invariant curve ε-close to this nullcline, where trajectories follow the slow timescale; such a curve is referred to as a slow manifold for the system. In this regime, we have that u t is positive below the u nullcline and negative above it; accordingly, the flow above a slow manifold (say, above and sufficiently separated from v = f (u)) is expected to move to the left (negative u direction) according to the fast timescale, and one has precisely the oppositely-directed fast motion below.
We focus on the parameter regime c < −1, for which it is well-known that the system is excitable [20,24]. Starting from a certain point (c, v 0 ), with v 0 sufficiently smaller than f (c), one can clearly construct a conjunction of four singular orbits consisting of σ 1 : a horizontal (fast) trajectory to the right from the point in question to the point (u 0 , v 0 ) satisfying f (u 0 ) = v 0 ; σ 2 : a critical trajectory within the curve v = f (u) from (u 0 , v 0 ) to the right critical point of f , say (u c , v c ); σ 3 : a horizontal (fast) trajectory to the left from the right critical point of f to the other point in v = f (u) satisfying v = v c ; σ 4 : a critical trajectory within the curve v = f (u) connecting the last point (u l , v c ) to the equilibrium (c, f (c)).
Accordingly, for ε ≈ 0, singular perturbation theory [25][26][27] assures the existence of a orbit γ • (t) of the isolated FHN starting from (c, v 0 ) that is asymptotic to (c, f (c)) and ε-close to the concatenation σ 1 ∪ σ 2 ∪ σ 3 ∪ σ 4 described above. For the purposes of this article, excitability is to be interpreted in the following sense: the system is at rest at (c, f (c)), but if presented with a sufficiently large instantaneous momentum transfer in the negative v direction, the system crosses threshold to the right (this constitutes an effective depolarization), and so resets to the rest state only after this threshold crossing; alternatively, if the kick is not sufficiently large, reset to the rest state happens directly without a depolarization. Strong excitation with discharge is thus possible, but not necessarily achieved.
Of note, the phase portrait for an isolated FHN cell contains a rich geometry that provides a mechanism for mixed-mode oscillations (MMOs). Recall that MMOs correspond to switching between small amplitude and large amplitude oscillations. Such patterns were first discovered in the Belousov-Zhabotinsky reaction [28,29] and have since been frequently observed in experiments and models of chemical and biological rhythms, particularly in the context of neuroscience [20,30,31]. There is a large number of works dealing with MMOs -both the theoretical and the applied aspects. We refer to [32] for a general review of various results on MMOs as well as [33] for a short summary. A remarkable way to obtain MMOs is through the canard phenomenon found in slow-fast systems. Canard solutions are phase space trajectories that follow a repulsive manifold for a certain amount of time. They were first discovered in the Van der Pol oscillator [34] with non-standard analysis techniques, and this led to numerous subsequent works, for example [27,[35][36][37]. In this context, small oscillations may correspond to canards that are repelled to a nearby attractive manifold, while large oscillations arise from canards which are repelled to a distant one. Combined with a return mechanism, this provides a setting for sustained MMOs [32,38,39].
The model we investigate also features the emergence of MMOs within an individual cell as well as their propagation through the feedforward chain. To draw a distinction between our case and the examples from the literature referenced above, we note that the first of these phenomena is based on an interplay between the impulsive periodic forcing and a phase space geometry that allows both relaxation oscillations and simple reset to equilibrium. Variations to the period of the forcing can lead to multiple consecutive relaxation oscillations (see fig. 20) that are occasionally followed by a reset to equilibrium. Some of the observed trajectories in this model are related to the canard phenomenon since they seem to follow the repulsive manifold for a non-trivial amount of time. Overall the dynamical mode is characterized by the property that a solution curve consists of multiple loops about the equilibrium (c, f (c)), some of which are sufficiently large to cross the depolarization threshold at u = 0, while others remain too small for this. This amounts to an evolution of the voltage variable u that features a combination of small-amplitude flutters and large-amplitude excursions, the latter resulting in the triggering of synaptic firing toward the neighboring cell. Since this complex phenomenology provides a mechanism for skipping some of the beats in a regular depolarization train, it also establishes a way for cells in this network to support a richer spike train encoding structure. We discuss this regime in some detail in section 4, and we explore the downstream propagation of MMOs in section 5.
We denote by ϕ(·) = (ϕ j (·)) N j=1 the unique solution to the version of eq. (FHN FFC) featuring N cells evolving dynamically from a prescribed initial value ϕ(0) = (ϕ · (0)) ∈ (R 2 ) N . We will be particularly interested in the initial prescription ϕ j (0) = (c, f (c)) for all j, corresponding to the scenario where the chain is initially kicked from its rest state. In accordance with standard terminology, a solution ϕ to such an initial value problem is said to be periodic of period τ if for all t ≥ 0. We refer to ϕ as a traveling wave solution (TW) provided that there exist a function ψ : [0, ∞) → mathbbR 2 and a j-independent time shift β ∈ R, such that (ϕ(t)) j = ϕ j (t) = ψ(t − j β), j ∈ 1, . . . , N , for all t ≥ 0. We note that, in the context of the treatment presented herein, a periodic TW in eq. (FHN FFC) corresponds to the presence of rhythmic and regular cascades of depolarization in the chain. Related phenomenology has been observed in the context of type I phase oscillator networks [40,41] under excitatory coupling, where a variety of rigorous and numerical results point to the existence and stability (both Lyapunov and structural) of traveling wave solutions in homogeneous feedforward chains. In addition to its connections with simple spiking network models, our model provides an excellent opportunity to address questions concerning wave propagation in the context of excitable cells, which differ from oscillatory neurons because of the attraction to a resting equilibrium in the absence of stimulus. Studying side-by-side oscillatory and excitatory models may reveal the extent to which the presence of stable regular waves of activity is a stimulus-driven phenomenon. If we suppose that cell S 1 receives a kick (by −A in the v-direction) sufficiently large to induce a simple fast trajectory across threshold resulting in depolarization, the second cell S 2 in the chain is also affected by the stimulus. Thus results a simple setting for exploring depolarization propagation in chains of excitable units coupled according to the equations above. This relates to the work [12,13] where traveling pulses are generated from a central generator and diffuse outwards, according to the FitzHugh-Nagumo reaction-diffusion model eq. (RD FHN). As in that prior work, of interest to the present study -will periodic excitation at S 1 propagate down the chain or be blocked? If it propagates, will that propagation be akin to a periodic wave with a well-defined velocity (like an effective diffusion in a homogeneous medium), or can the propagation be more complex?

Plan for the paper
The paper is organized as follows. In section 2 we give an overview of our results, with emphasis to distinguish rigorous work from deductions based solely on thorough numerical simulations; these results are summarized in fig. 2. In section 3, we establish a theoretical result guaranteeing the presence of a Lyapunov-stable traveling wave solution with periodicity that matches the forcing period. In section 4 we provide a rigorous result establishing a mechanism for especially simple mixed-mode oscillations at S 1 under periodic forcing with kicks. We also present numerical investigations of the variety of possibilities for more complex mixed-mode oscillations at S 1 . Then, in section 5.1, we relate our investigations of mixed-mode oscillation phenomenology at S 1 to downstream propagation, which can be significantly more complicated than the simple scenario established in section 3. We conclude with a survey of numerical experiments on propagation when dynamics at one or more of the upstream sites is complex and poorly understood; these numerics are described in section 5.2 and seem to indicate that complex behavior may be successively averaged / filtered out by the first few nodes of a chain, so as to obtain tame behavior downstream -such as a traveling wave -while the first few upstream sites behave in a more complicated manner. Section 6 provides some concluding remarks as well as directions for future work.

Low frequencies (rigorous)
If the frequency of incoming kicks is sufficiently low, S 1 is guaranteed to settle to within a small neighborhood of (c, f (c)) prior to any subsequent stimulation. For the purposes of this discussion, we treat this as having the forcing period α sufficiently large (say, there is some α 0 > 0 and the preceding remark holds Table 1 Values for the parameters used in the numerical simulations for α > α 0 ). Since we assume the kick size A to be sufficiently large for the vertical offset induced by a kick to clear the neighborhood of the critical manifold where the system can feature complex behavior (specifically, we have in mind avoiding canard trajectories; see, e.g., [36]), the resulting dynamics at S 1 will be time asymptotic to a periodic threshold traversal with period α, which guarantees S 2 receives a periodic sequence of kicks identical to that presented to S 1 . The argument outlined above clearly extends inductively because of the homogeneity in the network. The result is an asymptotic trend toward a solution close to an α-periodic traveling wave (TW) down the chain, a phenomenon resembling that studied for phase oscillator ensembles in [41] and [40]. Of note, there are some theoretical subtleties to be detailed during proofs. Nevertheless, the network supports traveling waves of arbitrary (sufficiently high) period, when the stimulus / kick size is sufficiently large in the sense presented above. We establish this rigorously in section 3.

Higher frequencies (analysis and numerics)
As the forcing period α is decreased, a complicated interplay between stimulus size and frequency begins to emerge. This interplay provides a setting for very diverse phenomenology that includes a variety of mixed-mode oscillations (MMOs) as well as obstruction to depolarization. Thinking of the stimulus/kick size as fixed going forward, exploring the interplay amounts to a parameter study by way of decreasing α. 1 See figs. 13 to 19.
Setting aside the observation that a full half-open interval (α 0 , ∞) of forcing periods features the tame behavior discussed above, we continue the α-parameter exploration. For moderately higher frequency, kicks are overwhelmingly likely to arrive while the dynamics is in a slow regime (and yet, before it reaches a small neighborhood of the equilibrium). There are two possibilities -a kick arrives either on the downswing toward the stable equilibrium (c, f (c)), or on the upswing slow branch to the right. Since the kicks always decrease the value of the recovery variable v, the former possibility can yield a fast subsequent trajectory that either crosses threshold for depolarization (provided that kick arrives sufficiently late) or wavers for a small amplitude oscillation around (c, f (c)), close to the left extremum of the critical manifold (see figs. 13, 14, 16 and 18); on the other hand the latter possibility yields various interesting complications, including a regime where the neuron is trapped forcing period α low forcing frequencies in the downstream Fig. 2 A schematic of the results discussed in this paper. Rigorous analysis at S 1 for forcing with α-periodic Dirac kicks reveals α-periodic response (at low frequency, with α > α 0 ) and 2α-periodic simple mixed-mode oscillation (at higher frequency, with α ∈ (α 2 , α 1 )); in turn, this analysis extends inductively down the chain, thus establishing the existence of α-periodic TW at low frequency (resp., 2α-periodic TW for an interval of higher frequencies). A detailed numerical investigation explores a bifurcation cascade as the period α is increased across the threshold α 1 and then dialed up to the threshold α 0 ; this analysis reveals successively more complex mixed-mode oscillations at S 1 , whose complexity decreases as they propagate downstream -that is, the chain filters out the complexity and decreases the number of large amplitudes before a small amplitude break (when examining the dynamics at S j , as j increases). A few numerical case studies of this filtering phenomenon (at S j , with j low) are also presented for illustration in the very high forcing frequency domain, where α < α 2 .
in small amplitude oscillations on the right (upswing) slow branch and never resets to the left slow branch, which drives the voltage u to +∞ and prevents it from undergoing future depolarizations.
At some point, provided that α is within an appropriate range, one simple mechanism for MMO at S 1 emerges as follows. For simplicity, assume the cell starts at its natural equilibrium (c, f (c)). The kick received to start is sufficiently large to induce a fast horizontal trajectory (assume f (c) − A is sufficiently smaller than f (−1) = −2) across threshold; thereafter, the subsequent stimulus received at t = α comes when the trajectory is in the refractory phase, but not yet sufficiently close to the equilibrium to guarantee clearing the threshold anew. As a consequence, the subsequent stimulus induces a much smaller oscillation consisting of a direct attraction to the equilibrium instead of a large oscillation. Thereafter, over the next α interval, the cell has sufficient time to reach the aforementioned neighborhood of the equilibrium that is essential to assuring that it clears the firing threshold at the subsequent stimulation; see figs. 9 to 11. In section 4 we prove that in this regime, the dynamics of S 1 features a 2α-periodic solution with this geometry. For definiteness, we define the range of α where these simple MMOs are supported as (α 2 , α 1 ). Accordingly, when α is between α 1 and α 0 , more complex mixedmode oscillations are present; see section 4 for details. On the other hand, the MMO property is absent when α is only slightly below α 2 . 2 The simple MMO at S 1 implies that the stimulus presented to S 2 amounts to a periodic sequence of kicks with period 2α. Of course, if 2α > α 0 , the overall effect is attraction to a solution close to a 2α-periodic traveling wave, with period equal to twice the forcing period for all S j with j ≥ 2 (again, this will be shown rigorously) -and this is the first instance in which we have a complete understanding of phenomenology beyond the simplest traveling wave for low frequency. But it is also possible that 2α < α 0 ; and in the latter case, we note there may yet be a 4α-periodic oscillation of simple mixed-mode type at S 2 , or a vastly more complicated situation; in the former case, the prior reasoning extends inductively; in the latter, we illustrate typically observed behavior through a few representative numerical simulations. In particular, these latter simulations speak to the manner in which MMOs are transmitted down the feedforward chain. Our numerics indicate that when they arise, complex MMOs induce less complex behavior (sometimes of MMO type, not always) down the chain; we investigate this attenuation further in section 5.2.

Existence of simple α-periodic TW
We begin with a detailed study of the behavior of S 1 when A > 0, assuming that the system begins at initial condition (c, f (c)). If it were to receive a single stimulus, this neuron would instantaneously find itself at (c, v 0 ), with v 0 = f (c) − A. Clearly, if A is sufficiently large, the system will feature a large trajectory orbit (including a depolarization event) that is asymptotic to (c, f (c)). Since is assumed to be small, our proofs are often inspired by the case = 0. In this case, the trajectory after kick is no more continuous, but evolves first on the right part of the critical manifold, and jumps to the left part when u reaches the value 1. There are however some differences between the case = 0 and > 0. For example, if = 0 one can prove existence and uniqueness of the periodic solution at S 1 . But the argument of uniqueness fails for > 0. Therefore, we consider below the cases = 0 and > 0 separately.
2 See our discussion of the state of affairs for α only slightly smaller than α2 in remark 10. The limit case = 0 Proposition 1. We assume = 0 and A > 0 fixed sufficiently large. If α is sufficiently large, S 1 admits a unique α-periodic trajectory which crosses the firing threshold u = K. Furthermore, this trajectory is stable.
Proof Existence Let F denote the discrete-time map that takes an initial condition on the left branch of the critical manifold v = f (u) to its image under kick + time-α flow. Notice that the curve under discussion can be parametrized by the v-coordinates of its points (1-1 correspondence). To prove existence of a periodic solution, it is sufficient to prove that F admits a fixed point. First note that F (f (c)) > f (c). This results from the fact that (c, f (c) is a stable stationary solution of the ODE. Now, we search for a point along the left branch of the critical manifold (above (c, f (c))) having the property that the trajectory ensuing from it returns with a smaller vcoordinate at time t = α. Such a point exists for sufficiently large α due to the attraction toward (c, f (c)). We denotev its v-coordinate. Then     F (v * ) = v * which corresponds to a periodic solution.
Uniqueness We still denote v * the v-coordinate of the kick-point on the periodic solution exhibited above. Uniqueness follows from the observation that all trajectories ensuing from points with v-coordinates between f (c) and the v-coordinate v * are mapped (after kick+time-α flow) to points with v-coordinates above v * . This results from the fact that any trajectory ensued below v * will spend more time on the right branch of the critical manifold than the trajectory ensued from v * . A symmetric argument can be applied from above, and together these observations ensure that the periodic solution is the only trajectory whose pre-kick v-coordinate remains invariant under the discrete dynamics. Local Stability We denote by (u 1 , v 1 )(t) the α-periodic solution. Let (u 2 , v 2 )(t) be another solution with initial conditions satisfying u 2 (0) > u 1 (0), u 2 (0) − u 1 (0) small. And let -u i (0 + ), i ∈ {1, 2} denote the values of u i on the right part of the critical manifold right after the kick; -t 1 denote the time at which (u 1 , v 1 )(t) reaches the fold point (1, f (1)) = (1, 2); -τ denote the time needed for u 2 to go from u 2 (0 + ) to u 1 (0 + ).
As a consequence, and also −u Remark 2. Figure 3 gives a representation of F for α = 10. It shows numerical evidence that |F (v * )| < 1 which implies local stability. Note that F can be computed almost explicitly. The calculation entails computing the times spent on the cubic -first in the right branch, then in the left one after the jump -and retrieving the value of v at time α; this leads to solving the equation where f −1 relates either to the right or left branch of the cubic.

The case > 0 (small)
We now proceed to the case ε > 0. The proof of the existence of a periodic solution is analogous to the case = 0. The proof of uniqueness does not extend to the case > 0. However, we show local stability, which in turn implies that the periodic solution is isolated. Stability is obtained by specific computations. Let F denote the mapping which maps an initial condition (u, v)(0)) to (u, v)(α) corresponding to the point on the trajectory after kick + time α under the flow.
Proof The existence of a periodic solution follows from the Brouwer theorem. Consider a small enough closed ball B around (c, f (c)); for α large enough, global asymptotic stability of (c, f (c)) for the ODE, ensures F (B) ⊂ B. We now consider the stability. We denote (u 1 , v 1 )(t) the α-periodic solution. Let (u 2 , v 2 )(t) be another solution with initial conditions close to (u 1 , v 1 )(0) . Let We have Now, we remark that u is controlled by the initial conditions, and the trajectories u 1 and u 2 are well known. Furthermore, a simple computation shows that the time for which f (u 1 (s)) > 0 is of O( ), while the time for which f (u 1 (s)) < 0 is uniformly greater than a positive constant. Altogether, assuming u(0) small enough leads to the conclusion that α 0 u 2 f (u 1 ) − 3u 1 u − u 2 dt < 0, which means that g(α) < g(0) and implies the result. Remark 4. The reader might think to apply the classical results for the stability of limit-cycles as in [42] p 216 (or [43] for the result with proof ). However, looking into the details shows that the jump brings some technical difficulties with respect to the local change of variables around the periodic solution, which prevent a direct application of these results. The general condition ensuring stability would write α 0 f (u 1 (s))ds < 0 which is finally very close to the argument provided here.
Remark 5. As mentioned above, the proof for uniqueness of the periodic solution provided for = 0 does not extend to the case ε > 0. The point is that when > 0 there is some time spent in the fast trajectories. A trajectory starting below another goes faster along the fast fibers but spend more time in the slow manifolds. It follows that for very close trajectories, arguments provided for = 0 do not apply. We will illustrate how this phenomenon may lead to the emergence of MMOs.
The existence of a travelling wave solution is stated now. Corollary 6. Assume the hypothesis are as in proposition 3. Then {S j } j≥1 has a α-periodic traveling wave solution.
Proof Since the chain is homogeneous in i, if every node starts in an initial condition corresponding to the periodic solution for S 1 , the trajectories will be the same for each i, with a delay between the node i and the node i + 1 corresponding to the time needed for the node i to reach the semi-axis u = 0, v < 0.
The next proposition deals with the stability of the traveling wave solution. Let (ϕ j (t)), 1 ≤ j ≤ N the traveling wave solution exhibited in the corollary 6. For any initial condition of the system eq. (FHN FFC), for each j ∈ {1, .., N }, we denote by (t n j ) n∈N the sequence of times at which the node S j receives a kick. The following proposition holds We conclude this discussion with a heuristic calculation of the velocity (and related quantities) for the traveling wave. Notice that in this simplest case, one can deduce an explicit estimate on the speed of the resultant α-periodic traveling wave, or equivalently, the time of propagation for pulses presented to the front node. For any fixed cell, the time lag between stimulus arrival and consequent excitation of its right neighbor is approximately equal to the amount of time that it takes for a trajectory to travel from u = c to u = K. Accordingly, we obtain Taking the additional simplifying assumption of starting from rest allows one to replace v with f (c) − A, and this yields the formula In addition, notice that the time needed to excite the whole chain is N × T e . We emphasize that this is the amount of time taken by a single pulse to travel through chain front to back. Given the integral units of the spacial coordinate, all of this allows us to deduce that the speed of the traveling wave is 1/T e . And lastly, in accordance with the definition of TW detailed in section 1.1, one can take the parameter β to be any element of T e + αZ; without loss of generality, we think of β as the smallest positive element of this set -namely, T e .

Mixed-mode oscillations at S 1
It should be noted that for α 1 < α < α 0 , prior to the simple MMO regime described in table 1 and established herein, continuous dependence on parameters implies a highly-complex transition from the simple α-periodic trajectory to complicated oscillations; recall fig. 2. For instance, this can be realized via canard trajectories that either follow the unstable branch of the critical manifold for some time prior to making a fast jump (a path that possibly includes threshold crossing) or otherwise canard trajectories that loop directly back toward the equilibrium without crossing threshold. It is equally worth noting that a possibly very complicated combination of small and large oscillations can ensue; see fig. 14. A detailed analysis of the bifurcation structure in this system will be the subject of subsequent investigations. For the present discussion, we establish the existence of simple MMO rigorously and then conclude this section with a few numerical investigations of the bifurcation cascade.
Proposition 8. Assume the condition on A as in the hypothesis of proposition 1. There is an interval of periods (α 2 , α 1 ) such that for each α in this interval, the dynamics of S 1 feature a 2α-periodic trajectory with one large loop that crosses threshold at u = 0 and one small simple loop that does not. (See fig. 9; this is the simplest mixed-mode oscillation scenario.) Remark 9. Altogether, this scenario results in an overall effect of doubling the depolarization period for the first cell. In this case, if 2α > α 0 , it follows   This trajectory then turns to the left and is simply attracted to the equilibrium (as t → ∞). In fact, there is a portion C 0 of the critical manifold consisting of points with v-coordinates slightly higher than the v-coordinate of the origination point of this special trajectory, all of which feature precisely this simple dynamics. Suppose B δ is a ball about (c, f (c)). It follows not only that all trajectories starting on C 0 will end up in B δ , but also that the time it takes before they enter (and never leave thereafter) is less than the time it takes for ψ(t) to enter it. This is because ψ(0) is the lowest (v-coordinate) point on C 0 .
We now deduce some necessary explicit time estimates for trajectories of the system. The time that it takes for a typical trajectory to change its v-coordinate from v 1 to v 2 can be computed by using the dynamics on the critical manifold and appealing to the formula in tandem with the change of variables v = f (u). Explicit calculations are straightforward for the example f (u) = 3u − u 3 and yield This formula is used in the calculations captured within fig. 12. We argue for ε = 0. Assume α is large enough but not too large, so that a trajectory that ensues from (c, f (c)) ends up on C 0 at time t = α. It follows that there is a range of α such that for all δ sufficiently small, the ball B δ about (c, f (c)) is mapped inside C 0 after one kick and time α flow of the ODE. We close by noting that there exists a range of δ for which the time needed for ψ(t) to enter B δ after one kick plus ODE flow is much smaller than smallest of the α in the aforementioned interval. (See fig. 12.) Altogether, this means that there is a closed ball B δ about (c, f (c)) that is mapped within itself by the sequence of two applications of kick plus time-α flow. By application of the Brouwer fixed point theorem, we are guaranteed the existence of a fixed point to the discrete time map consisting of the conjunction of kicks and flows just described. This property then assures the existence of the aforementioned 2α-periodic solution to the dynamics in S 1 under the α-periodic forcing.

Remark 10.
When α is only slightly smaller than α 2 , the observed phenomenology at S 1 is straightforward despite not being of MMO type. The arrival of the first kick is indeed too soon on the left branch of the critical manifold -in fact, so soon that the trajectory is not kicked below the equilibrium point (c, f (c)) on the graph of f ; accordingly, the solution undergoes a simple fast trajectory to the right and then continues down the left branch of the critical manifold until the second kick; the latter then causes a vertical drop sufficient for the trajectory to clear the critical manifold and cross threshold.
Just like in the case of simple MMO, this period-2α solution has immediate consequence for S 2 , which evidently sees a lower-frequency input. In fact, the chain {S j } for j ≥ 2 will entrain to the 2α-periodic TW that is guaranteed for that forcing period by the analysis in section 3, because the numerical values in table 2 indicate 2α > α 0 when α ≈ α 2 . Of course, for other choices of the system parameters (i.e., different from those declared in table 1), more complicated propagation scenarios might occur. In the interval between α 1 and α 0 , the system can feature longer sequences of regular depolarizations interspersed by occasional quiescent loops. This is an interesting and mathematically complicated parameter regime. We share here some numerical observations together with our proposed geometric mechanism for this behavior.
Remark 11. One of the remarkable phenomena occurring in the regime α 1 ≤ α ≤ α 0 is that, after each large-amplitude oscillation, the values of v at which the trajectory for S 1 receives a kick form an increasing sequence until attraction to the stable equilibrium causes a small oscillation. This gives rise to an effective periodicity with MMOs, that is, a rhythmic train consisting of sequences of consecutive large-amplitude oscillations interrupted by one small one.
Proposed mechanism. The reader should refer to figs. 16, 18 and 20. With reference especially to the third figure here, there is a region below the critical Fig. 13 Steady-state evolution of ϕ 1 = (u 1 , v 1 ) as a function of time. α = 8.3. Here α is in the interval (α 1 , α 0 ). After a kick and a first big loop the solution is not attracted to the stable point; instead, it goes for another large oscillation. However, the time needed to reach the right portion of the slow manifold is longer on this second loop because of the velocity dissipation when the trajectory passes close to the local minimal point (-1,-2) of the graph of f . When the next kick arrives, the solution is further up along the left branch of the slow manifold and consequently enters a neighborhood that is directly attracted to the stable point; in turn, this gives rise to an oscillation with a much smaller amplitude. As a result, the steady-state here appears to be a periodic (or nearly so) trajectory with two large and one small oscillations. manifold v = f (u) in which trajectories (these go mostly left-to-right) take progressively longer to go across from the left part of phase space over to the right branch of the manifold; this tendency agrees with the vertical ordering of these trajectories, with the higher ones taking a longer time. This phenomenon is most drastically observed when trajectories of the system track the critical manifold beyond the critical point at u = −1. In particular, given two such trajectories -one below the other -the left-to-right crossing time for the higher one is larger than the combination of crossing and ascent times for the     lower one, which is to say that if the system were to travel (say, from initial condition u(0) = −1) along these trajectories, it would take longer for it to reach the right jump point on the higher trajectory than on the lower one. Now, let us assume that the system is kicked from rest at (c, f (c)). After a large oscillation, provided that the next kick arrives when the trajectory is in the downward (slow) swing, the new (post-kick) v coordinate of the trajectory is bounded below (which is to say, above) v = f (c) − A. At this point, either the trajectory features a small oscillation with attraction to (c, f (c)), or it features another large oscillation (but one whose lower portion is at a level above the first one). It seems then that the MMO phenomenology observed in, e.g., fig. 18, is owed to the fact that when not attracted directly to (c, f (c)), the trajectory finds itself precisely in the region of phase space previously described. A sequence of large oscillations ensues, each one featuring a left-toright crossing above the prior, until the the kick arrives sufficiently early to induce a small loop with attraction to (c, f (c)). In effect, this last step renews the system and the process starts over. Fig. 20 A phase plot of a single trajectory of (u 1 , v 1 ) for α = 8.45, specifically, its successive passes after successive kicks (shown in various colors, in order) progressively closer to the left extremal point of the critical manifold. This is due to a speed dissipation effect that grows with proximity to the minimal point of the graph of f . This figure suggests the following mechanism: after six large oscillations (with depolarization) corresponding to six kicks, at the seventh kick, the trajectory is trapped toward the attractive equilibrium with a small oscillation. The time afterwards needed to reach a small neighborhood of the steady state is much shorter than α; consequently, at the next kick the trajectory is not distinguishable from the one following the first kick. Hence, this mechanism gives rise to a 6-large-1-small MMO.
We couple this observation with a description of the phenomenology at S 1 as the forcing period is decreased from α 0 to α 1 : In principle, there can be a very high number of large loops prior to a small resetting oscillation. This number is larger for α ≈ α 0 and steadily decreases, as α decreases to α 1 (see figs. 13, 14, 16 and 18). Notice that, as noted in proposition 8, for α in the interval (α 2 , α 1 ), there is just one large loop paired to a small one.

Remark 12.
Understanding the extent of this complex phenomenology constitutes an involved mathematical problem related to various prior work and deserves a specially dedicated investigation beyond the scope of our descriptive exposition here. A few examples of mathematical analysis related to the present problem are provided in the context of three-time-scale equations, see [44][45][46] and also [47]. In there, as in here, solution trajectories constrained to follow a canard (and thereby result in relaxation oscillations) will instead diverge from this expected path, thereafter tracking in the opposite direction. This behavior is captured nicely by the computable solutions of a 2d-integrable system. 5 Case studies of propagation and filtering of mixed-mode oscillations beyond S 1

A simple mechanism for filtration
Assume first that A is sufficiently large for proposition 8 to apply at S 1 ; furthermore, assume that α ∈ (α 2 , α 1 ). The dynamics of S 1 is then asymptotically attracted to a simple mixed-mode 2α-periodic oscillation that features depolarization only once in its period. It is immediate that S 2 receives forcing at half the frequency of the original input (the period is doubled to 2α). Thereafter, if 2α happens to exceed α 0 , the tail of the chain will support the corresponding stable 2α-periodic TW as opposed to a TW of period α. This scenario -both observed in our numerics and established rigorously -presents a simple mechanism for the failure of wave propagation, a property observed in a variety of models in neuroscience. We refer to this phenomenology as filtering, since it presents an effective downshift in spiking frequency from the perspective of signal processing. It is also worth noting that the presence of higher-complexity MMO for certain forcing periods (e.g., those between α 1 and α 0 ) provides a setting where this sort of failure of wave propagation (equivalently, wave mutation) can happen in a much more complicated manner. For instance, if 2α is in the doubling interval (α 2 , α 1 ), the dynamics at S 2 will be asymptotic to a 4αperiodic trajectory; in turn, S 3 sees a sequence of Dirac kicks with period 4α. On the other hand, in principle one could have 2α ∈ (α 1 , α 0 ), whereby S 2 supports an MMO with some complexity, no longer featuring uniform spacing between consecutive kicks, but instead skipping kicks every once in a while. Of course, some of these possibilities may not be achievable in practice given the model parameters fixed in our numerical investigations. Fig. 21 Traveling wave of period 2α for {S j } j≥2 . Represented by (u j ) 1≤j≤100 at time t = 100. α = 8. Notice the kick between j=1 and 2 which will not generate a traveling wave but rather die at j = 2 (compare with fig. 8); this is due to the fact that the wave profile starts at j = 2.

A few examples indicating even more-complex downstream filtering of mixed-mode oscillations
A cursory juxtaposition of figs. 16 and 22 reveals that S 2 presents a simpler MMO than S 1 at steady-state with α = 8.41 ∈ (α 1 , α 0 ). This is to say that the complexity of MMOs decreases down the chain. A similar observation follows from examination of figs. 24 and 25, albeit in the regime α < α 2 which corresponds to the very high forcing frequency setting. The voltage dynamics in time allows another perspective on the filtration properties of the network. Figure 26 illustrates an α-to-2α-to-4α doubling cascade through cells S 1 , S 2 , and S 3 for the instance in which α < α 2 , followed by 2α ∈ (α 2 , α 1 ) (notice that proposition 8 applies here), whereafter 4α > α 0 . Starting with the third cell, the tail of the chain {S j } j≥3 will feature a regular simple firing pattern with a 4α-periodic TW profile.
On the other hand, higher complexity of the MMOs at an upstream node (the situation where 2 j α is finds itself in (α 1 , α 0 ) during the course of the induction) has further interesting consequences. This is illustrated in fig. 27,   where the third cell in the chain sees a burst-like behavior with four cycles per period, three of them 'on' beats, and one of them a quiescent 'off' cycle. Of note, the second cell behaves almost like this, but features a genuine MMO with a small-amplitude oscillation instead of the quiet phase.

Conclusion / future work
We have investigated the response of a FitzHugh-Nagumo cell to periodic forcing with instantaneous Dirac impulse, when the latter is incorporated into the recovery variable. Furthermore, we have provided detailed studies on the propagation of such a response in the context of a spatially extended feedforward network of cells communicating via excitatory kicks. The results of our parameter exploration (in decreasing forcing period) point to transitions from simple regular depolarization cascades to various mixed-mode oscillations and their propagation / filtration through the network. This work points to increasing complexity as forcing frequency is increased, in competition with a Fig. 26 α = 4. Left: Voltage dynamics of S 1 , simple periodic oscillation with period 2α (twice the forcing due to the fact that every second stimulus arrives early, during the slow portion of the trajectory, so the dynamics continues toward the stable equilibrium after a brief rightward fast motion following that kick). Center: Voltage dynamics of S 2 , which receives periodic stimuli at period 2α ∈ (α 2 , α 1 ); this induces a simple MMO with one largeamplitude oscillation and one small-amplitude reset. The overall effective period is 4α Right: Voltage dynamics of S 3 shows a simple 4α-periodic oscillation with one depolarization in each of its periods. Left: Voltage dynamics of S 1 , simple periodic oscillation with period 2α (twice the forcing due to the fact that every second stimulus arrives early, during the slow portion of the trajectory, so the dynamics continues toward the stable equilibrium after a brief rightward fast motion following that kick). Center: Voltage dynamics of S 2 , which receives periodic stimuli at period 2α ∈ (α 1 , α 0 ); this induces a complex MMO with three large-amplitude oscillations and one small-amplitude reset. The overall effective period is 8α. The important point of emphasis here is that while the stimulus this cell presents to S 3 is periodic, it is not uniform. (Effectively, there are three consecutive beats and one small flutter phase in each period.) Right: Voltage dynamics of S 3 shows a simpler 8α-periodic oscillation with three depolarizations in each of its periods, interspersed by a quiescent phase distinct from the mixed-mode behavior in the second cell, in that the voltage ramps up instead of fluttering. This can be understood by noting the decomposition of the presented stimulus, which features three beats at 2α together with a fourth flutter that does not produce a stimulus. So instead of making a small recovery (flutter loop), S 3 in fact has time to recover directly to a neighborhood of the equilibrium prior to the process renewing. For reference, we have omitted the plot for S 4 but it is identical to the one for S 3 ; this suggests that the situation settles as j increases. Notice that the chain dynamics appears to converge to a TW that is periodic but qualitatively distinct from the simple regular TWs discussed in section 3. The latter TW features voltage dynamics at downstream sites that mix consecutive beats with single breaks in between, so they are not uniform.
strong regularization or damping effect as the excitation progresses through the medium.
Of note, worthwhile future directions for this research include the investigation of direct stimulation to the voltage variable (for comparison with work done on eq. (RD FHN)), as well as more detailed exploration of the correspondence between stimulus size, period of forcing, and steady-state values in the parameter regime where the natural equilibrium is a sink of spiral type (which may afford the possibility for mixed-mode oscillations with higher number of small-amplitude voltage resets). Another interesting avenue of exploration would be to formulate a probabilistic analysis of the corresponding model under stochastic forcing with, e.g., a Poisson spike train. The interplay between the randomness of kick arrival times, modulo well defined average frequency, and the slow-fast geometry in the FHN dynamics poses an interesting mathematical challenge; it also brings the model's applicability to the realm of random stimuli, an altogether different area of applications.
Lastly, our studies have revealed a particularly interesting, if somewhat pure, mathematical challenge. The possible complexity of MMOs present in a single FHN cell driven by regular Dirac kicks is evident in the numerical investigations. Our admittedly heuristic explanation for this phenomenon relies on the presence of canard trajectories that veer toward the unstable branch of the critical manifold despite starting out as "fast" trajectories. Within this realm, one can attempt to address the question of whether it is possible for the steady state solution to feature an arbitrarily high number of long (thresholdcrossing) canard trajectories prior to the small reset; that is, does this number increase without bound under adjustments to the period? or, is it bounded? Conversely, this question speaks to whether the length of time in between small-amplitude breaks can be arbitrarily large.