On numerical modeling of dispersive mechanical waves in lipid bi-layers

TL;DR

This paper numerically investigates how mechanical waves are generated and influenced by electrical nerve signals in lipid bi-layers, revealing key relationships between wave velocities and shapes in nerve pulse propagation.

Contribution

It introduces a coupled mathematical model to analyze the generation and characteristics of mechanical waves during nerve pulse propagation, highlighting their dependence on system parameters.

Findings

Mechanical wave velocity correlates with action potential velocity.

Wave shape depends on sound velocity in lipid bi-layers.

Mechanical waves influence the velocity and shape of action potentials.

Abstract

We investigate different mechanical effects which accompany the nerve pulse propagation by using mathematical modeling. The propagation process is composed by three connected phenomena: (i) the action potential (electrical signal) which is usually considered when nerve pulses are discussed, (ii) the mechanical wave propagating in the biomembrane and (iii) the pressure wave in the axoplasm inside the axon. The main goal of the present study is to investigate numerically how the mechanical wave is generated by the action potential and how the characteristics of the system are reflected in the emerging wave ensemble. %the existence and characteristics of the mechanical wave influence the parameters of an action potential part of the wave ensemble. The key characteristics for the coupled model system are: (i) the velocity of the peak of the mechanical pulse is associated with the velocity of the action potential regardless of the sound velocity value in the lipid bi-layer, (ii) the velocity of the front and the shape of the mechanical wave depends on sound velocity in the lipid bi-layer, (iii) the shape of the mechanical wave can have an effect on the velocity and shape of the action potential.

Equations

$$\begin{array}{rl}& Z_T=Z\left(Z-C_1-Z^2+C_{1}Z\right)-I+D{Z}_{XX},\\ & I_T=\epsilon \left(Z{C}_2-I\right),\end{array}$$

$$\begin{array}{rl}{U}_{TT}& =c^2{U}_{XX}+NU{U}_{XX}+M{U}^2{U}_{XX}+N{\left(U_X\right)}^2+\\ & +2MU{\left(U_X\right)}^2-H_1{U}_{XXXX}+H_2{U}_{XXTT}+g_3{I}_X,\end{array}$$

$$P_{TT}=c_f^2{P}_{XX}-\mu P_T+F(Z,J,U),$$

$$Z(X,0)=Z_o{\mathrm{sech}}^2{B}_{o}X;Z(X,T)=Z(X+2Km\pi ,T);m=1,2,\dots ,$$

$$\hat{Z}(k,T)=\mathrm{F}\left[Z\right]=\sum _{j=0}^{n-1}Z(j\Delta X,T)\exp \left(-\frac{2\pi \mathrm{i}jk}{n}\right),$$

$$\frac{{\partial }^{m}Z}{\partial X^m}={\mathrm{F}}^{-1}\left[(\mathrm{i}k{)}^m\mathrm{F}(Z)\right],$$

$$\begin{array}{rl}{U}_{TT}-H_2{U}_{XXTT}& =c^2{U}_{XX}+NU{U}_{XX}+M{U}^2{U}_{XX}+N{\left(U_X\right)}^2+\\ & +2MU{\left(U_X\right)}^2-H_1{U}_{XXXX}+g_3{I}_X\end{array}$$

$$U={\mathrm{F}}^{-1}\left[\frac{\mathrm{F}(\Phi )}{1+H_2{k}^2}\right],\frac{{\partial }^{m}U}{\partial X^m}={\mathrm{F}}^{-1}\left[\frac{(\mathrm{i}k{)}^m\mathrm{F}(\Phi )}{1+H_2{k}^2}\right].$$

$$\begin{array}{rl}{\Phi }_{TT}& =c^2{U}_{XX}+NU{U}_{XX}+M{U}^2{U}_{XX}+N{\left(U_X\right)}^2+2MU{\left(U_X\right)}^2-\\ & -H_1{U}_{XXXX}+g_3{I}_X.\end{array}$$

$$\omega =\sqrt{\frac{c^2{k}^2+H_1{k}^4}{1+H_2{k}^2}},$$

Full text

††journal: arXiv

1 Introduction

We investigate different effects associated with the nerve pulse propagation using mathematical modeling. The propagation process is composed of three connected phenomena: (i) the action potential (electrical signal) [15, 24] which is usually considered when nerve pulses are discussed, (ii) the mechanical wave propagating on the axon surface [18] and (iii) the pressure wave in the axoplasm inside the axon [32]. In this analysis the possible thermodynamical effects [5, 23] are left aside. Possibilities for accounting thermodynamic effects within the presented framework can be found in [30].

