On the complexity of signal propagation in nerve fibres

Biological systems are characterized by many interwoven processes over multiple scales where interactions between the various phenomena play a decisive role. These phenomena can be chemical, electrical, and/or mechanical, all embedded into a whole. This paper discusses the propagation of signals in nerve fibres, which exhibits clear signs of complexity: electrical signals (action potentials) are coupled to mechanical waves in the internal axoplasmic fluid and in the surrounding biomembrane. Here the underlying microstructure affects strongly the processes: the existence of ion currents changes the balance of ions within fibres, the opening of ion channels in the surrounding biomembrane is a crucial process, and the biomembrane itself has a microstructure composed of lipids. The whole process is governed by interactions, and the analysis of single processes has demonstrated the importance of nonlinearities. The main challenge is to build up a general model where the coupling of all related phenomena is taken into account. It is proposed that three processes – the propagation of an action potential and mechanical waves in the biomembrane and in the axoplasmatic fluid – be united into a general model with additional interaction forces for reflecting coupling. Such a model results in the emerging of a mutually interacting ensemble of waves. The preliminary numerical simulations cast light onto the possible validation of this general model reflecting the complexity of signal propagation in nerve fibres. The mathematically consistent modelling will allow not only the prediction of process characteristics but gives a possibility of understanding the role of governing factors in the whole complex process.


INTRODUCTION
Complex systems are characterized by interactions of their components. This is the main reason why a system behaves differently from the simple sum of the behaviours of its components. Such a behaviour was already known by Aristotle. In contemporary studies there are many examples of processes or phenomena where interactions bring about a new quality at the macrolevel. In this context, scaling and hierarchies constitute a framework while interactions rule the complex macrobehaviour. The present understanding of complexity started with ideas of L. van Bertalanffy and N. Wiener (mid-20th century) and developed in studies on self-organization, chaos theory, networks, etc. in the second half of the 20th century (see, for example Prigogine and Stengers [1], Bak [2], Nicolis and Nicolis [3],Érdi [4], etc.). From the classical studies of interacting fields such as in thermoelasticity, photoelasticity, etc., nowadays much attention is paid to real-life systems whether in biology or man-made systems (multi-agent systems, econophysics, just to name a few examples) where nonlinearity, heterogeneity, hierarchies, and non-stationarity describe the governing properties or mechanisms. Again, in many such studies the focus is on interacting fields or interacting waves [5]. In these studies an important tool is the mathematical modelling of interactions. Mathematically consistent modelling will allow not only the prediction of process characteristics but gives a possibility of understanding the role of governing factors in the whole complex process [6,7].
In the modelling of wave phenomena it is possible to grasp even within a single governing equation many complex effects such as stable solutions formed as a result of the balance between nonlinear and dispersive effects, wave trajectories may form a certain pattern like Korteweg-de Vries (KdV) soliton trajectories or chaotic solutions may emerge under an external impact. One of the crucial problems in the modelling of wave phenomena is the proper modelling of the medium, which as a rule is not homogeneous but is formed by many constituents not only on the atomistic scale but also on the mesoscale like microstructured materials.
The question whether the microstructure should be determined with the same accuracy as the macrostructure or some assumptions should be made has to be answered clearly.
One of the possibilities is to model the microstructure using the concept of internal variables that model the microstructure as an additional field [8,9]. The modelling of wave phenomena in biological systems faces actually similar problems of scales, hierarchical microstructures, interactions of different processes, etc. while in general, biological phenomena are complex [10].
The effects of nonlinearity, microstructure, interactions, and excitability must be taken into account, and the physical and biological systems bear many similarities [11].
In this paper, attention is paid to the propagation of signals in nerve fibres. Despite the excellent results in electrophysiology starting with Hodgkin and Huxley [12] and FitzHugh [13], the more recent studies have revealed more information about the processes in nerve fibres [14] and have challenged researchers to build more consistent models that could account for many interacting processes accompanying the propagation of signals in nerve fibres. In Section 2, general principles of the mathematical modelling of complex physical and biological processes are briefly analysed for explaining the framework of further analysis. Section 3 is devoted to modelling the single processes in nerves in order to prepare the constituents for a more general model. In Section 4, the united model of signal propagation is described including the interacting electrical signal (action potential), the mechanical wave in the surrounding biomembrane, and the wave in the axoplasmic fluid inside the nerve fibre. As a result, an ensemble of waves propagates in a fibre influenced by nonlinearities and interactions. In the final section, the discussion is presented with the results of a numerical simulation, open questions, and perspectives for further studies.

