Stability of coupled solitary wave in biomembranes and nerves

G. Fongang Achu and F. M. Moukam Kakmeni∗

Complex Systems and Theoretical Biology Group (CoSTBiG),
Laboratory of Research on Advanced Materials and Nonlinear Science (LaRAMaNS),
Department of Physics, Faculty of Science, University of Buea, PO Box 63, Buea, Cameroon
(Dated: November 15, 2019)

In this work, we consider the electromechanical density pulse as a coupled solitary waves repre-
sented by a longitudinal compression wave and an out-of-plane transversal wave (i.e., perpendicular
to the membrane surface). We analyzed using, the variational approach, the characteristics of the
coupled solitary waves in the presence of damping within the framework of coupled nonlinear Burger-
arXiv:1911.05993v1 [nlin.PS] 14 Nov 2019

Korteweg-de Vries-Benjamin-Bona-Mahony (BKdV-BBM) equation. It is shown that, the inertia

parameter increases the stability of coupled solitary waves while the damping parameter decreases
it. Moreover, the presence of damping term induces a discontinuity of stable regions in the inertia-
speed parameter space, appearing in he form of an island of points. Bell shape and solitary-shock
like wave profiles were obtained by varying the propagation speed and their linear stability spec-
trum computed. It is shown that bell shape solitary wave exhibit bound state eigenvalue spectrum,
therefore stable. On the other hand, the solitary-shock like wave profiles exhibit unbound state
eigenvalue spectrum and are therefore generally unstable.

I. INTRODUCTION in the form of a soliton, or sound wave, along with the

lipid bilayer [1, 2, 4, 11]. According to this theory, the
Nerve cells are encapsulated by a plasma membrane; compression of a lipid membrane will change its density
a thin quasi-two-dimensional layer consisting mainly of resulting to phase transitions from a liquid state to a
lipids and proteins. Lipid membrane is integral parts of gel state accompanied by the propagation of a density
every living cell; attributing an important role in signal- pulse in the gel state. Thus in the soliton theory, the
ing integration to the lipid membrane. They have been transmission, storage, and information processing are in-
reported to resist pressure, stretching, tension, and bend- trinsic properties of the lipid bilayer. The soliton model
ing. It has been proposed that these perturbations, in is based on the propagation of a localized density wave in
the form of a density or voltage pulses, may play a ma- the axonal membrane [4, 5]. The important requirement
jor role in inter or intracellular communications [1–3] as of the model is the empirically known lipid phase transi-
well as nerve pulse propagations [4, 5]. Signal transduc- tions slightly below the physiological temperatures. The
tion across the plasma membrane via various receptors soliton model has been successful in predicting the exact
and ion channels has received much attention. In a now- pulse propagation velocities in myelinated nerves. The
classic model proposed by Hodgkin and Huxley, the ionic propagation velocities are closely related to the lateral
hypothesis is used to explain the generation and propaga- sound velocities in the nerve membrane [4, 5]. More-
tion of action potential. According to the ionic hypothe- over, the soliton model explains the reversible temper-
sis, voltage-gated ions channels are seen as being respon- ature and heat exchanges observed in connection with
sible for action potential propagation along the plasma the nerve pulse. In the soliton theory, the appearance of
membrane of the neuronal axon [6–8]. The Hodgkin- a voltage pulse is explained to be a consequence of the
Huxley model is an electrical model based on two con- piezoelectric nature of partially charged and asymmetric
ductors (cytosol and extracellular space) separated by a cell membrane [4, 5].
capacitor (plasma membrane) with ions specific channels
(proteins). This model is the currently accepted model, A great effort is being devoted to the deciphering of
the mechanism of signal propagation in the axon, in par-
substantiated by the discovery of ion channels [9] and
their crystallization [10]. However, the ionic hypothe- ticular, the coupling between electrical and mechanical
waves in the form of the longitudinal compression waves
sis has been unable to explain the non-electric phenom-
ena observed during the propagation of action potential: [1–5, 11–19]. Sound propagates well in lipid membranes
in the form of the longitudinal compression wave. What
swelling of the axon, phase transition, shortening of the
axon, and adiabatic nature of action potential. A ther- is not is not obvious, though, is whether the out-of-plane
modynamic theory was therefore proposed by Heimburg transversal mode can propagates in disordered systems
and Jackson in 2005 to address some these issues [4]. A such as biomembrane and nerves. Many experiments
thermodynamic or soliton theory of nerve pulse propa- have been carried out to detect transverse wave prop-
gation hypothesizes that the action potential propagates agation in the lipid membrane. As early as 1988, Vo-
gel and Möbius investigated surface density fluctuations
in lipid monolayers at air/water interfaces by superpos-
ing two transverse capillary waves under resonance such
∗ Corresponding Author, Email: [email protected] that their two wave vectors are orthogonal [20]. They
 2

observed that the surface density fluctuations propagate II. THE COUPLED SOLITARY WAVES MODEL
longitudinally along the lipid monolayers. In addition, FOR NERVE IMPULSE TRANSMISSION.
the wavelength of the transverse mode of the transverse
capillary waves was found to be modulated by the lon- The propagation of action potential in the soliton
gitudinal mode of the wave at resonance. This results model for biomembranes and nerve is usually considered
clearly indicate the inherent coupling between the longi- as a longitudinal compression pulse. This is done by as-
tudinal and the transverse mode. In the soliton theory suming that the diameter of the axon is smaller than
for the action potential, a two-dimensional sound wave the length of the nerve pulse. However, recent studies
propagates along the membrane plane, described by lat- show that one need not only to consider the longitudi-
eral compressibility and a correlated change in density. nal compressional dislocations but also the out-of-plane
When applied to the axon, the action potential is consid- transversal mode (i.e., the mode perpendicular to the
ered as a one-dimensional represented by a longitudinal membrane surface) The propagation of an action poten-
compression pulse since the diameter of the axon is as- tial in the soliton model for biomembranes and nerve
sumed to be smaller than the length of the nerve pulse is usually considered as a one-dimensional represented
[4, 13–19]. In contrast, El Hardy and Machta consider by a longitudinal compression pulse; by assuming that
not only the longitudinal compressional dislocations but the diameter of the axon is smaller than the length of
also an out-of-plane transversal mode (i.e., perpendicular the nerve pulse. However, recent studies show that one
to the membrane surface) that Heimburg and Jackson did need not only to consider the longitudinal compressional
not considered. This means that El Hardy and Matcha dislocations but also the out-of-plane transversal mode
effectively consider a three-dimensional wave propagation (i.e., the mode perpendicular to the membrane surface).
along the axon. That is, two in a plane and one out of a By considering the nerve signal U as an electromechan-
plane and assuming the axon to display a circular cross- ical density pulse with two components, the improved
section as in the case of the soliton model, they get rid Heimburg-Jackson model equation is rewritten as
of one dimension [21]. Thus, the electromechanical den-
∂2U 2 
 
