Traveling wave solutions of Fitzhugh model with cross-diffusion

The Fitzhugh-Nagumo equations have been used as a caricature of the Hodgkin-Huxley equations of neuron firing to better understand the essential dynamics of the interaction of the membrane potential and the restoring force and to capture, qualitatively, the general properties of an excitable membrane. Even though its simplicity allows very valuable insight to be gained, the accuracy of reproducing real experimental results is limited. In this paper, we utilize a modified version of the Fitzhugh-Nagumo equations to model the spatial propagation of neuron firing; we assume that this propagation is (at least, partially) caused by the cross-diffusion connection between the potential and recovery variables. We show that the cross-diffusion version of the model, besides giving rise to the typical fast traveling wave solution exhibited in the original diffusion Fitzhugh-Nagumo equations, also gives rise to a slow traveling wave solution. We analyze all possible traveling wave solutions of the Fitzhugh-Nagumo equations with this cross-diffusion term and show that there exists a threshold of the cross-diffusion coefficient (the maximum value for a given speed of propagation), which bounds the area where normal impulse propagation is possible.


Introduction
Hodgkin, Huxley, and Katz in the 1940's explored mathematically and experimentally the nature of nerve impulses. Their work revealed that the electrical pulses across the membrane arise from the uneven distribution between the intracellular fluid and the extracellular fluid of potassium (K + ), sodium (Na + ) and protein anions (Sherwood 2001). When a neuron is not sending a signal, it is said to be at rest (and at approximately -70 mV) and the inside of a neuron is more negative relative to the outside. The change in the Na + and K + permeability allows for the movement of ions in and out of the cell by means of opening and closing of ion channels. The influx of Na + and efflux of K + results in electrical potential difference. Triggered by a stimulus the Na + channels open, the influx of Na + ions increases, the membrane depolarizes, and the potential voltage reaches a threshold level typically between -50 and -55 mV. At this time an explosive depolarization takes place, which rapidly moves the potential to a maximum of +30 mV (Sherwood (2001)). At a similar rate, the Na + and K + channels close and open, respectively, initiating membrane repolarization caused by an efflux of K + ions. The potential drops back to resting potential.
The intensity and magnitude of the forces behind repolarization are such that the membrane goes through a transient hyperpolarization of -10 mV below resting potential. This entire process of rapid change in potential from threshold to peak reversal and then back to the resting potential level is called action potential, impulse, or spike (see schematic diagram in Fig. 1).
This process by which a neuron fires was mathematically investigated by Hodgkin and Huxley, in 1952, with a four-variable model. In an effort to construct a simpler mathematical model of an excitable membrane FitzHugh (1961) proposed a two variable model. This model made it possible to illustrate the various physiological states involved in an action potential (such as resting, active, refractory, enhanced, and depressed) in the phase plane. One such portrait, which shows the spike-like behavior of the model, is given in Fig.2. The FitzHugh system, although a caricature of the Hodgkin-Huxley four equations, captures much of the same dynamical behavior. A more realistic model is one that depends on both space and time since electric currents cross the membrane of the cell and move along its axon lengthwise inside and outside. This mechanism makes it possible for electrical signals to be transmitted over long distance and thus propagate throughout the membrane without ever weakening or decreasing their initial strength. With this in mind a mathematical model of the diffusion of current potential was first proposed and studied by FitzHugh (1961FitzHugh ( , 1969 and Nagumo et al. (1962). Recent models have been proposed where the spatial solutions are conditioned by the effects of cross-diffusion "control" or "interactions" between components of the system (see, for example, Kuznetsov et.al. (1996), Othmer and Stevens (1997), Berezovskaya and Karev (2000), etc.).
Motivated by these works we modified the FitzHugh model to include a crossdiffusion connection between the potential and recovery variables. We suppose that, due to the semiconductor nature of the nerve membrane, the cross-diffusion regulation plays an important (perhaps, crucial) role in the spatial spreading of potential. This version of the model will provide an avenue for investigating successful propagation of an excitable neuron but also propagation failures, which are extremely important for many applications.
For example, the influence of certain drugs or external chemicals affect the rate at which sodium channels close and the rate at which potassium channels open, thus altering the normal dynamics of a firing potential membrane. A study conducted by Gubitosi-Klug and Gross (1996) showed that certain metabolites in ethanol accelerate the release of potassium ions from the brain cells. The increased potassium efflux in turn makes it difficult for cells to absorb enough calcium and thus inhibits the release of neurotransmitters (Highfield, 1996). The changes in the release of potassium ions result in the changes in the recovery phase of the excitable membrane. We include this effect of a generic drug by incorporating a cross-diffusion term in the original FitzHugh model.
In this work we explore the changes of the characteristics of the spatial propagation of nerve impulses brought by changes in the velocity of propagation and intensity of the cross-diffusion regulation. In particular, we are interested in the conditions of "normal" neuron firing propagation and investigate its possible violations.
The paper is organized in the following manner. Section 2 contains a brief description of local neuron dynamics within the framework of the FitzHugh model and the bifurcation portrait of the model. A cross-diffusion modification of the FitzHugh model aimed at providing an explanation of spatial modes like "traveling waves" is contained in Section 3. We show that fast and slow traveling waves can appear with respect to parameter values and follow their dependence by bifurcation analysis of corresponding wave systems; we show also that the "traveling spike" appears only before some threshold of the cross-diffusion coefficient. Section 4 contains the discussion of obtained results.
Proofs of the statements of section 3 are given in the Appendix.

The Fitzhugh equations as a local membrane model
The original Fitzhugh model (1961) describing dynamics of the physiological states of a nerve membrane contains membrane potential variable P and recovery variable Q (which plays the role of all other three variables in the Hodgkin-Huxley model, sodium activation/inactivation and potassium activation). The variable P shares the properties of both the membrane potential and excitability and thus describes the dynamics of the rising phase of neuron firing. The variable Q is responsible for accommodation and refractoriness and thus represents the dynamics of the falling phase of the action potential. The equations are given by where I , ρ, a and b are parameters of the system. The stimulating current is defined by the variable I; I<0 corresponds to a cathodal shock and I>0 corresponds to an anodal shock.
Additionally, it was shown (Khibnik et. al 1998) that a bifurcation of co-dimension 4 with symmetry "3-multiple neutral singular point with the degeneration" is realized in the vector field defined by system (3) in a vicinity of the parameter point M (k 1 =1, k 2 =0, ε=1). The main results of these works are summarized in this section below.
Let us consider the phase-parametric portrait (the bifurcation diagram) of the vector field (3), i.e., the partition of a vicinity of the parametric point M into all possible domains with topologically different phase portraits in a neighborhood of the point O (see Arnold (1993), Kuznetsov (1997)). For the sake of brevity, we call a cycle "small" if it contains a unique singular point, and "large" if it contains three singular points. (1) Surfaces , under conditions that 0<ε §1; (2) Lines SS: k 1 = 1, k 2 = 0, 1 Volokitin and Treskov (1994)  Let us emphasize that in the framework of the FitzHugh model the spike-regime (see 3 Cross-diffusion model of a transmembrane potential

Extending/Modifying Fitzhugh model
In this section, we propose an extension of the FitzHugh spatial model to include the implicit (hypothetical) cross-diffusion mechanism of the spatial propagation of the firing process of the neuron. As mentioned previously, a cross-diffusion term involving the recovery variable could be used for modeling the effect of a generic drug that affects the flow of sodium and/or potassium.
We first give a very brief review of the well-known FitzHugh-Nagumo model where t is time, x is a one-dimensional space variable and non-negative constants D, σ are the cross-diffusion and diffusion coefficients, respectively. The "diffusion" version of system (4), which corresponds to the case D=0, σ>0 is known as "model FitzHugh-Nagumo" (Nagumo et.al. (1964), FitzHugh (1969)). Many works were devoted to the study of its dynamics, and in particular, to the investigation of "traveling wave" solutions (see, e.g., Hastings (1976), Evans et al. (1982), etc.). A bifurcation approach applied to the study of traveling impulses and trains (see Kuznetsov (1997)) allowed to reveal "fast" and "slow" waves that can exist with the same values of "local parameters" ε, k 1.
To make clearer the role of the spatial distribution of the recovery variable (along the axon) and the meaning of the cross-diffusion term in the impulse propagation, we consider a cross-diffusion version of system (3): Mathematically, cross-diffusion equations possess special properties, which facilitate their research (Wei-Ming Ni, 1998, Berezovskaya, 1998. In what follows, we explore "traveling wave" solutions of system (5): where ξ = x + Ct and positive C is the velocity of the wave propagation. It can be checked that (p(ξ), q(ξ)) satisfy the following two-dimensional "wave system": For convenience, let us denote α = C 2 / (εC 2 -Dk 1 ) and change the independent variable h = ξ/C, then the wave system (6) becomes where β= sign(α). This system is defined for εC 2 ∫ Dk 1 .
Further, we study phase behaviors of (7) dependending on "local" parameters ε, k 1 , k 2 for arbitrary fixed values of constants D and C such that C 2 ∫Dk 1 /ε. In other words, the 8 domain of the parameters D, C under fixed values of ε, k 1 , k 2 is divided into two domains, in which the system exhibits qualitatively different behavior, by the parabola

Wave system of the cross-diffusion Fitzhugh model
System (7+) with C, D such that C 2 >Dk 1 /ε is called fast wave system, whereas system (7−) with C 2 <Dk 1 /ε is called slow wave system of model (5). For both cases system (7β), evidently, has from one up to three singular points (p*,q*) whose coordinates satisfy equations: Two singular points coincide and form two-multiple point with the parameter values belonging to boundaries S ≤ (see Fig.3 and Theorem 1). Thus, both the fast and slow wave systems have three singularities inside the curvilinear angle that is formed by S ≤ , and a single singular point outside the angle. Note that this unique singular point is a non-saddle (a node or spiral) for the fast wave system whereas it is a saddle for the slow one. Inside the parameter angle the fast wave system has two non-saddles and a saddle whereas the slow one has two saddles and a non-saddle (see Mathematical Appendix, Proposition 3). At (ii) Let 0<C 2 <Dk 1 /ε. There exist a neighborhood of the parameter point M (k 1 =1, k 2 = 0, ε=1) in which the vector field defined by system (7-) in a neighborhood of the phase point (p=0,q=0) has a bifurcation diagram, whose cut to the plane (k 1 , k 2 ) is topologically equivalent to the one presented in Fig.5a, 6 for arbitrary fixed positive ε<1 and in Fig.5a, 6 for arbitrary fixed ε>1.
The boundary surfaces in the parameter space correspond to the following bifurcations: Remark. The proofs of Theorems 1 and 2 (given in the Appendix) contain the analytical description of local bifurcations of the systems. Non-local bifurcations were investigated mostly by computation.

Traveling wave solutions of PDE and their profiles as solutions of wave ODE
The correspondence between traveling wave solutions (6) of model (5) and orbits of its wave system (7) is schematically given in Figures 7 through 9 (see, e.g., Volpert et al Berezovskaya andKarev (1999, 2000)).
Heteroclinic, homoclinic orbits, and limit cycles of a wave system correspond to wave fronts, pulses, and train solutions of the model, respectively. A heteroclinic curve connects two singular points; a homoclinic curve connects a single singular point to itself; and a phase limit cycle encloses an odd number of singular points. 4 see middle part of Fig.7a,b Thus, the description of all possible wave solutions of PDE system (5) is reduced to the analysis of phase curves and bifurcations in the wave ODE system (7) in which C is an "additional" parameter. We make this correspondence more formal with the following statement. iv) a wave train of the model corresponds to a limit cycle in the (u,v) phase plane of the wave system, see Fig.9.

Fast and slow traveling waves of the cross-diffusion system
We say that model (5) and three different trains in the areas 11 and 13a-b (see Fig.4); the impulses for the parameter points (k 1 , k 2 , ε) belonging to the boundaries P + , Pand R + , R -, see Fig.8a,b,c. the trains in the domains 4`a, 5` (see Fig.6); the impulses for the points (k 1 ,k 2 ,ε) belonging to the boundaries P + , P -(see Fig 8) -

Stationary distribution of the model
The stationary spatial distribution of model (5) is defined by the equations P t = Q t = 0. Then system (5) takes the form: and can be easily reduced to the second order equation with respect to P and x: Note that the wave system of model (5): coincides with system (8) if C=0. Traveling wave solution with C=0 is called a standing wave. One could write (9) as the Hamiltonian system whose Hamiltonian is of the form: H(P,W)=(2Dk 1 W 2 +2 P 2 (1-k 1 )+4Pk 2 -P 4 )/ (2Dk 1 ) + const.
System (10) has from one (a center) up to three (two saddles and a center) equilibria (P*, 0), where P* satisfies the equation The parameter portrait of (10) presented in Fig.10a is the "light version" of the portrait in  Fig. 10 a,b).

Proposition 2.
Depending on parameters k 1, k 2 (see parameter portrait in Fig.10, left) three types of standing waves are realized in model (5) for suitable initial values (p,q) : 12 i) trains corresponding to limit cycles (Fig.10a,b,c); ii) small impulses corresponding to small homoclinics (Fig.10a,b) and iii) fronts corresponding to heteroclinics (Fig. 10c) 4. Discussion and conclusion

