DYNAMIC TRANSITIONS OF THE FITZHUGH-NAGUMO EQUATIONS ON A FINITE DOMAIN

. The main objective of this article is to study the dynamic transitions of the FitzHugh-Nagumo equations on a ﬁnite domain with the Neumann boundary conditions and with uniformly injected current. We show that when certain parameter conditions are satisﬁed, the system undergoes a continuous dynamic transition to a limit cycle. A mixed type transition is also found when other conditions are imposed on the parameters. The main method used here is Ma & Wang’s dynamic transition theory, which can be used generally on diﬀerent set-ups for the FitzHugh-Nagumo equations.


1.
Introduction. Propagated signaling in axons is an important phenomenon in neural science, various models have been put forward to explain this behavior, the Hodgkin-Huxley (HH) equations [8] is a satisfactory and experimentally supported model. It is based on the assumptions that axon membranes can be regarded as capacitors, while different types of ion channels can be viewed as linear resistors, of which the conductance are measured from voltage-clamped data. The HH equations then can be constructed using an equivalent electric circuit.
Due to the difficulty of analytical mathematical study, several simplifications of the HH model have been made, including the FitzHugh-Nagumo (FN) equations, which have been extensively studied, for example in [4] [6][10] [11].
The FN equations read as follows: where 0 < a < 1, > 0, γ > 0. v corresponds to the voltage, y is associated with the recovering behavior of the opening of voltage-gated (v.g.) K + channels and the closing of N a + channels after voltage increases, a roughly represents the threshold of excitation for v, is a small number to fit the delayed response of y, and γ is 3936 YIQIU MAO another adjusting parameter. The FN model captures some behaviors of neurons, including action potentials, impulse trains, anodal break excitations and so on [4]. In a "space clamped" version of the FN model, the ∂ 2 v ∂x 2 term is removed and the dynamics of this equation can be explored using theories for finite dimensional dynamic systems. In particular, phase planes can be used to visualize the behavior of solutions. Finite dimensional methods can also be used for traveling wave solutions. Namely, if we set v and y to be functions of x + θt only, θ being a constant for traveling speed, then the equations can be converted to finite dimensional ODE systems. However, in a more realistic case when the axon length is regarded as finite and if proper boundary conditions are added, (1) becomes an infinite dimensional system.
In this paper, using the dynamical transition theory developed by Ma & Wang [14], we provide a complete characterization of the nonlinear stability and transition behavior of the system. One of the key ingredients is to reduce the underlying infinite dimensional system to a finite dimensional system generated by the unstable modes. The reduced system offers the complete information on the types and structures of the dynamic transitions.
Specifically, this paper will deal with the case where x ∈ [0, L], corresponds to an axon with length L, and with the Neumann boundary conditions, which amounts to saying that the axon is sealed at both ends, and with an uniform injection of current along the length of the axon. The result shows under appropriate restrictions on parameters, (among them, the smallness of ) when current I increases to a critical value, the system undergoes a Hopf bifurcation to a limit cycle, furthermore, the transition is continuous, which means that the size of the cycle gradually increases with current I. Under another restriction on parameters, however, a different type of transition is discussed, this transition is considered to be possibly irrevelant to real behavior of axons, but more out of the complications of the FN model.
The paper is arranged as follows: Section 2 introduces the FN model as a simplification of HH model and some possible equation setups, Section 3 provides necessary mathematical settings, Section 4 and 5 calculate the dynamic transitions, Section 6 explores some biological conclusions from those transitions.
2. The Fitzhugh-Nagumo model. The FN model can be regarded reasonably as the simplification of the original Hodgkin-Huxley equations: Here v again represents the voltage inside the membrane along the axon while the voltage outside is 0, I m represents the total ionic current goes through the axon membrane per unit distance, C is the capacitance of the membrane per unit distance, d is the diameter of the axon, R L is the specific resistivity of the fluid inside the membrane, barred quantities are constants. g N a and g K represent the conductance of v.g. (voltage-gated) sodium and potassium channels respectively. g l is the conductance of the resting channels.v N a ,v K andv l represent the equilibrium voltage if only v.g. N a + , K + and resting channels are open respectively. m, n represent the opening behavior of N a + and K + channels, while h measures the closing of N a + channels. τ h , τ m , τ n , m ∞ , h ∞ and n ∞ are determined experimentally as a function of voltage. Generally speaking, the steady state of h decreases with increasing voltage while n and m increases with voltage. And the time constants for m to reach steady state are much faster than that of h and n's. This model corresponds well with experiments [8].
Paper [1] offers an explanation of the connection between the HH equation and the FN equation. Roughly speaking, since m behaves much more rapidly than n and h, we can treat the opening behavior of N a + channels as simply a function of v. Allowing v = 0, I m = 0 to represent the rest state, we can modify I m as follows: where C 1 and C 2 are constants. The first term at right hand side corresponds to the opening behavior of N a + channels, second term z represents the collective behavior of the closing of N a + and the opening of the v.g. K + channels. The functions in I m are chosen to roughly fit the voltage-clamped data of axon membranes. After suitable choosing for the constants and non-dimensionalization, we acquire (1). There are several set-ups that can be used to explore the dynamic behavior. Here are a few: Space clamped: in this case a metal wire is added inside the cytoplasm along the length of the axon so the voltage remains same throughout the axon, ∂ 2 v ∂x 2 is removed and the equations become ODEs. Adding I at right hand side of the first equation would amount to adding a injected current. The dynamic behavior of this system of ODE is well studied in [4]. Several notable results are: 1: The absence of an all-or-none response to injected current, as opposed to the possible existence of a firing threshold observed in real neurons [9]. However a "quasi-threshold" does exist, see [3], [4]. 2: Anodal break excitation. In which the injected hyperpolarizing current is abruptly stopped and a transient spike is resulted. 3: A stable limit cycle occurs when the injected depolarizing current is big enough. Traveling wave solution: If v is a function of x + θt then the solution can be regarded as a steady propagating wave. In this case (1) can be regard equivalently as a 3-dimensional ODE. A stable homoclinic orbit has been found [6], [10], which corresponds to a single propagating pulse. Dynamics of this system have also been extensively studied for example in [11]. Finite domain with boundary conditions: In this paper we mainly focus on bounded domain, x ∈ [0, L], L > 0, then there are two possible choices for boundary conditions: This corresponds to the case when current I is injected uniformly into an axon with length L, both ends of which are sealed (no current going through). Equations of this form allow for simpler mathematical manipulations, and we will consider this setup in current paper.
This represents a more realistic case, exciting current I is injected at one end of an axon, while the other end being sealed (forming synapses). However the steady solution of this equation is hard to calculate without computation tools.
3. Mathematical setting and basic states. We will investigate equations (3) in this paper. For simplicity only constant steady states are considered. Those steady state solutions are given by: The solution near v = 0 can be written as: where v I is a function of I. It is chosen so that v 0 = 0 and v I is continuous with I in a neighborhood of I = 0. It's easy to see v I either keeps increasing with I or initially increases and then reaches a maximum at some value of I, whether v I reaches a maximum depends on the value of γ and a. Elementary computations show: • when a 2 −a+1 3 < 1 γ , v I continuously increases with I and can reach infinity.
where v m is the smaller root of the quadratic equation: After v m , if we continue increasing I, the only steady state solution will occur at a distance from v m . Before v m , there are multiple roots for (5) including the v I branch. In total, the set of solution v of equation (5) still contains (0, ∞).
Let new v equal to v −v I , new y equal to y −v I /γ, then the equations (3) become: We introduce a new parameter: Using Then, equation (7) can be characterized by: where u 0 = v0 y0 . The operators L λ : H 1 → H and G λ : Here the Hilbert spaces H 1 , H and H 1/2 are defined by: It can be verified that L is a sectorial operator (using for example Theorem 1.3.2 of [7]) and enjoys spectral properties as a linear completely continuous field (see [13]), G is a C ∞ bounded operator, and . The global existence and uniqueness of solutions to equation (8) are straightforward because of the negative leading term in G.

