Microtubules are the major part of the cytoskeleton. They are involved in nuclear and cell division and serve as a network for motor proteins. The first model that describes nonlinear dynamics of microtubules was introduced in 1993. Three nonlinear models are described in this chapter. They are longitudinal U-model, representing an improved version of the first model, radial φ -model and new general model. Also, two mathematical procedures are explained. These are continuum and semi-discrete approximations. Continuum approximation yields to either kink-type or bell-type solitons, while semi-discrete one predicts localized modulated waves moving along microtubules. Some possible improvements and suggestions for future research are discussed.
- partial and ordinary differential equations
- kink solitons
A cell is defined as eukaryotic if it has a membrane-bound nucleus. Such cells are generally larger and much more sophisticated than prokaryotic ones. Microtubules (MTs) are the basic components of cytoskeleton existing in eukaryotes . They are long structures that spread between a nucleus and a cell membrane. MTs play an essential role in the shaping and the maintenance of cells and are involved in cell division. Also, they represent a network for motor proteins. These proteins move with a velocity of  carrying a certain cargo such as mitochondrion.
All eukaryotic cells produce two kinds of tubulin proteins. These are and tubulins, or monomers, and they spontaneously arrange head to tail forming biologically functional subunit that we call a heterodimer, or a dimer for short. When intracellular conditions favor assembly, the dimers assemble into long structures called protofilaments (PFs). Microtubules are usually formed of 13 PFs, as shown in Figure 1.
Hence, MTs are long cylindrical polymers whose lengths vary from a few hundred nanometers up to meters in long nerve axons . Each dimer is an electric dipole whose mass and length are and , respectively . The component of its electric dipole moment in the direction of PF is . Consequently, MT as a whole appears to be a giant dipole with negatively charged end coinciding with biologically positive end (more active) and vice versa. This is the reason why an intrinsic electric field exists within MT.
MTs in non-neuronal cells are unstable structures. They exhibit dynamic instability behavior existing in phases of elongation (polymerization) or rapid shortening (depolymerization). This size fluctuation has been called as dynamic instability [7, 8]. Notice that the shrinkage rate is bigger than the growth rate (see Ref.  and references therein). MTs grow steadily at positive end, corresponding to the – subunit, and then shrink rapidly by loss of tubulin dimers at the negative end, corresponding to the –monomer. Many anticancer drugs, for example, taxol (paclitaxel), prevent growth and shrinkage of MTs and thus prevent cell proliferation .
MTs existing in neuronal cells are stable and, consequently, neurons, once formed, do not divide . This stability is crucial as there are evidences that neuronal MTs are responsible for processing, storage and transduction of biological information in a brain [4, 11].
It was mentioned that MTs represent the traffic road for motor proteins. Some more information can be found in Ref.  and in an exhaustive review paper . It suffices now to state that the cellular motors with dimensions of less than 100 nm convert chemical energy into useful work. These small machines have the fundamental role of dissipation in biological systems, which has been confirmed by both the theoretical and the experimental investigations . The molecular motors dissipate continuously and operate as irreversible systems .
It is clear that any molecular motor, to start moving, should obtain a certain signal. One of the promising dynamical mechanisms for intracellular signaling is solitary waves, which is explained in this chapter.
2. Mechanical models
MTs, as well as all biological systems, are nonlinear in nature. Strong covalent chemical bonds are usually modeled by linear “springs”, while weak chemical interactions, existing in all biological systems, are modeled by nonlinear “springs”. This means that expressions for energy of biological systems require nonlinear terms, which brings about nonlinear partial differential equations (PDEs) explaining nonlinear dynamics of these systems. This is the topic of the present chapter. We will see that, in case of MTs, the solutions of these nonlinear PDEs are solitary waves.
The word soliton was introduced in 1965 to designate solitary waves describing the propagation of excitations in continuous media with nonlinearity and dispersion . The first qualitative description of solitary waves dates back to 1834 when hydrodynamic engineer John Scott Russell observed them on a surface in a shallow channel . The wave was so stable that the engineer followed it about 1 or 2 miles. From then, there has been tremendous interest for various kinds of solitons in many branches of physics [15, 16, 17, 18, 19]. In this chapter, the terms soliton and solitary waves are treated as synonyms, which is commonly accepted in literature.
Solitons are localized waves possessing some interesting properties. The most important is their stability in a sense that they conserve their shape and energy after mutual interaction. In other words, they can pass through one another without annihilation. This was experimentally observed in neurons .
To model complex MT dynamics, we should introduce some simplifications. To the best of the author’s knowledge, all the models introduced so far have only one degree of freedom per dimer. Hence, for the models explained in this chapter, elementary subunits of PFs are dimers and they perform either longitudinal or angular oscillations and the appropriate models can be called as longitudinal or angular (radial), respectively.
The longitudinal contacts along PFs are much stronger than those between adjacent PFs [21, 22], which allows us to construct a simplified Hamiltonian of MT, which is, practically, Hamiltonian for a single PF only. However, the influence of the neighboring PFs is taken into consideration through the electric field. Namely, each dimer exists in the electric field coming from the dimers belonging to all PFs. Also, the nearest neighbor approximation is assumed.
The first model that describes nonlinear dynamics of MTs is a longitudinal one. It was introduced in 1993 by Satarić et al. . According to the model, the dimers perform angular oscillations but a coordinate , describing the dimer’s displacement, is a projection of the top of the dimer on the direction of PF. Therefore, the displacements are radial but the used coordinate is longitudinal. There is a real longitudinal model assuming longitudinal displacements of the dimers that we call as Z-model . Both U- and Z-models bring about equal crucial differential equations and the latter one will not be studied here.
Somewhat improved and more general version of the first nonlinear model is what we call as U-model  and this will be explained in the following paragraphs. Both models are based on the fact mentioned above that the dimers are electric dipoles and that the whole MT can be regarded as ferroelectric [23, 26], which means that the interaction between a single dimer and its surrounding can be modeled by W-potential [23, 27]. This yields to the following Hamiltonian for MT [23, 25, 28]
where dot means the first derivative with respect to time while the integer determines the position of the considered dimer in PF. The first term obviously represents a kinetic energy of the dimer of mass . The second one is interaction between the neighboring dimers belonging to the same PF in the nearest neighboring approximation and is an intradimer stiffness parameter. The next two terms represent the W-potential energy mentioned earlier, where the parameters and should be determined or, at least, estimated and are assumed to be positive. We should point out that the double-well potential is rather common in physics [27, 29, 30]. The very last term is coming from the fact that the dimer is the electric dipole existing in the field of all other dimers, where Q > 0 represents the excess charge within the dipole and E > 0 is internal electric field. The last three terms together can be regarded as unsymmetrical W-potential.
Our final goal is the function , describing nonlinear dynamics of MT. This function is a solution of so-called dynamical equation of motion, which can be obtained from Eq. (1). To derive it, we introduce generalized coordinates and defined as and . Using well-known Hamilton’s equations of motion and , we obtain the following discrete differential equation that should be solved
The last term is a viscosity force with being a viscosity coefficient . Therefore, nonlinear dynamics of MTs has been described by Eq. (2). Obviously, nonlinearity is coming from the fourth degree term in the W-potential.
It was explained earlier that we used some approximations to derive Eq. (2). However, we need one more to solve it. We now explain two mathematical methods for solving this equation. Practically, these two approaches are two approximations. They are continuum and semi-discrete approximations. We will see that the different mathematical procedures yield to different solutions. Therefore, the function depends not only on the physical system but also on the used mathematical method.
Let us explain the continuum approximation first. A question if MTs are discrete or continuum systems was studied in Ref. , where it was shown that the continuum approximation is valid. The continuum approximation means a transition , which allows a series expansion of the terms , that is, , where is the dimer’s length explained earlier. In fact, PF can be seen as one-dimensional crystal with being a period of the lattice. This straightforwardly brings about the following continuum dynamical equation of motion
This is PDE that cannot be easily solved. Hopefully, this equation can be transformed into an ordinary differential equation (ODE). It is well known that, for a given wave equation, a traveling wave is a solution which depends upon and only through a unified variable as , where and are constants. If we substitute the variables and by we straightforwardly transform Eq. (3) into the following ODE
Eq. (4) becomes the appropriate one in Ref.  for . Therefore, the U-model is more general than its predecessor introduced in Ref. . It is crucial that the parameter can be determined together with the function for known or estimated and . This is because has very important physical meaning. The first term in Eq. (3) is the inertial term and it is coming from the kinetic energy in Hamiltonian (1), while the second one is the elastic one. Therefore, positive means that the inertial term is bigger than the elastic one and vice versa.
Eq. (4) has already been solved using different mathematical procedures like standard procedure [23, 27, 29, 30] and method of factorization [31, 32]. There exists a group of procedures where the function is represented as a serious expansion over other known function like . The function is usually known and we plug into Eq. (4) and determine the coefficients . A common example for is a solution of Riccati equation, which is either tangent or tangent hyperbolic. As only the latter function may have physical meaning, we call the method as tangent hyperbolic function method (THFM) [25, 33, 34, 35] and extended or modified extended THFM . The function can also be one of Jacobian elliptic functions  and, even, unknown .
It is very likely that the most general procedure is the simplest equation method (SEM) [39, 40, 41] and its simplified version called as modified simplest equation method (MSEM) . According to SEM, the series expansion is [39, 40, 41].
where , and are coefficients that should be determined and represents the first derivative. In general, the function is known and represents a solution of a certain ODE of lower order than the equation that should be solved. A commonly used example is the Riccati equation 
To determine the positive integer in Eq. (6), we should plug , const, into Eq. (4) and concentrate our attention on the leading terms . One can easily show that for Eq. (4) as the leading terms are proportional to and .
In what follows, we assume .
Obviously, this is satisfied if all the coefficients are simultaneously equal to zero. This brings about a system of seven equations, which can be obtained using Mathematica or similar software . One of them can be written as
indicating two possible relationships between the parameters and . Hence, there are a few cases to be studied. They are as follows : (1) , ; (2) , ; (3) ; (4) , ; (5) , and (6) , .
It is obvious that the first case represents nothing but a simpler method called extended tanh-function method. The system mentioned earlier brings about 
The final result is 
where is the following three real solutions of the first of Eqs. (11)
Of course, these three solutions exist for . The case was discussed in Ref. .
All the three solutions are shown in Figure 2 for and . Of course, these solutions reproduce previously known results . Figure 2 shows that the solutions of Eq. (4) are kink and antikink solitons. More detailed analysis of their physical meaning is given in Ref. .
It was shown  that the second case is equal to the first one indicating that the value of is irrelevant if . In other words, we could have assumed a simpler version of the Riccati equation neglecting the term .
The third case is more interesting. It turns out that, instead of the three lines in Figure 2, that is, the three solutions, we obtain infinitely many lines corresponding to each of them . However, they represent three groups of parallel lines, which means that all these solutions are only shifted functions and, consequently, have equal physical meaning. Therefore, this case does not bring about any physically new result.
The final expression for is
This case yields to a new solution, which was not obtained using less general mathematical methods. However, it may be interesting from mathematical point of view only as diverges for .
Case 5 is a simplified version of SEM, explained in Ref. . The mentioned system brings about as well as
where a notation has been introduced to distinguish this parameter from , used in the previous cases. It is interesting to compare the polynomials for and , existing in Eqs. (11) and (17). We can see that.
The function is shown in Figure 3 for and for two values of the parameter . We notice very interesting result that is a bell-type soliton! This certainly demonstrates the advantage of SEM method over the less general ones.
It is important to study the physical meanings of the parameters and . Eq. (19) indicates that solitonic width is inversely proportional to and that does not affect maximum of the wave. Figure 3 shows that the amplitude of is a decreasing function with respect to .
Finally, the last case gives the solution
which is obviously divergent for , .
Therefore, all the cases are explained and we can see that the continuum approximation yields to both kink solitons and bell-type solitons. The latter may exist only if viscosity is neglected.
It was stated earlier that the coordinate was the projection of the top of the dimer on the direction of MT. A patient reader may ask how can be negative when this is the projection. This question is answered in Ref. .
It was mentioned earlier that there are two approximations that can be used to solve Eq. (2). Now we get back to Eq. (2) and study semi-discrete one [15, 28, 43]. A mathematical basis for the method is a multiple-scale method or a derivative-expansion method [16, 44]. We assume small oscillations
which straightforwardly transforms Eq. (2) into
where ω is the optical frequency of the linear approximation, q = 2π/λ > 0 is the wave number, cc represents complex conjugate terms and the function F0 is real. Of course, is the dimer’s length, as mentioned earlier. The function is continuous and represents an envelope, while exp(iθn), including discreteness, is a carrier component. Notice that the parameter ε exists in the function , but does not in exp(iθn). This is because the frequency of the carrier wave is much higher than the frequency of the envelope and we need two time scales, and , for those two functions. The same holds for the coordinate scales.
To simplify the problem, a continuum limit should be introduced as well as new transformations and . This allows a series expansion of , that is
where indexes and denote the first and the second derivative with respect to . Hence, the function becomes
This crucial expression represents a starting point for a series of important expressions. They can be obtained equating the coefficients for the various harmonics. For example, equating the coefficients for and neglecting all the terms with and one obtains the following expressions for the dispersion relation and the group velocity :
Also, the coefficients for and , respectively give 
which yields to
Eqs. (28) and (29) and new coordinates and , defined as and , allows us to simplify Eq. (27). An explanation for why the parameter exists in the time scaling but is absent in the space scaling is given in Refs. [45, 46]. If we consider the terms for again we obtain the well-known nonlinear Schrödinger equation (NLSE) for the function
where the dispersion coefficient and the coefficient of nonlinearity are
where parameters and represent envelope and carrier component velocities, while the amplitude and the soliton width have the forms
It is very difficult to deal with the parameters and as is completely unprecise statement. However, seems to be more practical. Hence, new parameters and have been introduced as , and . Finally, we can easily obtain the expression for the longitudinal displacement of the dimer at the position
One more parameter can be eliminated using the idea of coherent mode . This mode means that the envelope and the carrier wave velocities are equal. It follows from Eq. (35) that , which yields to the function . This means that the wave is the one phase function, preserving its shape in time.
To plot the function or, equivalently, we should know or estimate the values of a couple of the parameters. Of course, if 2D plot is chosen, can be presented as either a function of at a certain position or as a function of for chosen . Very detailed analysis of the parameter selection was done in Ref. . One example for is shown in Figure 4. Obviously, this is a localized modulated wave usually called as breather. We can see that its width is about 200 nm, which means that it covers about 25 dimers.
As a conclusion, we can state that the two mathematical procedures bring about even three results, that is, three different solitons. These are kinks, bell-type solitons and breathers. They may be signals for the motor proteins to start moving, as explained in Introduction.
Obviously, viscosity has been neglected. This will be explained in the following section, within the -model.
A weak point of the U-model is the last term in Eq. (1). A scalar product would be better choice for the potential energy, where is the distance between the centers of positive and negative charges within the dipole. This potential indicates the angle as a coordinate instead of the projection and the Hamiltonian for the radial model, which we call as -model, is [50, 51]
where is a moment of inertia of the dimer at the position . Notice that the W-potential does not exist in Eq. (38) even though the terms including and appear as a result of a series expansion of the cosine function. Instead of viscosity force introduced in the previous section, we introduce viscosity momentum , where is the viscosity coefficient [51, 52, 53]. Following the procedure explained earlier, we obtain
Like above, the solutions are kink solitons .
where is the soliton velocity, while and are corresponding sound velocities. According to Eq. (11), we can see that the U-model predicts as is positive. On the other hand,  means that, according to the -model, the kink belongs to the class of supersonic solitons. We will return to this issue in the next section.
Now, we switch to the semi-discrete approximation within the -model to solve the dynamical equation of motion, which is 
The final solution is the same as except that and are different. Therefore, both the U- and the -models predict the breather waves moving through MT.
5. General model of MTs
It was mentioned earlier that the weak point of the U-model is the last term in Eq. (1). Also, it is better to use the radial coordinate than the longitudinal one as we assume angular oscillations of the dimers. The scalar product , existing in the -model, solved these problems but the W-potential has been missing. In fact, a series expansion of gives and terms but with opposite signs from those in the U-model. These two terms are, practically, a potential that looks like W in a mirror having only one minimum surrounded by two maxima and, due to its shape, can be called as M-potential . This potential brings about , which is disputable result.
Therefore, we want to solve the mentioned problem regarding the U-model but to keep the W-potential, the coordinate and, probably, , that is, . This suggests the following Hamiltonian
where , and has the same meaning as in the -model. Let us call the model as general one (GM). The procedure mentioned earlier brings about
where, of course, . If we consider Eqs. (3) and (47), we can see that the last two terms in Eq. (47) may be the first derivatives of either W- or M-potential, depending on the sign of the terms in the brackets. However, these brackets may have different signs or can be zero. Therefore, the possible cases are:
Case 1: , Case 2: ,
Case 3: , , Case 4: , .
All of them are studied in Ref.  and they will be explained here briefly.
Case 1 straightforwardly yields to
and the final solution is 
One of the advantages of the GM over the -model is the value of amplitude. Namely, the amplitude of the kink soliton, according to the -model, is , coming from Eq. (40). This is unrealistic, too big value. Instead of , the appropriate factor, existing in the GM, is K, given by Eq. (49).
If viscosity is neglected, the GM brings about
Case 2 straightforwardly yields to
and to the final results
It is obvious that these results do not have physical meaning as is complex, while may diverge.
Case 3 brings about
as well as , which certainly means that Eq. (55) does not have any solution having physical sense.
The remaining Case 4 linearizes Eq. (47) and will not be studied here.
Therefore, the GM yields to the kink solitons as the previous two models do. However, this is the radial model and the problems with both the last term in Eq. (1) and the huge amplitude in the case of the -model have been solved. We should study one more issue. It was mentioned earlier that the U-model predicts the subsonic kink soliton, while the -model predicts the supersonic wave. How about the GM? It was shown that Case 1 yields to the solutions having physical sense and that . According to Eq. (49), we easily reach the final conclusion:
If and , then and the function is subsonic soliton, kink for the positive in Eq. (50) and antikink otherwise.
If and , then and the function is supersonic soliton, antikink for the positive in Eq. (50) and kink otherwise.
All this certainly suggests the advantages of the GM with respect to the previous two.
6. Conclusion and future research
In this chapter, the three models describing nonlinear dynamics of MTs are shown. The first one, the U-model, is the improved version of the first nonlinear model and it predicts subsonic kink solitons moving along MT. The second one, the radial -model, predicts the supersonic kinks. Finally, the GM is explained. This is the radial model which yields both possibilities regarding the kink’s speed. If we assume that the kink soliton is subsonic wave then we know the minimum value of the parameter , that is, , as explained earlier.
Two mathematical procedures are explained, continuum and semi-discrete approximations. It is very interesting that the final result depends not only on the physical system but on the mathematical methods as well. These solutions are the kink soliton and the breather. The question is which one, if any, really moves along MTs. This is not known in the moment and cannot be without experimental results.
It was demonstrated that the GM is better than the previous two models. However, this does not mean that it should not be improved. For example, there has been an attempt to improve the model introducing Morse potential instead of the harmonic one . The harmonic potential energy assumes that attractive and repulsive forces are equal. Morse potential is not symmetric and is good for both strong and weak interactions.
In this chapter, the dimers are considered as elementary units. However, their structure is more complicated and they include tubulin tales (TTs). Consequently, nonlinear dynamics of TTs should also be studied and some results already exist [56, 57].
The W-potential has two minima which means that it assumes existence of the two angles between the dimer and the direction of PF around which the dimer oscillates. One of the future tasks should be measuring these angles. First of all, such experiment would check if the W-potential is correct or not. If it is, then our knowledge of their values would improve the theory a lot.
One of the future research goals should be two-component model. This may mean that we should construct the model assuming two degrees of freedom. However, one of these degrees can be an internal one, which means that oscillations if monomers within the dimer should be taken into consideration. Notice that the two-component model may be the one studying electro-acoustic wave excitations .
This work was supported by funds from Serbian Ministry of Education, Sciences and Technological Development (grant No. III45010).