COMPLEXITY
Here we very briefly list the main features of complexity and modelling used later in the analysis of signal propagation in nerve fibres. The signatures of complexity in physical systems are described in many monographs, see for example [3,4]. Starting from simple nonlinear cases, many important phenomena characterize the behaviour of complex systems and much can be learned from them. It is even surprising that very simple nonlinear systems such as the logistic equation or the three-body system display rich dynamics that help to understand more complicated cases. The main effects for understanding complexity are (i) non-additivity and nonlinear interactions; (ii) sensitivity to initial conditions; (iii) there are several typical phenomena characterizing the behaviour of nonlinear systems like bifurcations, emergence of patterns or ensembles, thresholds, coherent states, etc; (iv) deterministic unpredictability; (v) despite the variety of nonlinear motions, there are several rules that govern the processes. For more information one should consult the Encyclopedia of Nonlinear Science [15]. Also the following should be stressed. The usual understanding (common sense) is that nonlinear models are just a little bit corrected linear models. The world around us, however, is deeply nonlinear and the linear models, as a rule, are simplifications.
The analysis of wave motion in physical systems has revealed many signatures of complexity. Nonlinearity is an important property of solids and fluids influencing wave motion, and it may be balanced with other physical effects such as dispersion and dissipation [16]. As a result, solitons (coherent structures) or shock waves may be formed [16,17], ensembles of solitons may emerge [18,19], interaction of 2D solitons may lead to amplification [20], patterns of soliton trajectories demonstrate the character of motion in time [18,19], interaction with the embedded microstructure reveals the emergence of wave ensembles -solitary trains [21], seismic waves may undergo amplification [22,23], etc. Interaction of waves and fields such as in photoelasticity or acoustoelasticity shows explicitly how new wave structures may be formed and how informative they are (see overview by Engelbrecht [5]).
In biological systems the situation is even more complicated for several reasons [11,24,25]: (i) there is a need for energy exchange with the surrounding environment; (ii) many chemical reactions and transfer mechanisms that are often characterized on the molecular level are involved; (iii) different time scales, adaptivity, and hierarchies should be taken into account; (iv) in mathematical terms, the biological systems are described by different types of mathematical equations, which may cause problems in solving them; (v) nonlinearities, diffusive effects, activity, excitability, spatiotemporal coupling, etc. play usually a significant role [26][27][28][29][30][31].
The mathematically consistent modelling will allow not only the prediction of process characteristics but gives a possibility of understanding the role of the governing factors in the whole complex process and also causality. But the knowledge in other fields such as the theory of continua, thermodynamics, chemistry, computational methods, etc. should be interwoven into a whole on the basis of biological functions and physiology [24]. Whatever model is derived, it should always be validated by experiments and if necessary, modified and improved on the basis of validation.

