1 Introduction

The development of the theory of the modulation instability (MI) as a well-known nonlinear phenomenon has attracted considerable attention and started almost simultaneously and occurred in parallel in hydrodynamics by Benjamin and Feir [1, 2] and in electrodynamics by Ostrovskii et al. [3, 4]. It can be found in many fields of physics, such as fluid dynamics [5], nonlinear optics. MI, an effect well known in the latter context, is the process by which a constant-wave background becomes unstable to sinusoidal modulations because of the presence of a focusing nonlinearity, leading to a pulsed wave, called a bright soliton train [6, 7]. MI in optical fibers was originally studied by Hasegawa and Brinkman [6] and their theoretical prediction was observed experimentally by Tai et al. [7].

Nonlinear excitations [8,9,10,11,12,13,14,15,16,17,18,19] (solitons, discrete breathers, intrinsic localized modes) have been drawing increasing attention over recent years and are widely believed to be responsible for several effects in molecular chains, such as charge and thermal conductivity, energy transfer and localization. A particularly interesting discrete system that support solitons and localized modes is deoxyribonucleic acid, or DNA. In this system, localization of energy has been suggested as a precursor of the transcription bubble [10], and moving localized oscillations as a method of transport of information along the double strand [15]. For this reason Fialko and Lakhno have studied the transport of charge and hole along the short DNA molecule by using the Peyrard-Bishop-Holstein (PBH) model [20, 21]. These authors investigate the impact of long-range transfer of charges through the DNA molecule. Many authors have shown that the most standard mechanism through which bright solitons or solitary wave structures appear is through the activation of MI of plane waves. In the above-mentioned contexts, MI has been suggested to be responsible for energy localization mechanisms leading to the formation of large amplitude nonlinear excitations in hydrogen-bonded crystals or DNA molecules.

Peyrard and Bishop [10] studied the process of denaturation in which only the transverse motion of bases along the hydrogen bond was taken into account. One of its improved version, proposed by Dauxois et al. [8], takes helicity into consideration and has been recently shown, by Zdravkovic et al. [11, 12] and then by Tabi et al. [13, 14], to better describe the extremely high amplitude stretching of the hydrogen bonds. It is therefore obvious that all the above-cited models describe either the radial opening of strands or the unwinding of the double helix, but not both of them simultaneously. In more recent years, due to experiments on single molecules of DNA, models with two and more degrees of freedom have been introduced with emphasis on the radial and torsional aspects. As an example, a model which has two degrees of freedom per base pair: one radial variable related to the opening of the hydrogen bonds and an angular one related to the twisting of each base-pair responsible for the helicoidal structure of the molecule, has been built by Barbi et al. [15], and later improved by Cocco and Monasson [16]. Such a model provides an extension of the Peyrard-Bishop approach towards a more realistic description of biological processes. It can also be useful for more general studies of the interaction between geometrical conformation and dynamical properties of the molecule referring to the recent mechanical experiments on DNA [17]. Nonlinear interactions between atoms in DNA can give rise to intrinsically localized breather-like vibration modes [18]. Such localized modes, being large amplitude vibrations of a few (2 or 3) particles, can facilitate the disruption of base pairs and in this way initiate conformational transitions in DNA. These modes can occur as a result of the modulational instability of continuum-like nonlinear modes [19], which is created by energy exchange mechanisms between the nonlinear excitations. The latter favors the growth of large excitations [22]. But, as far as we know, no work concerning analysis of MI through the Peyrard-Bishop model of DNA with rotational motion for the nucleotides has been presented in the literature. So, the aim of the present work is to show that the nonlinear dynamics of DNA can be described by coupled discrete equations. The paper is organized as follows. In Section 2, we propose the model and we derive the equations of motion for exciton and phonon excitations. In Section 3, linear analysis is investigated and predictions of some localized structure formations are made. The validity of this analysis is checked by numerical simulations in Section 4. Section 5 concludes the paper.

2 Model and equations of motion