While the action potential of the nerve pulse propagation has been thoroughly studied, the other two effects are not so well understood. For the action potential the Hodgkin–Huxley model [15] and its modifications are mostly used but there also exists a number of simplified models which are still capable of capturing the key characteristics like, for example, the FitzHugh–Nagumo model (FHN) [24] or the models based on the telegraph equation [6, 23]. In the present paper we use the FHN model for modeling the action potential although it must be stressed that in principle any of the existing action potential models can be used. Here the attention is on coupling effects regarding accompanying phenomena. The main criterion for the used model in the present framework is getting a correctly shaped action potential out of it.

For the mechanical wave we use the improved Heimburg-Jackson model [7, 8, 13, 25] which is, in essence, a Boussinesq-type model with nonlinear terms which are of the displacement-type and two dispersive terms (where the second dispersive term is the improvement to the original model [7]). In the solid mechanics the Boussinesq-type equations are used for investigating nonlinear microstructured materials [4]. The influence of internal structure of solids on wave propagation in such media clearly manifests itself in the wave dispersion. The modeling of such effects can be traced to the Born–von Karman lattice model [1]. Following such an approach, the basic idea is the continualization of discrete systems which brings the higher-order terms in the governing equations into account. An excellent overview on corresponding mathematical models is given by Maugin [22].

For the pressure wave it makes sense to use the classical wave equation as a starting point. Pressure gradients observed experimentally are small [32] allowing one to neglect nonlinear effects. However, axon dimensions are typically in range where viscous effects can be significant enough so the classical wave equation should be improved by taking viscosity into account. Like with other constituents of the present framework the particular model used for describing the phenomenon is not as important as getting a signal shape which is qualitatively similar to experimental observations. Obviously the individual models can be as detailed as needed and what is proposed here should be taken as a starting point and a proof of concept.

In a nutshell our idea is simple. Take established single models describing all effects related to the nerve pulse propagation that should be considered, propose coupling terms between processes that influence each other, and derive a coupled mathematical model describing an ensemble of waves related to the nerve pulse propagation. As noted, doing all this would be only a first step on a long road ahead as once the single models are collected into one framework, much work remains to be done for the experimental verification and for the proper calibration of the models in such a framework.

2 Model equations and statement of the problem

For modeling the action potential we start with the FHN model [24]

[TABLE]

where $Z$ is trans-membrane potential, $I$ is the combined ion (recovery) current, $D$ is a coefficient, $C_{i}=a_{i}+b_{i}$ where $a_{i}$ are the electrical activation coefficients, $b_{i}$ are the ‘mechanical’ activation coefficients and $X,T$ are the dimensionless spatial and time coordinates, respectively. Subscripts $X,T$ here and further denote partial derivatives with respect to the indicated coordinate. The coefficients $b_{i}$ could be taken as $b_{1}=-g_{1}U$ and $b_{2}=-g_{2}U$ where $g_{i}$ are constants and $U$ is from Eq. (2) – the feedback from the mechanical wave. Briefly, the underlying idea behind the ‘mechanical’ activation coefficients is that if the density of the lipid bi-layer (biomembrane) increases then it could make it harder for the ions to propagate through the membrane and if it decreases then ions could more easily penetrate through the lipid bi-layer.

The mechanical wave is modeled by the improved Heimburg - Jackson model [7, 8, 13]

[TABLE]

where $U=\Delta\rho$ is the dimensionless longitudinal density change in the lipid bi-layer, $N,M$ are the nonlinear coefficients, $H_{1},H_{2}$ are the dispersion coefficients, $c$ is the velocity in the unperturbed state, $g_{3}$ is the coupling coefficient and $I_{X}$ is the gradient of $I$ from Eq. (1). It should be noted that while recovery current $I$ is a positive pulse-type quantity, its gradient $I_{X}$ is a bi-polar quantity which is important for the stability of the solution of the Eq. (2).

The coupling term $g_{3}I_{X}$ in Eq. (2) can be further elaborated by taking into account the pressure and displacement effects from the ion current(s). As noted, an improved wave equation equation could be sufficient as a starting point for describing the pressure. As a rough approximation one can assume that the local transverse displacement is proportional to the local pressure change because the transverse displacement observed in the experiments is small [18] compared to the typical diameter of the axon. Here the lipid bi-layer is assumed to be ideally elastic in the presence of such small displacements. Then we could use