Comparing of bifurcation diagrams
We have shown that the parameter portrait of the local FH-model (3)  We emphasize that the parameter portrait in Fig.3 of the local FH-model also corresponds to the wave system of the cross-diffusion model with a small cross-diffusion coefficient D (under fixed propagation speed C) while the parametric portrait in Fig.5 corresponds to the same model with "large" cross-diffusion coefficient D. Hence, these portraits describe the model behavior before and after the threshold D= εC 2 /k 1 accordingly.

Scenarios of appearance and transformations of the traveling waves
The problem of our interest is the appearance and transformations of the traveling wave solutions depending on the model parameters D and C that characterize the axon abilities for the firing propagation. One could assume that these characteristics may change as a result of influence of certain drugs or external chemicals. We now trace the transformation of the traveling wave solutions by varying these parameters under the supposition that parameters k 1 , k 2 , ε have arbitrary fixed values close to k 1 =1, k 2 =0, ε=1.
Let the (positive) value of the speed propagation C be fixed and suppose the positive crossdiffusion coefficient increases. For D=0 the wave system of the model coincides with the local Fitzhugh model. This model demonstrates a spike (shown in Fig.1) if (k 1 *, k 2 *, ε * ) 13 belongs to the boundary R + . The wave system describes "pseudo-waves" and, in reality, there is no firing propagation.
On the contrary, if D> ε/k 1 C 2 then a traveling spike does not exist. The only possible traveling waves are "small" trains, impulses or fronts with the amplitudes less than is degenerate and describes the curve q = -p 3 + p(1-k 1 ) + k 2 that smoothly joins the points p 1 , p 2 , p 3 . Further increase of C leads to a "transformation" of the slow wave system into the fast one, and thus to the appearance of waves similar to the spike spreading along an axon. A behavior of the model under critical values D= k 1 C 2 /ε, evidently, cannot be studied in the framework of the two-dimensional model (5).