The Peyrard-Bishop (PB) model was introduced in 1989 to explain the fluctuations and the nonlinear phenomena observed in the DNA molecule [10]. Initially for this model, only one degree of freedom per base pair was recognized, the opening along the hydrogen bond of base pair n. Thereafter, the model with several degrees of freedom per base pair was elaborated by taking into account simultaneously the transverse and longitudinal motions characterized by u j and v j starting from the equilibrium position and the rotational motions characterized by 𝜃 j and ϕ j [23]. Each nucleotide is represented by a disc of mass m, the radius r and d represent in the equilibrium position between the discs. On each base pair, forces of interactions with the other base pair are exerted, which are represented by a nonlinear potential V (u j ,v j ,𝜃 j ,ϕ j ) [23]. The displacement of the H bond between two adjacent discs j used in this present work, can be written as [23].

$$ y_{j}=(u_{j}-v_{j}+d+2r)-r(\cos\theta_{j}+\cos\phi_{j}). $$
(1)

The lattice Hamiltonian for the system can be written as [23]

$$\begin{array}{@{}rcl@{}} H_{1}&=&\sum\limits_{j}\left\{\frac{p_{w_{j}}^{2}}{2m}+\frac{p_{\lambda_{j}}^{2}}{2m} +\frac{p_{\phi_{j}}^{2}}{2I}+\frac{p_{\theta_{j}}^{2}}{2I}\right\}+ \sum\limits_{j}\left\{\frac{k}{2}\left[(w_{j}-w_{j-1})^{2}+(\lambda_{j}-\lambda_{j-1})^{2}\right]\right\}\\ &&+ \sum\limits_{j}\left\{\frac{\xi}{2}\left[(\phi_{j}-\phi_{j-1})^{2}+(\theta_{j}-\theta_{j-1})^{2}\right]\right\}+\sum\limits_{j} U(x_{j}) \end{array} $$
(2)

where

$$\begin{array}{@{}rcl@{}} w_{j}&=&\frac{u_{j}+v_{j}}{\sqrt{2}},\\ \lambda_{j}&=&\frac{u_{j}-v_{j}}{\sqrt{2}}, \end{array} $$
(3)

In this equation, \( p_{w_{j}}\) (or \(p_{\lambda _{j}}\)) and \( p_{\phi _{j}}\) (or \(p_{\theta _{j}}\)) are the linear and angular moments, respectively. We consider that the mass m, the moment of inertia I are equal for all nucleotides [23]. Also, the strength constants k and ξ are the same for the system. In this present work we consider only the artificial and the homogeneous DNA molecule. Consequently all base pairs are identical. All parameters used in the work have been taken from in the references [24,25,26]. The interaction between two adjacent bases is given by the Morse potential [25]

$$ U\left( {y_{j}} \right) = D\left[ {\exp \left( {- a(y_{j}-y_{0}) } \right) - 1} \right]^{2}= D\left[ {\exp \left( {- a\sqrt{2}\lambda_{j}-g(\phi_{j},\theta_{j}) } \right) - 1} \right]^{2} $$
(4)

where D is the energy of dissociation of the base pair, is a parameter with dimension of inverse length, y 0 = 2r + d is the equilibrium point (distance between the centers of the discs), and the function g(ϕ j ,𝜃 j ) is given by

$$ g(\phi_{j},\theta_{j})=r(\cos\phi_{j}+\cos\theta_{j}). $$
(5)

In the following, the above extended PB model will be coupled to the excitons. Thus, the Hamiltonian of excitons for the system can be written as

$$ H_{2} = - {\sum}_{j} {V\left( {B_{j}^ + B_{j + 1} + B_{j}^ + B_{j - 1} } \right)}. $$
(6)

This Hamiltonian describes the transfer of charges between the pairs of bases in the adjacent stacking of the molecular orbitals which superimpose along the DNA molecule. \(B_{j}^{+}\) and B j are creation and annihilation operators, respectively, for charge carrier at the j th base pair of the double strand, and V represents the coupling transfer integral between π orbital at adjacent base pairs. For the charge-lattice interaction, we use on-site, Holstein-type coupling [26, 27]. Thus, the Hamiltonian interaction of the system can be expressed as:

$$ H_{3} = {\sum}_{j} {\chi \lambda_{j} } B_{j}^ + B_{j}, $$
(7)

where χ is the coupling vibrational and rotational constant. The total Hamiltonian of our system can be written as H = H 1 + H 2 + H 3. Now, we use semiclassical equations of motion [26], i.e, we treat the quantum charge mechanically and the vibrational and rotational motion classically. We obtain the following Schrödinger equation:

$$ i\hbar \frac{d}{dt}\varphi_{j} = -V\left( {\varphi_{j + 1} + \varphi_{j - 1} } \right) + \chi \lambda_{j} \varphi_{j} $$
(8)

where φ j is the probability amplitude for the charge carrier located at the j th base pair. Newton’s equations of motion for the strectching w j , λ j , and the rotational 𝜃 j , ϕ j motions become:

$$ m{ \ddot{w}_{j}}=k(w_{j+1}+w_{j-1}-2w_{j}), $$
(9)
$$\begin{array}{@{}rcl@{}} m\ddot{\lambda}_{j}&=&k(\lambda_{j+1}+\lambda_{j-1}-2\lambda_{j})+2\sqrt{2}aD(\exp(-a\sqrt{2}\lambda_{j}-g(\theta_{j}, \phi_{j}))-1)\\ &&\times (\exp(-a\sqrt{2}\lambda_{j}-g(\theta_{j},\phi_{j})))-\chi\lambda_{j}|\varphi_{j}|^{2}, \end{array} $$
(10)
$$\begin{array}{@{}rcl@{}} I\ddot{\phi}_{j}&=&\xi(\phi_{j+1}+\phi_{j-1}-2\phi_{j})-2raD\sin\phi_{j}(\exp(-a\sqrt{2}\lambda_{j}-g(\theta_{j}, \phi_{j}))-1)\\ &&\times (\exp(-a\sqrt{2}\lambda_{j}-g(\theta_{j},\phi_{j}))), \end{array} $$
(11)
$$\begin{array}{@{}rcl@{}} I\ddot{\theta}_{j}&=&\xi(\theta_{j+1}+\theta_{j-1}-2\theta_{j})-2raD\sin\theta_{j}(\exp(-a\sqrt{2}\lambda_{j}- g(\theta_{j},\phi_{j}))-1)\\ &&\times (\exp(-a\sqrt{2}\lambda_{j}-g(\theta_{j},\phi_{j}))). \end{array} $$
(12)

Equations (8)–(12) are equations of motion of our system. Let:

$$\begin{array}{@{}rcl@{}} X_{j}=\frac{\theta_{j}+\phi_{j}}{\sqrt{2}}, \\ {\Psi}_{j}=\frac{\theta_{j}-\phi_{j}}{\sqrt{2}}, \end{array} $$
(13)
$$\begin{array}{@{}rcl@{}} x_{j}=\frac{X_{j}}{\sqrt{2}},\\ \psi_{j}=\frac{{\Psi}_{j}}{\sqrt{2}}. \end{array} $$
(14)

Our priority is to transform (11)–(12) into one equation. Substituting (13)–(14) into (11)–(12), we obtain the following equation:

$$\begin{array}{@{}rcl@{}} I\ddot{\psi}_{j}&=&\xi(\psi_{j+1}+\psi_{j-1}-2\psi_{j})-2raD\sin\psi_{j}\cos x_{j}(\exp(-a\sqrt{2}\lambda_{j}-g(x_{j},\psi_{j}))-1)\\ &&\times(\exp(-a\sqrt{2}\lambda_{j}-g(x_{j},\psi_{j}))) \end{array} $$
(15)

