Extended Heimburg–Jackson Model

• The extended Heimburg–Jackson model is a nonlinear dispersive wave framework that represents nerve pulses as mechanical density waves with quadratic and cubic nonlinearities.
• It rigorously characterizes solitary-wave existence, polarity, and stability through explicit parameter regimes and an algebraic test based on energy moment derivatives.
• The model is extended to include coupling with myelin sheath dynamics, higher-order nonlinearities, and weak dissipation, yielding novel soliton and Lambert W-kink solutions.

The extended Heimburg–Jackson model is a class of nonlinear dispersive continuum models for nerve-pulse propagation in biomembranes and axons, derived from the Heimburg–Jackson proposal that nerve impulses can be represented as propagating mechanical density waves rather than solely as Hodgkin–Huxley electrical spikes. In the analytically most developed formulation, the model is the generalized Boussinesq equation

$v_{tt}+\bigl(-v+a v^2+b v^3\bigr)_{xx}+v_{xxxx}=0,$

with quadratic and cubic nonlinearities controlled by $a$ and $b$. Subsequent extensions retain the membrane-wave perspective while adding dissipation, mixed inertial–dispersive terms, coupling to electrical and pressure variables, a myelin-sheath internal field, or higher-order constitutive nonlinearities that generate generalized Duffing/Liénard reductions and Lambert $W$-kink solutions (Freistühler et al., 2013).

1. Canonical equation and traveling-wave reduction

The extended Boussinesq form studied in the stability analysis is

$$v_{tt} + (p(v))_{xx} + v_{xxxx}=0,\qquad p(v)=-v+a v^2+b v^3.$$

Within this formulation, $a$ controls the quadratic nonlinearity and $b$ controls the cubic nonlinearity. Their sign and relative magnitude determine the existence of positive or negative solitary waves, linear well-posedness at constant states, and the stability or instability of solitary waves. A structurally important point is that the mixed case $ab\neq 0$ cannot be reduced by a simple scaling, unlike the pure quadratic or pure cubic limits (Freistühler et al., 2013).

Solitary waves are introduced through

$v(x,t)=V(x-ct),\qquad V(\pm\infty)=0,$

which yields the profile equation

$V''=(1-c^2)V-aV^2-bV^3.$

The associated potential is

$a$0

so that the profile ODE is the Newton equation $a$1, with first integral

$a$2

A solitary wave must approach $a$3 as $a$4, so $a$5 must be a saddle point of the phase portrait. This requires

$a$6

This canonical reduction fixes the central mathematical structure of the model: a dispersive nonlinear wave equation whose traveling waves are homoclinic orbits of an effective Hamiltonian system. In the nerve-pulse interpretation, the field $a$7 represents a longitudinal density perturbation propagating in the membrane.

2. Solitary-wave existence and polarity

The existence theory for the generalized Boussinesq equation is complete in the sense that the admissible combinations of sign, speed, and nonlinearity are explicitly characterized. Positive solitary waves of speed $a$8 exist if and only if either $a$9, $b$0, and $b$1, or $b$2, $b$3, and

$b$4

Negative solitary waves exist if and only if either $b$5, $b$6, and $b$7, or $b$8, $b$9, and

$W$0

Thus the sign of $W$1 selects the admissible polarity when $W$2, while $W$3 permits both positive and negative solitary waves across the full subsonic interval (Freistühler et al., 2013).

The extrema of the solitary-wave profile are explicit: $W$4

$W$5

A positive solitary wave has maximum $W$6, and a negative solitary wave has minimum $W$7.

These formulas show that the extended Heimburg–Jackson model does not merely support a single pulse family. It supports distinct positive and negative branches, with admissibility constrained by the nonlinear constitutive coefficients. This is significant because later extensions exploit precisely this sensitivity of waveform class to constitutive structure, either by adding new fields or by adding higher-order nonlinear terms.

3. Stability theory, Heimburg–Jackson pulses, and instability regimes

