- Email: [email protected]
Coupled energy patterns in zigzag molecular chains
Accepted Manuscript Coupled energy patterns in zigzag molecular chains Bertrand Sadjo, Conrad B. Tabi, Hervais Edongue, Alidou Mohamadou PII: DOI: Reference:
S0165-2125(17)30052-5 http://dx.doi.org/10.1016/j.wavemoti.2017.04.008 WAMOT 2157
To appear in:
Wave Motion
Received date: 28 July 2016 Accepted date: 7 April 2017 Please cite this article as: B. Sadjo, C.B. Tabi, H. Edongue, A. Mohamadou, Coupled energy patterns in zigzag molecular chains, Wave Motion (2017), http://dx.doi.org/10.1016/j.wavemoti.2017.04.008 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Coupled Energy Patterns in Zigzag Molecular Chains Bertrand Sadjo1∗, Conrad B. Tabi1,2,3†, Hervais Edongue4‡, and Alidou Mohamadou5§ 1
Laboratory of Biophysics, Department of Physics, Faculty of Science, University of Yaound´e I, P.O. Box 812, Yaound´e, Cameroon.
2 3 4
Botswana International University of Science and Technology, Private Bag 16 Palapye, Botswana
The African Institute for Mathematical Sciences, 6-8 Melrose Rd, Muizenberg 7945, South Africa.
Laboratory of Material Sciences, Department of Physics, Faculty of Science, University of Yaound´e I, P.O. Box 812, Yaound´e, Cameroon.
5
Department of Physics, Faculty of Science, University of Maroua, P.O. Box 46, Maroua, Cameroon.
April 15, 2017
Abstract The contribution of both longitudinal and transversal nonlinear oscillations to energy localization is investigated in a zigzag molecular chain, which include simultaneously nearest- and next-nearest neighbor interactions. Coupled amplitude equations are found in the form of discrete nonlinear Schr¨odinger equations, whose plane wave solutions are found to be subjected to some instabilities. They are shown to be very sensitive to transverse and longitudinal couplings, which is confirmed via direct numerical simulations. The two available modes are found to be alternatively responsible for energy localization and transport. Thermal fluctuations effects bring about highly localized modes, along with narrow structures for efficient energy transport. [email protected] (B. Sadjo) Corresponding author (C. B. Tabi): [email protected]. Tel: +237-6788-735-50 ‡ [email protected] (H. Edongue) § [email protected] (A. Mohamadou)
∗
†
1
PACS number(s): 87.15.Aa, 05.45.Yv, 05.45.-a, 42.65.Tg Keywords: Zigzag chain, Modulational instability, Energy patterns.
Highlights: • The nonlinear dynamics of complex zigzag molecules is addressed. • Transverse and longitudinal dynamics are shown to both contribute to energy localization via modulational instability • Thermal fluctuations are responsible for highly localized energy patterns.
1
Introduction
The functional motion of proteins is strongly related to their structure. Usually, proteins are pictured as one-dimensional entities, generally involving only interactions among nearest-neighbors lattice sites [1, 3, 4]. During the last decades, considerable efforts have been made to show that beyond that simplified formulation of the problem, additional degrees of freedom can lead to a more suitable description of the way energy is transported and stored in biomolecules for specific purposes [1, 2]. In a few biological molecules like DNA, the obvious contribution of both longitudinal and transversal oscillations were already brought forth, showing how important their concomitant effects could be in the key processes of replication and transcription [5, 6, 7, 8]. On the other side, from the Davydov’s model [9, 10, 11], more complex models have been proposed, that show that for energy to be efficiently transported, the presence of at least three spines is important [12, 13, 14, 15, 16]. This has also been addressed recently [15, 17], where the strong impact
2
of radial and longitudinal displacements of peptides groups were confirmed to bring about new features in vibrational energy localization among hydrogen bonding spines. Albeit these considerable efforts to suitably capture and understand the energy transport phenomenon in biomolecules, it remains clear that much work still has to be done, which considers other aspects of their geometry, specifically, the mutual effect of their first and secondary structures. According to Ref. [18, 19], the simplest way to achieve this is to take into consideration first- and second-neighbor interactions, standing respectively for the deformation of valence bonds and the deformation of valence and torsion angles. In describing the interactions between first neighbors as transverse and those between second neighbors as longitudinal, one might be capable of describing the dynamics of more complex biological polymers by adopting their simplest two-dimensional structure in the form of zigzag chains [20], as represented in Fig. 1. The zigzag configuration is therefore a more generalized description of a broad range of molecular entities, whose the dynamics involves more than one degree of freedom, a good example being α−helix proteins, which include both rigid valence bonds (first-neighbor interactions) and soft hydrogen bonds (second-neighbor interactions). Compression soliton solutions, their stability and condition of existence were studied in Refs. [18, 19]. Indubitably, the existence of solitons implies the possibility of energy localization with reference to the pioneering work of Fermi, Pasta and Ulam [21]. One of the direct mechanisms leading to this phenomenon is modulational instability (MI), which to the best of our knowledge has not yet been used in the context of molecular zigzag chains. This constitutes the main objective of the present paper. MI is in fact the a process whereby a plane wave, under the competitive effects of dispersion and nonlinearity, breaks up into trains of soliton-like structures [22, 23, 24, 25, 26, 27, 28, 29]. In several occasions, it has also been shown to be responsible for the localization of energy in some biological lattice such as DNA [8, 30] and proteins [3, 15, 16, 17], with the obvious contribution of their structure and the involved interactions. In the rest of this work, we first present the zigzag model and we show that its transversal and longitudinal dynamics can be fully represented by
3
Figure 1: Schematic representation of a zigzag-shaped molecular chain, where U and V represent the intermolecular interaction potentials between first and second neighbors, respectively. xn and yn are respectively the longitudinal and transversal displacements, K1 and K2 are the stiffness constants for the first and second-neighbor forces. coupled discrete nonlinear Schr¨odinger (DNLS) equations, involving first and second dispersive and nonlinear terms. The linear stability analysis is thereafter performed and the existence of unstable patterns of energy is discussed. We also present some numerical applications of our predictions and discuss the features of the subsequent energy patterns, along with their behaviors under thermal fluctuations. Some concluding remarks are finally given.
2
Model and coupled amplitude equations
The planar zigzag model addressed in this work comprises two degrees of freedom xn (t) and yn (t) which represent the longitudinal and transversal displacements of the lattice site n with respect to their equilibrium positions. The total Hamiltonian that describes the dynamics of the whole lattice is written as H=
X n
(
1 M 2
"
dxn dt
2
+
dyn dt
2 #
)
+ Kl2 [U (rn ) + V (qn )] + ε0 Z(un , vn ) ,
4
(1)
where the first term stands for the kinetic energy, with M being the reduced mass of a chain particle. K is the stiffness constant, while intermolecular interactions between first and second neighbors are respectively described by the dimensionless potentials U (rn ) and V (qn ), and where the deviations from the equilibrium length rn and qn are such that s 2 2 1 xn+1 − xn yn + yn+1 + + h− −b rn = 2 l l s (2) 2 2 xn+1 − xn−1 yn+1 − yn−1 qn = 1+ + − 1. l l q b = h2 + 14 is the dimensionless equilibrium length of the rigid bonds and h is the thickness of the zigzag backbone, in units of the longitudinal lattice spacing l. The potentials U (r) and V (q) are normalized by U (0) = V (0) = 0 and U ”(0) = κ1 , V ”(0) = κ2 , where κ1 = K1 /K and κ2 = K2 /K are the dimensionless stiffness constants for the first- and second-neighbors forces, respectively. While ε0 describes the height of the potential Z(u, v), the latter brings about the zigzag features of the chain. The potential Z(u, v) is in fact made of two periodic functions f (u) and g(u), so that 0 ≤ f ≤ 1 and 0 ≤ g ≤ h. These two functions also have the periodicity properties f u + 21 = f (u)
and g(u + 1) = g(u), and their more explicit forms, as proposed in Ref. [18], are f (u) = sin2 (2πu) and g(u) = h sin2 (πu).
(3)
We assume f 00 ( n2 ) = Ω20 , with Ω0 being the dimensionless characteristic frequency of small oscillations of particles at the minima of the substrate potential Z(u, v), whose the complete expression is proposed in the form [18] 2 1 Z(u, v) = f (u) + Ω0 2 v − g(u) . 2
(4)
The features of the above potential are depicted in Fig. 2, versus the displacements u and v, for different values of the frequency Ω0 . In the vertices of the zigzag structure of the potential, the intermolecular bonds remain undistorted. The structure of this configuration does not change 5
Figure 2: Plots of the 2D potential versus u and v for h = 0.5 and: (a1)-(a2) Ω20 = 10, (b1)-(b2) Ω20 = 15 and (c1)-(c2) Ω20 = 20. with increasing Ω0 , but the limits between them and the maxima appear clearly. Furthermore, we q K introduce the dimensionless parameters t → M t, un (t) → xnl(t) , vn (t) → ynl(t) and we consider
the simplest harmonic form for the potentials U and V , such that U (r) = r2 /2 and V (q) = q 2 /2. Considering all this, Hamilton’s equations can be applied to the dimensionless version of (1), leading to the following coupled dynamical equations in un (t) and vn (t): 1 1 + un+1 − un + un − un−1 d2 un 2 2 =κ1 un+1 − 2un + un−1 − − dt2 1 + rbn 1 + rn−1 b 1 + un+2 − un 1 + un − un−2 ∂ − −ϑ Z(un , vn ) + κ2 un+2 − 2un + un−2 − 1 + qn+1 1 + qn−1 ∂un d2 v n h − vn−1 − vn h − vn − vn+1 =κ1 2h − vn−1 − 2vn − vn+1 − + dt2 1 + rn−1 1 + rn b b vn+2 − vn vn − vn−2 ∂ + κ2 vn+2 − 2vn + vn−2 − − −ϑ Z(un , vn ), 1 + qn+1 1 + qn−1 ∂vn
6
(5a)
(5b)
with ϑ =
ε0 . Kl2
After setting ϑ = 1 and expanding
∂Z ∂un
and
∂Z ∂vn
up to the third order, it is possible to
reduce Eqs.(5) to their amplitude equations. In performing such an expansion, one should keep in mind that the system should keep the balance between nonlinear and dispersive effects. Therefore, one of the techniques that usually work is the Fourrier series expansion method which assumes un and vn to have the trial forms [7, 8] un (t) =
+∞ X
ipωb t a(p) and vn (t) = n e
p=−∞
+∞ X
ipωb t b(p) , n e
(6)
p=−∞
where ωb is close to some linear oscillation frequencies. The Fourier coefficients are slowly varying (p)
with time, as they are functions of τ = 2 t, and they exponentially decay in p so that an ∼ p ∗ (p) (0) (0) (p) (−p) and bn ∼ p for p > 0, with an ∼ 2 (bn ∼ 2 , respectively). Moreover, an = an and ∗ (p) (−p) bn = bn . By applying trial solutions (6) to the expanded Eqs.(5), and considering the (1)
harmonics p = 1, one obtains a set of coupled equations whose time-dependent solutions an (t) and (1)
bn (t) are assumed to be in the form r 2ωb 8π 2 − ωb2 + α0 (3 + 4h2 ) + 2(κ2 − κ1 ) (1) an (t) = νn (t) exp i t p0 + 4α0 2ωb r 2ωb 2α0 (1 + 2h2 ) − 2κ1 − ωb2 + Ω0 2 (1) n bn (t) = (−1) ψn (t) exp i t , κ1 2ωb 4
where p0 = 2Ω0 2 π 4 h2 − 64π and α0 = 3
κ1 . 2b
(7a)
(7b)
Further mathematical expansions, using (7a), lead to the
following set of nonlinearly coupled equations iν˙n − P1 (νn+1 + νn−1 ) − P2 (νn+2 + νn−2 ) p0 |νn |2 νn + α0 (|νn+1 − νn |2 (νn+1 − νn ) − |νn − νn−1 |2 (νn − νn−1 )) − Q1 − κ (|ν 2 2 2 n+2 − νn | (νn+2 − νn ) − |νn − νn−2 | (νn − νn−2 )) + Q2 (|ψn − ψn+1 |2 (νn+1 − νn ) − |ψn − ψn−1 |2 (νn − νn−1 ))
− Q3 (|ψn+2 − ψn |2 (νn+2 − νn ) − |ψn − ψn−2 |2 (νn − νn−2 )) = 0
7
(8a)
iψ˙n − P3 (ψn+1 + ψn−1 ) + Q1
α0 (|ψn − ψn−1 |2 (ψn − ψn−1 ) + |ψn − ψn+1 |2 (ψn − ψn+1 ))
− κ (|ψ 2 2 2 n+2 − ψn | (ψn+2 − ψn ) − |ψn − ψn−2 | (ψn − ψn−2 ))
+ Q2 (|νn − νn−1 |2 (ψn − ψn−1 ) + |νn+1 − νn |2 (ψn − ψn+1 ))
(8b)
− Q4 (|νn+2 − νn |2 (ψn+2 − ψn ) − |νn − νn−2 |2 (ψn − ψn−2 )) = 0, with the different coefficients P1 = and P3 =
(1+h2 )α2 −κ1 , 2ωb
Q4 =
κ2 Q1 . 2
α0 (3+4h2 ) , 2ωb
P2 =
κ2 , 2ωb
Q1 = (p0 + 4α0 )−1 , Q2 = α0 /κ2 , Q3 =
κ2 , 4κ1
The two equations (8) are coupled DNLS equations which include
not only interactions between first neighbors, but also between second neighbors. Eq.(8a) describes the amplitude of longitudinal vibrations, while Eq.(8b) gives information regarding the transverse nonlinear oscillations. The two equations are nonlinearly coupled via the coefficients Q3 and Q4 which are functions of κ2 . Otherwise, they will completely be decoupled if κ2 → 0. Therefore, as initially described, each of the equations will consider only interactions between first neighbors. A single DNLS equation was already obtained in the context of charge transport in DNA [31] and coupled DNLS equations derived to describe nonlinear oscillations of DNA strands [7, 8, 32]. Their exact solutions were found in Ref. [7] to picture the way helicity and radial vibrations contribute to the unzipping and opening of the DNA duplex for the genetic code to be read. In the context of the present work, we intend to use the ones found here to perform the linear stability analysis of their plane wave solutions, in order to find regions where they are unstable, or more precisely, where solitonic and nonlinear energy patterns can be observed. However, direct numerical verifications will be performed on the generic Eqs. (5).
3
Linear stability analysis
The initial step to detect regions of parameters where modulated waves can be observed is to assume plane wave solutions for the system equations. For the specific case at hand, such solutions are supposed to be in the form νn = ν0 ei(k1 n−ω1 t) and ψn = ψ0 ei(k2 n−ω2 t) , where the frequencies ωj and
8
the wavenumbers qj , (j = 1, 2), satisfy the dispersion relations k1 4 4 2 ω1 = 2P1 cos(k1 ) + 2P2 cos(2k1 ) + Q1 |ν0 | p0 + 16κ2 sin (k1 ) − 16α0 sin 2 k1 k2 sin2 − Q3 sin2 (k1 ) sin2 (k2 ) + 16|ψ0 |2 Q2 sin2 2 2 k2 2 4 4 − k2 sin (k2 ) ω2 = 2P3 cos(k2 ) − 16Q1 |ψ0 | α0 sin 2 k2 k1 2 2 2 2 2 sin + Q4 sin (k1 ) sin (k2 ) . − 16|ν0 | Q2 sin 2 2
(9)
To further proceed, we assume ω1 = ω2 = ω and k1 = k2 = k, and one can get a relation between the two amplitudes ψ0 and ν0 as follows: (P3 − P1 ) cos(k) − P2 cos(2k) 8 (Q2 + α0 Q1 ) sin4 k2 (Q1 κ2 − Q3 ) sin4 (k) 16(α0 Q1 − Q2 ) sin4 k2 − Q1 p0 − 16(κ2 Q1 + Q4 ) sin4 (k) + |ν0 |2 . 16 (Q2 + α0 Q1 ) sin4 k2 + (Q1 κ2 − Q3 ) sin4 (k)
|ψ0 |2 =
(10)
To study the stability of the previous trial solutions, their behavior under small perturbations should be observed. These solutions can therefore be perturbed in amplitude or in phase, or both of them. However, when we consider the amplitude perturbations, i.e., νn = ν0 [1 + ηn (t)]ei(kn−ωt) and ψn = ψ0 [1 + ξn (t)]ei(kn−ωt) , we find the perturbations ηn and ξn to be governed by the coupled
9
set of linearized equations h i iη˙n + ωηn − P1 (ηn+1 + ηn−1 ) cos(k) + i(ηn+1 − ηn−1 ) sin(k) h i − P2 (ηn+2 + ηn−2 ) cos(2k) + i(ηn+2 − ηn−2 ) sin(2k) ∗ ∗ p0 (2ηn + ηn∗ ) + α0 [2(4 cos(k) − 3)ηn + (2 cos k − 2)(ηn+1 + ηn−1 + ηn∗ ) + (4 cos(k) − 2)(ηn+1 + ηn−1 ) + 4i(ηn+1 − ηn−1 ) sin(k)] − κ2 [(8 cos 2k − 6)ηn − Q1 |ν0 |2 ∗ ∗ + (2 cos(2k) − 2)(ηn+2 + ηn−2 + ηn∗ ) + (4 cos(2k) − 2)(ηn+2 + ηn−2 ) + 4i(ηn+2 − ηn−2 ) sin(2k)] ∗ ∗ 2 cos 2kξ + (4 cos(k) − 2)(η + ξ ) n n n ∗ ∗ 2 + 2(cos(k) − 1)(2ξn + ξn+1 + ξn−1 ) + (2 cos(k) − 1)(ξn+1 + ξn−1 ) + Q2 |ν0 | + 2(cos(k) − 1)(ηn+1 + ηn−1 ) + 2i(ηn+1 − ηn−1 ) sin k + 2i(ξn+1 − ξn−1 ) sin(k) ∗ ∗ 2 cos 4kξn∗ + (4 cos(2k) − 2)(ηn + ξn∗ ) + (2 cos(2k) − 2)(2ξn + ξn+2 + ξn−2 ) 2 + (2 cos(2k) − 1)(ξn+2 + ξn−2 ) + (2 cos(2k) − 2)(ηn+2 + ηn−2 ) =0 − Q3 |ψ0 | + 2i(ηn+2 − ηn−2 ) sin(2k) + 2i(ξn+2 − ξn−2 ) sin(2k) (11a)
10
h i ˙ iξn + ωξn − P3 (ξn+1 + ξn−1 ) cos(k) + i(ξn+1 − ξn−1 )sin(k)) ∗ ∗ α0 [(ξn+1 + ξn−1 )(2 − 2 cos k) − 4i(ξn+1 − ξn−1 ) sin(k) + (ξn+1 + ξn−1 )(2 − 2 cos k) ∗ ∗ + 8(1 − cos(k))ξ + (2 − 4 cos k)ξ + 2 cos 2kξ + 2(ξ + ξ ) cos 2k n n+1 n−1 n n 2 + 2i(ξn+1 − ξn−1 ) sin 2k] + κ2 [(ξn+2 + ξn−2 )(−2 + 2 cos(2k)) + 4i(ξn+2 − ξn−2 ) sin(2k) − Q1|ψ0 | ∗ ∗ ∗ + (ξn+2 + ξn−2 )(−2 + 2 cos 2k) + 8(−1 + cos(2k))ξn + (−2 + 4 cos(2k))ξn ∗ − 2 cos(4k)ξn − 2(ξn+2 + ξn−2 ) cos(4k) − 2i(ξn+2 − ξn−2 ) sin(4k)] (ξ + ξ )(1 − 2 cos(k)) + 2i(ξ − ξ ) sin(k) + (η + η )(1 − 2 cos(k)) n+1 n−1 n+1 n−1 n+1 n−1 ∗ ∗ − 2i(ηn+1 − ηn−1 ) sin(k) + (ηn+1 + ηn−1 )(2 − 2 cos(k)) + 4(ξn + ηn )(1 − cos k) 2 − Q2 |ν0 | + 2(1 − 2 cos k)ηn∗ − 2 cos(2k)ηn∗ + (ξn+1 + ξn−1 ) cos(2k) + i(ξn+1 − ξn−1 ) sin(2k) + (ηn+1 + ηn−1 ) cos(2k) + i(ηn+1 − ηn−1 ) sin(2k) (ξn+2 + ξn−2 )(−1 + 2 cos 2k) + 2i(ξn+2 − ξn−2 ) sin 2k + (ηn+2 + ηn−2 )(−1 + 2 cos 2k) ∗ ∗ + 2i(ηn+2 − ηn−2 ) sin 2k + (ηn+2 + ηn−2 )(2 cos 2k − 2) + 4(ξn + ηn )(−1 + cos 2k) =0 − Q4 |ν0 |2 ∗ ∗ + 2(−1 + 2 cos 2k)η + 2 cos 4kη − (ξ + ξ ) cos 4k − i(ξ − ξ ) sin 4k n+2 n−2 n+2 n−2 n n − (ηn+2 + ηn−2 ) cos 4k − i(ηn+2 − ηn−2 ) sin 4k (11b) ∗ t)
whose solutions can be suggested in the form ηn (t) = B1 ei(Kn−Ωt) + B2∗ e−i(Kn−Ω
and ξn (t) =
∗
B3 ei(Kn−Ωt) +B4∗ e−i(Kn−Ω t) , with K and Ω being the wavenumber and frequency of the perturbation, respectively. Reporting this into Eqs.(11) and making use of the dispersion relations (9), one obtains a set of four homogeneous equations for B1 , B2 , B3 and B4 , i.e., M × (B1 , B2 , B3 , B4 )T = 0, where the system matrix is given by a11 + Ω a12 a13 a14 a21 a22 − Ω a23 a24 M = a31 a32 a33 + Ω a34 a41 a42 a43 a44 − Ω 11
(12)
(13)
Figure 3: MI gains along with their corresponding stability/instability diagrams, versus the wavenumbers k and K for Ω20 = 10 and: (a1)-(a2) κ2 = 0.005κ1 , (b1)-(b2) κ2 = 0.01κ1 and (c1)-(c2) κ2 = 0.5κ1 , with κ1 = 0.85. with the matrix elements given in the Appendix. The above system admits non-trivial solutions given that the determinant of its matrix is null. This results into a fourth-order polynomial in Ω, whose solutions are either real or complex. However, we are interested in complex solutions Ω = Ωr + iΩi . Nevertheless, finding then is quite long and cumbersome, which requires some symbolic computations in order to extract such solutions whose the imaginary part gives the MI gain. In fact, the latter is the largest value among those corresponding to the various branches of the dispersion relation. In what follows, the gain is computed through the formula G(k, K) = |Im(Ω)|. The corresponding results are displayed in Figs. 3 and 4. The value of κ1 has been fixed and one plays both on the values of the frequency Ω0 and the force constant κ2 between second neighbors. 12
Figure 4: MI gains along with their corresponding stability/instability diagrams, versus the wavenumbers k and K for Ω20 = 15 and: (a1)-(a2) κ2 = 0.005κ1 , (b1)-(b2) κ2 = 0.01κ1 and (c1)-(c2) κ2 = 0.5κ1 , with κ1 = 0.85.
13
For example, in Fig. 3, we have fixed Ω20 = 10, and the gain shows only one region of instability for κ2 = 0.005κ1 . In this case, one clearly understands that only strong effects from the firstneighbor coupling contribute to the emergence of the instability zone. When κ2 increases to 0.01κ1 (Fig. 3(b1)-(b2)), the MI gain clearly shows the presence of another instability zone, while the primary one gets reduced. Following that, it therefore appears that increasing the value of κ2 is responsible for the explosion of the instability domain, as this is indeed confirmed in Figs. 3(c1)(c2), where Ω20 = 15. Obviously, from the instability diagram, the primary zone of instability has reduced, while there appear other two regions of instability. To remind, in regions of instability, the plane wave is supposed to be disintegrated into trains of solitary waveforms, or it is said to be stable when parameter are taken from the stability region. In Fig. 4, we have fixed Ω20 = 15, and we equally play on the value of κ2 , which appears to importantly influence the instability features. Contrarily o the previous case, for κ2 = 0.005κ1 (Fig. 4 (a1)-(a2)), there are three regions of instability. Again, this case shows the strong impact of the first-neighbor coupling, where the competition between the two directions of vibration is minimized. However, increasing κ2 to 0.01κ1 , the merged instability region is delocalized, therfore completely changing the region of parameters where modulated patterns are expected. Finally, when κ2 = 0.5κ1 , there is only one large region of instability as clearly shown by Fig. 4 (c1)-(c2). The behaviors induced by large and small values of Ω0 are a bit different, except that the two cases gives rise to instability regions. Also, there is an overlap of many instability points, which show that to some extent, modulated waves and nonlinear patterns of energy might be observed, with almost the same behaviors depending on the strength of the two coupling parameters. On the side, this cannot be confirmed by this analysis, which is just a procedure that gives information on the onset of MI.
14
Figure 5: Panel (a) shows the energy density as a function of time (increasing from bottom to top) and space (increasing from left to right) for κ2 = 0.005κ1 , k = 0.4π, K = 0.25π. (b1)-(b2) show the spatial expansion of longitudinal and transversal manifestations of MI at the respective instants t = 200 and t = 800. The time series for un and vn are also shown at the lattice points n = 50 and n = 150 (panels (c1) and (c2)).
4
Patterns of energy and coupled wave effects
In this section, we perform direct numerical simulations on the model equations (5) using the fourth-order Rung-Kutta computational scheme with periodic boundary conditions and time-step ∆t = 10−4 . We adopt the initial conditions un (t = 0) = u0 [1 + ζ cos(Kn)] cos(kn), u˙ n (t = 0) = u0 ω[1 + ζ cos(Kn)] sin(kn) (14) vn (t = 0) = v0 [1 + ζ cos(Kn)] cos(kn), v˙ n (t = 0) = v0 ω[1 + ζ cos(Kn)] sin(kn), with ζ being the perturbation strength fixed as ζ = 10−2 . Numerical simulations are performed on a lattice of 200 points and one extracts the information related to the contribution of both the longitudinal and transversal vibrations to the localization of energy. For this purpose, the energy is computed via its dimensionless expression " # X 1 dun 2 dvn 2 κ1 κ1 E= + + rn2 + qn2 + Z(un , vn ). 2 dt dt 2 2 n 15
(15)
For the wavenumbers k = 0.4π and K = 0.25π, with u0 = 0.5 and v0 = 0, an example of such simulations is given by Fig. 5, where we have fixed κ2 = 0.1κ1 . In Fig. 5 (a), we observe some coherent structures due to the activation of MI. The corresponding patterns are localized at specific sequences of lattice points and form rows over the time. This comportment is known, as it has been reported in numbers of models related to DNA [8, 32] and α−helix dynamics [16, 17]. However, what makes them special here is that they are the consequence of the competition between longitudinal and transversal oscillations. The value of κ2 , which is smaller than κ1 is expected not to really affect the pattern of energy, but at the same time, in Figs. 5 (b1) and (b2), one can notice that oscillations, from the two mode of vibrations, display train of waves, although transverse ones are more localized for t = 200. However, for t = 800, both of them almost display the same amplitude, which show that over long time, they all contribute to energy localization. The latter is also confirmed by Figs. 5 (c1) and (c2), where the time series of the solitonic structure are depicted at positions n = 50 and n = 150 and where the emerging phenomena is the out-of-phase oscillations of the two modes. This suggest that, each of the modes contributes to maintain energy in the system, making it available not only for specific metabolisms, but also to entertain the intrinsic dynamics of the molecule, as recently reported in the case of complex proteins with interspine coupling [15, 17, 33]. It was in fact concluded that, when two or more spines are coupled via hydrogen bonds, two of them might serve as energy reservoir, while the proper energy transport is ensured by the other. The same is also true for microtubules, as addressed in Ref. [34], where energy is always available to maintain intrinsically the molecule. Importantly, when κ2 increases to 0.5κ1 , the energy density increases and rows of energy get more and more restricted to small sequences of lattice points, as shown in Fig. 6(a). The main reason for this result is the features of transversal and longitudinal oscillations displayed in Fig. 6 (b1) and (b2), where trains of solitons are observed in both cases. Nevertheless, as only un has been excited, the behavior that appears in Fig. 6 (b1), at t = 200, seems to be the one expected. However, over long time evolution, i.e., t = 800 (Fig. 6 (b2)), transverse oscillations
16
Figure 6: Energy patterns (panel (a)) and solitonic structures due to MI (panels (b1)-(b2)) for κ2 = 0.5κ1 , k = 0.4π, K = 0.25π. In (a) time goes from bottom to top and space increases from left to right. Panels (c1) and (c2) show the time series of un and vn at the points n = 50 and n = 150. take advantage and contribute apparently to effectively localize energy. The time series of these waves still maintain the fact that the dynamics of the two modes remains out-of-phase, but one notices the increase in temporal oscillations at both n = 50 and n = 150. Also, the amplitude of such waves is not uniform as there are periods of high-amplitude oscillations. Nevertheless, the conclusion on the contribution of the two modes to energy localization remains the same as in Fig. 5. Energy patterns in proteins, and in biomolecules in general, are very sensitive to thermal fluctuations [16, 35, 36, 37]. This should however be consider, not only to evaluate their effects, but also because their inherent to biological media. For example, the fluctuational opening of DNA involves thermal fluctuations, and this enhances the creation of the bubbles that expose the bases out of the stack. As a whole, such mechanisms involves energy, which is normally brought by the hydrolysis of ATP, and conveyed to specific sites by specific enzymes such as RNA-polymerase. The latter surfs the long DNA molecule in search of the initiation sequence, and invests that energy for the replications process to be initiated [8, 32, 36]. It is evidenced, from a recent contribution, that 17
Figure 7: Noise effect on MI wave patterns for κ2 = 0.005κ1 , k = 0.4π, K = 0.25π. The energy density is represented with time increasing from bottom to top and space increasing from left to right (panel (a)). Panels (b1) and (b2) picture the spatial evolution of un and vn at instants t = 200 and t = 800, respectively. Panels (c1) and (c2) show the temporal evolution of un and vn at the lattice points n = 50 and n = 150. thermal fluctuations and dissipation can enhance high energy localization [8, 16, 32], resulting to the change in the vibration modes of the molecules. In the present case, to confirm that, we apply the random force Γn (t), with the statistical properties hΓn (t)i = 0 and hΓn (t)Γ0n0 (t)i = 2T γδn,n0 δ(t − t0 ),
(16)
to the two equations (5), where T is the temperature and γ the dimensionless friction coefficient taken to be equal to 10−2 . With in mind the study of the same effects as previously, we make numerical calculations at the temperature T = 300K, and we fix κ2 = 0.005κ1 . The corresponding results are depicted in Fig. 7, where for weak couplings, the coherency of the oscillations is destroyed, due partly to the increase in frequency of the oscillations. However, from Fig. 7(a), there are regions where energy is highly localized, due to the activation of MI. The time series shown in Figs. 7 (c1) and (c2) display complex oscillations at positions n = 50 and n = 150. Again, the amplitudes of 18
Figure 8: Noise effect on MI wave patterns for κ2 = 0.5κ1 , k = 0.4π, K = 0.25π. The energy density is plotted in panel (a) versus time (increasing from bottom to top) and space (increasing from left to right). Panels (b1) and (b2) show un and vn versus the lattice index n, at instants t = 200 and t = 800, respectively. The two panels (c1) and (c2) depict the time series of un and vn at positions n = 50 and n = 150. oscillations at those two positions are different, but with the same waveform. These features change for κ2 = 0.5κ1 , where the energy is more localized than in Fig. 6(a). There is only one pulse of energy, which has a solitonic shape (Fig. 8(a)), probably due to MI mechanism. This is the case of molecules whose energy remains concentrated in specific regions. That energy is due to the pulse solitons observed in Fig. 8(b1) and (b2) at different instants. However, the time series of the lattice point n = 50 shows that transverse and longitudinal oscillations take place, even where energy is not highly localized, which of course confirms our arguments on the permanent availability of energy in biological molecules. At position n = 150, one also confirms that the two modes display the same dynamics, except that their amplitudes are again not equal, while it seems noisy with slight modulations. In the case of α−helix proteins interactions with second neighbors are related to the deformation of the hydrogen bonds, and their highly localized oscillations, as shown in Figs. 7(c2) are in agreement with the observation made for a set of three coupled spines. Therefore, under 19
thermal fluctuations, they are favorable to the formation of breathers and to important energy outcomes.
5
Conclusion
Energy localization in zigzag molecular chains has been addressed in this paper. After reducing the whole dynamics of the complex structure to a set nonlinearly coupled DNLS equations, MI of plane wave solutions has been addressed and regions of instability have been found to be very sensitive to the coupling between longitudinal and transversal dynamics. This has been confirmed by direct numerical simulations on Eqs.(5), with specific interest in the way the subsequent nonlinear oscillations taking place in the two modes contribute to energy localization and transport in the chain. We have observed that in the absence of thermal fluctuations, such a dynamics is outof-phase, and we have evidenced that the two modes alternatively entertain permanent energy availability, not only for specific metabolisms to be achieved, but also for the intrinsic dynamics of the molecule. However, under thermal fluctuations, energy has been found to be highly localized due to the collective dynamics of the two modes. Importantly, under strong longitudinal coupling, energy has been found to display soliton-like behaviors following the individual dynamics of the two modes.
20
Appendix: Elements of the matrix (13) a11 = 2P1 cos(k) + 2(Q1 κ2 (|ψ0 |2 + 4|ν0 |2 ) + (4Q3 + Q2 )|ψ0 |2 + Q4 |ν0 |2 ) cos(2k) − 2(Q1 κ2 |ν0 |2 + Q3 |ψ0 |2 ) cos(4k) + 2|ψ0 |2 (Q2 − Q3 ) − Q1 |ν0 |2 (6κ2 + 6α0 − p0 ) − 4(Q3 − Q2 )|ψ0 |2 − 2P1 cos(K + k) − P2 cos(2(K + k)) − 2p0 Q1 |ν0 |2 h i 2 − α0 Q1 |ν0 | 2(4 cos(k) − 3 − 2 cos(K)) + 8 cos(K + k) h i − Q2 |ψ0 |2 2(2 cos(k) − 1 − 2 cos(K)) + 4 cos(K + k) h i 2 − Q1 κ2 |ν0 | 2(4 cos(2k) − 3 − 2 cos(2K)) + 8 cos(2(K + k)) h i − Q3 |ψ0 |2 (4 cos(2k) − 2 − 4 cos(2K)) + 4 cos(2(K + k))
a12 = −|ν0 |2 [p0 Q1 − α0 Q1 |ν0 |2 (2 cos(k) − 2)(2 cos(K) + 1) + κ2 Q1 (2 cos(2k) − 2)(2 cos(2K) + 1)] h i a13 = −2Q2 |ψ0 |2 (2 cos(k) − 2 − cos(K)) + 2 cos(K + k) h i 2 − Q3 |ψ0 | 2(2 cos(2k) − 2 − cos(2K)) + 2 cos(2(K + k)) h i a14 = −Q2 |ψ0 |2 2 cos(2k) + 4 cos(k) − 2 + (4 cos(k) − 4) cos(K) h i 2 − Q3 |ψ0 | 2 cos(4k) + 4 cos(2k) − 2 + (4 cos(2k) − 4) cos(2K) , a21 = a12 h i a24 = −2Q2 |ψ0 |2 (2 cos(k) − 2 − cos(K)) + 2 cos(K − k) h i 2 − Q4 |ψ0 | 2(2 cos(2k) − 2 − cos(2K)) + 2 cos(2(K − k))
21
a22 = 2P1 cos(k) + 2(Q1 κ2 (|ψ0 |2 + 4|ν0 |2 ) + (4Q3 + Q2 )|ψ0 |2 + Q4 |ν0 |2 ) cos(2k) − 2(Q1 κ2 |ν0 |2 + Q3 |ψ0 |2 ) cos(4k) + 2|ψ0 |2 (Q2 − Q3 ) − Q1 |ν0 |2 (6κ2 + 6α0 − p0 ) − 4(Q3 − Q2 )|ψ0 |2 − 2P1 cos(K − k) − P2 cos(2(K − k)) − 2p0 Q1 |ν0 |2 h i − α0 Q1 |ν0 |2 2(4 cos(k) − 3 − 2 cos(K)) + 8 cos(K − k) h i − Q2 |ψ0 |2 2(2 cos(k) − 1 − 2 cos(K)) + 4 cos(K − k) h i − Q1 κ2 |ν0 |2 2(4 cos(2k) − 3 − 2 cos(2K)) + 8 cos(2(K − k)) h i − Q3 |ψ0 |2 (4 cos(2k) − 2 − 4 cos(2K)) + 4 cos(2(K − k))
a31 = −Q2 |ν0 |2 [2 cos(K) − 4(cos(k) − 1) − 4 cos(K + k) + 2 cos(K + 2k)] − Q4 |ν0 |2 [4(cos(k) − 1) − 2 cos(K) + 4 cos(2K + k) − 2 cos(2(K + 2k))] a32 = −Q2 |ν0 |2 [(4 − 4 cos(k)) cos(K) + 2(1 − 2 cos(k)) − 2 cos(2k)] − Q4 |ν0 |2 [(−4 + 4 cos(2k)) cos(2K) + 2(−1 + 2 cos(2k)) + 2 cos(4k)] a34 = −α0 Q1 ψ02 [(4 − 4 cos(k) cos(K) + 2(1 − 2 cos(k)) + 2 cos(2k)] − κ2 Q1 |ψ0 |2 [(−4 + 4 cos(k)) cos(2K) + 2(−1 + 2 cos(2k)) − 2 cos(4k)] a23 = a14 ,
a43 = a34 ,
a41 = a32
a42 = −Q2 |ν0 |2 [4(cos(2k) − 1) − 2 cos(2K) − 4 cos(K + k) + 2 cos((K − 2k))] − Q4 |ν0 |2 [4(cos(2k) − 1) − 2 cos(2K) − 4 cos(K − k) − 2 cos(2(K − 2k))]
22
a33 = 2P1 cos(k) + 2(Q1 κ2 (|ψ0 |2 + 4|ν0 |2 ) + (4Q3 + Q2 )|ψ0 |2 + Q4 |ν0 |2 ) cos(2k) − 2(Q1 κ2 |ν0 |2 + Q3 |ψ0 |2 ) cos(4k) + 2|ψ0 |2 (Q2 − Q3 ) − Q1 |ν0 |2 (6κ2 + 6α0 − p0 ) − 4(Q3 − Q2 )|ψ0 |2 − 2P3 cos(K + k) − Q2 |ν0 |2 [2(2 cos(K) − 2 cos(k)) − 4 cos(K + k) + 2 cos(K + 2k)] − α0 Q1 |ψ0 |2 [4 cos(K) + 8(1 − cosk) − 8 cos(K + k) + 4 cos(K + 2k)] − Q4 |ν0 |2 [4(cos(2k) − 1) − 2 cos(2K) + 4 cos(2(K + k)) − 2 cos(2(K + 2k))] − κ2 Q1 |ψ0 |2 [8(cos(2k) − 1) − 4 cos(K) + 8 cos(2(K + k)) − 4 cos(2(K + 2k))] a44 = 2P1 cos(k) + 2(Q1 κ2 (|ψ0 |2 + 4|ν0 |2 ) + (4Q3 + Q2 )|ψ0 |2 + Q4 |ν0 |2 ) cos(2k) − 2(Q1 κ2 |ν0 |2 + Q3 |ψ0 |2 ) cos(4k) + 2|ψ0 |2 (Q2 − Q3 ) − Q1 |ν0 |2 (6κ2 + 6α0 − p0 ) − 4(Q3 − Q2 )|ψ0 |2 − 2P3 cos(K − k) − Q2 |ν0 |2 [2(2 cos(K) − 2 cos(k)) − 4 cos(K − k) + 2 cos(K − 2k)] − α0 Q1 |ψ0 |2 [4 cos(K) + 8(1 − cosk) − 8 cos(K − k) + 4 cos(K − 2k)] − Q4 |ν0 |2 [4(cos(2k) − 1) − 2 cos(2K) + 4 cos(2(K − k)) − 2 cos(2(K − 2k))] − κ2 Q1 |ψ0 |2 [8(cos(2k) − 1) − 4 cos(K) + 8 cos(2(K − k)) − 4 cos(2(K − 2k))]
23
References [1] A. C. Scott, Phys. Rep. 217 (1992) 1. [2] A. C. Scott, Phys. Rev. A 26 (1982) 578. [3] R.Y. Ondoua, C. B. Tabi, H. P. Ekobena, A. Mohamadou and T. C. Kofan´e, Eur. Phys. J. B 85 (2012) 318. [4] J. D. Tchinang Tchameu, C. Tchawoua and A. B. M. Togueu, Wave Motion 65 (2016) 112. [5] C. B. Tabi, A. Mohamadou and T. C. Kofan´e, Eur. Phys. J. D 86 (2008) 374. [6] A. K. Dang, C. B. Tabi, H. P. F. Ekobena and T. C. Kofan´e, Int. J. Quant. Chem. 115 (2015) 34. [7] C. B. Tabi, A. Mohamadou and T. C. Kofan´e , Phys. Lett. A 373 (2009) 2476. [8] C. B. Tabi, J. Phys.: Cond. Matter 22 (2010) 414107. [9] A. S. Davydov, Solitons in Molecular Systems, (Kluwer Academic Publishers, Dordrecht, 1981). [10] A. S. Davydov and A. D. Suprun , Ukr. J. Phys. 19 (1974) 44. [11] A. S. Davydov and N. I. Kislukha, Phys. Stat. Sol. (b) 39 (1973) 465. [12] M. Daniel and M. M. Latha, Physica A 298 (2001) 351. [13] E. A. Bartnik, J. A. Tuszynski and D. Sept, Phys. Lett. A 204 (1995) 263. [14] M. Daniel and M. M. Latha, Physica A 240 (1997) 526. [15] C. B. Tabi, J. C. Mimshe, H. P. Ekobena, A. Mohamadou and T. C. Kofan´e, Eur. Phys. J. B 86 (2013) 374. 24
[16] H. P. F. Ekobena, C. B. Tabi, A. Mohamadou and T. C. Kofan´e, J. Phys.: Condens. Matter 23 (2011) 375104. [17] C. B. Tabi, R.Y. Ondoua, H. P. Ekobena, A. Mohamadou and T. C. Kofan´e, Phys. Lett. A 380 (2016) 2374. [18] P. L. Christiansen, A. V. Savin and A. V. Zolotaryuk, J. Comput. Phys. 134 (1997) 1o8. [19] A. V. Savin, L. I. Manevich, P. L. Christiansen and A. V. Zolotaryuk, Phys. -Usp. 42 (1999) 245. [20] A. V. Savin, I. P. Kikot, M. A. Mazo and A. V. Onufriev, PNAS 110 (2013) 2816. [21] E. Fermi, J. Pasta and S. Ulam, in Collected Papers of Enrico Fermi, Vol. 2 (Univ. of Chicago Press, Chicago, 1965), p. 978 [22] C. B. Tabi, A. K. Dang, R. O. Doko, H. P. F. Ekobena and T. C. Kofan´e, Physica A 442 (2016) 498. [23] A. Mohamadou, A. J. Kenfack and T. C. Kofan´e, Phys. Rev. E 72 (2005) 036220 [24] K. Matsuoka and T. Kawahara, Wave Motion 38 (2003) 117. [25] I. Daumont, T. Dauxois and M. Peyrard, Nonlinearity 10 (1997) 617. [26] E.J. Parkes, Wave Motion 13 (1991) 261. [27] C. B. Tabi, A. S. Et´em´e, and A. Mohamadou, Physica A 474 (2017) 186. [28] I. Ma¨ına, C. B. Tabi, A. Mohamadou, H. P. F. Ekobena and T. C. Kofan´e, Chaos 25 (2015) 043118. [29] A. S. Et´em´e, C. B. Tabi and A. Mohamadou, Commun. Nonl. Sci. Num. Simul. 43 (2017) 211. 25
[30] C. B. Tabi, A. Mohamadou and T. C. Kofan´e, J. Phys.: Condes. Matter 20 (2008) 415104. [31] A. K. Dang, C. B. Tabi, H. P. F. Ekobena and T. C. Kofan´e, Chaos 115 (2012) 34. [32] C. B. Tabi, G. Bineli and A. Mohamadou, J. Biol. Phys. 41 (2015) 391. [33] J. C. F. Mimshe, C. B. Tabi, H. Edongue, H. P. Ekobena, A. Mohamadou and T. C. Kofan´e, Phys. Scr. 87 (2013) 025801. [34] S. Zdravkovi´c, A. N. Bugay, G. F. Aru and A. Maluckov, Chaos 24 (2014) 023139. [35] M. Peyrard and A.R. Bishop, Phys. Rev. Lett. 62 (1989) 2755. [36] T. Dauxois, M. Peyrard and A. R. Bishop, Phys. Rev. E 47 (1993) R44. [37] T. Dauxois and M. Peyrard, Phys. Rev. E. 51 (1995) 4027.
26
Highlights: • The nonlinear dynamics of complex zigzag molecules is addressed. • Transverse and longitudinal dynamics are shown to both contribute to energy localization via modulational instability • Thermal fluctuations are responsible for highly localized energy patterns.
S0165-2125(17)30052-5 http://dx.doi.org/10.1016/j.wavemoti.2017.04.008 WAMOT 2157
To appear in:
Wave Motion
Received date: 28 July 2016 Accepted date: 7 April 2017 Please cite this article as: B. Sadjo, C.B. Tabi, H. Edongue, A. Mohamadou, Coupled energy patterns in zigzag molecular chains, Wave Motion (2017), http://dx.doi.org/10.1016/j.wavemoti.2017.04.008 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Coupled Energy Patterns in Zigzag Molecular Chains Bertrand Sadjo1∗, Conrad B. Tabi1,2,3†, Hervais Edongue4‡, and Alidou Mohamadou5§ 1
Laboratory of Biophysics, Department of Physics, Faculty of Science, University of Yaound´e I, P.O. Box 812, Yaound´e, Cameroon.
2 3 4
Botswana International University of Science and Technology, Private Bag 16 Palapye, Botswana
The African Institute for Mathematical Sciences, 6-8 Melrose Rd, Muizenberg 7945, South Africa.
Laboratory of Material Sciences, Department of Physics, Faculty of Science, University of Yaound´e I, P.O. Box 812, Yaound´e, Cameroon.
5
Department of Physics, Faculty of Science, University of Maroua, P.O. Box 46, Maroua, Cameroon.
April 15, 2017
Abstract The contribution of both longitudinal and transversal nonlinear oscillations to energy localization is investigated in a zigzag molecular chain, which include simultaneously nearest- and next-nearest neighbor interactions. Coupled amplitude equations are found in the form of discrete nonlinear Schr¨odinger equations, whose plane wave solutions are found to be subjected to some instabilities. They are shown to be very sensitive to transverse and longitudinal couplings, which is confirmed via direct numerical simulations. The two available modes are found to be alternatively responsible for energy localization and transport. Thermal fluctuations effects bring about highly localized modes, along with narrow structures for efficient energy transport. [email protected] (B. Sadjo) Corresponding author (C. B. Tabi): [email protected]. Tel: +237-6788-735-50 ‡ [email protected] (H. Edongue) § [email protected] (A. Mohamadou)
∗
†
1
PACS number(s): 87.15.Aa, 05.45.Yv, 05.45.-a, 42.65.Tg Keywords: Zigzag chain, Modulational instability, Energy patterns.
Highlights: • The nonlinear dynamics of complex zigzag molecules is addressed. • Transverse and longitudinal dynamics are shown to both contribute to energy localization via modulational instability • Thermal fluctuations are responsible for highly localized energy patterns.
1
Introduction
The functional motion of proteins is strongly related to their structure. Usually, proteins are pictured as one-dimensional entities, generally involving only interactions among nearest-neighbors lattice sites [1, 3, 4]. During the last decades, considerable efforts have been made to show that beyond that simplified formulation of the problem, additional degrees of freedom can lead to a more suitable description of the way energy is transported and stored in biomolecules for specific purposes [1, 2]. In a few biological molecules like DNA, the obvious contribution of both longitudinal and transversal oscillations were already brought forth, showing how important their concomitant effects could be in the key processes of replication and transcription [5, 6, 7, 8]. On the other side, from the Davydov’s model [9, 10, 11], more complex models have been proposed, that show that for energy to be efficiently transported, the presence of at least three spines is important [12, 13, 14, 15, 16]. This has also been addressed recently [15, 17], where the strong impact
2
of radial and longitudinal displacements of peptides groups were confirmed to bring about new features in vibrational energy localization among hydrogen bonding spines. Albeit these considerable efforts to suitably capture and understand the energy transport phenomenon in biomolecules, it remains clear that much work still has to be done, which considers other aspects of their geometry, specifically, the mutual effect of their first and secondary structures. According to Ref. [18, 19], the simplest way to achieve this is to take into consideration first- and second-neighbor interactions, standing respectively for the deformation of valence bonds and the deformation of valence and torsion angles. In describing the interactions between first neighbors as transverse and those between second neighbors as longitudinal, one might be capable of describing the dynamics of more complex biological polymers by adopting their simplest two-dimensional structure in the form of zigzag chains [20], as represented in Fig. 1. The zigzag configuration is therefore a more generalized description of a broad range of molecular entities, whose the dynamics involves more than one degree of freedom, a good example being α−helix proteins, which include both rigid valence bonds (first-neighbor interactions) and soft hydrogen bonds (second-neighbor interactions). Compression soliton solutions, their stability and condition of existence were studied in Refs. [18, 19]. Indubitably, the existence of solitons implies the possibility of energy localization with reference to the pioneering work of Fermi, Pasta and Ulam [21]. One of the direct mechanisms leading to this phenomenon is modulational instability (MI), which to the best of our knowledge has not yet been used in the context of molecular zigzag chains. This constitutes the main objective of the present paper. MI is in fact the a process whereby a plane wave, under the competitive effects of dispersion and nonlinearity, breaks up into trains of soliton-like structures [22, 23, 24, 25, 26, 27, 28, 29]. In several occasions, it has also been shown to be responsible for the localization of energy in some biological lattice such as DNA [8, 30] and proteins [3, 15, 16, 17], with the obvious contribution of their structure and the involved interactions. In the rest of this work, we first present the zigzag model and we show that its transversal and longitudinal dynamics can be fully represented by
3
Figure 1: Schematic representation of a zigzag-shaped molecular chain, where U and V represent the intermolecular interaction potentials between first and second neighbors, respectively. xn and yn are respectively the longitudinal and transversal displacements, K1 and K2 are the stiffness constants for the first and second-neighbor forces. coupled discrete nonlinear Schr¨odinger (DNLS) equations, involving first and second dispersive and nonlinear terms. The linear stability analysis is thereafter performed and the existence of unstable patterns of energy is discussed. We also present some numerical applications of our predictions and discuss the features of the subsequent energy patterns, along with their behaviors under thermal fluctuations. Some concluding remarks are finally given.
2
Model and coupled amplitude equations
The planar zigzag model addressed in this work comprises two degrees of freedom xn (t) and yn (t) which represent the longitudinal and transversal displacements of the lattice site n with respect to their equilibrium positions. The total Hamiltonian that describes the dynamics of the whole lattice is written as H=
X n
(
1 M 2
"
dxn dt
2
+
dyn dt
2 #
)
+ Kl2 [U (rn ) + V (qn )] + ε0 Z(un , vn ) ,
4
(1)
where the first term stands for the kinetic energy, with M being the reduced mass of a chain particle. K is the stiffness constant, while intermolecular interactions between first and second neighbors are respectively described by the dimensionless potentials U (rn ) and V (qn ), and where the deviations from the equilibrium length rn and qn are such that s 2 2 1 xn+1 − xn yn + yn+1 + + h− −b rn = 2 l l s (2) 2 2 xn+1 − xn−1 yn+1 − yn−1 qn = 1+ + − 1. l l q b = h2 + 14 is the dimensionless equilibrium length of the rigid bonds and h is the thickness of the zigzag backbone, in units of the longitudinal lattice spacing l. The potentials U (r) and V (q) are normalized by U (0) = V (0) = 0 and U ”(0) = κ1 , V ”(0) = κ2 , where κ1 = K1 /K and κ2 = K2 /K are the dimensionless stiffness constants for the first- and second-neighbors forces, respectively. While ε0 describes the height of the potential Z(u, v), the latter brings about the zigzag features of the chain. The potential Z(u, v) is in fact made of two periodic functions f (u) and g(u), so that 0 ≤ f ≤ 1 and 0 ≤ g ≤ h. These two functions also have the periodicity properties f u + 21 = f (u)
and g(u + 1) = g(u), and their more explicit forms, as proposed in Ref. [18], are f (u) = sin2 (2πu) and g(u) = h sin2 (πu).
(3)
We assume f 00 ( n2 ) = Ω20 , with Ω0 being the dimensionless characteristic frequency of small oscillations of particles at the minima of the substrate potential Z(u, v), whose the complete expression is proposed in the form [18] 2 1 Z(u, v) = f (u) + Ω0 2 v − g(u) . 2
(4)
The features of the above potential are depicted in Fig. 2, versus the displacements u and v, for different values of the frequency Ω0 . In the vertices of the zigzag structure of the potential, the intermolecular bonds remain undistorted. The structure of this configuration does not change 5
Figure 2: Plots of the 2D potential versus u and v for h = 0.5 and: (a1)-(a2) Ω20 = 10, (b1)-(b2) Ω20 = 15 and (c1)-(c2) Ω20 = 20. with increasing Ω0 , but the limits between them and the maxima appear clearly. Furthermore, we q K introduce the dimensionless parameters t → M t, un (t) → xnl(t) , vn (t) → ynl(t) and we consider
the simplest harmonic form for the potentials U and V , such that U (r) = r2 /2 and V (q) = q 2 /2. Considering all this, Hamilton’s equations can be applied to the dimensionless version of (1), leading to the following coupled dynamical equations in un (t) and vn (t): 1 1 + un+1 − un + un − un−1 d2 un 2 2 =κ1 un+1 − 2un + un−1 − − dt2 1 + rbn 1 + rn−1 b 1 + un+2 − un 1 + un − un−2 ∂ − −ϑ Z(un , vn ) + κ2 un+2 − 2un + un−2 − 1 + qn+1 1 + qn−1 ∂un d2 v n h − vn−1 − vn h − vn − vn+1 =κ1 2h − vn−1 − 2vn − vn+1 − + dt2 1 + rn−1 1 + rn b b vn+2 − vn vn − vn−2 ∂ + κ2 vn+2 − 2vn + vn−2 − − −ϑ Z(un , vn ), 1 + qn+1 1 + qn−1 ∂vn
6
(5a)
(5b)
with ϑ =
ε0 . Kl2
After setting ϑ = 1 and expanding
∂Z ∂un
and
∂Z ∂vn
up to the third order, it is possible to
reduce Eqs.(5) to their amplitude equations. In performing such an expansion, one should keep in mind that the system should keep the balance between nonlinear and dispersive effects. Therefore, one of the techniques that usually work is the Fourrier series expansion method which assumes un and vn to have the trial forms [7, 8] un (t) =
+∞ X
ipωb t a(p) and vn (t) = n e
p=−∞
+∞ X
ipωb t b(p) , n e
(6)
p=−∞
where ωb is close to some linear oscillation frequencies. The Fourier coefficients are slowly varying (p)
with time, as they are functions of τ = 2 t, and they exponentially decay in p so that an ∼ p ∗ (p) (0) (0) (p) (−p) and bn ∼ p for p > 0, with an ∼ 2 (bn ∼ 2 , respectively). Moreover, an = an and ∗ (p) (−p) bn = bn . By applying trial solutions (6) to the expanded Eqs.(5), and considering the (1)
harmonics p = 1, one obtains a set of coupled equations whose time-dependent solutions an (t) and (1)
bn (t) are assumed to be in the form r 2ωb 8π 2 − ωb2 + α0 (3 + 4h2 ) + 2(κ2 − κ1 ) (1) an (t) = νn (t) exp i t p0 + 4α0 2ωb r 2ωb 2α0 (1 + 2h2 ) − 2κ1 − ωb2 + Ω0 2 (1) n bn (t) = (−1) ψn (t) exp i t , κ1 2ωb 4
where p0 = 2Ω0 2 π 4 h2 − 64π and α0 = 3
κ1 . 2b
(7a)
(7b)
Further mathematical expansions, using (7a), lead to the
following set of nonlinearly coupled equations iν˙n − P1 (νn+1 + νn−1 ) − P2 (νn+2 + νn−2 ) p0 |νn |2 νn + α0 (|νn+1 − νn |2 (νn+1 − νn ) − |νn − νn−1 |2 (νn − νn−1 )) − Q1 − κ (|ν 2 2 2 n+2 − νn | (νn+2 − νn ) − |νn − νn−2 | (νn − νn−2 )) + Q2 (|ψn − ψn+1 |2 (νn+1 − νn ) − |ψn − ψn−1 |2 (νn − νn−1 ))
− Q3 (|ψn+2 − ψn |2 (νn+2 − νn ) − |ψn − ψn−2 |2 (νn − νn−2 )) = 0
7
(8a)
iψ˙n − P3 (ψn+1 + ψn−1 ) + Q1
α0 (|ψn − ψn−1 |2 (ψn − ψn−1 ) + |ψn − ψn+1 |2 (ψn − ψn+1 ))
− κ (|ψ 2 2 2 n+2 − ψn | (ψn+2 − ψn ) − |ψn − ψn−2 | (ψn − ψn−2 ))
+ Q2 (|νn − νn−1 |2 (ψn − ψn−1 ) + |νn+1 − νn |2 (ψn − ψn+1 ))
(8b)
− Q4 (|νn+2 − νn |2 (ψn+2 − ψn ) − |νn − νn−2 |2 (ψn − ψn−2 )) = 0, with the different coefficients P1 = and P3 =
(1+h2 )α2 −κ1 , 2ωb
Q4 =
κ2 Q1 . 2
α0 (3+4h2 ) , 2ωb
P2 =
κ2 , 2ωb
Q1 = (p0 + 4α0 )−1 , Q2 = α0 /κ2 , Q3 =
κ2 , 4κ1
The two equations (8) are coupled DNLS equations which include
not only interactions between first neighbors, but also between second neighbors. Eq.(8a) describes the amplitude of longitudinal vibrations, while Eq.(8b) gives information regarding the transverse nonlinear oscillations. The two equations are nonlinearly coupled via the coefficients Q3 and Q4 which are functions of κ2 . Otherwise, they will completely be decoupled if κ2 → 0. Therefore, as initially described, each of the equations will consider only interactions between first neighbors. A single DNLS equation was already obtained in the context of charge transport in DNA [31] and coupled DNLS equations derived to describe nonlinear oscillations of DNA strands [7, 8, 32]. Their exact solutions were found in Ref. [7] to picture the way helicity and radial vibrations contribute to the unzipping and opening of the DNA duplex for the genetic code to be read. In the context of the present work, we intend to use the ones found here to perform the linear stability analysis of their plane wave solutions, in order to find regions where they are unstable, or more precisely, where solitonic and nonlinear energy patterns can be observed. However, direct numerical verifications will be performed on the generic Eqs. (5).
3
Linear stability analysis
The initial step to detect regions of parameters where modulated waves can be observed is to assume plane wave solutions for the system equations. For the specific case at hand, such solutions are supposed to be in the form νn = ν0 ei(k1 n−ω1 t) and ψn = ψ0 ei(k2 n−ω2 t) , where the frequencies ωj and
8
the wavenumbers qj , (j = 1, 2), satisfy the dispersion relations k1 4 4 2 ω1 = 2P1 cos(k1 ) + 2P2 cos(2k1 ) + Q1 |ν0 | p0 + 16κ2 sin (k1 ) − 16α0 sin 2 k1 k2 sin2 − Q3 sin2 (k1 ) sin2 (k2 ) + 16|ψ0 |2 Q2 sin2 2 2 k2 2 4 4 − k2 sin (k2 ) ω2 = 2P3 cos(k2 ) − 16Q1 |ψ0 | α0 sin 2 k2 k1 2 2 2 2 2 sin + Q4 sin (k1 ) sin (k2 ) . − 16|ν0 | Q2 sin 2 2
(9)
To further proceed, we assume ω1 = ω2 = ω and k1 = k2 = k, and one can get a relation between the two amplitudes ψ0 and ν0 as follows: (P3 − P1 ) cos(k) − P2 cos(2k) 8 (Q2 + α0 Q1 ) sin4 k2 (Q1 κ2 − Q3 ) sin4 (k) 16(α0 Q1 − Q2 ) sin4 k2 − Q1 p0 − 16(κ2 Q1 + Q4 ) sin4 (k) + |ν0 |2 . 16 (Q2 + α0 Q1 ) sin4 k2 + (Q1 κ2 − Q3 ) sin4 (k)
|ψ0 |2 =
(10)
To study the stability of the previous trial solutions, their behavior under small perturbations should be observed. These solutions can therefore be perturbed in amplitude or in phase, or both of them. However, when we consider the amplitude perturbations, i.e., νn = ν0 [1 + ηn (t)]ei(kn−ωt) and ψn = ψ0 [1 + ξn (t)]ei(kn−ωt) , we find the perturbations ηn and ξn to be governed by the coupled
9
set of linearized equations h i iη˙n + ωηn − P1 (ηn+1 + ηn−1 ) cos(k) + i(ηn+1 − ηn−1 ) sin(k) h i − P2 (ηn+2 + ηn−2 ) cos(2k) + i(ηn+2 − ηn−2 ) sin(2k) ∗ ∗ p0 (2ηn + ηn∗ ) + α0 [2(4 cos(k) − 3)ηn + (2 cos k − 2)(ηn+1 + ηn−1 + ηn∗ ) + (4 cos(k) − 2)(ηn+1 + ηn−1 ) + 4i(ηn+1 − ηn−1 ) sin(k)] − κ2 [(8 cos 2k − 6)ηn − Q1 |ν0 |2 ∗ ∗ + (2 cos(2k) − 2)(ηn+2 + ηn−2 + ηn∗ ) + (4 cos(2k) − 2)(ηn+2 + ηn−2 ) + 4i(ηn+2 − ηn−2 ) sin(2k)] ∗ ∗ 2 cos 2kξ + (4 cos(k) − 2)(η + ξ ) n n n ∗ ∗ 2 + 2(cos(k) − 1)(2ξn + ξn+1 + ξn−1 ) + (2 cos(k) − 1)(ξn+1 + ξn−1 ) + Q2 |ν0 | + 2(cos(k) − 1)(ηn+1 + ηn−1 ) + 2i(ηn+1 − ηn−1 ) sin k + 2i(ξn+1 − ξn−1 ) sin(k) ∗ ∗ 2 cos 4kξn∗ + (4 cos(2k) − 2)(ηn + ξn∗ ) + (2 cos(2k) − 2)(2ξn + ξn+2 + ξn−2 ) 2 + (2 cos(2k) − 1)(ξn+2 + ξn−2 ) + (2 cos(2k) − 2)(ηn+2 + ηn−2 ) =0 − Q3 |ψ0 | + 2i(ηn+2 − ηn−2 ) sin(2k) + 2i(ξn+2 − ξn−2 ) sin(2k) (11a)
10
h i ˙ iξn + ωξn − P3 (ξn+1 + ξn−1 ) cos(k) + i(ξn+1 − ξn−1 )sin(k)) ∗ ∗ α0 [(ξn+1 + ξn−1 )(2 − 2 cos k) − 4i(ξn+1 − ξn−1 ) sin(k) + (ξn+1 + ξn−1 )(2 − 2 cos k) ∗ ∗ + 8(1 − cos(k))ξ + (2 − 4 cos k)ξ + 2 cos 2kξ + 2(ξ + ξ ) cos 2k n n+1 n−1 n n 2 + 2i(ξn+1 − ξn−1 ) sin 2k] + κ2 [(ξn+2 + ξn−2 )(−2 + 2 cos(2k)) + 4i(ξn+2 − ξn−2 ) sin(2k) − Q1|ψ0 | ∗ ∗ ∗ + (ξn+2 + ξn−2 )(−2 + 2 cos 2k) + 8(−1 + cos(2k))ξn + (−2 + 4 cos(2k))ξn ∗ − 2 cos(4k)ξn − 2(ξn+2 + ξn−2 ) cos(4k) − 2i(ξn+2 − ξn−2 ) sin(4k)] (ξ + ξ )(1 − 2 cos(k)) + 2i(ξ − ξ ) sin(k) + (η + η )(1 − 2 cos(k)) n+1 n−1 n+1 n−1 n+1 n−1 ∗ ∗ − 2i(ηn+1 − ηn−1 ) sin(k) + (ηn+1 + ηn−1 )(2 − 2 cos(k)) + 4(ξn + ηn )(1 − cos k) 2 − Q2 |ν0 | + 2(1 − 2 cos k)ηn∗ − 2 cos(2k)ηn∗ + (ξn+1 + ξn−1 ) cos(2k) + i(ξn+1 − ξn−1 ) sin(2k) + (ηn+1 + ηn−1 ) cos(2k) + i(ηn+1 − ηn−1 ) sin(2k) (ξn+2 + ξn−2 )(−1 + 2 cos 2k) + 2i(ξn+2 − ξn−2 ) sin 2k + (ηn+2 + ηn−2 )(−1 + 2 cos 2k) ∗ ∗ + 2i(ηn+2 − ηn−2 ) sin 2k + (ηn+2 + ηn−2 )(2 cos 2k − 2) + 4(ξn + ηn )(−1 + cos 2k) =0 − Q4 |ν0 |2 ∗ ∗ + 2(−1 + 2 cos 2k)η + 2 cos 4kη − (ξ + ξ ) cos 4k − i(ξ − ξ ) sin 4k n+2 n−2 n+2 n−2 n n − (ηn+2 + ηn−2 ) cos 4k − i(ηn+2 − ηn−2 ) sin 4k (11b) ∗ t)
whose solutions can be suggested in the form ηn (t) = B1 ei(Kn−Ωt) + B2∗ e−i(Kn−Ω
and ξn (t) =
∗
B3 ei(Kn−Ωt) +B4∗ e−i(Kn−Ω t) , with K and Ω being the wavenumber and frequency of the perturbation, respectively. Reporting this into Eqs.(11) and making use of the dispersion relations (9), one obtains a set of four homogeneous equations for B1 , B2 , B3 and B4 , i.e., M × (B1 , B2 , B3 , B4 )T = 0, where the system matrix is given by a11 + Ω a12 a13 a14 a21 a22 − Ω a23 a24 M = a31 a32 a33 + Ω a34 a41 a42 a43 a44 − Ω 11
(12)
(13)
Figure 3: MI gains along with their corresponding stability/instability diagrams, versus the wavenumbers k and K for Ω20 = 10 and: (a1)-(a2) κ2 = 0.005κ1 , (b1)-(b2) κ2 = 0.01κ1 and (c1)-(c2) κ2 = 0.5κ1 , with κ1 = 0.85. with the matrix elements given in the Appendix. The above system admits non-trivial solutions given that the determinant of its matrix is null. This results into a fourth-order polynomial in Ω, whose solutions are either real or complex. However, we are interested in complex solutions Ω = Ωr + iΩi . Nevertheless, finding then is quite long and cumbersome, which requires some symbolic computations in order to extract such solutions whose the imaginary part gives the MI gain. In fact, the latter is the largest value among those corresponding to the various branches of the dispersion relation. In what follows, the gain is computed through the formula G(k, K) = |Im(Ω)|. The corresponding results are displayed in Figs. 3 and 4. The value of κ1 has been fixed and one plays both on the values of the frequency Ω0 and the force constant κ2 between second neighbors. 12
Figure 4: MI gains along with their corresponding stability/instability diagrams, versus the wavenumbers k and K for Ω20 = 15 and: (a1)-(a2) κ2 = 0.005κ1 , (b1)-(b2) κ2 = 0.01κ1 and (c1)-(c2) κ2 = 0.5κ1 , with κ1 = 0.85.
13
For example, in Fig. 3, we have fixed Ω20 = 10, and the gain shows only one region of instability for κ2 = 0.005κ1 . In this case, one clearly understands that only strong effects from the firstneighbor coupling contribute to the emergence of the instability zone. When κ2 increases to 0.01κ1 (Fig. 3(b1)-(b2)), the MI gain clearly shows the presence of another instability zone, while the primary one gets reduced. Following that, it therefore appears that increasing the value of κ2 is responsible for the explosion of the instability domain, as this is indeed confirmed in Figs. 3(c1)(c2), where Ω20 = 15. Obviously, from the instability diagram, the primary zone of instability has reduced, while there appear other two regions of instability. To remind, in regions of instability, the plane wave is supposed to be disintegrated into trains of solitary waveforms, or it is said to be stable when parameter are taken from the stability region. In Fig. 4, we have fixed Ω20 = 15, and we equally play on the value of κ2 , which appears to importantly influence the instability features. Contrarily o the previous case, for κ2 = 0.005κ1 (Fig. 4 (a1)-(a2)), there are three regions of instability. Again, this case shows the strong impact of the first-neighbor coupling, where the competition between the two directions of vibration is minimized. However, increasing κ2 to 0.01κ1 , the merged instability region is delocalized, therfore completely changing the region of parameters where modulated patterns are expected. Finally, when κ2 = 0.5κ1 , there is only one large region of instability as clearly shown by Fig. 4 (c1)-(c2). The behaviors induced by large and small values of Ω0 are a bit different, except that the two cases gives rise to instability regions. Also, there is an overlap of many instability points, which show that to some extent, modulated waves and nonlinear patterns of energy might be observed, with almost the same behaviors depending on the strength of the two coupling parameters. On the side, this cannot be confirmed by this analysis, which is just a procedure that gives information on the onset of MI.
14
Figure 5: Panel (a) shows the energy density as a function of time (increasing from bottom to top) and space (increasing from left to right) for κ2 = 0.005κ1 , k = 0.4π, K = 0.25π. (b1)-(b2) show the spatial expansion of longitudinal and transversal manifestations of MI at the respective instants t = 200 and t = 800. The time series for un and vn are also shown at the lattice points n = 50 and n = 150 (panels (c1) and (c2)).
4
Patterns of energy and coupled wave effects
In this section, we perform direct numerical simulations on the model equations (5) using the fourth-order Rung-Kutta computational scheme with periodic boundary conditions and time-step ∆t = 10−4 . We adopt the initial conditions un (t = 0) = u0 [1 + ζ cos(Kn)] cos(kn), u˙ n (t = 0) = u0 ω[1 + ζ cos(Kn)] sin(kn) (14) vn (t = 0) = v0 [1 + ζ cos(Kn)] cos(kn), v˙ n (t = 0) = v0 ω[1 + ζ cos(Kn)] sin(kn), with ζ being the perturbation strength fixed as ζ = 10−2 . Numerical simulations are performed on a lattice of 200 points and one extracts the information related to the contribution of both the longitudinal and transversal vibrations to the localization of energy. For this purpose, the energy is computed via its dimensionless expression " # X 1 dun 2 dvn 2 κ1 κ1 E= + + rn2 + qn2 + Z(un , vn ). 2 dt dt 2 2 n 15
(15)
For the wavenumbers k = 0.4π and K = 0.25π, with u0 = 0.5 and v0 = 0, an example of such simulations is given by Fig. 5, where we have fixed κ2 = 0.1κ1 . In Fig. 5 (a), we observe some coherent structures due to the activation of MI. The corresponding patterns are localized at specific sequences of lattice points and form rows over the time. This comportment is known, as it has been reported in numbers of models related to DNA [8, 32] and α−helix dynamics [16, 17]. However, what makes them special here is that they are the consequence of the competition between longitudinal and transversal oscillations. The value of κ2 , which is smaller than κ1 is expected not to really affect the pattern of energy, but at the same time, in Figs. 5 (b1) and (b2), one can notice that oscillations, from the two mode of vibrations, display train of waves, although transverse ones are more localized for t = 200. However, for t = 800, both of them almost display the same amplitude, which show that over long time, they all contribute to energy localization. The latter is also confirmed by Figs. 5 (c1) and (c2), where the time series of the solitonic structure are depicted at positions n = 50 and n = 150 and where the emerging phenomena is the out-of-phase oscillations of the two modes. This suggest that, each of the modes contributes to maintain energy in the system, making it available not only for specific metabolisms, but also to entertain the intrinsic dynamics of the molecule, as recently reported in the case of complex proteins with interspine coupling [15, 17, 33]. It was in fact concluded that, when two or more spines are coupled via hydrogen bonds, two of them might serve as energy reservoir, while the proper energy transport is ensured by the other. The same is also true for microtubules, as addressed in Ref. [34], where energy is always available to maintain intrinsically the molecule. Importantly, when κ2 increases to 0.5κ1 , the energy density increases and rows of energy get more and more restricted to small sequences of lattice points, as shown in Fig. 6(a). The main reason for this result is the features of transversal and longitudinal oscillations displayed in Fig. 6 (b1) and (b2), where trains of solitons are observed in both cases. Nevertheless, as only un has been excited, the behavior that appears in Fig. 6 (b1), at t = 200, seems to be the one expected. However, over long time evolution, i.e., t = 800 (Fig. 6 (b2)), transverse oscillations
16
Figure 6: Energy patterns (panel (a)) and solitonic structures due to MI (panels (b1)-(b2)) for κ2 = 0.5κ1 , k = 0.4π, K = 0.25π. In (a) time goes from bottom to top and space increases from left to right. Panels (c1) and (c2) show the time series of un and vn at the points n = 50 and n = 150. take advantage and contribute apparently to effectively localize energy. The time series of these waves still maintain the fact that the dynamics of the two modes remains out-of-phase, but one notices the increase in temporal oscillations at both n = 50 and n = 150. Also, the amplitude of such waves is not uniform as there are periods of high-amplitude oscillations. Nevertheless, the conclusion on the contribution of the two modes to energy localization remains the same as in Fig. 5. Energy patterns in proteins, and in biomolecules in general, are very sensitive to thermal fluctuations [16, 35, 36, 37]. This should however be consider, not only to evaluate their effects, but also because their inherent to biological media. For example, the fluctuational opening of DNA involves thermal fluctuations, and this enhances the creation of the bubbles that expose the bases out of the stack. As a whole, such mechanisms involves energy, which is normally brought by the hydrolysis of ATP, and conveyed to specific sites by specific enzymes such as RNA-polymerase. The latter surfs the long DNA molecule in search of the initiation sequence, and invests that energy for the replications process to be initiated [8, 32, 36]. It is evidenced, from a recent contribution, that 17
Figure 7: Noise effect on MI wave patterns for κ2 = 0.005κ1 , k = 0.4π, K = 0.25π. The energy density is represented with time increasing from bottom to top and space increasing from left to right (panel (a)). Panels (b1) and (b2) picture the spatial evolution of un and vn at instants t = 200 and t = 800, respectively. Panels (c1) and (c2) show the temporal evolution of un and vn at the lattice points n = 50 and n = 150. thermal fluctuations and dissipation can enhance high energy localization [8, 16, 32], resulting to the change in the vibration modes of the molecules. In the present case, to confirm that, we apply the random force Γn (t), with the statistical properties hΓn (t)i = 0 and hΓn (t)Γ0n0 (t)i = 2T γδn,n0 δ(t − t0 ),
(16)
to the two equations (5), where T is the temperature and γ the dimensionless friction coefficient taken to be equal to 10−2 . With in mind the study of the same effects as previously, we make numerical calculations at the temperature T = 300K, and we fix κ2 = 0.005κ1 . The corresponding results are depicted in Fig. 7, where for weak couplings, the coherency of the oscillations is destroyed, due partly to the increase in frequency of the oscillations. However, from Fig. 7(a), there are regions where energy is highly localized, due to the activation of MI. The time series shown in Figs. 7 (c1) and (c2) display complex oscillations at positions n = 50 and n = 150. Again, the amplitudes of 18
Figure 8: Noise effect on MI wave patterns for κ2 = 0.5κ1 , k = 0.4π, K = 0.25π. The energy density is plotted in panel (a) versus time (increasing from bottom to top) and space (increasing from left to right). Panels (b1) and (b2) show un and vn versus the lattice index n, at instants t = 200 and t = 800, respectively. The two panels (c1) and (c2) depict the time series of un and vn at positions n = 50 and n = 150. oscillations at those two positions are different, but with the same waveform. These features change for κ2 = 0.5κ1 , where the energy is more localized than in Fig. 6(a). There is only one pulse of energy, which has a solitonic shape (Fig. 8(a)), probably due to MI mechanism. This is the case of molecules whose energy remains concentrated in specific regions. That energy is due to the pulse solitons observed in Fig. 8(b1) and (b2) at different instants. However, the time series of the lattice point n = 50 shows that transverse and longitudinal oscillations take place, even where energy is not highly localized, which of course confirms our arguments on the permanent availability of energy in biological molecules. At position n = 150, one also confirms that the two modes display the same dynamics, except that their amplitudes are again not equal, while it seems noisy with slight modulations. In the case of α−helix proteins interactions with second neighbors are related to the deformation of the hydrogen bonds, and their highly localized oscillations, as shown in Figs. 7(c2) are in agreement with the observation made for a set of three coupled spines. Therefore, under 19
thermal fluctuations, they are favorable to the formation of breathers and to important energy outcomes.
5
Conclusion
Energy localization in zigzag molecular chains has been addressed in this paper. After reducing the whole dynamics of the complex structure to a set nonlinearly coupled DNLS equations, MI of plane wave solutions has been addressed and regions of instability have been found to be very sensitive to the coupling between longitudinal and transversal dynamics. This has been confirmed by direct numerical simulations on Eqs.(5), with specific interest in the way the subsequent nonlinear oscillations taking place in the two modes contribute to energy localization and transport in the chain. We have observed that in the absence of thermal fluctuations, such a dynamics is outof-phase, and we have evidenced that the two modes alternatively entertain permanent energy availability, not only for specific metabolisms to be achieved, but also for the intrinsic dynamics of the molecule. However, under thermal fluctuations, energy has been found to be highly localized due to the collective dynamics of the two modes. Importantly, under strong longitudinal coupling, energy has been found to display soliton-like behaviors following the individual dynamics of the two modes.
20
Appendix: Elements of the matrix (13) a11 = 2P1 cos(k) + 2(Q1 κ2 (|ψ0 |2 + 4|ν0 |2 ) + (4Q3 + Q2 )|ψ0 |2 + Q4 |ν0 |2 ) cos(2k) − 2(Q1 κ2 |ν0 |2 + Q3 |ψ0 |2 ) cos(4k) + 2|ψ0 |2 (Q2 − Q3 ) − Q1 |ν0 |2 (6κ2 + 6α0 − p0 ) − 4(Q3 − Q2 )|ψ0 |2 − 2P1 cos(K + k) − P2 cos(2(K + k)) − 2p0 Q1 |ν0 |2 h i 2 − α0 Q1 |ν0 | 2(4 cos(k) − 3 − 2 cos(K)) + 8 cos(K + k) h i − Q2 |ψ0 |2 2(2 cos(k) − 1 − 2 cos(K)) + 4 cos(K + k) h i 2 − Q1 κ2 |ν0 | 2(4 cos(2k) − 3 − 2 cos(2K)) + 8 cos(2(K + k)) h i − Q3 |ψ0 |2 (4 cos(2k) − 2 − 4 cos(2K)) + 4 cos(2(K + k))
a12 = −|ν0 |2 [p0 Q1 − α0 Q1 |ν0 |2 (2 cos(k) − 2)(2 cos(K) + 1) + κ2 Q1 (2 cos(2k) − 2)(2 cos(2K) + 1)] h i a13 = −2Q2 |ψ0 |2 (2 cos(k) − 2 − cos(K)) + 2 cos(K + k) h i 2 − Q3 |ψ0 | 2(2 cos(2k) − 2 − cos(2K)) + 2 cos(2(K + k)) h i a14 = −Q2 |ψ0 |2 2 cos(2k) + 4 cos(k) − 2 + (4 cos(k) − 4) cos(K) h i 2 − Q3 |ψ0 | 2 cos(4k) + 4 cos(2k) − 2 + (4 cos(2k) − 4) cos(2K) , a21 = a12 h i a24 = −2Q2 |ψ0 |2 (2 cos(k) − 2 − cos(K)) + 2 cos(K − k) h i 2 − Q4 |ψ0 | 2(2 cos(2k) − 2 − cos(2K)) + 2 cos(2(K − k))
21
a22 = 2P1 cos(k) + 2(Q1 κ2 (|ψ0 |2 + 4|ν0 |2 ) + (4Q3 + Q2 )|ψ0 |2 + Q4 |ν0 |2 ) cos(2k) − 2(Q1 κ2 |ν0 |2 + Q3 |ψ0 |2 ) cos(4k) + 2|ψ0 |2 (Q2 − Q3 ) − Q1 |ν0 |2 (6κ2 + 6α0 − p0 ) − 4(Q3 − Q2 )|ψ0 |2 − 2P1 cos(K − k) − P2 cos(2(K − k)) − 2p0 Q1 |ν0 |2 h i − α0 Q1 |ν0 |2 2(4 cos(k) − 3 − 2 cos(K)) + 8 cos(K − k) h i − Q2 |ψ0 |2 2(2 cos(k) − 1 − 2 cos(K)) + 4 cos(K − k) h i − Q1 κ2 |ν0 |2 2(4 cos(2k) − 3 − 2 cos(2K)) + 8 cos(2(K − k)) h i − Q3 |ψ0 |2 (4 cos(2k) − 2 − 4 cos(2K)) + 4 cos(2(K − k))
a31 = −Q2 |ν0 |2 [2 cos(K) − 4(cos(k) − 1) − 4 cos(K + k) + 2 cos(K + 2k)] − Q4 |ν0 |2 [4(cos(k) − 1) − 2 cos(K) + 4 cos(2K + k) − 2 cos(2(K + 2k))] a32 = −Q2 |ν0 |2 [(4 − 4 cos(k)) cos(K) + 2(1 − 2 cos(k)) − 2 cos(2k)] − Q4 |ν0 |2 [(−4 + 4 cos(2k)) cos(2K) + 2(−1 + 2 cos(2k)) + 2 cos(4k)] a34 = −α0 Q1 ψ02 [(4 − 4 cos(k) cos(K) + 2(1 − 2 cos(k)) + 2 cos(2k)] − κ2 Q1 |ψ0 |2 [(−4 + 4 cos(k)) cos(2K) + 2(−1 + 2 cos(2k)) − 2 cos(4k)] a23 = a14 ,
a43 = a34 ,
a41 = a32
a42 = −Q2 |ν0 |2 [4(cos(2k) − 1) − 2 cos(2K) − 4 cos(K + k) + 2 cos((K − 2k))] − Q4 |ν0 |2 [4(cos(2k) − 1) − 2 cos(2K) − 4 cos(K − k) − 2 cos(2(K − 2k))]
22
a33 = 2P1 cos(k) + 2(Q1 κ2 (|ψ0 |2 + 4|ν0 |2 ) + (4Q3 + Q2 )|ψ0 |2 + Q4 |ν0 |2 ) cos(2k) − 2(Q1 κ2 |ν0 |2 + Q3 |ψ0 |2 ) cos(4k) + 2|ψ0 |2 (Q2 − Q3 ) − Q1 |ν0 |2 (6κ2 + 6α0 − p0 ) − 4(Q3 − Q2 )|ψ0 |2 − 2P3 cos(K + k) − Q2 |ν0 |2 [2(2 cos(K) − 2 cos(k)) − 4 cos(K + k) + 2 cos(K + 2k)] − α0 Q1 |ψ0 |2 [4 cos(K) + 8(1 − cosk) − 8 cos(K + k) + 4 cos(K + 2k)] − Q4 |ν0 |2 [4(cos(2k) − 1) − 2 cos(2K) + 4 cos(2(K + k)) − 2 cos(2(K + 2k))] − κ2 Q1 |ψ0 |2 [8(cos(2k) − 1) − 4 cos(K) + 8 cos(2(K + k)) − 4 cos(2(K + 2k))] a44 = 2P1 cos(k) + 2(Q1 κ2 (|ψ0 |2 + 4|ν0 |2 ) + (4Q3 + Q2 )|ψ0 |2 + Q4 |ν0 |2 ) cos(2k) − 2(Q1 κ2 |ν0 |2 + Q3 |ψ0 |2 ) cos(4k) + 2|ψ0 |2 (Q2 − Q3 ) − Q1 |ν0 |2 (6κ2 + 6α0 − p0 ) − 4(Q3 − Q2 )|ψ0 |2 − 2P3 cos(K − k) − Q2 |ν0 |2 [2(2 cos(K) − 2 cos(k)) − 4 cos(K − k) + 2 cos(K − 2k)] − α0 Q1 |ψ0 |2 [4 cos(K) + 8(1 − cosk) − 8 cos(K − k) + 4 cos(K − 2k)] − Q4 |ν0 |2 [4(cos(2k) − 1) − 2 cos(2K) + 4 cos(2(K − k)) − 2 cos(2(K − 2k))] − κ2 Q1 |ψ0 |2 [8(cos(2k) − 1) − 4 cos(K) + 8 cos(2(K − k)) − 4 cos(2(K − 2k))]
23
References [1] A. C. Scott, Phys. Rep. 217 (1992) 1. [2] A. C. Scott, Phys. Rev. A 26 (1982) 578. [3] R.Y. Ondoua, C. B. Tabi, H. P. Ekobena, A. Mohamadou and T. C. Kofan´e, Eur. Phys. J. B 85 (2012) 318. [4] J. D. Tchinang Tchameu, C. Tchawoua and A. B. M. Togueu, Wave Motion 65 (2016) 112. [5] C. B. Tabi, A. Mohamadou and T. C. Kofan´e, Eur. Phys. J. D 86 (2008) 374. [6] A. K. Dang, C. B. Tabi, H. P. F. Ekobena and T. C. Kofan´e, Int. J. Quant. Chem. 115 (2015) 34. [7] C. B. Tabi, A. Mohamadou and T. C. Kofan´e , Phys. Lett. A 373 (2009) 2476. [8] C. B. Tabi, J. Phys.: Cond. Matter 22 (2010) 414107. [9] A. S. Davydov, Solitons in Molecular Systems, (Kluwer Academic Publishers, Dordrecht, 1981). [10] A. S. Davydov and A. D. Suprun , Ukr. J. Phys. 19 (1974) 44. [11] A. S. Davydov and N. I. Kislukha, Phys. Stat. Sol. (b) 39 (1973) 465. [12] M. Daniel and M. M. Latha, Physica A 298 (2001) 351. [13] E. A. Bartnik, J. A. Tuszynski and D. Sept, Phys. Lett. A 204 (1995) 263. [14] M. Daniel and M. M. Latha, Physica A 240 (1997) 526. [15] C. B. Tabi, J. C. Mimshe, H. P. Ekobena, A. Mohamadou and T. C. Kofan´e, Eur. Phys. J. B 86 (2013) 374. 24
[16] H. P. F. Ekobena, C. B. Tabi, A. Mohamadou and T. C. Kofan´e, J. Phys.: Condens. Matter 23 (2011) 375104. [17] C. B. Tabi, R.Y. Ondoua, H. P. Ekobena, A. Mohamadou and T. C. Kofan´e, Phys. Lett. A 380 (2016) 2374. [18] P. L. Christiansen, A. V. Savin and A. V. Zolotaryuk, J. Comput. Phys. 134 (1997) 1o8. [19] A. V. Savin, L. I. Manevich, P. L. Christiansen and A. V. Zolotaryuk, Phys. -Usp. 42 (1999) 245. [20] A. V. Savin, I. P. Kikot, M. A. Mazo and A. V. Onufriev, PNAS 110 (2013) 2816. [21] E. Fermi, J. Pasta and S. Ulam, in Collected Papers of Enrico Fermi, Vol. 2 (Univ. of Chicago Press, Chicago, 1965), p. 978 [22] C. B. Tabi, A. K. Dang, R. O. Doko, H. P. F. Ekobena and T. C. Kofan´e, Physica A 442 (2016) 498. [23] A. Mohamadou, A. J. Kenfack and T. C. Kofan´e, Phys. Rev. E 72 (2005) 036220 [24] K. Matsuoka and T. Kawahara, Wave Motion 38 (2003) 117. [25] I. Daumont, T. Dauxois and M. Peyrard, Nonlinearity 10 (1997) 617. [26] E.J. Parkes, Wave Motion 13 (1991) 261. [27] C. B. Tabi, A. S. Et´em´e, and A. Mohamadou, Physica A 474 (2017) 186. [28] I. Ma¨ına, C. B. Tabi, A. Mohamadou, H. P. F. Ekobena and T. C. Kofan´e, Chaos 25 (2015) 043118. [29] A. S. Et´em´e, C. B. Tabi and A. Mohamadou, Commun. Nonl. Sci. Num. Simul. 43 (2017) 211. 25
[30] C. B. Tabi, A. Mohamadou and T. C. Kofan´e, J. Phys.: Condes. Matter 20 (2008) 415104. [31] A. K. Dang, C. B. Tabi, H. P. F. Ekobena and T. C. Kofan´e, Chaos 115 (2012) 34. [32] C. B. Tabi, G. Bineli and A. Mohamadou, J. Biol. Phys. 41 (2015) 391. [33] J. C. F. Mimshe, C. B. Tabi, H. Edongue, H. P. Ekobena, A. Mohamadou and T. C. Kofan´e, Phys. Scr. 87 (2013) 025801. [34] S. Zdravkovi´c, A. N. Bugay, G. F. Aru and A. Maluckov, Chaos 24 (2014) 023139. [35] M. Peyrard and A.R. Bishop, Phys. Rev. Lett. 62 (1989) 2755. [36] T. Dauxois, M. Peyrard and A. R. Bishop, Phys. Rev. E 47 (1993) R44. [37] T. Dauxois and M. Peyrard, Phys. Rev. E. 51 (1995) 4027.
26
Highlights: • The nonlinear dynamics of complex zigzag molecules is addressed. • Transverse and longitudinal dynamics are shown to both contribute to energy localization via modulational instability • Thermal fluctuations are responsible for highly localized energy patterns.