Dynamics of Nerve Pulse Propagation in a Weakly Dissipative Myelinated Axon

We analytically derived the complex Ginzburg-Landau equation from the Liénard form of the discrete FitzHugh Nagumo model by employing the multiple scale expansions in the semidiscrete approximation. The complex Ginzburg-Landau equation now governs the dynamics of a pulse propagation along a myelinated nerve fiber where the wave dispersion relation is used to explain the famous phenomena of propagation failure and saltatory conduction. Stability analysis of the pulse soliton solution that mimics the action potential fulfills the Benjamin-Feir criteria for plane wave solutions. Finally, results from our numerical simulations show that as the dissipation along the myelinated axon increases, the nerve impulse broadens and finally degenerates to front solutions.


Introduction
Wave propagation and pattern formation in a variety of excitable media can be effectively described by reaction-diffusion equations.The FitzHugh Nagumo (FHN) model [1] [2] is one of the simple examples of a two-dimensional excitable dynamics that governs such systems.This model have been used as a simplified version of the Hodgkin-Huxley model [3] of neuron firing.The FHN model which is characterized by a recovery mechanism furnishes us with a better understanding of the essential dynamics of the interaction of membrane potential and qualitatively captures 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.
Myelinated axons modelled by discrete FHN equations incorporate a richer dynamics than its continuous counterparts.Also, the mathematical study of spatially discrete models is challenging because of special and poorly understood phenomena occurring in them that are absent if the continuum limit of these models is taken.For example, in the discrete FHN model, there exists a coupling threshold for its propagation while the continuous system sustains propagation of all coupling strengths [4]- [7].Enormous efforts to understand propagation failure in FHN system were made by Booth et al. [6], in which they considered slow recovery and very special limiting values of the parameters characterizing the bistable source and the spatial diffusivity in the FHN system.By imposing some special boundary conditions, they studied the evolution of localized front from one cell to another until it failed to propagate.Also, by using the Morris-Lecar model for myelinated nerve fiber, Hastings et al. [8] demonstrated that travelling waves could be observed even though they had difficulties to extend their results to the FHN system.Our main focus is to capture the various profiles of nerve impulses along a myelinated axon when the dissipative effects is reduced to its minimum value.This will be experimentally very difficult to be realized, but we hope to theoretically investigate this effect by using perturbation techniques.
The effects of dissipation is prominent in most physical and biological systems [9]- [11], such as the myelinated FHN fiber which is a discrete system made of periodic structures called Ranvier nodes [12].Dissipative systems are spatiotemporally organized patterns that are far from a steady state configuration because of the dynamic imbalance between spatial interactions and temporal instabilities.This dissipative effect spread thoughout a neural network because of the spatially localized connectivity between adjacent group of neurons, leading to the diffusive transmissions to neighboring cells.In fact, from experimental and theoretical standpoint, waves travelling over dissipative neural network mostly vanish due to collision [13]- [15] and lack of an external energy source to take care of attenuation [16].However, we still observe localized short excitations of nerve impulses [17], ultrashort pulses from passively mode-locked lasers, travelling waves in cortical networks, Bose-Einstein condensates in cold atoms [18], in various autonomous dissipative media, suggesting that more stable solitonic profiles can be observed provided measures are taken to minimize dissipation.
The control of pulse propagation in dissipative media can also be achieved by subjecting the system to high frequency periodic perturbations.For a discrete system like myelinated axons, this method has the advantage that the external frequency can either suppress or enhance the pulse propagation [19].The purpose of this study is to examine the evolution of a nerve impulse, responsible for carrying electrical signals along a weakly dissipative myelinated axon.This issue is brought into sharp focus because neurons effectively participate in a collective spatiotemporal sharing of vital information.For some conservative media like in optical fibers, such crucial signals are conveyed by trains of solitons [20]- [22] which are robust to external perturbations.We derived the Complex Ginzburg-Landau (CGL) equation from a weakly dissipative FHN model, where the pulse solution clearly depicts the action potential typical in neural networks.The linear stability analysis of the action potential generally shows that the nerve impulse is relatively unstable to small pertubations.
The rest of the work is organized as follows; in Section 2, after a brief introduction of the FHN model, we derived the model CGL equation in a weak dissipative medium using the multiple scale expansion in the semi-discrete approach.The analytic solution of the CGL equation is obtained in Section 3 following the method highlighted in [23], to obtain the Pereira and Stenflo [24] pulse soliton solution.A linear stability analysis is carried out in Section 4 to check the robustness of the nerve impulse when subjected to a minimal perturbation.Consequently, the Benjamin-Feir criteria is satisfied from the qualitative analysis of the modulational instability.Numerical simulations of the CGL equation depict the evolution of the pulse soliton whose width increases as the dissipation increases.In the absence of dissipation, we observe a stable profile of the nerve impulse.Finally Section 5 gives a summary of the entire work.
Here n u and n v are the membrane potential and the recovery variable respectively at the th n excitable membrane site known as Ranvier node.The recovery variable n v depicts the slow dynamics because the parameters , , 1. a b c  For Equation (1a), the first term, n u , accounts for the positive feedback, where depo- larization enhances more depolarization through the voltage-gated sodium channel.The second term, 3 3 , is a rapid negative feedback loop, corresponding to the auto-inactivation of the sodium channel.The third term n v − represents a recovery process which may be physically responsible for regulating the outward potassium currents that oppose depolarization.The last term is the discrete diffusive term with coupling strength K, that is proportional to the difference in internodal currents through a given Ranvier node.The first term of (1b) i.e. a mainly measures the potassium leakage current while n bu captures the activation of the voltage-gated po- tassium channel by the membrane potential n u , which increases the magnitude of n v .Lastly, n cv − controls the pumping of potassium ions out of the neuron.
Differentiating Equation (1a) with respect to time and substituting the value of n v  from (1b) yields where 0 0 1 2 0 , , , ,D γ γ γ Ω and 1 D are all constants.System (2) is the Liénard form of the FHN model of a myelinated nerve fiber which may be considered as the modified form of the van der Pol equation.For 0 2) becomes the cubic Liénard equation with linear damping and it is sometimes regarded as a generalization of damped oscillations.Furthermore, within the linear regime and when the dissipation is neglected, we obtain the well known linear harmonic oscillator that finds numerous applications in both classical and quantum physics.The Liénard type of equations have been intensively investigated from both mathematical and physical perspectives, and their study remains an active field of research in mathematical physics [25]- [28].
Equation ( 2) also mimics a one dimensional chain of atoms with unit mass, harmonically coupled to their nearest neighbors, characterized by dissipation and subjected to nonlinear on-site potential.A lot of difficulties is encountered to analytically solve this equation, however we will use a perturbation technique to minimize the effects of dissipation and also attempts to find appropriate solutions.In this light, all the dissipative coefficients 2 1 , D γ are perturbed to the order 2 ε , where 1 ε  is a slow variable parameter.Also, we greatly minimize the rate at which the leakage potassium ions are discharged from the neuron by perturbing 0 γ to order 3 ε .
Consequently keeping just terms up to order 2 ε , Equation ( 2) is rewritten as There are always frequency limits within which normal propagation of nerve impulse signals are observed across the FHN myelinated axon.In order to define the appropriate frequency range, we employ a perturbative technique where a suitable solution of the FHN model containing ε parameter is used.Upon substitution of this solution into the diffusive FHN model (3), the dispersion relation is obtained with terms at order of 0 e n iθ ε .Hence, we consider a simplified solution of our FHN model (3) of the form; , e , e , with n qn t θ ω = − where q is the normal mode wave number and ω is the angular frequency.In the semidiscrete approximation, ( ) ψ is supposedly independent of the "fast" variables t and n.Instead, it depend on the "slow" variables defined by mation is then made with the wave amplitudes while the discrete nature of the phase in maintained, details of this method is given in [25] [29].
Collecting terms at order 0 e n iθ ε gives the dispersion relation which is plotted in Figure 1.From Equation ( 5), the linear spectrum has a gap min 0 ω = Ω and it is limited by the cut-off frequency ( ) The variation of the diffuseness of the plasma membrane through the dispersion coefficient 0 D , physically