[TABLE]

where $P$ is pressure, $c_{f}$ is the sound velocity in axoplasm, $\mu$ is viscosity coefficient. The $F$ is the coupling term accounting for the possible influence from the action potential (ion currents) and mechanical wave in biomembrane. Should the improved wave equation be insufficient it is an option to use some modification, like, for example, the 2D formulation of wave propagation in elastic tube [21]. The existence and characteristics of a pressure wave accompanying the action potential has been measured by Terakawa [32]. An illustrative scheme and a block diagram for the proposed combined model are shown in Fig. 1.

However, in the present paper a two component model containing the action potential (Eq. (1)) and the mechanical wave (Eq. (2)) is used. Simulations including a pressure wave (Eq. (3)) in the numerical scheme can be found in [9].

The main goal of the present study is to investigate numerically how the mechanical wave is generated by the action potential and how the characteristics of the system are reflected in the emerging wave ensemble.

3 Initial and boundary conditions, numerical method

A $\operatorname{sech}^{2}$-type localized initial condition with an initial amplitude $Z_{o}$ is applied to Eq. (1) and we make use of the periodic boundary conditions

[TABLE]

where $K=160$, i.e., the total length of the spatial period is $320\pi$. For Eq. (2) we take initial excitation to be zero and make use of the periodic boundary conditions. The solution representing the mechanical wave described by Eq. (2) is generated over time as a result of coupling between the Eqs (1) and (2).

For numerical solving of the Eqs (1) and (2) we use the discrete Fourier transform (DFT) based pseudospectral method (PSM) [11, 29]. Variable $Z$ is represented in the Fourier space as

[TABLE]

where $n$ is the number of space-grid points ($n=2^{12}$ in the present paper), $\Delta X=2\pi/n$ is the space step, $k=0,\pm 1,\pm 2,\ldots,\pm(n/2-1),-n/2$; $\mathrm{i}$ is the imaginary unit, $\mathrm{F}$ denotes the DFT and $\mathrm{F}^{-1}$ denotes the inverse DFT. Basically, the idea of the PSM is to approximate space derivatives by making use of the DFT

[TABLE]

reducing therefore the partial differential equation (PDE) to an ordinary differential equation (ODE) and then to use standard ODE solvers for integration with respect to time.

The regular PSM algorithm is derived for $u_{t}=\Phi(u,u_{x},u_{2x},\ldots,u_{mx})$ type equations. In our case, however, we have a mixed partial derivative term $H_{2}U_{XXTT}$ in the Eq. (2) and thus the standard PSM has to be modified [16, 17, 29]. Therefore we rewrite the Eq. (2) so that all partial derivatives with respect to time are at the left-hand side of the equation

[TABLE]

and introduce a new variable $\Phi=U-H_{2}U_{XX}.$ After that, making use of properties of the DFT, one can express the variable $U$ and its spatial derivatives in terms of the new variable $\Phi$:

[TABLE]

Finally, Eq. (2) can be rewritten in terms of the variable $\Phi$ as

[TABLE]

In Eq. (9) all partial derivatives of $U$ with respect to $X$ are calculated in terms of $\Phi$ by using expression (8) and therefore one can apply the PSM for numerical integration of Eq. (9). The FHN model (1) is already written in the form of two first-order PDE’s which can be solved by the standard PSM without any further modifications.

The calculations are carried out with the Python package SciPy [19], using the FFTW library [12] for the DFT and the F2PY [27] generated Python interface to the ODEPACK FORTRAN code [14] for the ODE solver. We used ‘vode’ with options nsteps$=10^{7}$, rtol$=10^{-10}$, atol$=10^{-12}$ and $\Delta T=0.1$.

4 Model parameters

