Applied Mathematics and Computation 242 (2014) 353–360
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Are microtubules discrete or continuum systems? S. Zdravkovic´ a,⇑, A. Maluckov a, M. Ðekic´ b, S. Kuzmanovic´ c, M.V. Sataric´ d a
Institut za nuklearne nauke Vincˇa, Univerzitet u Beogradu, Poštanski fah 522, 11001 Beograd, Serbia Institut za multidisciplinarna istrazˇivanja, Univerzitet u Beogradu, Ulica kneza Višeslava 1, 11030 Beograd, Serbia Prirodno-matematicˇki fakultet, Univerzitet u Prištini, Kosovska Mitrovica, Serbia d Fakultet tehnicˇkih nauka, Univerzitet u Novom Sadu, 21000 Novi Sad, Serbia b c
a r t i c l e
i n f o
Keywords: Microtubule Kink soliton Method of factorization Continuum approximation
a b s t r a c t In this paper we study nonlinear dynamics of microtubules (MTs) relying on so-called u-model. A crucial discrete differential equation is transformed into a partial differential equation using a continuum approximation. Both the continuum and the discrete equations are solved and an excellent agreement of the results shows that MTs can be viewed as continuum systems. Thus we proved that the continuum approximation is applicable in research of nonlinear dynamics of MTs. Ó 2014 Published by Elsevier Inc.
1. Introduction Microtubules (MTs) are crucial polymerized structures comprising a cytoskeleton. They also serve as a network for motor proteins. Their structure is known and described in many books and papers [1–6]. In the present work we offer only some essential information relevant for this paper. MTs are hollow cylinders whose outer and inner radii are about 25 nm and 17 nm, respectively. Their length may span dimensions from the order of micrometer to the order of millimetre. The cylindrical walls are formed usually by 13 parallel protofilaments (PFs). Each PF represents a series of proteins called tubulin dimers, with the length of l ¼ 8 nm. They are electric dipoles whose longitudinal component of the electric dipole moment is p ¼ 337 Debye [7]. This means that MTs are ferroelectrics, which is crucial point for the choice of an appropriate model describing their dynamics. The oldest nonlinear model describing complicated dynamics of MTs is probably the one introduced about 20 years ago [8]. The model assumes only one longitudinal degree of freedom per dimer. In fact, according to the model, the dimer rotates but the used coordinate un is a projection of the top of the dimer at the position n on the direction of PF. Its improved and more general version, which we call u-model, was introduced recently [9]. We want to mention two more models. They also assume one degree of freedom per dimer. This coordinate is either an angle [10] or a longitudinal shift of the dimer [6]. We call them as u-model and z-model, respectively. An important simplification, used in all the models mentioned above, is based on a fact that the interaction between the neighboring dimers belonging to the same PF is much stronger than the interaction between the dimers belonging to neighboring PFs [11,12]. This means that the Hamiltonian of MT is nothing but the Hamiltonian of a single PF. However, electric field strength existing in the Hamiltonian comes from all the dimers, which means that the influence of the dimers from the neighboring PFs is taken into consideration through this field.
⇑ Corresponding author. E-mail address:
[email protected] (S. Zdravkovic´). http://dx.doi.org/10.1016/j.amc.2014.05.068 0096-3003/Ó 2014 Published by Elsevier Inc.
S. Zdravkovic´ et al. / Applied Mathematics and Computation 242 (2014) 353–360
354
In Section 2, we very briefly describe the u-model, which is important for this paper. We use a continuum approximation. A crucial differential equation is solved in Section 3 both analytically, using a method of factorization, and numerically. Almost perfect agreement between the results of these two approaches is demonstrated. A main goal of this paper is to study if the continuum approximation can be applied in the MT research. This is a crucial question for researchers working in theoretical biophysics. We solved an initial discrete partial differential equation in Section 4 and showed that its solutions are very close to the analytical results, which proves our expectations regarding the continuum approximation. The paper is closed by concluding remarks in Section 5. 2. u-Model of MT dynamics If un is the longitudinal coordinate of the dimer at the position n then the Hamiltonian for MT is
H¼
Xm n
k u_ 2n þ ðunþ1 un Þ2 qEun þ Vðun Þ ; 2 2
ð1Þ
where dot means the first derivative with respect to time, m is a mass of the dimer and k is an intra-dimer stiffness parameter [8,9]. The first three terms respectively represent a kinetic energy of the dimer, a potential energy of the chemical interaction between the neighboring dimers belonging to the same PF and a potential energy of the electric dipole existing in the electric field of all other dimers of the MT. Obviously, E is the magnitude of the intrinsic electric field at the site n and q represents the excess charge within the dipole. It is assumed that q > 0 and E > 0. Finally, the last term in the Hamiltonian is a well known double-well potential [8]
1 1 Vðun Þ ¼ Au2n þ Bu4n ; 2 4
ð2Þ
which is rather common in physics and used to study dynamics of ferroelectric systems [13–15]. Positive parameters A and B should be estimated [8]. Using Hamilton’s equations, generalized coordinates qn ¼ un , pn ¼ mu_ n and Eqs. (1) and (2) we easily obtain a following dynamical equation of motion
€ n k ðunþ1 þ un1 2un Þ qE Aun þ Bu3n þ c u_ n ¼ 0: mu
ð3Þ
This is a difference–differential equation where a dot denotes a first derivative with respect to time. The last term is a vis_ where c is a viscosity coefficient [8]. cosity force F v ¼ c u, Using a continuum approximation un ðtÞ ! uðx; tÞ and a series expansion
un1 ! u
@u 1 @2u 2 lþ l : @x 2 @x2
ð4Þ
Eq. (3) can be straightforwardly transformed into the following dynamical equation of motion
m
2 @2u @u 2@ u kl qE Au þ Bu3 þ c ¼ 0: 2 2 @x @t @t
ð5Þ
Obviously, Eq. (5) is a nonlinear partial differential equation (PDE). It is well known that PDE can be transformed into an ordinary differential equation (ODE). This is a traveling wave uðnÞ solution which depends upon x and t only through a unified variable n
n jx xt; where
ð6Þ
j and x are constants. Substitution of x and t by n straightforwardly transforms (5) into the following ODE 2
d w 2
dn
q dw 1 3 r þ ðw wÞ ¼ ; a dn a a
ð7Þ
where a dimensionless function w was introduced by the relation
u¼
rffiffiffi A w: B
ð8Þ
The expressions for the parameters existing in (7) are 2
a¼
mx2 kl A
j2
;
qE
r ¼ qffiffi ; q ¼ A
A B
cx A
:
ð9Þ
We treat the parameters q and r as input and determine w and a. A combined potential that consists of the last two terms in the Hamiltonian (1) can be written as
VðwÞ ¼
A2 f ðwÞ; B
1 1 f ðwÞ ¼ rw w2 þ w4 : 2 4
ð10Þ
S. Zdravkovic´ et al. / Applied Mathematics and Computation 242 (2014) 353–360
355
This is obviously a nonsymmetrical double well potential. Positions of its two minima (right and left) and maximum are solutions of the equation
w3 w r ¼ 0:
ð11Þ
The corresponding values are [9]
2 wR ¼ pffiffiffi cos F; 3
ð12Þ
pffiffiffi 1 wM ¼ pffiffiffi cos F þ 3 sin F ; 3
ð13Þ
pffiffiffi 1 wL ¼ pffiffiffi cos F þ 3 sin F ; 3
ð14Þ
where
F¼
1 arccos 3
r 2 ; r0 ¼ pffiffiffi : r0 3 3
ð15Þ
There are a few methods for solving Eq. (7) [8,9,14,16]. Among them, the modified extended tanh-function method is probably the simplest one. The method is based on premise that the solution of (7) is in the form of tangent hyperbolic trial function. In fact, we look for the values of certain parameters that allow the trial function to satisfy Eq. (7). In what follows, we explain how a method of factorization [17] can be used for this purpose. Essentially, this procedure is more general as no trial function is required. 3. Method of factorization One of the goals of this paper is to show how the method of factorization can be used in nonlinear dynamics of MTs, in particular for solving Eq. (7). A first step is to get rid of the constant term in the right hand side of (7). This can be done using a new function UðnÞ defined as
wðnÞ ¼ UðnÞ þ d;
ð16Þ
where d is a constant that should be determined. From (7) and (16) we easily obtain 2
d U dn2
q dU 1 þ PðUÞ ¼ 0; a dn a
ð17Þ
where
PðUÞ ¼
1
a
1
UðU2 þ 3dU þ 3d2 1Þ UðU U ÞðU Uþ Þ;
a
ð18Þ
if
d3 d r ¼ 0:
ð19Þ
According to Eqs. (11) and (19) we see that the solutions of (19) are given by Eqs. (12)–(14). Notice that
wR > 0;
wM < 0;
wL < 0;
ð20Þ
as well as
wR þ wM þ wL ¼ 0:
ð21Þ
The zeros of the polynomial PðUÞ are
U ¼
3d
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4 3d2 ; 2
Uþ ¼
3d þ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4 3d2 : 2
ð22Þ
One can easily check that these expressions for the three values for d are
pffiffiffi
UðRÞ ¼ 3 cos F sin F ¼ wL wR UðRÞ þ
UðMÞ ¼ 2 sin F ¼ wL wM UðMÞ þ
)
pffiffiffi ¼ 3 cos F þ sin F ¼ wM wR pffiffiffi ¼ 3 cos F sin F ¼ wR wM
;
ð23Þ
) ;
ð24Þ
S. Zdravkovic´ et al. / Applied Mathematics and Computation 242 (2014) 353–360
356
and
)
UðLÞ ¼ 2 sin F ¼ wM wL UðLÞ þ ¼
pffiffiffi 3 cos F þ sin F ¼ wR wL
ð25Þ
:
We are now ready for solving Eq. (17). We rewrite it in a form [17]
d d f2 ðUÞ f1 ðUÞ UðnÞ ¼ 0; dn dn
ð26Þ
where f1 ðUÞ and f2 ðUÞ are real functions that should be determined. If we compare Eqs. (26) and (17) we straightforwardly obtain
f1 ðUÞf2 ðUÞ ¼
PðUÞ
U
;
ð27Þ
df1 ðUÞ q ¼ : a dU
ð28Þ
and
f1 ðUÞ þ f2 ðUÞ þ U
Eqs. (18) and (27) indicate an appropriate choice for the functions f1 ðUÞ and f2 ðUÞ. A key question is a sign of a. For a < 0 a good choice is a ffi f1 f1 ðUÞ ¼ pffiffiffiffi ðU U Þ a
) ð29Þ
;
ffi ðU Uþ Þ f2 f2 ðUÞ ¼ ap1ffiffiffiffi a
where the parameter a should be determined. If we plug (29) into (28) we obtain
a2 ¼
1 ; 2
ð30Þ
and
1 q aU Uþ ¼ pffiffiffiffiffiffiffi : a a
ð31Þ pffiffiffi
If we had assumed a positive a, i.e. a in (29), we would have obtained an imaginary a instead of (30). Hence, we conclude that a < 0, which was obtained previously using different procedures [9,16]. According to Eqs. (30) and (31) we can conclude that a and U 2Uþ have opposite signs. Also, Eqs. (23)–(25) show that
U 2Uþ ¼ 3wM
for d ¼ wR and a ¼ p1ffiffi2
U 2Uþ ¼ 3wR for d ¼ wL ;
)
d ¼ wM and a ¼ p1ffiffi2
;
ð32Þ
which brings about 2
2q a aM ¼ 9w 2
M
2
2q a aR ¼ 9w 2
for d ¼ wL ;
R
9 > = : ; and a ¼ p1ffiffi2 >
for d ¼ wR and a ¼ p1ffiffi2 d ¼ wM
ð33Þ
Before we proceed we want to give a few comments. Note that the first and the second expressions in (32) and (33) were obtained for negative and positive values of a, respectively. Also, the minus sign could be included in the second expression in (29) instead of in the first one. However, this yields the same result. Finally, there is one more choice corresponding to (29). This is a ffi f10 f10 ðUÞ ¼ pffiffiffiffi ðU Uþ Þ a
)
ffi ðU U Þ f20 f20 ðUÞ ¼ ap1ffiffiffiffi a
ð34Þ
:
Using the procedure explained above we obtain (30) again as well as 2
2q a aM ¼ 9w 2
M
2
a aL ¼ 29wq2 L
for d ¼ wL for d ¼ wR
9 > = ; ; and d ¼ wM >
ð35Þ
pffiffiffi instead of (33). All the three values for a in Eq. (35) are obtained for a ¼ 1= 2 while the positive a does not give any solution. To obtain the function w we should solve the equation [17]
dU ¼ f1 ðUÞUðnÞ; dn
ð36Þ
S. Zdravkovic´ et al. / Applied Mathematics and Computation 242 (2014) 353–360
357
which can be seen from (26), and use Eq. (16). There are two choices for f1 , given by (29) andp(34), and three possible values ffiffiffi for d. As an example, we pick up (29) and d ¼ wR . Eqs. (32a), (33a), (34a), (36) and a ¼ 1= 2 give
dU 3w ¼ L dn: 2q UðU Uþ Þ
ð37Þ
Notice that wL < 0, while Uþ is given by Eq. (23b). For U > 0 and U < Uþ we obtain diverging solutions that, of course, do not have physical meaning. However, for Uþ < U < 0 we obtain
U¼
Uþ 3wL Uþ ðnn0 Þ 2q
1 þ e
ð38Þ
;
where n0 is a constant of integration. Finally, Eqs. (16), (23b) and (38), for n0 ¼ 0, yield
9 wR Þ expðv1 nÞ = w1 ¼ wR þ ðwM1þexpð v1 nÞ : ; v1 ¼ 3wL ðw2RqwM Þ
ð39Þ
To understand physical meaning of the function w1 we should see its asymptotic values. As
w1 ð1Þ ¼ wR ;
w1 ðþ1Þ ¼ wM ;
v1 > 0 these values are: ð40Þ
which means that the solution w1 describes a transition from the right minimum to the maximum. pffiffiffi Therefore, the solution given by (39) is a result of the combination d ¼ wR and a ¼p1= 2. Using the same procedure as ffiffiffi above we can check that the same result would be obtained for d ¼ wM and a ¼ 1= 2. The remaining four combinations give two more solutions. These are
9 wR Þ expðv2 nÞ = w2 ¼ wR þ ðwL1þexpð v nÞ 2
ð41Þ
;
v2 ¼ 3wM ðw2qR wL Þ and
9 wM Þ expðv3 nÞ = w3 ¼ wM þ ðwL 1þexpð v3 nÞ : ; v ¼ 3wR ðwM wL Þ 3
ð42Þ
2q
pffiffiffi pffiffiffi pffiffiffi For thepsolution w2 we need to assume either d ¼ wL , a ¼ 1= 2 or d ¼ wR , a ¼ 1= 2, while d ¼ wM , a ¼ 1= 2 and d ¼ wL , ffiffiffi a ¼ 1= 2 give the solution w3 . The appropriate asymptotic values are
w2 ð1Þ ¼ wR ;
w2 ðþ1Þ ¼ wL ;
ð43Þ
w3 ð1Þ ¼ wL ;
w3 ðþ1Þ ¼ wM :
ð44Þ
and
All the three solutions are shown in Fig. 1 (solid lines) for crucial in nonlinear dynamics of MTs.
r ¼ 0:85r0 and q ¼ 1. We see that kink and antikink solitons are
analyt.
1.0
numer.
ψ1 ( αL ) ψ2 (αM)
ψ3 (αR)
ψ
0.5
0.0
-0.5
-1.0
-6
-4
-2
0
2
4
6
ξ Fig. 1. Solutions w1 , w2 and w3 for
r ¼ 0:85r0 and q ¼ 1. Solid and dashed lines denote analytically and numerically derived solutions, respectively.
S. Zdravkovic´ et al. / Applied Mathematics and Computation 242 (2014) 353–360
358
We would like to point out that the same results were obtained previously using different mathematical procedures. It turned out that the modified extended tanh-function method and Jacoby elliptic functions are useful tools but the procedure explained here is the simplest one. Also, as stated above, this method is the most general one as it does not require any trial function. 4. Numerical research In order to prove the analytical findings, the numerical calculations of the model Eq. (7) are performed. The 2nd order ODE (7) is firstly transformed into the following system of the two 1st order ODEs by defining two variables YðnÞ ¼ wðnÞ and ZðnÞ ¼ dYðnÞ=dn
)
dY dn
¼Z
dZ dn
¼ qa Y a1 ðY 3 YÞ þ ra
ð45Þ
:
The next task was establishment of the initial (asymptotic [10]) conditions for both numerical variables: values of Y are chosen according to the analytical estimations given by Eqs. (40), (43), and (44), while the values of Z were chosen to be close to zero in all three cases. Finally the system (45) is solved by adopting the Runge–Kutta numerical procedure properly adjusting the step size [10]. Obtained steady state solutions for each of the initial conditions are compared with the corresponding analytical ones in Fig. 1. A perfect agreement between the analytical and the numerical solutions is obvious. Numerical procedure shows the possibility to obtain two different steady state solutions for each value of the parameter a depending on the initial system preparation. One of them fits the analytically derived solution as presented in Fig. 1. New
two solution types numerically obtained
1.0
αL αL αM αM αR αR
ψ
0.5
0.0
-0.5
-1.0
-6
-4
-2
0
2
4
6
2
4
6
ξ
(a) ψa (αM ) ψ3 (αR )
-0.4
ψ
ψb (αL)
-0.6
-0.8
-6
-4
-2
0
ξ
(b) Fig. 2. (a) Six solutions of Eq. (7) obtained numerically. (b) The enlarged part of (a).
S. Zdravkovic´ et al. / Applied Mathematics and Computation 242 (2014) 353–360
ψ1
1.6
ψ2
ψ3
359
discr. eq. (num.) cont. eq. (anal.)
ψ
0.8
0.0
-0.8 -6
-4
-2
0
ξ
2
4
6
Fig. 3. Comparative presentation of the solutions w1;2;3 derived analytically under the continuum approximation and numerically by solving the discrete system (3).
solutions are noted by dashed green line for aL , dashed violet line for aM and solid cyan line for aR , respectively. All this is shown in Fig. 2(a). The part (b) of Fig. 2 is enlarged in order to clearly distinguish different solutions. The main goal of this paper is to study either the continuum approximation can be used to study nonlinear dynamics of MTs. Validity of the continuum approximation used above is shown by numerically solving the discrete model (3). Numerical relaxation method [18] is applied to the system of N discrete nonlinear equations and obtained results are presented in Fig. 3. An excellent agreement between the curves coming from analytical and numerical investigations show that MT can be treated as the continuum system, i.e. that the continuum approximation is applicable. 5. Concluding remarks As was stated above, a key question in modeling the nonlinear dynamics of MTs is possibility to apply continuum approximation. Studying both continuum and discrete model equations of MTs we proved that this approximation can be used, which significantly simplifies research of MT systems. It was mentioned in Introduction that MT served as a network for motor proteins. To start moving a motor protein needs a signal. We postulated that such a signal is nothing but the kink soliton studied in this paper [6]. It was shown that kink solitons move, so that they are responsible for energy transfer, along MTs. However, our preliminary results indicate that breather-type solitons, that are localized modulated solitonic waves, may also play an important role in MT dynamics. This topic will be considered in detail in the forthcoming publication. It is obvious that we have relied on classical physics in this paper.It is widely accepted that the dynamics of proteins in the context of their tertiary structure is adequately described by classical laws [19]. The quantum treatment should be applied for secondary structure protein elements, primarily a-helix. This is because the peptide groups within a-helix have masses small enough for quantum mechanical treatment. However, we here deal with mesoscopic polymer with elementary mass units represented by protein tubulin dimers with the mass of about 1:1 105 atomic units (Daltons). Acknowledgments This research was supported by funds from Serbian Ministry of Education and Sciences, Grants III45010 and OI171009. References [1] P. Dustin, Microtubules, Springer, Berlin, 1984. [2] J.A. Tuszyn´sky, S. Hameroff, M.V. Sataric´, B. Trpisová, M.L.A. Nip, Ferroelectric behavior in microtubule dipole lattices: implications for information processing, signaling and assembly/disassembly, J. Theor. Biol. 174 (1995) 371–380. [3] M.V. Sataric´, J.A. Tuszyn´sky, Nonlinear dynamics of microtubules: biophysical implications, J. Biol. Phys. 31 (2005) 487–500. [4] M. Cifra, J. Pokorny´, D. Havelka, O. Kucˇera, Electric field generated by axial longitudinal vibration modes of microtubule, BioSystems 100 (2010) 122–131. [5] D. Havelka, M. Cifra, O. Kucˇera, J. Pokorny´, J. Vrba, High-frequency electric field and radiation characteristics of cellular microtubule network, J. Theor. Biol. 286 (2011) 31–40. [6] S. Zdravkovic´, M.V. Sataric´, S. Zekovic´, Nonlinear dynamics of microtubules – a longitudinal model, Europhys. Lett. 102 (2013) 38002. [7] J.E. Schoutens, Dipole–dipole interactions in microtubules, J. Biol. Phys. 31 (2005) 35–55. [8] M.V. Sataric´, J.A. Tuszyn´ski, R.B. Zˇakula, Kink like excitations as an energy-transfer mechanism in microtubules, Phys. Rev. E 48 (1993) 589–597. [9] S. Zdravkovic´, L. Kavitha, M.V. Sataric´, S. Zekovic´, J. Petrovic´, Modified extended tanh-function method and nonlinear dynamics of microtubules, Chaos Solitons Fractals 45 (2012) 1378–1386.
360
S. Zdravkovic´ et al. / Applied Mathematics and Computation 242 (2014) 353–360
[10] S. Zdravkovic´, M.V. Sataric´, A. Maluckov, A. Balazˇ, A nonlinear model of the dynamics of radial dislocations in microtubules, Appl. Math. Comput. 237 (2014) 227–237. [11] P. Drabik, S. Gusarov, A. Kovalenko, Microtubule stability studied by three-dimensional molecular theory of solvation, Biophys. J. 92 (2007) 394–403. [12] E. Nogales, M. Whittaker, R.A. Milligan, K.H. Downing, High-resolution model of the microtubule, Cell 96 (1999) 79–88. [13] M.A. Collins, A. Blumen, J.F. Currie, J. Ross, Dynamics of domain walls in ferrodistortive materials. I. Theory, Phys. Rev. B 19 (1979) 3630–3644. [14] A. Gordon, Nonlinear mechanism for proton transfer in hydrogen-bonded solids, Phys. B 146 (1987) 373–378. [15] A. Gordon, Kink dynamics in hydrogen-bounded solids, Phys. B 151 (1988) 453–461. [16] S. Zekovic´, S. Zdravkovic´, L. Kavitha, A. Muniyappan, Employment of Jacobian elliptic functions for solving problems in nonlinear dynamics of microtubules, Chin. Phys. B 23 (2014) 020504. [17] O. Cornejo-Pérez, J. Negro, L.M. Nieto, C. Rosu, Traveling-wave solutions for Korteweg-de Vries–Burgers equations through factorizations, Found. Phys. 36 (2006) 1587–1599. [18] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes in Fortan, second ed., Cambridge University Press, 2000. [19] K. Huang, Lectures on Statistical Physics and Protein Folding, World Scientific Publishing Co., 2005. Chapter 16.
