On mechanisms of electromechanophysiological interactions between the components of nerve signals in axons

TL;DR

This paper reviews experimental evidence and proposes a mathematical model for electromechanophysiological interactions in axons, highlighting three physical mechanisms that couple electrical, mechanical, and thermal effects during nerve signal propagation.

Contribution

It introduces a modified mathematical model incorporating physical mechanisms of electromechanophysiological coupling in nerve signals, based on experimental insights.

Findings

Identification of three key physical coupling mechanisms.

Development of a mathematical model describing wave interactions.

Analysis of variable influences on nerve signal interactions.

Abstract

Recent studies have revealed the complex structure of nerve signals in axons. There is experimental evidence that the propagation of an electrical signal (action potential) is accompanied by mechanical and thermal effects. In this paper, first an overview is presented on experimental results and possible mechanisms of electromechanophysiological couplings which govern the signal formation in axons. This forms a basis for building up a mathematical model describing an ensemble of waves. Three physical mechanisms responsible for coupling are (i) electric-lipid bi-layer interaction resulting in the mechanical wave in biomembrane; (ii) electric-fluid interaction resulting in the mechanical wave in the axoplasm; (iii) electric-fluid interaction resulting in the temperature change in axoplasm. The influence of possible changes in variables which could have a role for interactions are analysed and the concept of internal variables introduced for describing the endothermic processes. The previously proposed mathematical model is modified reflecting the possible physical explanation of these interactions.

Equations

$$\begin{array}{rl}{Z}_T& =D{Z}_{XX}-J+Z\left(Z-C_1-Z^2+C_{1}Z\right),\\ J_T& ={ϵ}_1\left(C_{2}Z-J\right),\end{array}$$

$$P_{TT}=c_2^2{P}_{XX}-{\mu }_2{P}_T+F_2(Z,J),$$

$$\begin{array}{rl}{U}_{TT}& =c_3^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}-{\mu }_3{U}_T+F_3(Z,J,P),\end{array}$$

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

$$F_2={\eta }_1{Z}_X+{\eta }_2{J}_T+{\eta }_3{Z}_T,$$

$$F_3={\gamma }_1{P}_T+{\gamma }_2{J}_T-{\gamma }_3{Z}_T,$$

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

$$F_4={\tau }_1{J}^2+{\tau }_2\left(P_T+{\phi }_2(P)\right)+{\tau }_3\left(U_T+{\phi }_3(U)\right)-{\tau }_{4}K,$$

$$F_4={\tau }_1{Z}^2+{\tau }_2\left(P_T+{\phi }_2(P)\right)+{\tau }_3\left(U_T+{\phi }_3(U)\right)-{\tau }_{4}K.$$

$$K_T+{ϵ}_{4}K=\zeta J.$$

$$K_T={\phi }_4(J)-\frac{K-K_0}{{\tau }_K},\text{where}{K}_0=0,{\tau }_K=\frac{1}{{ϵ}_4},$$

$${\phi }_2(P)={\lambda }_2\int P_{T}dT\text{and}{\phi }_3(U)={\lambda }_3\int U_{T}dT,$$

$${\phi }_4(J)=\zeta \int JdT$$

Full text

††journal: arXiv

1 Introduction

Although the studies into the propagation of signals in nerves have a long history [1, 2], the research is going on. It is accepted that beside the electrical signal (action potential), the accompanying non-electrical manifestations like mechanical effects together with local temperature changes should also be taken into account in order to get a full picture on this fascinating phenomenon. There is strong experimental evidence [3, 4, 5, 6, 7, 8, 9] on accompanying effects and there are several mathematical models [10, 11, 12, 13, 14] proposed for describing a signal with such effects. However, there is no general consensus about the possible physical mechanisms governing the processes of propagation and coupling. A need for such a consensus is stressed by many researchers, from Hodgkin [15] to recent studies [16, 17, etc.]. In this paper, a brief overview on experimental results and proposed physical mechanisms of coupling is given and analysed within the framework of a robust mathematical model [18, 19]. As a result, on the basis of known experimental studies, the coupling forces are specified in more detail in order to reflect better the physics of the process. It is concluded that there are several mechanisms responsible for emerging an ensemble of waves in nerve fibres. Based on experimental results, these mechanisms are cast into the mathematical descriptions.

In Section 2 a brief overview of experimental studies is given. The possible physical mechanisms are analysed in Section 3. On the basis of this analysis, possible mathematical descriptions and coupling forces are proposed in Section 4. The mathematical model with coupling forces related to physical effects is described in Section 5 which enlarges our previous proposals [18, 19, 11]. Finally, the conclusions are presented in Section 6.

2 Brief overview of experimental studies

The contemporary understanding of processes in nerve fibres is based on experimental studies in axons. In simplified physical terms, axons are cylindrical tubes embedded into an extracellular fluid. The wall of the tube has a bilayered lipid structure called biomembrane formed by phospholipids (a head group and a fatty acid tail) which altogether is about 3-4 nm thick. Inside the tube is axoplasm (intracellular fluid) with cytoskeletal elements where the action potential (AP) propagates. The squid giant axon which is the classical experimental object has the diameter up to 1 mm but the usual mammalian axons have the diameter around 20 $\mathrm{\SIUnitSymbolMicro m}$. A more detailed description of axon morphology is given, for example, by Debanne et al. [20] and of biomembranes – by Mueller and Tyler [21].

