A Partial Averaging of Lagrangian Equations is Equivalent to Partial Averaging of its Hamiltonian for the Reduced Model

Control of a Model of DNA Division via Parametric Resonance1

Abstract

We study the internal resonance, energy transfer, activation mechanism, and control of a model of DNA division via parametric resonance. While the system is robust to noise, this study shows that it is sensitive to specific fine scale modes and frequencies that could be targeted by low intensity electro-magnetic fields for triggering and controlling the division. The DNA model is a chain of pendula in a Morse potential. While the (possibly parametrically excited) system has a large number of degrees of freedom and a large number of intrinsic time scales, global and slow variables can be identified by (i) first reducing its dynamic to two modes exchanging energy between each other and (ii) averaging the dynamic of the reduced system with respect to the phase of the fastest mode. Surprisingly the global and slow dynamic of the system remains Hamiltonian (despite the parametric excitation) and the study of its associated effective potential shows how parametric excitation can turn the unstable open state into a stable one. Numerical experiments support the accuracy of the time-averaged reduced Hamiltonian in capturing the global and slow dynamic of the full system.

In this paper we study the internal resonance, energy transfer, activation mechanism, and control of a model of DNA division via parametric resonance. While DNA macro-molecules are robust to noise, our study shows that they are sensitive to specific fine scale modes and frequencies that could be targeted by low intensity electro-magnetic fields for triggering and controlling the division. The suggested method of control is supported not only by the observation that DNA vibrations induced by electric-fields or microwave absorbtion are an experimental reality but also by the fact that electric field-induced molecular vibrations have already been used as a noninvasive cell transfection protocol. Our study also raises the question on whether enzymes are using the proposed mechanism to initiate the opening of DNA strands. This question is to put into correspondence with increasing theoretical and experimental evidence that low-frequency vibrations do exist and play significant biological functions in proteins, DNA molecules, and other bio-macromolecules.

1 Introduction

The model considered in this paper is a chain of pendula in a Morse potential, with torsional springs between pendula, that mimic real DNA characteristics [24, 38]. Previous studies [24, 8, 9], mainly numerical, showed that this model exhibits an intriguing phenomenon of structured activations observed in many bio-molecules: while the system is robust to noise, it is sensitive to certain specific fine scale modes that can trigger the division. Below, we will describe briefly the results of our analytical study on this intriguing phenomenon and our effort in the control of this DNA model via parametric resonance.

By using the Fourier modal coordinates [8, 12, 23], this model can be seen as a small nonlinear perturbation of harmonic oscillators with frequencies . The coarse variable of the (approximate) th mode, which is the average angle of the pendula, corresponds to the angle of the frame defined by the radius of gyration in our molecular studies [40, 41, 42]. This variable is a global and slow variable and it plays an important role as the reactive coordinate for our DNA work. Moreover, this reactive mode forms a nearly resonance with any other mode, each of which has an frequency. This fact leads to small denominators or coupling terms in the corresponding averaged equations or normal forms [26, 27, 28, 29, 18, 14, 35, 20, 2]. Since other modal frequencies are not rationally commensurate or have significant time scale separation, we do not expect strong resonance among them. Extensive simulations confirmed our expectation. We observed: (i) the energy transfers mainly from an excited mode to the reactive mode, triggering the division, (ii) only an extremely small amount of energy transfers from the excited mode to one or two other modes via near resonance. This observation, together with a rigorous error estimate [8], show that two-mode or three-mode truncation, i.e., a truncated system that includes the reactive mode, the excited mode, and perhaps a mildly affected mode should provide an adequate reduced model for our analytical study on the activation mechanism.

By applying the method of partial averaging [12, 1, 26, 28, 18, 14, 35] for nearly resonance, we obtained the averaged equations for the reduced models of this chain of Morse oscillators up to nonlinear terms of very high degree. These averaged reduced equations not only reveal the coupling between the action of the excited modes and the dynamics of the reactive mode, they also shed lights on the phase space structure of the activation mechanism. Moreover, they allow us to estimate analytically the minimum activation energy for each excited mode. These analytical estimates not only match well with those obtained from simulations of the full DNA model, but also uncover a relationship between the frequency of the excited mode and its corresponding minimum activation energy. These findings also provide an analytical and deeper understanding of the internal mechanism that is responsible for the phenomenon of structured activations.

Based on our understanding of its internal dynamics, we introduce a method for controlling the division of this DNA model via parametric excitation, that is in resonance with its internal trigger modes. By choosing appropriate external excitations and frictions, we uncover a class of trajectories that show how the parametric resonance can be used to drive the averaged reduced model from its (almost) equilibrium state to its open state. The identification of an effective Hamiltonian (after reduction and partial averaging) opens the possibility of studying the global phase space structures of our averaged reduced model and sheds lights on the significant trajectories mentioned above. Moreover, we extend the results for the averaged reduced model with parametric resonance and friction to the reduced as well as the full model with parametric resonance and friction. These findings support the conjectures that (i) low intensity electro-magnetic fields can be used with parametric resonance to inject energy into the trigger modes and destabilize the DNA chain and keep it near the open state for replication and transcription, (ii) enzymes may use similar method to initiate open state dynamics for DNA replication and transcription. Although the issues of inhomogeneity, helicity, and environmental effects (such as noises) will be ignored for now, they will be taken into consideration in our future work.