∂  c20 + αU + βU.U . ∂U + ν ∂  ∂U 

sity pulses in the nerve can effectively be considered as

=
a vector soliton with two components. That is, the lon-
 
∂t2 ∂x ∂x ∂x2 ∂t
gitudinal component corresponding relative height field, ∂4U ∂4U
and the transverse components corresponding to lateral −η1 4 + η2 2 2 . (1)
∂x ∂x ∂t
stretch field, and have been observed during nerve pulse
propagation in garfish olfactory nerve, squid giant axon, In Eq. (1), x is the spatial coordinate along the mem-
and hippocampal neuron [21]. brane cylinder and t is time. U describes a two compo-
nent electromechanical density pulse define as
Stimulation of peripheral nerves can be used for the U = u(x, t)er + v(x, t)eθ , (2)
treatment of the dysfunction of the lower urinary tract,
chronic pain, epilepsy, and other neurological disorders where, u(x, t) and v(x, t) correspond to the longitudinal
[22]. It is worth noting in this connection that the com- and transverse components respectively. er is the unit
position of the phospholipid bilayer of neurons in terms vector in the direction of u(x, t) while eθ is the unit vector
of its constitutive lipid has been strongly implicated in in the direction of v(x, t). er and eθ satisfy the orthonor-
the the onset of neurobiological and psychiatric disor- mal properties er .eθ = 0, eθ .eθ = 1, and er .er = 1. ν is
ders [23]. It will be of great interest to understand the the damping coefficient while η2 and η1 are the dispersion
characteristics of signal propagation in a lipid bilayer; a parameters.
“device” for memory storage and information process- The second dispersion coefficient η2 is physically re-
ing, as it might play a role in the functional success lated to the inertia effects of the membrane structure
of neurotherapy. As a consequence, a right mathemat- and have been shown to influence the width of the soli-
ical model of pulse propagation in the lipid membrane tary pulse. In this improved model, η1 determines the
and nerve is therefore required. In this work, we have limiting velocity at high frequencies and η2 governs how
considered the coupled solitary wave model for action this limit is reached [16, 17] and in [19] it was estab-
potential derived from the improved soliton model of lished that η2 increases the stability the nerve signals. A
nerve. We investigated in section II, using the varia- great deal of work has been done on the effects of damp-
tional approach the characteristics of a coupled solitary ing coefficient ν on the dynamics of nerve pulse. It has
waves in the presence of damping within the framework of been established that ν causes modulation instability in
coupled nonlinear Burger-Korteweg-de Vries-Benjamin- the nerve, and possible generation of localized periodic
Bona-Mahony equations (BKdV-BBM) derived from the wave trains [18]. Such localized wave trains were ana-
vector Boussinesq equation describing the dynamics of lytically and numerically obtained in [19]. The parame-
two components nerve impulse in the soliton model for ters α and β describe the nonlinear elastic properties of
biomembranes and nerves. In the last section, the stabil- membranes. At a temperature slightly above the melt-
ity of the coupled solitary wave solution is studied both ing transition, the lipid membrane has negative values for
analytically and numerically. the parameter α and positive values for the parameter β
 3

and in low-frequency limit for dipalmitoyl phosphatidyl- where ǫ and ǫ3 are chosen in such a way to balance the
choline (DPPC) membranes, the initial speed of sound effects of nonlinearity and damping, and substituting into
c2 c2
c0 = 176.6 m/s, α = −16.6 ρA0 , and β = 79.5 (ρA0)2 with Eq. (1), terms of order ǫ4 give
0 0
initial density of the membrane ρA0 = 4.035 × 10
−3
gm−2 .
If we assume a small oscillation of the density pulse,
U → ǫU and the system parameters co → ǫc0 , η2 →
ǫ−2 η2 , ν → ǫν, and the transformation

y = ǫ(x − t), s = ǫ3 t, (3)

∂u 1 2 ∂u η1 ∂ 3 u ∂3u µ ∂ 2 u β 2 ∂u
+ (c0 + βu2 ) − 3
− η2 2 − + v = 0, (4)
∂s 2 ∂y 2 ∂y ∂y ∂s 2 ∂y 2 2 ∂y

∂v 1 2 ∂v η1 ∂ 3 v ∂3v µ ∂2v β ∂v
+ (c0 + βv 2 ) − 3
− η2 2 − 2
+ u2 = 0. (5)
∂s 2 ∂y 2 ∂y ∂y ∂s 2 ∂y 2 ∂y