where g(X j ,Ψ j ) = g(x j ,ψ j ) = 2r(cosx j cosψ j ). Equation (9) represents a pure harmonic lattice with plane solutions. Our attention will be focused on the nonlinear (10)–(12). For small angular motion we have cosx j ≃ 1, sinψ j ≃ ψ j and \(\exp (-a\sqrt {2}\lambda _{j}-g(x_{j},\psi _{j}))\simeq \exp (-a\sqrt {2}\lambda _{j})\) because g(x j ,ψ j ) ≃ 2r and exp(−2r) → 0 when r →∞ Using this approximation, and (15) becomes:

$$ I\ddot{\psi}_{j}=\xi(\psi_{j+1}+\psi_{j-1}-2\psi_{j})-2raD\psi_{j}(\exp(-a\sqrt{2}\lambda_{j})-1) (\exp(-a\sqrt{2}\lambda_{j})). $$
(16)

Finally, after using all cited approximations, the equations of motion becomes three discrete equations for three variables,

$$ i\hbar \frac{d}{dt}\varphi_{j} = -V\left( {\varphi_{j + 1} + \varphi_{j - 1} } \right) + \chi \lambda_{j} \varphi_{j}, $$
(17)
$$ I\ddot{\psi}_{j}=\xi(\psi_{j+1}+\psi_{j-1}-2\psi_{j})-2raD\psi_{j}(\exp(-a\sqrt{2}\lambda_{j})-1) (\exp(-a\sqrt{2}\lambda_{j})), $$
(18)
$$\begin{array}{@{}rcl@{}} m\ddot{\lambda}_{j}&=&k(\lambda_{j+1}+\lambda_{j-1}-2\lambda_{j})+2\sqrt{2}aD(\exp(-a\sqrt{2}\lambda_{j}-g(\theta_{j},\phi_{j}))-1)\\ &&\times (\exp(-a\sqrt{2}\lambda_{j}-g(\theta_{j},\phi_{j})))-\chi\lambda_{j}|\varphi_{j}|^{2}. \end{array} $$
(19)

After expanding the terms in the exponential up to the third order and using the continuum approximation, we obtain the following equations

$$ \ddot{\lambda}+c_{0}\frac{\partial^{2}\lambda} {\partial x^{2}}+W_{g}(\lambda-\alpha_{0}\lambda^{2}+\gamma_{2}\lambda^{3})+c_{2}|\varphi|^{2}=0, $$
(20)
$$ \ddot{\psi}+c_{1}\frac{\partial^{2}\psi} {\partial x^{2}}-\beta(\psi\lambda-\alpha_{0}\psi\lambda^{2}+\gamma_{2}\psi\lambda^{3})=0, $$
(21)
$$ i\frac{d}{dt}\varphi= P_{1}\frac{\partial^{2}\varphi}{\partial x^{2}}+Q_{1}\varphi+Q_{2}\lambda\varphi, $$
(22)

where \(W_{g}=\frac {4a^{2}D}{m}\), \(\alpha _{0}=\frac {3a\sqrt {2}}{2}\), \(c_{0}=\frac {-k}{m}\), \(c_{1}=\frac {-\xi }{I}\), \(\gamma _{2}=\frac {7a^{2}}{3}\), \(\beta =\frac { m r W_{g}\sqrt {2}}{2I}\), \(Q_{1}=-\frac {V}{\hbar }\), \(Q_{2}=\frac {\chi }{\hbar }\), \(P_{1}=-\frac {2V}{\hbar }\) and \(c_{2}=\frac {\chi }{m}\).

3 Linear stability analysis

In order to perform the linear stability analysis of system (20) and (21), we assume that:

$$\begin{array}{@{}rcl@{}} \varphi=\varphi_{0}\exp i(k x-w_{0}t),\\ \lambda=\lambda_{0} ,\\ \psi=\psi_{0}\exp i(k x-w_{0}t), \end{array} $$
(23)

with real constants w 0, λ 0, ψ 0 and φ 0 complex constant and where and are the wave number and frequency, respectively, of the system without perturbation. Introducing the above relation into (20) and (21), we get:

$$ w_{0}= \eta (\frac{A-k^{2}}{B+k^{2}}), $$
(24)

where

$$\begin{array}{@{}rcl@{}} A=\frac{\beta c_{2}}{W_{g} c_{1}}, \\ B=\frac{Q_{1}+Q_{2}\lambda_{0}}{P_{1}}, \\ \eta =\frac{c_{1}}{P_{1}} \end{array} $$
(25)

and

$$ k^{4}-\mu k^{2}+\mu_{1}=0 $$
(26)

where

$$ \mu=-\frac{1}{{P_{1}^{2}}}(2P_{1}(Q_{1}+Q_{2}\lambda_{0})+c_{1}) $$
(27)
$$ \mu_{1}=- \beta c_{2} \frac{|\varphi_{0}|^{2}}{W_{g}}+(Q_{1}+Q_{2}\lambda_{0})^{2}. $$
(28)

Equation (26) have solutions given by:

$$ k_{+}^{2}=\frac{\mu+(\sqrt{\mu^{2}-4\mu_{1})}}{2}, $$
(29)
$$ k_{-}^{2}=\frac{\mu-(\sqrt{\mu^{2}-4\mu_{1})}}{2}. $$
(30)

The major remark is that all the particles of the molecular chain oscillate in all the directions (see Fig. 1), creating nonlinear effects in the molecule. Figure 1 plots the coefficients of nonlinear dispersion relation versus the wave number K 1 and we note different behaviors in the dynamics of the base pairs. We see that the dispersion curve oscillates. This behavior is probably induced by the competition between the longitudinal and the rotational motions. Adding a small perturbation in above the equilibrium state, that is

$$\begin{array}{@{}rcl@{}} \varphi=(\varphi_{0}+\epsilon\varphi_{1})\exp i(K x+w_{0}t),\\ \lambda=\lambda_{0}+\epsilon\lambda_{1},\\ \psi=(\psi_{0}+\epsilon\psi_{1})\exp i(K x+w_{0}t), \end{array} $$
(31)

and using it in (20) and (21), with the help of condition (24), we write φ 0 = a 1 + i a 2, φ 1 = u + i v, ψ 0 = b 1 + i b 2, ψ 1 = u 1 + i v 2 where we are taking a 1 = a 2 for the sake of simplicity, and we separate imaginary and real parts as follows:

$$ P_{1}\frac{\partial^{2}u}{\partial x^{2}}-\frac{\partial v}{\partial t}+Q_{1}u+Q_{2}a_{1}\lambda_{1}=0, $$
(32)
$$ P_{1}\frac{\partial^{2}v}{\partial x^{2}}+\frac{\partial u}{\partial t}+Q_{1}v+Q_{2}a_{2}\lambda_{1}=0, $$
(33)
$$ c_{0}\frac{\partial^{2}\lambda_{1}}{\partial x^{2}}+\frac{\partial^{2} \lambda_{1}}{\partial t^{2}}+w_{g}(1-2\alpha_{0}\lambda_{0}+3\gamma_{2}{\lambda_{0}^{2}})\lambda_{1}+2c_{2}(ua_{1}+a_{2}v)=0, $$
(34)
$$\begin{array}{@{}rcl@{}} \frac{\partial^{2}u_{1}}{\partial t^{2}}&+&2w_{0}\frac{\partial v_{2}}{\partial t}-{w_{0}^{2}}u_{1}+c_{1}\left( \frac{\partial^{2} u_{1}}{\partial x^{2}}-2K\frac{\partial v_{2}}{\partial x}-K^{2}u_{1}\right)\\ &-&\beta((b_{1}\lambda_{1}+\lambda_{0} u_{1})-2\alpha_{0}\lambda_{0}b_{1}\lambda_{1}- 3b_{1}{\lambda_{0}^{2}}\gamma_{2}\lambda_{1}-\gamma_{2}{\lambda_{0}^{3}}u_{1})=0, \end{array} $$
(35)
$$\begin{array}{@{}rcl@{}} \frac{\partial^{2}v_{2}}{\partial t^{2}}&-&2w_{0}\frac{\partial u_{1}}{\partial t}-{w_{0}^{2}}v_{2}+c_{1}\left( \frac{\partial^{2} v_{2}}{\partial x^{2}}+2K\frac{\partial u_{1}}{\partial x}-K^{2}v_{2}\right)\\ &-&\beta((b_{2}\lambda_{1}+\lambda_{0} v_{2})-2\alpha_{0}\lambda_{0}b_{2}\lambda_{1}- 3b_{2}{\lambda_{0}^{2}} \gamma_{2}\lambda_{1}-\gamma_{2}{\lambda_{0}^{3}}v_{2})=0. \end{array} $$
(36)

