On solutions of a Boussinesqtype equation with amplitudedependent nonlinearities: the case of biomembranes
Abstract
Boussinesqtype wave equations involve nonlinearities and dispersion. In this paper a Boussinesqtype equation with amplitudedependent nonlinearities is presented. Such a model was proposed by Heimburg and Jackson (2005) for describing longitudinal waves in biomembranes and later improved by Engelbrecht et al. (2015) taking into account the microinertia of a biomembrane. The steady solution in the form of a solitary wave is derived and the influence of nonlinear and dispersive terms over a large range of possible sets of coefficients demonstrated. The solutions emerging from arbitrary initial inputs are found using the numerical simulation. The properties of emerging trains of solitary waves waves are analysed. Finally, the interaction of solitary waves which satisfy the governing equation is studied. The interaction process is not fully elastic and after several interactions radiation effects may be significant. This means that for the present case the solitary waves are not solitons in the strict mathematical sense. However, like in other cases known in solid mechanics, such solutions may be conditionally called solitons.
keywords:
nonlinearities, dispersion, solitary waves, solitons, emergence, interactions1 Introduction
The celebrated wave equation which is based on the conservation of momentum, models the motion with a finite speed. In order to account for accompanying physical phenomena, the wave equation must be modified. For conservative systems, the Boussinesqtype equations are widely used. The original Boussinesq equation was derived for surface waves on a fluid layer (Boussinesq, 1871; Rayleigh, 1876) but nowadays such equations are used also in solid mechanics (Christov et al., 2007). The main features of Boussinesqtype equations are: (i) bidirectionality (d’Alembert operator); (ii) nonlinearity of any order; (iii) dispersion of any order (presence of space and time derivatives of the fourth order or higher (Christov et al., 2007, etc). Beside fluid mechanics, there are many studies of such equations derived using various physical assumptions (Maugin, 1995, 1999; Porubov, 2003; Christov et al., 2007; Berezovski et al., 2013; Engelbrecht et al., 2015, etc). In solid mechanics, nonlinearity is caused by the nonlinear stressstrain relation and nonlinear strain tensor, i.e., physical and geometrical nonlinearities are involved (see, for example, Engelbrecht (1997)). The governing equations involve then type terms (), i.e., the displacement gradients enter the model. For example, the simple 1D equation reads
(1) 
where , , is the velocity in the unperturbed state and is the nonlinear parameter. Here and further, the indices and denote the differentiation with respect to the indicated variable. One could say that actually the velocity is calculated like
(2) 
The dispersive effects in solids are due to the geometry (Porubov, 2003) or due to the microstructure (Mindlin, 1964; Berezovski et al., 2013, etc). Then terms like , etc. appear in governing equations. The combined action of nonlinear and dispersive effects may give rise to solitary waves (Christov et al., 2007; Maugin, 2011; Engelbrecht et al., 2011, etc).
During the last decade the interest to mechanical waves in biomembranes has been growing fast (Heimburg and Jackson, 2005; Andersen et al., 2009; Appali et al., 2012, etc). The biomembranes have a special structure, made of lipids (Heimburg and Jackson, 2005; Mueller and Tyler, 2014) and in this case nonlinear effects are different from that in solids. Based on experimental results, the nonlinearity in biomembranes can be accounted in the velocity like (Heimburg and Jackson, 2005)
(3) 
where and are coefficients and is the density change along the axis of the biomembrane. This means that contrary to the gradienttype nonlinearity, the displacementtype nonlinearity appears in governing equations for waves in biomembranes. The HeimburgJackson model (Heimburg and Jackson, 2005), improved by Engelbrecht et al. (2015) takes such nonlinearities into account together with dispersive term(s). The governing equation is then of the Boussinesqtype and may lead to the emergence of solitary waves.
In this paper, the improved HeimburgJackson model (Engelbrecht et al., 2015) is systematically studied in detail needed for describing the possible emergence of solitary waves. After describing the derivation of the governing Boussinesqtype equation (Section 2), the following questions are analysed: (i) deriving the steady solutions to the governing equation (Section 3); (ii) finding the solutions for an arbitrary input (Section 4); (iii) studying the interaction of waves (Section 5). In this way, the existence of solitary waves is shown, the emergence of trains of solitary waves is demonstrated, and finally, the interaction of solitary waves shows whether the solitary waves are solitons in the classical sense. As it is well known, solitons interact with each other elastically without losses like elementary particles and only the phase shifts show the interaction effects (Zabusky and Kruskal, 1965; Drazin and Johnson, 1989; Salupere et al., 2002). In many physical systems the interaction is accompanied by radiation, i.e., the process is not fully elastic. In this case the solitary waves can only conditionally be called solitons. The final remarks are presented in Section 6 where the special features of solutions to this Boussinesqtype equation with displacementdependent nonlinearities are summarised. The analysis is wider than only the case of biomembranes and includes many combinations of governing parameters.
2 Derivation of the governing equation
The signal propagation in a nerve fibre is a complicated phenomenon. The nerve fibre itself can be modelled as a tube filled with axoplasm and surrounded by the extracellular fluid. The wall of the tube is made of a biomembrane (Debanne et al., 2011). The biomembrane is a very special biological structure made of phospholipids with hydrophobic tails directed to inside of the membrane, i.e., away from the intra and extracellular fluid (Mueller and Tyler, 2014). In general words, the lipid membrane represents a special biological microstructure with complicated properties. The concentration of ions within and outside of a fibre is different but the ion change can occur through the ion channels. These channels are closed at the rest but can be opened under electrical or mechanical impact (Mueller and Tyler, 2014).
The electrophysiological model describing the propagation of an electrical signal called the action potential was derived by Hodgkin and Huxley (1952) and is based on telegraph equations and on opening and closing the ion channels under the electrical impact. However, this model cannot explain all the complex effects in the nerve fibres. Experiments by Iwasa et al. (1980) and Tasaki (1988) have clearly demonstrated the swelling of the surrounding biomembrane and the accompanying heat exchange. This means that an action potential is accompanied also by a mechanical wave in the fibre wall. A mathematical model governing such a wave is proposed by Heimburg and Jackson (2005, 2007). Their model is based on the wave equation, i.e., on the balance of momentum and written in terms of density change in the longitudinal direction:
(4) 
Two essential assumptions are made. First, it is assumed that the velocity of a wave in a circular biomembrane is related to the compressibility of the lipid structure and can be taken as
(5) 
where and are coefficients and is the velocity of the small amplitude sound wave (Heimburg and Jackson, 2005).
The second assumption is to add ad hoc higher order term to the governing equation responsible for dispersion. The governing equation is then
(6) 
where is a constant. Equation (6) is a Boussinesqtype equation (see, for example, Christov et al., 2007). Heimburg and Jackson (2005) have demonstrated that Eq. (6) possesses a solitary pulsetype solution. There are several further studies analysing such solutions (Heimburg and Jackson, 2007; Andersen et al., 2009; Appali et al., 2012, etc). Equation (6) has been improved by Engelbrecht et al. (2015) in order to remove the discrepancy that at higher frequencies the velocities are unbounded. Following the ideas from the solid mechanics (Mindlin, 1964; Berezovski et al., 2013) and supported by the Lagrangian formalism, the inertial term is added to the governing equation:
(7) 
where and are dispersion coefficients.
The importance of the additional dispersion term can be explained with dispersion analysis. It has been shown (Engelbrecht et al., 2015) that in case of only one dispersion term (Eq. (6)), the phase velocity is expressed as and it tends to infinity as the wave number grows. In case of the second fourth order mixed dispersion term the propagation velocity is bounded as it can be seen in Fig. 1. The bounding velocity for high frequency harmonics is defined by the ratio of the dispersion coefficients () and the coefficient is related to the rate of change of the velocity from low frequency to the high frequency domain. Higher valued coefficient means that the transition from the low frequency speeds to the higher frequency speed is more rapid (see Fig. 1).
From the viewpoint of solid mechanics the importance of the fourth order mixed derivative is not surprising as it is well known that the presence of only spatial derivatives in the governing equation can lead to instabilities (Maugin, 1999). Moreover, the mixed fourth order derivative is related to the inertia of the microstructure and it is shown by Maurin and Spadoni (2016a) that both dispersive terms arise naturally when all effects involved in wave propagation in solids are considered and this has also been demonstrated experimentally (Maurin and Spadoni, 2016b).
The focal point of this paper is the full analysis of Eq. (7). Further it is convenient to use the dimensionless form of Eq. (7), which will take the form
(8) 
where , , and , . Here is a certain length, for example, the fibre diameter.
Equation (8) must be solved under initial and boundary conditions formulated in the dependent variable .
3 Steady solutions
In this section we focus our analysis on undistorted travelling waves in the form
(9) 
where is some function and is dimensionless wave velocity (Ablowitz, 2011; Drazin and Johnson, 1989). Substituting this into Eq. (8) we get
(10) 
Integrating Eq. (10) twice we get after some rearranging
(11) 
where and are constants of integration. Since we are looking for solitary wave solutions, then we may add boundary conditions that as and therefore (Ablowitz, 2011; Drazin and Johnson, 1989). Now the Eq. (11) is multiplied by and integrated to get
(12) 
which can be rewritten as
(13) 
where
(14) 
is a fourthorder ‘pseudopotential’. Note that for the classical KdV equation the ‘pseudopotential’ is of the third order (Drazin and Johnson, 1989).
The existence of solitary waves can be analysed by either investigating the behaviour of the ‘pseudopotential’ (14) or the phase portrait of Eq. (12). In case of the ‘pseudopotential’ (14) also applies for the HeimburgJackson model (6) and has been analysed by Lautrup et al. (2011) for a particular set of parameters that were determined experimentally and are relevant for the solitary wave propagation in biomembranes (, ). Here the analysis is more general and the signs of the parameters and are not fixed.
The four zeros of the polynomial (14) are
(15) 
Double zero at indicates the saddle point, which is minimal requirement for the existence of solitary waves (Ablowitz, 2011; Drazin and Johnson, 1989). The following analysis can be divided into two parts: the cases of , which has also been analysed previously (Peets et al., 2016) and . Attention is paid to the signs of and which govern the structure of solutions.
(i)
The case of . For this case the analysis is pretty straightforward. It can be deduced from aforementioned restrictions and from Eq. (15) that in this case the additional condition for the velocity is
(16) 
which means that the in case of and the solitary waves governed by the Eq. (8) will always travel slower than the low frequency sound. This is in good agreement with the actual pulse propagation in biomembranes (Heimburg and Jackson, 2005; Lautrup et al., 2011).
The ‘pseudopotential’ (14) and the phase portrait (12) for this case have been plotted in Fig. 2 for (left column) and for (right column), respectively. The ‘pseudopotential’ (14) has been plotted in the top row and the phase portraits for the case in the middle and for the case is plotted in the bottom row for reference.
The existence of solitary wave solution requires that has a local minimum at with at least one local maximum next to it (Figs 2a,b). Alternatively one can study the phase portrait (Figs 2c,d): solitary wave solutions exist when a saddle point and a homoclinic orbit exists (shown in blue). The amplitude of the solitary wave in both cases is determined by . It is clear that while the magnitude of the amplitude of a solitary wave depends on the ratio of the parameters and together with the velocity , the sign of the amplitude is determined only by the parameter : in case of positive solitary wave emerges and in case of the amplitude will be negative. It can also be shown that higher values of result in lower amplitudes meaning that the lower amplitude solitary waves travel faster as it has been shown earlier (Peets and Tamm, 2015; Tamm and Peets, 2015).
The case of is shown in Fig. 3, where it can be seen that the behaviour of the ‘pseudopotential’ and the phase portrait is significantly different from the case of , . Although solitary wave solutions are possible in the region where , only amplitude is realised, because the local maximum between and is closer to the saddle point (cf. Dauxois and Peyrard, 2006). Also, in this case the speed of the solitary wave can have any value between zero and meaning that also in this case the solitary wave travels slower than the speed of the low frequency sound. As in case of , here also the magnitude and the sign of the amplitude is determined by the parameters , and and the higher velocities result in lower amplitudes.
Recalling that the analytical solution of Eq. (8) is (Peets et al., 2015, 2016)
(17) 
where and is the velocity of the solitary wave, the solitary wave solutions for the given cases are plotted in Fig. 4 where the solid lines represent the case of and the dashed represents the case of , which are the original HeimburgJackson equation (6).
Although solitary waves in case of and only exist when , periodic solutions to Eq. (12) also exist when
(18) 
This case is demonstrated in Fig. 5 where it can be seen that in this case ‘pseudopotential’ (14) is positive between the points and and a stable orbit exists (shown in blue) which means an existence of a periodic solution. What is interesting is that the phaseportrait in this case looks similar to Fig. 3c only it has been slightly shifted to the right and the higher amplitude part is realised.
It can also be seen in Figs 2, 3 and 4 that in the case of the second dispersion coefficient more localised solution is obtained: the greater value of the quantity means the steeper slope (and hence the smaller width) of the solitary wave. The effect of the dispersive term on the width of a solitary wave is demonstrated in Fig. 6, where it can be seen that higher values of result in more localised solutions.
(ii)
Since Eq. (12) can be thought of as conservation of ‘pseudoenergy’, then if then also has to be negative. In case of the condition means that periodic solutions emerge even in the case of as it is seen in Fig. 7. Here also the periodic solution oscillates between the values and (shown in blue in Fig. 7), which corresponds to the region where as it can be seen in Fig. 2a. Similar result is obtained when , only with negative amplitude.
In case of and the ‘pseudopotential’ will only have regions and in case of a solitary wave with amplitude exist (Fig. 8a,c). In case of and , a solitary wave with amplitude exists (Fig. 8b,d). Like in previous cases the structure of the phase portrait depends on the sign of the coefficient and the sign of (amplitude of the solution) depends on the sign of the coefficient . Unlike in case of where smaller amplitude solitary waves travel faster, in case of the higher amplitude waves travel faster. Also recall that if then periodic solution exists with the same coefficients(see Fig. 5).
Last case we mention is , and when no solutions exist.
4 Solutions emerging from arbitrary initial conditions
In the present paper we use the pseudospectral method (PSM) to solve the governing equation (8) under localised initial conditions demonstrating the influence of parameters and on the evolution of solutions. The PSM is a well established method which is used for solving PDE’s and ODE’s on regular basis. The advantages and disadvantages of the PSM are well explored in the literature (Fornberg and Sloan, 1994; Fornberg, 1998). Here two points are worth of highlighting: (i) the PSM requires one to use periodic boundary conditions, (ii) the governing equations have to be in a suitable form for applying the PSM with time derivatives on the left hand side and spatial derivatives on the right hand side of the equation. The first point is not a problem, however, taking a look at Eq. (8) it is evident that we have a mixed partial derivative term . We use a change of variables for transforming the governing equation (8) for allowing the application of the PSM (Engelbrecht et al., 2015; Peets and Tamm, 2015; Tamm and Peets, 2015). The basic idea of the PSM is to find the spatial derivatives by making use of the properties of the Fourier transform and then solve the resulting ODE with respect to time derivate by making use of the commonly available schemes for numerical solving of the ODE’s.
For initial and boundary conditions for the systematic analysis we use a pulse–type localised initial condition in the form of type profile: where , i.e., the total length of the spatial period is . Here the amplitude and the width of the initial pulse are and . The initial phase velocity is meaning that the initial condition splits into two pulses propagating in the opposite directions. Some examples are provided using differentt combination of parameters in which case the used parameters are noted separately.
The calculations are carried out with the Python package SciPy (Jones et al., 2007) with Python interface to the ODEPACK FORTRAN code (Hindmarsh, 1983) for the ODE solver.
The typical solution of Eq. (8) can be seen in Fig. 9. Under the arbitrary initial condition the initial disturbance splits into two equal wave structures propagating in opposing directions (initial velocity was zero). In the balanced dispersion case the emerging waveprofiles are of the solitary wave type and very stable even through multiple interactions.
In addition to the formation of solitary waves a number of different waveprofile regimes exist for the solutions of the governing equations (8) depending on the parameters but also on the initial conditions. To name the ones investigated previously:
(i) solitary waves (single or as a part of solitary wave train, see Figs 10,11);
(ii) Airy or reverse Airy like oscillatory structures (see Fig. 10);
(iii) hybrid solution where part of the initial pulse evolves into a train of solitary waves and remainder of the initial pulse forms an oscillatory structure (Peets et al., 2015; Peets and Tamm, 2015; Tamm and Peets, 2015).
From the viewpoint of nerve pulse propagation the most interesting one is the solitary wave solution, however, the rest of the solution types can not be ignored either as these might be relevant somehow for either nerve pulse propagation or some kind of pathologies. It should be emphasised that not only the equation parameters are important in determining what kind of solution evolves from the initial excitation but also the character of said initial excitation is important. As an example see Fig. 10 where some solutions corresponding to the different parameters and initial condition amplitudes are presented. Depending on the dispersion type the initial excitation sign determines if the emergining wave structure is composed of solitary pulses or Airy or reverseAiry type oscillatory packet under the parameter combinaion used in Fig. 10. In essence, a parameter combination which in Fig. 10 leads to solitary wave train would be of Airytype if the initial excitation is with opposite sign and vice versa, which previously was of Airytype solution would be a train of solitary waves. Another interesting phenomenon which must be mentioned is the case where smaller amplitude solitary waves can travel faster than the high amplitude ones as can be seen in Fig. 11. Under the suitable parameter combinations it is possible to observe a solutions where both negative and positive amplitude solitary waves can exist simultaneously and if the smaller amplitude solitary waves travel faster then the larger negative amplitude solitary waves travel even faster.
However, this is not an universal symmetry but needs the right ratio of parameters. The most common solution type seems to be the oscillatory structure with few solitary pulses where some part of the initial pulse energy is sufficient to form one or more solitary pulses and the remainder forms an oscillatory trail either in front or behind (depending of the dispersion type) of the propagating solitary waves. As the system is conservative and we have periodic boundary conditions then the resulting wave profiles keep interacting until after sufficiently long evolution the radiation from not fully elastic interactions causes the solitary waves to merge into the oscillatory structure where they can no longer be detected separately. For all practical purposes the integration time needed for that to happen is so long that before this scenario becomes relevant one would have to question if the decision to not take dissipation into account is still relevant from the viewpoint of physical interpretation of the results.
In Figs. 12, 13 and 14 one can see selected example countour plots with isoline interval of from to for the amplitude.
In addition we are tracking waveprofile peak trajectories by finding the exact local maxima of the wave profiles by making use of the properties of the Fourier transform (Salupere, 2009) (reconstructing the wave profile from the Fourier spectrum to minimise inaccuracies from using the discrete grid) for finding the exact spatial coordinates of the pulse peaks at each time step for the purposes of finding the waveprofile velocities. Following observations can be highlighted from Figs. 12, 13 and 14:
(i) the dispersion parameters have a strong effect on the evolution of the wave profiles – the main pulse velocities are clearly different depending on the dispersion parameters and in addition the dispersion type determines on which side (relative to the propagation direction) the secondary wave structures emerge from the main pulse. Increasing the parameter increases the main pulse propagation velocity as predicted by dispersion analysis (Peets and Tamm, 2015);
(ii) in the balanced dispersion case the waveprofile remains stable through several interactions even under relatively strong nonlinear parameters (see Fig. 14) in addition, it should be emphasised that this is not a classical dispersionless case as both dispersive terms are, in fact, nonzero but are just balanced against each other resulting in a situation where for majority of the frequencies the group and phase speed are equal to each other.
(iii) The nonlinear parameters and have some influence on the waveprofile propagation velocities (see, for example, Fig. 12 vs. Fig. 14).
In the case of the normal dispersion increasing the nonlinearity leads to the slower propagation velocity for the wave profiles. In the case of the anomalous dispersion the main pulse velocity remains almost the same, however, where the effect is more prominent is the secondary oscillatory structures where in Fig. 14 the secondary structure starts to separate at approximately while in Fig. 12 the secondary pulses start to separate at approximately meaning that in the case of higher nonlinear parameters the secondary structures have propagated at higher velocity. This is in agreement with previous results where it has been demonstrated that due to the uncommon (in the context of Boussinesq type equation) nonlinear terms there can exist parameter combinations where the smaller amplitude solitary waves propagate faster than the higher amplitude ones (Tamm and Peets, 2015).
Next, let us take a more detailed look how the the nonlinear and dispersive parameters influence the observable quantities of the wave profiles under the positive and negative initial conditions. We observe the speed of the peak of the main pulse. We track the coordinates of the peak of the main pulse by reconstructing the waveprofile shape from the full Fourier spectrum at each time step. Parameters and change from to with the step size of and for dispersion related parameters and three combinations are recorded – an normal dispersion case (), balanced dispersion case () and anomalous dispersion case ().
(i) , :
Normal dispersion case. The negative amplitude initial condition results in the faster propagation velocity of the pulse than the positive amplitude initial condition. Increasing leads to increase of the observed main pulse velocity and increasing of the main pulse amplitude. On the other hand, increasing the nonlinear parameter towards zero leads to decrease in main pulse velocity for the negative amplitude initial condition. The positive initial amplitude case has velocity gains for the main pulse. As far as the amplitude goes, if the initial condition is positive then the main pulse amplitude increases and if negative then decreases. The oscillatory structures amplitudes remain roughly the same when nonlinear parameter changes.
Balanced dispersion case. For the negative initial amplitude case the main pulse velocity is greater than in the case of the positive initial amplitude. Increasing the parameter leads to increasing main pulse velocity and the main pulse amplitude. Increasing the parameter towards zero leads to the main pulse velocity moving closer to one (the normalised speed of sound in the present context). If the initial condition is with negative amplitude then increasing parameter results in the reduction of velocity for the main pulse and if the initial amplitude is positive then the main pulse propagation velocity is increased. The main pulse amplitude is close to half of that of the initial condition and the oscillatory structures are practically nonexisting.
Anomalous dispersion case. In the anomalous dispersion case the negative amplitude initial condition results in the faster main pulse velocity in the case of positive initial amplitude than in the case if the initial amplitude is negative. What is different in this case is that increasing leads to decrease of the observed main pulse velocity in the case of negative initial condition while in the case of positive amplitude initial condition this leads to increase of the observed main pulse velocity. Increasing the parameter leads to small increase of the main pulse velocity under both used initial condition signs.
(ii) , :
Normal dispersion case. The negative amplitude initial condition leads to a greater main pulse velocity than the positive amplitude initial condition. Under the both initial condition signs decreasing the nonlinear parameter (towards the zero) leads to a small decrease of the main pulse velocity. In the case of the negative amplitude initial condition the main pulse amplitude is greater than in the case of the positive amplitude initial condition and the observed oscillations are larger for the case with positive initial amplitude than in the case with the negative initial amplitude. Increasing parameter leads to decrease in the main pulse velocity in the case of the negative amplitude initial condition while in the case of the positive initial amplitude the main pulse velocity remains the same. Increasing parameter towards zero leads to marginally greater amplitude for the main pulse in the case of negative amplitude initial condition while in the case of positive initial amplitude the main pulse amplitude is unaffected by the changes in the nonlinear parameter . The oscillatory structure magnitude is unaffected in the normal dispersion case.
Balanced dispersion case. The main pulses propagate with the velocity close to one for the both considered initial condition signs. Decreasing parameter leads to a decrease in the observed main pulse velocity and main pulse amplitude under both of the initial condition signs. Increasing parameter towards the zero leads to main pulse speeds closer to one under both considered initial condition signs. It should be mentioned that the amplitudes of the main pulses are greater than half of the initial pulse height which is due to the nonlinear effects which are combined with relatively weak dispersion under the used parameters combination.
Anomalous dispersion case. The main pulses propagate with velocity greater than one under both of the considered initial condition signs, however, the main pulse amplitudes and associated oscillatory structures are different. Increasing the parameter leaves the observed propagation speed the same but decreases the observed main pulse amplitude and leaves the observed oscillatory structures about the same. Increasing the parameter does not affect the main pulse velocity significantly in the considered dispersion case regardless of the sign of the initial amplitude. However, increasing the nonlinear parameter decreases the main pulse amplitude and increases the amplitude of the oscillatory structures under the both considered initial condition signs.
(iii) , :
Normal dispersion case. The negative amplitude initial condition leads to a smaller main pulse velocity than the initial condition with positive amplitude. Regardless of the velocity difference the amplitudes of the pulses are comparable and the same can be observed for the oscillatory tails. Increasing the parameter does not affect noticeably the solution corresponding to the negative initial amplitude while in the case of positive initial amplitude the main pulse velocity and amplitude increase with increasing the parameter towards the zero. Increasing the parameter leads to a small decrease for the main pulse velocity for the negative initial amplitude case and to a increase of the main pulse velocity in the case of positive amplitude initial condition. In addition, increasing nonlinear parameter decreases the main pulse amplitude in the case of negative initial condition and in the case of positive amplitude initial condition the amplitude of the main pulse is increased when parameter increases.
Balanced dispersion case. The main pulses tend to propagate at velocity close to one and maintain amplitude which is close to the half of that of the initial condition. Increasing the parameter towards zero leads, in general, to the increase of the main pulse velocity and amplitude. Increasing the parameter leads to decrease of the main pulse velocity and amplitude in the case of the negative initial condition amplitude and to the increase of the main pulse velocity and decrease of the amplitude in the case of positive amplitude initial condition.
Anomalous dispersion case. The main pulse propagation velocity is greater than one and the influence of the nonlinear parameters is small as far as the main pulse evolution is concerned. Increasing the parameter leaves the main pulse velocity the same but reduces the main pulse amplitude by a small amount in the case of the negative initial amplitude. If we have the positive initial amplitude then increasing the parameter leads to a reduction of the velocity of the main pulse. Increasing the parameter leads to a decrease in the main pulse propagation velocity and amplitude under both initial condition signs.
(iv) , :
Normal dispersion case. For the considered nonlinear parameters signs the decreasing the parameter leads in the case of negative initial amplitude to the decrease of the main pulse and oscillatory structure amplitude while the velocity remains practically the same. In the case of the positive initial condition amplitude there is small decrease in the main pulse velocity if parameter decreases while the amplitude and oscillatory structure remain practically the same for the amplitude of the main pulse and for the oscillatory structure. Decreasing the nonlinear parameter leads to increase of velocity and decrease of the main pulse amplitude and the oscillatory structure amplitude in the case of negative initial condition amplitude. In the case of positive initial amplitude decreasing leads to the decrease of the main pulse velocity and amplitude and the increase of the oscillatory structure amplitude.
Balanced dispersion case. In the balanced dispersion case the main pulse propagation velocities are close to one and the main pulse amplitudes close to the half of the initial condition amplitude. Decreasing parameter in the case of negative initial condition amplitude leads to decrease of the main pulse velocity and the amplitude of the oscillatory structure. In the case of positive initial condition amplitude decreasing parameter leads to marginally smaller main pulse amplitude and marginally larger oscillatory structure amplitude with minor drop also in the main pulse propagation velocity. Decreasing the parameter leads in the case of negative initial condition amplitude to a greater propagation velocity of the main pulse as well as to a increased main pulse amplitude and decrease of the oscillatory structure amplitude close to zero. Decreasing the parameter in the case of positive initial amplitude leads to a decrease of the main pulse velocity while the amplitude of the pulse is increased and the oscillatory structure amplitude is suppressed.
Anomalous dispersion case. Decreasing the parameter leads to a marginal decrease of the main pulse and oscillatory structure amplitude if the initial condition is with negative amplitude and to a increase of the main pulse velocity and amplitude and to a decrease of the oscillatory structure amplitude if the initial condition has a positive amplitude. Decreasing the parameter leads to a small decrease in the main pulse velocity for the amplitude and increase of the oscillatory structure amplitude if the initial condition is with negative amplitude. If the initial condition has a positive amplitude then decreasing the parameter leads to a increase in the main pulse velocity and amplitude and to a decrease of the oscillatory structure amplitude.
5 Interaction of solitons
In general, a soliton can be described as a stable particlelike state of a nonlinear system (Dodd et al., 1982). Another way of describing the phenomenon we call soliton is through its properties. A soliton is a wave in the nonlinear environment that (1) has a stable form, (2) is localized in space and (3) restores its speed and structure after interaction with another soliton (Drazin and Johnson, 1989; Engelbrecht, 1995). Solitons emerge when there is a balance in the system between dispersive and nonlinear effects. In essence it can be said that solitons are nonlinear waves that behave between interactions like linear waves. A solitary wave is usually a wave in the nonlinear environment where all the key properties of solitons are not strictly fulfilled. For example, if the interaction between two waves is not entirely elastic (or it is not possible to observe the interaction) or if the form of the wave is not sufficiently stable in time, then the wave is often called a solitary wave to distinguish it from the soliton.
In Fig. 15 one can see the HJ model (6) solitary wave prpagation (top) and interaction (bottom). The parameters are the same as in Section 3 except . From Fig. 15 it is clear that while the single HJ pulse is stable it is a solitary wave, not a soliton, because the interaction with another such wave is not elastic as there is significant radiation even during the first interaction event and the shape of the waveprofile is not properly restored after the interaction. However, it should be noted that the parameter combinations can exist where the solitary wave solutions can be relatively stable with almost no radiation.
The interaction of single solitary waves depends on the parameters of he model and consequently, on the dispersion type (normal, anomalous). We start with the following set , , , , , , and . In Fig. 16 one can see the interactions if the parameter and in the bottom if – the interactions are remarkably similar and non disruptive with the main difference being that the solitonic waveprofiles are more localized if . In this case the interactions have almos no radiation (negligible radiation two orders of magnitude smaller than the main pulse amplitude at ). Amplitude isolines are separated by from 0.001 to 0.013 in Fig. 16.
In Fig. 17 waveprofiles and corresponding phase plots after the five interaction events () are presented. The solitonic waveprofiles corresponding to higher values of are more localised as expected (the waveprofile in Fig. 17 is propagating to the left). The small distortions to the waveprofiles are easier to spot in phase plots (right), in particular the small radiation close to zero which is two orders of magnitude smaller, as noted, than the main pulse under the used parameter combination.
Let us return to a parameter set presented in Section 3 for the analytical solution. It turns out that there is also a scenario possible where the solitonic solutions with the additional dispersive term are more stable through interactions than the solitonic solutions if parameter . In Fig. 18 the case is presented at the top and the case in the bottom. The amplitude isolines are separated by 0.02 from 0.02 to 0.3. It is clear that under the parameter set used in Figs 15 and 18 the solitonic waves corresponding to have greater amount of radiation than the case which is relatively stable in comparison throughout interactions. Neither of the cases can be considered solitons in the strict mathematical case (Drazin and Johnson, 1989) as in both cases there is significant enough radiation after only three interactions. While unrelated to the mechanics of soliton interactions it is interesting to remark that the used numerical algorithm performs approximately three times faster if .
In addition it should be noted the the parameter set used for the systematic analysis in Section 4 has nonneglible amounts of radiation during the interactions as well. However, this is not easy to see in Figs. 12, 13 and 14 at the given scale. For that parameter set the amplitude loss due to the interaction events is less than 10% over dozen interaction events.
6 Final remarks
The systematic analysis of solutions to the special Boussinesqtype equation with the displacementdependent nonlinearities has revealed several interesting phenomena. The analysis is focused on Eq. (7) (or its dimensionless form (8)) which is the improved HeimburgJackson model for describing the longitudinal wave process in biomembranes. Like every wave equation it describes the process generated by initial and/or boundary conditions expressed in terms of the dependent variable. Here the variable under consideration is the change of the density in the longitudinal direction. In terms of this variable the existence of solitary solutions is demonstrated, the emergence of trains of solitary pulses is shown and the properties of emergence analysed, and the interaction of single solitary waves and trains studied. The governing nonlinear wave equation is actually a novel mathematical model compared to the conventional models in continuum theory where the nonlinearities as a rule are of the deformationdependent.
The analysis can be summarised with following conclusions:

The additional dispersive term with the coefficient (or in the dimensionless form) in addition to the ad hoc dispersive term (Heimburg and Jackson, 2005) describes actually the influence of the inertiality of the microconstituents (lipids) of the biomembrane. This corresponds to the understandings of continuum mechanics of microstructured solids (Mindlin, 1964) and is demonstrated also experimentally (Maurin, 2015). This term regulates the width of the solitary pulse (see Fig. 6) and such an effect can be used for determining the value of from experiments. It also determines how fast the transition from the low frequency speeds to the high frequency speeds occurs (see Fig. 1);

Soliton trains can be emerged from an arbitrary initial condition. These results were obtained by numerical simulation by using the pseudospectral method (Salupere, 2009). Depending on the signs of coefficients and , the nonlinear effects start to influence the emergence either from the front or from the back of the propagating pulse (see Fig 11). For the case of a biomembrane one has , and the train emerging from a positive input starts with smaller solitons which travel faster than the bigger ones. This is different from the conventional case of nonlinear evolution equations (the KdV equation, for example). In the case of a negative input, the train is headed by bigger solitons which travel faster (see Fig. 11). It has been shown that there are several wave types possible: solitary waves (Fig. 10), oscillatory (Airytype) waves (Fig. 10), and hybrid solutions.

The interaction of solitary waves is not fully elastic (see Figs 16, 18) which shows that these solitary waves are not solitons in the strict sense (Drazin and Johnson, 1989). However, like in other Boussinesqtype equations (Christov et al., 2007; Engelbrecht et al., 2011), the radiation effects accompanying every interaction start cumulating rather slowly and the interacting solitons keep their shape for a rather long time. It gives the ground to call emerging solitary waves modelled by Eq. (7) (or Eq. (8)) solitons like it is done in other physical cases (Maugin, 2011).
Biological structures as a rule have high complexity because the macrobehaviour is strongly influenced by the embedded microbehaviour. Mathematical modelling is a tool not only for describing biological processes and performing experiments in silico. The behaviour of biomembrane is an excellent example how the microstructure (lipids) of a membrane has a direct impact on wave phenomena along the membrane. The analysis of the governing equation (7) (or Eq. (8)) presented above demonstrates the richness of the model from the viewpoint of mathematical physics and opens the ways for physiological experiments concerning the properties of biomembranes.
Acknowledgements
This research was supported by the European Union through the European Regional Development Fund (Estonian Programme TK 124) and by the Estonian Research Council (projects IUT 3324, PUT 434).
References
Footnotes
 journal: XYZ
References
 Ablowitz, M. J., 2011. Nonlinear Dispersive Waves. Asymptotic Analysis and Solitons. Cambridge Univ Press, Cambridge.
 Andersen, S. S. L., Jackson, A. D., Heimburg, T., 2009. Towards a thermodynamic theory of nerve pulse propagation. Prog. Neurobiol. 88 (2), 104–13.
 Appali, R., Van Rienen, U., Heimburg, T., 2012. A comparison of the HodgkinHuxley model and the soliton theory for the action potential in nerves. In: Iglic, A. (Ed.), Adv. Planar Lipid Bilayers Liposomes, Vol. 16. Academic Press, Ch. 9, pp. 275–299.
 Berezovski, A., Engelbrecht, J., Salupere, A., Tamm, K., Peets, T., Berezovski, M., 2013. Dispersive waves in microstructured solids. Int. J. Solids Struct. 50 (1112), 1981–1990.
 Boussinesq, J., 1871. Théorie de l’intumescence liquide, applelée onde solitaire ou de translation, se propageant dans un canal rectangulaire. Comptes Rendus l’Academie des Sci. 72, 755–759.
 Christov, C. I., Maugin, G. A., Porubov, A. V., 2007. On Boussinesq’s paradigm in nonlinear wave propagation. Comptes Rendus Mécanique 335 (910), 521–535.
 Dauxois, T., Peyrard, M., 2006. Physics of Solitons. Cambridge University Press, Cambridge.
 Debanne, D., Campanac, E., Bialowas, A., Carlier, E., Alcaraz, G., 2011. Axon physiology. Physiol. Rev. 91 (2), 555–602.
 Dodd, R., Eilbeck, J., Gibbon, J., Morris, H., 1982. Solitons and nonlinear wave equations. Academic Press Inc. LTD., London.
 Drazin, P., Johnson, R., 1989. Solitons: and Introduction. Cambridge University Press, Cambridge.
 Engelbrecht, J., 1995. Beautiful dynamics. Proc. Estonian Acad. Sci. Physics/Mathematics 1 (44), 108–119.
 Engelbrecht, J., 1997. Nonlinear Wave Dynamics. Complexity and Simplicity. Kluwer, Dordrecht.
 Engelbrecht, J., Salupere, A., Tamm, K., 2011. Waves in microstructured solids and the Boussinesq paradigm. Wave Motion 48 (8), 717–726.
 Engelbrecht, J., Tamm, K., Peets, T., 2015. On mathematical modelling of solitary pulses in cylindrical biomembranes. Biomech. Model. Mechanobiol. 14, 159–167.
 Fornberg, B., 1998. A Practical Guide to Pseudospectral Methods. Cambridge University Press, Cambridge.
 Fornberg, B., Sloan, D., 1994. A review of pseudospectral methods for solving partial differential equations. Acta Numer. 3, 203–267.
 Heimburg, T., Jackson, A., 2007. On the action potential as a propagating density pulse and the role of anesthetics. Biophys. Rev. Lett. 2, 57–78.
 Heimburg, T., Jackson, A. D., 2005. On soliton propagation in biomembranes and nerves. Proc. Natl. Acad. Sci. U. S. A. 102 (28), 9790–5.
 Hindmarsh, A., 1983. ODEPACK, a Systematized Collection of ODE Solvers. Vol. 1. NorthHolland, Amsterdam.
 Hodgkin, A. L., Huxley, A. F., 1952. A quantitative description of membrane current and its application to conduction and excitation in nerve. J. Physiol. 117 (4), 500–544.
 Iwasa, K., Tasaki, I., Gibbons, R., 1980. Swelling of nerve fibers associated with action potentials. Science 210 (4467), 338–339.
 Jones, E., Oliphant, T., Peterson, P., 2007. SciPy: open source scientific tools for Python.
 Lautrup, B., Appali, R., Jackson, A. D., Heimburg, T., 2011. The stability of solitons in biomembranes and nerves. Eur. Phys. J. E. Soft Matter 34 (6), 1–9.
 Maugin, G. A., 1995. On some generalizations of Boussinesq and KdV systems. Proc. Estonian Acad. Sci. Phys. Math. 44, 40–55.
 Maugin, G. A., 1999. Nonlinear Waves in Elastic Crystals. Oxford University Press, Oxford.
 Maugin, G. A., 2011. Solitons in elastic solids (19382010). Mech. Res. Commun. 38 (5), 341–349.
 Maurin, F., 2015. Wave propagation in periodic buckled beams. Ph.D. thesis, Ècole Polytechnique Fédérale de Lausanne.
 Maurin, F., Spadoni, A., 2016a. Wave propagation in periodic buckled beams. Part I: Analytical models and numerical simulations. Wave Motion.
 Maurin, F., Spadoni, A., may 2016b. Wave propagation in periodic buckled beams. Part II: Experiments. Wave Motion.
 Mindlin, R. D., 1964. Microstructure in linear elasticity. Arch. Ration. Mech. Anal. 16 (1), 51–78.
 Mueller, J. K., Tyler, W. J., 2014. A quantitative overview of biophysical forces impinging on neural function. Phys. Biol. 11 (5), 051001.
 Peets, T., Tamm, K., 2015. On mechanical aspects of nerve pulse propagation and the Boussinesq paradigm. Proc. Estonian Acad. Sci. 64 (3S), 331–337.
 Peets, T., Tamm, K., Engelbrecht, J., 2015. Numerical investigation of mechanical waves in biomembranes. In: Elgeti, S., Simon, J.W. (Eds.), Conf. Proc. YIC GACM 2015 3rd ECCOMAS Young Investig. Conf. 6th GACM Colloquium, July 2023, 2015, Aachen, Ger. pp. 1–4.
 Peets, T., Tamm, K., Engelbrecht, J., 2016. On the role of nonlinearities in the Boussinesqtype wave equations. Wave Motion (1), 1–7.
 Porubov, A. V., 2003. Amplification of Nonlinear Strain Waves in Solids. World Scientific, Singapore.
 Rayleigh, L., 1876. On waves. Philos. Mag. 1 (1), 257–279.
 Salupere, A., 2009. The pseudospectral method and discrete spectral analysis. In: Quak, E., Soomere, T. (Eds.), Appl. Wave Math. Springer Berlin Heidelberg, Berlin, pp. 301—334.
 Salupere, A., Peterson, P., Engelbrecht, J., 2002. Longtime behaviour of soliton ensembles . Part I ââ Emergence of ensembles. Chaos, Solitons & Fractals 14, 1413–1424.
 Tamm, K., Peets, T., 2015. On solitary waves in case of amplitudedependent nonlinearity. Chaos, Solitons & Fractals 73, 108–114.
 Tasaki, I., 1988. A macromolecular approach to excitation phenomena: mechanical and thermal changes in nerve during excitation. Physiol. Chem. Phys. Med. NMR 20, 251–268.
 Zabusky, N., Kruskal, M., 1965. Interaction of “solitons” in a collisionless plasma and the recurrence of initial states. Phys. Rev. Lett. 15 (6), 240–243.