For stability analysis, the equation is rewritten as the first-order system

$W$8

A traveling wave $W$9 is orbitally stable in the standard $$v_{tt} + (p(v))_{xx} + v_{xxxx}=0,\qquad p(v)=-v+a v^2+b v^3.$$0 sense if initial closeness implies global existence and closeness modulo spatial translation. The central variational quantity is the moment of instability

$$v_{tt} + (p(v))_{xx} + v_{xxxx}=0,\qquad p(v)=-v+a v^2+b v^3.$$1

The Grillakis–Shatah–Strauss / Bona–Sachs criterion used in the analysis states that a solitary wave is stable if and only if $$v_{tt} + (p(v))_{xx} + v_{xxxx}=0,\qquad p(v)=-v+a v^2+b v^3.$$2, and unstable if $$v_{tt} + (p(v))_{xx} + v_{xxxx}=0,\qquad p(v)=-v+a v^2+b v^3.$$3 (Freistühler et al., 2013).

The paper defines Heimburg–Jackson pulses as solitary waves satisfying

$$v_{tt} + (p(v))_{xx} + v_{xxxx}=0,\qquad p(v)=-v+a v^2+b v^3.$$4

This is the parameter regime singled out in the original Heimburg–Jackson biological proposal. The main theorem is that all Heimburg–Jackson pulses are stable. For positive waves,

$$v_{tt} + (p(v))_{xx} + v_{xxxx}=0,\qquad p(v)=-v+a v^2+b v^3.$$5

and the proof reduces the sign of $$v_{tt} + (p(v))_{xx} + v_{xxxx}=0,\qquad p(v)=-v+a v^2+b v^3.$$6 to the positivity of

$$v_{tt} + (p(v))_{xx} + v_{xxxx}=0,\qquad p(v)=-v+a v^2+b v^3.$$7

In the relevant regime,

$$v_{tt} + (p(v))_{xx} + v_{xxxx}=0,\qquad p(v)=-v+a v^2+b v^3.$$8

which forces $$v_{tt} + (p(v))_{xx} + v_{xxxx}=0,\qquad p(v)=-v+a v^2+b v^3.$$9 on the integration interval and hence $a$0.

A notable structural result is that linear well-posedness at constant states is governed by exactly the same inequality. Linearizing about a constant state $a$1 gives

$a$2

and Fourier analysis yields linear well-posedness if and only if $a$3 for all $a$4, which is equivalent to

$a$5

Thus, in this family, the equation is linearly well-posed at all constant states if and only if the Heimburg–Jackson condition holds, and in the same regime all solitary waves are stable.

Outside the Heimburg–Jackson regime, the stability picture is mixed rather than uniform. If $a$6 and

$a$7

there exist thresholds $a$8 such that positive waves with $a$9 are stable and positive waves with $b$0 are unstable. There also exist $b$1 such that negative waves with $b$2 are unstable and negative waves with $b$3 are unstable. The same conclusion holds with positive and negative waves interchanged if $b$4. The paper also gives an explicit closed-form discriminant $b$5 whose sign matches the sign of $b$6, so $b$7 implies stability and $b$8 implies instability. This establishes an exact algebraic stability test and yields a fast/slow wave transition resembling the transition known from FitzHugh–Nagumo theory (Freistühler et al., 2013).

4. Improved Heimburg–Jackson dynamics and myelin-coupled microstructure

A distinct extension studies the deformation of the myelinated axon wall by retaining the improved Heimburg–Jackson membrane equation and coupling it to an additional internal field for the myelin sheath. The improved Heimburg–Jackson equation is

$b$9

where $ab\neq 0$0 is the longitudinal density change in the biomembrane, $ab\neq 0$1 and $ab\neq 0$2 are dispersive coefficients, $ab\neq 0$3 is a dissipation coefficient, and $ab\neq 0$4 couples the mechanical wave to the action potential $ab\neq 0$5, ionic currents $ab\neq 0$6, and axoplasmic pressure $ab\neq 0$7 (Tamm et al., 2021).

