# DNA Breathing Dynamics in the Presence of a Terahertz Field

B. S. Alexandrov  Theoretical Division and Center for Nonlinear Studies, Los
Alamos National Laboratory, Los Alamos, New Mexico 87545  V. Gelev  Harvard
Medical School, Boston, Massachusetts 02215  A. R. Bishop  Theoretical
Division and Center for Nonlinear Studies, Los Alamos National Laboratory, Los
Alamos, New Mexico 87545  A. Usheva  Harvard Medical School, Boston,
Massachusetts 02215  K. Ø. Rasmussen  Theoretical Division and Center for
Nonlinear Studies, Los Alamos National Laboratory, Los Alamos, New Mexico
87545

###### Abstract

We consider the influence of a terahertz field on the breathing dynamics of
double-stranded DNA. We model the spontaneous formation of spatially localized
openings of a damped and driven DNA chain, and find that linear instabilities
lead to dynamic dimerization, while true local strand separations require a
threshold amplitude mechanism. Based on our results we argue that a specific
terahertz radiation exposure may significantly affect the natural dynamics of
DNA, and thereby influence intricate molecular processes involved in gene
expression and DNA replication.

Spectroscopic techniques in the terahertz range are currently emerging as new
tools for the investigation of biological macromolecules Markelz . In spite of
the experimental difficulties (caused by the water’s giant Debye dipole
moment, which leads to a substantial dielectric relaxational loss in the THz
range) a notable amount of research effort is being devoted to the development
of sophisticated THz bio-imaging Nature Photonics . Nevertheless, very little
is known about THz-radiation’s influence on biological systems, and the
mechanisms that govern this influence. The possibility that low frequency
electromagnetic radiation may affect genetic material, enzymatic reactions,
etc. was introduced long ago Frohlich , and since then has been a subject of
constant debate. The energy of such radiation is too low to directly disrupt
any chemical bonds or cause electronic transitions. Only a resonance-type
interaction might lead to an appreciable, biological effect. In biomolecules
such interactions are possible through the ubiquitous hydrogen bonds that have
energies in the THz range. Numerous in vivo and in vitro experiments have been
conducted to clarify low frequency radiation’s ability to cause biological
effects, such as chromosomal aberration, genetic damage etc. The experimental
studies have been conducted under a variety of conditions, but mostly at
frequencies below 0.01 THz, power below 1 $mW/cm^{2}$, and short exposure
times. The data collected in these conditions led to mixed conclusions: some
studies reported significant genetic damages while others, although similar,
showed none Mutation Research . The major international research project,
”THz-bridge” THz-Bridge , which was specifically concerned with THz radiation
genotoxicity concluded that: under some specific conditions of exposure,
change in membrane permeability of liposomes was detected and an induction of
genotoxicity was observed to occur in lymphocytes. Hence, this project
confirmed the existence of THz genotoxicity, but it remains unclear under
which specific conditions such effects occur.

Recent measurements confirm that only extended (6 hours) exposure to a weak
THz field can cause genomic instability in human lymphocytes Korenstein .
Independently, it was reported that neurons exposed in vitro to powerful THz
radiation (over 30 $mW/cm^{2}$ ) cause infringement of the morphology of the
cellular membranes and intracellular structures Russ1 . The same work also
showed that at decreasing power and/or at different frequencies the
morphological changes do not occur. Recently, it was further pointed out that
exposure to a low level THz radiation can interfere with the protein-
recognition processes Korenstein II . It was also shown that exposure of mice
to 3.6 THz high power (15 $mW/cm^{2}$) radiation, for 30 min caused behavioral
changes Russ2 , while under short (5 min) exposure the changes could not be
detected.

Figure 1: Spontaneous formation of a localized permanent opening in dsDNA in
the presence of a THz field with amplitude $A=144$ pN, and $\Omega=2$ and a
spatial fluctuation with high amplitude (see text).

Thus the available experimental data strongly suggest that THz-radiation can
affect biological function, but only under specific conditions, viz. high
power, or/and extended exposure, or/and specific THz frequency.

The present work is devoted to developing a qualitative description of the
existence of such THz effects. We argue that the appropriate conditions
correspond to irradiation parameters that assure the existence of linear and
nonlinear resonances between DNA conformational dynamics and THz field Fig.1.