The nonlinear coupled Burger-KdV-Benjamin-Bona- ence as fluid dynamics [25, 26], plasma physics [27, 28],
Mahony equations (4) and (5) describe the dynamics and in the recent wave of activity on the dynamics of
of two coupled solitary waves in the improved soliton trapped Bose condensates, as well as in the the general
model of nerves. It can be used to model the dynamics context of the Hamiltonian systems [29]. Here however,
of the mechanical fields (i.e., relative height field and the we are interested in the solitary wave solution of the
lateral stretch field) that has been observed in a neuron. coupled system (4) and (5) via variational method. To
It is instructive to note that, a linearly coupled KdV proceed we first applied traveling wave ansatz method
equations have been used to model two layer settings in and obtained a set of coupled modified Lienard’s equa-
various physical systems, such as stratified fluids with tions. It is worth noting that, the modified coupled Lien-
a superimposed shear flow, and in dual-core optical ards equation has been used to model the dynamics of
waveguides carrying ultrashort pulses. In such a system, two coupled bursting neurons by considering a set cou-
two types of vector soliton solutions (i.e., symmetric pled Hindmarsh-Rose (HR) neuron model subjected to
and asymmetric solitons) are usually obtained using an external periodic excitation. Numerical simulations
the variation approximation method [31]. In addition, revealed the existence of some bifurcation structures in-
the coupled KdV-Benjamin-Bona-Mahony equation cluding saddle-nodes, symmetry breaking and period-
have been used to model the motion of small-amplitude doubling route to chaos [30]. It is imperative to note that
long waves on the surface of an ideal fluid under the the addition of the frictional term to the classical soliton
gravity force, and in situations where the motion model gives an additional term to the Lagrangian of the
is sensibly two dimensional. In such a system, one system; the so-called frictional potential term. In order
of the component models the elevation of the fluid to obtain the Lagrangian density with the damping term,
surface from the equilibrium position while the other the conventional Lagrangian formalism is inconsistent.
component models the horizontal velocity in the flow [32]. Therefore we employed the differential approach which
consists of multiplying the conventional Lagrangian with
an exponential factor [31, 40–44].
Proceeding as described above, we first seek travelling
wave solution of the coupled system (4) and (5) by us-
III. VARIATIONAL FORMALISM AND
ANALYTICAL SOLUTION TO THE COUPLED
ing the traveling-wave ansatz, u = u(y − ξs) = u(z) and
MODEL EQUATION v = v(y − ξs) = v(z) (where ξ is the solitary wave veloc-
ity). Substituting into Eqs. (4) and (5), we obtained the
coupled Lienard’s equation
The understanding of physical, and, in particular, non-
linear, phenomena through conservation laws and varia-
tional formalism is probably the most fundamental and ∂2u γ ∂u c2o − 2ξ β 3 β
+ + u+ u + uv 2 = 0, (6)
universal approaches to theoretical physics [25]. In this ∂z 2 γo ∂z 2γo 6γo 2γo
context the natural question is what kind of information
can we potentially get from knowing the energy carried
by the solitary waves in the nerve fiber?. This approach ∂2v γ ∂v c2o − 2ξ β 3 β
is widely acknowledged in such branches of nonlinear sci- + + v+ v + vu2 = 0, (7)
∂z 2 γo ∂z 2γo 6γo 2γo
 4

2ξη2 −η1
where γ = −ν 2 and γo = 2 . In the variational Lagrangian L, given by
Lagrangian approach, Eqs. (6) and (7) are restated as a Z ∞
′ ′
variational problem in terms of the Lagrangian L given L= L dz, (10a)
by −∞

 3γ − 2 3γo2 (2ξ − c2o ) − γ 2 

 
c2o − 2ξ  2 = +  (A2 + B 2 )
 
′ 1 2 γ
3γ 6γγo2 a2

2
 e γo z
2 

L =  uz + vz − u +v 


2 4γo  
γπ γπ

β  4 β 2 2  γγ z
 × cosec
u + v4 + 2γo 2γo a

−  u v e o . (8)
24γo 2γo  2
4γo2 

β  γ
γ −  −  (A4 + B 4 + 12A2 B 2 )

In Eq. (8), the exponential term e γo z denotes the damp- 72γo3 a4 a2

ing characteristics of the coupled solitary waves. In order
 
γπ γπ
to analyze the dynamics of the coupled solitary waves, we × cosec . (10b)
2γo 2γo a
use a reduced variational principle and assume solitary
wave ansatz as a trial function. Thus we choose an ansatz It should be noted that, the values of the amplitudes
of the form A and B, and the inverse width a, for the soliton so-
lutions are real, and can be gotten in terms of the sys-
(u, v) = (A, B)sech(az), (9) tem parameters by using the Euler-Lagrangian equations
∂L ∂L ∂L
∂A = ∂B = ∂a = 0. Thus, we obtain the system of equa-
and substituting this into Eq. (8), we obtain the effective tions

3γ − 2 3γ02 (2ξ − c20 ) − γ 2

 2
4γo2  2

β γ
A + 6B 2 = 0,

+ 2 2
− 3 4
− 2 (11a)
3γ 6γγ0 a 36γ0 a a
3γ − 2 3γ02 (2ξ − c20 ) − γ 2
 2
4γo2  2

β γ
B + 6A2 = 0,

+ 2 2
− 3 4
− 2 (11b)
3γ 6γγ0 a 36γ a a
 2 2 2
  2 0 2
3γ − 2 γ − 3γ 0 (2ξ − c )
0   2 β 3γ 4γ
 A + B2 − − 2o A4 + B 4 + 12A2 B 2 = 0.
 
+ (11c)



3γ 2
6γγ0 a 2 72γ0 3 a 4 a

Note that, in order to obtain Eq. (11c), we assumed a Now let turn to the analysis of the amplitudes (A, B)
small
n oscillation in the inverse pulse width a, such that of the soliton as depicted by Eq. (12). Since (A, B) and
a2 must be real for a pulse soliton, then from Eq. (12) ,
o
cot 2γγπ0 a ≈ 2γγπ0 a . Solving simultaneously for A and B
it easy to deduce that a2 is bounded within the region
in Eqs. (11a) and (11b) we obtain ν2 2 c2o −6ξ
s 4(2η2 ξ−η1 ) < a < 4+3ν . Since ν > 0, it follows that
6γo a2 (2γo2 a2 {3γ − 2} + 3γo {2ξ − c2o }) the propagation speed ξ is also bounded in the region
(A, B) = ± . (12) η1 c2o 2
7βγ {γ 2 − 4γo2 a2 } 2η2 < ξ < 6 . Figure 1 shows the evolution a against