Action potential
The famous Hodgkin-Huxley (HH) model is a basic description of the electrophysiological signal in nerve fibres [12]. It is based on the telegraph equation where the inductance is neglected and the ion currents added. In the original notations it reads: where v is the voltage; C is the capacitance; r s designates the resistance determined as 1/(πa 2 σ 1 ), where a is the radius of an axon and σ 1 is the conductivity; and I i is the ion current. Standard notations x = distance along the fibre and t = time are used. A phenomenological expression for the ion current is Here G K and G Na are the maximum potassium and sodium conductances, respectively; G L is a constant leakage conductance, C m is the membrane capacitance per unit area, v R is the resting potential, and v K , v Na , v L are the corresponding equilibrium potentials. The phenomenological (hidden in terms of continuum mechanics) variables n, m, and h govern the 'turning on' and 'turning off' of the membrane conductances. The ion current is a nonlinear function of three phenomenological (hidden) variables that govern the 'turning on' and 'turning off' of the membrane conductances. These phenomenological variables are described by the kinetic equations where the equilibrium values and relaxation times are described by complicated expressions based on experiments.
The HH model is a benchmark for action potential models, but further studies have improved it considerably by explaining the ion pumps in more details [32][33][34].
The simplified FitzHugh-Nagumo (FHN) model takes into account only one phenomenological variable and can be presented in the following form [13,15]: where z is the scaled voltage, j is the recovery current, D 1 is a coefficient, ε is the time-scale difference, the activation parameters satisfy conditions 0 < a < 1 and b > 0, and x and t are dimensionless space and time, respectively. Equations (4) can also be presented as one third-order partial differential equation [35]. Following the ideas of evolution equations (see for example [15,16], it is possible to derive also an evolution equation for the action potential [36]: where z 1 is the scaled voltage and ξ is a moving frame ξ = c 0 t − x. Following the FHN model with one phenomenological variable, f (z 1 ) is a quadratic function and g(z 1 ) is a linear function. Note that the moving frame involves velocity c 0 determined from the full telegraph equation. The final velocity c of a nerve pulse is different from c 0 . In fact, the ionic mechanisms are more complicated and include several pump, exchange, and background currents [37]. However, all the models based on telegraph type (electrical cable) theories do not take into account any other thermodynamical variables as well as possible mechanical effects [38]. Nevertheless, the nonlinear models described above are able to reflect the main features of measured action potentials: (i) the existence of a threshold for an input, (ii) the all-or-none phenomenon of the generated signal, (iii) the existence of an asymmetric localized pulse-type signal with an overshoot, (iv) the existence of a refraction length due to the overshoot, and (v) the possible annihilation of the counter-propagating signals [35,39]. However, being a solitary pulse, the action potential represents a wave in an active medium and is not a soliton [11]. The evolution equation (5) is a one-wave equation and is not able to describe the counter-propagating waves.
From the viewpoint of complexity, one has to stress the existence of a critical threshold, the formation of a coherent state, a special form of interactions, and the significant role of nonlinearities governing the process. In addition, the existence of phenomenological variables brings these models closer to microstructured materials where such hidden variables play a role in wave propagation [8,9].

Pressure waves in an axoplasmic fluid
The axoplasm within a nerve fibre is conmposed of 87% water held together by cytoskeleton [40]. The pressure waves in the axoplasm accompanying the propagation of an action potential were recorded by Terakawa [41].
In modelling, such a complicated structure is derived as a pseudoplastic fluid [42] or as a viscous compressible fluid [43,44]. This means that a pressure wave can be described by Navier-Stokes equations and in the first approximation by the 1D model while the biomembrane can be treated as an elastic tube. In this sense there is a similarity between the blood flow in the aorta [45] and pressure waves in the intersticial fluid.
The governing equation for 1D (along the x-axis) pressure waves is a momentum balance [46]. In terms of longitudinal velocity v x it reads: where ρ is the density,p is the pressure, µ is the viscosity, and F is the body force.
In terms of pressure in the 2D setting for waves in a fluid surrounded by a cylindrical tube (shell) the governing equations are [47]: , where x and r are cylindrical coordinates, c f is the velocity, and U x and U r are longitudinal and transverse displacements, respectively.

Mechanical waves in the surrounding biomembrane
Biomembranes are able to resist pressure, tension, stretch, and bending [48]. The deformation can be induced by a mechanical or an electrical impact [49][50][51]. The theoretical model of Griesbauer et al. [50] for pulses in lipid monolayers at the air-water interface is a second-order wave equation that is viscously coupled (by the first derivative) to the liquid (water) underneath. The existence of 2D sound waves along a lipid monolayer was demonstrated experimentally by Shrivastava and Schneider [51]. They showed that the localized selfsupporting pulses have a threshold amplitude and allor-none nature. The need to take nonlinearity in compressibility together with dispersive effects into account is stressed [51]. Such a nonlinear model with dispersion is proposed by Heimburg and Jackson [14], Andersen et al. [52], etc. The governing equation for a 1D pulse in a cylindrical biomembrane is proposed by using the energy density e [14] which results in the following equation of motion: where the last term is added in order to model dispersive effects [14]. Here u is the density change (u = ∆ρ 0 ), c 0 is the initial sound velocity, and p, q, h 1 are constants.
This is a Boussinesq-type equation where nonlinearity is caused by the compressibility, which has an impact on the velocity: In other words, the influence of the stiffness of the biomembrane caused by the repulsion potentials between the heads and tails of lipid molecules [53] is taken into account.
Due to the existence of nonlinearities and dispersion in Eq. (9), its solution may have a solitary character [14]. In order to improve dispersive properties of such a model, the use of more realistic dispersion terms like it is done for modelling waves in microstructured solids has been proposed [54]: The second dispersive term with the coefficient h 2 reflects the inertial properties of the lipid structure and it influences the width of the solitary pulse. In the improved model h 1 determines the limiting velocity at high frequencies and h 2 governs how fast this limit is reached.
In order to compare these theoretical results with experiments [55], one has to understand that Eqs (9) and (11) describe longitudinal waves and that the transverse displacements w measured by Tasaki [55] can be found from the derivative of the longitudinal profile [44,54,56]: where r is the radius of the axon and k is a constant. In the theory of rods k is the Poisson ratio. The nonlinear wave equations (9) and (11) describe longitudinal density waves in biomembranes and are able to model localized pulses due to the balance of nonlinear and dispersive effects. Several authors propose to call such models 'soliton models' for signals in nerve fibres [14,38,52,57]. As shown by Engelbrecht et al. [58], in the strict mathematical sense [17] these solitary pulses are not solitons because their interaction is not fully elastic. However, like in other Boussinesq-type models, it is customary to use the notion of solitons despite this discrepancy [59]. A full description of solutions of Eq. (11) over a large set of its parameters including the formation of solitary trains from an arbitrary initial input and the analysis of interactions is given in [58].

A COUPLED MODEL
The ideas of the coupling of physical effects in signal propagation in nerve fibres were already emphasized a long time ago [39,60] and more recently by Bennett [61] and Andersen et al. [52]. It is known, for example, that the ion channels may be voltagegated and also mechanically sensitive [48] and that an action potential generates an axoplasmic pressure wave in the intersticial fluid [43,44]. Both can in principle activate the ion channels in the surrounding biomembrane. Experiments by Terakawa [41] and Tasaki [55] demonstrated coupling effects. In general, one can say that the signal transmittance in neural systems is a complex physiological process [62,63]. It all calls for framing 'a theory that incorporates all observed phenomena in one coherent and predictive theory of nerve signal propagation' [52]. In terms of complexity, the goal is to formulate a model that will be able to describe an ensemble of waves of different physical origin (electrical and mechanical) in which the nonlinearities play a decisive role.
The simplest idea is to base the 'coherent theory' on existing mathematical models (action potential, the pressure wave in the axoplasmic fluid, and the mechanical wave in the surrounding biomembrane, see above) but introducing the coupling forces like it is done in continuum mechanics [64]. The following assumptions are made: (i) electrical signals are the carriers of information [62] and trigger all the other processes; (ii) the axoplasm in a fibre can be modelled as a fluid where a pressure wave is generated due to an 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 a pressure wave in the fluid [40,42,44]; (iii) the biomembrane is able to deform (stretch, bend) under a mechanical impact [50,51]; (iv) the channels in biomembranes can be opened and closed under the influence of electrical signals as well as of the mechanical input; this means that the tension of a membrane leads to an increase of the transmembranal ion flow and the intracellular actin filaments may influence the motions at the membrane [48,65]; (v) there is strong experimental evidence on electrical or chemical transmittance of signals from one neurone to another [61,66] although the role of mechanical transmission is also discussed [67]. The proposed simplified set of the coupled model is the following: The initial condition at t = 0 is where z is an electrical pulse. The action potential is governed by the FHN-type model (4) in the form of two coupled equations: where a = (a 1 + b 1 ) and b = (a 2 + b 2 ). Here a 1 , a 2 control 'electrical' activation and b 1 , b 2 control 'mechanical' activation. Indeed, as shown by Morris [68], the ion channels may be mechanosensitive. Compared with the original FHN model (4), the modification is introduced by differentiating electrical and mechanical activation. The pressure wave in the axoplasm is governed by a 1D Navier-Stokes model: where v is the velocity, ρ is the density, µ ν is the viscosity, and F 1 (z) is a force from the action potential.
In the biomembrane, the longitudinal wave is governed by where the notations follow Eqs (9) and (11), and F 2 (z, v) is a force from the processes in the axoplasm. In the biomembrane the transverse wave is governed by Eq. (12): where r is the axon radius and k is the coefficient which in the theory of rods is the Poisson ratio but in the present case, due to the complicated structure of a biomembrane, needs to be determined from experiments. The flowchart of such a model is depicted in Fig. 1. Clearly, the process depends on the microstructure of a nerve fibre: first, there are two structural elementsthe axoplasmic fluid and the surrounding biomembrane; second, the biomembrane has a distinct microstructure consisting of lipids. In order to describe the dynamics in such a complicated structure, the physical background has to be well understood. Physics of single elements is described from the viewpoint of the propagation of the action potential starting from studies of Hodgkin and Huxley [12] up to the contemporary understanding on axon physiology [62] with a lot of details. The physics of biomembranes is intensively studied not only because they are the elements of nerve fibres but because they are the building blocks of cells [48,69,70]. Due to the special structure of biomembranes made of phospholipids, it is important to understand their properties from the viewpoint of the mechanics [40,57,71] as well as from the viewpoint of the opening of the ion channels [48,72,73]. The next problem is the modelling of waves, which up to now has been focused mostly on single waves. The propagation of action potentials is widely studied [11][12][13]62]. The waves in biomembranes are described using the soliton model of Heimburg and Jackson [14] and their coworkers [38,74, etc.]. The full analysis of the governing equation (11) has revealed several types of solutions [58]. The possible pressure wave in the interstitial fluid is modelled theoretically [43,44] and measured by Terakawa [41].
The proposed simplified model (13)- (17) is an attempt to unite the possible components into a coupled set of equations of motion, resulting in the emergence of an ensemble of waves, which is a signature of complexity. Such ensembles of waves are known, for example, for the KdV solitons [75] and for the Schrödinger solitons [76]; the interaction of fields may also result in ensembles of waves like in thermoelasticity [77]. Note that the governing equations (13)-(15) are nonlinear and the coupling forces model the interaction. The possible ensemble of waves is schematically presented in Fig. 2.
The model satisfies directly assumptions (i)-(iv). The last assumption (v) stresses that the electrical signal has the main role in transmittance, and mechanical waves accompany the main signal.
From the viewpoint of modelling complex biological systems [24], it is obvious that the effects such as the energy exchange and spatio-temporal coupling between the single waves are taken into account, the ionic transport (the influence of the microstructure) is crucial to keep the action potential propagating, and nonlinearities play a decisive role. A threshold for triggering the process exists predicted by the HH [12] or by the FHN [35] models.

