Modified extended tanh-function method and nonlinear dynamics of microtubules

Modified extended tanh-function method and nonlinear dynamics of microtubules Slobodan Zdravkovića,*, Louis Kavithab,c, Miljko V. Satarićd, Slobodan Zekovića, Jovana Petrovića a Institut za nuklearne nauke Vinča, Univerzitet u Beogradu, Poštanski fah 522, 11001 Beograd, Serbia b Department of Physics, Periyar University, Salem-636 011, India c The Abdus Salam International Centre for Theoretical Physics, Trieste, Italy d Fakultet tehničkih nauka, Univerzitet u Novom Sadu, 21000 Novi Sad, Serbia ABSTRACT We here present a model of nonlinear dynamics of microtubules (MT) in the context of modified extended tanh-function (METHF) method. We rely on the ferroelectric model of MTs published earlier by Satarić et al [1] where the motion of MT subunits is reduced to a single longitudinal degree of freedom per dimer. It is shown that such nonlinear model can lead to existence of kink solitons moving along the MTs. An analytical solution of the basic equation, describing MT dynamics, was compared with the numerical one and a perfect agreement was demonstrated. It is now clearer how the values of the basic parameters of the model, proportional to viscosity and internal electric field, impact MT dynamics. Finally, we offer a possible scenario of how living cells utilize these kinks as signaling tools for regulation of cellular traffic as well as MT depolymerisation. 1. Introduction Microtubules are major cytoskeletal proteins. They are hollow cylinders formed by protofilaments (PF) representing series of proteins known as tubulin dimers [2,3]. There are usually 13 longitudinal PFs covering the cylindrical walls of MTs. The inner and the outer diameters of the cylinder are 15nm and 25nm , while its length may span dimensions from the order of micrometer to the order of millimetre. Each dimer is an electric dipole whose length and longitudinal component of the electric dipole moment are l = 8nm [2-4] and p = 337 Debye [5], respectively. The constituent parts of the dimers are α and β tubulins, corresponding to positively and negatively charged sides, respectively [2-4]. * Corresponding author. E-mail addresses: szdjidji@vinca.rs (S.Zdr.), louiskavitha@yahoo.co.in (L.K.), bomisat@neobee.net (M.V.S.), zekovic@vinca.rs (S.Zek.), jovanap@vinca.rs 1 In this paper we demonstrate how METHF method [6-10] can be used in the study of nonlinear dynamics of MTs. The paper is organized as follows. In Section 2 we explain the well known model for MTs we rely on [1]. The modification of the model presented in this paper is a generalization of the original one and will be referred to as u-model. The model brings about a crucial nonlinear differential equation, describing nonlinear dynamics of MTs. In Section 3 we briefly describe METHF method. Then we solve the basic nonlinear differential equation, mentioned above. We show that its solution is a kink-like solitonic wave. This result is compared with numerical solutions in Section 4. Finally, in Section 5, we give concluding remarks. In particular, we emphasize the biological importance of the studied kink-like solitons. 2. U-model of MTs The model we rely on assumes only one degree of freedom of dimers motion within the PF [1]. This is a longitudinal displacement of a dimer at a position n denoted as u n and thus we call the model as u-model. The overall effect of the surrounding dimers on a dipole at a chosen site n can be described by a double-well potential [1] Vd (u n ) = − 1 1 Au n2 + Bu n4 2 4 (1) where A and B are positive parameters that should be estimated. As an electrical dipole, a dimer in the intrinsic electric field of the MT acquires the additional potential energy given by [1] Vel (u n ) = −Cu n , C = qE , (2) where E is the magnitude of the intrinsic electric field at the site n, as the dimer n exists in the electric field of all other dimers, and q represents the excess charge within the dipole. It is assumed that q > 0 and E > 0 . The Hamiltonian for one PF is represented as follows ⎤ ⎡m 2 k 2 H = ∑ ⎢ u& n + (u n +1 − u n ) + V (u n )⎥ , 2 ⎦ n ⎣2 (3) where dot means the first derivative with respect to time, m is a mass of the dimer, k is an intra-dimer stiffness parameter and the integer n determines the position of the considered dimer in the PF [1]. The first term represents a kinetic energy of the dimer, the second one is a potential energy of the chemical interaction between the neighbouring dimers belonging to the same PF and the last term is the combined potential V (u n ) = −Cu n − 1 1 Au n2 + Bu n4 . 2 4 (4) 2 It is obvious that the nearest neighbour approximation is used. However, this does not mean that the influence of the neighbouring PFs is completely ignored as the value of the electric field E depends also on the dipoles belonging to the neighbouring PFs. By using the generalized coordinates q n and p n , defined as q n = u n and p n = mu& n , applying a continuum approximation u n (t ) → u ( x, t ) and making a series expansion u n ±1 → u ± 1 ∂ 2u 2 ∂u l l+ 2 ∂x 2 ∂x (5) we can straightforwardly obtain an appropriate dynamical equation of motion. In order to derive a realistic equation, the viscosity of the solvent should also be taken into consideration. This can be achieved by introducing a viscosity force Fv = −γ u& into the obtained dynamical equation of motion, where γ is a viscosity coefficient [1]. All this brings about the following nonlinear partial differential equation m 2 ∂u ∂ 2u 2 ∂ u − = 0. − qE − Au + Bu 3 + γ k l 2 2 ∂t ∂t ∂x (6) It is well known that, for a given wave equation, a travelling wave u (ξ ) is a solution which depends upon x and t only through a unified variable ξ ξ ≡ κ x −ωt , (7) where κ and ω are constants. Substitution of x and t by ξ transforms Eq. (6) into the following ordinary differential equation (ODE) (mω 2 ) − k l 2κ 2 u ′′ − γω u ′ − Au + Bu 3 − qE = 0 . (8) By introducing a dimensionless function ψ through the relation u= A ψ, B (9) a much more convenient equation can be obtained. This final ODE reads αψ ′′ − ρψ ′ − ψ + ψ 3 − σ = 0 , (10) and contains the following new parameters: mω 2 − kl 2κ 2 α= , A (11) 3 σ= ρ= qE A A B , γω (12) (13) A du . dξ It was already mentioned that this approach represented a certain improvement of the original model, explained in [1]. If we compare Eq. (10) with the appropriate one in [1] we can see that they are equal for α = −1 . Therefore, this approach is more general. We treat the parameters ρ and σ as an input and will determine values of dinamical parameters of the system. We will see that the final result depends on ρ and σ only, i.e. on the parameters that determine their values. The crucial equation (10) will be solved in the next section. Before we proceed we want to discuss the potential energy V (u ) , defined by Eq. (4). This step is very important to understand the physics behind Eq. (10) and its solutions. Using the procedure explained above we can easily obtain the following convenient expression for this potential and u ′ ≡ V (ψ ) = A2 f (ψ ) , B (14) where 1 1 f (ψ ) = −σψ − ψ 2 + ψ 4 . 2 4 (15) The function f (ψ ) is shown in Fig. 1 for two values of the parameter σ . For σ = 0 the function f (ψ ) and, consequently, the potential V (ψ ) , is symmetric (curve a) while for the increasing σ the right minimum becomes deeper and the left one becomes shallower and elevated. To find the values of ψ for which f (ψ ) reaches a maximum and minima we should solve the equation f ′(ψ ) = −σ − ψ + ψ 3 = 0 . (16) 4 0.0 -0.2 a a f (ψ ) b -0.4 -0.6 -1.0 -0.5 0.0 ψ 0.5 1.0 Fig. 1. The function f (ψ ) for: (a) σ = 0 and (b) σ = 0.3 . According to the procedure explained in Appendix we can easily obtain the following roots of Eq. (16): ψR = 2 3 cos F , ψ max = 1 ψL = − 1 3 3 (− cos F + (cos F + (17) ) 3 sin F , ) 3 sin F , (18) (19) where ⎛σ ⎞ 1 F = arccos⎜⎜ ⎟⎟ . 3 ⎝σ0 ⎠ (20) The functions ψ R and ψ L correspond to the right and left minimum of the function f (ψ ) and a critical value of the parameter σ has the value 5 σ0 = 2 3 3 . (21) This means that the three real roots of Eq. (16) exist for σ < σ 0 . The dependence ψ ex (σ ) is shown in Fig. 2, where ψ ex stands for ψ R , ψ max and ψ L . We see that ψ R is a slowly increasing function of σ while ψ L and ψ max approach each other for the increasing σ . This means that the right minimum moves to the right for the higher σ while the maximum and the left minimum of the function f (ψ ) become closer. All this can be noticed in Fig. 1. For σ = σ 0 the left minimum and the maximum disappear, coalescing in the saddle point. ΨR 1.0 0.5 Ψex(σ) σ0 0.0 Ψmax -0.5 ΨL -1.0 0.0 0.1 0.2 σ 0.3 0.4 Fig. 2. The values of ψ ex corresponding to the extrema of the potential V (ψ ) as a function of the parameter σ . The values ψ R and ψ max for a couple of values of σ are shown in Table 1. Of course, ψ max does not exist for σ = σ 0 . 6 Table 1 σ 0 0.1 0.2 0.3 ψR 1 1.047 1.088 1.125 1.155 0 -0.101 -0.209 -0.339 na ψ max σ0 All this suggests that the higher σ , i.e. the bigger value of qE , increases stability of MTs dynamics around the right minimum which becomes deeper. This requires further research and should be checked by stability analysis which will be a topic of a separate publication. 3. Modified extended tanh-function method In what follows we briefly outline the basic features of METHF method. The method will be applied for solving Eq. (10). According to this procedure we look for the possible solution in the form [6,7] ψ = a 0 + ∑ (ai Φ i + bi Φ −i ) , M (22) i =1 where the function Φ = Φ (ξ ) is a solution of the well known Riccati equation Φ′ = b + Φ 2 (23) and Φ ′ is the first derivative [6,7]. The parameters a 0 , ai , bi and b are real constants that should be determined as well as an integer M . The possible solutions of (23) depend on the parameter b as follows: ( ) ( ) 1) If b > 0 then Φ = b tan b ξ , or Φ = − b cot b ξ , 1 2) If b = 0 then Φ = − , ξ ( ) (24) (25) ( ) 3) If b < 0 then Φ = − − b tanh − b ξ , or Φ = − − b coth − b ξ . (26) Balancing the orders of ψ ′′ and ψ 3 with respect to the new function Φ we find M = 1 . Namely, the highest orders of the function Φ in the expressions for ψ ′′ and ψ 3 are Φ M + 2 and Φ 3M respectively and they are equal for M = 1 . Eqs. (22) and (23) yield the expressions for ψ ′ , ψ ′′ and ψ 3 as functions of Φ . For example, the second derivative of ψ is 7 ψ ′′ = 2a1bΦ + 2a1Φ 3 + 2b1b 2 Φ −3 + 2b1bΦ −1 . (27) Eq. (10) in terms of the expressions for ψ , ψ ′ , ψ ′′ and ψ 3 brings about the crucial equation A1Φ + B1Φ −1 + A2 Φ 2 + B2 Φ −2 + A3 Φ 3 + B3 Φ −3 + A0 = 0 , (28) where the following set of abbreviations is introduced: A0 = −a0 + a0 + 6a0 a1b1 − ρ a1b + ρ b1 − σ , (29) A1 = −a1 + 3a0 a1 + 2α a1b + 3a1 b1 , (30) B1 = −b1 + 3a0 b1 + 2α bb1 + 3a1b1 , (31) A2 = 3a0 a1 − ρ a1 , (32) B2 = 3a0 b1 + ρ bb1 , (33) A3 = 2α a1 + a1 (34) 3 2 2 2 2 2 2 3 and B3 = 2α b 2 b1 + b1 . 3 (35) Of course, Eq. (28) is satisfied if all these coefficients are simultaneously equal to zero which renders a system of seven equations. Before embarking on solving such a system, we examine its behaviour under the conditions given by Eqs. (24)-(26). The solutions expressed through tangents and cotangents cannot be biophysically tractable as these functions diverge. The same can be said for hyperbolic cotangent. This means that we are looking for the acceptable solutions for which b < 0 , a1 ≠ 0 and b1 = 0 , which reduces the system mentioned above. Hence, the system can be reduced to the following system of four equations − a 0 + a 0 − ρ a1b − σ = 0,⎫ ⎪ 2 − 1 + 3a 0 + 2α b = 0, ⎪ ⎬ ρ = 3a 0 a1 , ⎪ 2 ⎪ 2α = −a1 . ⎭ 3 (36) Its solutions, i.e. the values of the parameters b , a 0 , a1 and α , are given through 8 8a0 − 2a0 + σ = 0 , 3 a1 = ρ , 3a 0 (37) (38) 2 a α =− 1 2 (39) and 3a 0 − 1 2 b= a1 2 . (40) Notice that a 0 a1 > 0 and α < 0 . Also, b < 0 holds for the inequality 1 . 3 a0 < 2 (41) It is easy to check that the requirement (41) is equivalent to σ < σ 0 , which was discussed earlier. It was mentioned above that this extended version of the model reduces to the model from [1] if α = −1 . To calculate the value for α we have to know the values of the parameters determining σ and ρ . The values of these parameters have not been determined yet. It suffices now to discuss a negative sign of α , as can be seen from (39). This parameter comes from the first two terms in Eq. (6). These are inertial and elastic terms. A question is which contribution is higher and the negative value of α suggests that the elastic term in (6) is bigger than the inertial one. Using the expression for wave speed v= ω κ (42) Eq. (11) can be written as α= κ2 (mv A 2 − kl 2 ) . (43) The fact that α < 0 may indicate a small velocity of the wave and/or a big k , i.e. strong chemical bond between the neighbouring dimers in PF, as can be seen from (3). The importance of α and our intention to get as much information about it as possible is reason that we did not perform rescaling of Eq. (10), which would remove α from the 9 equation. Notice that we cannot a priori exclude the value α = 0 , which prevents rescaling. Eq. (37) looks like (16) and Appendix should be used again. The requirement for the existence of the three real solutions is σ < σ 0 again. Hence, the roots of (37) are a 01 = a 02 = 1 2 3 1 2 3 a 03 = − 1 3 (cos F + 3 sin F , ) (44) (cos F − 3 sin F , ) (45) cos F (46) where the expressions for F and σ 0 are given by (20) and (21). Of course, all these three values depend on σ and they are shown in Fig. 3. One can see that (37) has three real solutions for σ < σ 0 and only one ( a 03 ) for σ > σ 0 . To obtain Eqs. (44)-(46) a relationship cos −1 ( β ) + cos −1 (− β ) = π (47) was used. A next step is to construct the function ψ (ξ ) . Suppose that ψ 1 describes the system in the state corresponding to the right minimum in Fig. 1. Using Eqs. (26), (38) and (40) we easily obtain ( ) Φ = − − b tanh − b ξ = − 3a01 ρ ⎞ ⎛ 3a 2 2 1 − 3a01 tanh⎜⎜ 01 1 − 3a 01 ξ ⎟⎟ . ⎠ ⎝ ρ (48) Finally, Eqs. (22) with M = 1 and (48) yield the following expression for the function ψ1 : ⎛ 3a ⎞ ψ 1 (ξ ) = a 01 − 1 − 3a 01 2 tanh⎜⎜ 01 1 − 3a 01 2 ξ ⎟⎟ , ⎠ ⎝ ρ (49) where a 01 is determined by (44), (20) and (21). We can see that ψ 1 (ξ ) depends on ρ and σ only as a 01 is a function of σ . Notice that the jump of the function ψ 1 (ξ ) from − ∞ to + ∞ depends on σ only while the solitonic width, i.e. slope, depends on both ρ and σ . It is obvious that the solitonic width is proportional to viscosity. 10 1.0 a01 0.5 a0(σ) a02 σ0 0.0 -0.5 a03 -1.0 0.0 0.1 0.2 σ 0.3 0.4 Fig. 3. The values of a 0 as a function of the parameter σ . The function (49) is shown in Fig. 4 (analytic curves) for the two values of the parameter σ and for ρ = 1 . Obviously, this is an anti-kink soliton. The parameter ρ affects its slope, i.e. the width of the soliton, but not the character of the wave. According to Eq. (49) we can calculate the initial and the final value of ψ 1 that is ψ 1 (−∞) = a 01 + 1 − 3a01 2 (50) and ψ 1 (+∞) = a 01 − 1 − 3a01 2 . (51) Using Eq. (44) we can straightforwardly prove that ψ 1 (−∞) = ψ R , ψ 1 (+∞) = ψ max (52) as sin F − 3 cos F is negative. 11 Fig. 4. The function ψ 1 (ξ ) for ρ = 1 and: (a) σ = 0 and (b) σ = 0.3 . Therefore, nonlinear dynamics of MTs is described by antikink solitons. In other words, the system exists at the right minimum of Fig. 1 and the value of ψ (ξ ) is ψ 1 (−∞) = ψ R . When some energy is supplied, released by hydrolysis of guanosine triphosphate (GTP), the system jumps to the maximum and ψ becomes ψ max , as explained by (52). The way how this transition is performed is explained by the solitonic wave (49). Of course, sooner or later, the system spontaneously returns to its minimum and so on. These transitions represent displacements of the dimer, described by u , which is defined by Eq. (9). Obviously, for a quantitative treatment we should know the values of the parameters. Some estimates have already been done [1] but this is still an open question. We should keep in mind that the model described in this paper is somewhat different from the original one, introduced in [1]. For example, Eqs. (11), (38) and (39) represent a relationship between some parameters and will help in their estimations. Anyway, the parameter selection requires further research and is not the topic of this work. A careful reader may ask why only the function ψ 1 (ξ ) has been studied in this paper. There are three solutions of (37) and, consequently, three possible functions ψ (ξ ) . However, it does not seem that the functions ψ 2 (ξ ) and ψ 3 (ξ ) , corresponding to a 02 and a 03 , are physically relevant. For example, one can easily show that ψ 3 (−∞) = ψ L , ψ 3 (+∞) = ψ max , (53) 12 which could be expected. Obviously, the soliton ψ 1 (ξ ) , corresponding to the deeper energy minimum, is more relevant. Also, we can show that ψ 2 (−∞) = ψ R , ψ 2 (+∞) = ψ L , (54) which may be physically relevant only for σ < σ 0 . Therefore, we study in some more details only one out of the three possible solutions as this one is physically important. Also, we believe that the solution ψ 1 (ξ ) is relevant for triggering mechanism, which will be explained in the last section. We do not study the case when only one real solution of (37) exists. Such the case corresponds to the positive b , which does not have physical meaning, as was explained above. Due to the large uncertainty in the values of the parameters A and B of the potential (4), we concentrate on the characteristics of the soliton that do not depend on these parameters. An important example is the wave velocity. After a simple algebra with Eqs. (11), (13), (36) and (42), it can be expressed as kl 2 v2 = m+ γ2 . (55) 2 18a 0 A We consider the solution a 0 = a 01 given by Eqs. (44), (20) and (21). According to the function ψ 1 (ξ ) given by (49) we can write 3a 0 ρ 1 − 3a 0 ξ = 2 3a 0 ρ 1 − 3a 0 κ ( x − vt ) ≡ χ ( x − vt ) . 2 (56) Hence, χ= 3a 0κA 1 − 3a 0 γω 2 . (57) Obviously, χ determines the solitonic width Λ through [11] χ= 2π . Λ (58) It is convenient to write Λ as an integer of dimer’s length l , i.e. Λ = Nl . (59) All this brings about the following expression for v : 13 v= 3a 0 ANl 1 − 3a 0 2 . 2πγ (60) We can eliminate A from (55) and (60) and obtain mv 2 + C 0 v − kl 2 = 0 (61) where C0 = N lγ 1 − 3a0 2 12π a 0 . (62) Hence, the expression for v reads v= C 0 ⎛⎜ 4mkl 2 −1+ 1+ 2 2m ⎜⎝ C0 ⎞ ⎟. ⎟ ⎠ (63) It is obvious that the speed v does not depend on A and B but does depend on some other parameters. Using (62) one can see that v is the decreasing functon of viscosity γ . To estimate v we can use [1]: m = 11 ⋅ 10 −23 kg and γ = 5.3 ⋅ 10 −11 kg s but the values of N and k are not known. Also, of special importance is the intensity of electric field E as this determines a 0 through σ . Namely, Eq. (12) can be written as σ= qE A A B = qdE A dA B = pE A dA B , (64) where p ≡ qd is the longitudinal component of the electric dipole moment and the distance between the centres of positive and negative charges in any dimer is d ≈ 4nm [5]. Therefore, the estimates of E , A and B are of crucial importance. An attempt to estimate these parameters is in progress even though some estimates were done in [1]. Some values of the velocity v for N = 30 and four values of k are given in Table 2 for different values of σ (different strengths of the intrinsic electric field E ). 14 Table 2 σ v(m v(m v(m v(m 0 s) s) s) s) 18.9 92.1 179.2 341.4 0.1 0.2 0.3 σ0 15.6 76.6 150.2 290.0 12.8 63.4 125.1 243.9 10.2 50.5 100.1 197.0 6.3 31.5 62.8 124.8 k k k k = 0.1 N m = 0.5 N m = 1N m = 2N m 4. Numerical solutions In order to test the described analytical procedure and to study behaviour of the soliton for those parameters of the MT for which analytical solutions do not exist, we have solved Eq. (10) numerically. The standard shooting method with the 4th order RungeKutta integrator was used. Particular solutions are determined by their asymptotic behaviour and are centred at x = 0 . The numerical step was 10 −5 . Numerical solutions were first compared to the analytical solutions obtained by METHF method. Excellent agreement shown in Fig. 4 validates our analytical approach. We further generated numerical solutions of Eq. (10) for those parameters α for which our analytical procedure does not give solutions. To enable a direct comparison, we kept the parameters ρ and σ constant and again considered dimers within a MT with the intrinsic electric field zero ( σ = 0 ) and ρ = 1 , and with a non-zero electric field ( σ = 0.3 ) and ρ = 1 . The results are shown in Fig. 5 and Fig. 6, respectively. The values for a 01 were calculated according to (37). As the parameter α multiplies only a derivative in Eq. (10) it does not influence asymptotic solutions. These are determined by σ only. For small absolute values of α the soliton is narrow and has a shape similar to tanh function, meaning that the transition between the initial and the final position of the dimer is fast and smooth. On the other hand, for big α the dimmer experiences oscillation before it stabilizes in the final value. The period and amplitude of oscillations increase with α . These results indicate that there is an optimal parameter α for a given process to occur. 15 Fig. 5. The function ψ 1 (ξ ) for ρ = 1 and σ = 0 for different values of α Fig. 6. The function ψ 1 (ξ ) for ρ = 1 and σ = 0.3 for different values of α 16 5. Conclusions In this paper we thoroughly revisited an earlier work [1] dedicated to nonlinear dynamics of tubulin dimers within the same PF of a stable cellular MT. The basis of the model is the quartic nonlinear potential describing the mean field of chemical bonds between neighbouring tubulin dimers. This potential does dominantly depend on longitudinal displacements u n of the tubulin dimers along MT axis, as shown by Eq. (1). Inclusion of dipole-dipole interaction between tubulin dimers as well as the impact of viscous damping of cytosol has led to dimensionless nonlinear ODE (10). The effective potential (15) was discussed respecting the parameter σ which is proportional to the magnitude of intrinsic electric field created by MT itself. To solve Eq. (10) we have relied on METHF method [6-10]. The physically tractable solutions have the shape of anti-kink soliton (49), well known in nonlinear ferroelectrics and ferromagnetics. These solutions for the two values of σ are compared with the numerical solutions and a perfect agreement was demonstrated. It might be important to point out that the Hamiltonian is only what was taken from Ref. [1]. In addition, the way how the viscosity term was introduced in the dynamical equation of motion is widely accepted and can be found in many references. The procedure explained in this paper is completely different from the one in Ref. [1]. It turned out that the METHF method is both elegant and useful. It enables us to estimate the maximum of electric field E . Namely, the highest value of σ is σ 0 , given by (21). According to Eq. (64), we see that this is proportional to E. As the values of p , d , A and B can be found in literature we can calculate the maximum of E . This is not the topic of this paper as we believe that the values of the parameters are still an open problem, requiring further research. It was mentioned above that the crucial ODE for ψ in this paper and in Ref. [1] are equal only if α = −1 , which certainly means that this approach is more general. One can say that the final result of the both papers, this one and Ref. [1], is the function ψ (ξ ) . We showed only the function ψ 1 (ξ ) describing dynamics of MT when the system is in the deeper minimum of the potential well. This is shown in Figs. 4-6. There is one very important advantage of the method applied here over the one used in Ref. [1]. If we look at Eq. (3.11) of Ref. [1] we see that ψ (ξ ) does not depend neither on ρ nor on σ , which is a big problem. In fact, Eq. (3.18) was an attempt to introduce the dependence of electric field but it turned out that the term with σ was negligible. As for ψ 1 (ξ ) , existing in this paper, it depends on both parameters as a 01 depends on σ . This is the big advantage as ρ and σ , being proportional to viscosity and internal electric field, are internal parameters and we certainly expect that the final result, i.e. dynamics of MT, should depend on them. In addition, the expression for ψ 1 (ξ ) shows that the solitonic width is proportional to viscosity, which is physically plausible. Also, Figs. 4-6 show how σ affects the final result. Kink velocities are estimated in both in this paper and in Ref. [1]. We can see that these estimations are different. The estimated values of solitonic speed here are in the range from about 10m s to a few hundreds of m s . This interval may be comparable to 17 the speed of nerve axon potential, spanning to 30m s , but is much greater than the speed of MT motor protein which is, in a case of kinesin, around 1µm s , or ionic waves within the cell (a few µm s ). We expect that these stable solitons could be elicited for example by the energy released in hydrolysis of tubulin-nucleotide GTP, or by the interaction with other MT associated proteins. We here offer the ideas of how these solitons could be utilized by a cell for some important mechanisms underlying its functional dynamics. First, it is plausible to expect that anti-kink soliton reaching the MT tip could cause the detachment of last tubulin dimer. This is one of triggering mechanisms responsible for depolymerisation of MTs. Second, it is well known that MTs serve as a “road network” for motor proteins (kinesin and dynein) dragging different “cargos” such as vesicles and mitochondria to different sub-cellular locations. Kinesin motors move along MT towards MT plus-end and dynein motors move in opposite direction. Both classes of motors are present on the same “cargo” [12] and it is still unclear which mechanism decides which motor will take the action and drag the cargo in the proper direction. We argue that the cell’s compartment which needs the specific cargo will launch the soliton of above kind sending it as a signal along the closest MT. The soliton will activate proper motors being close to MT along which the soliton propagates. For example, if the need for mitochondria is expressed in the regions close to the cell membrane, the solitons created in these MTs will propagate down MTs towards cell nucleus. These solitons will turn on kinesin motors in order to attach to the same MTs and to carry mitochondria towards MT plus-ends which are close to cell membrane. It might be useful to mention that some similar ideas were elaborated earlier [13,14]. Appendix A polynomial x 3 + px + q = 0 , p<0 (A.1) has three real roots for D < 0 or one real and two complex conjugate roots for D > 0 , where [15] D= q2 p3 + . 4 27 (A.2) For the negative D the solutions of (A.1) are [15] ⎛ ϕ + 2kπ ⎞ x k = 2 3 r cos⎜ ⎟, 3 ⎠ ⎝ k = 0,1, 2 , (A.3) where 18 r= − p3 27 (A.4) and cos ϕ = − q . 2r (A.5) For the positive D the solution of (A.1) is [15] p 3 x=− sin( 2ϕ ) − (A.6) where ⎛ w⎞ tan ϕ = 3 tan⎜ ⎟ ⎝2⎠ (A.7) and 3 ⎤ ⎡ 2 2 p ⎛ ⎞ ⎢ w = arcsin ⎜ − ⎟ ⎥ . ⎢q ⎝ 3 ⎠ ⎥ ⎦ ⎣ (A.8) Acknowledgements This research was supported by funds from Serbian Ministry of Sciences, grants III45010 and OI171009. L. Kavitha gratefully acknowledges the financial support from UGC in the form of major research project, BRNS, India in the form of Young Scientist Research Award and ICTP, Italy in the form of Junior Associateship. S. Zdravković gratefully acknowledges the hospitality of the Abdus Salam International Centre for Theoretical Physics, Trieste, Italy, where a big part of this work was done. References [1] [2] Satarić MV, Tuszyńsky JA, Žakula RB. Kink like excitations as an energytransfer mechanism in microtubules. Phys. Rev. E 1993; 48:589-597. Dustin P. Microtubules. Berlin: Springer; 1984. 19 [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] Tuszyńsky JA, Hameroff S, Satarić MV, Trpisová B, Nip MLA. Ferroelectric Behavior in Microtubule Dipole Lattices: Implications for Information Processing, Signaling and Assembly/Disassembly. J. Theor. Biol. 1995; 174:371380. Satarić MV, Tuszyńsky JA. Nonlinear Dynamics of Microtubules: Biophysical Implications. J. Biol. Phys. 2005; 31:487-500. Schoutens JE. Dipole–Dipole Interactions in Microtubules. J. Biol. Phys. 2005; 31:35–55. Ali AHA. The modified extended tanh-function method for solving coupled MKdV and coupled Hirota–Satsuma coupled KdV equations. Phys. Lett. A 2007; 363:420-425. El-Wakil SA, Abdou MA. New exact traveling wave solutions using modified extended tanh-function method. Chaos, Solitons Fractals 2007; 31:840-852. Kavitha L, Prabhu A, Gopi D. New exact shape changing solitary solutions of a generalized Hirota equation with nonlinear inhomogeneities. Chaos, Solitons Fractals 2009; 42:2322–2329. Kavitha L, Akila N, Prabhu A, Kuzmanovska-Barandovska O, Gopi D. Exact solitary solutions of an inhomogeneous modified nonlinear Schrödinger equation with competing nonlinearities. Math. Comput. Modelling. 2011; 53:1095–1110. Kavitha L, Srividya B, Gopi D. Effect of nonlinear inhomogeneity on the creation and annihilation of magnetic soliton. J. Magn. Magn. Mater. 2010; 322:1793– 1810. Zdravković S, Satarić MV. Transverse interaction in DNA molecule. BioSystems 2011; 105:10–13. Gross SP. Hither and yon: a review of bi-directional microtubule-based transport. Phys. Biol. 2004; 1:R1-R11. Satarić MV, Budimski-Petković L, Lončarević I, Tuszyńsky JA. Modeling the role of intrinsic electric fields in microtubules as an additional control mechanism of bi-directonal intracellular transport. Cell, Biochem. Biophys. 2008; 52:113124. Satarić MV, Kozmidis-Luburić U, Budimski-Petković L, Lončarević I. Intrinsic electric field as a control mechanism of intracellular transport along microtubules. J. Comput. Theor. Nanosci. 2009; 6:721-731. Smirnov VI. Kurs vishey matematiki, tom 1. Moscow: Nauka 1965. (In Russian). 20