The way to understanding the propagation of APs is paved by experiments of Hodgkin and Huxley [22]. They have measured the AP in the non-myelinated squid axon and explained the role of ionic (Na+ and K+) currents. The absolute amplitude of the AP was about 100 mV, duration (without the overshoot) about 1 ms and the velocity 18.9 ms*-1*. The contemporary measurements have demonstrated that the velocities of electric signals in nerves vary in a large interval (from ca 1 ms*-1* to ca 100 ms*-1*). The measurements of the AP together with the quantitative measurements of ionic currents [22] form the basis for the corresponding Hodgkin-Huxley (HH) model [15] which is nowadays called also the Hodgkin-Huxley paradigm.

Many experiments have shown that the propagation of an AP in a nerve fibre is accompanied by transverse displacements of the biomembrane which mean changes in the axon diameter [3, 4, 5]. These local changes called also swelling are small being in the range of 1-2 nm and compared with the diameter of fibres are of several orders smaller. This is confirmed by recent experiments by Yang et al. [6]. In addition to the deformation of the biomembrane, the pressure wave in axoplasm has been measured by Terakawa [7]. In a squid axon the amplitude of a pressure wave, measured simultaneously with the AP, was about 1 to 10 mPa. The temperature change during the passage of an AP is also measured, it is of the range about 20-30 $\mathrm{\SIUnitSymbolMicro}$K for garfish [4] and much less for bullfrog [23]. The earlier findings on mechanical and thermal effects are summarised by Watanabe [24], and more recently by Andersen et al. [16].

The transverse displacement $W$ of the cylindrical biomembrane is associated to the longitudinal displacement $U$ of the biomembrane. This effect – $W$ is proportional to the gradient of $U$ – is well understood in mechanics [25] for the theory of rods. It means that the bipolar $W$ measured by Tasaki [4] corresponds to the unipolar $U$ (and vice versa). The possible deformation of a biomembrane under loading is studied by measuring the transverse displacement [26, 27] and interpreted then as an accompanying mechanical wave along the biomembrane [12, 26]. This longitudinal wave may have a soliton-type shape [12], i.e., is unipolar. It must be noted that the excitable plant cells (Chara braunii) behave in a similar way: the electrical signal is coupled with a mechanical effect [28].

To sum up, there is strong experimental evidence about several effects accompanying the propagation of APs. We shall analyse next the proposed mechanisms of coupling.

3 Possible physical mechanisms of coupling

The whole ensemble of waves and effects is the following: (i) action potential AP and the corresponding voltage $Z$ associated with the ion current(s), here denoted by $J$; (ii) longitudinal wave LW in the biomembrane with the amplitude $U$; (iii) pressure wave PW in the axoplasm with the amplitude $P$; (iv) transverse displacement TW with the amplitude $W$; (v) temperature change $\Theta$. Note that above just one ion current is listed as a variable. The crucial question is: what are the physical mechanisms which link these signal components into a whole.

The starting point of the analysis is related to causality. The classical understanding in axon physiology is the HH paradigm: the whole process is electricity-centred and starts with generating the AP. However, according to the original HH model, the accompanying effects (see above) are not considered. Recently much attention is also paid to another paradigm: the whole process starts with generating the mechanical disturbance (LW) as proposed by Heimburg and Jackson [12]. The governing equation of the LW has a soliton-type solution describing the region of higher density in the biomembrane and that is why this approach is called soliton theory. The pros and contras of both approaches are analysed by Appali et al. [29] and Meissner [30].

Later in Section 4 we follow the HH paradigm as the fundamental approach in contemporary axon physiology [31, 20] but try to link it to the accompanying effects. This means that one should pay attention to electrical-to-mechanical and mechanical-to-electrical couplings and to possible heat generation.

Qualitative observations from experiments

The following observations can be noted from the published experiments: -The measured TW is bi-polar for the squid giant axon and its peak coincides with the peak of the AP [4], both have approximately the same duration.

-The measured TW is close to uni-polar for the garfish olfactory nerve and its peak coincides with the peak of the AP [5].

-The peak of the force developed at the axon surface coincides fairly accurately with the peak of the AP [4].

-The peak of the pressure PW lags behind the peak of the AP for the squid giant axon [7].

• Mechanical and thermal signals are in phase with the voltage changes [26].

• The shape and the width of a TW are similar to those of the measured AP (without the overshoot) for the rat neuron [6].

• The experiments with Chara braunii cells have demonstrated that the mechanical pulse (out-of-plane displacement of the cell surface) propagates with the same velocity as the electrical pulse and is (in most cases) of bi-polar nature [28].

• The AP and temperature $\Theta$ for the garfish olfactory nerve are almost in phase and the duration of the positive phase of heat production is very close to the duration of the depolarizing phase of the AP [4].

• The AP is narrower than the temperature change and the thermal response cannot be directly proportional to the change of voltage [8].

• Good correlation exists between the initial positive heat and the potassium (K+) leakage [9].

