Mechanical Models of Microtubules

Abstract

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.

1. Introduction

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 [1]. 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 $0.1-2\mathrm{\mu m}/\mathrm{s}$ [2] carrying a certain cargo such as mitochondrion.

All eukaryotic cells produce two kinds of tubulin proteins. These are $\alpha$ and $\beta$ 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.

Figure 1.

A tubulin dimer, a protofilament and a microtubule [3].

Hence, MTs are long cylindrical polymers whose lengths vary from a few hundred nanometers up to meters in long nerve axons [4]. Each dimer is an electric dipole whose mass and length are $m=1.8\times {10}^{-22} \mathrm{kg}$ and $l=8 \mathrm{nm}$ , respectively [5]. The component of its electric dipole moment in the direction of PF is $p=337\text{Debye}=1.13\times {10}^{‐27} \mathrm{Cm}$ [6]. 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. [9] and references therein). MTs grow steadily at positive end, corresponding to the $\beta$ – subunit, and then shrink rapidly by loss of tubulin dimers at the negative end, corresponding to the $\alpha$ –monomer. Many anticancer drugs, for example, taxol (paclitaxel), prevent growth and shrinkage of MTs and thus prevent cell proliferation [10].

MTs existing in neuronal cells are stable and, consequently, neurons, once formed, do not divide [4]. 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. [9] and in an exhaustive review paper [12]. 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 [13]. The molecular motors dissipate continuously and operate as irreversible systems [13].

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 [14]. 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 [15]. 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 [20].

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.

3. U-model

The first model that describes nonlinear dynamics of MTs is a longitudinal one. It was introduced in 1993 by Satarić et al. [23]. According to the model, the dimers perform angular oscillations but a coordinate $u$ , 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 [24]. 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 [25] 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]

$$H_u=\sum _n\left[\frac{m}{2}{{\dot{u}}_n}^2+\frac{k_u}{2}{\left(u_{n+1}-u_n\right)}^2-\frac{1}{2}{\mathit{Au}}_n^2+\frac{1}{4}{\mathit{Bu}}_n^4-{\mathit{QEu}}_n\right]$$ E1

where dot means the first derivative with respect to time while the integer $n$ determines the position of the considered dimer in PF. The first term obviously represents a kinetic energy of the dimer of mass $m$ . The second one is interaction between the neighboring dimers belonging to the same PF in the nearest neighboring approximation and $k$ is an intradimer stiffness parameter. The next two terms represent the W-potential energy mentioned earlier, where the parameters $A$ and $B$ 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 $u_n\left(t\right)$ , 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 $q_n$ and $p_n$ defined as $q_n=u_n$ and $p_n=m\mathrm{d}{u}_n/\mathrm{d}t$ . Using well-known Hamilton’s equations of motion $\mathrm{d}{p}_n/\mathrm{d}t=-\mathrm{d}H/\mathrm{d}{q}_{\mathrm{n}}$ and $\mathrm{d}{q}_n/\mathrm{d}t=\mathrm{d}H/\mathrm{d}{p}_{\mathrm{n}}$ _, we obtain the following discrete differential equation that should be solved

$$m{\ddot{u}}_n=k_u\left(u_{n+1}+u_{n-1}-2u_n\right)+{\mathit{Au}}_n-{{\mathit{Bu}}_n}^3+\mathit{QE}-\gamma {\dot{u}}_n$$ E2

The last term is a viscosity force with $\gamma$ being a viscosity coefficient [23]. 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 $u_n\left(t\right)$ 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. [31], where it was shown that the continuum approximation is valid. The continuum approximation means a transition $u_n\left(t\right)\to u\left(x,t\right)$ , which allows a series expansion of the terms $u_{n\pm 1}$ , that is, $u_{n\pm 1}\to u\pm \frac{\partial u}{\partial x}l+\frac{1}{2}\frac{{\partial }^{2}u}{\partial x^2}{l}^2$ , where $l$ is the dimer’s length explained earlier. In fact, PF can be seen as one-dimensional crystal with $l$ being a period of the lattice. This straightforwardly brings about the following continuum dynamical equation of motion

$$m\frac{{\partial }^{2}u}{\partial t^2}-k_u{l}^2\frac{{\partial }^{2}u}{\partial x^2}-\mathit{QE}-\mathit{Au}+{\mathit{Bu}}^3+\gamma \frac{\partial u}{\partial t}=0$$ E3

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 $u\left(\xi \right)$ is a solution which depends upon $x$ and $t$ only through a unified variable $\xi$ as $\xi =\mathit{\kappa x}-\mathit{\omega t}$ , where $\kappa$ and $\omega$ are constants. If we substitute the variables $x$ and $t$ by $\xi$ we straightforwardly transform Eq. (3) into the following ODE

