Modelling of processes in nerve fibres at the interface of physiology and mathematics

TL;DR

This paper develops a coupled mathematical model integrating electrical, mechanical, and thermal effects in nerve fibers, bridging physiology and physics to better understand nerve pulse propagation.

Contribution

It introduces a novel coupled system of equations based on physical principles that models energy transduction and wave interactions in nerve fibers.

Findings

The model generates combined electrical, mechanical, and thermal wave ensembles.

It provides a framework for analyzing energy transduction in nerve fibers.

The approach unifies physics, physiology, and mathematics in nerve modeling.

Abstract

The in silico simulations are widely used in contemporary systems biology including the analysis of nerve pulse propagation. As known from numerous experiments, the propagation of an action potential is accompanied by mechanical and thermal effects. This calls for an analysis at the interface of physics, physiology and mathematics. In this paper, the background of the model equations governing the effects in nerve fibres is analysed from a physical viewpoint and then discussed how to unite them into a system by using the coupling forces. The leading hypothesis associates the coupling to the changes of variables, not to their values or amplitudes. This hypothesis models actually the physiological mechanisms of energy transductions in a fibre. The general assumptions in modelling the processes and the properties of the coupled system of equations are described. The dimensionless mathematical model which couples the action potential with mechanical waves together with temperature effects is presented in the Appendix. This model generates an ensemble of waves including the electrical signal and mechanical and thermal effects.

Equations

$$\begin{array}{rl}{Z}_T& =D{Z}_{XX}-J+Z\left(Z-G_1-Z^2+G_{1}Z\right);\\ J_T& =\epsilon \left(G_{2}Z-J\right),\end{array}$$

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

$$\begin{array}{rl}{U}_{TT}& =c^2{U}_{XX}+NU{U}_{XX}+M{U}^2{U}_{XX}+N{U}_X^2+\\ & +2MU{U}_X^2-H_1{U}_{XXXX}+H_2{U}_{XXTT}+F_2(Z,J,P),\end{array}$$

$$W\propto U_X,$$

$${\Theta }_T=\alpha {\Theta }_{XX}+F_3(Z,J,U),$$

$$\begin{array}{rl}Z(X,0)& =A_z{\mathrm{sech}}^2{B}_{o}X;Z_T(X,0)=0;\\ Z(X,T)& =Z(X+2Lm\pi ,T);m=1,2,\dots ,\end{array}$$

Full text

††journal: arXiv

1 Introduction

Signal propagation in nervous systems is extremely important for all multicellular animals including human beings. Nerve signals control motion, behaviour and consciousness in many respects. The rich history of studies into the function of nerves is described in many overviews [7, 34], etc and a concise history of neuroscience spiced with mathematics is given by Scott [38]. Contemporary studies are characterized by attempts to describe interactions between the biological and physiological functions of a nervous system across multiple scales [36, 22]. The better knowledge about functional behaviour of normal signalling in nerves helps to understand also the axon dysfunction, neuronal communication, synaptic efficacy, etc needed for understanding neuronal activity [12]. For this purpose, in silico studies of biological processes are gaining more and more attention [35]. It means developing interdisciplinary approaches combining ideas of physiology, physics and mathematics in studies on nervous function.

In this paper, a possible approach is envisaged for mathematical modelling of the signal propagation in axons. In Section 2 the general background of basic model equations is analysed followed in Section 3 by physiological description of signals in nerve fibres. The assumptions made in the modelling are presented and analysed in Section 4. Further, in Section 5, a concrete mathematical model based on those assumptions is described. The discussion in Section 6 summarizes the analysis and gives the hints for further studies. The mathematical model in the dimensionless form is given in the Appendix.

2 Mathematical modelling and basic equations

Complexity of biological processes involving structural and functional hierarchies needs interdisciplinary approaches in order to model phenomena under consideration. This is why systems biology uses methods and tools from physics and mathematics. Sometimes such an approach is called bio-mathematical modelling [22], sometimes a physical scientist’s perspective [11], sometimes simply working at the interface of physics and biology [6]. The mathematical modelling of biological processes means casting the physiological phenomena and functional behaviour into the language of mathematics. One could recall the saying of Galileo Galilei: ”The Book of Nature is written in the language of mathematics” but using this idea is far away of being simple. Indeed, biological systems need energy exchange with the surrounding environment and represent the systems far from the thermodynamic equilibrium; they involve many transfer mechanisms that should be analyzed on the molecular level; the process may involve different time and space scales; different types of mathematical equations are needed for describing the system; different physical effects like nonlinearities, dissipation, excitability, spatiotemporal coupling, etc must be taken into account [43]. However, mathematical modelling is not only a tool for describing biological processes but also extremely useful for performing in silico experiments. It is said that “simple conceptual models can be used to uncover new mechanisms that experimental science has not yet encountered“ [44]. In addition, a properly chosen mathematical model has a predictive power and helps to understand the causality of the process. Without any doubt, all mathematical models need experimental validation.