We solve Eqs (1) and (2) under localized initial conditions and periodic boundary conditions (4). The parameters for the model equations are $D=1$, $\varepsilon=0.01$ or $\varepsilon=0.05$, $a_{1}=0.2$, $a_{2}=0.2$, $N=0.05$, $M=0.02$, $H_{1}=0.5$, $H_{2}=0.75$, $g_{1}=0.05$, $g_{2}=0.05$, $g_{3}=0.02$. Note that as the model Eqs (1) and (2) are dimensionless then also the model parameters are here dimensionless. The parameter relations to the quantities with dimensions can be found, for example, in [7]. As noted, the length of the spatial period is $320\pi$ and the localized initial condition is given in the middle of the spatial period in the form of sech2-type pulse. The initial condition parameters are $Z_{0}=2$ and the pulse width parameter $B_{0}=1$. It should be noted that if the time-scale parameter is taken as $\varepsilon=0.01$ then such an initial condition is above the threshold leading to emergence of the propagating AP. However, if $\varepsilon=0.05$ then the initial pulse is not sufficient for generating the propagating AP. Such an initial condition is, in essence, a very narrow high amplitude ‘spark’ for the FHN model (1). If the ‘spark’ is above the threshold then the normal FHN pulse is generated in the middle of the spatial period which starts to propagate from there to the positive and negative coordinate $X$ directions. Due to the coupling between the FHN model (1) and the improved Heimburg-Jackson model (2) the mechanical wave is generated in the biomembrane. The initial condition for the Eq. (2) is that the system is at rest (zero initial conditions). If the initial condition for the Eq. (1) is below the threshold then the initial electrical signal vanishes over a short time but a small amplitude mechanical wave is still generated as a function of recovery current gradient until the initial excitation vanishes.

It should be noted that the improved Heimburg-Jackson equation is conservative, i.e., it does not generate nor lose energy, however, due to the added coupling with the FHN equation (1) with the improved Heimburg-Jackson equation it is no longer conservative as it can change energy with the Eq. (1).

The improved Heimburg-Jackson model has two dispersive terms controlled by dispersion parameters $H_{1}$ and $H_{2}$. The dispersion relation for Eq. (2) can be written as

[TABLE]

where $k$ is the wave number and $\omega$ is the dimensionless frequency. The ratio of dispersion parameters $H_{1}^{2}/H_{2}^{2}$ determines the bounding velocity in the system while parameter $H_{2}$ determines the rate at which the bounding velocity is achieved as the frequency of the signal increases. If we have the normal dispersion ($c^{2}>H_{1}^{2}/H_{2}^{2}$) then this means that group speed $c_{gr}$ (which usually is the speed of energy propagation in the system) is smaller than the phase speed $c_{ph}$ and if we have an anomalous dispersion ($c^{2}<H_{1}^{2}/H_{2}^{2}$) then it is the opposite, i.e., $c_{gr}>c_{ph}$ (for more details see [2, 3, 7, 8, 25]). In essence the second dispersive term $U_{XXTT}$ removes the ability of sufficiently high frequencies to travel at near infinite velocities (i.e., makes the model causal). Since the mixed derivative term $U_{XXTT}$ reflects the inertial properties then the parameter $H_{2}$ is related to the inertial properties of the biomembrane [25]. Based on experimental considerations [13] the dispersion type for the lipid bi-layer should be anomalous which is the case under the considered model parameters. The phase velocity curves for the Eq. (2) under the considered parameters are shown in Fig. 2.

5 Results and discussion

The solutions for the coupled model equations (1) and (2) at $T=400$ are shown in Figs 3 and 4. As noted previously, the initial ‘spark’ either leads to a formation of the propagating FHN action potentials which then leads to formation of mechanical waves through coupling (Fig. 3) or if the initial value of the ‘spark’ is below the threshold for the FHN model (1) then it vanishes over a short time leading to a formation of only a small amplitude mechanical wave which proceeds to propagate without the FHN action potential (Fig. 4). Note the difference in the scale of Figs 3 and 4. Here only the waves propagating to the left are shown. There exist similar wave-profiles propagating to the right which are not shown.

In Fig. 3 we compare action potentials $Z$ and mechanical waves $U$ at dimensionless time $T=400$ provided the velocity of the low frequency sound is different for the lipid bi-layer. A number of essential characteristics for the coupled model system can be observed: (i) the velocity of the peak of the mechanical pulse is similar to the velocity of the action potential regardless of the sound velocity value in the lipid bi-layer within the considered parameter range, (ii) the velocity of the front and the shape of the mechanical wave depends on sound velocity in the lipid bi-layer for the Eq. (2), (iii) the shape of the mechanical wave can have an effect on the velocity and shape of the action potential (if $c^{2}=0.25$ for the mechanical wave then the corresponding action potential $Z$ has propagated further compared with the case $c^{2}=0.125$). The mechanisms for changing the sound velocity in the lipid bi-layer and different scenarios depending on the mechanical disturbance velocity when compared against the velocity of an electrical pulse have also been studied by El Hady and Machta using different theoretical considerations [5].