Possible mechanisms

The well-studied HH model describes the AP dependence on ion currents [15]. Based on experiments with the squid giant axons, the role of opening and closing of ion channels has been demonstrated resulting in the asymmetric shape of the AP with an overshoot, the existence of a threshold for generating the pulse and the refractory period. The HH model is analysed in many details as the fundamental model of axon physiology [20]. It has been shown that in addition to Na+ and K+ ion currents there may be many more ion currents regulating the AP properties, especially for human atrial APs [32, 33]. However, the HH model does not describe accompanying effects, although Hodgkin said [15]: “In thinking about the physical basis of the action potential perhaps the most important thing to do at the present moment is to consider whether there are any unexplained observations which have been neglected in an attempt to make the experiments fit into a tidy pattern”.

The physical mechanism governing the HH model is based on the flow of ions through the biomembrane upon the change of voltage in axoplasm. Contemporary understanding is that the ion channels may be voltage-gated like in the HH model but also mechanically-sensitive [21, 34]. This must be taken into account in building up a fully coupled model.

Gross et al. [35] have analysed electromechanical transductions in nerves. For electrical-to-mechanical transductions the mechanisms of electrostriction and piezoelectricity are analysed and argued that both mechanisms could predict the swelling effects. For mechanical-to-electrical transductions it is proposed that the stress-induced changes due to surface charges influence the intracellular electric field.

A promising mechanism for coupling the electrical and mechanical signals is the flexoelectric effect which is manifested in the deformation of the biomembrane curvature under an imposed electric field [36]. The flexoelectric effect is used by Chen et al. [13] for modelling the coupling of the AP and mechanical wave (TW). The classical HH model combined with cable theory includes density change in the biomembrane induced by the flexoelectricity. The biomembrane is modelled as an elastic (or viscoelastic) tube where the flexoelectric force is included into the governing equation. This force depends on the local change in the membrane potential. The changes in the axon diameter are taken into account and the system is reciprocal – it can be triggered either by an electrical pulse resulting in an AP or by a mechanical stimulus.

A coupled model based on the primary AP which generates all other effects is proposed by El Hady and Machta [14]. A Gaussian profile (a pulse) for an AP is taken as a basic without modelling. The assumption is made that in the fibre the membrane has potential energy, and axoplasmic fluid has kinetic energy. The idea is that a surface wave (meaning the surface of the fibre) is generated in the membrane and in the bulk field within the axon the linearised Navier-Stokes equations are used for calculating the pressure. The ensemble includes voltage (pulse, not a typical AP), radial membrane displacement TW and the lateral displacement inside the axon (i.e., the PW). The heat is assumed to be produced as additional release of mechanical energy, summing transverse changes in the diameter and lateral stretch. So the sequence: AP – mechanical waves – heat (temperature) is followed.

Rvachev [37] has proposed that the axoplasmic pressure pulse (PW) triggers all the process. The PW triggers the Na+ channels and the local HH voltage spike develops which in its turn opens the Ca2+ channels. Free intracellular Ca2+ activates then the contraction of filaments in the axoplasm which gives rise to the radial contraction in the lipid bilayer. The similar idea that a PW could cause the excitement is proposed by Barz et al. [38].

Based on experimental results, Terakawa [7] has suggested that the pressure PW arises either from a change in electrostriction across the axoplasm or from a change in charge-dependent tension along the axoplasm. He states that pressure response is correlated with membrane potential and not with the membrane current. A slight influence of electro-osmotic water flow to pressure response is detected.

Abbott et al. [39] discuss the heat generation in Maia nerves. They give three possible reasons for heat production: (i) the positive heat is derived from the energy released during the rising phase of the AP and the negative heat due to the absorption of energy during the falling phase of the AP; (ii) the positive heat is due to the interchange of Na+ and K+ ions and the negative heat represents the partial reversal of this interchange; (iii) the heat production is related to exothermic and endothermic chemical reactions. Richie and Keynes [8] have supported similar experimental data like Abbott et al. [39]. They also stated that energy of the membrane capacitor is proportional to the voltage square and that the thermal response cannot be directly proportional to the change of voltage.

Tasaki and Byrne [23] have analysed the heat production in bullfrog myelinated nerve fibres. They estimated theoretically the relation of voltage $V$ to temperature $\Theta$ in dependence of time constant $RC$ (capacity x resistance in operational amplifier). For $RC$ longer, they used $V$ related to $\Theta$ and for $RC$ shorter $V$ related to $d\Theta/dt$.

To sum up, there is no consensus about the coupling and transduction of energy. Some studies concentrate upon the electromechanical transductions, some – upon the coupling of an electrical signal and temperature. This way or another, the qualitative experimental observations serve as guides in modelling. According to the HH paradigm, the process is triggered by an electrical stimulus which generates the AP but it is also proposed that the stimulus could be of a mechanical character. This could be the LW [12] or the PW [37]. The model by Chen et al. [13] based on using the flexoelectric effect involves reciprocity of electrical and mechanical components of the process. The generation of the temperature during the process is associated either with the electrical signal [39, 23] or to mechanical effects [14]. However, the role of chemical reactions in producing temperature changes is also under discussion [39].