We will base our analysis on the Peyrard-Bishop-Dauxois (PBD) model of dsDNA 2
, which is arguably the most successful available model for describing the
local DNA pairing/unpairing (breathing) dynamics 6 . Of particular importance
for the THz effect is our previous demonstration of strong correlations
between regulatory activity, such as protein-DNA binding and transcription,
and the equilibrium propensity of dsDNA for local strand separation -8 ; 3 ; 4
; 5 . With this background it is natural to consider the influence of THz
radiation on dsDNA dynamics within the PBD modeling framework. One
complication is that the specific physical nature of the interactions between
DNA and the THz electromagnetic field is not known in detail. However, given
the sensitivity of THz radiation to the strand pairing state of DNA 15 , we
will here simply augment the PBD to include a drive in the THz frequency
range, without specifying the precise nature of the underlying physical
coupling. Since the wavelength of the THz field is larger than 100 $\mu m$ and
the characteristic size the investigated dsDNA sequences is less than 20 $nm$
it is reasonable to consider an uniform THz field. In this fashion the
breathing dynamics (represented by the normalized separation $y_{n}$ of
complementary bases) of the $n$’th base pair in the presence of a
monochromatic spatially homogeneous external field is described by

$\displaystyle m\ddot{y}_{n}=$ $\displaystyle-$ $\displaystyle
U^{\prime}(y_{n})-W^{\prime}(y_{n+1},y_{n})-W^{\prime}(y_{n},y_{n-1})$ (1)
$\displaystyle-$ $\displaystyle m\gamma\dot{y}_{n}+A\cos{(\Omega t)},$

where, the Morse potential $U(y_{n})=D_{n}(\exp{(-a_{n}y_{n})}-1)^{2}$
represents the hydrogen bonding of the complementary bases. The parameters
$D_{n}$ and $a_{n}$ depend on the type of the base pair (A-T or G-C).
Similarly,
$W(y_{n},y_{n-1})=\frac{1}{2}\chi[y_{n},y_{n-1}](y_{n}-y_{n-1})^{2}$ with
$\chi[y_{n},y_{n-1}]=k(1+\rho\exp{(-\beta(y_{n}-y_{n-1}))})$ represents the
stacking energy between consecutive base pairs. The term $m\gamma\dot{y}_{n}$
is the drag caused by the solvent while $A\cos{(\Omega t)}$ is the (THz)
drive. For simplicity, we consider here a homogeneous poly(A)
note_on_parameters  DNA molecule with 64 base pairs. Our numerical simulations
of this set of nonlinear coupled equations (using periodic boundary
conditions) showed that the primary response of the system is spatially
homogeneous (Fig.1 the first 40 $ps$). However, we observed that a presence of
fluctuations may lead to the creation of persistent spatially localized
openings of the double-stranded molecule, Fig.1. Because, such openings
(bubbles) are known to functionally affect dsDNA we focus here on
characterizing the conditions under which these states appear.

Since, the primary response of of the driven system is spatially uniform it
can be understood in terms of a single classical damped and driven Morse
oscillator. Such systems have been studied quite extensively in many other
contexts 7 , and are known to display rich dynamical behavior including limit
cycles, period-doubling bifurcations etc. An illustration of this is shown in
Fig.2, A, for our

Figure 2: A) Bifurcation diagram for a poly(A) molecule in an external THz
field with amplitude $A=144$ pN. B) The instability curves for wave vector
$k=0$ (solid curve) and $k=\pi$ (dashed curve) are plotted. On the vertical
axis is the amplitude of the external field, and on the horizontal axis is the
frequency $\Omega$ of the THz field.