The mathematical models of dynamical processes have a physical background and are derived from basic principles: conservation of mass, motion and energy together with thermodynamical conditions. Dealing with motion, the wave equation is the backbone of all the models – one of the most important equations of mathematical physics [39]. Another important equation is the diffusion equation governing the distribution of heat (or variation in temperature). Both of these classical equations of mathematical physics are well studied but too simple for describing most practical cases and need modifications in order to grasp effects which are characteristic to phenomena under consideration [43, 15]. So the mathematical models of biological processes are usually derived from those simple equations with many specific additional terms and conditions in order to better reflect realities but certainly related to basic principles. And one idea more: “simple conceptual models can be used to gain insight, develop intuition, and understand “how something works”” [44]. Below such a modelling is demonstrated for the case of action potential in a nerve fibre.

3 Signals in nerve fibres, physiology and physics

The basic component in the process of signal propagation in nerve fibres is the electrical signal called action potential (AP). Since the pioneering studies of Hodgkin and Huxley [25], the electrophysiological model has been developed in many details [10, 9, 12] etc. The process is physically complicated due to special structure of nerve fibres. Bearing in mind unmyelinated axons, a nerve fibre in mechanical terms is a tube in a certain environment. The wall of the tube has a bilayered lipid structure – biomembrane – which can be considered as a microstructured medium. Inside the tube is the axoplasmatic fluid (axoplasm), called also intracellular fluid. In mechanical terms axoplasm is a viscoelastic medium containing also cytosceletal filaments. The biomembrane has voltage-sensitive (also mechano-sensitive) ion channels through which the ions can pass from the intracellular fluid to the environment (extracellular fluid) and vice versa. The ion concentration has an important role in emerging of an AP. The biomembrane as the cornerstone in dynamics of cells, plays an important role in biology in general, not only for nerve fibres [32].

A signal in such a complex structure has many components and there are many experimental studies demonstrating that an AP is accompanied by mechanical effects and temperature changes [41, 42, 1, 27] etc. Kaufmann [30] has stated that electrical action potentials are inseparable from force, displacement, temperature, entropy and other membrane variables. Indeed, the AP generates also the pressure wave (PW) in the axoplasm and a longitudinal wave (LW) in the biomembrane. The LW means a local longitudinal compression which changes the diameter of an axon (swelling) reflected by a transverse displacement (TW). This process is characterized also by heat production and temperature ($\Theta$) changes. It is a real challenge as said by Andersen et al [2]: to frame a theory that incorporates all observed phenomena in one coherent and predictive theory of nerve signal propagation. Such a theory depends on many physical properties of an axon and extracellular fluid. The AP depends on the properties of the axoplasm and ion currents, the PW depends on the properties of the axoplasm. The LW and the TW depend on the elastic properties of the biomembrane, while the $\Theta$ depends on the thermodynamical parameters of the whole system. All these quantities are coupled into a whole which reflects the complexity of the process of signal propagation in an axon.

Albert Einstein has said: “Everything must be made as simple as possible but not simpler”. In case of signal propagation in nerves, the main question is what is “simpler”? There are many physiological properties of axons [9, 12], etc which influence the process. First, all the structural elements have a certain microstructure (the lipid molecules of the biomembrane and the ion channels, the filaments and proteins in the axoplasm, etc) which play a role in signal formation and propagation. The coupling forces between the AP and the accompanying effects should reflect the physics of electro-mechanical transduction and be determined with a suitable accuracy. The thermodynamical properties of all the system should be taken into account. An overview by Drukarch et al [13] describes the present state-of-art of analysis of this fascinating process.

Although molecular mechanisms (ion currents) influence the process, the values of main variables and axon parameters are far from the molecular range. In terms of continuum theories this means the micro- or mesoscale [21] which justifies the usage of models derived from basic conservation laws of continua. Starting from the basics (see Section 2), the wave equations are the cornerstones of all the processes which propagate in media with certain finite velocities. The diffusion equations describe the processes related to heat production. Given the complicated structure of an axon, these model equations should be modified.

4 Assumptions in modelling

In what follows, an attempt is described following the ideas of Andersen et al [2] to build up a mathematical model which could account for main effects in signal propagation in nerve fibres. Based on classical understandings of axon physiology [9, 12], the following basic assumptions are made:

(i) electrical signals are the carriers of information [12] and trigger all the other processes;

(ii) the axoplasm in a fibre can be modelled as fluid where a pressure wave is generated due to electrical signal; here, for example, the actin filaments in the axoplasm may influence the opening of channels in the surrounding biomembrane but do not influence the generation of pressure wave in the fluid [2];

(iii) the biomembrane is able to deform (stretch, bending) under mechanical impact [24];

(iv) the channels in biomembranes can be opened and closed under the influence of electrical signals as well as of the mechanical input; it means that tension of a membrane leads to the increase of transmembranal ion flow and the intracellular actin filaments may influence the motions at the membrane [32, 4];

(v) there is strong experimental evidence on electrical or chemical transmittance of signals from one neuron to another [5, 26] although the role of mechanical transmission is also discussed [3].

A different scheme of signal generation is proposed by Heimburg and Jackson [24] where the mechanical wave in the biomembrane (LW) has the leading role.

The trivial idea for constructing a mathematical model is to collect all the single equations describing the processes in a nerve fibre but these should be coupled into a joint system. The question how the coupling forces between the equations (i.e., signal components) are modelled is crucial. It is not only a mathematical construct but the forces reflect the electro-mechanical (or mechano-electrical) transduction mechanisms.

We introduce the main hypothesis [19]: all mechanical waves in axoplasm and surrounding biomembrane together with the heat production are generated due to changes in electrical signals (AP or ion currents) that dictate the functional shape of coupling forces. The seconding hypothesis is: the changes in the pressure wave may also influence the waves in biomembrane.

It must be noted that already in the 19th century the German physiologist Emil Du Bois-Reymond has noticed that “the variation of current density, and not the absolute value of the current density at any given time, acts as a stimulus to a muscle or motor nerve” [23]. This statement is called the Du Bois-Reymond law.

Two essential remarks are following: (i) the changes of variables mean mathematically either the space or time derivatives; (ii) the pulse-type profiles of electrical signals mean that the derivatives have a bi-polar shape which is energetically balanced. Based on these remarks, the functional shapes of coupling forces are proposed in the form of first-order polynomials of gradients or time derivatives of variables [19, 16].

Due to reciprocity, the mechanical effects could also influence the electrical properties of the membrane and consequently the behaviour of the AP (see Chen et al, [8]). In mathematical terms it means that the model should have the variable-dependent parameters. In the present model the feedback from mechanical effects to the AP is proposed through introducing the dependence on LW into the ion current. It is possible to include some scale factors into the coupling forces but such a proposal needs more studies in order to prove a need for such an assumption. As the dynamical changes in the axon diameter (the amplitude of the TW) are many orders less than the diameter itself [41], the electrical properties of the fibre are taken constant during all the process.

In addition it is assumed that the TW is derived from the LW like in mechanics of rods [37] – the transverse displacement is a space derivative of the longitudinal displacement. Concerning the posssible temperature changes, there are several possibilities to build up a mathematical model [40].

It must be stressed that the assumptions described above are logically consistent within the present framework of axon electrophysiology [9, 12] and experimental studies on accompanying effects. What is important – all the components of a wave ensemble including the AP are coupled. The model serves the contemporary trends in computational biology and could help experimentalists to check the possible transduction mechanism [44].

5 Model equations

The mathematical analysis of the AP is usually based on the celebrated Hodgkin-Huxley (HH) equation which is derived from a telegraph equation neglecting inductivity [25]. The wave character of this benchmark model is preserved by adding the influence of ion current, so the result is a reaction-diffusion type equation. Three phenomenological variables govern the “turning on” and “turning off” the membrane conductance and are described by relaxation equations. The values of all the parameters needed for calculations are determined by ingenious experiments. As a result the profile of an AP has a typical asymmetric shape with an overshoot responsible for the relaxation length between pulses. An AP can be generated by an input over a certain threshold value. Contemporary studies [9, 12] have demonstrated that this model, complicated as it is, should be made even more complicated for describing effects like neuronal timing and synaptic efficacy, etc. However, in order to describe a correct functional shape of an AP, a simpler model was proposed by Nagumo et al [33] which has taken only one ion current into account. The main properties of the AP described above, are preserved and that is why this model called after FitzHugh-Nagumo (FHN) is used in the robust model proved by many calculations like for example those by Engelbrecht et al [19, 16].

