On the Hybrid Model of Nerve Pulse: Mathematical Analysis and Numerical Results ()
1. Introduction
The mechanisms and properties for Nerve pulse generation and propagation form one of the most investigated physiological activities in Neuroscience [1] - [11] . The modelling of this process is historically related to the analysis of action potentials (pulse like voltage waves that carry information along a nerve fiber) whose measurements by electrophysiological techniques have well documented references. The Hodgkin-Huxley model [1] was one of the earliest mathematical models that described and captured the propagation of an asymmetric pulse in the nerve fiber. Hodgkin and Huxley showed how measurements of the conductive parameters of a nerve fiber can be used to directly calculate both the shape and velocity of an action potential on a Squid’s giant axon. According to their model, the nerve pulse is a self-regenerative wave associated with the electrochemical activity of the nerve cell and due to the flow of Sodium and Potassium ions through the protein pores on the nerve membrane. It is now well established that during the generation and transmission of the nerve impulse, the leading edge of the depolarization region triggers adjacent membranes to depolarize causing a self-propagation of the excitation, related to the transmembrane potential, down the nerve fiber [1] [12] [13] . As suggested by Hodgkin-Huxley, a convenient way to describe the nerve impulse propagation is to think of the nerve fiber as an electrical transmission line. Thus, in its most conventional formulation, the Hodgkin-Huxley electrical model assumes currents in intracellular and extracellular fluids to be ohmic so that the net transmembrane current is the sum of ionic (due to the flow of Sodium and potassium ions) and capacitive currents (due to the charge stored on the membrane capacitor). The conservation law for current flow across the membrane was modelled by the equation [1] :
(1)
where V is the transmembrane voltage (action potential), Cm the membrane capacitance, D the diffusion coefficient and F is a contribution of ion currents.
Besides the well-established electrical activity of the nerve membrane, there are experimental evidences [14] - [19] of existing mechanical forces acting on the membrane [21] that are also contributing to the nerve impulse generation [20] . To this last point, it is assumed that the fluids surrounding the nerve cytoplasm induce a sound-like wave when crossing the membrane (Osmotic flow across lipid membranes) and the feedback process governing this flow is represented by a nonlinear dependence of the sound on the net density difference of fluids across the membrane [7] . To describe the formation of density pulses associated with mechanical constraints acting on the nerve membrane, there have been proposals to include nonlinearity in the membrane compressibility together with dispersive effects related to the ladder structure of the nerve fiber [20] . In this respect a nonlinear mathematical model including dispersions was proposed by Heimburg and Jackson [7] , and later on extended by several authors (see e.g. ref. [11] ). The underlying equation, considering a one-dimensional density pulse propagating along a biomembrane of a cylindrical shape, has been modelled by [7] [11] :
(2)
where u is the density difference between fluids across the membrane, c0 is the speed of sound in the homogeneous regime (i.e. in the absence of pressure waves), h and p are constants [20] (hereafter we shall set
for simplicity, [11] ).
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 a full account of these two activities.
In [22] , 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 H-H cable equation with a feedback from ion flows across the nerve membrane represented by a density-dependent membrane capacitance was proposed. Properties such as the steady-state regime of the model were considered and exact soliton solutions determined.
In the current paper we consider the non-steady state regime of the model proposed in [22] from a more analytic and mathematical stand point. We establish an appropriate mathematical framework to handle such complex equations and, formulate and prove conditions under which such an initial value problem will admit a unique solution. We then use the results from the framework to verify that the modified model in the non steady state will have a unique solution. Finally, a numerical experiment is designed to estimate the solutions in the non steady state regime and check the stability of these solutions as they propagate in the x-t plane.
2. The Model
In [22] , Mengnjo, Dikandé and Ngwa formulated a model of the nerve impulse taking into account the mechanosensory processes of the nerve cell. In their formulation, the influence of mechanotransduction processes on the generic mechanism of the action potential is investigated analytically by translating to a mathematical model that consists of two coupled nonlinear partial differential equations. The coupling equations are: the Korteweg de Vries equation governing the spatio-temporal evolution of the density difference between intracellular and extracellular fluids across the nerve membrane, and a modified Hodgkin-Huxley cable equation for the transmembrane voltage with a diode-type membrane capacitance. The membrane capacitance is assumed to vary with the difference in density of ion carrying intracellular and extracellular fluids thus ensuring an electromechanical feedback mechanism and consequently an effective coupling of the two nonlinear equations. The equations in [22] are:
(3)
(4)
(5)
With the help of the exact one-soliton solution,
, of the density difference equation, Equation (4), the cable equation, Equation (3), becomes
(6)
where V is the transmembrane voltage. Taking
, the first three bound states of the model (3)-(5) are [22] :
(7)
(8)
(9)
Here, we address the issue of propagation of these bound states by solving, numerically, the spatio-temporal Equation (6). To attempt such a numerical experiment, we have to ensure that the equation is mathematically well-posed (i.e. the solution exists, is unique and depends continuously on the initial data
). We remark that the nerve impulse lasts for a finite time T and the length, L, of the axon is also finite. Consequently the time and the space domains will be considered as bounded intervals.
3. Mathematical Analysis of the Model
In this section, we state and prove conditions under which an initial value problem of the form
(10)
(11)
with homogenous boundary conditions
will have a unique solution. Non-homogeneous boundary conditions can be reduced to homogenous ones by homogenization depending on their smoothness.
3.1. Sobolev Spaces of Solutions
We now introduce some of the function spaces that form the basis of the mathematical framework for the analysis of an initial value problem of the form (10-11) [23] .
1) Given an interval I and the Hilbert space X,
will denote the space of bounded k-continuously differentiable functions of the form
equipped with the norm
For example the space
is the space of bounded continuous functions of the form
with the norm
2) The space
will denote the space of measurable functions u such that
and
with the norm
3)The space
equipped with the norm
. It is also the closure of the space
of infinitely continuously differentiable functions with compact support in
.
4) The space
is the dual space of
in the sense of distributions.
5) The space
is the completion of
with respect to the norm
6) The space
with the norm
where
7) The space
of infinitely continuously differentiable functions from
to
which have a compact support on
.
3.2. Existence and Uniqueness of Solutions
Let us consider the operator
defined by
The system (10-11) becomes
(12)
(13)
(14)
Lemma 1. Define
on
. Then
satisfies the coercivity condition
(15)
where
are positive constants.
Proof. Let
. Since
and
we see that f and g are positive and bounded (as continuous functions on bounded domains). Then there exists
such that
and
for any
. We have that
(16)
(17)
(18)
(19)
From the Poincare inequality,
where C is a positive constant. Consequently
where
and
.
To continue our way to the existence and uniqueness, we need to establish a relation between the space of solutions and the space of continuous functions.
Lemma 2. Suppose
, then
.
Proof
Let
,
then
(20)
(21)
(22)
Consequently
So if we choose
in such a way that
is the mean value of
in
, we have
Then
Hence
We can now prove the existence and uniqueness of solutions of problem (10-11).
Theorem 1. Assume that the function
and
are given. Then (10-11) has a unique solution
Proof.
Set
where
is the constant defined in inequality (15). Then
and the Equation (10) becomes
Hence, the system (10-11) is equivalent to
(23)
(24)
(25)
We will show that (23-25) has a unique solution and conclude the existence and uniqueness of the solutions of (10-11). Let us set
Then from the proof of Lemma 1,
and so
Consequently, from the coercivity condition (15), Poincare inequality and integration by parts,
(26)
Proof of uniqueness
Let u be a solution of (23-25) then we have
and
Consequently, using (26)
and
since
is compactly embedded in
. Therefore,
So if
and
are two solutions of (23) then
then
.
Proof of existence
is separable and so has a countable dense subset
. From the Gram-Schmidt orthogonalization process we can assume that the
are linearly independent such that
Let
be the orthogonal projection from
on to
and let
be the solution of
(27)
The system (27) is a system of linear ODEs for the coefficients
and so has a unique solution as a consequence of Picard-Lindelof Theorem [23] .
Using a similar proof as in Lemma 1 one has
and
(28)
since (by the linearity and continuity of
,
). Then
is bounded in
and consequently has a subsequence still denoted
that converges weakly in
to a limit u. Let
be of the form
(29)
for some
where
. For
, one has
so
Taking the limit one has
(30)
Since the functions
are dense in
, (30) holds in the sense of distributions. This implies
From (30) we can also see that
. Hence
and
from Lemma 1.
Let us now prove that
. Let
such that
. Functions of the form (29) are dense in
. So if
has the form (29) and
, one has in (27) (using integration by parts over t)
Taking the limit, one has
(31)
Multiplying (23) by
and taking the integral over
yields
(32)
Equations (31) and (32) imply that
.
Remark. Considering the model (6), we observe that the function
and each of the bound states
Hence, by Theorem 1, the initial value problem (6) has a unique solution
4. Numerical Experiments
This section investigates numerically the behavior of the model (6) using the bound states
in (7-9) as the initial states of the system, the two points boundary value problem is
(33)
(34)
(35)
The smoothness of the domain and the data allow us to use the Forward and Backward Euler Finite difference methods.
The Finite Difference Method (FDM) aims at approximating straight forwardly the value of the solution of a given problem at a chosen set of discrete points of the domain. Suppose
is partitionned into
subintervals
,
of length
with
and similarly
is partitionned into
subintervals
,
of length
with
. Thus
is approximated by the set of points
,
.
4.1. The Forward Euler Finite Difference Method
Setting
and
,
, using the difference quotients
and
, Problem (33-35) can be approximated by the discrete Problem
(36)
The convergence and accuracy of the method require its stability and consistency. The stability is proved in the subsequent lemmas.
Lemma 3. Let U and V be arbitrarily discrete functions defined on G. Define
. If
for any
,
, then
Proof. The proof is done by induction on j. For
,
.
Suppose for
fixed,
. Let us show that
Hence by induction, for any
,
,
.
Define the operator
where
.
The first Equation (33) can be written in the form
in G. It follows using the existence and uniqueness of the solution that
is invertible. Knowing that the fuctions f and g are smooth and bounded, we assume the following conditions for stability.
(37)
(38)
The step size
is suffciently small so that
. (39)
Based on Assumptions (37-39), the following Lemma holds.
Lemma 4. Suppose that Assumptions (37-39) hold and consider the operator
where Id is the identity operator from
to
. There exist
such that
.
Proof. Let V be in
, using the definition of
one has
from Assumptions (37-38)
from Assumption (39)
Hence
and we can apply Lemma 3 to see that
Hence
.
Lemma 4 shows that the sequence
is uniformly bounded, consequently the FDM is stable. It remains to show that the FDM is consistent. We can see from the Taylor’s formula that for any function V twice differentiable with respect to x and differentiable with respect to t
We have the following convergence theorem.
Theorem 2. Let V be the solution of (33-35) and W be the solution of (36). Then
as
and
.
Proof. From Lemma 3 one has
The Forward Euler FDM is consistent of order 1 and stable with the stability conditions (37 - 39). The stability and consistency imply the convergence of the FDM with an accuracy
.
4.2. The Backward Euler Finite Difference Method
Setting
and
,
, using the difference quotients
and
, problem (33-35) can be approximated by the discrete Problem
(40)
As in the previous subsection, it will be shown that the Backward Euler FDM is stable and consistent then convergent. Define the operator
where
. We suppose the following assumptions:
(41)
The step size
is suffciently small so that
(42)
Based on these assumptions one uses the expression of to obtain
then
.
Setting
, it follows that
and then
. Hence
where
and
. It can be concluded that
and that the sequence
is uniformly bounded. Consequently the FDM is stable. One can still use the Taylor formula to show that the method is consistent of order
. The following theorem show the convergence.
Theorem 3. Let V be the solution of (33-35) and W be the solution of (40). Then
as
and
.
Proof.
The Backward Euler FDM is consistent of order 1
and converges under the conditions (41-42).
4.3. Example 1
This example investigates numerically the behavior of the model (6) using the bound states
in (33-35) as the initial states of the system for
. The stability conditions (37-42) are respected by the respective FDM used.
Table 1 and Table 2 present the error estimates for the Forward and the backward Euler FDM using the maximum error of two successive approximated solutions
and of two successive refinements
and
using the formula
The consistency order
is computed using the formula
The following figures are obtained for a space mesh spacing
and a time mesh spacing
to respect the stability condition. Results in Figures 1-3 depicted below showing that the bound states do actually propagate in the x-t plane.
4.4. Example 2
We consider the same problem (33-35) with
varying in the set
using a backward Euler Finite difference method, the following figures are obtained for a space mesh spacing
and a time mesh spacing
. Results in Figures 5-7 depicted below are still showing that the bound states do actually propagate in the x-t plane. Figure 4 shows the behaviour of the approximated solutions at
.
Figure 1. Sketches of the first bound mode for the action potential, in the steady state and its evolution in the x-t plane.
Figure 2. Sketches of the second bound mode for the action potential, in the steady state and its evolution in the x-t plane.
Table 1. Error estimates of the Forward Euler FDM.
Table 2. Error estimates of the Backward Euler FDM.
Figure 3. Sketches of the third bound mode for the action potential, in the steady state and its evolution in the x-t plane.
Figure 4. Bound mode: V1, Bound mode: V2 Bound mode: V3.
Figure 5. Bound mode: V1,
Bound mode: V1,
Bound mode: V1,
.
Figure 6. Bound mode: V2,
Bound mode: V2,
Bound mode: V2,
.
Figure 7. Bound mode: V3,
Bound mode: V3,
Bound mode: V3,
.
5. Conclusions
We have considered the hybrid model for nerve pulse generation and propagation proposed by Mengnjo, Dikandé and Ngwa. This model consists of a system of coupled partial differential equations. By considering the steady states of the system, we formulated an initial value problem that governs the evolution of the action potential. This is a linear parabolic equation with variable coefficients and an initial profile. To analyse the system, we formulated and proved, in general, the conditions under which such an initial value problem will have a unique solution and verified that the initial value problem governing the evolution of the action potential satisfies these conditions and hence has a unique solution. Due to the complex nature of the equation (variable coefficients in both independent variables), there are no known analytical techniques that can be used to solve the equation. Because of this difficulty we designed a numerical experiment to simulate the approximate solutions and thus gain insight in to their long term behaviour. The numerical experiments reveal that the solutions are bounded pulses propagating in the x-t plane with constant velocity and shape similar to that of the initial profile. We interpret this to mean that the initial profile which is a steady state solution to the MDN model evolves by translation without attenuation. This observation falls in line with what should be expected of a healthy nerve cell: Generate a pulse and propagate it without distortion.
In our model, we have combined both electrical and mechanical processes to generate the nerve pulse. We have demonstrated analytically that the propagation mechanism of the action potential is actually an interplay between the electrical and thermodynamic processes of the nerve cell. Considering the fact that the thermodynamic processes can be perturbed by both internal and external factors, we can subject our model to such factors and observe its behaviour. This opens up possible applications in the modelling of pathologies of nervous pathways. This could lead to a breakthrough in understanding some nervous disorders and more effective treatments designed. We thus see a major application of our work in the field of medicine.