system parameters, in terms of the Poincaré section $\dot{y}_{n}(t=mT)$, where
$T=2\pi/\Omega$, versus the driving frequency $\Omega$ for several integer
values of $m$. Feigenbaum’s period-doubling route to chaos is clearly seen in
this figure, which demonstrates that this system’s main instability mechanism
is through period-doubling. As described earlier, our interest is to
investigate whether the presences of a THz field can lead to the creation of
spatially localized unbinding of the DNA double-strand. Such states are most
easily created in connection with the period-doubling events, as a result of a
spatially inhomogeneous perturbation. In order to investigate this phenomenon
we adopt a method developed in Ref. 9 . In this work the linear stability
analysis (based on the Floquet theorem) of a single Morse oscillator was
performed, and it was shown that in a rotating-wave approximation Bessel  the
dynamics ensuing from a given perturbation is governed by the extended Hill
equation. The stability analysis of the Hill equation has been developed by
Ince 11  based on his extension of the Whittaker’s method of solving Mathieu’s
equation 20 . Since, our system is spatially extended, a perturbation may be
spatially inhomogeneous. In this case the response of the system is
characterized by $y_{n}(t)=y^{0}(t)+z_{n}(t)$, where $y^{0}(t)$ is the
spatially homogeneous and temporally periodic solution described above (see
Fig. 2), A, and $z_{n}(t)$ represents the perturbations in whose dynamics we
are interested. In the most general form the perturbation is given in terms of
Fourier modes: $z_{n}(t)=\sum_{k}\exp[ikn]\xi_{k}(t)$ with wavenumber, $k$,
and amplitude $\xi_{k}(t)$. As with the single oscillator, the dynamics of the
Fourier amplitudes $\xi_{k}(t)$ of the perturbation are governed by set of
uncoupled extended Hill equations. However, the coefficients of these
equations depend on the wave-number $k$. The Ince stability analysis allows us
to analytically determine note1  the Floquet exponents (i.e. the stability
criteria) as functions of the spatial profile as given by $k$. For more
technical details on the method see Ref. 12 . The result of our analysis are
the stability curves dividing the parameter space $(A/m,\Omega,k)$ into
stability and instability zones. Figure 2, panel B shows the stability
boundaries of a spatially homogeneous perturbation ($k=0$, solid curve) and of
a staggered ($k=\pi$, dashed curve) perturbation. Perturbations represented by
other wave-numbers have stability boundaries in-between those shown in Fig. 2,
B. Below these curves the spatially homogeneous states perform stable temporal
oscillations at the frequency of the external drive and its higher harmonics,
while above the curve the perturbation leads to period-doubled oscillations.
It is clearly seen that for driving frequencies below $\Omega\sim 2.5$ THz
spatially homogeneous ($k=0$) perturbations lead to the strongest instability
(occurring at the lowest value of the driving amplitude). However, above
$\Omega\sim 2.5$ THz the staggered ($k=\pi$) perturbation becomes the primary
instability mode. In the region above the stability boundary the period-
doubled state will therefore tend to adopt the spatial profile of the most
unstable Fourier-mode (homogeneous or staggered). The presence and location of
the

Figure 3: The $k=\pi$ spatial pattern, at $\Omega=2.52$ THz, is presented. The
upper panel shows evolution of the system. The lower panel shows the power
spectrum of the breather motion. Figure 4: Breather formation at $\Omega=2.00$
THz. The upper panel shows the evolution of the breather. The lower panel
shows the power spectrum of the breather motion.