Substituting Eqs. (12) into Eq. (11c), and solving the ξ. In Fig. 1(i), solitary pulse are stable for ξ > 800, and
resulting equation for a2 , we obtain two roots given by unstable for ξ < 800 for the given membrane parameters.
s On the other hand, Fig.1 (ii) and Fig. 1(iii) shows the
a21 − 4a2 a0 effect of the inertia parameter η2 and damping parameter
a2 = −a1 + , (13a) ν on the stability of the soliton pulse. The plot clearly
a0
s shows that increasing η2 increases the stability of the
a21 − 4a2 a0 soliton pulse while increasing ν decreases its stability.
a2 = −a1 − , (13b) As a consequence, stable soliton pulse solutions describe
a0
by Eq. (9) is possible if there is a balance between
where damping and inertia effects of the lipid membrane.
ao = −8γo4 (3γ − 2), (14a)
In order to further understand the dynamics of the
a1 = −2γo2γ 2 (3γ
− 2) + 24γo4 (2ξ − γo2 ) − 4γ 2 , (14b) solutions with respect to the pulse width, we turn our
and a2 = 4γ 4 − 15γo2 γ 2 (2ξ − c2o ). (14c) attention to the parameter space plot. First, we find
the bifurcation points by setting a2 (η2 , ξ, ν) = 0 and we
qNote that the quantity a1 > 0, and that
obtain the fixed points in terms of the parameters η2
a21 −4a2 a0
a0 << a1 . As a consequence, the two roots of and ν by solving for ξ. We then plot the 3D surface
2
a are approximately equal. a20 = F (η2 , ν) as a function of two variables η2 and ν for
 5

tisymmetric case while the dark-dark, and bright-bright

solitary wave solutions corresponds to the symmetric
case. It should noted here that the bright-dark mode
has been recorded experimentally in the garfish olfac-
tory nerve, squid giant axon, and hippocampal neuron
[1–4, 21].

FIG. 1. Evolution of the square of the soliton width as a func-

tion of the propagation speed ξ, as predicted by Eq. (13b).
(i) The order parameters are ν = 0.05, η1 = 0.05, η2 = 0.05.
For a2 < 0, that is complex pulse width (which has no phys-
ical significance), Eq. (9) describe periodic solutions for real
amplitudes A and B. While for a2 > 0, that is real pulse
width Eq. (9) describe pulse solitons. (ii) The effect of iner- FIG. 2. The evolution of soliton amplitude with propagation
tia parameter on the dynamics of the soliton. Increasing the speed ξ and initial parameter η2 for different values of the
inertia parameter leads to an increase in the stability of the damping coefficient ν and their corresponding contour plot in
soliton pulse. (ii) The effect of damping parameter on the the η2 − ξ space of Eq. (12): (i)-(iv) ν = 0.05, (ii)-(v) ν =
dynamics of the soliton pulse.(iv) The level-zero contour plot 0.08 and (iii)-(vi) ν = 0.5. The other model parameters are
of Eq. (13a). The other model parameters areare η1 = 0.05, c2
ρA
0 = 4.035×10
−3
, β = 79.5 (ρA0)2 . The amplitude response as
c2
c0 = 176.6, ρA
0 = 4.035 × 10
−3
, β = 79.5 (ρA
0
2 . The level-zero 0
a function of the propagation speed and the initial parameter
0 )
contour depicts all the solution to the equation a20 = F (η2 , ν) shown the threshold and saturation effect. The evolution of
in the contemplated parameter range. Increasing the damp- amplitude as a function of damping parameter ν show the
ing parameter leads to a decrease in the stability of the soliton effect of discontinuity in the propagation speed.
pulse.
Now consider the behavior of the soliton amplitudes
with the parameters of the system. Analysis of the soliton
the bifurcation point ξ = ξ0 . Figure 2(iv) shows the 3D amplitudes can be used to predict nonlinear phenomenon
surface a20 = F (η2 , ν) as introduced in Eq. (13a). This such as symmetry breaking or discontinuity. It should
surface depends on two parameters η2 and ν. The actual be noted that in [31], symmetry breaking in a linearly
colours of any spot (η2 , ν, a20 ) ǫ ℜ3 on the surface depends coupled Korteweg-de Vries system was observed. It was
on the height or level of a20 above the η2 −ν plane. In this shown that the linear KdV equation exhibit symmetry
particular graph, the “height” of a20 ranges from 0−10000 breaking depending on the propagation velocity. Also in
in the surface plot. Ref. [33], the authors investigated numerically the exis-
So far we have analyzed the solutions obtained using tence, and stability of solitons in parity-time-symmetric
the square of the pulse width described by Eq. (13a). If optical media characterized by a generic complex hy-
we considered that the amplitudes A and B are real func- perbolic refractive index distribution and fourth-order
tions and that a2 > 0, then Eq.(9) is said to describe soli- diffraction. They showed that the fourth-order diffrac-
tary wave solution which can be symmetric (i.e., A = B) tion coefficient greatly alters parity-time-breaking points.
or antisymmetric (i.e., A = −B ) [31]. Hence, from Eqs. In lipid membrane, discontinuity or non-linearity in state
(12) and Eqs.(9), the two-component soliton model ad- diagram has been studied both theoretically and experi-
mits three possible types of solutions: dark-dark, bright- mentally [1–5, 34], and it is thought to results from the
bright, and bright-dark solitary wave solutions. The liquid-expanded-liquid condensed phase change in lipid
bright-dark solitary wave solution correspond to the an- monolayers. In fact, it is reported that the velocity of
 6

0.05 0.08 per bound on the maximum amplitude (all-or-none) and