$${\alpha }_u{\Psi }^{″}-{\rho }_u{\Psi }^{\prime }-\Psi +{\Psi }^3-\sigma =0$$ E4

where $u^{\prime }\equiv \mathrm{d}u/\mathrm{d}\xi$ and

$$u=\sqrt{\frac{A}{B}}\Psi ,{\alpha }_u=\frac{m{\omega }^2-k_u{l}^2{\kappa }^2}{A},\sigma =\frac{\mathit{qE}}{A\sqrt{A/B}},{\rho }_u=\frac{\mathit{\gamma \omega }}{A}$$ E5

Eq. (4) becomes the appropriate one in Ref. [23] for ${\alpha }_u=-1$ . Therefore, the U-model is more general than its predecessor introduced in Ref. [23]. It is crucial that the parameter ${\alpha }_u$ can be determined together with the function $\Psi$ for known or estimated $\sigma$ and ${\rho }_u$ . This is because ${\alpha }_u$ 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 ${\alpha }_u$ 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 $\Psi$ is represented as a serious expansion over other known function like $\Psi =\sum _{k=0}^N{A}_k{\mathrm{\Phi }}^k$ . The function $\mathrm{\Phi }$ is usually known and we plug $\Psi$ into Eq. (4) and determine the coefficients $A_k$ . A common example for $\mathrm{\Phi }$ 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 [36]. The function $\mathrm{\Phi }$ can also be one of Jacobian elliptic functions [37] and, even, unknown [38].

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) [42]. According to SEM, the series expansion is [39, 40, 41].

$$\Psi =A_0+\sum _{k=1}^N\left(A_k{\mathrm{\Phi }}^k+B_k{\left(\frac{{\mathrm{\Phi }}^{\prime }}{\mathrm{\Phi }}\right)}^k\right)$$ E6

where $A_0$ , $A_k$ and $B_k$ are coefficients that should be determined and ${\mathrm{\Phi }}^{\prime }$ represents the first derivative. In general, the function $\mathrm{\Phi }=\mathrm{\Phi }\left(\xi \right)$ 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 [40]

$${\mathrm{\Phi }}^{\prime }+{\mathrm{\Phi }}^2-2a\mathrm{\Phi }-b=0,a,b=\text{const}$$ E7

To determine the positive integer $N$ in Eq. (6), we should plug $\Psi =c/{\xi }^p$ , $c=$ const, into Eq. (4) and concentrate our attention on the leading terms [42]. One can easily show that $N=1$ for Eq. (4) as the leading terms are proportional to ${\xi }^{-\left(p+2\right)}$ and ${\xi }^{-3p}$ .

The well-known general solution of Eq. (7) is [39, 40]

$$\mathrm{\Phi }=a+\sqrt{a^2+b}\text{tanh}\left[\sqrt{a^2+b}\left(\xi -{\xi }_0\right)\right]$$ E8

In what follows, we assume ${\xi }_0=0$ .

Our next step is determination of the parameters $A_0$ , $A_1$ , $B_1$ , $a$ , $b$ and ${\alpha }_u$ . According to Eqs. (6) for $N=1$ and (7), we obtain the expressions for ${\Psi }^{\prime }$ , ${\Psi }^{″}$ and ${\Psi }^3$ as required by Eq. (4), which yields to the following expression:

$$K_3{\mathrm{\Phi }}^3+{K_3}^{\text{'}}{\mathrm{\Phi }}^{-3}+K_2{\mathrm{\Phi }}^2+{K_2}^{\text{'}}{\mathrm{\Phi }}^{-2}+K_1\mathrm{\Phi }+{K_1}^{\text{'}}{\mathrm{\Phi }}^{-1}+K_0=0$$ E9

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 [39]. One of them can be written as

$$K_3\equiv \left(A_1-B_1\right)\times \left[2{\alpha }_u+{\left(A_1-B_1\right)}^2\right]=0$$ E10

indicating two possible relationships between the parameters $A_1$ and $B_1$ . Hence, there are a few cases to be studied. They are as follows [39]: (1) $B_1=0$ , $a=0$ ; (2) $B_1=0$ , $a\ne 0$ ; (3) $A_1=B_1$ ; (4) $2{\alpha }_u=-{\left(A_1-B_1\right)}^2$ , $A_1{B}_1\ne 0$ ; (5) $A_1=0$ , $a\ne 0$ and (6) $A_1=0$ , $a=0$ .

It is obvious that the first case represents nothing but a simpler method called extended tanh-function method. The system mentioned earlier brings about [39]

$$8{A_0}^3-2A_0+\sigma =0,{\alpha }_u=-\frac{{A_1}^2}{2},A_1=-\frac{{\rho }_u}{3A_0},b=\frac{1-3{A_0}^2}{{A_1}^2}$$ E11

The final result is [39]

$${\Psi }_i=A_{0i}+A_{1i}{\mathrm{\Phi }}_i,{\mathrm{\Phi }}_i=\sqrt{b_i}\text{tanh}\left(\sqrt{b_i}\xi \right)$$ E12

where $A_{0i}$ is the following three real solutions of the first of Eqs. (11)

$$A_{01}=\frac{1}{2\sqrt{3}}\left(\text{cos}F+\sqrt{3}\text{sin}F\right),A_{02}=\frac{1}{2\sqrt{3}}\left(\text{cos}F-\sqrt{3}\text{sin}F\right)$$ E13

$$A_{03}=-\frac{1}{\sqrt{3}}\text{cos}F,F=\frac{1}{3}\text{arccos}\left(\frac{\sigma }{{\sigma }_0}\right),{\sigma }_0=\frac{2}{3\sqrt{3}}$$ E14

Of course, these three solutions exist for $\sigma <{\sigma }_0$ . The case $\sigma >{\sigma }_0$ was discussed in Ref. [25].

All the three solutions are shown in Figure 2 for $\sigma \approx 0.9{\sigma }_0$ and ${\rho }_u=1$ . Of course, these solutions reproduce previously known results [25]. 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. [25].

Figure 2.

The functions $\Psi \left(\xi \right)$ for ${\rho }_u=1$ and $\sigma =0.34$ .

It was shown [42] that the second case is equal to the first one indicating that the value of $a$ is irrelevant if $B_1=0$ . In other words, we could have assumed a simpler version of the Riccati equation neglecting the term $2a\mathrm{\Phi }$ .

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 [39]. 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.

Case 4 is suggested by Eq. (10). The system of seven equations, mentioned earlier, gives the first and the last term in Eq. (11) as well as

$$a=0,A_1=2B_1,{\alpha }_u=-\frac{{B_1}^2}{2},B_1=-\frac{{\rho }_u}{3A_0}$$ E15

The final expression for $\Psi$ is

$$\Psi \left(\xi \right)=A_0-\sqrt{1-3{A_0}^2}\left(\text{tanh}\hspace{0.17em}y+\frac{1}{\text{sinh}y}\right),y=\frac{3A_0}{2{\rho }_u}\sqrt{1-3{A_0}^2}\xi$$ E16

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 $\Psi$ diverges for $\xi =0$ .

Case 5 is a simplified version of SEM, explained in Ref. [42]. The mentioned system brings about ${\rho }_u=0$ as well as

$$B_1=-\frac{{\mathrm{A}}_0}{a},{\alpha }_u=-\frac{{B_1}^2}{2},{{\mathrm{A}}_0}^3-{\mathrm{A}}_0-\sigma =0,b=\frac{a^2\left({{\mathrm{A}}_0}^2-1\right)}{2{{\mathrm{A}}_0}^2}$$ E17

where a notation ${\mathrm{A}}_0$ has been introduced to distinguish this parameter from $A_0$ , used in the previous cases. It is interesting to compare the polynomials for $A_0$ and ${\mathrm{A}}_0$ , existing in Eqs. (11) and (17). We can see that.

$${\mathrm{A}}_{0i}=-2A_{0i},i=1,2,3$$ E18

which means that the values for ${\mathrm{A}}_{0i}$ are given by Eqs. (13), (14) and (18). We can easily show that the final solution for $\Psi$ is [39]

$$\Psi ={\mathrm{A}}_0-\frac{{\mathrm{A}}_0{K}^2}{{\text{cosh}}^2\left(\mathit{aK\xi }\right)\left[1+K\text{tanh}\left(\mathit{aK\xi }\right)\right]},K=\frac{\sqrt{3{{\mathrm{A}}_0}^2-1}}{{\mathrm{A}}_0\sqrt{2}}$$ E19

Obviously, this function cannot diverge for any value of $\mathit{a\xi }$ but only for $-1<K<1$ . Also, $K$ should be real and these two requirements eliminate ${\Psi }_2$ and ${\Psi }_3$ [39], which means that $\Psi$ and ${\mathrm{A}}_0$ in Eq. (19) are ${\Psi }_1$ and ${\mathrm{A}}_{01}$ .

The function ${\Psi }_1\left(\xi \right)$ is shown in Figure 3 for $a=0.1$ and for two values of the parameter $\sigma$ . We notice very interesting result that is a bell-type soliton! This certainly demonstrates the advantage of SEM method over the less general ones.

Figure 3.

A bell-type soliton for a = 0.1 and σ = 0.34 (a) and σ = 0.1 (b).

It is important to study the physical meanings of the parameters $a$ and $\sigma$ . Eq. (19) indicates that solitonic width is inversely proportional to $\left|a\right|$ and that $a$ does not affect maximum of the wave. Figure 3 shows that the amplitude of ${\Psi }_1$ is a decreasing function with respect to $\sigma$ .

Finally, the last case gives the solution

$$\Psi =\pm \frac{\sqrt{2}}{\text{sin}\left(2\sqrt{-b}\xi \right)},b<0$$ E20

which is obviously divergent for $\sqrt{-b}\xi =\mathit{k\pi }$ , $k=0,\pm 1,\pm 2,\dots$ .

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 $u$ was the projection of the top of the dimer on the direction of MT. A patient reader may ask how $u$ can be negative when this is the projection. This question is answered in Ref. [28].

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

$$u_n\left(t\right)=\epsilon \hspace{0.17em}{\mathrm{\Phi }}_n\left(t\right),\epsilon <<1$$ E21

which straightforwardly transforms Eq. (2) into

$$\epsilon m{\ddot{\mathrm{\Phi }}}_n=\epsilon k\left({\mathrm{\Phi }}_{n+1}+{\mathrm{\Phi }}_{n-1}-2{\mathrm{\Phi }}_n\right)+\epsilon A{\mathrm{\Phi }}_n-{\epsilon }^{3}B{{\mathrm{\Phi }}_n}^3+\mathit{qE}+\mathrm{O}\left({\epsilon }^4\right)$$ E22

According to the semi-discrete approximation, we look for wave solution which is a modulated wave, that is [28, 45]

$${\mathrm{\Phi }}_n\left(t\right)=F\left(\xi \right)e^{i{\theta }_n}+\epsilon F_0\left(\xi \right)+\epsilon F_2\left(\xi \right)e^{i2{\theta }_n}+\mathrm{cc}+\mathrm{O}\left({\epsilon }^2\right)$$ E23

$$\xi =\left(\epsilon nl,\epsilon t\right),{\theta }_n=nql-\omega t$$ E24

where ω is the optical frequency of the linear approximation, q = 2π/λ > 0 is the wave number, cc represents complex conjugate terms and the function F₀ is real. Of course, $l$ is the dimer’s length, as mentioned earlier. The function $F$ is continuous and represents an envelope, while exp(iθ_n), including discreteness, is a carrier component. Notice that the parameter ε exists in the function $F$ , 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, $t$ and $\epsilon t$ , for those two functions. The same holds for the coordinate scales.

To simplify the problem, a continuum limit $\mathit{nl}\to z$ should be introduced as well as new transformations $Z=\epsilon z$ and $T=\epsilon t$ . This allows a series expansion of $F\left(\xi \right)$ , that is

$$F\left(\epsilon \left(n\pm 1\right)l,\epsilon t\right)\to F\left(Z,T\right)\pm F_Z\left(Z,T\right)\epsilon l+\frac{1}{2}{F}_{\mathit{ZZ}}\left(Z,T\right){\epsilon }^2{l}^2$$ E25

where indexes $Z$ and $\mathit{ZZ}$ denote the first and the second derivative with respect to $Z$ . Hence, the function ${\mathrm{\Phi }}_n\left(t\right)$ becomes

$${\mathrm{\Phi }}_n\left(t\right)=F{e}^{\mathit{i\theta }}+F^{\ast }{e}^{-\mathit{i\theta }}+\epsilon F_0+F_2{e}^{i2\theta }+{F_2}^{\ast }{e}^{-i2\theta }$$ E26

where ${}^{\ast }$ stands for complex conjugate and $F\equiv F\left(Z,T\right)$ . All this allows us to obtain the expressions existing in Eq. (22), such as ${\mathrm{\Phi }}_{n+1}+{\mathrm{\Phi }}_{n-1}-2{\mathrm{\Phi }}_n$ , ${\mathrm{\Phi }}_n$ and ${{\mathrm{\Phi }}_n}^3$ , and Eq. (22) becomes [28]

$$\begin{array}{l}\left({\epsilon }^3{F}_{\mathit{TT}}-2i{\epsilon }^2\omega F_T-{\mathit{\epsilon \omega }}^{2}F\right)e^{\mathit{i\theta }}-\left(4i{\epsilon }^3\omega F_{2T}+4{\epsilon }^2{\omega }^2{F}_2\right)e^{i2\theta }+\mathrm{cc}=\\ =\epsilon \frac{k}{m}\left\{2F\left[\text{cos}\left(\mathit{ql}\right)-1\right]+2\mathit{i\epsilon }{\mathit{lF}}_Z\text{sin}\left(\mathit{ql}\right)+{\epsilon }^2{l}^2{F}_{\mathit{ZZ}}\text{cos}\left(\mathit{ql}\right)\right\}{e}^{\mathit{i\theta }}\\ +\epsilon \frac{k}{m}\left\{2\epsilon F_2\left[\text{cos}\left(2\mathit{ql}\right)-1\right]+2i{\epsilon }^2{\mathit{lF}}_{2Z}\text{sin}\left(2\mathit{ql}\right)\right\}{e}^{i2\theta }+\frac{C}{m}\\ +\epsilon \frac{A}{m}\left(F{e}^{\mathit{i\theta }}+\epsilon F_0+\epsilon F_2{e}^{i2\theta }\right)-{\epsilon }^3\frac{B}{m}\left(3{\left|F\right|}^{2}F{e}^{\mathit{i\theta }}+F^3{e}^{i3\theta }\right)+\mathrm{cc}+\mathrm{O}\left({\epsilon }^4\right)\end{array}$$ E27

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 $e^{\mathit{i\theta }}$ and neglecting all the terms with ${\epsilon }^2$ and ${\epsilon }^3$ one obtains the following expressions for the dispersion relation $\omega =\omega \left(q\right)$ and the group velocity $\mathrm{d}\omega /\mathrm{d}q$ :

$${\omega }^2=\frac{2k_u}{m}\left[1-\text{cos}\left(\mathit{ql}\right)\right]-\frac{A}{m},V_g=\frac{l{k}_u}{\mathit{m\omega }}\text{sin}\left(\mathit{ql}\right)$$ E28

Also, the coefficients for $e^{i0}=1$ and $e^{i2\theta }$ _, respectively give [28]

$${\epsilon }^2{F}_0=-\frac{C}{A},F_2=0$$ E29

which yields to

$$u_n=\epsilon {\mathit{Fe}}^{i{\theta }_n}-\frac{C}{A}+\mathrm{cc}$$ E30

Eqs. (28) and (29) and new coordinates $S$ and $\tau$ , defined as $S=Z-V_{g}T$ and $\tau =\epsilon T$ , allows us to simplify Eq. (27). An explanation for why the parameter $\epsilon$ exists in the time scaling but is absent in the space scaling is given in Refs. [45, 46]. If we consider the terms for $e^{\mathit{i\theta }}$ again we obtain the well-known nonlinear Schrödinger equation (NLSE) for the function $F$

$${\mathit{\text{iF}}}_{\tau }+P{F}_{\mathit{SS}}+Q{\left|F\right|}^{2}F=0$$ E31

where the dispersion coefficient $P$ and the coefficient of nonlinearity $Q$ are

$$P=\frac{1}{2\omega }\left[\frac{k_u{l}^2}{m}\text{cos}\left(\mathit{ql}\right)-{V_g}^2\right],Q=-\frac{3B}{2\mathit{m\omega }}$$ E32

Even though Eq. (31) is PDE, its solution exists. This well-known solution, existing for $\mathit{PQ}>0$ , is [15, 47, 48]

$$F\left(S,\tau \right)=A_e\text{sech}\left(\frac{S-u_e\tau }{L_e}\right)\text{exp}\frac{{\mathit{iu}}_e\left(S-u_c\tau \right)}{2P},u_e>2u_c$$ E33

where parameters $u_e$ and $u_c$ represent envelope and carrier component velocities, while the amplitude $A_e$ and the soliton width $L_e$ have the forms

$$A_e=\sqrt{\frac{{u_e}^2-2u_e{u}_c}{2\mathit{PQ}}},L_e=\frac{2P}{\sqrt{{u_e}^2-2u_e{u}_c}}$$ E34

It is very difficult to deal with the parameters $u_e$ and $u_c$ as $u_e>2u_c$ is completely unprecise statement. However, $u_c/u_e<0.5$ seems to be more practical. Hence, new parameters $U_e$ and $\eta$ have been introduced as $U_e=\epsilon u_e$ , $\eta =u_c/u_e$ and $0\le \eta <0.5$ [45]. Finally, we can easily obtain the expression for the longitudinal displacement of the dimer at the position $n$

$$u_n\left(t\right)=A_0\text{sech}\left(\frac{\mathit{nl}-V_{e}t}{L}\right)\text{cos}\left(\mathrm{\Theta }\mathit{nl}-\mathrm{\Omega }t\right)-\frac{C}{A}\equiv U_n\left(t\right)-\frac{C}{A}$$ E35

where

$$A_0\equiv 2\epsilon A_e=2\left|U_e\right|\sqrt{\frac{1-2\eta }{2\mathit{PQ}}},L\equiv \frac{L_e}{\epsilon }=\frac{2P}{\left|U_e\right|\sqrt{1-2\eta }}$$ E36

$$V_e=V_g+U_e,\mathrm{\Theta }=q+\frac{U_e}{2P},\mathrm{\Omega }=\omega +\frac{\left({\mathrm{V}}_g+\eta U_e\right)U_e}{2P}$$ E37

One more parameter can be eliminated using the idea of coherent mode [49]. This mode means that the envelope and the carrier wave velocities are equal. It follows from Eq. (35) that $V_e=\mathrm{\Omega }/\mathrm{\Theta }$ , which yields to the function $U_e\left(\eta \right)$ . This means that the wave $u_n\left(t\right)$ is the one phase function, preserving its shape in time.

To plot the function $u_n\left(t\right)$ or, equivalently, $U_n\left(t\right)$ we should know or estimate the values of a couple of the parameters. Of course, if 2D plot is chosen, $U_n\left(t\right)$ can be presented as either a function of $t$ at a certain position $n$ or as a function of $n$ for chosen $t$ . Very detailed analysis of the parameter selection was done in Ref. [28]. One example for $q=2\pi /\mathit{Nl}$ 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.

Figure 4.

Function $U\left(n\right)$ for $t=50\mathrm{ns}$ , $A=2.9\times {10}^{-3}\mathrm{N}/\mathrm{m}$ , $B=1.7\times {10}^{14}\mathrm{N}/{\mathrm{m}}^3$ , $k_u=150A$ , $N=15$ and $\eta =0.43$ .

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 $\phi$ -model.

4. $\phi$ -model

A weak point of the U-model is the last term in Eq. (1). A scalar product $-\overrightarrow{p}\cdot \overrightarrow{E}=-\mathit{QdE}\text{cos}{\phi }_n$ would be better choice for the potential energy, where $d$ 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 $u$ and the Hamiltonian for the radial model, which we call as $\phi$ -model, is [50, 51]

$$H_{\phi }=\sum _n\left[\frac{I}{2}{{\dot{\phi }}_n}^2+\frac{k_{\phi }}{2}{\left({\phi }_{n+1}-{\phi }_n\right)}^2-\mathit{pE}\text{cos}{\phi }_n\right]$$ E38

where $I$ is a moment of inertia of the dimer at the position $n$ . Notice that the W-potential does not exist in Eq. (38) even though the terms including ${{\phi }_n}^2$ and ${{\phi }_n}^4$ 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 $M_v=-\mathrm{\Gamma }\dot{\phi }$ , where $\mathrm{\Gamma }$ is the viscosity coefficient [51, 52, 53]. Following the procedure explained earlier, we obtain

$${\alpha }_{\phi }{\Psi }^{″}-{\rho }_{\phi }{\Psi }^{\prime }+\Psi -{\Psi }^3=0$$ E39

where

$$\phi =\Psi \sqrt{6},{\alpha }_{\phi }=\frac{I{\omega }^2-k_{\phi }{l}^2{\kappa }^2}{\mathit{pE}},{\rho }_{\phi }=\frac{\omega \mathrm{\Gamma }}{\mathit{pE}}$$ E40

Like above, the solutions are kink solitons [50].

It is very interesting to compare the expressions for ${\alpha }_u$ and ${\alpha }_{\phi }$ , given by Eqs. (5) and (40). They can be written as

$${\alpha }_u=\frac{m{\kappa }^2}{A}\left(v^2-{c_u}^2\right),{c_u}^2=\frac{k_u{l}^2}{m}$$ E41

and

$${\alpha }_{\phi }=\frac{I{\kappa }^2}{\mathit{pE}}\left(v^2-{c_{\phi }}^2\right),{c_{\phi }}^2=\frac{k_{\phi }{l}^2}{I}$$ E42

where $v=\omega /\kappa$ is the soliton velocity, while $c_u$ and $c_{\phi }$ are corresponding sound velocities. According to Eq. (11), we can see that the U-model predicts $c_u>v$ as $A$ is positive. On the other hand, ${\alpha }_{\phi }=2{{\rho }_{\phi }}^2/9>0$ [50] means that, according to the $\phi$ -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 $\phi$ -model to solve the dynamical equation of motion, which is [51]

$$\epsilon I{\ddot{\mathrm{\Phi }}}_n=\epsilon k_{\phi }\left({\mathrm{\Phi }}_{n+1}+{\mathrm{\Phi }}_{n-1}-2{\mathrm{\Phi }}_n\right)-\mathit{\epsilon pE}{\mathrm{\Phi }}_n+{\epsilon }^3\mathit{pE}{{\mathrm{\Phi }}_n}^3+\mathrm{O}\left({\epsilon }^4\right)$$ E43

where $\phi =\epsilon \mathrm{\Phi }$ has been used. Of course, Eq. (43) is analog to Eq. (22). Following the procedure explained in the previous section, we straightforwardly obtain $F_0=0$ and $F_2\left(\xi \right)=0$ , as well as

$${\omega }^2={\omega }_0^2+\frac{4k}{I}{\text{sin}}^2\left(\mathit{ql}/2\right),{\omega }_0=\sqrt{\mathit{pE}/I},V_g=\frac{lk}{\mathit{I\omega }}\text{sin}\left(\mathit{ql}\right)$$ E44

where ${\omega }_0$ is the lowest frequency of the oscillations [59]. Also, we easily obtain NLSE (31), where

$$P=\frac{1}{2\omega }\left[\frac{l^{2}k}{I}\text{cos}\left(\mathit{ql}\right)-{V_g}^2\right],Q=\frac{3\mathit{pE}}{2\mathit{I\omega }}.$$ E45

The final solution ${\phi }_n\left(t\right)$ is the same as $U_n\left(t\right)$ except that $P$ and $Q$ are different. Therefore, both the U- and the $\phi$ -models predict the breather waves moving through MT.

Finally, viscosity should be introduced in the semi-discrete approximation [51]. Due to viscosity momentum $M_v=-\mathrm{\Gamma }{\dot{\mathrm{\Phi }}}_n$ _, the final result ${\phi }_n\left(t\right)$ includes the expected exponential term $e^{-\beta t}$ , where $\beta =\mathrm{\Gamma }/2I$ [51].

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 $\phi$ than the longitudinal one as we assume angular oscillations of the dimers. The scalar product $-\overrightarrow{p}\cdot \overrightarrow{E}=-\mathit{QdE}\text{cos}{\phi }_n$ , existing in the $\phi$ -model, solved these problems but the W-potential has been missing. In fact, a series expansion of $\text{cos}{\phi }_n$ gives ${{\phi }_n}^2$ and ${{\phi }_n}^4$ 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 [54]. This potential brings about ${\alpha }_{\phi }>0$ , which is disputable result.

Therefore, we want to solve the mentioned problem regarding the U-model but to keep the W-potential, the coordinate $\phi$ and, probably, $I{\omega }^2<{\mathit{kl}}^2{\kappa }^2$ , that is, $\alpha <0$ . This suggests the following Hamiltonian

$$H=\sum _n\left[\frac{I}{2}{{\dot{\phi }}_n}^2+\frac{k}{2}{\left({\phi }_{n+1}-{\phi }_n\right)}^2-\frac{A}{2}{{\phi }_n}^2+\frac{B}{4}{{\phi }_n}^4-\mathit{pE}\text{cos}{\phi }_n\right]$$ E46

where $A>0$ , $B>0$ and ${\phi }_n$ has the same meaning as in the $\phi$ -model. Let us call the model as general one (GM). The procedure mentioned earlier brings about

$$\left(I{\omega }^2-{\mathit{kl}}^2{\kappa }^2\right){\phi }^{″}-\mathrm{\Gamma }\omega {\phi }^{\prime }-\left(A-\mathit{pE}\right)\phi +\left(B-\frac{\mathit{pE}}{6}\right){\phi }^3=0$$ E47

where, of course, $\phi \equiv \phi \left(\xi \right)$ . 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: $\left(A-\mathit{pE}\right)\left(B-\frac{\mathit{pE}}{6}\right)>0$ ,      Case 2: $\left(A-\mathit{pE}\right)\left(B-\frac{\mathit{pE}}{6}\right)<0$ ,

Case 3: $A=\mathit{pE}$ , $B\ne \frac{\mathit{pE}}{6}$ ,                  Case 4: $A\ne \mathit{pE}$ , $B=\mathit{pE}/6$ .

All of them are studied in Ref. [54] and they will be explained here briefly.

Case 1 straightforwardly yields to

$${\alpha }_1{\Psi }^{″}-{\rho }_1{\Psi }^{\prime }+\Psi -{\Psi }^3=0$$ E48

where

$$\phi =\pm \sqrt{\frac{A-\mathit{pE}}{B-\mathit{pE}/6}}\Psi \equiv \mathit{K\Psi },{\alpha }_1=\frac{I{\omega }^2-{\mathit{kl}}^2{\kappa }^2}{\mathit{pE}-A},{\rho }_1=\frac{\mathrm{\Gamma }\omega }{\mathit{pE}-A}$$ E49

and the final solution is [54]

$$\phi \left(\xi \right)=\frac{K}{2}\left[1+\text{tanh}\left(\frac{3}{4{\rho }_1}\xi \right)\right],{\alpha }_1>0$$ E50

Eq. (50) holds for both positive and negative ${\rho }_1$ . Therefore, $\phi \left(\xi \right)$ represents kink soliton if ${\rho }_1>0$ and antikink one for negative ${\rho }_1$ , which is shown in Figure 5.

Figure 5.

Kink soliton for ${\rho }_1=-1$ (a) and antikink soliton for ${\rho }_1=1$ (b) for $K=1$ .

One of the advantages of the GM over the $\phi$ -model is the value of amplitude. Namely, the amplitude of the kink soliton, according to the $\phi$ -model, is $\sqrt{6}$ , coming from Eq. (40). This is unrealistic, too big value. Instead of $\sqrt{6}$ , the appropriate factor, existing in the GM, is K, given by Eq. (49).

If viscosity is neglected, the GM brings about

$${\phi }_0=K\text{tanh}\left(\xi /a\right),a^2=2{\alpha }_1$$ E51

Case 2 straightforwardly yields to

$${\alpha }_2{\Psi }^{″}-{\rho }_2{\Psi }^{\prime }+\Psi +{\Psi }^3=0,\phi =\pm \sqrt{\frac{\mathit{pE}-A}{B-\mathit{pE}/6}}\Psi \equiv K^{\prime }\Psi$$ E52

and to the final results

$${\phi }_2=i\frac{K^{\prime }}{2}\left[1+\text{tanh}\left(\frac{3}{4{\rho }_2}\xi \right)\right],{\rho }_2\ne 0$$ E53

and

$${\phi }_{20}=K^{\prime }\text{tan}\left(\xi /a\right),a^2=-2\alpha ,{\rho }_2=0$$ E54

It is obvious that these results do not have physical meaning as ${\phi }_2$ is complex, while ${\phi }_{20}$ may diverge.

Case 3 brings about

$${\alpha }_3{\Psi }^{″}-{\rho }_3{\Psi }^{\prime }+{\Psi }^3=0,{\alpha }_3=\frac{I{\omega }^2-{\mathit{kl}}^2{\kappa }^2}{B-\mathit{pE}/6},{\rho }_3=\frac{\mathrm{\Gamma }\omega }{B-\mathit{pE}/6}$$ E55

as well as $a_0=a={\alpha }_3={\rho }_3=0$ , 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 $\phi$ -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 $\phi$ -model predicts the supersonic wave. How about the GM? It was shown that Case 1 yields to the solutions having physical sense and that ${\alpha }_1>0$ . According to Eq. (49), we easily reach the final conclusion:

1. If $A>\mathit{pE}$ and $B>\mathit{pE}/6$ _, then ${\rho }_1<0$ and the function $\phi \left(x,t\right)$ is subsonic soliton, kink for the positive $K$ in Eq. (50) and antikink otherwise.

2. If $A<\mathit{pE}$ and $B<\mathit{pE}/6$ _, then ${\rho }_1>0$ and the function $\phi \left(x,t\right)$ is supersonic soliton, antikink for the positive $K$ 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 $\phi$ -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 $A$ , that is, $A>\mathit{pE}$ , 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 [55]. 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 [58].

Finally, we should bear in mind the cytological and medical applications of the research explained in this chapter [59, 60, 61].

Acknowledgments

This work was supported by funds from Serbian Ministry of Education, Sciences and Technological Development (grant No. III45010).