This way or another, in most of studies (both experimental and theoretical), the coupling is associated to local changes of fields which cause changes in the whole system.

4 Modelling of coupled signals and coupling forces

In order to overcome (or to unite) the possible differences in proposing the various mechanisms of interactions described in Section 3, the best way is to return to basics. This means starting from physical considerations. Indeed, the signal propagation in nerves is a dynamical process which in terms of continuum theories is far from the molecular range (although some molecular mechanisms may be of influence) and the amplitudes of the components (AP, PW, TW, $\Theta$) are in micro or meso-scale. The question is then what is the role of wave equations or diffusion equations which usually govern the physical dynamical processes, in nerve pulse propagation. The celebrated HH model is based on the telegraph equation where the inductivity is neglected but ion currents added resulting in a reaction-diffusion equation which predicts a finite velocity of the AP. Leaving aside the physiological details, the longitudinal wave in the biomembrane and the pressure wave in axoplasm are waves and the role of wave equations describing the process is obvious. And in physical systems where temperature effects are measured, the Fourier law and the diffusion equation are the cornerstones of the possible description. So the conservation laws, well-known in continuum theories could be used. In general formulation they all include the possible forces which correspond to given physical situations (see, for example, [40]) and may also include more variables for which some thermodynamical constraints must be satisfied [41]. The next step is to describe the forces in correspondence with governing mechanisms. Although for the modelling of signals in nerves, several mechanisms (see Section 3) have been proposed, it is a challenge to find the proper physical description of possible forcing. Certain flexibility is desirable because it has been mentioned in several studies (for example, in [39]) that there might be several mechanisms involved simultaneously in the signal propagation.

Consequently, the crucial issue is how to construct the forces and how to model the coupling. It is proposed that the main hypothesis for constructing these forces could be [18]: 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 hypotheses are: the changes in the pressure wave may also influence the waves in biomembrane and mechanical waves may influence the AP and ion currents. This means reciprocity between the signal components. In many studies referred to in Section 3, the changes are mentioned [35, 13, 7, 23] as reasons for interactions. And back to the history: the German physiologist Emil Du Bois-Reymond has noticed in the 19th century 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” [42]. This statement is called the Du Bois-Reymond law.

What is change? In mathematical terms changes mean either space ($X$) or time ($T$) derivatives of variables. This gives the glue for proposing the functional shapes of forces which at the first approximation could be described in the form of first-order polynomials of gradients or time derivatives of variables ($Z_{X}$, $J_{X}$, $U_{X}$, $P_{X}$, and $Z_{T}$, $J_{T}$, $U_{T}$, $P_{T}$). Here and further subscripts $X$ and $T$ denote partial derivatives with respect to space and time, respectively. Such an approach involves also certain flexibility in choosing the model.

Further on, a possible approach in modelling is envisaged based on ideas described above [18, 19, 11]. The processes which compose the leading effects in signal propagation can each be described by single model equations. In the coupled model these single equations are united into a system by coupling forces. For a proof of concept, the coupling could be modelled by a simpler approach: from the AP and ion currents to all the other effects. It has been proposed that the process could be divided into primary and secondary components [43]. The primary components are characterised by corresponding velocities and their mathematical models are derived from wave equation(s). These components are the AP, LW and PW. The secondary components are either derived from the primary components like TW or their models are derived from basic equations which do not possess velocities like temperature $\Theta$. In the latter case the diffusion-type equation could serve as a basic mathematical model.

On the basis of experimental studies (Section 3) it seems plausible that there are several physical mechanisms of coupling. This concerns electrical-mechanical (AP to PW and AP to LW) and electrical-thermal (AP to $\Theta$) transduction. In essence, the mechanisms should also include feedback and coupling between all the components of the ensemble. However, these effects could be of more importance in pathological situations (axon dysfunction) or for detailed understanding of neural communication and neural activity in general.

As stated above, the coupling forces could be described by changes in variables. Note that gradients (space derivatives) act along the axon and time derivatives across the membrane. Based on thermodynamics (see Richie and Keynes [8], Tamm et al, [11]), for temperature changes one should consider also the possible effects of $Z$ or $Z^{2}$ (alternatively $J$ or $J^{2}$). In general terms the basic mechanisms might be the following:

(i) electric-biomembrane interaction resulting in a mechanical response (LW, variable $U$);

(ii) electric-fluid (axoplasm) interaction resulting in a mechanical response (PW, variable $P$);

(iii) electric-fluid (axoplasm) interaction resulting in a thermal response ($\Theta$).

However, as said above, the feedback from mechanical waves to other components of the whole ensemble might also influence the process. In more detail these interactions can be characterised:

(a) influence from the AP:

Pressure change in axoplasm – $Z_{X}$, i.e., AP gradient along the axon axis can influence the pressure as a result of charged particles present in axoplasm reacting to the potential gradient along the axon; 2. 2.

Pressure change in axoplasm – $Z_{T}$, i.e., potential changes across the biomembrane can lead to a pressure change proportional to the potential change through electrically motivated membrane tension changes (also taken into account for LW, and this is accounting for the influence from the mechanical displacement which is taken proportional to $Z_{T}$); 3. 3.

