Model of the Nerve Impulse with Account of Mechanosensory Processes: Stationary Solutions

Mechanotransduction refers to a physiological process by which mechanical forces, such as pressures exerted by ionized ﬂuids on cell membranes and tissues, can trigger excitations of electrical na-tures that play important role in the control of various sensory ( i.e. stimuli-responsive) organs and homeostasis of living organisms. In this work, the inﬂuence of mechanotransduction processes on the generic mechanism of the action potential is investigated analytically, by considering a mathematical model that consists of two coupled nonlinear partial diﬀerential equations. One of these two equations is the Korteweg-de Vries equation governing the spatio-temporal evolution of the density diﬀerence between intracellular and extracellular ﬂuids across the nerve membrane, and the other is Hodgkin-Huxley cable equation for the transmembrane voltage with a self-regulatory ( i.e. diode-type) membrane capacitance. The self-regulatory feature here refers to the assumption that membrane capacitance varies with the diﬀerence in density of ion-carrying intracellular and extracellular ﬂuids, thus ensuring an electromechanical feedback mechanism and consequently an eﬀective coupling of the two nonlinear equations. The exact one-soliton solution to the density-diﬀerence equation is obtained in terms of a pulse excitation. With the help of this exact pulse solution the Hodgkin-Huxley cable equation is shown to transform, in steady state, to a linear eigenvalue problem some bound states of which can be obtained exactly. Few of such bound-state solutions are found analytically.


Introduction
The mechanism by which the action potential is generated and transmitted along the neuronal network, is a complex phenomenon that has been actively investigated in Mathematical Physics [1][2][3][4][5][6][7][8]. A pioneer model of the nerve-impulse generation was proposed by Hodgkin and Huxley (hereafter referred to as HH model) [1]; the model rests on the picture of the nerve impulse as an electric voltage propagating in form of pulse wave along the nerve fiber. Thus according to the HH model the nerve impulse is a self-regenerative wave associated with the electrochemical activity of the nerve cells, due to flow of ion currents (Na and K) through specific ion channels [1]. This wave is assumed to propagate with a constant shape following a mechanism which can be summarized as follow: During the generation and transmission of the nerve impulse, the leading edge of the depolarization region of the nerve triggers adjacent membranes to depolarize. This causes a selfpropagation of the excitation related to the transmembrane potential down the nerve fiber [1,9,10].
Hodgkin and Huxley suggested that a convenient way to describe the propagation of this transmembrane potential, is to view the nerve fiber as an electrical transmission line. Hence in its most conventional formulation, the HH electrical model assumes currents in intracellular and extracellular fluids to be ohmic such that the net transmembrane current is the sum of ionic and capacitive currents. In this context, the conservation law for current flow across the membrane leads to the following reaction-diffusion equation [1]: where V is the transmembrane voltage (also referred to as action potential), C m is the membrane capacitance, D is the diffusion coefficient and F (V ) is the contribution from ion currents to the transmembrane voltage generation. However, besides the well-established electrical activity of the nerve membrane, there are experimental evidences [11][12][13][14][15][16] of existing mechanical forces acting on the membrane [8] also contributing to the nerve impulse generation [17]. To this last point, it is well known that at equilibrium the nerve membrane is in a liquid phase; however under mechanical forces the membrane undergoes a structural change to a gel phase resulting in an increase of its compressibility, reminiscent of structural properties of anharmonic solids [17,18]. To describe the formation of density pulses associated with mechanical constraints acting on the nerve membrane, there has been proposals to include nonlinearity in the membrane compressibility together with dispersive effects related to the ladder structure of the nerve fiber [17]. In this respect a nonlinear mathematical model including dispersions was proposed by Heimburg and Jackson [4], and later on extended by several authors (see e.g. ref. [7]). The underlying equation, considering a one-dimensional (1D) density pulse propagating along a biomembrane of a cylindrical shape, is given by [4,7]: where u is the density change across the membrane due to mechanical contraints on the membrane, c 0 is the sound speed in the homogeneous regime (i.e. in the absence of pressure waves), h and p are constants [17] (hereafter we shall set p = −1 for simplicity [7]).
Given the evident simultaneous contributions of both electrical and mechanical activities of the nerve in the generation of the action potential, a better description of the generic mechanism of the action potential requires full account of these two activities. In this work we propose a mathematical model which combines the two nonlinear partial differential equations, namely the one describing the spatio-temporal evolution of the density pulse across the membrane, and the HH cable equation with a feedback from ion flows across the nerve membrane representede by a density-dependent membrane capacitance. The full model being analytically untractable, in the present context we focus mainly on the generic mechanism of the action potential. In this regard we consider the steady-state regime of the model, where exact soliton solutions can be found. In a future work numerical simulations will be carried out on the spatio-temporal model, to check the stability of these exact soliton solutions as they propagate along the nerve.