2 A Model of DNA Division and Structured Activations

Our analytical study has been inspired by the work of Mezić and Eisenhower at UC Santa Barbara and Marsden and Du Toit at Caltech. Their studies [24, 8, 9, 11], mainly numerical, showed that this DNA model exhibits an intriguing phenomenon of structured activations. While Eisenhower and Mezić did apply Arnold’s method of partial averaging [1] to a truncated Hamiltonian in their analytical study of a chain of Duffing oscillators [10, 11, 12], they did not extend it to the DNA model. The reason that they gave in [12] for studying the Duffing case was that “The exponential form of the Morse potential makes analytical progress difficult,…” Since we have not seen this kind of simplification in the DNA literature, we will, in this study keep the Morse potential (and the added difficulty).

2.1 A Model of DNA Division

The model was first introduced in [24]. It is a chain of equivalent pendula that rotate about the axis of a fixed backbone with an angle measured from the upward position. The pendula interact with nearest neighbors along the backbone through harmonic torsional coupling, and with pendula on the opposing immobilized strand through a Morse potential that has two stable equilibria and a saddle. See Figure 1.

Figure 1: (a) A model of DNA division with 30 pendula. (b) The phase space of a single pendulum in a Morse potential without coupling. It has two stable equilibria and a saddle. The black curve is the homoclinic orbit that separates two types of motion, namely, the libration near the equilibria , and the flipping across the saddle .

The Hamiltonian that describes the motion of these coupled pendula is given by

(1)

with periodic boundary condition, . The first term is the kinetic energy terms of -pendula. The second term is the torsional coupling terms. The third term is the Morse potential terms, which model the hydrogen bonds of the respective DNA base pairs. In Eq. (1), and represent the mass and the length of each pendulum, and is the generalized momenta conjugate to . The parameter determines the strength of the nearest neighbor coupling, while the parameter determines the strength of the Morse potential. The parameter determines the equilibrium distance of the Morse potential, while the parameter determines the width of the Morse potential. All the parameter values were chosen to best represent the typical values for the opening and closing dynamics of DNA division [4, 22] as follows: amu (typical mass of a DNA base (nucleotide)), nm (typical radius of DNA), eV, eV, nm, and nm.

After dividing the both sides of Eq. (1) by , the Hamiltonian can be non-dimensionalized as

(2)

where is the dimensionless momentum defined with respect to the dimensionless time . Thus, in the present study, one unit time () corresponds to 0.272 ps. In Eq. (2), the dimensionless amplitude of Morse potential is a small parameter and is equal to . We have also introduced the dimensionless decaying coefficient of Morse potential , and the dimensionless equilibrium distance .

For our analytical study, we use, as in previous numerical studies [24, 11], a Hamiltonian system composed of 30 coupled pendula. Figure 1(a) is a model of 30 coupled pendula. Before studying its dynamics, it is instructive to look at a single pendulum in a Morse potential without coupling. Figure 1(b) shows the phase space of such single pendulum. It has two stable equilibria and a saddle. The black curve is the homoclinic orbit that separates two types of motion, namely, the oscillation near the equilibria , and the flipping across the saddle . The n-coupled pendula have similar but much more complicated behaviors. First, the system has two global stable equilibria, achieved when all the pendula have the same angular displacements (thus nullifying the nearest neighbors coupling) and each is positioned at the equilibria of a single pendulum. It also has a rank one saddle at where all pendula are at the upward position. For small energy, the pendula are liberating near one of the global stable equilibria where all angles are the same and equal to . For large enough energy, it has been observed that a local activation can cause the pendula to move collectively from one energy basin to the other and to flip across the rank one saddle.

2.2 Phenomenon of Structured Activations

Previous studies [24, 8, 9, 11], mainly numerical, showed that this model exhibits an intriguing phenomenon of structured activations observed in many bio-molecules: while the system is robust to noise, it is sensitive to certain specific fine scale modes that can trigger the division. Figure 2 provides the numerical data for such claim. The figure shows the relationship between the initial amount of energy injected for various types and shapes of activation and the time for DNA division.

Figure 2: Figure shows the relationship between the initial amount of energy injected for various types and shapes of activation and the time for DNA division. Reprinted with permission from B. Eisenhower, PhD Dissertation, University of California Santa Barbara, (2009) [11].

In order to appreciate this figure and the phenomenon of structured activations, we need to introduce the Fourier modal coordinates which relate to the original system coordinates through the following linear transformation ():

(3)

where and . Here, we have assume that the number of pendula is even. If is odd, we merely need to have the middle column corresponding to removed and the column altered accordingly.

These modal coordinates help to reveal the natural dynamics of the system by diagonalising the linear coupling terms and rewriting the Hamiltonian as follows:

(4)

where is the Morse potential for a single pendulum. More specifically, the matrix is nothing but the matrix of orthonormal eigenvectors used for the diagonalization and are their corresponding eigenvalues. By using this coordinate system, the model can be seen as a small perturbation of harmonic oscillators with frequencies , This can also be seen clearly if we write the equation of motion in the Lagrangian form:

(5)

where

(6)

is the Morse potential term, and are the partial derivatives of with respect to respectively.

Observe that the coordinate of the (approximate) th Fourier mode, given as follows