Principle of exchange of stabilities.
To study the dynamic transitions of this problem, first consider the linearized eigenvalue problem of (8): or: To solve this equation, we use separation of variables: where C 1 and C 2 are constants. It's straightforward to calculate the solutions as: To avoid the case when β is a multiple root, which will cause unnecessary trouble in the discussion follows, we require the ∆ of (13) to be nonzero: Assume (14) is satisfied, then the two roots β k1 and β k2 of equation (13) are different. If both of them are real, the eigenvectors of (8) can be written as: (without loss of regularity, let β k1 > β k2 ) and L λ Ψ ki = β ki Ψ ki , i = 1 or 2. For other k's, the eigenvectors can be written (assume Imβ k1 > 0): and L λ Ψ k2 = Imβ k1 Ψ k1 + Reβ k1 Ψ k2 . It's evident this is a complete set of eigenvectors for H.
From (13), the following Theorem for principle of exchange of stabilities can be proved: Theorem 4.1. For the linearized equation (12), we have the following: (1) With conditions (14) and the following holds true: the following holds true: where λ 0 = γ.
Proof. The first case: When λ = λ 0 , if k = 0, equation (13) will have two roots 0 and 1/γ − γ, then the second inequality shows β 02 < 0. For other k's, elementary quadratic equation analysis shows their eigenvalue β ki 's must have negative real parts for λ 0 . For β 01 , similar analysis shows in a neighborhood of λ 0 , the greater root of (13) with k = 0 increases with λ. The second case: When λ = λ 0 , if k = 0, the quadratic equation (13) will have two pure imaginary roots: ± γ(1/γ − γ)i. They are non-zero numbers according the second inequality. For other k's, since the symmetric axis lies to the left of the y-axis, corresponding eigenvalues who are non-real numbers will only have negative real parts. If there is a real non-negative eigenvalue, then calculation shows λ = 1/γ + k 2 π 2 /L 2 , hence at λ = λ 0 there is no k for this equality to hold because λ 0 < 1/γ. As with β 01 , increasing λ will move symmetric axis to the right hence the real part will be positive when λ > λ 0 and negative when λ < λ 0 .
If all the βs have negative real parts, the steady state 0 of (8) is stable though attractor basin can vary. Our purpose is to discuss what will occur if the change of λ causes eigenvalues to cross the imaginary axis. For the next section, we investigate two possibilities: the first crossing occurs at real and at complex eigenvalues.