The Model
The nerve fiber can be pictured as a cylinder with walls made from the cell membrane surrounded by intracellular and extracellular fluids [6]. The intracellular fluid stands for a conductive liquid with a high concentration of potassium ions but low concentration of sodium and chlorine ions, while the cell membrane acts like a barrier preventing ions from intracellular liquid from mixing with extracellular liquid. Due to the difference in ion concentrations between intracellular and extracellular fluids, a resting potential is expected to set up across the membrane. If the nerve is depolarized, or more precisely in the presence of a stimulus, the membrane becomes selectively permeable to ion currents flowing rapidly into the cell, reversing the polarity of the action potential.
In the picture of electric cable, the selective regulation of transmenbrane ions by the cytoplasm suggests a "management" (or selfregulation) of the charges stored in the nerve membrane capacitance. Consequently we can expect changes in the electrostatic potential of the membrane during a propagating density pulse, indicating a possibility for an electromechanical coupling between the membrane density and the electrostatic potential. Such electromechanical coupling has been widely reported in neurophysiological experients [19][20][21][22], and can be linked with changes in the membrane capacitance as a result of variation of the fuid density across the membrane.
To include the electromechanical coupling in the mathematical model for the nerve impulse generation, we reformulate the HH equation for the action potential which leads to the following new mathematical model for the nerve-impulse generation and propagation (neglecting the contribution of F (V ), for simplicity): Characteristic parameters in the above set of coupled equations are defined as follow: • Q m (u) is the total charge stored in the membrane capacitance at time t, • D is the diffusion coefficient, • c 0 is the sound speed along the membrane, • h is the fourth-order dispersion coefficient [4].
By definition, the term in the left of Equation (3) represents the current I c across the membrane capacitance C m (u) = ∂Q m (u)/∂V . This means that I c is related to the instantaneous charge in the membrane i.e. Q m = C m (u)V (assuming that the membrane capacitance is not explicitely dependent on the transmembrane voltage), through the relation I c = ∂Q m /∂t. The term in the right-hand side of Equation (3) is the diffusion current related to charge flow across the membrane (with D = a/2R, where a is the radius of the cylindric nerve and R is the membrane resistance).
We consider that at some given time t the membrane capacitance C m is determined by the difference u in densities of the fluids crossing the membrane, and propose a law of dependence of C m on u of the following form: where κ is real and constant. With Equation (5) we can transform Equation (3) into a reaction-diffusion problem with a linear source term namely: Since formula (5) suggests that the membrane capacitance C m (x, t) is fixed by the net difference u(x, t) in densities of the ionic fluids crossing the membrane at time t and position x, in the next section we solve analytically Equation (4) in order to obtain the anaytical expression of u(x, t).

Pulse Solution to the Density-Difference Equation
Equation (4) describes the motion of waves in one space dimension. Let us assume that the solutions we seek are localized coherent structures obeying: lim and propagating in the direction of positive x. Thus, using the coordinate transformations: we can rewrite Equation (4) as: where α and β are constants depending on c 0 , h and c. Note that parameters α and β can be set to any value by coordinate transformation, however the values we shall use in this work are α = 6 and β = −1 for which Equation (9)  Vries equation [23]. As we are interested in exact analytical solutions to the density-wave equation, let us first explicitely set the associate mathematical problem. Starting with an initial density pulse u 0 (x), we are interested in the initial-value problem: u(x, 0) = u 0 (x).
We look for the existence of special solutions called solitons, which are travelling waves with permanent shape obtained by means of the Inverse-Scattering Transform [23].
As a useful recall, within the framework of the Inverse-Scattering Transform any exactly integrable nonlinear evolution equation for a function u = u(x, t) can be mapped onto a pair of spectral problems so-called Lax pair [24], namely: where L and A are two linear operators the coefficients of which depend on the function u and its derivatives. Suppose: and then the associated evolution equation; turns to the following Korteweg-de Vries equation: The first step in the Inverse-Scattering Transform is to solve the forward scattering problem namely: where the spectral parameter λ ∈ R because the operator L is selfadjoint. If the potential u in (12) decays sufficiently quickly as |x| → ∞, then the linear problem (16) can be written: There are two possible kinds of eigenfunctions for problem (17): 1. bound states: If λ < 0, then we may define k = √ −λ to have: The values of k for which d = 0, so that Ψ(x) is bounded as x → +∞, are the discrete eigenvalues of L, denoted k n . The corresponding functions Ψ n = c n e −knx are bound states. It is convenient to normalize Ψ n so that: with this normalization condition the constant c n is the normalization constant.
The constant a is called the transmission coefficient and b the reflection coefficient.
The eigenvalues k, k n together with the constants c n , b and a constitute the scattering data. The time dependence of these scattering data ensues from the second auxiliary equation i.e.: in this case we have the following relations: c n (t) = c n (0)e 8k 2 n t , k n (t) = k n (0).
The inverse scattering problem is the problem of reconstructing the wave potential u on the basis of the scattering data. Given the normalization constant c, the eigenvalues −k 2 n for the bound states, and the reflection and transmission coefficients for the continuous spectra, the potential u(x, t) is recovered from: where the implicitly time-dependent function K(x, y) satisfies the Gelfand-Levitan-Marchenko equation [25,26]: in which: We now apply this procedure to the initial value problem (10), where the initial profile is taken to be: We note the following: 1. Since u(x, 0) ∈ L 2 (R), the operator L is self-adjoint and thus has a real spectrum.
With the initial potential (24), the scattering problem (17) at t = 0 becomes: where τ = tanhx. Recall that the Associated Legendre equation is of the general form [27]: where n and k are integers. The corresponding bounded solutions are given by: where: are Associated Legendre polynomials. Hence the bounded solution of Equation (25), corresponding to λ = −k 2 1 , exists for k = k 1 = 1 and the discrete eigenfunction is therefore: where: Since R sech 2 x dx = 2, the normalized eigenfunction becomes: The asymptotic behavior of this solution yields: so that: Consequently c 1 (t) = √ 2e 4t . This information is sufficient for the reconstruction of u(x, t), given that we have chosen an initial profile for which b(k) = 0 for all k. It follows that Equation (23) becomes: and the integral Equation (22) is therefore: This implies that K(x, z; t) = L(x, t)e −z , where L is such that:  Equation (31) can be solved directly to yield: Thus, which is a pulse soliton. The spatial profile of the solution to the density-difference equation is sketched in Figure 1, in dimensionless units.