obtained instability curves (Fig. 2, B) were verified by exploring the
parameter-space $(A/m,\Omega,k)$ by direct simulations of Eq. (1). The
dynamics of small spatially random (uniformly distributed in the interval
$[-\epsilon,\epsilon]$ and thus containing all wave-numbers) perturbation
$z_{n}$ introduced on top of the spatially homogeneous and temporally periodic
solution $y^{0}(t)$ was followed numerically. As this perturbation injects
energy into all wave-numbers, the ensuing dynamics will emphasize the mode
with the strongest instability. For point A (Fig. 2, B), the results are shown
in Figure 3. The upper panel shows that after the introduction of the
perturbation (at 40 ps) a spatially staggered structure forms. The power
spectrum of this dynamics (lower panel in Fig. 3), shows that the frequency of
the vibrations is mainly $\Omega/2$, i.e. a period-doubling transition has
occurred as expected from the above analysis. The power spectrum (lower panel)
shows that the signature of the drive persists at $\Omega$, in addition to a
component at zero frequency, which results from the asymmetry of the base-pair
potential $U(y_{n})$. A similar scenario is observed for driving frequencies
below $\Omega\sim 2.5$ THz, but in these cases the solutions remain spatially
homogeneous after the introduction of the perturbation. The fact that the
instability curves for $0<k<\pi$ lie between the two curves in Fig. 2, B shows
that it is, in general, impossible for any spatially localized feature with a
finite ($2~{}\mbox{base pairs}<\lambda<\infty$) length scale $\lambda$ to
emerge as a result of a linear instability. This is in contrast with many
other nonlinear condensed matter models, where spatially localized and
temporally periodic modes (intrinsic localized modes or breathers) commonly
arise through linear modulational instabilities. For our model of dsDNA a
perturbation of finite amplitude $\epsilon$ is required to create localized
unbindings of the DNA double strand through nonlinear mechanisms. An
illustration of this is given in Fig. 4 for $\Omega=2.0$ THz, where the
dynamics caused by a random perturbation with $\epsilon>0.14$ is shown to
result in a localized opening of the DNA double strand. Again the power
spectrum, given in the lower panel of Fig. 4 shows, as expected, that the
localized state vibrates at the frequency $\Omega/2$ in contrast to the
uniform background, which remains synchronized to the external drive at
$\Omega$. A more detailed study of the nonlinear mechanisms underlying the
creation of such localized unbinding states reveals that a localized injection
of a certain amount of energy is required. Due to the asymmetry of the base
pair potential, an effective way to achieve this is a localized compression of
the double strand. Using this insight we found the perturbation
$z_{n}=0.42\cos{(\frac{\pi}{4}(n-n_{0}))}$, for $-4\leq n-n_{o}\leq 4$
($z_{n}=0$ otherwise) to be an effective perturbation for the creation of a
localized unbinding state at $n_{0}$, Fig.1. The phase portrait (not shown) of
the base pair $n_{0}$, at the transient time period just before the breather
is established, reveals the breather formation mechanism. First, the locally
injected energy from the compression must be sufficient for a few consecutive
base pairs to undergo a temporary local melting transition, i.e. reach the
plateau of the Morse potential. Subsequently, the stacking interactions and
friction dissipate the energy to gradually reach a steady state. At given
driving amplitude and frequency, the dsDNA molecule has two steady
conformational states with different average energies and amplitudes. Only in
these states is the externally pumped energy exactly compensated by the energy
dissipated by the friction. The first state is vibrating predominantly at the
driving frequency $\Omega$, while the second mainly oscillates at $\Omega/2$.
Hence, the breather formation proceeds in three stages: 1) all the base pairs
are in the first state then; 2) a few of them undergo a temporary melting
transition as a result of the perturbational energy injection. Finally, 3) the
perturbed base pairs reach the second steady conformational state, after
dissipating their additional energy.

In summary, we have found that a THz field may cause dynamical separations of
the DNA double-strand. In the presence of weak perturbations e.g. thermal
fluctuations, small amplitude response occurs (at half the driving frequency)
either in a spatially uniform manner or, at higher-frequencies, in an unusual
spatially dimerized form. In the latter case, neighboring base pairs oscillate
in an out-of-phase fashion. However, large localized openings (bubbles) in the
DNA double strand can only occur via a nonlinear mechanism requiring a spatial
perturbation above a certain amplitude threshold that is determined by the
intensity and the frequency of the THz field.

We previously showed experimentally that the introduction of an artificial
permanent bubble in dsDNA, via a short mismatched segment, is in fact
sufficient for transcription even in the absence of any auxiliary factors 5 .
Our present finding of large localized openings resulting from an external THz
field therefore underscores the importance of including the interactions of
genomic DNA with the surrounding environment in DNA models. The amplitudes of
the resulting openings observed in the presence of a THz field are significant
compared to the dynamic signatures of protein binding sites and transcription
start sites 4 . This suggests that THz radiation may significantly interfere
with the naturally occurring local strand separation dynamics of double-
stranded DNA, and consequently, with DNA function.

Based on the model results present here, we believe that the main effect of
THz radiation is to resonantly influence the dynamical stability of the dsDNA
system. Hence, our instability curves define the parameters space (i.e. the
amplitude of the field (or the power), and frequency) in which THz radiation
can have an immediate effect. In contrast, nonlinear instability may occur at
any point of the instability diagram but it requires fluctuations with
significant amplitudes. In biological systems, such fluctuations are generated
thermally. Hence, the occurrence of a fluctuation with sufficiently large
amplitude is very rare, and therefore extended exposure is required for the
THz effect to take place, via nonlinear instabilities, especially if the power
is small. In this framework, it is than natural that the character of THz
genotoxic effects are probabilistic rather than deterministic.