(7)

is the average amplitude of the pendula except for a constant factor of . It plays a special role as the collective variable, the reactive coordinate, and the slow variable for the system. Other (n-1) modal coordinates are the fine scale variables, the bath coordinates, and the fast variables.

Figure 3: Figure shows a sequence of six snapshots of the evolution of 30 pendula from one equilibrium to the other equilibrium . Because the coupling is much stronger than the nonlinearity, the transition is collective, closely follow the mean, the thick blue line.

The role of approximate th mode (or reactive) coordinate can be seen from the sequence of six snapshots of the evolution of 30 pendula from one equilibrium to the other equilibrium of the Morse potential for a single pendulum. For clarity of illustration, only one pendulum is perturbed as the initial activation. Because the coupling is much stronger than the nonlinearity, the transition is collective, closely follow the mean, the thick blue line. Therefore, the average angle of pendula can be used to monitor and mark the time of DNA division.

Now we are ready to appreciate Figure 2. The figure shows the relation between the initial amount of energy injected for various Fourier modes of activation and the time for DNA division. The curves have been constructed by choosing the initial activation as a single Fourier mode. Take the magenta curve of the 4th mode as an example. Assuming that the initial position of the system is at its equilibrium point , certain amount of the momentum of the 4th mode is injected for the initial activation at and its amplitude could be modified to vary the amount of activation energy . After the system evolves for awhile, the time of division could be determined as the time when the average angle first crosses the position where . In this way, a curve of is obtained that shows the amount of activation energy vs the time to DNA division for each Fourier mode. The integration time is set for 3000 units which is approximately equal to . Notice that there are 14 such curves from left to right representing those from 0th mode to 13 mode respectively. Each curve has an asymptote at the low energy limit and will be named the minimum activation energy for each excited mode. White “”s show the data for the random noise. From the figure, we can see that the minimum activation energy depends on the way this energy is injected into the system. While the system is robust to noise, it is sensitive to certain specific fine scale modes that can trigger the division. There may also exist a relation between the minimum activation energy of each Fourier mode and its modal frequency. In this paper, we would like to develop the analytical methods to reveal the activation mechanism and to compute the minimum activation energy for each mode. Moreover, we want to develop the techniques for controlling the real DNA division via low intensity electro-magnetic fields and to reveal how enzymes initiate the DNA opening dynamics.

3 Analytical Study of Structured Activations

Let us start our analytical study of the phenomenon of structured activations.

3.1 Nearly Resonance and Partial Averaging

Recall that the equations of motion of our system are given by Eq. (2.2) which can be seen as a nonlinear perturbation of harmonic oscillators with frequencies Note that besides , varies from 0.2091 to 2. Hence, the reactive mode forms a nearly resonance with any other mode, each of which has an frequency:

(8)

While such a relationship may not be seen as a resonance in the classical “engineering sense”, it certainly is in the mathematical sense. This fact leads to small denominators and coupling terms in the corresponding averaged equations or normal form. Other modal frequencies, from 1st to 14th mode (and from 15th to 29th mode), are not rationally commensurate and do not have significant time scale separation. We do not expect strong resonance among them. Hence, we believe that nearly resonance should be the main focus of our study.

Nayfeh et al. [26, 28], Haller [18], Feng and Liew [14], Tuwankotta and Verhulst [35], Langford and Zhan [20] and Broer et al. [2] have studied such a degenerate resonance. Both Langford and Zhan and Broer et al. used the method of normal form. Nayfeh et al., Haller, and Feng and Liew applied a modified averaging method directly to a two mode truncation of a simple mechanical system with parametric or external excitation. Tukwankotta and Verhulst applied a similar method to the study of nonlinear wave equations.

While Eisenhower and Mezić did not mentioned the term, nearly resonance, in their papers, they did apply Arnold’s method of partial averaging [1] to a truncated Hamiltonian in the study of a chain of Duffing oscillators. In [12], the reason that Eisenhower and Mezić gave for studying the Duffing case was that “The exponential form of the Morse potential makes analytical progress difficult,…” Since we have not seen this kind of simplification in the DNA literature, we will, in this study keep the Morse potential (and the added difficulty).

3.2 Two Mode Truncation Is Adequate

In order to carry out the analytical work for such a high dimensional system, certain reduced model is needed. The success of Eisenhower and Mezić [12], Du Toit et al. [8], Nayfeh et al. [26, 28], and other researchers [23, 18, 35] shows that the method of modal truncation is a viable method if it is used with care. Our understanding of the nearly resonance in our DNA model and the results of our extensive numerical simulation have convinced us that two mode truncation is adequate for our analytical study of the phenomenon of structured activations.

Figure 4: The figure shows the projections of a sample trajectory on the phase spaces of the first 15 modes. is the amount of initial kinetic energy injected into the 6th mode.

Figure 4 shows the projections of a sample trajectory on the phase spaces of the first 15 modes. The activation has been chosen to be a single 6th mode. From the numerical simulation which use and as the initial condition (where is the amount of energy injected into the 6th mode), we observe that (i) the energy transfer takes place over-whelmingly from the excited mode to the reactive mode, inducing the division, (ii) only an extremely small amount of energy transfers from the excited mode to one or two other modes via near resonance. Hence, we expect that the two-mode or three-mode truncations i.e., a truncated system that includes the reactive mode, the excited mode, and perhaps a mildly affected mode, should provide an adequate reduced model for our analytical study on the nearly resonance, the energy transfer, and the activation mechanism.