Pressure change in axoplasm – $J_{T}$, i.e., the axoplasm volume change from ion currents in and out of axon through the biomembrane plus the effect of possible osmosis which is assumed to be proportional to the ionic flows; 4. 4.

Density change in biomembrane – $Z_{T}$, i.e., electrically induced membrane tension change or flexoelectricity; 5. 5.

Density change in biomembrane – $J_{T}$, i.e., membrane deformation as a result of an ionic flow through the membrane (ion channels deforming surrounding biomembrane when opening/closing); note that ion channels are not modelled here in the FHN model explicitly; 6. 6.

Temperature change – $Z^{2}$ or $J^{2}$, i.e., temperature increase from current flowing through the environment (power); this effect is related to Joule heating; 7. 7.

Temperature change – endothermic term is dependent on the integral of $J$ which is taken to be proportional to concentration of reactants which decays exponentially in time after the signal passage; needs a new kinetic equation to be added.

(b) influence from the pressure P:

Density change in biomembrane – $P_{T}$, i.e., membrane deformation (displacement) from the local pressure changes inside the axon; 2. 2.

Temperature change – $P_{T}$, i.e., reversible local temperature change: when pressure increases then temperature increases proportionally and when pressure decreases then temperature decreases proportionally, this happens at the same timescale as pressure changes; 3. 3.

Temperature – the irreversible local temperature increase from energy consumed by viscosity (friction) – this is actually time integral of $P_{T}$.

(c) influence from the mechanical wave U in biomembrane:

Action potential change – $U_{T}$, i.e., ion currents are suppressed when membrane density is increased and amplified when membrane density is decreased; 2. 2.

Temperature change – $U_{T}$, i.e., the reversible local temperature change, when density increases the local temperature is increased proportionally and when density decreases the local temperature decreases proportionally, this happens at the same timescale as the density changes in biomembrane; 3. 3.

Temperature change – the irreversible local temperature increase from energy consumed by the added friction/viscosity term in longitudinal density change model – this is actually time integral of $U_{T}$ which is proportional to $U$.

These are the plausible processes which can influence the dynamics of nerve pulse propagation. However, not all of the listed effects might be energetically at the same scale as far as dynamics of the nerve pulses are concerned. At this stage, the possible influence of temperature changes might have on properties of fibres are disregarded as experimental observations are showing temperature changes low enough (up to about 50 $\mu$K) to be probably negligible for these processes [8, 23]. In addition, the possible changes in the fibre diameter which could change the physical properties are also disregarded being of many scales lower compared with leading effects. In general, the significance of the listed effects must be determined by experimental and theoretical studies in the future.

5 Mathematical model with specified coupling forces

We have proposed a mathematical model of nerve pulse propagation including the accompanying effects [18, 44]. This model served as a proof of concept and was based on using the basic equations of mathematical physics which were modified in order to reflect the physiological effects. The following models have been used: the reaction-diffusion model for describing the AP (the FHN model) and accompanying ion current; the modified wave equation for describing the LW in the biomembrane (the iHJ model); the wave equation with dissipation for the PW; the diffusion equation for the $\Theta$. The coupling forces [18] have been proposed to model the interaction between the components of the whole ensemble. The criticism on this model [45] is based on the seeming inconsistency of the coupling of the AP and LW described by indicated models. The main point of the criticism is that in the HH model the capacitance of the neural membrane is taken constant [45]. However, it must be noted that in the coupled model one has used the modelling of an AP only for obtaining a correct (measured by numerous experiments) asymmetric shape of an AP with an overshoot and a corresponding ion current. Note that the overshoot has a definite role in signal propagation [16, 20]. Certainly, other models could also be used for this purpose. For example, El Hady and Machta have used just a symmetrical Gaussian pulse without an overshoot for the AP [14]. The present understanding of electrophysiology give strong evidence of the importance of electrical signals [31, 20] in axons with an asymmetric shape.

In mathematical terms, such “simple conceptual models can be used to gain insight, develop intuition, and understand “how something works”” [46]. In physical terms, much must be understood before one could realise a complete unifying model. Following the HH paradigm [47], the causality relies on the propagation of an electrical signal in the axon. It must be noted that, alternatively, following Heimburg and Jackson [48] or Mussel and Schneider [49] the mechanical wave in a biomembrane can generate voltage changes. The manifestations of physical processes described briefly in Section 3 are suggesting the causality related to the AP but indicate the need to modify both the model equations as well as coupling forces for describing the accompanying effects. The main reason for the modification is the possible plurality of dynamical processes and coupling effects. Changes compared with the previous model [18, 19, 11] are following: dissipation is added to the iHJ equation; the coupling force for the temperature equation is specified accounting for multiple mechanisms which include also endothermic processes. The latter means proposing an additional kinetic equation involving the ion current. It means that the whole system involves the dissipative effects and the endothermic processes involve internal variables. Internal variables are not observable and characterise the internal structure or processes [50]. It has been used for describing some microscopic effect leading to local structural rearrangements [50] which have no inertial influence. Such a concept is widely used in contemporary continuum mechanics [52] and relies on thermodynamical description of dissipative effects. The first attempt to use internal variables in nerve pulse dynamics is described by Maugin and Engelbrecht [53] for describing the role of ion currents. Here such a formalism [54] is used for describing the possible chemical processes [39] which influence the temperature changes. The mechanicsm of such changes is proposed to be similar to the behaviour of phenomenological variables in the HH model – the change from one level to another.

