Continuum Rod Modeling of DNA
In order to apply a continuum rod model to DNA, one must
extract appropriate continuum parameters (e.g. intrinsic shape, stiffnesses)
from existing experimental or computational data:
X-ray crystal structures
molecular mechanics energy-minimization computations
molecular dynamics computations
minimum-energy shapes within the
wedge-angle model
In all of these cases, the model underlying the data is not a continuum rod, but
rather some discrete model, either an all-atom model (for the first three) or
a base-pair-level model (for the last). The challenge is to design
algorithms to translate information from these discrete
models to the continuum model, where we can take advantage of
rod theory and efficient computations.
The next section gives an example of this parameter extraction, in the
case where the underlying model is the wedge-angle model
The wedge-angle model
The wedge-angle model describes each DNA base-pair as a rigid body,
or equivalently a 3D axis system or frame
(d1,d2,d3).
The origin of the frame is at the center of the base-pair, the d3 axis points to the
next base-pair center, and the d1 axis tracks the rotation of the DNA sugar-phosphate
double helix.
DNA intrinsic shape is incorporated into the wedge-angle model by allowing
the minimum-energy stacking orientation of one frame on the previous
frame to depend on the base-pairs involved. These minimum-energy
stacking orientations have been computed by various researchers,
using a mixture of real experiments and atomic simulations (although
there is still considerable disagreement among these results).
Smoothing the wedge-angle intrinsic shape
Given a wedge-angle minimum-energy shape, we convert it to a
continuous shape by a smoothing procedure
involving a combination of data filtering and least-squares fitting.
Many of the mathematical techniques in this smoothing are standard, but the
particular combination required to produce a continuum shape
with the same global curvatures as the wedge-angle shape but
smoother local curvatures is specialized to this DNA problem.
The natural frame
In addition, there is an important transformation required to
accommodate the rapid intrinsic twist of the DNA double-helix
(which makes a full turn every 10.5 base-pairs); if not removed,
this intrinsic twist causes undesirable rapid variations
in the continuum rod parameters. Instead of using
the real DNA d1 vector field of the intrinsic shape, we replace it by
the natural, or parallel-transport, or zero-twist, frame D1.
This natural frame D1(s) is uniquely defined (up to choice of
D1(0)) by the centerline.
For each value of the arclength parameter s, let w(s) denote the
angle between d1(s) and D1(s) in the
rod intrinsic shape.
The remarkable fact is that if the continuum rod is uniform
and has circular cross-section (i.e.
K1(s)=K1=K2(s)
in the energy functional),
then a static equilibrium
of the rod with the real frame can be recovered exactly
from a static equilibrium of the rod with the natural frame
simply by rotating d1(s) (and d2(s)) by w(s).
Said another way, the real-to-natural-frame transformation and
the continuum equilibrium computation commute, so that
we may adopt the following strategy:
(1) Transform from real frame to natural frame
(2) Compute equilibrium configuration
(3) Transform from natural frame to real frame
Real DNA, by virtue of its anisotropic base-pairs, probably does not
have equal bending stiffnesses K1=K2 as assumed
in the above argument. However, since DNA has large intrinsic twist,
this anisotropy should be effectively averaged out over the length scale
of several base-pairs, so that the above assumption is accurate (see
Kehrbaum & Maddocks (1998)
for proofs of results of this type).
Verification
The accuracy of the wedge-angle-to-continuum procedure has been verified
as follows.
Given a DNA sequence, construct its wedge-angle intrinsic shape.
Apply the above procedure to produce a continuum intrinsic shape.
Compute the
lowest-energy ring equilibrium
with d1(0)=d1(1)
within the continuum model for that intrinsic shape.
Take this continuum solution as an initial guess
for the determination of the lowest-energy ring equilbrium
with d1(0)=d1(1) within the wedge-angle model
(this requires a high-dimensional
constrained minimization not described here).
By this procedure,
we directly measure the error introduced in our
conversion to the continuum problem, which is quite small, both
in energy (less than 0.5%) and in configurations (see the
superpositions of centerlines and helices below):
Modelling DNA cyclization
Cyclization as an experimental probe of DNA mechanical properties
Many biophysicists are trying to understand the role of
DNA intrinsic shape and flexibility in phenomena
such as transcription. However, the direct measurement of
these shapes and stiffnesses is often difficult, so instead
we turn to indirect observations such as cyclization
(the joining of two ends of a DNA molecule to form a looped molecule).
Cyclization is a promising probe for DNA shape and flexibility
because the energy of a cyclized molecule (and hence its probability of
forming) depends quite sensitively on its intrinsic shape (e.g. C-shaped
molecules form cycles much more readily than S-shaped ones)
and stiffnesses (e.g. a very flexible sequence in just the right location
can lead to a very low cyclization energy).
We thus are faced with an inverse problem:
given experimental cyclization energies, estimate parameters
of shape and flexibility of the DNA by coupling the experiments to a
model such as the continuum rod model.
Test problem: continuum predictions of cyclization probabilities
(cf. Manning, Maddocks, and Kahn (1996) )
As a first test, we conducted a study
with Jason Kahn
(Dept. of Chemistry and Biochemistry, University of Maryland College Park)
to solve the forward problem:
given currently accepted intrinsic-shape and stiffness parameters
for "normal" DNA (i.e. DNA containing no sequence thought to have properties
in disagreement with the parameters),
can the continuum computations recreate experimental cyclization
probabilities? The test problem involved
11 molecules, of lengths between 150 and 160 base-pairs, each
containing an A-tract (intrinsically bent by about
90o) as well as a CAP binding site (intrinsically bent
by about 10o). Within the family of 11 molecules, the
length of the region between the two bends is varied so that the
spectrum between "C-shaped" and "S-shaped" molecules is covered.
Within the continuum model, a cyclized molecule is
exactly a twisted ring equilibrium:
As described elsewhere, we compute the set of stable
twisted ring equilibria as 
is varied. The solution branches for three of the 11 molecules (branch color
indicates the molecule) are superimposed in the figure below.
The energy versus m3 of each ring equilibrium
is plotted. In addition, the actual cyclized equilibria
(those with = 0)
are indicated with circles. Stable equilibria are indicated with
solid lines, and unstable equilibria with dashed lines.
For each molecule, the lowest-energy stable cyclized equilibrium
should be the one seen experimentally. (If two stable cyclized
equilibria are close enough in energy, they should both be seen
experimentally. Further, in this case, the two equilibria are
likely to have links differing by 1. Both of these predictions are
confirmed by experimental results).
The energies of these lowest-energy stable cyclized equilibria
are compared to the experimental results for the 11 DNA molecules
in the figure below. The continuum energies are plotted
as a function of the ratio of twisting stiffness to bending stiffness
K3/K1, since the appropriate value
of this parameter for DNA is under some debate. The table at the right
shows the experimental results.
The continuum model does not contain entropy effects,
so its energies are a constant shift 
down from the experimental
free energies. The ordering and relative spacing of the energies
match experiment to within experimental error, and if we include
the entropic shift as a free parameter, the continuum results match
experiment to within experimental error.
Computational cost
The entire process of starting from a base-pair sequence, constructing
a continuum intrinsic shape, and computing from it cyclization equilibria,
takes about 5 minutes for the above molecules. Further, once this process
is completed, it takes only a few seconds to compute the
cyclization equilibria at a different value of a parameter like
K3/K1. Since the central computation in the
interpretation of cyclization experiments is a search among DNA shape
and flexibility parameters to best
match experimental cyclization energies, the appeal of these continuum
computations is clear.
Conclusions from test study
This test case thus demonstrates that apart from an entropic term,
the continuum computations can solve the forward problem of
predicting experimental cyclization energies given "correct"
rod parameters. Even without the ability to compute this
entropy shift, the continuum model offers
a promising approach to solving the inverse problem,
using the following strategy. Given experimental cyclization
energies, use the continuum model to search within a space
of DNA shape and stiffness parameters plus an entropic
shift parameter to find the best fit to the experimental
results. Then, as a verification, perform a single
Monte Carlo simulation (or some other method that does account
for entropy). Since the continuum computations
are so much faster than Monte Carlo or similar computations (seconds
per parameter step versus tens of hours), this strategy should be
quite efficient.
Shape and Flexibility of the TATA box
(with Jason Kahn)
Having established the value of the continuum approach for the
forward problem of predicting cyclization energies, we come
to the real goal: to use the continuum approach in inverse
problems to extract mechanical properties of special DNA
regions of biological interest.
We are currently pursuing
this question, in collaboration with Jason Kahn, to study
the TATA-box sequence (a key player in the initiation of transcription),
both with and without the presence of the TATA-box-binding-protein.
In each case, cyclization
experiments are well different from what would be predicted
from "standard" DNA shape and stiffness parameters, so we are
seeking to use the continuum model to determine a shape and/or
stiffness profile of the TATA-box region consistent with the
experimental results.
We have studied the relationship
between intrinsic DNA shape and multiplicities of stable
nicked equilibria and stable cyclized equilibria.
Nicked DNA are DNA loops
in which only one of the sugar-phosphate strands is closed.
We adopted the model for nicked DNA of Katritch and
Vologodskii (1997) in which r(0)=r(1) and d3(0)=d3(1) for
the rod, but no constraint is placed
on the relative orientation of d1(0) and d1(1).
Equivalently, we seek points in
E vs. m3 cyclization bifurcation diagrams
which are both stable and local minima in E.
It is clear that real nicked DNA can not spin freely about its
tangent vector as this model suggests, but we have adopted it
for comparisons with the results of Katritch and Vologodskii.
When no intrinsic curvature is present, then because of the
symmetries of the rod, there
is a degenerate circle of (neutrally) stable nicked equilibria
(corresponding to the bottom-most point in the
perfect diagram), as well as
degenerate circles of (neutrally) stable cyclized equilibria (which precise
links are stable depends on the length of the DNA and the value
of K3/K1).
Perturbation results
show that for infinitesimal perturbations,
these circles of degenerate stable equilibria are expected to perturb
to a pair of equilibria, one stable and one unstable, unless
two special integrals
involving the rod intrinsic shape vanish.
Thus, for the macroscopic intrinsic curvatures seen in real DNA,
the likelihood of multiple stable equilibria should correlate well
with the smallness of these two special integrals. This hypothesis
has been verified via tests on large sets of DNA sequences.
Return to my home page
rmanning@haverford.edu