What follow are the equations for any two-mode truncation (the reactive mode and the th mode)

(9)

where

(10)

is the reduced Morse potential term, and are its partial derivatives with respect to respectively. Note: for notational simplicity, will be used later for whenever there is no ambiguity.

Error Estimates for the Two Mode Truncation.

Here we show that the solution trajectories of our two mode truncation Eq. (3.2) are within of the solution trajectories of the original full system described in Eq. (2.2) for at least times.

The proof is an extension of the one in Du Toit et al. [8], where a general system that has the same form as Eq. (2.2) was studied. In that paper, the authors (i) proposed a general technique for obtaining an degree of freedom reduced system and (ii) were able to provide a rigorous error estimate for their procedure. First, they introduced an approximation by replacing Eq. 2.2(b) with the analytical solutions of the unperturbed linear system, as defined by

(11)

where are the initial conditions. Then, they obtain their reduced equation

(12)

via simple substitutions. Hence, the information contained in the other modes persists in the reduced equation of the reactive mode via the initial conditions. The error estimate that they provided claims that the solution trajectories of Eq. (12) are within of the solution trajectories of the original full system described in Eq. (2.2) for at least times. This result is shown by applying a standard error analysis technique: substituting a formal expansion of the solutions, using the Lipschitz continuity of , and then applying the Gronwall lemma as is done in the proof of Theorem 9.1 in [37] (see also [34] for a more detailed analysis on why removing small nonlinear perturbation to harmonic oscillations incurs small error).

Notice that the equations of our two mode truncation Eq. (3.2) can be obtained by simply setting all the initial conditions for Eq. (2.2) equal to zero except if . Hence, the solution trajectories of our two mode truncation Eq. (3.2) should also be within of the solution trajectories of the original full system described in Eq. (2.2) for times .

Moreover, this error estimate is valid for an arbitrary forcing function. But for our DNA model, the exponential decay of the Morse potential with the distance implies that when the pendula escapes the immediate vicinity of the opposing pendula, the Morse potential and its consequent perturbation are effectively zero. Hence, our error bound is very loose because it utilizes only the fact that is small, but not the specific feature of our DNA model, which is that the forcing term is also small in a large region of phase space. This is why we numerically observed that the two mode truncation remains accurate over a timescale much longer than .

Figure 5: Figure shows the amount of activation energy vs the time of DNA division for the 6th mode. The blue “o”s are the data points for the full equations of 30 modes; the red “*”s are correspondent data points for the two mode truncation. Both have almost the same minimum activation energy .

Figure 5 show the amount of activation energy vs the time of DNA division for the 6th mode. The blue “o”s are the data points for the full equations of 30 modes; the red “*”s are correspondent data points for a two mode truncation (the reactive mode and the 6th mode). Both have almost the same minimum activation energy . Even their time of division differs very little if the actuation energy is slightly larger than the minimum activation energy (away from the asymptote). The above observation also holds true for all the other modes.

Now we are ready to do the truncation and apply the method of partial averaging to obtain the averaged equations for the reduced models of this chain of coupled Morse oscillators, and use them to study the activation mechanism and compute the minimum activation energy for each mode.

3.3 Averaged Reduced Hamiltonian

Since in our case, the partial averaging of Euler-Lagrangian equations is equivalent to the partial averaging of its Hamiltonian and our main concern is on energy transfer, we preferred to use the Hamiltonian formulation in this study. The equivalency was shown in §A of the Appendix.

Recall that the Hamiltonian of the two mode reduced model is given by

(13)

The averaged reduced Hamiltonian can be obtained as follows. First, we approximate (via Taylor expansion) the Morse potential , which involves the exponential function, by a polynomial of degree 26 at

(14)

where are the coefficients of the expansion. As pointed out in [12], the Morse potential is indeed a difficult function to work with. Since the geometry of the Morse potential requires a polynomial of degree at least for an accurate approximation, we decide to use such a high degree expansion for our computation of the averaged reduced Hamiltonian.

Then, we use the angle-action coordinates defined as follows

(15)

and rewrite the reduced Hamiltonian as . Notice that, besides , the action is also a slow variable. Hence, the averaged reduced Hamiltonian can be obtained by averaging the only fast variable :

(16)

After renaming variables, the averaged reduced Hamiltonian is given by

(17)

where are the Cartesian coordinates of the reactive mode; are the action and the constant frequency of the other mode; is the energy value at the saddle; are polynomials in . For example, is a 13 degree polynomial in given by

(18)

where are numerical coefficients. Also, for the simplicity of notation, the letter is dropped from the averaged reduced Hamiltonian .

Error Estimates for the Averaging.

By applying the standard averaging theory for differential equations (see for instance [30]; we again recall that averaging the Hamiltonian is equivalent to averaging the equations), we can conclude that the solution trajectories of the averaged reduced equations are within of the solution trajectories of the reduced equations for at least times.

3.4 Error Estimates of the Entire Treatment.