For sake of completeness, the modified model is presented here. The FHN equations describing the AP are

[TABLE]

where $C_{i}=a_{i}+b_{i}$ and $b_{i}=-\beta_{i}U$. Here $Z$ is action potential, $J$ is ion current, $\epsilon_{1}$ is the time-scales difference parameter, $a_{i}$ is the “electrical” activation coefficient, $b_{i}$ is the “mechanical” activation coefficient and $U$ is a longitudinal density change from lipid bi-layer density model.

For the PW a modified wave equation with is used

[TABLE]

where $P$ is pressure, $\mu_{2}$ is a viscosity coefficient, $F_{2}$ is the coupling term accounting for the possible influence from the AP and TW.

Longitudinal wave in lipid bilayer is modelled by the improved HJ model

[TABLE]

where $U=\Delta\rho$ is the longitudinal density change, $c_{3}$ is the sound velocity in the unperturbed state, $N$,$M$ are nonlinear coefficients, $H_{i}$ are dispersion coefficients and $\mu_{3}$ is the dissipation coefficient. Note that $H_{1}$ accounts for the elastic properties of the bi-layer and $H_{2}$ the inertial properties. The $F_{3}$ is the coupling term accounting for the possible influence from the AP and PW. The transverse displacement (TW) is $W\propto U_{X}$.

The local temperature change $\Theta$ is governed by

[TABLE]

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

The coupling term $F_{2}$ is

[TABLE]

where $Z_{X}$ accounts for the presence of charged particles under the influence of the potential gradient (along the axon), $J_{T}$ accounts for the ionic flows into and out of an axon (across the membrane) and $Z_{T}$ accounts for the possible pressure change as a result of membrane tension changes from the electrical field. Contact term $F_{3}$

[TABLE]

where $P_{T}$ accounts for possible membrane deformation because of pressure changes (pressure to TW to LW), $J_{T}$ accounts for possible membrane deformation as a result of ionic flows through ion channels and $Z_{T}$ accounts for possible electrically induced membrane tension change. Note the sign, assuming that density decreases with the increasing tension.

For $F_{4}$ the following values were used in [11]: $F_{4}=\tau_{1}Z$ or $F_{4}=\tau_{2}Z^{2}$ or $F_{4}=\tau_{3}J$ or $F_{4}=\tau_{4}J^{2}$ or $F_{4}=\tau_{5}U$ or $F_{4}=\tau_{6}U^{2}$ or $F_{4}=\tau_{7}Z_{T}+\tau_{8}J_{T}$ or $F_{4}=\tau_{9}J_{T}+\tau_{10}U_{X},$ etc. These terms were considered from the viewpoint of possible mathematical description. However, it is possible that Eq. (4) needs to be improved in order to match better experimental observations when the time-scale for the temperature changes is larger than the time-scale for AP changes. Note that Abbott et al. [39] have argued about the possible combination of many reasons for heat production. Consequently, for modelling the temperature changes one could use

[TABLE]

where

[TABLE]

or

[TABLE]

Here $K$ is an internal variable, modelling all endothermic processes. In principle, it describes the change from one level to another like the phenomenological variables describing ion currents in the HH model [22]. The corresponding evolution equation for $K$ in Eqs (8) and (9) is [52, 53]

[TABLE]

An alternative way of representing the evolution of $K$ could be

[TABLE]

where $K_{0}$ is equilibrium level (in our setting equilibrium has been taken as zero level) and $\tau_{K}$ is the relaxation time. Clearly, $K$ as an internal variable has a certain relaxation time before reaching its equilibrium value. Following the discussion of Abbott et al. [39], the thermal influence of chemical reactions might initially be endothermic and in the later stages of the recovery – exothermic. In principle, this internal variable can describe both endo– and exothermic chemical influence, which means that the value of K has to change the sign (overshooting equilibrium value). However, this might need several ion currents to be taken into account.

In Eqs (8) and (9)

[TABLE]

where integrals in (12) characterise energy transfer from mechanical waves into the temperature increase as a result of friction–type terms in Eqs (2) and (3) and $\lambda_{i}$ are coefficients. In dimensionless case $\lambda_{i}=\mu_{i}$.

In Eqs (8) – (10) $\tau_{i},\lambda_{i},\zeta,\epsilon_{i}$ are coefficients. The term $J^{2}$ in Eq. (8) (or $Z^{2}$ in Eq. (9)) is either related to $\mathcal{P}=J^{2}r$ or $\mathcal{P}=Z^{2}/r$ where $r$ is resistance and $\mathcal{P}$ is power accounts for the Joule heating from the (electrical) current flowing through the environment. The term $P_{T}$ accounts for the local temperature increase if the pressure increases and decrease when pressure decreases (reversible), term $\int P_{T}dT$ accounts for the temperature increase from the energy lost to viscosity in Eq. (2), term $U_{T}$ accounts for the local temperature increase when the local density is increased and decrease if the local density decreases (reversible), term $\int U_{T}dT$ accounts for the temperature increase from the energy lost to viscosity in Eq. (3) and finally, as noted, the term $K$ takes into account the temperature decrease from endothermic processes. In Eq. (10), term $J$ is used on assumption that endothermic process intensity is proportional to