To include myelin, the model is extended to

$ab\neq 0$8

Here $ab\neq 0$9 is an internal field describing the mechanical influence of the myelin sheath, $v(x,t)=V(x-ct),\qquad V(\pm\infty)=0,$0 and $v(x,t)=V(x-ct),\qquad V(\pm\infty)=0,$1 are characteristic wave speeds of the membrane and myelin components, $v(x,t)=V(x-ct),\qquad V(\pm\infty)=0,$2 is the characteristic frequency scale of the myelin internal mode, and $v(x,t)=V(x-ct),\qquad V(\pm\infty)=0,$3 are coupling coefficients. The authors stress that $v(x,t)=V(x-ct),\qquad V(\pm\infty)=0,$4 is best viewed as an internal variable or extra degree of freedom, not merely an auxiliary correction.

In the linearized model, eliminating $v(x,t)=V(x-ct),\qquad V(\pm\infty)=0,$5 yields an equivalent higher-order PDE for $v(x,t)=V(x-ct),\qquad V(\pm\infty)=0,$6. The elimination makes two consequences explicit. First, the effective low-frequency wave speed is renormalized as

$v(x,t)=V(x-ct),\qquad V(\pm\infty)=0,$7

Second, the reduction generates additional higher-order terms; the paper remarks that a slaving-principle elimination would produce a single higher-order equation with $v(x,t)=V(x-ct),\qquad V(\pm\infty)=0,$8th-order dispersive terms. The characteristic short-wave speed is

$v(x,t)=V(x-ct),\qquad V(\pm\infty)=0,$9

The linear dispersion relation of the coupled system is

$V''=(1-c^2)V-aV^2-bV^3.$0

It supports two distinct branches: an acoustic branch and an optical branch. In the long-wave limit, the reduced acoustic speed is

$V''=(1-c^2)V-aV^2-bV^3.$1

while

$V''=(1-c^2)V-aV^2-bV^3.$2

The optical branch has zero group velocity and infinite phase velocity as $V''=(1-c^2)V-aV^2-bV^3.$3. Larger $V''=(1-c^2)V-aV^2-bV^3.$4 weakens the influence of $V''=(1-c^2)V-aV^2-bV^3.$5, so the correction to $V''=(1-c^2)V-aV^2-bV^3.$6 and the higher-order effects vanish in the limit of large $V''=(1-c^2)V-aV^2-bV^3.$7.

Dissipation enters through $V''=(1-c^2)V-aV^2-bV^3.$8 and, in the effective one-equation representation, generates additional imaginary terms in the dispersion relation: $V''=(1-c^2)V-aV^2-bV^3.$9 Numerically, the paper uses a pseudospectral method on a periodic domain of spatial period $a$00 with $a$01 grid points and initial data

$a$02

with $a$03, $a$04, and $a$05. This initial condition splits into two equal counterpropagating pulses of amplitude $a$06. For reference parameters $a$07, $a$08, $a$09, $a$10, $a$11, with no dissipation and no external forcing, the pulse narrows as it travels, increases in amplitude, and eventually breaks into an oscillatory wave packet. Increasing $a$12 or $a$13 slightly increases speed and amplitude, whereas increasing $a$14 or $a$15 reduces both. Under the tested parameters, inclusion of the myelin sheath tends to slow the mechanical wave (Tamm et al., 2021).

5. Higher-order nonlinearities and Lambert $a$16-kink solitons

A further extension augments the membrane constitutive law by third- and fourth-order polynomial terms. The governing equation for the longitudinal density change $a$17 is

$a$18

The parameters $a$19 are interpreted respectively as the sound speed in the fluid phase, empirical nonlinear elastic coefficients, an elastic/dispersion coefficient, an inertial coefficient from lipid-molecule inertia, and viscous damping from the surrounding axoplasmic fluid (Mendoza-Millán et al., 23 Jul 2025).