There are three approximations that we made in our treatment: Taylor approximations of the potential, truncation into two modes, and averaging. We found that the truncation and averaging respectively induces and errors for at least and times. In addition, as long as the solution remains bounded, Taylor approximations result in a small error in the nonlinear forces, which again by Gronwall only induces error in the solution till at least times. Put together, our entire treatment induces at most an error till at least time.

Moreover, this error estimate is valid for an arbitrary forcing function. But for our DNA model, the exponential decay of the Morse potential with the distance implies that when the pendula escapes the immediate vicinity of the opposing pendula, the Morse potential and its consequent perturbation are effectively zero. Hence, our error bound is very loose because it utilizes only the fact that is small, but not the specific feature of our DNA model, which is that the forcing term is also small in a large region of phase space. This is why we numerically observed that our approximation remains accurate over a timescale much longer than .

3.5 Phase Space of Averaged Reduced Equations

Given the averaged reduced Hamiltonian, we obtain the averaged reduced Hamiltonian equations as follow

(19)

Notice that the action of the excited mode is a constant of motion. Hence, the averaged equations of the reactive mode, i.e., the equations, and their phase space structures are parametrized by .

Figure 6: (a) The contour plots for the averaged reduced Hamiltonian when . This phase space has a separatrix which separates two types of motion: liberation (also referred to as “breathing” in [4]) and flipping. (b) Three such phase spaces stacked up in increasing in the space. Notice that the phase space and the separatrix “shrinks” towards the saddle as increases. Together, these separatrices form a homoclinic manifold and can be used in studying the minimum activation energy for each excited mode.

Therefore, the averaged reduced Hamiltonian, Eq. (17), can be used to study the phase spaces of the reactive mode of the averaged reduced equations that are parametrized by . Figure 6(a) shows the contour plots of the averaged reduced Hamiltonian and the phase space for the reactive mode of the averaged reduced equations, when . Notice that this phase space has a separatrix which separates two types of motion: liberation and flipping. Figure 6(b) shows three of such phase spaces stacked up in increasing (action of excited mode). In the space, the vertical axis can also be seen as the axis of increasing activation energy by scaling with , . Notice that the phase space and the separatrix “shrinks” towards the saddle as increases. Together, these separatrices form a homoclinic manifold and can be used in studying the minimum activation energy for each excited mode.

3.6 Analytical Study of Minimum Actuation Energy

Recall that for large enough energy in the excited mode, the energy transfered to the reactive mode will surpass those at the saddle and the system will be driven from its initial equilibrium across the separatrix, causing the flipping. See Figure 4. But for the averaged reduced system, this process manifests itself via the changes in the phase space and the separatrix of the reactive mode parametrized by . See Figure 6(b) for illustration. For activation energy slightly larger than the minimum activation energy, ), the initial stable equilibrium at , now parametrized by , will cross the separatrix parametrized by , move into the flipping region of phase space parametrized by and induce the DNA division. Hence, can be found by the condition that the point is on the separatrix that passes through . Since the separatrix is the energy curve defined by

(20)

where is the averaged reduced Hamiltonian, can be found by solving the following equation for

(21)

where is a saddle at and is the initial equilibrium point at .

As pointed out earlier, of Eq. (17) has been computed using a Taylor expansion of the Morse potential at . Strictly speaking, it may be better to denote it as . Even though we have already tried to find an excellent Taylor expansion that can take care of the approximation of at both the neighborhood near the saddle and the basins around the stable equilibria , the accuracy of and at the latter is still not as good as those at the former. For example, while the actual stable equilibria should be at , the approximated ones will be at if is used. Since any significant error in the approximation of the Morse potential at the stable equilibria will cause a large error in the estimation of the minimum activation energy, another averaged Hamiltonian has also been computed using an expansion of the Morse potential at one of the stable equilibria :

(22)

Here, the letter stands for the word “equilibria”. Now, The minimum activation energy can be analytically estimated by solving the following equation

(23)

Straightforward substitution will give us

(24)

After simplification, the minimum activation energy can be analytically estimated by solving

(25)

For the two-mode truncation of the reactive mode and the 6th mode, . Therefore, , as compared with from numerical simulation.

Figure 7: (a) Analytical estimation of minimum activation energies vs numerical simulation obtained from two mode reduced models. The data for magenta “*”s are from analytical computation. The data for cyan “o”s are from numerical simulation obtained from the two-mode truncations. They match very well. (b) Numerical simulation of two mode reduced models vs numerical simulation obtained from the full model. The blue “o”s are the data from numerical simulation of the full equations (30 modes). The red “*” are from simulation of the two-mode truncations.

Figure 7(a) compares the analytical estimation of minimum activation energy with the numerical simulation obtained from two-mode reduced model. The data for cyan “o”s are from the numerical simulation obtained from the two-mode truncation. For example, the 6th data point (for the 6th mode) with needs minimum activation energy . The data for magenta “*”s are from analytical computation. For example, the 6th data point with has Clearly, the data from numerical simulation and the data from analytical computation match very well. Moreover, (except the 15th mode) can be approximated as a parabola

(26)

Here, we would like to remark that the use of two different Taylor expansions of the Morse potential function , one at and another at , not only allows us to deal with the difficulty posed by the exponential form of the Morse potential function but also improves the accuracy of our analytical estimation of the minimum activation energy for each excited mode.