Then, we insert

$$ u=u_{01}\exp i(K_{1}x-{\Omega} t)+c.c., $$
(37)
$$ v=v_{01}\exp i(K_{1}x-{\Omega} t)+c.c., $$
(38)
$$ \lambda_{1}=\lambda_{01}\exp i(K_{1} x-{\Omega} t)+c.c., $$
(39)
$$ u_{1}=u_{10}\exp i(K_{1} x-{\Omega} t)+c.c., $$
(40)
$$ v_{2}=v_{20}\exp i(K_{1} x-{\Omega} t)+c.c., $$
(41)

into (32)–(36), where K 1 and Ω are the perturbation wave number and frequency, respectively. cc is the complex conjugate. We arrive to a five coupled linear equations

$$ M(u_{01},v_{01},\lambda_{01},u_{10}, v_{20})^{T}=0 $$
(42)

with

$$M=\left( \begin{array}{lllll} m_{11} & i{\Omega} & m_{13} & 0 & 0 \\ -i{\Omega} & m_{11} & m_{23} & 0 & 0 \\ m_{31} & m_{32} & -{\Omega}^{2}+m_{33} & 0 & 0 \\ 0 & 0 & m_{43} & -{\Omega}^{2}+m_{44} & -i(2w_{0}{\Omega}+m_{45})\\ 0 & 0 & m_{43} & i(2w_{0}{\Omega}+m_{45}) & -{\Omega}^{2}+m_{44} \end{array}\right) $$

where

$$\begin{array}{@{}rcl@{}} m_{13}=Q_{2}a_{1}, \\ m_{23}=Q_{2}a_{2}, \\ m_{32}=2c_{2}a_{2}, \\ m_{31}=2c_{2}a_{1}, \\ m_{11}=-P_{1}{K_{1}^{2}}+Q_{1}, \\ m_{45}=-c_{1}KK_{1}, \\ m_{33}=-c_{0}{K_{1}^{2}}+W_{g}(1+2\alpha_{0}\lambda_{0}+3{\lambda_{0}^{2}}\gamma_{2}), \\ m_{43}=-\beta b_{1}+2\alpha_{0}\lambda_{0}b_{2}+3{\lambda_{0}^{2}}\gamma_{2}b_{2}, \\ m_{44}=-{w_{0}^{2}}-c_{1}{K_{1}^{2}}-c_{1}K^{2}-\beta\lambda_{0}+{\beta\lambda_{0}^{2}}\gamma_{2}. \end{array} $$
(43)
Fig. 1
figure 1

Parameters of the nonlinear dispersion relation (44) for m = 300a m u, a = 4.41Å −1,χ = 1.2e VÅ −1,D = 0.05e V,V = 0.005e V,k = 0.04,r = 0.3Å −1,I = 300a m u,φ 0 = 10−3Å −1

Full size image

The condition that (42) has a nontrivial solution requires its determinant to be zero. This gives the eigenvalue equation:

$$ {\Omega}^{8}-P{\Omega}^{6}-T{\Omega}^{5}+M{\Omega}^{4}+N{\Omega}^{3}+R{\Omega}^{2}+Z{\Omega}+Q=0 $$
(44)