(i) (ii) are self-sustaining. In [35], Shrivastava further proposed
0.04
0.06
that the observed solitary shock waves at lipid interface
0.03
0.04
also provide evidence for the detonation of shock waves
at such interfaces. That is, the self-sustaining shock wave
u,v

u,v
0.02
0.02 utilizes the latent heat of phase transition of the lipids
0.01 (i.e., chemical energy stored in lipids interface) reinforc-
0
0 ing it in the process. One needs to stress here that the
general form of the nonlinear evolution we have consid-
-0.01 -0.02
-200 0 200 -200 0 200 ered allows for other possible sources of chemical energy
0.5
x 0.8
x that can reinforce the propagating shock wave. Further-
(iii) (iv) more, shock waves have been reported to be signature
0.4
0.6 for traumatic brain injury [36]. It is well known that
0.3
exposure of biological cells to shock waves causes dam-
0.4
0.2 age to the cell membrane, however, the mechanisms by
u,v

u,v

0.1
0.2
which damage is caused and how it depends on physical
0 parameters of the shock waves such as shock waves veloc-
-0.1
0 ity, shock waves duration, shock waves shape, and am-
plitudes is poorly understood. Figure 5 shows generation
-0.2
-200 0 200 -0.2 of solitary shock-like waves at particular propagation ve-
-200 0 200
x x locity, with different amplitudes and shock waves profiles.
Sliozberg et al. demonstrated numerically using coarse-
grained model of lipids vesicle the principle of damage
FIG. 3. The evolution of pulse shape as a function of propa- induced by shock waves by direct passage through the
gation speed ξ within one of the one of stable island describe cranium. The results show that the structural integrity
in figure 4(iv): (i) ξ = 165.8, (ii) ξ = 165.5,(iii) ξ = 165.0 and of the lipid vesicles is altered as pores are formed in the
ξ = 164.9. The other model parameters are ǫ = 0.001, t = 0, lipid membrane. As a result, the membrane becomes
c2
η1 = 0.045, ν = 0.05, ρA 0 = 4.035 × 10
−3
and β = 79.5 (ρA0)2 . permeable to sodium, potassium and calcium ions and
0
The solitary wave evolve into a solitary oscillatory shock- therefore a possible source of chemical energy for the
like waves as the propagation speed is gradually reduced. shock waves [36]. It should be stressed that, a number
Waves with smaller propagation speed turns to propagates of scientists have shown that membranes can display ion
with higher amplitudes. channel- like current traces even in the complete absence
of proteins provided only that the membrane is close to a
melting transition [37–39]. These events are called “lipid
channels” and consists of small pores or defects in the
sound decreases discontinuity as the system goes from
lipid membranes. We proposed that the observed soli-
LE phase to LE-LC coexistence phase in liquid monolay- tary shock-like wave in the coupled-model in might be
ers. Figures 3 (iv) and 3 (v) exhibits in parameter space
responsible for the defects in the lipid membrane near
(η2 , ξ) the stable and unstable regions. One observes the phase transitions.
discontinuity of the stable region, which appear like the
Consider the interaction of the two coupled solitary
island of points when the damping parameter is varied.
waves. From Eqs. (8), (9), and (10a), and noting that
It should be noted that discontinuity in speed implies A2 = B 2 , we obtain the effective potential energy func-
a discontinuity in the compressibility. Discontinuity in
tion Us of the coupled solitary waves
biomembranes indicates an abrupt-phase transition, and
γ 2 γπ(c2o − 2ξ) + γ 3 π
 
near phase transitions, stable solitary waves can propa- γπ
Us = o cosec A2
gate along the nerve fiber. In order to investigate stable 2γγo3 a2 2γo a
solitary waves near the region of discontinuity, we plotted βπγ γ 2

4γo2
 
γπ

the soliton profile within the stable region in Fig.3 (iv) for + − cosec A4 . (15)
different values of the propagation speed (see Fig.3). It 72γo4 a4 a2 2γo a
is worth noting that as the speed of the wave is gradually In order to conveniently analyze the characteristic of
reduce, the single pulse created evolved into an oscilla- energy function, we plot the curve of Us against ξ. In
tory solitary shock-like wave. Figure 2 (vi), exhibits con- view of Fig. 4, Us is positive in modeling region and
tinuity in the stable region within the parameter space therefore the potential function of interaction is always
(η2 , ξ) resulting from an increase in the damping coeffi- negative. Thus, the interaction between the two solitary
cient. In Refs. [34, 35], Shrivastava et al. have recently waves is attractive; this agrees with the fact of neural
found that a two-dimensional solitary shock wave rem- coupling oscillation.
iniscent of propagating spikes in nerves can be induced In this section, we analyzed the characteristics of two
in monolayer lipid membrane near a phase transition. coupled solitary waves in nerve impulse transmission. We
These waves have a threshold for excitation and an up- employ an approximate approach using the variational
 7

IV. LINEAR STABILITY ANALYSIS

(a) (b) In non-linear waves, solitons and their linear stabil-

15000 15000 ity properties are crucial. In integrable systems, solitons
admit analytical expressions and they are generally sta-
10000 ble against small perturbations. In non-integrable sys-
10000
tems, solitons generally do not admit analytical expres-
5000 sions and can be either stable or unstable. The under-
5000 standing of many properties of these nonlinear waves can
Us

0 be achieved through the calculation of their spectra. In-