Figure 7(b) shows the comparison between numerical simulation obtained from the two mode truncations with the numerical simulation obtained from the full model. The blue “o”s are the data from numerical simulation of the full equations (30 modes). For example, the 6th data point (for the 6th mode) with needs minimum activation enesrgy . The red “*” are from simulations of the two-mode truncations . Besides the 14th mode (with error less than ), all other data values of blue “o”s and red “*”s have differences less than 1%. Therefore, for the study of the minimum actuation energy, the analytical estimation provides accurate prediction for the full system. .

3.7 Summary of Analytical Results on Structured Activations

By applying the method of partial averaging, we have obtained the averaged reduced equations Eq. (3.5) for a chain of coupled Morse oscillators. (i) These equations reveal the coupling between the action and energy of the excited mode and the dynamics of the reactive mode, as well as the phase space structure of the activation mechanism. (ii) They allow us to estimate analytically the minimum activation energy for each excited mode and discover a relation between the frequency of the excited mode and their corresponding minimum activation energy. (iii) These estimates match very well with the numerical simulations obtained from the reduced and the full model. The results show that the nearly internal resonance is responsible for the structured activation of our DNA model.

3.8 Remark on Pitchfork Bifurcation

We observe that the reactive mode of the averaged reduced model has a pitchfork bifurcation at . This can be seen clearly by studying the following averaged equations of the reactive mode parametrized by :

(27)
Figure 8: Figure (a) shows that the reactive mode of the averaged reduced system has a pitchfork bifurcation at the parameter value . Figure (b) shows a graph of which has a simple zero at .

Notice that the equilibria of this system are given by and where are the solutions of

(28)

See the blue curve in Figure 8(a). Clearly, are stable equilibria, except at . In fact, they are surrounded by liberation contours as shown in Figure 6. As for the equilibria for each , their stability are determined by the sign of the coefficient of the linear term in , namely, . For , they are saddles. For , they are stable equilibria. Figure 8(b) shows a graph of which has a simple zero at . Hence, the system has a pitchfork bifurcation as shown in Figure 8(a). As will be shown in the following section, this bifurcation will play an important role in our study on the control of this DNA model via parametric resonance.

4 Control via Parametric Resonance

Building on our understanding of its internal dynamics, we want to control the division of this DNA model via parametric excitation, that is in resonance with its internal trigger modes. This effort is guided by two observations and two conjectures: The two observations are (i) the averaged reduced model has a pitchfork bifurcation at —- the physical interpretation of this is that if enough energy is injected into the trigger mode with a value larger than its bifurcation value, the chain of pendula will remain near its open state; (ii) parametric resonance is an efficient way for energy transfer from an external source. Hence, our two conjectures are (i) low intensity electro-magnetic fields can be used with the mechanism of parametric resonance to pump energy into the trigger modes and keep the real DNA chain near the opening state for replication and transcription. Here we note that electric field-induced molecular vibrations have already led to the development of a noninvasive cell transfection protocol that enables foreign DNA molecules to cross cell membranes and penetrate into the cytoplasm by eliciting vigorous vibration between molecules and cells [44, 33]; (ii) enzymes may use similar method to initiate the opening dynamics for DNA replication and transcription. For discussions and evidence of important relations between low-frequency vibrations in DNA molecules (and other bio-macromolecules) and significant biological functions, we refer to [4], [5], [3], [31], [32], [13], [46], [21].

Our three main results in this section are as follows: (i) We identify excitation parameters (as a function of friction) and a class of trajectories that show how parametric resonance can be used to drive the averaged reduced model from its almost equilibrium state to its open state. (ii) By identifying the averaged reduced effective Hamiltonian (after parametric excitation), we uncover the global phase space structure and analyze the characteristic trajectories mentioned above. (iii) We extend the results for the averaged reduced model with parametric excitation and friction to the reduced and to the full model with parametric resonance and friction. The findings uncover a method for controlling DNA division via parametric resonance. Although the issues of inhomogeneity, helicity, and environmental effects (such as noises) are ignored for now, they will be taken into consideration in our future work.

4.1 Equations of Motion with Parametric Excitation

The equations of motion of our full DNA model with parametric excitation and frictions can be written as follows:

(29)

where , , is the Morse potential function and is its derivative. Notice that the left hand side is the original equations of motions, without parametric excitation or frictions. See Eq. 2. Moreover, are the amplitude and the frequency of the parametric excitation, respectively; is in a nearly 1:2 parametric resonance with a chosen internal trigger mode ; is the frictional coefficient.

Note the friction terms represents energy loss caused by interaction with surrounding molecules. Noise terms (ignored here) would represent energy gain caused by (thermal) interaction with surrounding molecules [25]. The parametric excitation terms represent interactions with electro-magnetic fields, we note that DNA vibrations induced by electric-fields or microwave absorbtion are an experimental reality [44, 33].

As pointed out in Ref. [44], the bases of DNA have dipole moments [7], which could couple with an external electromagnetic field. In principle, this coupling can induce a periodic wave-like forcing over DNA bases. Since an electromagnetic field has a polarization, it is reasonable to assume that the coupling between an electromagnetic field and a DNA base is dependent on the orientation of the base, . This motivates us to introduce the parametric excitation term in our model equation (the first term on the right-hand side of Eq. (29) above).

After using the Fourier modal coordinates, the equations of motion for a corresponding two mode reduced model are given by

