EMERGENCE OF TRAVELLING WAVES IN SMOOTH NERVE FIBRES

. An approximate analytical solution characterizing initial condi- tions leading to action potential ﬁring in smooth nerve ﬁbres is determined, using the bistable equation. In the ﬁrst place, we present a non-trivial sta- tionary solution wave. Then, we extract the main features of this solution to obtain a frontier condition between the initiation of the travelling waves and a decay to the resting state. This frontier corresponds to a separatrix in the projected dynamics diagram depending on the width and the amplitude of the stationary wave.

1. Introduction. A major problem in mathematical neuroscience and applied science is to understand and control the excitability of a nerve fiber so that external stimulations induce the firing of action potentials or a decay to the resting state. Depending on the power of an external stimulus, a fiber can indeed stay at its resting state or its membrane can be charged so that it induces a local change in the transmembrane potential. This external stimulation, which can be realized by a small local electrode, can be seen as an initial condition corresponding to a spatially distributed transmembrane voltage along the nerve fiber. In 1937, Rushton introduced a concept of the liminal length. According to this hypothesis, there is a minimum length of a fiber which must be raised above the excitation threshold in order to initiate the propagating wavefront [1]. The emergence of travelling waves is not only function of the amplitude but also of the spatial distribution of the stimulation. In order to quantify the effects of these parameters, we use the bistable equation (corresponding to the FitzHugh-Nagumo (FHN) model [2] without recovery variable) which is widely used to investigate the mechanism of the propagation [3,4,5] in nerve fibres. McKean and Moll [6], using this model with a piecewise linear representation of the sodium current flowing through the membrane ((f (V ) in eq.(1)), demonstrated that for the initial data corresponding to a local stimulation with a single electrode, there exists a threshold surface in the space of initial data that separates subcritical initial conditions, which decay to zero, from the supercritical initial conditions, which expand into a pair of propagating wavefronts, sometimes called trigger waves. Aronson and Weinberger [7] used a maximum principle to prove that, for the bistable equation, sub-and super-critical initial data were bounded by a threshold surface. The qualitative characterization of this threshold surface has been performed by  [8] by using the gradient flow structure of the bistable equation. More recently, Neu et al. [9] proposed to rewrite this problem with a variational approach by introducing a gradient flow of the energy defined for a steady-state solution of the system. Using a projection of this gradient onto an approximate solution space, this method leads to obtain the quantitative conditions on the width and amplitude of initial conditions which enable the emergence of a travelling wave. They focused on an asymptotic case where the threshold parameter corresponding to the activation of the sodium current is very small (α 1 in eq.(2)). Using the same methodology, we propose an approximate determination of the conditions needed to fire a nerve membrane whilst introducing this threshold parameter α as a detuning parameter. This is motivated by the fact that this parameter can take a wide range of values usually not very small, depending on the kind of neurons. In this article, we consider the bistable equation which includes the stable and the unstable excited states due to a cubic polynomial function as a description of the sodium current. In section 2, the analytical shape of the stationary wave is determined, corresponding to initial conditions. In section 3, a parametric form of the stationary solution is used to define an analytical shape of the separatrix in the projected dynamics diagram. This separatrix corresponds to a stable manifold and provides the conditions which allow the emergence of travelling waves. In section 4, we propose a determination of the stable manifold, using the stable manifold theorem and the Poincaré's linearization [10], that is the frontier between an emerging action potential and a decay to the resting state of the membrane.

Fife and McLoad
2. Analytical stationary wave. In the first place, let us consider the bistable equation in an infinite continuous one-dimensional medium such as In this model of the smooth nerve fibre, V denotes the voltage accross the membrane, D represents the diffusion coefficient, t is the time and x is the longitudinal axis of the fibre. The nonlinear function f (V ) corresponds to a cubic polynomial function so that and represents the sodium current flowing through the membrane. The parameter α acts as a threshold between the passive and the active role of the sodium conductance (0 < α < 1), so that the steady state V = 0 is the resting state and the steady state V = 1 is the excited one.
The system is completed by Neumann conditions, so that lim x→±∞ dV dx = 0. In this context, we look for non uniform stationary waves where V (x = 0) = 0 and lim x→±∞ V (x) = 0, so that the nerve fibre is at its resting state when |x| → ∞. Due to symmetry, we can reduce this problem to a semi-infinite medium by set- Multiplying both sides by dV dx and recognizing that eq. (3) can be rewritten such as Note that this equation could also be obtained by introducting a potential function, as in [11]. Taking into account the Neumann conditions requirements, eq. (5) yields, after integration, where has three real roots and can be expressed by V = 0 and V = 1 are the minima of F (V ), corresponding to stable states, and V = α is a local maximum corresponding to an unstable state. To integrate eq. (6), F (V ) must be positive, therefore excursion of V must lie in the range [0, The left-hand side of eq. (8) corresponds to the integral form of an elliptic function [12] yielding Equation (9) corresponds to the stationary wave whose width depends on the diffusion parameter D and on the threshold parameter α. The amplitude of V s is maximum when x = 0 so that V s (0) = V 1 and is only a function of α. Because of symmetry, this result can be extended to an infinite medium, as illustrated in Fig. 1.
Note that the existence of this stationary wave (in accordance with the results obtained in [11]) is assured only when α ∈]0, 1/2[, which corresponds to the realistic biological case [3]. 3. Projected dynamics diagram. The shape of the stationary wave expressed by eq. (9) is rewritten into a more general form The parameters (a, γ, k) of eq.(10) depend on the value of the nonlinear threshold α, while k also depends on the diffusion coefficient D. Therefore a variation of these parameters from their specified values destroys the stationary behaviour of the system. In order to understand the influence of such a change, we use a variational approach leading to capture the different evolutions of an initially distributed voltage defined by eq. (10). In the following, k and a are defined as variable parameters, while γ remains fixed by the threshold α (this last condition avoids untractable analytical developments). Therefore a is linked to the amplitude of the wave while k is linked to its width. The diffusive coefficient D is only present in the definition of k. In order to simplify the following analytical study, we consider here the case where D = 1. We project the dynamics of the stationary wave in the (a, k) space to determine the projected dynamics diagram. Considering the wave V (x, t) with V (x, t) → 0 as |x| → ∞, let its energy E be defined as where F (V ) is the antiderivative function defined in eq. (6). The energy E can also be expressed by considering a Lagrangian function L, so that with V (x) null on the boundaries and Ω as the domain on which E is defined.
Respecting the Euler-Lagrange conditions [13], if E has an extremum, its variational derivative must be null, therefore Recognizing that L(V (x), Consequently, The gradient flow of eq. (10) is projected onto a two-dimensional space formed by the parameter a and the width k of a pulse. The resulting projected dynamical system is a pair of ODEs describing the temporal evolution of the parameter a and the width of the pulse, and is investigated using the methodology detailed by Neu et al. [9]. Let us recall, the main steps leading to the projected gradient flow: The wave (eq. (10)) is considered as a parametric representation, so that where β(t) = (a(t), k(t)) is a time-dependent vector of parameters. Considering the time derivative of V (eq. (16)), the gradient flow (eq. (11)) implies where δE δV is a linear combination of ∂ j V , j = 1, ..., N . To extract the ODEs for β j (t), the eq. (17) is multiplied by ∂ i V dx and then integrated  The energy is a function of V , i.e function of a vector β of parameters, so that the functional derivative is Therefore, equation (18) Using linear algebraic notation, (eq. (20)) is expressed such as where the gradient operator is taken with respect to β and M is a N × N symmetric matrix. The recalled steps are applied on eq. (10), leading to with m ij determined in the appendix in function of γ. Applying eqs. (21-22) and determining the gradient operator as proposed in the appendix, the projected dynamical system can be expressed such as which describes the evolution of an initial condition defined by eq. (10).
Calculating the nullclines of the system (23), four equilibrium points R(0, 0), P 1 (a P 1 , 0), P 2 (a P 2 , 0) and N (a N , k N ) are found so that R(0, 0) is an attractor point and it characterizes the return of a system to a resting state; P 1 (a P 1 , 0) corresponds to a source; P 2 (a P2 , 0) is a sink and characterizes the excited state and N (a N , k N ) is a saddle, as illustrated in Figure 2.
Depending on values of a and k, we can deduce four domains in the projected dynamics diagram (see Fig. 2). Domains (1) and (2) characterize the stimulation, i.e. initial conditions defined by (a, k) leading to excite a nerve fibre. Domains (3) and (4) are so that initial conditions defined by (a, k) lead to a return to the resting state. Therefore, the stable manifold of the saddle node (a N , k N ) corresponds to a frontier between the emergence of a travelling wave and a decay to the resting state. Determination of this separatrix leads to obtain the relationship between a and k, thus between the amplitude and the width of the initial condition giving the frontier. An approximate solution can be reached by using the Poincaré's linearization theorem. 4. Determination of the stable manifold. Let (x, y) be a basis obtained from the linear basis so that (0, 0) corresponds to N. The x-axis corresponds to the unstable linearized eigenvector of the saddle node N, while the y-axis corresponds to its stable linearized eigenvector. With λ 1 (resp. λ 2 ), the stable (resp. unstable) eigenvalue of N, we can rewrite the where x ∈ R n and y ∈ R n . g i (x, y), i = 1, 2, contains the nonlinear parts of the equation which vanish, as their first derivatives at the origin. The stable manifold theorem is used to realize the Poincaré's linearization [10] of eqs. (24) and (25).
The stable manifold, denoted S(x), is such as y = S(x) and ∂S ∂x (0) = 0, leading to Replacing all terms in y by S(x) in eq. (25) yields Furthermore, The stable manifold of order seven has been determined and compared to the frontier obtained by numerical simulations from eq. (10), as illustrated in Fig. 2. There is a good correlation between the numerical result and the stable manifold of order seven, compared with the first order linearization of the separatrix. Nevertheless, even if the accuracy is increased for small k, it diverges when k exceeds a certain value, which is due to the fact that the eigenvalues of the saddle are resonant, at least numerically. This divergence becomes more and more pronounced as the order n of the approximation is increased, reaching the saddle asymptotically. However, this method can be applied to determine the (a, k) relationship in the range 0 ≤ k ≤ k N to whatever desired accuracy. Otherwise, a compromise has to be made between accuracy and the possible variation of k.
In this study, the criterion for the initiation of the travelling waves has been investigated, based on the existence of stationary waves. We have also determined analytically the existence of a frontier between the initiation of a travelling wave and its return to the resting state in a projected dynamics diagram. We suggest that our results could be applied to determine the minimum punctual electrode stimulation or the minimum size of the electrode needed to generate an action potential. It would also be interesting to apply this method in the bidimensional excitable media case modelling cardiac tissue.