[TABLE]

which is taken as the time integral of ion currents at the location. This concentration-like quantity is decaying exponentially in time which is accounted by the term $-\epsilon_{4}K$ in Eq. (10). In other words, we are considering some kind of endothermic chemical reaction. Some other processes are outlined in the literature [51, 49] but not considered here, as possibly relevant, like, the phase change of the lipid bi-layer. In the present framework the phase change, if relevant, could be accounted by the term $U_{T}$ when using an alternate physical interpretation where change in the membrane density is a result of the phase change. However, unless noted otherwise the interpretation used in the present framework for the $U$ is that this is the ‘normal’ density change as a result of mechanical wave propagating in elastic environment.

A block diagram of the presented model is shown in Fig. 1, where the main components of the model framework and interactions between previously established individual component models are highlighted. The proposed framework is highly flexible and any of the components can be replaced with an improved model with the desired descriptive detail, like, for example, one could use an HH model instead of the FHN model without having to change the framework significantly.

Example results for the cases where only a single source of thermal energy is present is shown in Fig. 2. The profiles for the $Z,P,U$ are the same but which is different is the thermal response $\Theta$ curves. Qualitatively what is relevant here is the location of the thermal response peak compared to the wave ensemble and the shape of the initial thermal response curve.

Two distinct thermal response cases are shown in Fig. 3. It must be noted that the PW and LW profiles are also significantly different. A case with full $F_{3}$ in the form of $\tau_{1}J^{2}+\tau_{2}\left(P_{T}+\varphi_{2}(P)\right)+\tau_{3}\left(U_{T}+\varphi_{3}(U)\right)-\tau_{4}K$ and the viscosity in $P$ (2) and $U$ (3) is shown in left panel in Fig. 3. In the right panel in Fig. 3 the irreversible thermal processes have been significantly suppressed by setting coefficient $\tau_{1}$ to zero (related to Joule heating) and reducing the viscosity by two orders of magnitude compared to the case in the left panel. This means that the thermal response in Fig. 3 right panel is dominated by the reversible processes (temperature increase when density increases in axoplasm and membrane and temperature decrease if the density decreases with friction related terms suppressed). It can also be noted that viscosity or lack of it has a significant effect on the shape of the PW and LW profiles under the used parameters. The governing equations for the PW and LW are conservative before adding the contact forces $F_{1}$, $F_{2}$ and the dissipative terms. Without accounting for some kind of friction-like term in PW and LW the different sound velocities in the corresponding environments have an opportunity to spread the wave-profiles out as the driving force from the AP is taken with a different velocity than the sound speed in axoplasm and biomembrane. While not demonstrated in the figures, it should be noted the thermal signal drop off rate is dependent on the relaxation time $\tau_{K}$. If relaxation is taken fast enough it is possible to model thermal response curves which are almost perfectly in-phase with the propagating wave ensemble. Some endothermic chemical processes in the context of the nerve pulse propagation have been discussed by Abbott et al. [39]. Alternatively thermal response in-phase with the propagating wave ensemble can emerge if the reversible thermal changes are dominating (Fig. 3 right panel).

Figure 4 demonstrates the effect of Joule heating as a function of $Z^{2}$. In the left panel the thermal source and thermal sink have been balanced $\tau_{1}Z^{2}\approx K$ so that the thermal response settles eventually back to its starting value. In the middle panel the balance between thermal source and thermal sink terms is varied, demonstrating that depending on the balance between these terms there is a different residual level of heat increase after the nerve pulse wave ensemble has passed. In the right panel the thermal source and thermal sink are balanced and the relaxation time is varied. Demonstrating that depending on the relaxation time related parameters, the equilibrium or residual level can be reached at different rates. In terms of qualitative thermal response curve characteristics there is no practical difference between formulating the Joule heating in terms of $Z^{2}$ or $J^{2}$.