(30)

where and are the coordinates for the reactive mode and the -mode, respectively; is the reduced Morse potential term and are its partial derivatives with respect to respectively. See Eq. (3.2) and Eq. (10) for detail.

After applying the method of partial averaging to Eq. (4.1) using the angle-action variables

(31)

renaming variables, and setting , we obtain the averaged reduced equation of motion

(32)

where are the action and the phase of the chosen internal trigger mode; is a detuning parameter, defined by ; is the averaged Morse term, and are its partial derivatives with respect to respectively. See Eq. (3.5) for comparison.

The derivation is similar to what has been done in [28, 18] and is parallel to the method used in §A of the Appendix with only minor modification. We first rewrite Eq. (4.1) in the following first order form:

(33)

where

(34)

By using the angle-action variables defined by Eq. (31), we can transform Eq. (4.1) into the Lagrangian Standard Form

(35)

where and

(36)

Hence, we can apply the standard averaging theory (by averaging from to ) and obtain the averaged reduced equations

(37)

where . Then, after carrying the averaging (that is similar to §A), renaming variables, and setting , we obtain Eq. (4.1).

Moreover, the averaged reduced system (without friction) has an effective Hamiltonian

(38)

that can provide insights on the global phase space structure of this averaged reduced model. Notice that since angle-action variables are used in our derivation, the Hamiltonian that we have obtained is canonical. See reference [18] for comparison.

4.2 Merging Local Bifurcation Analysis with Global Geometry of the Effective Hamiltonian

Detailed local bifurcation analysis of the averaged reduced system can be used to reveal the ranges of and other parameters, , where the desired dynamics may be available [26, 28]. For example, for where , the frequency response curves for the averaged reduced system can be depicted as in Figure 10. These frequency response curves can be obtained by studying the fixed points of Eq. (4.1) and their stability that are parametrized by the detuning parameter . Notice first that for all fixed points, . Moreover, there are two types of fixed points, or .

Fixed Points with .

There exist two sub-cases: (i) and the fixed point is , (ii) and the fixed points are . These two types of solutions can be obtained by solving the following equation

(39)

with . Notice that in our study of the fixed points and their stability with , the following Hamiltonian polar coordinates has been used:

(40)

That is why when .

For case (i), the eigenvalues of the fixed point are given by

(41)

where . For example, as decreases from to , this fixed point goes through two bifurcations, one at and another at . In the interval , it is a rank two saddle (saddlesaddle), marked with red dashed lines. In the intervals and , it is a saddlestable foci, marked with magenta dashed and dotted lines. See Figure 10(a)(b).

For case (ii), the eigenvalues of the fixed points are given by

(42)

where is the second partial derivative of the averaged Morse potential term with respect to and is the first partial derivative of with respect to . In the interval , these two fixed points are stable (stable focistable foci), marked with blue lines. See Figure 10(c)(d). Even though they do go through bifurcations at and , we do not include their analysis in this paper because they are far outside of the region where we have found the desired dynamics.

Fixed Points with

For these fixed points, . Again, there exist two sub-cases. (iii) x=0 and the fixed points are where is obtained by solving the following equation for

(43)

with . See Figures 9(a)(b) for the curves of these fixed points. (iv) and the fixed points are where are obtained by solving the two nonlinear equations Eq. (39) and Eq. (43) simultaneously. See Figures 9(c)(d) for the curves of these fixed points.

Figure 9: Figures (a)(b) and Figures (c)(d) show the curves of fixed points for and respectively. Figure (a): (1) the fixed points of the upper left branch of are stable (stable focistable foci) for , marked with the blue line ; (2) the fixed points of lower left branch are saddlesaddle, marked with the red dashed line; (3) the fixed points of right branch are saddlestable foci, marked with the magenta dashed and dotted line. Figure 9(b): for simplicity of presentation, only part of the solution that corresponds to upper left branch in Figure (a) has been marked with the blue line (stable foci stable foci). Figures (c)(d): (1) the fixed points of the left branches are saddlestable foci, marked with the magenta dashed and dotted lines; (2) the fixed points of the right branches are stable (stable focistable foci), marked with the blue lines.

For case (iii), the eigenvalues of the fixed points are given by

(44)

Numerical computation shows that (1) the fixed points of the upper left branch of in Figure 9(a) are stable (stable focistable foci) for ; (2) the fixed points of lower left branch are saddlesaddle; (3) the fixed points of right branch are saddlestable foci. Note: for simplicity of presentation, only part of the solution in Figure 9(b) that corresponds to upper left branch in Figure 9(a) has been marked with blue line (stable foci stable foci).

For case (iv), the characteristic equation of the eigenvalues of the fixed points is given by

(45)

where

Numerical computation shows that (1) the fixed points of the left branches of Figures 9(c)(d) are saddlestable foci; (2) the fixed points of the right branches are stable (stable focistable foci).

Figure 10: Frequency response curves. Figure shows how the fixed points and their stability change as the detuning parameter is varied. Since for all the fixed points, and , only (a) vs and (b) vs are plotted. For example, in the neighborhood of , the system has four types of fixed points for each : (1) the red dashed curves where denote saddle saddle; (2) the solid blue curves where denote stable focistable foci (only positive is drawn); (3) the solid blue curves where denote stable focistable foci; (4) the magenta dashed dotted curves where denote saddle stable foci.