where

$$\begin{array}{@{}rcl@{}} P=m_{11}^{2}+m_{33}+2m_{44}-4{w_{0}^{2}}, \\ T=4w_{0}m_{45}, \end{array} $$
(45)
$$ N=w_{0}(4m_{33}m_{45}-4m_{11}^{2}m_{45}), $$
(46)
$$ M=2m_{33}m_{44}-2m_{31}m_{11}m_{13}+m_{11}^{2}(2m_{44}+m_{33}-4{w_{0}^{2}})+m_{44}^{2}+4m_{33}{w_{0}^{2}}- m_{45}^{2}, $$
(47)
$$\begin{array}{@{}rcl@{}} R&=&2m_{11}^{2}m_{33}m_{44}+4m_{23}m_{44}m_{32}m_{11}-m_{33}m_{44}^{2}-(m_{11}m_{44})^{2}+8m_{11}m_{31} m_{13}{w_{0}^{2}}\\ &&-m_{33}m_{45}^{2} -(m_{11}m_{45})^{2}-4{w_{0}^{2}}m_{11}^{2}m_{33}, \end{array} $$
(48)
$$ Z=w_{0}m_{45}(8m_{11}m_{13}m_{31}-4m_{11}^{2}m_{33}), $$
(49)
$$ Q=2m_{45}^{2}m_{11}m_{13}m_{31}-m_{11}^{2}m_{33}m_{45}^{2}-2m_{11}m_{13}m_{31}m_{44}^{2}+m_{11}^{2}m_{33} m_{45}^{2}. $$
(50)

In order to solve for Ω, we need to find the roots of (44), which is a eighth order polynomial equation. It is not easy to obtain analytical results. The coefficients P, T, M, N, R, Z and Q add to this difficulty, because a full listing requires several pages. In this framework, we plotted in Fig. 1, these seven quantities with respect to the wave number K 1, and the following results have been observed: To make sure that our system is stable we have plotted the parameters of (44), as function of wave number K 1. We observe two regions. The higher zone corresponding to stability and the lower zone corresponding to instability (see Fig. 1a, b and c). One obtains the same result with Tabi et al. [24]. We also plot the growth rate of instability and one finds a similar result with that obtained by Tabi et al. [24]. But it is necessary to note here a change of concavity obtained in Fig. 1a and e due to the presence in our system of another degree of freedom. This change of concavity explains the denaturation observed in an experimental way on the dynamics of DNA molecule. All particles of the molecular chain can vibrate in opposite direction, which creates the oscillation motions of double helix. Consequently the bases pair separates completely from the double helix. The regions of instability are regions where the initial plane wave is supposed to break into trains of soliton like objects. We have plotted on the Fig. 2, the behavior of the growth rate of instability according to the size of the nucleotide. When it is high the peak of instability decreases. The nucleotide size clearly fluctuates with instability.

Fig. 2
figure 2

Growth rate versus the wavenumber of the perturbation K 1 for m = 300a m u, a = 4.41Å −1, χ = 0.3e VÅ −1, D = 0.04e V, V = 0.0015e V, I = 300a m u and φ 0 = 10−3Å −1

Full size image

4 Numerical experiment