Us
η 2= 0.110 deed, spectral analysis of the equations appearing after
0 ν = 0.55
-5000 η 2= 0.112 linearization of the governing non-linear equations near,
ν = 0.85
η 2= 0.114 ν = 1.50 the solitons has become widespread tools. Bound eigen-
-5000 states, which do not grow with propagation or in time,
50 100 150 50 100 150
ξ ξ of the linear problem, are often called internal modes.
Internal modes supported by the single mode soliton in
FIG. 4. The energy of the two solitary waves Us as a function the soliton model of nerve has been studied in [19]. In
of propagation speed ξ. (a) The effect of the inertia parameter this section, we compute the linear-stability spectrum of
η2 for ν = 0.55. (b) The effect of the damping parameter ν for the solitary wave obtained; which consists of eigenval-
η2 = 0.110. The energy of the two solitary waves decrease as ues of the linear-stability operator of the solitary wave
the damping coefficient and inertia coefficient increases. The [48]. The linear-stability spectrum contains valuable in-
other model parameters are η1 = 0.045, c0 = 176.6, ρA 0 = formation on the solitary wave. For instance, if this spec-
c2
4.035 × 10−3 and β = 79.5 (ρA0)2 . trum contains eigenvalues with positive real parts, then
0
the solitary wave is linearly unstable. In this case, the
largest real part of the eigenvalues gives the maximal
principle and obtained a solution for the coupled solitary growth rate of perturbations. If the spectrum contains
wave within the framework of the soliton model described purely imaginary discrete eigenvalues, these eigenvalues
by coupled KdV-BBM equations in the presence of damp- are the internal modes [48, 49]. In particular, persistent
ing. We showed that, a small friction potential present oscillations of the soliton width and position are due to
in the coupled solitary wave model induced discontinuity these modes.
(or phase transitions) in the amplitude (density change)- To proceed, we consider a small perturbation in u and
velocity parameter space, a reminiscent for propagation v given by u = u0 + ϕu1 and v = v0 + ϕv1 where ϕ << 1
of solitary shock-waves and a single soliton pulse. The and u0 and v0 is the solution to the Eqs. (4) and (5) given
relevant biological implication is also discussed. In the by Eq. (9). Substituting u = u0 + ϕu1 and v = v0 + ϕv1
next section we will study linear stability analysis of the into Eqs. (4) and (5) and linearising about ϕ the linearise
solitary waves solution obtained in this section. system

∂u1 1 2 ∂u1 ∂u0 η1 ∂ 3 u1 ∂ 3 u1 ν ∂ 2 u1

+ c0 + β(u20 + v02 ) + β(u0 u1 + v0 v1 ) − − η2 − = 0, (16a)
∂s 2 ∂y ∂y 2 ∂y 3 ∂y 2 ∂s 2 ∂y 2
∂v1 1 2 ∂v1 ∂v0 η1 ∂ 3 v1 ∂ 3 v1 ν ∂ 2 v1
+ c0 + β(u20 + v02 ) + β(u0 u1 + v0 v1 ) − − η2 − = 0. (16b)
∂s 2 ∂y ∂y 2 ∂y 3 ∂y 2 ∂s 2 ∂y 2

Since the linearisation is near a solitary wave, i.e, Eq. (9), problem:
then Eqs. (16a) and (16b), reduces into an eigenvalue
QV = λV. (17)

In Eq. (17), V(y) = (u1 , v1 )T is the eigenfunction, y =

(y1 , ..., yN ), where N is the number of spatial dimensions.
λ is the eigenvalue while Q is the linearization operator
(usually non-Hermitian) and is given by

(1 − η2 ∂y2 )−1 (G1 ∂y − η1 ∂y3 − ν∂y2 + G0 ) (η2 ∂y2 − 1)−1 G2

 
1
Q= 2 (18)
2 −1
(η2 ∂y − 1) G2 (1 − η2 ∂y ) (G1 ∂y − η1 ∂y3 − ν∂y2 + G0 )
2 −1
 8

2 2
G0 = −2βA sech (ay)tanh(ay),
 |Im(λ)| > 0, then the solution (u, v) will grow exponen-
G1 = c20 + 2βA2 sech2 (ay), (19) tially with y (i.e., it is unstable), otherwise, the solution
G = 2aA2 sech2 (ay) tanh(ay).
 (u, v) is stable. In the existence range of the approxima-
2
tion described above, the dependence of the perturbation
The continuous spectrum of the system is determined growth rate (the most unstable growth rate) on the pa-
analytically by examining the matrix Q as y → ∞. Pro- rameter η1 and η2 are exhibited Fig. 5 (a) and 5(c-d)
ceeding with this approximation, it is easy to show that respectively. As it can clearly be observed, η1 and η2
the continuous eigenvalues are reduces the instability growth rate of the solitons.
In order to confirm the stability or instability of the
1 c20 k + η1 k 3 νk 2 system, one needs to solve the eigenvalue equation (17)
 
λ=± i + . (20) numerically. One of the methods used to solved such a
2 1 + η2 k 2 1 + η2 k 2
problem is the finite difference discretizations method.
where k is the perturbation wave number. Discrete eigen- However, the accuracy of this method is quite low. In
values do exist, however, they can only be obtained nu- addition, this method has been reported to give spuri-
merically. The local growth rate of the solitons or the ous eigenfunction even when the eigenvalues obtain are
gain is given by Im(λ). The maximum growth rate approximately correct. A more accurate method is the
of the soliton can be obtained by solving the equation Fourier collocation method (FCM) [48]. By this method,
dIm(λ)
dk = 0. Applying this definition, it is easy to show the eigenfunction V is expanded into a Fourier series
that, and Eq. (17) is turned into a matrix eigenvalue prob-
p lem for the Fourier coefficients of the eigenfunction V.
2 η2 c20 − 3η1 ± (3η1 − η2 c20 )2 − 4η1 η2 c20 The structure of the resulting matrix depends heavily on
k± = . (21)
2η1 η2 the structure of the spatial dimension of Eqs. (4) and
(5). A detail analysis of the Fourier collection method
can be found in [48–50]. The results from the numerical
×105
simulation is illustrated in Fig. ( 6).
2 Figure 6 (a1, b1, c1) illustrate the unperturbed soliton
(a) η =0.01 (b)
2 profiles for different values of the propagation speed ξ
1.5 η =0.02
2 100 k used in our numerical simulation while Fig. 6 (a2, b2,
Im(λ)

η 2=0.03 -
1 k+ c3) corresponds to the numerically computed eigenvalue
k