DISCUSSION
The proposed model (13)-(17) is a robust one and based on the causality of the general process generated by an electrical input (13) following the present understanding of axon physiology [62]. Quite certainly the model could be improved by using a more sophisticated HHtype model for the action potential and results of further studies, especially after the nature of coupling forces is better understood. The main question is related to the character of the coupling of an electrical signal with the deformation of the surrounding biomembrane, i.e. to the electro-mechanical coupling [48]. As the opening and closing of the ion channels are decisive for signal propagation, this question is closely related to the voltage sensitivity and/or mechanosensitivity of channel formation. Although the opening of ion channels may be caused besides the electrical signals also by mechanical signals [48], it might be possible to ignore this influence at the first step of model calibration (i.e. by taking b 1 = b 2 = 0). The force exerted on the biomembrane by an electrical signal could be related to the ion current through the pressure wave in the axoplasm. It is argued that this force is related to electrostriction [41,78] and can be taken as a function of the square of the voltage [44,48,79].
The propagation of the action potential is possible only because of the outflux and/or the influx of ions from the interstitial fluid. This influence is modelled either by adjusting nonlinear parameters a i , b i in Eq. (14) or by an additional external force F 3 (u, v) that acts on the initial system. Such a possibility was already shown by FitzHugh [13], who demonstrated the influence of the 'stimulus intensity' on the solution. In the case of the present model Eq. (14) this assumption means that a i + b i = const and an additional force appears on the right-hand side of Eq. (14). The present model describes the coupling of various components in an ensemble at the given moment, but the coupling could also be time-dependent taking into account the history of the process. There are several ways to model this complex process. El Hady and Machta [44] elaborated a mechanism of the coupling of electrical and mechanical signals based on the assumption that the potential energy of the process is stored in the surrounding biomembrane and kinetic energy in the axoplasmic fluid resulting in mechanical surface waves. The HH model is used for describing the action potential, and the force exerted on the biomembrane is taken proportional to the square of the voltage. The axoplasmic fluid is described by the linearized Navier-Stokes equation.
Another coupled model of electrical and mechanical signals based on the spring-dampers (dashpots) system is proposed by Jérusalem et al. [80]. The ion currents are again calculated using the HH model and calibrated for a guinea pig's spinal cord white matter. It is stressed that the strain and the strain rate have a decisive role linking electrophysiology and structural damage.
There are but a few experiments on simultaneous measurements of electrical and mechanical signals in nerve fibres. Terakawa [41] measured the propagation of the intracellular pressure in the axoplasm accompanying the action potential in squid giant axons. He showed that the pressure is correlated with the membrane potential but 'the pressure response was slower'. Experiments by Gonzales-Perez et al. [81] also demonstrated measured electrical and mechanical signals.
These signals, however, were shown not in real time but were only temporarily aligned. This indicates also the possibility of a delay in the deformation of the biomembrane due to the opening of ion channels. One should note also the measurements of coupled electrical signals and mechanical pulses in excitable plant cells (Chara braunii), which have shown synchronization of velocities of both pulses [82].
Our model, as presented in Section 4, is a system of differential equations, which actually means the assumption of continuous media for all structural constituents. There is a bulk of experimental values characterizing ion currents and membrane properties [37, 73,80, etc.] that can be used for evaluating the parameters of the model. As described above, two important problems are (i) quantification of the mechanisms of opening the ion channels and (ii) determination of the coupling forces.
As a first approximation the coupling forces should be related to the movement of ions through the membrane during the propagation of the nerve pulse. At the beginning, when the action potential forms, sodium ions move into the axon through the membrane increasing the pressure in the axoplasm locally, which can cause displacement of the lipid bilayer and result in the longitudinal density change. Similarly, when the potassium ions start moving out from the axon, the pressure in the axoplasm will change. The equilibrium is restored by the ion pumps, and another nerve pulse can be propagated. One can speculate that the density changes in the lipid bilayer could interfere with the activity of ion channels and their opening/closing time constants. This is just one possibility of constructing these coupling forces. The reality could be much more complex. For example, the local osmolarity might have a significant effect on the process and should be considered.
The calculations of single motions under an initial input, i.e.
solving Eqs (14)- (16) in their dimensionless form without coupling, demonstrate expected characteristics. The solution of the FHN equation (14) for the action potential AP demonstrates all the basic characteristics: the existence of a threshold, the formation of a stable asymmetric profile, and the existence of the refraction length [35]. The solution pressure wave of the 1D Navier-Stokes equation (15) is well known and can be compared with the flow in blood vessels [83,84, etc.]. The initial pulse-like excitation travels with changes due to viscosity and wall deformation if this is taken into account (not in our case).
The longitudinal mechanical wave in the biomembrane described by Eq. (16) can be found in the form of a solitary wave [14,58]. Note that the FHN and Heimburg-Jackson equations are actually two-wave equations [16] and from an initial input two waves are generatedone to the right and the other to the left. In terms of physiology this means that ortho-and antidromic waves can be generated from a localized excitation [79].
The question of synchronization of the velocities in a coupled model (cf. Fillafer et al. [82]) is of utmost importance. The velocity of action potential has been experimentally determined in many studies [39,85,86], and it depends significantly on the character of nerve fibres (diameter, myelin sheath, etc.). Early experiments established the velocity in squid axons in the range of 20 m s −1 but later the range was determined from 3 m s −1 to 100 m s −1 . The velocity of a mechanical pulse in a biomembrane is determined in the range of 170 m s −1 [14]. In plant cells both velocities are significantly lower but synchronized in the range of ca 10 m s −1 [82].
A numerical experiment by using the pseudospectral method [54,58] was carried out for solving Eqs (14) and (16), leaving aside the pressure wave and taking the coupling force as F 2 (z). Its functional shape is related to the changes of the ion current, which is coupled to the changes of voltage on the biomembrane. In the simple FHN model only one ion current is taken into account without specifying its character [35]. Note that in the HH model the sodium and potassium ions are specified with their complicated kinetics [39]. Here we assume that the mechanical wave is generated by F 2 (z) related to the gradient of the ion current. As far as the action potential (and the supporting ion current) is pulse-like in the moving frame, the corresponding bipolar gradient is energetically balanced. This is an essential property and means also satisfying the conservation laws. In numerical calculations the dimensionless variables Z and J in Eq. (14) and U in Eq. (16) are used, while n is the discrete dimensionless space (X = n i ∆X).The time and space scales match real values in large axons but as said above, the synchronous action of velocities needs further studies.
The generated ion current calculated from the FHN model (14) and its gradient are shown in Fig. 3. This bipolar gradient is used for the determining of the dimensionless coupling force F 2 (Z) = γ 1 J X , where γ 1 is a constant. If the coupling force is of a localized pulse type, the energetical balance will be distorted. At this stage we do not specify more the physical process of coupling but our assumption follows the discussion by Barz and Barz [67], El Hady and Machta [44], Gonzales-Perez et al. [81], and others. Note that using the gradient of the ion current for determining the driving force opens the door for possible modifications. Namely, the various ion currents like in the HH models can be differently influenced by various anaesthetic drugs [33].  Fig. 3. Normalized generated ion current J and its gradient J X calculated from the FHN Eq. (14). Here n is the discrete dimensionless space (X = n i ∆X).
The initial condition (13) is taken as a narrow sech 2 -type pulse with an amplitude above the threshold, which leads to a stable action potential. Figure 4 demonstrates a snapshot of the coupled process. The parameters in Eqs (14) and (16) are chosen such that lead the action potential and the mechanical wave in the biomembrane to propagate in the phase as an ensemble [44]. This corresponds to the situation in the experiments by Terakawa [41]. Note that this case corresponds to the anomalous dispersion of the mechanical wave, which explains the small vibration ahead of the main pulse [21]. The profile of the transverse wave can be determined from the calculated profile of the mechanical wave. Its functional bipolar shape is similar to the experimental results [41,55]. All amplitudes are normalized and profiles in Fig. 4 reflect the functional shape of the waves (AP, LW) and the accompanying ion current J in the ensemble. It is possible that if the velocity of the mechanical wave exceeds the velocity of the action potential, the mechanical wave in the biomembrane will affect the ion channels [48]. On the other hand, if the velocity of an action potential exceeds the velocity of the mechanical wave, the mechanical wave may influence the recovery process [41].
Numerical results demonstrate the possibility of building a coupled model for the complex process of signal propagation in nerve fibres. Even the robust model composed by Eqs (14) and (16) shows the propagation of an ensemble. Further studies are in progress in order to match the theoretical model with experimental data and to specify the character of the forcing terms together with the role of a pressure wave in the axoplasm. The possible synchronization of velocities of waves in the ensemble needs also further analysis. The earlier experiments [41,55] showed that the waves in the ensemble propagate in phase. However, the later estimations of velocities of single waves in nerve fibres differ considerably (cf. for example Hodgkin and Huxley [12] and Heimburg and Jackson [14]). The recent measurements of velocities in plant cells demonstrate clearly synchronization [82]. These contradicting results warrant further studies. Experimental verification of the model should follow studies of Terakawa [41], who measured all the components of a signal simultaneously. The regulation of the width of the longitudinal mechanical wave in a biomembrane permits the estimation of the coefficient h 2 as a measure of inertial properties of a lipid structure.
In conclusion, the signals in nerve fibres are composed of an ensemble of waves that are interacting with each other. The joint 'electromechanical' signal [14,52] exhibits clear signatures of complexity in dynamical processes [5]: nonlinearity and interaction between the components resulting in a unique biological unit [87]. Moreover, as a result of propagation, temporary changes occur in the structure of the biomembrane [14]. In addition, the process will be generated only if an initial input outvalues a critical threshold.
Many recent studies have tried to combine the physiological effects characteristic of signal propagation in nerve fibres into a more general model based on the explanation of possible interactive effects. In this paper, an attempt is made to build up a robust mathematical model based on coupled equations of motion as demonstrated in the analysis of many other physical problems [5]. To the best knowledge of the authors, this is the first attempt to describe the whole process in terms of differential equations. We certainly are aware that there are effects not considered in the proposed model such as heat production, the influence of the anaesthetic drugs, the role of proteins, etc. While mathematically clear and consistent (causality and reciprocity involved), the model can serve as a backbone inviting for further clarification and modification from the viewpoint of physiology. Characteristically of complex systems, the ensemble of waves resulting from the proposed model is more than single waves.