Figure 10 is the combined result of the case studies. The frequency response curves in this figure show how the fixed points and their stability change as the detuning parameter (and hence the frequency ) of parametric excitation is varied. By analyzing the relationship between all these curves, we can see that the averaged reduced model may have the desired dynamics in the neighborhood of . For example, at , the phase space has four types of fixed points as follows:

  1. is a rank two saddle.

  2. where are stable foci. Notice that these fixed points locate at the bottom of the potential well of the averaged reduced system.

  3. where are stable foci. Notice that the coordinates of these fixed points mark the DNA division.

  4. where are stable foci saddle. If the frictional coefficient is small. These fixed points are essentially rank one saddles.

According to the theory of tube dynamics [19, 15, 16], these modified rank one saddles may provide a low energy pathway from the neighborhood of the bottom of the well to the region which marks the DNA division.

Global Geometry of the Effective Hamiltonian

While the local bifurcation analysis does provide many basic ingredients for our study, it does not by itself give a clear and global picture of the dynamics of the averaged reduced system. Hence, the effective Hamiltonian is needed to fill in this gap. Notice that in the Hamiltonian polar coordinates Eq. (40), corresponds to the case where . Therefore,

(46)

provides an effective potential for the averaged reduced system.

Figure 11: (a) Contour plots of the effective potential energy in the space. (b) An example of a trajectory that shows how the parametric resonance drives the averaged reduced system from its almost equilibrium state to its open state in the space.

Figure 11(a) shows the energy contours of this effective potential when . It provides us with the insights for the global phase space structure of the averaged reduced system. From the figure, we can see clearly how the four types of fixed points obtained previously fit together within the global geometry of the effective Hamiltonian. Moreover,

  • The parametric excitation represented by the parameters has turned , which marks the DNA division, into a sink,

  • it also creates two low barrier rank one saddles that are close to the two DNA equilibrium states .

Hence, the addition of parametric excitation to the averaged reduce system should allow certain trajectories with a little energy in the trigger mode to move from an almost equilibrium state, over the saddle, navigate down the energy contours, and reach the sink. All these insights drawn from the local bifurcation analysis and the effective Hamiltonian enable us to generate a class of trajectories that show how the parametric resonance drive the averaged reduced system from its almost equilibrium state to its open state. Figure 11(b) shows an example of this kind of trajectories in the space.

4.3 Parametric Resonance Drives DNA to Division

The data for this trajectory are given as follow. For the system parameter, we have . For the initial condition, we have . The integration is done using the averaged reduced equations, Eq. (4.1). Notice that this trajectory in Figure 11(b) starts at the equilibrium position of the reactive mode but with certain small amount of energy in the trigger mode (th mode). Without the parametric excitation, the system will liberate near the equilibrium state if there is no friction or die down if the friction exists. However, if the parametric excitation is turned on at with the data provided above, the parametric resonance will inject energy into the trigger mode, increase the value of , make the trajectory to reach the region that marks the DNA division () but with large energies in the trigger mode. See Figure 14 below for a physical interpretation of this class of trajectories when it is near the DNA open state.

Here, we would like to make a remark on the amount of initial energy in the trigger mode. Note that if we increase the amplitude of the parametric excitation, the magenta dashed dotted curve of Figure 10(a) will shift downward. Similarly, the rank one saddles in Figure 11(a) will also shift downward. These mean that the amount of initial energy needed in the trigger mode, namely, the value of , can be lower for large . Numerical simulations of the full system confirm this observation.

4.4 Extend the Results to the Reduced and the Full Models

Figure 12 shows two corresponding trajectories, one for the reduced model and another for the full model. The trajectory in Figure 12(a) is generated with the following data. For the system parameters, we set as before. For the initial condition, we set and integrate the trajectory using the reduced equations, Eq. (4.1) where

(47)
Figure 12: (a) An example of a trajectory that shows how the parametric resonance drive the reduced system from its almost equilibrium state to its open state. (b) A corresponding trajectory for the full model.

As for the trajectory in Figure 12(b), it is generated with the following data. For the system parameters, we set as before. For the initial condition, we set and integrate the trajectory using the full equations, Eq. (29) where is obtained by Eq. (47) and are determined by the Fourier modal transformation, Eq. (3).

Figure 13: Figure shows the projection of the same trajectory of the full model on the phase space of its first 15 modes. The parametric excitation injects energy into the 6th mode, some of which transfers to the reactive mode and drive the full system from the almost equilibrium position to the region which marks the DNA division.

Remarks on This Class of Special Trajectories

Here, we would like to make a few remarks:

  • Despited its simplicity (when compared to the full system), the averaged reduced model is surprisingly accurate as illustrated by the fact that the three trajectories in Figure 11(b), Figure 12(a), and Figure 12(b) are very similar. Without a careful study of the averaged reduced equations, it may be difficult to guess that such a class of trajectoris will exist in the averaged reduced model, let alone in the reduced and the full model.

  • It is also interesting to point out that the initial conditions for the reduced and the full models are essentially the same. Hence, if a small amount of initial energy is in a single Fourier mode, the dynamics of the reduced model of a two mode truncation looks very similar to the dynamics of the full equations of our DNA model. Figure