After nondimensionalization and the traveling-wave ansatz $a$20, $a$21, the equation reduces, after two integrations and with $a$22, to

$a$23

This is written as a damped nonlinear oscillator

$a$24

described as a Liénard-type equation, reducing to a Duffing oscillator in the special case $a$25.

The analytical tool is the factorization ansatz

$a$26

with compatibility conditions

$a$27

For $a$28, the reduced cubic equation is treated as a FitzHugh–Nagumo-type equation and yields classical kink/antikink profiles. When the quartic and quintic terms are retained, the factorization produces a new explicit solution class,

$a$29

where $a$30 is the Lambert $a$31 function defined by $a$32. The auxiliary coefficient $a$33 must satisfy a cubic equation, and the paper states that $a$34 should be real and positive for physical admissibility.

The paper identifies the biomembrane regime

$a$35

for both the classical kink section and the Lambert $a$36 section, now allowing $a$37 and $a$38 to contribute in the latter. It also gives the effective-potential relation for $a$39,

$a$40

and interprets solitary-wave existence in terms of a local minimum at zero together with at least one adjacent local maximum. The paper connects classical kink/antikink solutions to step-like transitions between membrane states and interprets Lambert $a$41-kinks as asymmetric wavefronts that may model nonlinear recovery or hysteresis. It further argues that different Lambert $a$42 branches may correspond to bistability or threshold excitation, and it advances the broader claim that the extended Heimburg–Jackson model may admit supersymmetric soliton pairs with identical wavefront speed but different governing equations (Mendoza-Millán et al., 23 Jul 2025).

6. Weakly dissipative envelope dynamics, axoplasmic fluid, and impulse collision

Another extension begins from the Heimburg–Jackson density-wave equation with nonlinear sound speed

$a$43

so that

$a$44

To model dissipative coupling to the axoplasmic fluid, the equation is augmented by

$a$45

yielding a weakly dissipative nonlinear wave equation after nondimensionalization (Pavón-Torres et al., 2024).

For low-amplitude nonlinear excitations, the coefficients are scaled with a small parameter $a$46, and the long-wave stretched variables

$a$47

lead to a Burgers–KdV-type evolution equation involving nonlinearity, dispersion, and weak dissipation. A multiple-scale expansion with leading harmonic form then produces the damped nonlinear Schrödinger equation

$a$48

with

$a$49

In the absence of damping, this reduces to the standard focusing NLSE admitting bright solitons.

The perturbation analysis uses the bright-soliton ansatz

$a$50

where $a$51 represent amplitude, center, phase, and velocity/frequency. Karpman–Solov’ev–Maslov perturbation theory yields adiabatic evolution equations for these parameters. For two colliding pulses,

$a$52

with each $a$53 a bright soliton, the interaction terms include

$a$54

Using symmetric variables $a$55 and the interaction regime $a$56, $a$57, $a$58, the core adiabatic equations become

$a$59

$a$60

$a$61

The collision analysis distinguishes two regimes. With no axoplasmic damping, $a$62, the amplitude is constant and the solutions show periodic exchange-like behavior. With gain or loss, $a$63, and the interaction changes accordingly. In the symmetric case, $a$64 is interpreted as increasing the oscillation period and causing repulsion-like separation, while $a$65 decreases the period and causes attraction-like compression. In the antisymmetric case, the solution gives monotonic separation or repulsion, stronger under absorption and slower under gain.

This framework is explicitly contrasted with Hodgkin–Huxley annihilation. The paper states that, unlike Hodgkin–Huxley predictions that colliding action potentials annihilate, the Heimburg–Jackson / soliton picture permits orthodromic and antidromic impulses to penetrate each other and preserve their soliton character, consistent with the invertebrate collision experiments of Gonzalez-Perez et al. (2014). In that interpretation, the axoplasmic fluid does not merely damp the pulse; it modifies soliton parameters through weakly dissipative adiabatic evolution (Pavón-Torres et al., 2024).