50 spectrum. As can clearly seen, Fig. 6 (a2) describe a

0.5 bound state or internal mode. In particular, persistent
oscillations of the soliton width, and position are due to
0 0
0 50 100 6 7 8 these modes. On the other hand, Fig. 6 (a3) illustrate
×105 k ×106
η2 ×10-5 to the stable propagation of the perturbed solution while
5 8
(c) (d)
Fig. 6 (b3, c3) illustrate unstable propagation of the
4 η 1=0.05 η =0.05 intensity of the perturb solution.
max (Im(λ))

max(Im(λ))

6 1
η =0.08 η 1=0.08
3 1
η =0.10 4 η =0.10
1 1
2

1 2

0 0
0 0.5 1 0 0.5 1
-4
η2 ×10 η2 ×10-4

FIG. 5. (a)The perturbation growth rate of solitons versus k

with different values of η2 , η1 = 0.05 (b) Bifurcation diagram
of soliton states in the two components model. 0 < η2 <
5.32 × 10−5 , the wave number k = 0, k = 5.1 × 10−5 , there
is a growth in the wave number k and at k = 5.32 × 10−5 ,
the soliton split into two soliton states with wave number k−
and k+ . Note that their antisymmetric components −k− and
−k+ are not shown on the diagram and η1 = 0.55 (c) and
(d) corresponds to the maximum perturbation growth rate of
solitons states k− and k+ respectively versus η2 with different
values of η1 . The order parameter c0 = 176.6.

The stability of the perturbed soliton u, v is related

to the imaginary parts Im(λ) of all eigenvalues λ. If
 9

ied by a a wide variety of applied scientists; and, partly

as a result of recent interest shown by applied mathe-
maticians in the nonlinear diffusion equation, single pulse
dynamics is now a rather well understood physical phe-
nomenon [13–19]. One might suppose that this is the end
of the story; on the contrary, it is only the beginning. Re-
cent studied have shown that the electromechanical den-
sity pulses in the nerve can effectively be considered as a
vector soliton with two components. That is, the longitu-
dinal component corresponding relative height field and
the transverse components corresponding to the lateral
stretch field and have been observed experimentally in
garfish olfactory nerve, squid giant axon, and hippocam-
pal neuron [21].
In this work, we have considered the electromechanical
density pulse as a two coupled solitary waves represented
by longitudinal compression wave and an out-of-plane
transversal pulse (i.e., perpendicular to the membrane
surface) and analyzed using the variational approach the
characteristics of the coupled solitary waves in the pres-
ence of damping within the framework of coupled nonlin-
ear Burger-Korteweg-de Vries-Benjamin- Bona-Mahony
equations(BKdV-BBM) derived from the vector solito
FIG. 6. (a1, b1, c1) Initial solitary wave profile as given by
model equation for biomembranes and nerves. In par-
Eq. (9). (a2, b2, c2) Numerically computed linear stability ticular, we have shown that, there must be a balanced
spectra. (a3) Stable and (b2, c3) unstable propagations of between damping and inertia effects for a stable coupled
perturbed solution u = u0 + ϕu1 and v = v0 + ϕv1 . (a1, solitary waves to propagates within the axon. Further-
a2, a3) ξ = 165.8, (b1, b2, b3) ξ = 165.5, and (c1, c2, c3) more, we have shown that the presence of damping co-
ξ = 165.0. The other model parameters are ǫ = 0.001, η1 = efficient causes a discontinuity in the (η2 , ξ) parameter
0.045, ν = 0.05, ρA 0 = 4.035 × 10
−3
, ϕ = 10−33 , η2 = 0.11, space. In particular, we observed a discontinuity in sta-
c2
and β = 79.5 (ρA )2 .
0 ble regions stable region which appears like the island of
0
points. Analysis of the solitary wave energy shows that
both η2 and ν causes a decrease the energy of the cou-
pled solitary waves. In addition, the potential function
V. DISCUSSION AND CONCLUSION of interaction is shown to be always negative. Thus, the
interaction between the two solitary waves is attractive;
which is agreement with the fact of neural coupling os-
Electromechanical solitary wave measurements on cillation. We performed linear stability analysis and the
nerve fibers began in the early 1980s with the measuring results shows that both η2 and ν decreases the perturba-
of various non-electrical components of action potential tion growth rate of the coupled solitary waves. Numeri-
and the investigation of the chemistry of phase transi- cal simulation of the linearized equation shows that the
tions in nerve fibers, and its importance for nerve pulse bell shape solitary wave profile yields a bound eigenvalue
propagation by Tasaki [51] and by Kaufmann proposal spectrum while the solitary shock-like solitary waves pro-
of sound waves as a physical basis for action potential duce an unbound eigenvalue spectrum. Thus the bell
propagation in the nerve. This mechanical aspect of ac- shape profile is generally stable for a small perturbation
tion potential gained more attention when Heimburg and while the solitary-shock like profile is generally unstable.
Jackson proposed the soliton model for nerve pulse prop- The instability of these shock-like solitary waves might
agation in the early 2000s [4]. Since that time, the dy- be responsible for traumatic brain injury and damage to
namics of a pulse, propagation has been extensively stud- the cell membrane [36].

[1] J. Griesbauer, S. Bossinger, A. Wixforth, and M. Schnei- [4] T. Heimburg and A. D. Jackson, PNAS 102 9790 (2005).
der, Phy. Rev. Lett. 108, 198103 (2012). [5] S. S. L. Andersen, A. D. Jackson and T. Heimburg, Prog.
[2] J. Griesbauer, M. F. Schneider, Biophy. J. 97, 2710 Neurobiol. 88 104 (2009).
(2009). [6] A. L. Hodgkin and A. F. Huxley, J. Physiol. London 117,
[3] J. Griesbauer, S. Bossinger, A. Wixforth, and M. Schnei- 500(1952).
der, Phy. Rev. E 86, 061909 (2012). [7] A. L. Hodgkin and A. F. Huxley, J. Physiol. 104, 176
 10