Ω =
. There is a lower cutoff mode 0 q = with frequency 0 0 ω = Ω and upper cutoff mode π q = with frequency ( ) accounts for the alteration of ions movement across pumps and ion channels at the Ranvier nodes.Clearly for normal propagation of action potential across the myelinated axon, there exists a frequency range 0 max ( ) (upper cutoff mode π q = ).In this frequency range, the voltage-gated sodium ions channels are highly concentrated in the Ranvier nodes which are myelin free.When an action potential is generated at the axon hillock, the influx of sodium ions causes the adjacent Ranvier nodes to depolarize, resulting in an action potential at the node.This also triggers depolarization of the next Ranvier node and the eventual initiation of an action potential.Action potentials are successively generated at neighboring Ranvier nodes; therefore, the action potential in a myelinated axon appears to jump from one node to the next, a process called saltatory conduction.For 0 f f < , the action potential initially generated from one Ranvier node jumps to the next in a discontinuous manner.This process allows the speed of conduction to be greatly increased and completely distorts switching between adjacent Ranvier nodes.Consequently, this leads to propagation failure because of the inability of the nerve impulse to move across the axon.Lastly, for max f f > , the nerve impulse can not propagate because its amplitude exponentially decays to zero.However considering the nonlinearity of the medium, it may be possible to observe normal propagation along the axon provided the amplitude exceeds a certain threshold value.This mode of propagation called supratransmission has been realized in many physical systems [30]- [33].
Terms of order 1 e n iθ ε where 0 sin is the group velocity as depicted in Figure 2.
Finally, we collect terms proportional to 2 e n iθ ε to have By considering the reference mobile frame i  the wave, we obtain ( ) where P, is the real dispersion coefficient, and r i Q Q are respectively the real and imaginary parts of the nonlinear coefficient, while 0 R < is the dissipative coefficient that causes the attenuation of action potential because of lost of energy during propagation.These parameters are given by; Equation ( 8) is the CGL equation and generally speaking, it represents one of the most-studied nonlinear equations in the physics world today.This is because it gives a qualitative and quantitative description of a myriad of physical activities [11] [25] [34].In neural networks, the propagation of modulated nerve impulses observed is governed by CGL equation, which clearly demonstrates how neurons participate in processing and sharing of information [25].In this study which we focus on the propagation of action potential along a myelinated axon, it is incumbent on us to obtain appropriate solutions that depicts the asymmetric structural features of action potential.Physically in most conservative media [35], solitary waves are understood as carriers of energy, however in this our context where dissipation is still prominent we need to balance the energy outflow in order to observe stable nerve impulses.