Remark 1.
There is a slight complications of parameters, namely, for the nonlinear term in the full equation (7) to be defined, a 2 −a+1 3 ≥ λ must be satisfied. Thus for the PES to take place in the full equation, one of the following two parameter relations must be satisfied for each of the real and complex transition case: When neither of these inequalities are satisfied, there will be no eigenvalue crossing of the imaginary axis, and the constant steady state solutions of (7) remain locally stable.

Dynamic transitions.
Thanks to the principle of exchange of stabilities, we can make use of Ma & Wang's dynamic transition theory to explore the dynamics of equation (8) near the origin when λ crosses λ 0 . Their theory classifies all transitions generated by the unstable modes into three categories: Type-I (Continuous), Type-II (Jump) and Type-III (Random). For Type-I transition, the transition states are characterized by a local attractor, which attracts a neighborhood of the basic solution. This means that as the control parameter crosses a critical value, the transition states stay in the close neighborhood of the basic state. The onset of the Rayleigh-Benard convection [12] is an example for this type. For Type-II transition, the transition states are represented by some local attractors away from the basic solution and the system undergoes a more drastic change. For Type-III transition, the transition states are represented by two local attractors, with one as in Type-I transition and the other as in Type-II transition. In this case the fluctuations of the basic state can be divided into two regions such that fluctuations in one of the regions lead to continuous transitions and those in the other region lead to jump transitions. We refer to [14] for details. For the equation (8), we have the following transition theorem: (14) is satisfied for a neighborhood of λ = λ 0 we have the following for equation (8): 1 Real transition case. Assume condition (17) and (19) is satisfied, then: (1) Equation (8) has a type III (random) transition from 0 at λ 0 = 1/γ. More precisely, there exists a neighborhood U ⊂ H of u = 0 such that U is separated into two disjoint open sets U λ 1 and U λ 2 by the stable manifold Γ λ of u = 0 such that the following hold: a. U = U λ 1 + U λ 2 + Γ λ ; b. the transition in U λ 1 is a jump transition; c. the transition in U λ 2 is continuous. (2) Equation (8) bifurcates in U λ 2 to a unique singular point v λ on λ > λ 0 (which correspond to I < I 0 ) that is an attractor such that for every 2 Complex transition case. Assume condition (18) and (20) is satisfied, then the equation (8) undergoes a type I (continuous) transition to a periodic orbit when the parameter λ across γ, which correspond to the value of I 0 with v I being the first root of λ = γ (hence reachable by initially increasing I). The bifurcated periodic solution can be expressed as: x 1 (t) = Reβ 01 |b| 1/2 cos((Imβ 01 )t) + o(|Reβ 01 |), Proof. First we will discuss the real transition case: Denote the linear space spanned by Ψ 01 as E λ 1 , E λ 2 is spanned by rest of the eigenvalues so that H 1 = E λ 1 ⊕ E λ 2 . let P i : H → E λ i be the canonical projection, let Φ be the center manifold function, u = xΨ 01 + y, then projecting the equation (8) to E λ 1 will yield: dx dt = β 01 (λ)x + P 1 G(x + y, I), dy dt = L λ y + P 2 G(x + y, I).
By the existence of center manifold, near the origin the first equation can be written as: Here because elementary linear algebra shows: By Theorem A.1.1 of [14], Hence Then it follows Putting the center manifold function in (21), we obtain At critical crossing, β 01 = 0, hence the equation becomes: Since β 02 < 0, we will have Then, according to Theorem 2.3.2 of [14], we derive the conclusion for the type III transition stated in the theorem.
By the existence of center manifold, near the origin first equation can be written as: By Theorem A.1.1.
(2) of [14], here x = x1 x2 . After putting the center manifold functions into (24), we will obtain At critical crossing, Reβ 01 = 0, hence the equations become: Thanks to (18) and (20), this is always a negative number. Then, according to Theorem 2.3.7 of [14], the assertions for complex transitions are proved.

