Nonlinear dynamics of DNA systems with inhomogeneity effects
Abstract
We investigate the nonlinear dynamics of the PeyrardBishop DNA model taking into account site dependent inhomogeneities. By means of the multiplescale expansion in the semidiscrete approximation, the dynamics is governed by the perturbed nonlinear Schrödinger equation. We carry out a multiplescale soliton perturbation analysis to find the effects of the variety of nonlinear inhomogeneities on the breatherlike soliton solution. During the crossing of the inhomogeneities, the coherent structure of the soliton is found stable. The global shape of the inhomogeneous molecule is merged with the shape of the homogeneous molecule. However, the velocity, the wavenumber and the angular frequency undergo a timedependent correction that is proportional to initial width of the soliton and depends on the nature of the inhomogeneities.
keywords:
DNA, Nonlinear inhomogeneities, Breather soliton, Perturbation method.1 Introduction
DNA plays an important role in the carrier, protection, transmission, suppression, replication, transcription, recombination, repair and mutation of genetic information in biological systems. Several processes in the cell start with the binding of an protein enzyme at a promoter site of the DNA, this binding is known to change the conformation of DNA by generating a nonlinear localized excitation which causes the opening of few base pairs in the molecular chain sty . This excitation explains the transition conformation, the regulation of transcription, the denaturation and charge transport in terms of polarons and bubbles kal . Many theoretical models based on the longitudinal and transversal motions (as well as bending, stretching and rotations) have been proposed to describe the dynamics of DNA double helix Englander ; kong ; oka2 . One of the most popular of these models is the PeyrardBishop (PB) model in which DNA is considered as a helicoidal structure with a collection of particles connected with springs PB ; dau .
The biological processes which are executed in nature are carried out in a setting of inhomogeneities such as the protein enzymes which catalyse billion chemical reactions occurring anytime in a biological system. These processes are taking place not as isolated entities, but they do so in a molecular crowded environment. The particular biological functions of DNA imply the presence of different sites along the strands such as promotor (), coding (), several regulatory regions, (, , ), terminator (), which contain a specific sequence of base pairs and naturally make the strands inhomogeneous. The inhomogeneities in the DNA molecular chain can also be due to the presence of abasic sitelike nonpolar mimic of thymine or external molecules in the sequence or to the fact that the two different base pairs of the real DNA, and , combine in different ways constituting the genetic code lad . Along the same line, in their study on the dynamics of DNA with periodic and localized inhomogeneities both in stacking and in hydrogen bonds, using the plane base rotator model, Daniel et al. have shown that the inhomogeneity is found to modulate the width and velocity with which the open state configuration travels along the double helical chain. They also demonstrate that the inhomogeneity introduces fluctuations in the open state configuration represented by the kink/antikinktype soliton danD ; danA . Agüero et al., using the generalized coherent states approach to averaging the quasispin Hamiltonian for DNA, have found that, in a weakly saturating approximation, the formation of compacton/anticompacton pairs for the hydrogen bond displacements is strongly influenced by the inhomogeneity due to the action of an external agent on a specific site of the DNAagu .
The results obtained in the above studies demonstrate that, the interplay between nonlinearity and disorder in the DNA to biological processes is yet clearly far from being understood and rather encourage us to extend the study in another case: the case where the inhomogeneities are supposed to be due to the presence of additional molecules such as drugs, mutants, carcinogens or others in specific sites of DNA sequences along the molecular chain and cause damages or mutations on it. The mutations which are the accidental changes observed in the genetic code contains, occur constantly and are due to the actions of endogenous or external agents such as redoxcycling events involving environmental toxic agents and Fenton reactions mediated by heavy metals. Also reactive oxygen and nitrogen compounds produced by macrophages and neutrophils at sites of inflammation and infection arising as byproducts from oxidative respiration. The above chemical agents can attack DNA, leading to adducts that impair base pairing or block DNA replication and transcription processes, base loss, or DNA singlestrand breaks (SSB). Furthermore, when the DNAreplication apparatus encounters a SSB or certain other lesions, doublestrand breaks (DSB) are formed val . Cancer usually results from a series of mutations within a single cell. Often, a faulty, damaged, or missing gene is to blame. The gene provides proteins that stop mutated cells from dividing. Without this protein, cells divide unchecked and become tumours bie . The DNA lesions listed above can be gathered in three groups: the benefit mutations which create genetic diversity and keep population healthy, the silent mutations which have no effect at all (such diseases do not manifest) and the mutations which lead to diseases (base loss, SSB or DSB). Such DNA lesions are extremely toxic and difficult to repair. Then, the most promising directions in biophysics is the study of inhomogeneous nonlinear models of DNA, because this can give new interesting relations between the physical nonlinear properties of DNA and its biological functioning. Studies can lead to the discovery of the new mechanisms of regulation of fundamental biological functions of DNA, what can ”bridge” the nonlinear physics of DNA and medicine yaku1 . For instance, breatherimpurity interactions and its scattering in the PB model was studied in detail by Kyle et al. kyl , they show that the impurity can act as a catalyst and generates a larger excitation during the fusion of the breather mode due to the nonlinearity of the system and the one due to the mass inhomogeneity. With this in mind, in this work, we aim to study the nonlinear dynamics of the DNA double helix taking into account site dependence inhomogeneities, by considering the PB model. In this case, the DNA dynamics is governed by a perturbed nonlinear Schrodinger equation (PNLSE) and thus the problem boils down to solve this equation.
The paper is organized as follows. In Section 2, we propose the model Hamiltonian and derive the discrete equations of motion (EOM) for the inphase and outofphase motions respectively, using the multiplescale expansion in the semidiscrete approximation, the lattice equation is reduced to the PNLSE. In Section 3, the solitonic parameters and the first order soliton solution are obtained through the soliton perturbation technique. The effect of the inhomogeneity on the velocity, wavenumber, angular frequency, shape and position of the soliton solution is discussed. Section 4 concludes the paper.
2 Model Hamiltonian of DNA dynamics and equations of motion.
We consider the PB model PB for DNA denaturation where the degrees of freedom and associated to each base pair correspond to the displacements of the bases from their equilibrium positions along the direction of the hydrogen bonds that connect the two bases in a pair. A coupling between the base pairs due to the presence of the phosphate groups along the DNA strands is assumed to be inhomogeneous gra , so that the Hamiltonian for the model is given by:
(1) 
where and are the nucleotide mass and the elastic coupling constant in the same strand. The quantity represents the inhomogeneity site dependent character introduced in the transfer of the stacking energy between and base pairs. It shows that the stacking energy between neighbouring base pairs is site dependent function. As mentioned above, the inhomogeneity represents the intercalation of the compounds (drugs, mutants or carcinogens) between neighbouring base pairs along the DNA molecular chain without any distortion of the strands. The interactions between two bases in a pair are done by the hydrogen bonds which are modelled by the Morse potential given by:
(2) 
where is the depth of the Morse potential, is the width of the well. The values of parameters used to describe the motions of the two strands are those from the dynamical and denaturation properties of DNA. They are oka ; pey : amu, eV/Å, and Å. Our system of units (amu, Å, eV) defines a time unit () equal to seconds. The variables and are introduced:
(3) 
Taking into account Eqs. (2) and (3), the Hamiltonian of the system becomes
(4) 
where and represent the inphase and the outofphase motions respectively. By using the Hamiltonian of Eq. (4), the EOM of the system are given by:
(5) 
(6) 
The equation of variable describes the linear waves (phonons) while the one of variable describes the nonlinear waves (solitons). Hence, we restrict our attention on the second EOM (Eq. (6)) in which the small amplitude oscillation of the nucleotides is assumed around the bottom of the Morse potential allowing the following transformation pey .
(7) 
Replacing defined below into Eq. (6) and considering the system slightly inhomogeneous, our investigations will be limited to the analysis of the dynamical behaviour of the stretching motion of each base pairs, represented by the solution of Eq. (6). We assume , where is a sitedependent function which measures the inhomogeneity in the stacking. is treated perturbatively by assuming that, it contributes little enough to the whole DNA dynamics which is dominated by the first and second terms of the following equation obtained up to the third order of the Morse potential,
(8) 
where , , and . Equation (8) is the equation describing the dynamics of the outofphase motion of a weakly inhomogeneous DNA model. It can be solved directly using the series expansion unknown function method zdr or the exponential rational function method tala . In this paper, we are looking for the envelope soliton in the small amplitude approximation as a perturbed plane wave solution. For this purpose, it has been shown that the multiplescale expansion in the semidiscrete approximation is the most adapted technique oka2 ; pey ; rem . This approximation is a perturbation technique in which the amplitude is treated in the continuum limit while the carrier waves are kept discrete. The technique allows the study of the modulation of the wave, thus the soliton solution is looking in the form
(9) 
where C.C. stands for complex conjugate and . , and represent the optical frequency of the linear approximation of the base pair vibrations, the distance of neighbouring bases in the same strand and the wavenumber respectively. Nonlinear terms in Eq. (8) incite one to predict that, through frequencies superpositions, the first harmonics of the wave will contain terms in as well as terms without any exponential dependence. The amplitudes , and will be equally considered to change slowly in space and time. For these purposes, the continuum limit approximations and the multiplescale expansion will be applied on those amplitudes. Thus the amplitudes are treated as a function of variables according to the new space and time scales and respectively. Hence the solution which depends on these new sets of variables is found as a perturbation series of functions. We will consider here that . By Taylor expansion and up to the second order in , we obtain for spatial derivatives
(10) 
and their temporal derivatives
(11) 
The same procedure is used for and . Using Eqs. (9), (10) and (11) together with Eq. (8) and collecting the coefficients for the different powers of (, ), one obtains the angular frequency and the group velocity of the wave
(12) 
and
(13) 
respectively. The functions and in Eq. (9) can be expressed through as
(14) 
with
(15) 
where , while is a solution of the following PNLSE
(16) 
After changing to the frame moving at the group velocity of the carrier wave by defining and , we get:
(17) 
where
(18) 
In the following, we rescale the variable as and define a new function such as . Inserting the above considerations in Eq. (17), we get the PNLSE in the form
(19) 
where is assumed to be the perturbation parameter and the amplitude of the inhomogeneity.
It should be noted that Eq. (19) is the PNLSE, regulating the dynamics of the envelope soliton which represents the outofphase motion of DNA in the presence of inhomogeneity. In the case of a system without inhomogeneity in the lattice () or when the wave vibrates at the frequency close to (), the above equation is reduced to the following standard Nonlinear Schrödinger equation (NLSE)
(20) 
with the well known soliton solution jia ; danE ; karp
(21) 
where , , , are four real parameters which represent the height (as well as the width), the velocity, the initial position and initial phase of the propagating soliton respectively.
3 Effect of site dependence inhomogeneity on the open state
In this section, the soliton perturbation technique (SPT) is used to construct the firstorder perturbed soliton solution of Eq. (19) and to find the modifications on its parameters due to the inhomogeneity. The technique allows to consider the inhomogeneity as a perturbation term due to the fact that, is a small positive constant measuring the weakness of the perturbation.
Following the work done by jia ; danE ; karp , Eq. (19) is linearized by transforming the independent variable into several variables
(22) 
where the subscripts stand for partial differentiation with respect to the time .
It is more convenient to represent everything in the coordinate system moving with the soliton. Then, we use as a new space independent variable in place of . Under this assumption of quasistationary and due to the new time scale introduced below, the soliton parameters , , , and are now supposed to be function of the slow time variables . and are independent of jia ; danE ; karp . The one soliton Eq. (21) solution of the standard NLSE is rewritten here for convenience in the form:
(23) 
with
(24) 
By assuming that the wave propagates at the velocity close to the quarter of the group velocity, Eq. (24) in Eq. (19) leads,
(25) 
where .
We assume the solution to be on the form:
(26) 
where is the amplitude of the envelope soliton. Hence, Eq. (26) in Eq. (25) gives the linearized PNLSE in the form:
(27) 
with
(28) 
We use the poincaré expansion type
(29) 
Inserting Eq. (29) into Eq. (27) and equating the coefficients of each power of , we obtain the following equations at different orders of : At the order , we get the following equation known as stationary NLSE
(30) 
where ;
At the order , we get the stationary PNLSE in the form:
(31) 
where means the complex conjugate.
The term given in Eq. (29) is assumed to be on the form:
(32) 
Introducing Eq. (32) into Eq. (31) gives:
(33) 
with
(34) 
and
(35) 
and are two selfadjoint operators.
3.1 Variation of the solitonic parameters
In this part, the variation of the parameters of the soliton is evaluated by assuming that, the soliton propagate with an amplitude (as well as a width) and a velocity when the perturbation is switched off. To evaluate the modifications of the solitonic parameters, it is necessary to solve the homogeneous parts of Eq. (33) which admit and as solutions respectively. The nonsecularity conditions jia ; karp give us the following four important formulas which determine how the solitonic parameters (i.e. the amplitude, the position, the velocity and the angular frequency) are modified by the inhomogeneity.
(36) 
(37) 
(38) 
(39) 
with
(40) 
where . and are the real and imaginary parts of given in Eq. (25).
To find the variation of the soliton parameters explicitly, we have to evaluate the integrals found in the righthand sides of Eqs. [(36)(39)] which can be carried out only on supplying the specific forms of . Hence, we consider the localized inhomogeneity in the form of the hyperbolic tangent function and the exponential inhomogeneity in the form of exponential function separately. The results give the time dependence of the parameters of the soliton as:
(41) 
The first equation of the Eq. (41) shows that the amplitude (as well as the width) of the soliton remains constant, showing that the number of base pairs which take part in the opening process remains constant during the propagation. We observe in the second equation that the position of the soliton is not affected by the inhomogeneities. However, from the third and fourth equations, we found that the velocity and the angular frequency of the soliton can get a correction depending of the type of the inhomogeneities in the lattice.
The localized inhomogeneity can represent the intercalation of a compound between neighbouring base pairs or the presence of a defect on an abasic site in the DNA chain. In order to understand the effects of this type of inhomogeneity in the lattice, we substitute , in Eq. (41), and after integration we obtain and , which can be written in terms of the original time variable by using the expression and as:
(42) 
Equation (42) demonstrates that, the localized inhomogeneity in the above form affects the velocity and the angular frequency of the soliton, they get a correction during the crossing of this inhomogeneity. The nature of the correction depends of the nature of the inhomogeneity represented by the sign of its amplitude which can be either positive or negative. When is positive (), the inhomogeneity corresponds to an energetic barrier and on the other hand, when is negative (), the inhomogeneity behaves as a potential well danA .