Solution of the Equation of Motion
Since we are dealing with a CGL equation, it is incumbent on us to look for a propagating wave solution of the form: which upon substitution into Equation ( 8) to obtain Let us assume that ( ) ( ) where d is the chirp parameter and 0 φ is an arbitrary phase which we set to zero (i.e.0 0 φ = ) without loss of generality.Equation ( 15) is obviously a constraint imposed on ( ) φ ξ because the chirp could have a more general functional dependence on 1 ξ .However this restriction allows us to find all possible analytic pulse-like solutions of the modified CGL Equation (8).Consequently, Equation ( 14) becomes We now have two second order ordinary differential equations (ODE) relative to the same dependent variable, ( ) a ξ .To have a common solution, the two equations must be compatible which generally in most systems is not the case.However, for this particular system, they can be made compatible by a proper choice of the parameters.
The following procedure is employed in order to find the compatibility conditions of the two equations; initially we eliminate the first derivatives from Equation ( 16) to have which upon integration and setting the integration constant to zero yields On the other hand, we eliminate the second derivative from Equation ( 16), obtaining ( ) Equations ( 18) and ( 19) must coincide, consequently leading to the following conditions: ( ) In order for us to obtain pulse solutions with real amplitudes, we assume that 0, 0.50, 0, , sech e , 2 0.5 with the amplitude of the solution plotted in Figure 3.
We now obtain the exact analytic solution of the nerve impulse propagating along the myelinated axon by substituting ( 21) into the FHN model solution (4) to have;  where Figure 4 depicts the spatial profile of this nerve impulse which may be considered as an electrical signal that propagates over a long distance without a change in amplitude.Also, when the parameter ε is varied, the form of the wave is greatly distorted as in Figure 5. Finally, Figure 6 depicts the spatiotemporal