The second important component in the signal complex is the LW in the biomembrane. Due to the complicated internal structure (microstructure) of a biomembrane made of lipids, a simple wave equation is not sufficient for describing the wave motion but nonlinearity and dispersion must be also taken into account. Heimburg and Jackson [24] have derived a wave equation (HJ equation for short) where the compressibility of the lipid molecules has an impact on velocity and an ad hoc dispersive term describes dispersion. Engelbrecht et al [17] have shown that based on mechanics of microstructured media dispersive effects should be described by two terms – one like proposed by Heimburg and Jackson [24] reflecting elasticity and another term, reflecting inertia of lipid molecules. The governing equation for the LW in biomembranes is then described by the improved Heimburg-Jackson (iHJ) equation. The possible solutions of this Boussinesq-type iHJ equation are analysed by Engelbrecht et al [18]. The remarkable property of the HJ and iHJ equations is that they possess solitonic solutions [24, 17]. However, the influence of the inertia of lipids makes the soliton narrower. It must be stressed that the LW is generated by an AP and the formation of a soliton from such a special input takes time and depends on the energy of the input.

The third component in the signal complex is the PW in the axoplasm. As far as axoplasm can be treated as a viscoelastic fluid, the classical Navier-Stokes equations can be used. However, given the very small amplitudes of the PW [41, 42], nonlinear effects can be neglected and that is why the classical wave equation with a viscous term added can effectively be used.

These three components in the signal complex are all wave-like and possess definite finite velocities. In such a way these waves are the primary components in the ensemble of waves which are coupled into a whole by coupling forces modelled following the main hypothesis (see Section 4).

The remaining components (the TW and temperature $\Theta$) have no finite velocities of themselves and are derived from primary ones. This is the reason why these components in an ensemble are secondary [20].

The TW is derived from the LW by derivation following the theory of rods [37]. The derivation of a suitable model for temperature $\Theta$ generation is more complicated because there is no widely accepted consensus about the mechanism of the heat generation. The simulations show [40] that there are several mathematical models that could match various experimental results. It seems, however, that a diffusion equation with a source term based on the conservation of energy, is best suited for modelling the temperature change accompanying the AP. The source term like in the case of primary components, depends on changes of other variables.

The collection of all equations needed for describing an ensemble of waves in an axon involving primary and secondary components is presented in the Appendix.

6 Discussion

What is described above, is a straightforward modelling of the complex process in nerve fibres. Without any doubt, this is a robust model and serves first of all as a proof of concept. However, after demonstrating that the basic equations of mathematical physics in their modified forms can be used to build up the backbone of nerve pulse dynamics, the further studies could enlarge significantly the description of physiological effects. Here we have been working on the integrative level [36] trying to unite the different functional elements into a whole.

The main points in the proposed model can be listed as follows:

(i) the proposed model is a robust one but describes qualitatively correctly the profiles of essential signal components in a nerve fibre; as a matter of fact, an ensemble of waves propagates in the fibre;

(ii) the AP is described by the FHN model which takes only one ion current into account which is assumed to be similar to the sodium current [33]; a possibility to distinguish between the voltage-sensitive and mechano-sensitive ion channels [32] is accounted by different parameters governing the ion current;

(iii) the model at this step can be taken as a proof of concept and the possible next step in modifying the analysis is to use the HH model instead of the FHN model;

(iv) the important novel hypothesis for constructing the joint model is related to the coupling forces between the single components of the full signal;

(v) it is assumed that the coupling forces depend on changes of coupled signals (pulses) i.e., on their derivatives, not on their amplitudes;

(vi) by generating the AP from an initial input, the other waves are generated due to the coupling forces;

(vii) the properties of the TW which has been measured in several experiments depend on the mechanical properties of the biomembrane through the LW;

(ix) there are several possible ways to model the temperature changes accompanying the AP demonstrated in experiments.

It seems that the crucial problem in all mathematical models is related to the transduction of energy from one component in the ensemble to another. Several possibilities to model the coupling have been proposed. El Hady and Machta [14] have elaborated a mechanism of electro-mechanical coupling based on the assumption that the potential energy of the process is stored in the biomembrane and kinetic energy in the axoplasmic fluid resulting in mechanical surface waves. The HH model is used for describing the AP and the force exerted on the biomembrane is taken proportional to the square of the voltage while the axoplasmic fluid is described by the linearized Navier-Stokes equation. Another coupled model of electrical and mechanical signals based on spring-dampers (dashpots) system is proposed by Jerusalem et al [28]. Johnson and Winslow [29] have discussed on the basis of physiological effects, how to unite the HH model on the macroscopic level and the soliton model [24] on the microscopic level. In their approach the Na*+* ion current is responsible electro-mechanical interaction. Chen et al [8] have modelled the coupling of electrical and mechanical (transverse displacement) effects by using the idea of flexoelectricity. They have used the HH model for the AP and the linear viscoelastic Maxwell-type model for the biomembrane where due to the flexoelectric effect the AP generates a strain gradient bending the membrane upwards. This means that the AP is accompanied by a transverse displacement. The flexoelectric force depends on the change of the membrane potential. In this analysis the biomembrane is a homogeneous tube with elastic and viscous properties. The numerical calculations by means of the finite element method permit to solve direct (a mechanical effect generated by an AP) and reverse (an AP generated by the mechanical effect) problems. This means that reciprocity of processes is possible. The processes in both myelinated and unmyelinated axons were analyzed.