Conclusion
We utilized a modified version of the Fitzhugh-Nagumo equations to model the spatial propagation of neuron firing; we assumed that this propagation is essentially caused by the cross-diffusion connection between the potential and recovery variables. This modification, which includes the implicit (although hypothetical) cross-diffusion mechanism could help to explore the effect of a generic drug in the neuron firing process. The incorporation of the cross-diffusion mechanism is motivated by the experimental observation that certain the local model (3) is transformed to the generalized Lienard form: Note, that model (5) after transformation (12) reads as the cross-diffusion modification of (13) with the coefficient Dk 1 : A traveling wave solution of system (14-15) is defined as a pair of bounded functions where C>0 is a velocity of propagation. One can verify that after introducing an independent variable h = ξ/C the functions {u(h), z(h)} satisfy the wave system Here α = C 2 /(εC 2 -Dk 1 ) if C 2 ∫ Dk 1 /ε, β= sign(α) and F(u,z;α) = αf(u) +αzg 1 (u) + αz 2 G(u,z) with f(u), g 1 (u) and G(u,z) given by (14). Note, that α = 1/ε if D=0 and ε∫0. It means that the vector field defined by system (17+) coincides with the vector field defined by system (13).
Let's now replace the capital letters in (12) by small letters, reduce p and q via p = (z+u)/k 1 , q = u-k 2 (k 1 ∫ 0), and substitute into system (7β). Dividing equation for u ξ by C∫ 0 we get system (17β).
For any α, 0<α < ¶ vector field (A1, A2) has at least one singular point (u 0 , 0) where u 0 is a root of the cubic polynomial: For arbitrary δ 1 ∫ 0 f(u) may have two additional roots. An appearance of these roots corresponds to the appearance of two additional singular points of the vector field.
Generally, these additional points are a saddle and a node. They are seen to come into existence when -au 3 + uδ 1 + δ 2 =0, -3au 2 + δ 1 =0 and one can easily find the boundaries of the curvilinear angle in (δ 1 ,δ 2 )-space that correspond to existence of two-multiple roots u 0 of polynomial f(u) .