The Action Potential
Having obtained the analytical expression for the density pulse u(x, t), let us now look for possible solutions to Equation (6). In this purpose we assume that the axon is a cylinder of length L, and that the nerve impulse lasts for T units of time. With these considerations, the model turns to the following boundary-values system: Equation (34) needs to be solved numerically, to gain a consistent understanding of the spatio-temporal evolution of the action potential along the axon. Nonetheless there exists a regime in which exact analytical solutions can be obtained. Focusing on this specific regime, namely the steady-state regime where the density pulse u(x, t) ≡ u 0 (x), the action potential can be reduced to its stationary form v(x) and as a sequel Equation (34) turns to: This is a linear-operator equation with zero eigenvalue in which µ = 2κ D . Equation (38) is actually similar to Equation (25), and for this reason we can apply the same variable change to find an Associated Legendre equation but now with n(n + 1) = µ and k = 0. For this case, bounded solutions assuming arbitrary values of the integer n are the Legendre polynomials: v n (τ ) = 1 2 n n! d n dτ n (τ 2 − 1) n , n = 1, 2, 3, · · · . (39) Exact analytical expressions for the first three bound states are given below: The three bound states are sketched in (Figure 2). Note that these bound states are not the only possible solutions to the steady-state equation for the action potential, other possible solutions can be obtained for different values of the parameter µ = 2ε D linking the diffusion coefficient D with the coefficient κ coupling the two equations.
It is worthwhile to stress here that the steady-state solutions obtained analytically, are exact initial profiles of the action potential and hence are expected to play a significant role in numerical simulations of Equation (34), in order to find exact spatio-temporal profiles for the propagating action potential V (x, t). Owing to the extremely rich variety of possible steady-state profiles for the action potential, represented by distinct Legendre polynomials v n (τ ) for different values of µ = 2ε D = n(n + 1), to gain a consistent insight from numerical simulations it will be useful to consider a reasonable number of bound states. Therefore numerical simulations of the model deserve a specific framework of study and will be considered in a future work.

Conclusion
We have introduced a mathematical model to describe the influence of mechanosensory processes, on the generation and propagation of the nerve impulse. The model rests on a picture of the axon having the form of a long cylinder, with the membrane regarded as a capacitive diode for which the capacitance is constantly adjusted by the difference in densities of intracellular and extracellular ionic fluids, on either side of the membrane. Mathematically the model consists of a Korteweg-de Vries equation for the ionic density exchange across the nerve membrane, coupled to the HH cable equation with a variable capacitance. While the full model turns out to be analytically untractable, we obtained that its stationary regime represents a very rich source of various possible exact soliton solutions. Namely we obtained a single-pulse soliton for the density difference, and several distinct possible solitary-wave profiles for the action potential including a pulse, a kink, a pulse-in-kink structure and so on. Given their large numbers, the issue of their propagation requires solving numerically the corresponding spatio-temporal equation. This numerical study must take into consideration several of these profiles, which requires a specific framework of study. This last aspect will be considered in a separate work.