The numerical scheme used (with minor changes) for the example solutions is described in [44, 11] and references therein. The parameters used in the present numerical examples are the following. The FHN equation (1): $D=1$, $a_{i}=0.2$, $\beta_{i}=0.025$, $\epsilon_{1}=0.018$. The pressure equation (2): $c_{2}^{2}=0.09$, $\mu_{2}=0.05$ and coupling force $F_{2}$ (5): $\eta_{1}=10^{-3}$, $\eta_{2}=10^{-2}$, $\eta_{3}=10^{-3}$. The iHJ equation (3): $N=-0.05$, $M=0.02$, $H_{1}=0.2$, $H_{2}=0.99$, $c_{3}^{2}=0.1$ and coupling force $F_{3}$ (6): $\gamma_{1}=10^{-2}$, $\gamma_{2}=10^{-3}$, $\gamma_{3}=10^{-5}$. Parameter $\mu_{3}$ in (3) is taken as $\mu_{3}=0.05$ unless noted otherwise. The improved heat equation (7): $\alpha=0.05$ and thermal source/sink function $F_{4}$ (8) (or $\eqref{f3eqz}$): $\tau_{4}=\epsilon_{4}/2=0.005$ and $\tau_{1,2,3}$ is variable parameter. In addition, in (12) $\lambda_{2}=\mu_{2}=0.05$ and $\lambda_{3}=\mu_{3}=0.05$ unless noted otherwise. In (10) $\zeta=0.005$ and $\epsilon_{4}=0.01$ unless noted otherwise. The parameters for the numerical scheme and initial conditions: $n=2^{11}$ (number of spatial grid points), $L=128$ (the length of spatial period is $L\cdot 2\pi$), $\Delta T=1$ and $T_{f}=800$ (integration time in dimensionless units). For the initial condition $A_{z}=1.2$ (amplitude of the initial excitation for $Z$) and $B_{z}=1$ (width of initial sech2-type excitation for $Z$) while all other initial conditions are taken zero ($P,U,\Theta$ at rest).

6 Discussion

What has been described above, is the modelling of a complex process of propagation of signals in nerve fibres at the interface of physics, mathematics and physiology [55]. The physical background of wave motion is taken as a basic phenomenon, then processes formulated in the mathematical language and finally based on physiological understandings on the single effects coupled into a whole. In some sense we have followed A. Einstein [56] who said: “…mathematical construction enables us to discover the concepts and laws connecting them which give us the key to the understanding of the phenomena of Nature”.

An overview on experimental physical effects permits to describe the possible mechanisms of interaction between electrical and mechanical signals accompanied by temperature changes. Based on our hypothesis on the role of changes in electrical signals as sources for accompanying effects [18], the full list of possible factors are described and the coupling forces formulated. Compared with our previous studies, the number of possible coupling effects is enlarged and a new internal variable is introduced in order to account for possible endothermic reactions which according to experiments [39] could influence the temperature changes. Like the phenomenological variables in the HH model, the change of this internal variable involves a certain equilibrium level and the corresponding relaxation time. In order to justify the energy transduction, the dissipative term is included into the governing equation of the LW. Although we have followed the HH paradigm by starting the process by an electrical signal, the possible feedback is taken into account by coupling forces.

So in the full coupled model the primary governing equations which are related to wave motion have a physical basis derived from basic principles of continuum theory. The non-wave-like phenomena are governed either by diffusion equation (temperature change) or derived from the primary variables (this concerns the TW). The numerical results described in Section 5, are quantitatively similar to experimental profiles of changes (see references in Section 5).

The presented model is certainly robust and could be improved or modified in many aspects. As explained in Section 5, the FHN model has been used just for obtaining a qualitatively correct shape of an AP with an overshoot and the corresponding ion current. If the HH model with many physical constants would be used then there would be a possibility to specify the role of various ion currents. This could modify the generation of accompanying effects and may explain the retardation of refraction and also the influence of anaesthetics. In the proposed model the mechanical activation of ion channels in the FHN model is described phenomenologically by means of various coefficients. Such a proposal might also need refinements by specifying activation forces from lipids (influence from the biomembrane) or filaments (influence from the axoplasm) [57]. The LW is governed by iHJ equation which is an excellent description of a mechanical wave in a biomembrane [12]. However, one should note that although this equation has a soliton-type solution, its generation from an input needs time [58] and the time-scale might not be large enough for a real soliton to be emerged during the AP propagation. In addition, the existence of ion channels in biomembranes is a clear sign for inhomogeneities in the lipid bi-layer which could also influence the mechanical wave. This is again a challenge for further improvement of the model. The molecular effects related to ion movement, the influence of membrane proteins and the existence of filaments in the axoplasm may have a considerable effect to coupling forces. At this stage, the possible role of chemical reactions accompanying the whole process is taken into account by an internal variable related to possible temperature changes.

It is clear that every theoretical model should be validated by experiments. The mathematical modelling described above may be used to uncover new mechanisms that experimental science has not yet encountered [46]. Such an approach is justified from the viewpoint of physical considerations which must be correlated with physiological phenomena. The possible influence of various coupling forces depending on changes of variables demonstrated above (Section 5) indicates the possibilities for interpretation of experiments.

The explanations on coupling forces given above actually represent a rather full theoretical framework of possible interactions. Although there is a doubt whether it is possible to build up “general unifying models” [45], one cannot overlook the mathematical description of the physical background because the signal propagation in nerves is a complex electro-mechano-thermo-physiological process spiced with molecular effects resulting in an ensemble of waves. Surely the different models used so far for the description of single effects may have different levels of idealization but this is not a fundamental obstacle. Every attempt to understand this complexity is welcome and needs a close collaboration of theoreticians and experimentalists. In other words, it is a challenge to put the pieces back to the whole. The proposed model is still far from being ideal but forms a backbone in modelling based on physical considerations of wave motion expressed in the mathematical language.

Acknowledgment

This research was supported by the European Union through the European Regional Development Fund (Estonian Programme TK 124) and by the Estonian Research Council (IUT 33-24).