Beginning of proofs of Theorems 1,2
The following statement is given in [Khibnik et al., 1998].
Proof of Corollary 2. For any δ 3 <0 the divergence of vector field (A1, A2), divJ=a(g 1 (u)+2zG(u,z)+z 2 G z (u,z))= a(( δ 3 -3au 2 ) -3az(2u+z)), keeps its sign in some neighborhood of zero phase point containing inside all singular points of the field: divJ<0 for a >0 and divJ>0 for a <0. Thus, the vector field has no closed orbits in this neighborhood; any non-saddle singular point (if exists) is a sink for a >0 and a source for a <0. It proves the Corollary in the case a >0 and shows that for a<0 the field (A1, A2) has no limit cycles, homoclinics, and heteroclinics pairs that join the same pair of saddles.
Existence of boundaries L + for δ 1 >0 and Lfor δ 1 <0 for a <0 in the parameter portrait was checked numerically.

5.5.Local bifurcations of the vector-fields (A1),(A2) depending on α
Let us define in the space of parameters {δ 1 , δ 2 , δ 3 , α} of the vector field J the following sets: Coming back to the initial parameters k 1 , k 2 ,ε (see Proposition 2) we get the formulas given in Theorems 1,2.
The bifurcation "two-multiple cycles" is realized on the surface C, which touches the surfaces H ± by lines CH ± ; The bifurcation "a small loop composed by one of separatrix pairs of the saddle point" is realized on the surfaces P + , Pand "a big loop composed by one of separatrix pairs of the saddle point" is realized on the surfaces R + , R -. Surfaces S + , H + and P + have common lines SH + (see Fig. 3a) as well as S -, Hand Phave common lines SH - (Bogdanov, 1973).
Surfaces R + , Rhave common lines of touching with surfaces S + , S -. They also intersect bifurcation surfaces R + , R -(see Fig. 3) such that four separatrixes compose "8"at the phase plane; this bifurcations was studied by D. Turaev (1985).
The bifurcation "a small loop composed by one of separatrix pairs of the saddle point" is realized on the surfaces P + , P -. Surfaces S + , H + and P + have common lines SH + as well as S -, Hand Phave common lines SH -.
The bifurcation "upper, lower (respectively) heteroclinics of saddle singular points" is realized on the surfaces L +, Las for δ 3 >0 so for δ 3 <0. For δ 3 >0 L +, Land P +, Phave common line of intersection such that two heteroclinics simultaneously join saddle points.