According to results obtained by many authors about the modulational instability, it is known that linear stability analysis can determine the instability field in space parameter and envisage qualitatively how the amplitude of a modulation sideband evolves at the onset of the instability. However, such an analysis is based on the linearization around the undisturbed wave carrier. Clearly, the linear approximation should fail at large time scales as the amplitude of an unstable sideband develops exponentially. Linear stability analysis therefore cannot indicate the long time evolution of a modulated nonlinear plane wave. In order to check the validity of our analytical approach and to explore the formation of localized modes, we carried out numerical simulations of our model by using the standard Fourier Transform method, with an integration time step of 0.055s. In our numerical simulations, the initial conditions are at time t = 0, are coherently modulated plane waves of the form: λ n = ψ n = λ 0[1 + 0.01cos(K n)]cos(K 0 n) and φ n = φ 0(1 + 0.01exp(i K n))exp(i K 0 n), where K = π/65, and K 0 = 3π/65 are the wavenumbers for the perturbation and the carrier waves , respectively. For our study, we use λ 0 = φ 0 = 0.002. Our aim in this numerical analysis is to point out the impact of the tight-binding hopping parameter on the occurrence of soliton-like objects induced by MI. In Fig. 3, we see that the initial condition tends to disintegrate during the propagation, leading to break-up of the wave into a pattern of wave trains. In fact, nonlinear interactions can give rise to very stable excitations, called solitons, which can travel without changing their shape. These excitations are very robust and important in the coherent transfer of energy and charges. If, in the first case, the patterns are constant, the second case shows how increasing V influences the distribution of wave patterns. One can conclude that MI is the most standard mechanism through which bright solitons or solitary structures can appear in physical systems. Dauxois et al. [22,23,24,25,26,27,28,29] have suggested that, such localized oscillations can be precursor of bubbles that appear in the thermal denaturation of DNA. They also showed that those structures could be used to describe the leading phenomena known as replication and transcription.

Fig. 3
figure 3

Spatiotemporal evolution of the amplitude of the initial plane wave which breaks into a wave train having the shape of a soliton-like object in DNA molecule, as predicted by the analytical predictions for m = 300a m u, a = 4.41Å −1, χ = 0.3e V Å −1, D = 0.04e V, V = 10−3 e V see Fig. 3a, and V = 15−3 e V see Fig. 3b, r = 3 Å −1, I = 300a m u and φ 0 = 10−3Å −1

Full size image

We observe in Fig. 4 the formation of localized structures thus translating the emergence of charges. These charges can migrate in all the directions of the wave propagation. The charges can occupy the pairs of bases, or can meet in a site given on the molecule of DNA, and give place to a collision (Fig. 4a–b). The charges are trapped within the molecule. For χ = 0.6 and χ = 0.8 one observes an electronic distribution of charges which migrates within the molecule. But starting from χ = 1.6 − 1.8, this charge distribution decreases and one observes small radiations coming from the sites. The charge localizes oneself and tends to form thin rows on fewer lattices. This result was suggested by Berashevich et al. [30] which considered the Peyrard-Bishop-Holstein model in the presence of an electric field. As a whole, increasing the coupling constant of charges reduces the instability domains and better enhances spreading, in terms of charge pattern. We note that there is explosion of domains of instabilities when the value of χ decreases; therefore, it is possible to obtain a chaotic situation when all the molecule is saturated by charges. The charge density decreases when χ increases and the number of areas of instability falls when χ increases.

Fig. 4
figure 4

Manifestations of modulational instability of charge transport in the DNA model under the influence of charge-vibrational and rotational coupling constant for m = 300a m u, a = 4.41 Å −1, D = 0.04e V, V = 0.1e V, r = 0.3Å −1, I = 300a m u, m = 300a m u, φ 0 = 10−1Å −1, with a χ = 0.6e VÅ −1, b χ = 0.8e VÅ −1, c χ = 1.4e VÅ −1, d χ = 1.8e VÅ −1

Full size image

5 Conclusion

We have investigated load spreading through MI in a DNA model. The linear stability analysis has shown that the model could be subjected to MI, as indicated by the representation of growth rate of instability. These analytical predictions have been verified numerically, where patterns of charge have been displayed. In this case, we have observed the region of stability/instability due to the nucleotide size. The impact of the coupling constant χ has been pointed out, and it has been found that increasing the value of χ better reduces the instability domains, and enhances charge spreading, in terms of charge pattern. We also showed the effect on size of nucleotides and the coupling transfer constant V. But other factors are not treated here, such as noise and thermal effect or internal vibrational excitations can also have a non-negligible influence on charge propagation.

  • The authors declare that they have no competing financial interests.

  • There are no competing interests related to this work.

  • There are no known conflicts of interest associated with this work.