We acknowledge Dr. Voulgarakis for initial discussions regarding the subject
of this work.This research was carried out under the auspices of the U.S.
Department of Energy at Los Alamos National Laboratory under Contract No. DE-
AC52-06NA25396. and it was supported by the National Institutes of Health (R01
GM073911 to A.U.).

## References

  * (1) A. G Markelz, IEEE Journal of Selected Topics in Quantum Electronics, 14, 180 (2008).
  * (2) L. Ho, M. Pepper, P. Taday, Nature Photonics, $\bf{2}$, 541 (2008). A. P. George et al. Opt. Exp. $\bf{16}$, 1578 (2008).
  * (3) H. Frohlich, Proc. Nat. Acad. Set. USA, $\bf{72}$, 4211 (1975).
  * (4) L. Verschaeve, Mutation Research, $\bf{681}$, 259 (2009).
  * (5) http://www.frascati.enea.it/THz-BRIDGE/, (2004).
  * (6) A. Korenstein-Ilan et al. Rad. Res. $\bf{170}$, 224 (2008).
  * (7) J. S. Olshevskaya et al., Proc. IEEE Sibircon, 210 (2008);J. S. Olshevskaya et al., Zhurn. Vyss. Nervn. Deyat. Im. P. Pavlova, $\bf{59}$(3), 353 (2009).
  * (8) A. Homenko et al.,Bioelectromagnetics $\bf{30}$,167 (2009).
  * (9) N. P. Bondar et al., Bull. of Exp. Biol. and Med., $\bf{145}$, 4 (2008).
  * (10) M. Peyrard, and A.R. Bishop, Phys. Rev. Lett. $\bf{62}$, 2755, (1989);T. Dauxois, M. Peyrard, A.R. Bishop, Phys. Rev. E $\bf{47}$, R44 (1993).
  * (11) B. S. Alexandrov et al., J. Phys. Cond. Matter $\bf{21}$, 034107 (2009), and references therein.
  * (12) B. S. Alexandrov et al., Nucleic Acids Res. 37, 2405 (2009).
  * (13) C.H. Chu et al.,Nucleic Acids Res. $\bf{32}$, 1584 (2004).
  * (14) C.H. Chu et al., B. Phys. J. $\bf{95}$, 597 (2008).
  * (15) B. S. Alexandrov et al., PLoS Comp. Biol. $\bf{5(3)}$, e1000313 (2009); A. Usheva, T. Shenk, Proc Natl Acad Sci US $\bf{93}$, 13571, (1996).
  * (16) M. Brucherseifer et al., Appl. Phys. Lett. $\bf{77}$, 11 (2000).
  * (17) The specific parameter values are: $D_{n}$ =0.05 eV, $a_{n}$ = 4.2 $\AA^{-1}$, k=0.025 eV$\AA^{-2}$, $\beta$=0.35 $\AA^{-1}$, and $\rho$=2. The linear resonance frequency is $\omega=$1.22 THz. For the solvent drag we use $\gamma=1$ps-1, consistent with a typical relaxation time of water.
  * (18) C. Lie George and Yuan Jian-Min J. Chem. Phys. $\bf{84}$ (10) (1986).
  * (19) A.M. Samson et al., Radiofizika $\bf{33}$, 49 (1990).
  * (20) In this specific rotating-wave approximation the Morse potential is expanded using the equality: $\exp{(r\cos(\phi))}=I_{0}(r)+2\sum^{\infty}_{n=1}{I_{n}(r)\cos(n\phi)}$ where $I_{n}(r)$ are modified Bessel functions of the first kind.
  * (21) E. L. Ince, Roy. Astronom. Soc. $\bf{75}$, 436 (1915).
  * (22) E. T. Whittaker, Proc. Edin. Math. Soc. $\bf{xxxii}$, 76 (1913).
  * (23) For the analytical work we used a harmonic stacking potential, $\chi=const$.
  * (24) A. Vanossi et al., Phys. Rev. E $\bf{62}$, 7353 (2000).