(1945). Clark, Phys. Rev. Lett. 86, 564 (2001).

[8] A. L. Hodgkin, A.F. Huxley, and B. Katz, J. Physiol. [30] E. M. Yamakou, E. M. Inack, and F. M. Moukam Kak-
116, 424(1952). meni, Nonlin. Dyn. 83, 541 (2015).
[9] E. Neher and B. Sakmann, Nature 260, 799 (1976). [31] A. Espinosa-Ceron, B. A. Malomed, J. Fujioka, and R.
[10] A. D. Doyle, J. M. Pfuetzner, A. Kuo, S. L. Cohen, B. F. Rodriguez, CHAOS 22, 033145 (2012).
T. Chait, and R. Mackinnon, Science, 280, 69 (1998). [32] J. L. Bona, M. Chen, and J. C. Saut, J. Nonlinear Sci.
[11] G. H. Kim, P. Korstein, L. Obeid, and B. M. Salzberg, 12, 283 (2002).
Biophys. J. 92, 3122(2007). [33] C. G. L. Tiofack, F. Ndezana, A. Mohamadou, and T. C.
[12] S. Schrivastava and M. F. Schneider, J. R. Soc. Interface Kofane, Phys. Rev. E. 97, 032204(2018)
11, 0098 (2014). [34] S. Shrivastava, K. H. kang, and M. F. Schneider, Phys.
[13] F. Contreras, F. Ongay, O. Pavon,and M. Aguero, Int. J. Rev. E 91, 012715(2015)
Mod. Nonlin. Theory Appl. 2, 7 (2013). [35] S. Shrivastava, Proc. Mtgs. Acoust. 34, 045034 (2018)
[14] F. Contreras, H. Cervantes, M. Aguero, and Ma. de L. [36] Y. Sliozberg and T. Chantawansri, J. Chem. Phys. 141,
Najera, J. Nonlin. Dyn. 2014, 710152 (2014) 184904 (2014).
[15] E. V. Vargas , A. Ludu, R. Hustert, Peter Gumrich, A. D. [37] A. Blicher , K. Wodzinska , M. Fidorra , M. Winterhalter,
Jackson, and Thomas Heimburg, Biophys. Chem. 153, and T. Heimburg, Biophys. J. 96, 4581(2009).
159 (2011). [38] K. R. Laub, K. Witschas, A. Blicher, S. B. Madsen, A.
[16] J. Engelbretch, T. Peets, K. Tamm, M. Laasmaa, and Luckhoff , and T. Heimburg, Biochim. Biophys. Acta
M. Vendelin, Proc. Estonian Academy of Sciences 67, 28 1818, 1123 (2012).
(2018) [39] A. Blicher and T. Heimburg, PloS One 8, e65707 (2013).
[17] J. Engelbretch, K. Tamm, and T. Peets, Biomech. Chem. [40] D. Anderson, Phys. Rev. A 27, 3135 (1983).
Phys. Med. NMR 14, 159 (2015). [41] G. Iooss and D. D. Joseph, Elementary stability and Bi-
[18] G. Fongang Achu, F. M. Moukam Kakmeni, and A. M. furcation Theory (Springer, New York, 1980).
Dikande, Phy. Rev. E 97, 012211 (2018). [42] D. Schuch, J. Phys.: Conf. Ser. 380, 20120(2012)
[19] G. Fongang Achu, S. E. Mkam Tchouobiap, F. M. [43] P. Caldirola, Nuovo Cimento 18, 393 (1941)
Moukam Kakmeni, and C. Tchawoua, Phy. Rev. 98, [44] E. kanai, Progr. Theor. Phys. 3, 537 (1948).
022216 (2018). [45] F. G. Mertens, N. R. Quintero, and A. R. Bishop, Phys.
[20] V. Vogel and D. Möbius, Langmuir 1989, 129 (1988). Rev. E 81, 016608 (2010).
[21] A. E. Hady and B. B. Machta, Nat. Commun. 6, 6697 [46] F. G. Mertens, N. R. Quintero, I. V. Barashenkov, and
(2005). A. R. Bishop, Phys. Rev. E 84, 026614 (2011)
[22] W. M. Grill, S.E. Norman, and R. V. Bellamkonda, [47] N. R. Quintero, F. G. Mertens, and A. R. Bishop, Phys.
Annu. Rev. Biomed. Eng. 11,1(2009). Rev. E 91, 012905 (2015)
[23] L. Tonello and M. Cocchi, ‘ NeuroQuantology 8, 1 (2010). [48] J. Yang, Nonlinear Waves in Integrable and Noninte-
[24] J. Engelbrecht, K. Tamm, and T. Peets, Biomech. Model. grable Systems (Society for Industrial and Applied Math-
Mechanobiol. 14, 159 ( 2015). ematics, Philadelphia, USA, 2010)
[25] P. J. Morrison, Rev. Mod. Phys. 70, 467(1998). [49] D. Gottlieb and S.A. Orszag, Numerical Analysis of Spec-
[26] R. S. MacKay and P. G. Saffman, Proc. R. Soc. London, tral Methods: Theory and Applications (Society for In-
Ser. A 406, 115 (1986). dustrial and Applied Mathematics, Philadelphia, USA,
[27] A. Hasegawa, Plasma Instabilities and Nonlinear Effects 1977 ).
(Springer-Verlag, Berlin, 1975) [50] M. Grillakis, Comm. Pure Appl. Math. 41, 747(1988)
[28] A. J. Brizard, J. J. Morehead, and A. N. Kaufman, Phys. [51] I. Tasaki, Physiol. Chem. Phys. Med. NMR 20,
Rev. Lett. 77, 1500 (1996). 251(1988).
[29] D. L. Feder, A. A. Svidzinsky, A. L. Fetter, and C. W.