In Fig. 4 one can see the wave-profiles of the mechanical wave at $T=400$ when the initial excitation for the Eq. (1) is below the threshold at three different values for the $c^{2}$. It can be noted that while the front of the wave-profile propagates at the same velocities as in Fig. 3, the peak of the packet travels faster than the corresponding wave-profile in Fig. 3. While not shown here it should be noted that the system can support also solitonic solutions, for example, should we take $c^{2}=0.755$ under the considered parameters (when $\varepsilon=0.05$) we would get a solution where the initial pulse is separated into a train of solitonic pulses while in the shown cases we have an oscillatory structure instead of the solitonic pulses. The solitonic solutions for the Eq. (2) have been studied previously in [8, 25, 26] and for the original model with the simplified dispersion in [8, 13, 20, 33].

The numerical results demonstrate a number of characteristics for the combined nerve pulse model which are qualitatively in line with the observations from the experiments. The moving action potential generates a mechanical wave similar to certain extent to results of El Hady and Machta in [5] although through a different mechanism. It must be stressed that presented profiles are for the longitudinal density change while in experiments usually an transverse displacement is measured. However, the connection between longitudinal density change and transverse displacement can be established in a straightforward manner by drawing inspiration from theory of rods [26, 28]. The transverse displacement in the theory of rods is proportional to the gradient of the longitudinal density change. In Fig. 5 a wave-profile of the longitudinal density change under the considered parameters at $c^{2}=0.25$ (left) and its gradient (right) are shown. The shape of the gradient is similar to the transverse displacement profiles measured in the experiments [18, 31] and as expected, its amplitude is rather small. The shape of transverse displacement derived from the generated mechanical wave (see Fig. 5) is in a good agreement with various experimental observations [18, 31].

The numerical results demonstrate that under the proposed model the mechanical wave can influence the characteristics of the electrical signal, like, for example, its shape or propagation velocity. However, under the considered parameters the influence of the mechanical wave is relatively small on the electrical pulse. Then there is a question of velocity synchronization between the different signals in the ensemble. For example, Heimburg and Jackson have measured the unperturbed (low frequency) sound velocity in the lipid layer to be $171.4\;\mathrm{m/s}$ at $37\mathrm{C^{o}}$ and $176.6\;\mathrm{m/s}$ at $45\mathrm{C^{o}}$ and the minimum possible velocity for their solitonic solution is assessed to be approximately $115\;\mathrm{m/s}$ [13] while typically the propagation velocity for the action potential can be much lower than these values. On a related note, in principle, sound velocity in axoplasm which is composed of approximately 80% of water could be higher than the typical propagation velocity of an action potential as well. However, in the experiment by Terakawa [32] the observed pressure changes were propagating with similar velocities to the propagating action potentials. The numerical results capture the ‘synchronization’ of velocities between mechanical and electrical pulses when tracking the pulse peaks. However, it is also evident that under the considered parameters the wave front of the mechanical wave propagates at higher velocity (the bounding velocity in Fig. 2) under the presented parameter combinations than the action potential. The mechanical effects related to the action potential propagation in axons and lipid bi-layers have been studied in a number of experiments [10, 18, 31] and the numerical results are qualitatively similar, although, as noted earlier, the results presented here are in the dimensionless form.

The model can also include directly the pressure waves in the axoplasm [32]. In this case the improved Heimburg-Jackson model for the mechanical wave would contain an additional coupling term as a function of the local pressure [9]. The framework can be further extended to include thermal effects [30].

Finally, the combined model presented here results in a wave ensemble. However, the model has a number of properties which invite further investigation and analysis. The key phenomena observed in the coupled model are: (i) the velocity of the peak of the mechanical pulse is similar to the velocity of the action potential regardless of the sound velocity value in the lipid bi-layer, (ii) the velocity of the front and the shape of the mechanical wave depends on sound velocity in the lipid bi-layer, (iii) the shape of the mechanical wave can have an effect on the velocity and shape of the action potential.

Acknowledgements

This research was supported by the European Union through the European Regional Development Fund (Estonian Programme TK 124) and by the Estonian Research Council (projects IUT 33-24, PUT 434). The paper reflects the ideas of a talk at The Tenth IMACS International Conference on Nonlinear Evolution Equations and Wave Phenomena: Computation and Theory, March 29 – April 01, 2017. Georgia, USA.