Remark 2.
As is stated in Theorem 5.1, if γ, , a satisfy (18) and (20), the complex transition of (8) with parameter λ also happens in (7) with parameter I. However, if γ, , a satisfy (17)and (19), then the random type transition will not occur by increasing of parameter I from I = 0 in equation (7). It's because, as is discussed in (6), v I will reach a maximum when λ = 1/γ, hence if we increase I, λ will reach 1/γ but will not cross it. So effectively the transition will not occur by increasing of I, although (23) still captures the local behavior around 0 when λ reaches 1/γ.

6.
Conclusions. The complex transition shows that when γ is fixed, whenever is small enough, the system will undergo a continuous type dynamic transition from the basic state to a spatio-temporal oscillating pattern (Hopf bifurcation). This corresponds to the behavior of repetitive firing of neurons when there is a large enough persistent stimulating current. Though in experiment this repetitive firing is damped after several spikes [5], it was thought this damping was due to poor condition of the axon. Hence the continuous dynamic transition to a limit cycle provides a good description to the real axon behavior. What our model doesn't capture is the space variation of the bifurcated solutions, since all the non-constant functions will be eventually damped out as proved above. This might be due to the uniform injection setup and symmetric boundary condition. Nonetheless the method used in this paper provides an effective tool to reduce a PDE to finite dimensions and also gives a characterization of the dynamic behavior. With some computation tools, dynamics of (4) might provide more spatially variant details of the transition. The requirement for to be small shows that if the opening of voltage gated K + channels and the closing of N a + channels are not slow enough, the transition to periodic orbit won't occur.
One possible discrepancy with real axon is that, because of the continuous transition, limit cycle gets continuously larger as λ across γ. While in the real axon experiments in [5], data shows that the amplitudes of spikes are almost constant. This is possibly because of the same reason underlying the apparent paradox between the all-or-none behavior of real axons and the lack of a real threshold in HH models, namely, a small change of parameters can cause large changes in solutions even though the dependence on parameters is continuous.
The frequency of the resulting impulse trains, according to Theorem 5.1, is inversely proportional to Imβ 01 at least when λ is in a small neighborhood of λ 0 . It can be calculated that Imβ 01 decreases as λ increases from λ 0 , this corresponds well with the observed firing behavior of neurons. Namely, the larger the current, the more frequently neurons fire.
The real transition observed in (8) might be degenerated in the FN equations (3). If we denote the value of I when v I reaches v m as I m , then although no transition