As stressed by Coveney and Fowler [11], coupled models across several processes could provide also the route for calculating many unknown parameters. That explains why systems biology turns so much attention to the modelling of biological complexity involving molecular and continuum approaches. The in silico experiments permit to cover a large area of possible physical parameters in order to find suitable sets verified by experiments in vivo or in vitro. This is also the case of nerve signals where the molecular structure of the biomembrane affects the process and the conjectured coupling forces need further studies and quantization. Although much is known about the fascinating process of nerve signalling, the full picture needs further concerted efforts of experimentalists and theoreticians. A recent detailed overview [13] reflects the contemporary insights in this field that could lead to “the formulation of a more extensive and complete conception of the nerve impulse”.

In modelling described above, the path was from structures to dynamics trying to grasp the leading effects in mathematical terms. Such an approach is certainly robust which is widely used in systems biology [31]. The further modifications of the model should turn to details like accounting of the various ionic currents, membrane channels fluctuations, general oxygen consumption, CO2 output, the influence of carbohydrates and membrane proteins, etc together with elaborating the mechanisms of heat production. This will enhance the predictive power of the model which at this stage is related to qualitative estimations. These estimations correspond to known experimental results but need quantitative analysis.

The modelling at the interface of physiology and mathematics casts more light to the fascinating phenomenon of nerve pulse propagation. The structure of mathematical equations, especially of coupling forces helps to understand the causality of effects. Typically to biological processes, the coupling of various phenomena or fields must be taken into account. In this case one could modify the Du Bois-Reymond law: in dynamical processes every variation of fields acts as a stimulus to other fields. As explained above, mathematically the variation means the space or time derivatives of variables. Last but not least, as stressed by Konrad Kaufmann, one cannot forget physical background of biological processes. That is why for such a complex biological process as the propagation of signals in nerves the modelling associated with basic equations of mathematical physics is justified.

Acknowledgements

This research was supported by the Estonian Research Council (project IUT 33-24).

The mathematical model

The mathematical model in the dimensionless setting presented below has been elaborated and used for numerical simulations in earlier publications [17, 18, 19, 16, 20, 40]. Here and further subindex $X$ denotes partial derivative with respect to space and $T$ denotes partial derivative with respect to time. The primary (wave-like) components:

AP (the FHN equation):

[TABLE]

where $G_{i}=a_{i}+b_{i}$ and $b_{i}=-\beta_{i}U$. Here $Z$ – action potential, $J$ – ion current, $\varepsilon$ – the time-scales difference parameter, $a_{i}$ – “electrical” activation coefficient, $b_{i}$ – “mechanical” activation coefficient and $U$ – longitudinal density change from lipid bi-layer density model (3) and $\beta_{i}$ is an coupling coefficient.

PW (the wave equation with viscous and source terms):

[TABLE]

where $P$ – pressure, $\mu$ – viscosity coefficient. The $F_{1}$ is the coupling term accounting for the possible influence from the AP and TW.

LW (the iHJ equation with a source term):

[TABLE]

where $U=\Delta\rho$ is the longitudinal density change, $c$ is the sound velocity in the unperturbed state, $N,M$ are nonlinear coefficients, $H_{i}$ are dispersion coefficients. Where $H_{1}$ accounts for the elastic properties of the bi-layer and $H_{2}$ the inertial properties. The $F_{2}$ is the coupling term accounting for the possible influence from the AP and PW.

The secondary components:

TW (related to the LW):

[TABLE]

where $W$ – transverse displacement, $U$ – the LW.

$\mathbf{\Theta}$ (the diffusion equation with a source term):

[TABLE]

where $\Theta$ – temperature, $\alpha$ – thermal conductivity coefficient. The $F_{3}$ is the coupling term accounting for the possible influence from the AP and LW.

The initial condition is in terms of Z:

[TABLE]

while other variables have initially zero values and the same periodic boundary conditions. Here $A_{z}$ is the initial AP amplitude, $B_{o}$ is the pulse width and $L$ is the number of $2\pi$ sections in space.