Linear Stability and Numerical Analysis
Stability is a very important property of a wave profile in a neural network, since it determines whether patterns can be observed experimentally, or utilized for diagnostic purposes.Recall that the phenomenon of modulational instability results when a steady-state solution is subjected to a weak perturbation, which eventually leads to the exponential growth of its amplitude along the line of propagation due to the interplay between the nonlinearity and dispersive effects of the medium.Initially, we consider the stability of the trivial homogeneous solution , 0 ψ ξ τ = by substituting the perturbations of the form 1 2 e iqξ λτ + in the linearized part of Equation ( 8).This leads to the characteristic equation ( ) where q is a spatial wavenumber of the perturbation and λ determines the growth rate.The corresponding dispersion relation reads and when all the eigenvalues ( ) q λ have negative real parts, the homogeneous solution is asymptotically stable.Clearly the stability depends on the dissipative coefficient R and we conclude that the trivial solution is unstable since 0 R < .
We now study the existence and stability of plane wave solutions ( ) of Equation (8).We obtain the relation between the unknown amplitude 0 ψ , wavenumber q and frequency ω of the plane wave solution as follows; ( )  ( ) After substituting (28) into (29) and using Equation (25), we obtain an equation involving two linearly independent functions .Requiring that the coefficients of these functions are zero, we arrive at a system of linear equations for the unknowns ψ + and ψ − : with ( ) ( ) Solutions ( ) k λ can now be found explicitly and the maximum of their real parts determines the stability of plane waves.Consequently we solve Equation (31) to obtain ( ) ( ) where is the discriminant of the quadratic Equation (31) PQ , while they are stable for negative values.Consequently, one expects to find spatially localized nerve impulses along the myelinated axon, for any wave carrier whose wave vector is in the positive range of r PQ .The condition 0 r PQ > , is also the condition for the existence of solitons, the study of the long-term evolution of a plane wave injected as initial condition in the neural network shows that, after a period during which the waves start to self-modulate as predicted by the linear stability analysis we just performed, this change persists until the amplitude of the plane waves vanishes in some regions.Sometimes, this modulational instability leads to a train of pulses [36], but for our case, this effect generates a single nerve impulse produced by the superposition of continuous waves oscillating at a constant frequency shift.We now perform the numerical simulations of the CGL Equation ( 8) by using a Runge-Kutta scheme with fixed step size and initial conditions obtained from the analytic solution (21) of the CGL equation to check the long term effects of modulational instability.Figure 8 depicts the spatiotemporal evolution of the nerve impulse along the myelinated axon as the dissipative parameter is varied.
In the first case, i.e. Figure 8(a) where the dissipation is zero, we observe that the pulse conserves its shape, indicating that energy is not lost during propagation.This creates an ideal platform where vital information is transmitted along the myelinated axon with little or no distortion.As the dissipation increases as in Figure 8(b) & Figure 8(c), the nerve impulse becomes wider and flatter and may be termed a flat-top soliton.As time tends to infinity, the flat top solitons becomes unstable and decomposes into two fronts.
The corresponding contour plot of the spatiotemporal evolution of the nerve impulse in Figure 8 is given by Figure 9.This numerical results clearly confirms the prediction of our linear stability analysis by demonstrating that solutions for 0 r PQ > are unstable, and furnishes us with the details beyond the linear limits.This shows that the nerve impulse degenerates to fronts which are also solutions of the CGL equation [37].

Conclusions
In this study, we addressed the issue of nerve pulse propagation along a weakly dissipative myelinated axon modelled by the discrete FHN model.The effect of dissipation is always a big nuisance during the transmission of electrical signals across a neural network, consequently crucial information is always lost.We transformed the FHN model to its Liénard form and minimized the effects of dissipation by perturbing the appropriate parameters to higher order.The proposed solution for the Liénard equation carried the ε parameter, where the dispersion relation that governed the propagation of nerve impulse was obtained at the first order of terms in ε .The second and third order terms in ε gave the group velocity and the CGL equation respectively.The pulse soliton solution of the CGL equation was derived analytically, and clearly mimiced the action potential typical in neural networks.We carried out liner stability analysis in Section 4, with results showing that the pulse soliton was generated as a result of modulational instability of the plane wave solution.By numerically integrating the CGL equation and varying the dissipation, we demonstrated that nerve impulses conserved its shape in a dissipative free medium.As the myelinated axon becomes more and more dissipative, the nerve impulse  From this theoretical study carried on a dissipative myelinated axon, we believe it will furnish us with the appropriate knowledge of predicting the physiological state of real neurons.For instance, a sick neuron which is usually considered as dissipative can be clearly distinguished from a healthy neuron and consequently lead to

Figure 1 .
Figure 1.Linear dispersion curve of the nerve impulse for 0 0.04 D = and 2 0

v
being the group velocity of

Figure 7 Figure 7 .
Figure 7 shows the amplitude of the plane wave solution 0 ψ as a function of the dissipative parameter R − (since 0 R < ) for the determinant of the coefficient matrix vanishes:
unstable and degenerates to fronts.

Figure 9 .
Figure 9. Corresponding contour plot of the evolution of bright nerve impulse in Figure 8 when the dissipation is varied for 0.00064, 0.600, 0.50.
. It is a very useful component because it determines the nature of roots of the growth rate ( ) time resulting to the instability of the plane waves which tend to self-modulate with a wave vector k.The plane wave solutions of the CGL Equation (8) clearly manifest the the classical Benjamin-Feir scenario where plane waves are unstable for positive i Q >