When greater than zero, the correction is positive, what makes the velocity and the angular frequency increasing during the crossing of the inhomogeneity. The amount of these increments is proportional to the initial height (as well as the initial width) of the soliton. The higher the initial amplitude of the soliton, the greater the increments. The increasing in the velocity of the soliton helps to overcome the barrier due to the inhomogeneity and the soliton will propagate easily along the inhomogeneous DNA chain without formation of a bound state. Reporting these observations in Eq. (23), we notice that during the crossing of the inhomogeneity, the soliton propagates and vibrates quickly, and the DNA breathing mode becomes faster.

In the case of less than zero, the correction is negative what makes the frequency and the velocity decreasing. The decreasing of the velocity slows down the soliton. The soliton stops and vanishes when the original time satisfies the condition , where is a constant which depends on the inhomogeneity. The propagating time of the soliton depends on its initial velocity and width. The wider the soliton, the greater the propagating time.

When , the velocity and the angular frequency of the soliton remain constant.
In all the above cases, the soliton which is a coherent structure formed by involving a few base pairs, moves along the helical chain in the form of bubble without dissipation or any other form of deformation is found stable.
Next, we consider . This case corresponds to the physically interesting problem of the exponential distribution of similar molecules along the DNA helical chain. We have and . Written in the original time variable we get
(43) 
The same observations as in the previous type of inhomogeneity are also done here. We notice from the results that considering the localized and exponential inhomogeneities, the propagating velocity and frequency of the soliton undergo a correction during the crossing of the inhomogeneities. In Figures (1) and (2), the timeevolution of the velocity and angular frequency of the soliton are depicted as a function of the type and nature of the inhomogeneities in the lattice. The figures show that the correction terms can be positive or negative according to the nature of the inhomogeneities.
One comparing Eq. (42) and Eq. (43), we notice that, the absolute value of the correction terms is greater in the case of exponential inhomogeneity as can be seen in the above figures. This is because in this case, the inhomogeneities occur exponentially in the entire chain of the DNA molecule in terms of external agents.
3.2 Firstorder perturbed soliton solution
In this part of our work, we will seek for the breatherlike soliton solution for Eq. (25). As we know, once the seed solution is chosen as breather soliton solution via Eq. (23), the firstorder perturbed breatherlike soliton solution can be constructed using the SPT by solving Eq. (33) and report the solutions and in Eq. (32). We firstly solve the homogeneous part of the first equation in Eq. (33), which admits the following two particular solutions:
(44) 
and
(45) 
Then, the general solution can be obtained by using the following formula:
(46) 
where and are two arbitrary constants to be determined. We construct the solution , by substituting the expressions of , and in Eq. (46). After the evaluation of the integrals, we obtain:
(47) 
where and are the integrals depending on the form of the inhomogeneities and are given by:
(48) 
The term in is a secular term which makes the solution unbounded. The above term can be removed if we assume that the arbitrary constant . Using the boundary conditions and , we obtain and . By using the above results in Eq. (47), we get the general solution in the form:
(49) 
Now we give the solution of the second equation in Eq. (33), its homogeneous part admits the particular solutions given by:
(50) 
while the general form of is given by:
(51) 
where and are two arbitrary constants to be determined. We follow the same procedure by removing the secular terms (), we get also and applying the following boundary conditions , we have
(52) 
Taking into account the above results, the general form of the envelope soliton is
(53) 
where and depend on the type of inhomogeneity and are given in Eq. (48). The effects of the inhomogeneities on the shape of the first order perturbed soliton is therefore reduced to the integration of and .
We construct the first order perturbed soliton solution by substituting the corresponding values of , for respective given below and evaluating the integrals and , which involve very lengthy algebra. Thus, in the case , we use the values of and of Eq. (42). The Eq. (53) becomes,
(54) 
We then repeat the same procedure for constructing the first order perturbed soliton solution in the case of and we obtain the perturbed soliton solution as:
(55) 
Eq. (54) (and Eq. (55)) is the breatherlike soliton which represents the envelope of the soliton solution. Its describes the open state configuration in the individual strand of the DNA, which collectively represents a bubble moving along the inhomogeneous DNA molecule.
We have depicted in Figure 3 the schematic representation of the square of the envelope soliton as a function of the type of inhomogeneities in the lattice. Figure 3a presents the breather soliton moving in the homogeneous DNA chain. The case of localized inhomogeneity is depicted in Figure 3b while the case of exponential inhomogeneity is depicted in Figure 3c. We notice that, in both inhomogeneous cases, the robust nature of the breatherlike soliton is not modified as it propagates along the DNA molecule. The global shape of the molecule with inhomogeneities is merged with its shape without taking into account the inhomogeneities. Similar results were observed by Mvogo et al. in their investigations on the effects of localized inhomogeneity in the form of hyperbolic secant function in the helical proteins chain mvo . Also we found in Ref. saha that, using the numerical methods, the authors demonstrate that neither the opening of base pairs in DNA molecule nor the topological character of the breatherlike soliton are affected by the localized inhomogeneity in the form of hyperbolic secant function introduced in stacking and hydrogen bonding energies of DNA.
Considering Eq. (7) together with Eq. (9), we assume that the configuration of the molecule can be fully described by the following soliton solution with a rescaled time
(56) 
where
(57) 
and are the ”effective” wavenumber and angular frequency of the soliton solution respectively. , and represent the envelope soliton, its velocity and angular frequency respectively. They depend on the type of inhomogeneity and are given by [Eq. (42) and Eq. (54)] and [Eq. (43) and Eq. (55)] for the localized and the exponential inhomogeneities respectively. Eq. (56) represents the stretching of the base pairs also known as the breathing modes experimentally observed in DNA molecule putn . Since the wave velocity depends on the type of inhomogeneity in the lattice, from the first equation of Eq. (57) we notice that the wavenumber is a function of the inhomogeneities, it changes as time progressed due to the inhomogeneities.
The configuration of the molecule and the elongation of the outofphase motion vs. time are depicted as a function of the type and nature of the inhomogeneities present in the lattice. Figures 4 and 5 represent the configuration of the molecule while Figures 6 and 7 present the elongation of the outofphase motion vs. time. As predicted by Eqs. (42), (43) and (57), we observe that the wave velocity, the wavenumber and the angular frequency of the soliton solution get a correction which can be positive or negative depending on the nature of the inhomogeneities:

For , the inhomogeneity behaves as an energetic barrier. The above parameters increase as time progresses due to the inhomogeneities. Figure 4 shows the increasing of the wavenumber and wave velocity (Figure 4d) while Figure 6 shows the increasing of the angular frequency. The increasing of the velocity helps the soliton to overcome this energetic barrier. Therefore, the soliton can propagate easily without formation of a bound state.

For , the decreasing of the above wave parameters is observed in Figure 5 for the wavenumber and in Figure 7 for the angular frequency. In this case, the inhomogeneity behaves as a potential well in which the soliton is trapped. Thus the soliton slows down, stops and vanishes at where is equal to and for localized and exponential inhomogeneities respectively.

For the inhomogeneity is switched off. This case corresponds to the soliton travelling in a homogeneous media with the constant velocity, wavenumber and angular frequency.
Also, as in Figures 1 and 2, we observe that the absolute value of the corrections in the velocity and wavenumber (Figure 4 and Figure 5), and in the angular frequency (Figure 6 and Figure 7) of the soliton is greater in the case of exponential inhomogeneity. We also notice that the inhomogeneities do not affect the bubble sizes. The height and the width of the bubble remain constant during its propagation in the inhomogeneous DNA molecule.
It is clearly known that the transcription process is a very biological complex phenomenon yang . Thus the increasing of the velocity, the wavenumber and the angular frequency of the soliton during the crossing of the inhomogeneities makes difficult the execution of the biological functions of DNA such as reading, transcription and recombination of the genetic code to mention a few. These difficulties undoubtedly involve errors which can rise to mutations, damaging, missing of genes and biological diseases such as: sickle cell anemia, heart diseases, high blood pressure, Alzheimer’s disease, diabetes, cancer, obesity, eye disease, epilepsy or strokelike episodes which are due to mutations or disorder caused by a combination of environmental factors and mutations.
4 Conclusion
The wave dynamics of the PB model of DNA at the physiological temperature was studied. The classical PB model, assume the bases to be equal. In our model we include the differences of bases and sequences through site dependent springconstant. Using the multiplescale expansion method in the semidiscrete approximation, the outofphase motion has been described by the PNLSE. In the homogeneous limit, the dynamics is governed by the breather soliton of the integrable NLSE. Indeed, to understand the effect of inhomogeneities on the base pairs opening, we carried out a perturbation analysis using multiplescale SPT. For implementing this, we linearized the PNLSE with the aids of poincarétype asymptotic expansion. The breatherlike soliton solution which represents the opening of base pairs travelling along the inhomogeneous DNA chain in the form of bubble is constructed for different forms of the inhomogeneities (localized and exponential).
The results showed that the bubble profile is not affected by the inhomogeneities. However the velocity with which the base pairs are opening or closing, the wavenumber with which the soliton propagates and angular frequency with which the base pairs vibrate increase, decrease or remain uniform and even the soliton stops depending on the nature of the inhomogeneities. The high velocity, wavenumber and angular frequency of the soliton can create disorder in the execution of DNA biological functions and make coding or reading errors which can leads to a variety of biological diseases.
The inhomogeneous DNA models can explain DNA biological functions such as replication, transcription and recombination more viably than the homogeneous one, since it can explain and predict the biological diseases due to genetic mutations. The number of DNA damages in a single human cell exceeds every day klu and must be counteracted by special DNA repair processes. Despite the relevance of this work, the impact of inhomogeneities on the DNA biological processes is not yet clearly understood, since the nature generally selects inhomogeneous DNA and the gap between the nonlinear physics of DNA and medicine is still big.
Acknowledgements
J.B. Okaly is in debt to Dr. Ndzana Fabien II of the University of Maroua, Maroua Cameroon for some recommendations and fruitful discussions.
References
 (1)