Effects of Protein-Induced Local Bending and Sequence Dependence

Article pubs.acs.org/JCTC

Cite This: J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Effects of Protein-Induced Local Bending and Sequence Dependence on the Configurations of Supercoiled DNA Minicircles Yoav Y. Biton* Department of Mechanical Engineering, SCE, Shamoon College of Engineering, Beer Sheva 84100, Israel ABSTRACT: The statistical mechanics of a short circularized DNA molecule, a DNA minicircle, with a prescribed linking number depends heavily on the mechanical and geometrical properties of the DNA, which are known to be functions of the sequence. A description of a general numerical scheme used for the performed advanced simulation is presented with examples that reveal the effects of sequence dependence and local deformations, caused by protein binding, on the configurations in a canonical ensemble of a supercoiled minicircle. Using a realistic course grain model in which the sequence-dependent elasticity, the intramolecular electrostatic interactions, and the impenetrability of the DNA molecule are taken into account, the bifurcation of equilibria of supercoiled minicircle DNA with the induced bending as a parameter is shown. The unique Monte Carlo scheme for constrained DNA molecules was utilized in order to calculate site-to-site distance probabilities for each pair of sites in the molecule. The simulated examples show not only that the sequence alone can play a significant role in bringing two remote sites to a contact but also the possible existence of several competing regions of contact. The results suggest that the local deformation caused by protein binding can yield a global configurational change, dominated by slithering motion, which brings two (originally) remote sites to close proximity, and that the nature of such effect is related to the sequence architecture.

1. INTRODUCTION The sequence composition comprising a DNA molecule bears more than just the genomic information. A significant part of it is composed of interspersed segments that have no genetic functional role but serve as spacers between genes and have important structural functionality. Whether or not it comprises a gene, the sequence-dependence of both the intrinsic, i.e., stress-free, configuration of a DNA molecule (or a subsegment of it) and its local elastic resistance to modes of bending and twisting plays a crucial role in a diversity of processes including DNA packaging,1 genetic recombination,2 DNA looping,3,4 and kinetoplast DNA replication.5 As early as about three decades ago it was suggested that a DNA segment that is wrapped around a single histone octamer reaches a configuration of approximately minimal elastic energy by orienting its highly curved subsegments in a way that tends to minimize the amount of excess bending due to the confinement of the DNA to the cylindrical geometry of the histone octamer.1 DNA molecules of the same length but different sequence composition that are subject to identical external forces and geometric constraints may have very different mechanical responses, because their mechanical properties and intrinsic geometry depend on the sequence. These sequence-dependent properties of the DNA, i.e., the local elastic moduli and intrinsic geometry, are known to have significant departure from the mean values that are normally used when the molecule is treated as homogeneous.6 [See also the discussion in ref 7.] Thus, the mechanical response of a DNA molecule cannot be described accurately enough using models that do not take into account the sequence-dependent mechanical properties of the © XXXX American Chemical Society

DNA. On the other hand, the very important and nontrivial models, based on the regularity of the DNA structure, that treated the DNA as a continuous, homogeneous, and transversely isotropic elastic rod with values of elastic bending and twisting moduli appropriate to the mean observed persistence length of double stranded B-DNA8 can provide the general shape of equilibria of long enough DNA molecules. Other models treating the DNA as an inhomogeneous, continuous rod were presented in refs 9 and 10. [As was mentioned in ref 7, such treatments are based on the special Cosserat theory of rods, which was described for example in ref 11 and was generally restricted to inextensible and unshearable rods with no coupling between the modes of deformations.] When the impenetrability of a DNA molecule is accounted for [by proper treatment of contact points and curves between subsegment of the same molecule] and the molecule is assumed to be intrinsically straight and homogeneous, it is possible to analytically calculate equilibria of highly supercoiled molecules as functions of the topological state of the closed molecule.12−14 The works of Coleman and Swigon showed the existence of multiple equilibrium configurations for given end conditions and topological state.14 The use of these elastic rod models made it possible to study the minimum energy configurations of circularized or looped DNA molecules. However, when getting into the resolution of base pairs, it is necessary to distinguish between closed configurations that differ only by modes of eversion, i.e., Received: October 28, 2017 Published: March 20, 2018 A

DOI: 10.1021/acs.jctc.7b01090 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

intramolecular electrostatic forces exerted by its phosphate groups, and because these forces are not completely screened out by positively charged counterions, they can render the molecule’s equilibrium configurations sensitive to changes in the ionic strength of the aqueous medium. In research on DNA elasticity it has been common practice to assume that the effects of changes in the ionic strength can be accounted for by doing little more than allowing the persistence length and effective diameter of the DNA to be functions of the ionic strength, and then estimating the values of these functions under the further assumption that the DNA molecule is approximately straight. These assumptions, which are not correct for intrinsically curved DNA or even for intrinsically straight DNA that has been circularized into minicircles or supercoiled plasmids, were not made in the theory and calculations presented in refs 17−19. In ref 18 the molecule was treated as a tube-like structure with boundaries composed of spheres (with centers in coincidence with the assigned centers of the base pairs) and cylinders (with axes connecting adjacent barycenters), so that contact points were accounted for accurately using exact inequality constraints. This geometric approach was employed here for the excluded volume treatment in the Monte Carlo calculations. In order to be able to account for thermal fluctuations we have very recently developed a Monte Carlo scheme for direct generation of configurations confined to any prescribed end conditions and topological state, in which intramolecular electrostatic interaction, excluded volume effects, sequence dependence of intrinsic geometry, and elastic moduli are all taken into account,20 and used it to predict Lac Repressor mediated DNA looping probabilities. Here, for the first time, a thorough calculation of equilibrium configurations is combined with a Monte Carlo approach in order to investigate the effects of both the sequence composition and the protein-induced bending on the configurations of a short supercoiled minicircle. In general, when a protein binds to a DNA site of order 10− 20 base pairs in length, it renders a significant local deformation and stiffens the fully bound site. [See, e.g., the review of Harteis and Schneider.21] Such deformation results from each of the dimeric arms of the Lac repressor tetramer,22−24 by Heat Unstable (HU) proteins,25 or by the leucine zipper domain in regulatory proteins, which tend to create a local bend that is oriented toward the minor groove.26 Recent molecular simulations showed how proteins that bind to nonspecific DNA sites can affect the configuration and elastic response of DNA.27 The simulated effect of local bends, caused by the nonspecific binding of HU proteins to the DNA, on the Lac Repressor mediated looping propensities was reported in ref 28. I here present simulations concerning two possible effects on the favorable configurations of supercoiled DNA. The first is the effect of sequence-dependent intrinsic geometry and elastic moduli on the statistical distribution of configurations of supercoiled minicircles. In my virtual lab I chose, as an example, to investigate a 339 bp minicircle DNA with linking number 29, i.e., a molecule that was twisted about three complete turns prior to its cyclization. The production of this minicircle, with its specific sequence and topology, was performed successfully and reported in ref 29. By analyzing the complete statistics of site-to-site distances, I show how a given sequence can result in two competing states with different pairs of contact regions. The second effect is that of an induced local bending on the favorable configurations of supercoiled 339 bp minicircles that

modes varying the rotational orientation of the DNA molecule [Modes varying the rotational orientation were discussed in ref 15, where it was experimentally suggested that CRP induced local bending can obviate eversion modes of a short, (nearly) intrinsically straight, 284 base-pair (bp) minicircle DNA. I refer to modes varying the rotational orientation as eversion modes. Modes of eversion (or prevertion) preserve the axial curve connecting the assigned centers of the base pairs along the molecule and the magnitude of total bending and twisting between adjacent base pairs but change the two orthogonal components of the local curvature.] or slithering [A slithering mode of a closed molecule is referred to here as a mode characterizing motion that shifts the base pairs along the unchanged axial curve, like material points in a synchronous rubber belt that transmits rotational motion between a couple of wheels.], which for the case of a transversely isotropic and homogeneous molecule are indistinguishable, because they do not cause any change in the elastic energy. As a result the molecule may have infinitely many equilibria at that topological state. Examples of these two modes are shown in Figure 1 for

Figure 1. Schematic description of slithering (left) and eversion (right) modes in a supercoiled DNA minicircle shown as a closed tube. The left column shows how the motion shifts the blue subsegment from right (shown in the bottom of the column) to left (top). The right column shows how the rotational orientation of each cross section is varied. The material points marked with a red line along the molecule change their orientation with respect to the (unchanged) axial curve of the molecule from outward (bottom) to inward (top). To make this motion visible additional arrows are depicted. Each such arrow is in the cross section associated with its base pair. It originates at the center of its cross section and points in the direction of the material point marked in red on the external surface of the tube.

homogeneous and transversely isotropic mincircle DNA. It is clear from the figure, that for a specific inhomogeneous sequence each configuration would generally have different elastic energy. The evidence regarding the significance of the inhomogeneity of the DNA motivated the development of a naturally discrete models for the elastic response of DNA that can account for the sequence dependent intrinsic geometry and elastic moduli for any given function relating translational and angular deformations to the local elastic energy of the molecule.16 Furthermore, since a DNA molecule is subject to B


Article


(n+1)-th base pair with respect to the n-th, i.e., by a function of the components (with respect to the triad (dn1,dn2,dn3)) of the n n+1 vectors dn+1 − xn. The components Dnij j , j = 1, 2, 3, and r = x n+1 = dni ·dn+1 of the vectors d form a 3 × 3 orthogonal matrix j j with its entries functions of the three angles, (θn1, θn2, θn3), called the tilt, roll, and twist. [See the left-hand side of Figure 2.] The displacement variables, (ρn1,ρn2, ρn3), called shift, slide, and rise [See the right-hand side of Figure 2.] are related to the components rni = rn · dni of rn by a coordinate transformation of the type

are taken to be homogeneous. I show the remarkable influence of the bending angle on the distribution of configurations. I first present the theoretical model and then discuss the bifurcation analysis, the resulting equilibrium configurations, and the Monte Carlo results for each of the two cases.

2. THEORY In this theory, as in refs 16−18, base pairs are represented by rigid planar objects, and the energy of a duplex DNA molecule with N base-pair steps is determined when there is given, for each n, both the location of the center of the rigid object that represents the n-th base pair and a right-handed orthonormal triad (dn1,dn2,dn3) that is embedded in the base pair as shown in the center of Figure 2. The assigned center of a base pair is here

ρin = R̃ ij(θ1n , θ2n , θ3n)r jn

Rnij

R̃ ij(θn1,

(2)

θn2,

θn3)

The numbers = are the components of an orthogonal matrix, and the functions R̃ ij are independent of n. The elastic energy of the n-th base-pair step is thus given by a function of six kinematical variables ψ n = ψ̌ n(θ1n , θ2n , θ3n , ρ1n , ρ2n , ρ3n )

It is a matter of importance in the molecular biology of DNA that the function ψ̌ n depends on the nucleotide composition of the n-th and (n+1)-th base pairs. It is generally assumed that ψ̌ n, as a function of (θn1,θn2,θn3,ρn1,ρn2,ρn3) alone, is independent of the nucleotide composition of other base pairs, e.g., the (n−1)-th and (n+2)-th, from which it follows that a base-pair step can be one of ten different types. In the present calculations this assumption is made; however, the theory remains valid in more general cases in which the elastic energy function is assumed to be influenced by the composition of base pairs other than the nth and the (n+1)-th. The kinematical variables θni and ρni were defined with precision by El Hassan and Calladine30 in 1995. A detailed description of this parametrization and the implied relations between the components Rij and Dij and (θn1, θn2, θn3) is given in the appendix of ref 17. As in ref 18, a ”closed” DNA molecule, i.e., a molecule that is ”circularized” so that its ends are joined and aligned in such a way that they form a single base pair and the axial curve of the molecule has the topology of a circle, is considered. It is assumed that no external forces or moments act on a base pair other than those that result from self-contact or the long-range electrostatic interaction of negatively charged phosphates in the same polymeric molecule. It is further assumed, as an approximation, that the two negative charges associated with a base pair are located at the barycenter of the base pair, which implies that the electrostatic energy Φ of the DNA molecule depends on the configuration through an equation of the form

Figure 2. A schematic drawing of the two adjacent base pairs forming the n-th base-pair step of a DNA molecule is shown in the center. Each nucleotide base in the n-th base pair is covalently bonded at its darkened corner to one of the two sugar phosphate backbone chains. The 5′-3′ direction of the left chain is upward. The gray-shaded long edges represent the side that is in the minor groove of the DNA. Shown on the left- and right-hand sides are schematic representations of the kinematical variables that describe the relative orientation and displacement of consecutive base pairs. Each drawing illustrates a case in which one of the kinematical variables has a positive value and the others (with the exception of ρ3) are set equal to zero. The case in the center drawing depicts the combined effect on the configuration of the shown values of the tilt, roll, twist, shift, slide, and rise that were set (just for this nonrealistic illustration) to be equal to 25°, 26°, 36°, 6 Å, 4 Å, and 10 Å, respectively.

referred to as the barycenter (center of mass) of the base pair. In refs 30 and 31, this position was taken to be the midpoint of the line connecting the C8 of the purine base to the C6 of the pyrimidine base. To avoid early inconsistencies in this important issue, a rigorous definition of the position and orientation of reference frames attached to a base pair was introduced as a standard in 2001.32 The piecewise linear curve * composed of the line segments that connect the spatial points xn for n = 1, 2, ..., N represents here the axial curve of the molecule. The elastic energy Ψ of a configuration is taken to be the sum over n of the energy ψn of interaction between the n-th and the (n+1)-th base pairs, i.e.

1 Φ= 2

∑ ϕn (4)

n=1

̃n

in which ϕ = ϕ (x ,x ,...,x ) is the electrostatic energy associated with the n-th base pair. In the present work Manning’s theory of charge condensation33 is employed. With this, the total energy U of the DNA molecule is the sum of its elastic energy Ψ and its electrostatic energy Φ

U=Ψ+Φ

1 2

N

(5)

In each configuration Ψ is given by eqs 1 and 3. The function ψ̌ n is assumed, as in the usual practice, to be a quadratic form in the excess tilt Δθn1, the excess roll Δθn2, the excess twist Δθn3, the excess shift Δρn1, the excess slide Δρn2, and the excess rise Δρn3. These quantities are defined, for i = 1, 2, 3, by the relations

∑ ψn n=1

N+1

n

N

Ψ=

(3)

(1)

where ψ , the elastic energy of the n-th base-pair step, is given by a function of the relative orientation and displacement of the n

θin = oθin + Δθin C

(6) DOI: 10.1021/acs.jctc.7b01090 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article


n-th base pair and the charges associated with base pairs other than the adjacent neighbors [It is assumed here that the local interaction between adjacent neighbors is included in the elastic part of the energy.] is employed:

and n

ρi =

ρn o i

+ Δρi

n oθ i

n

(7)

n oρi

in which and are intrinsic values, i.e., values appropriate to the stress free state of the n-th base-pair step. The stress-free state of a single base-pair state is regarded here as the relative displacement and orientation between the two base pairs comprising the step when they are free from any form of external force. Hence, the intrinsic parameters have the values characterizing the global minimum energy state of an unconstraint molecule. The local elastic energy can therefore be written as 3 n

ψ =

3

∑∑ i=1 j=1 3

+

1 n n n Fij Δθi Δθj + 2

3

ϕn = Q

m≠n,n±1

3

3

∑∑ i=1 j=1

1 n n n Fij Δθi Δθj 2

is given by the (11)

in which c is the concentration of (monovalent) salt in moles per liter. The constant Q is given by

GijnΔθinΔρjn

i=1 j=1

Q=

q2 4π ϵ0ϵw

(12)

with ϵ0 representing the permittivity of free space, and ϵw representing the dielectric constant of water. In accord with Manning’s theory33,36 q is set equal to 24% of the charge of the two phosphate groups associated with each base pair, i.e., q = 2 × 0.24e−, where e− is the charge of an electron. Our model18 accounts for the impenetrability of the DNA molecule using geometric inequality constraints. This enables us to investigate the equilibria and canonical distribution of highly supercoiled mincircles as functions of the ionic strength in the media. In other practical approaches the dependence of the intramolecular electrostatic interactions on ionic strength is addressed by setting the value of an effective diameter and using it as the impenetrable diameter of the molecule in accord with the Stigter’s theoretical prediction.37 The dependence of such effective diameter on the ionic strength was also determined by fitting simulation data to experimental results.38 As a DNA configuration is specified (up to a rigid body motion) by the 6N kinematical variables, and as it is assumed here that the shift, slide, and rise of each base-pair step are fixed, it is permissible to identify a configuration with the set of 3N components of α

(8)

The elastic moduli Fijn, Gijn, and Hijn, like the intrinsic parameters oθni and oρni , are constants that depend on the nucleotide composition of the n-th and (n+1)-th base pairs. The intrinsic values appropriate to the model were estimated in ref 6 by calculating the mean values of the kinematical variables from X-ray crystal structures of DNA protein complexes. The same data were used to relate the covariance matrix associated with the fluctuations of the kinematical variables to the inverse of the symmetric matrix of elastic moduli that quantifies the elastic energy of a single base-pair step as a quadratic form with respect to variations of the kinematical variables. The estimated values of the elastic moduli were properly normalized to yield DNA persistence length in accord with experimentally observed values. In a recent theoretical model the elastic energy associated with a base-pair step was taken to be a function of the relative position and orientation, not only between the two adjacent base pairs comprising the step but also between each of the two complementary bases in a base pair.34 To do so they34 regarded the possibility of the molecule to be prestressed, which, together with the assumed symmetry of the two antiparallel strands, may yield nonlocal sequencedependence of the minimum energy state of a base-pair step. On average, the translational deformations, i.e., changes in shift, slide, and rise, have a relatively minor effect on the configuration of a DNA molecule in comparison with angular deformation. On the other hand, in ref 35 it was shown that a slide mode of deformation may have a major effect on a doublehelical structure when coupled with roll. This properly phased coupling facilitates the tight wrapping of the DNA around the core of histone proteins.35 As in the present work I do not wish to investigate the slide-roll coupling effects; for the calculations presented here I took the assumption that the values of ρni are fixed and equal to their intrinsic values. This assumption reduced the required computational effort significantly and allowed me to investigate both sequence-dependent effects and protein-induced bending without worrying about this coupling. Hence, the local elastic energy is reduced to ψn =

(10)

κ = 0.329 c

3

2

dnm = |x n − x m|

The Debye screening parameter in units of Å formula

∑ ∑ 1 HijnΔρin Δρjn i=1 j=1

e−κdnm dnm

−1

3

∑∑

∑

α = (α1 , α2 , α3 , α4 , ..., α3N ) = (θ11 , θ21 , θ31 , θ12 , ..., θ3N ) (13)

A closed DNA molecule with N base pairs is here treated, as though it is one with N+1 base pairs that has its two terminal base pairs in coincidence (the (N+1)-th does not carry any charge), and hence obeys the end conditions x1 = x N + 1,

d1i = d iN + 1,

i = 1, 2, 3

(14)

Although it is generally a complicated expression, with the kinematical description outlined above one can express the equations in eq 14 as six nonlinear scalar equations of the general form Ya(α) = 0,

a = 1, ..., 6

(15)

An equilibrium configuration is defined as a configuration for which the first variation of the total energy U vanishes for any admissible variation in the kinematical variables (or in the positions and orientations of the orthonormal triads), i.e., a variation obeying the end conditions and for which selfpenetration is not present. The form of the general equilibrium equations was derived from this condition in ref 18 in which an algorithm for solving the system was introduced. Such equilibrium configuration is stable if the second variation of

(9)

In accord with Manning’s theory of charge condensation33 the following expression for the energy resulting from the electrostatic interaction between the charge associated with the D


Article


basis available, one can randomly select a vector Δβ that lies in the null space of J and has a norm not greater than χ, which represents the maximum allowed move length. [The generation of this random move must be taken with care so that there is a uniform probability to move to each point in the hypersphere.] Thus, I have now 6 random components, ξa, a = 1, 2, ..., 6, of the move

the total energy is positive for any admissible variation. In ref 18 a detailed description of the second variation in terms of an analytically calculated symmetric Hessian matrix H was given. The smallest eigenvalue, λ, of H must be positive for stable equilibria. A complete description of the method used to find equilibrium configurations and to determine their stability is described in ref 18 and used here for the calculations of equilibrium configurations.

6

Δβ =

3. THE MONTE CARLO SCHEME: A GEOMETRICALLY CONSTRAINED METROPOLIS ALGORITHM To calculate a canonical ensemble of configurations obeying the Boltzmann distribution, I iteratively change a configuration using a move that is acceptable, i.e., in accord with (14), under strict detailed balance and ergodicity conditions39,40 [Detailed balance and ergodicity conditions ensure unbiased sampling and possible explorations of the complete configurational space.], calculate its total energy, and apply the acceptance criteria of the Metropolis scheme. Configurations with selfpenetration were properly rejected. To verify this excluded volume effect I used the values of the intramolecular distances dnm needed for the calculation of the electrostatic energy, to calculate analytically the surface of a tube with a diameter of 20 Å.18 Suppose that prior to a move one has a closed configuration, i.e., a configuration satisfying (15) with no self-penetration. Of course, a change in the values of the kinematical variables may result in a violation of these imposed constraints. In each iteration in the present scheme I start with a closed configuration and select randomly four base-pair steps, 1 ≤ n1 < n2 < n3 < n4 ≤ N. I then change only the 12 corresponding n n n n n angular variables, {θ11, θ21, θ31, θ12, ..., θ34}, while holding all the remaining variables fixed. One can now write the six equations in eq 15 as functions of only 12 unknowns, so that the subsegments between the chosen base-pair steps are rigid links for that move. [This makes the treated problem similar to the problem of forward kinematics of robotic manipulators with 12 degrees of freedom.] A linearization of the equations about a closed configuration yields Y(β + Δβ) = JΔβ + ( Δβ 2 )

a=1

Y(β + Δβ + JT σ ) = 0

β new = β + Δβ + JT σ

(21)

which satisfies the constraints (14) to within machine accuracy. Because the random move was performed in the Null space of J, not more than 3 iterations are needed for the calculations of the six components of σ in each complete Metropolis move. For a range of values of χ appropriate for the Metropolis Monte Carlo algorithm used here there exists a 1-1 mapping from the configuration β + Δβ to βnew. In order to make sure that this is the case, the complete move was confined to have norm smaller than χ, i.e.

(16)

β new − β < χ

(22)

This condition and the uniform distribution of Δβ imply that, for χ small enough, the probability p(β → βnew) to move from β to βnew equals p(βnew → β). [A schematic description of the scheme for a three-dimensional problem with two constraints can be seen in Figure 3 of ref 20 where the method was implemented on more general constraints for the calculation of looping probabilities.] The fact that the complete move can cover the 12-dimensional configuration space and that for each iteration 4 base-pair steps are randomly chosen implies ergodicity. The conditions (14) do not completely determine the topology of a closed molecule. For one can, in theory, nick a single strand of the molecule, twist by a complete turn the basepair on one side of the nick while holding fixed the base pair on the other side, and then join together the two ends of the nicked strand to form a closed molecule for which the constraints (14) still hold. Such a procedure changes by one the linking number, Lk, of the molecule. The linking number is a topological invariant, an integer. It is equal to the number of times that (either) one of the two (closed) strands is linked to the axial curve. [As the molecule can be closed in such a way

(17)

the components of Y are (Y1, Y2, ..., Y6), i.e., the values of the left-handed side of the constraints (15), and the 6 × 12 matrix J is the Jacobian matrix characterizing the linearization of (15), i.e. a = 1, ..., 6,

(20)

Since I have in hand the random move Δβ, the relation 20 is a set of six nonlinear equations for the six unknown components of σ. Implicit function theorem and the smoothness of the hypersurface characterized by (14) assures the uniqueness of σ for a small enough value of Δβ . In the calculations performed here I set that condition by taking χ to be χ ≤ 5 π/180. This restricts a single change in each angular variable to be of a magnitude smaller than 5 degrees. The system is solved using a modified Newton−Raphson scheme with the associated Jacobian calculated analytically for each iteration. With this I find a new configuration, namely

β = {β1 , β2 , β3 , β4 , ..., β12} = {θ1n1 , θ2n1 , θ3n1 , θ1n2 , ..., θ3n4}

∂Ya , ∂βm

(19)

The final move is necessary to eliminate the second order linearization error. This is done by moving along the space complement to the null space of J, i.e., the range of (JT), spanned by the 6 rows of J. This means that the correction move must satisfy the relation

where the components of the random move β are

Jam =

∑ ξa ba

m = 1, ..., 12 (18)

One now attempts to find a finite nonzero variation Δβ for which the right-hand side of (16) vanishes. A variation that satisfies the linear part of eq 16, ignoring for now the second order error, must be in the six-dimensional null space of J. Based on this, I perform a random move within the boundaries of a hypersphere of radius χ centered at the original point in the 12-dimensional configuration space. For this, I calculate in each iteration a basis of six orthonormal 12-dimensional vectors, ba, a = 1, 2, ..., 6, that span the null space of J (each vector in the basis is perpendicular to each of the rows of J). Having such a E


Article

Journal of Chemical Theory and Computation Table 1 DNA sequence TTTATACTAA 10 ACGCGCGAGG 60 CGGTGGAATA 110 GTACCGAGTT 160 GGTCACTAAT 210 GTTTTACAGT 260 ATTGATATTT 310

CTTGAGCGAA 20 CAGCTGTATG 70 TTTCGTTTCC 120 CCGACACTTT 170 ACCATCTAAG 220 ATTATGTAGT 270 ATATCATTTT 320

ACGGGAAGGG 30 GCGAAATGAA 80 TACTACGACT 130 CATTGAGAAA 180 TAGTTGATTC 230 CTGTTTTTTA 280 ACGTTTCTCG 330

TATCACCGAA 50 CCGGAAAACG 100 GGAAGCCTAT 150 CTCTGTTACA 200 CATATGTTGT 250 AATTTAATAT 300

defined to account for sequence-dependent interactions that stabilize the double helix was introduced in refs 50 and 51. An extension of the worm-like chain model that includes spontaneous kink formation was introduced in ref 52 in order to explain the experimental results of Cloutier and Widom.53 [Cloutier and Widom53 used DNA cyclization assays to show that even short DNA segments of 94 base pairs have an unexpectedly large cyclization J-factor because of spontaneous sharp bending of the DNA segments. These results are still considered controversial, see, e.g., ref 54.] Recent experimental results55 suggested that such spontaneous kinks form for loops of length 60 base pairs or less. My current model does not take into account the effect of spontaneous kink formation but can be extended to account for such local deformations.

that its axial curve has the topology of a knot, not only its linking number but also its knot type must be specified.] The linking number of a closed molecule is given by the relation [For general definitions see refs 41−43.]

L k = Wr + Tw

TTTTCACCGA 40 AGAGTTCTTC 90 ACTATCAGCC 140 GATGCCTCAG 190 ATAGTGACTG 240 TGCAAAATCT 290 TTCAGCTTT

(23)

in which the writhe, Wr, is a measure of the chiral distortion of the axial curve from planarity, and the total twist, Tw, is the number of times the one DNA strand is wound about the axial curve. The scheme described above for the generation of closed configurations would not in general preserve the linking number, as in some rare events the change in the kinematical parameters may result in a self-crossing of two subsegments of the molecule, and by that change the linking number. For this reason in the Monte Carlo simulation the linking number was numerically calculated once in about 100 moves. To do so the writhe, Wr, was calculated according to the scheme proposed in appendix B of ref 44 and the total twist (in turns), Tw, was calculated using the formula proposed by ref 45, which is proven to yield the exact value. Other methods suitable for studies of the statistical mechanics of DNA molecules of length in the order of a few hundreds to several thousands of base pairs using the Gaussian sampling approach were recently presented.46−49 In contrast to a Metropolis algorithm, the Gaussian sampling method has the ability to efficiently generate a very large canonical ensemble (up to 1014 configurations) of uncorrelated configurations. However, a significant drawback of this method is that, in order to simulate a molecule that is subject to geometric constraints such as (14), one needs to relax the geometric constraints and to make use of a small fraction of configurations obeying them. [The typically large canonical ensemble generated using Gaussian sampling permits one to estimate cyclization J-factors and looping propensities directly from the fraction of configurations obeying the relaxed constraints.] Furthermore, this method cannot be utilized when other contributions to the quadratic elastic energy, such as the intramolecular electrostatic energy of a DNA molecule in solution, or even nonquadratic terms in the elastic energy, are taken into account. The geometrically constrained Metropolis algorithm enables one to generate a large canonical ensemble of configurations subject to any possible set of geometric and topological constraints. A very interesting approach to an efficient calculation of ring closure probability in which a mesoscopic Hamiltonian is

4. SIMULATION AND RESULTS In this section I show the simulation results for two 339 bp minicircles. The first is constructed according to a given sequence but does not include protein-induced curvature. The second was taken to be intrinsically straight and homogeneous, but for which I mimic a general protein effect by taking the 11 base-pair step subsegment, with its center at step n = 170, to be rigid and to have a given curvature in its middle step (by controlling the value of the intrinsic roll, oθ170 2 ). 4.1. Effect of Inhomogeneity. In their research Fogg, Kolmakova, Rees, Magonov, Hansma, Perona, and Zechiedrich29 were able to generate and isolate large quantities of several distinct topoisomers, of a closed 339 bp DNA molecule labeled here A339 [DNA topoisomers are molecules with identical sequence that differ only in their linking number.]. The sequence in each of these 339 bp topoisomers can be expressed as a series of the four nucleotide bases A, T, C, G (Adenine, Thymine, Cytosine, Guanine) in one of the two Watson−Crick strands. In the present case that sequence is shown in Table 1. As the molecule is closed, the first and last base pairs form a base-pair step, and hence an additional subsegment of 6 T’s is formed. The topoisomers obtained in the experiments reported in ref 29 were examined using gel electrophoresis and atomic force microscopy (AFM). The results suggested the existence of topoisomers with linking numbers varying from Lk = 32 down to highly supercoiled minicircles with Lk = 26. A similar minicircle with 336 base pairs was used recently as a good representative of the supercoiled DNA loops found in nature.56 This length was found to be appropriate because it has about 10 different F


Article


Figure 3. Total energy U of the minicircle A339 at c = 1 × 10−1 M as a function of ζ3/2π. The value of ζ3/2π is essentially the excess link. The angle ζ is the relative rotation angle between the (virtual) (N+1)-th and the 1-st base pairs (that are constrained to lie on the same plane with coinciding barycenters).

Figure 4. Equilibrium configurations at c = 1 × 10−1 M of seven topoisomers of A339 with linking numbers Lk = 29, 30, ..., 35. Two perpendicular views of the intrinsic configuration are shown in the upper right corner of the figure. Each equilibrium configuration is drawn in two perpendicular views. Subsegments comprised of more than three adjacent A’s (Adenine) and T’s (Thymine) are drawn in blue and red. Subsegments that are in contact with sequentially remote base-pairs are shown in black. The arrows indicate the position of base pair 1 and the 5′-3′ direction (growing n). A lengthy search of other equilibrium configurations suggests that the shown configurations are the minimum energy configurations of A339 confined to the indicated linking numbers. The energy values are given in kBT, with kB representing the Boltzmann constant and T representing the absolute temperature which was set to 300 K.

distinguishable topoisomers.56 I calculated equilibrium configurations of the minicircle molecule treated in ref 29, for several topoisomers, i.e., minicircles with several values of Lk. [An analysis of Lac repressor-partitioned DNA rings, based on minimum elastic energy, of a shorter minicircle (182 bp) was

reported recently in ref 57.] The calculated configurations show a remarkable agreement with the representative configurations obtained in ref 56 using electron cryotomography (cryo-ET). In a very recent work the free energy of contact formation between base pairs of 336 bp minicircles (with the same G


Article


Figure 5. Four color maps showing the spatial J-factor associated with a couple (m,n) of remote base pairs (ordered in accord with the given sequence) in mole per liter for four different values of the radius d of the spherical volume. Representative configurations (with Lk = 29) of each of the two probable states are given. On the right the configurations correspond to the regions on the upper right of each color map, and on the left the three configurations correspond to the region on the lower left. In each of the configurations the first base pair is shown together with an arrow pointing to the direction of growing n along the minicircle. The couple (m,n) of base pairs that have the closest proximity are marked with their numbers and the distance between the associated barycenters. As in Figure 4, subsegments comprised of at least 4 adjacent A’s (Adenine) and T’s (Thymine) in a single strand are drawn in blue and red. In each configuration all base pairs that have other sequentially remote base pairs to within a distance of 30 Å (right-hand side) or 60 Å (left-hand side) are marked in yellow.

energy, U, on ζ3/2π along the equilibrium path, that, based on a thorough study, is very much likely to contain the minimum energy configurations for each value of ζ3/2π, is shown in Figure 3. [The shown equilibrium path contains at least one stable configuration for each value of ζ3/2π.] The graph shows that, as was suggested in 29, Lk = 32 for the ”most relaxed” closed configuration of A339. The (presumably) minimum energy configurations of A339 with seven linking numbers varying from 29 to 35 are shown in Figure 4. The configurations with Lk = 29 and Lk = 35 have points of self-contact associated with base-pair steps that are shown in black. Two sequentially remote sites are considered to be in contact when the distance of closest approach between the associated axial curves is exactly 20 Å. [For precise definitions of the contact points and the assumed geometry of the surface of the DNA in the present model see ref 18.] Contrary to the previous examples in which the treated molecules were assumed homogeneous and therefore the values of the total energy of their equilibrium configurations are practically indifferent to slithering modes, for the present example, the equilibrium configurations are highly sensitive to such modes. As a result, the proximity and relative orientation of sequentially remote sites (that might be in contact) is highly sequence-dependent. For the topoisomer with Lk = 29 there are

sequence for n = 1, 2, ..., 336) was calculated using molecular dynamic simulation for 10 different values of excess link.58 A proper use of a Morse potential was made in order to account for sequence-dependent effects, and the results appear to confirm that there are pairs of remote subsegment that have lower free energy of contact formation.58 All the results in this example were calculated at c = 1 × 10−1 M. The values of the intrinsic parameters and elastic moduli associated with each of the 10 unique base-pair steps are in accord with the database indicated in Zheng, Colasanti, Lu, and Olson59 which were derived by statistical analysis of X-ray crystallography data.6,60 Starting with end conditions that are close enough to those of the intrinsic configuration, a series of equilibrium configurations of open molecules were calculated by successively varying the end conditions such that the final calculated configuration in the series is closed, i.e., obeys eq 14. With a closed equilibrium configuration in hand, the twist angle between the first and (N+1)-th (virtual) base pairs, ζ3, was varied (as a bifurcation parameter) using standard continuation methods to follow a single equilibrium path. In such an equilibrium path, the number ζ3/2π, which may be regarded as the excess link, attains integral values that correspond to different linking numbers Lk. The dependence of the total H


Article

Journal of Chemical Theory and Computation two points of contact. The first is located between the central base-pair steps in the subsequences TTTCACCG and TAGTGACT that are 139 base pairs apart. The second point of contact is located between the central base-pair step in the subsequence GAGGCAGC and the base pair in the center of the subsequence AATACCA that is 151 base pairs apart. For the topoisomer with Lk = 35 there is only a single point of contact. It is located between the central base-pair steps in the subsequences TCACCGAT and TGCCTCAG that are 149 base pairs apart. The different number of contact points between the negatively and positively supercoiled minicircle is a consequence of the assumed intrinsic geometry of the sequence comprising the molecule. In Figure 4 subsegments of at least four consecutive A’s or T’s are marked in order to show the existence of A-tracts in the investigated sequence. A-tracts are believed to have a molecular structure that gives rise to global curvature of the DNA double helix.61 They have relatively high intrinsic bending together with high bending stiffness, and their elastic response and intrinsic geometry were found to be sensitive to modes of bending deformations.62 The calculated equilibria may imply that contact occurs with high likelihood between the two subsegments that are in close proximity in the minimum energy configuration; however, entropic effects may dominate. To check this, I performed a Monte Carlo simulation and calculated the canonical ensemble of 109 configurations for which Lk = 29 with the complete siteto-site statistical information. With the results it was possible to find the J-factor, J0, that expresses the concentration of one site in a sphere of radius d, centered at the other site. The calculation of J0 is in accord with the spatial part of the J-factor discussed in ref 63 J0 =

p(dmn < d) , NAvVd

Vd =

4 3 πd 3

Figure 6. Probability distribution function showing the probability of the couples (57,213) in solid black and (132,288) in dotted red to be in distance less than d.

(24) Figure 7. Probability density function showing the probability of the couples (57,213) in solid black and (132,288) in dotted red to be in distance within a small range around d.

where NAv is Avogadro’s number, and p(dmn < d) ≃

N (dmn < d) NTotal

(25)

roll. However, the intrinsic twist of a DNA in its B form makes a long enough molecule to respond mechanically as a nearly transversely isotropic molecule. To make the current model as real as possible, and yet homogeneous, the intrinsically straight, 339 base-pair minicircle, labeled here /339, was assumed to have bending rigidity for roll that is half that for tilt, Fn22 = 0.5Fn11, with the harmonic mean, FnH, of Fn11 and Fn22 equal to 0.0427 kBT/deg2 for all the base-pair steps. This yields a persistence length of 476 Å.64 Hence, for the molecule /339 treated here it was assumed that the diagonal elastic moduli are

the probability of the couple of base-pair steps m and n to be in distance dmn < d was calculated as the ratio between the number N(dmn < d) of configurations for which this condition is valid and NTotal, the total number of configurations in the canonical ensemble. The results together with representative configurations are shown in Figure 5. One can clearly see that there are two competing pairs of regions of contact. This suggests that DNA architecture can drive a supercoiled molecule to slither between multiple states bringing specific sites to close proximity. To have a closer look at two particular couples of base pairs the calculated probability distribution function P(dmn ≤ d) for the couples (57,213) and (132,288) is shown in Figure 6. According to these calculations the two base pairs in each couple have the highest probability to be in contact. Each point on the graph shows the probability of the base pairs n in each couple to be in a spherical volume of radius d centered at the barycenter of base pair m. As this radius increases, this probability approaches unity. The probability density function, i.e., the derivative of P(dmn ≤ d) with respect to d, shown in Figure 7, reveals two states with high probability for each of the couples. 4.2. Protein-Induced Curvature Effect. A DNA molecule has, in general, significantly higher resistance to bending deformations that distort the tilt, than to those that distort the

FHn = 0.0427kBT /deg 2,

n n F11 = 2F22 = 1.5FHn ,

n F33 = 1.05FHn

(26)

and that there is no coupling between the modes of deformations, i.e., Fij = 0 for i ≠ j. The intrinsic kinematical values of /339, for n = 1, 2, ·..., N, were set to have the following values: 2π n n n , o ρ1n = 0, oθ1 = 0, oθ2 = 0, oθ3 = 10.6 ρn o 2

= 0,

ρn o 3

= 3.4 Å

(27)

The molecule treated in the simulations was closed with Lk = 29, and the protein-induced bending was considered by I


Article


Figure 8. A bifurcation diagram showing the distance dcm between the centroid of the protein bound subsegment (with base pairs 165 ≤ n ≤ 176) and that of each calculated equilibrium configuration as a function of the bending angle β. Three enlargements of the framed rectangular regions are shown together with 12 representative examples of equilibrium configurations. The location of each configuration is marked with blue circles on the associated branch. The red dots show bifurcation points or folds. The bent subsegment is depicted in blue on each of the configurations and with a bold arc next to it.

shown in the enlarged boxes. These paths fold around β = 0 because a change in the sign of the bending angle gives rise to a significant mode of eversion. The calculated values of the total energy U, the distance dcm, and the induced bending angle, β, for each of the shown configurations, are indicated in Table 2.

assuming that the roll for n = 170 is fixed and equal to the bending angle β. All the other kinematical parameters for 165 ≤ n ≤ 175 were held fixed in their intrinsic values. Since the treated molecule is homogeneous, the sequential position of the bent subsegment was chosen arbitrarily. A thorough bifurcation analysis was done in order to reveal the dependence of stable equilibria on the protein-induced bending angle β. The results are depicted in Figure 8. For each calculated equilibrium configuration, the bifurcation diagram shows the distance, dcm, between the centroid [the arithmetic mean of the position vectors xn in the protein-bound subsegment] of the proteinbound subsegment (with base pairs 165 ≤ n ≤ 176) and the centroid of the complete molecule at the equilibrium configuration. The stability of each equilibrium was determined by calculating the minimal eigenvalue λ of the Hessian matrix H that characterizes the second variation of the total energy at equilibrium.18 The results show two branches containing stable equilibria, i.e., paths with positive λ. A branch of symmetric configurations, in which the bent region is located close to the centroid of the molecule’s configuration, contains stable configurations for bending angles β < 11.6° and has unstable equilibria in the range 11.6° < β < 56.8° (See the unstable configuration E in Figure 8.). For β > 56.8° the symmetric equilibria on this branch are again stable. The shown stable configurations K and H are in this branch. From the second bifurcation point of this branch, in the range −1.1° < β < 56.8° appears a second branch with equilibrium configurations that are stable. The equilibria on the branch show a shift of the bent site away from the contact region in the center (See the stable configurations A, B, and D on this branch.) This shift is more significant as β decreases in this range. All paths of stable equilibria appear in red and are marked with arrows between the relevant bifurcation points. Twelve configurations are shown in Figure 8 together with their associated points on the diagram. Complex regions on the diagram, around β = 0, are

Table 2. States and Values of the 12 Equilibrium Configurations Shown in Figure 8 conf

state

U [kBT]

dcm [Å]

β [deg]

A B C D E F G H I J K L

stable stable unstable stable unstable unstable unstable stable unstable unstable stable unstable

132.2 131.1 135.3 129.7 129.7 137.9 140.4 128.9 140.3 135.2 129.9 133.0

99.5 35.3 30.3 40.1 36.9 36.2 54.4 51.5 49.7 99.5 0.3 18.9

0 25 25 50 50 50 70 70 70 70 0 0

My Monte Carlo results for six different cases of the supercoiled molecule, protein free, zero bending (with rigid binding site), and bending angles β equal to 15, 30, 45, and 60 degrees, are shown in Figure 9. In each case the calculated J0 for d = 25 Å is shown. The results show that for β greater than 30 degrees the molecule fluctuates primarily around configurations in which the bent region is away from the contact region, and this becomes clearer as the bending angle is increased. Thus, while for β = 0 the protein induced stiff subsegment tends to be in the contact region, for greater bending angles, around β = 30°, the bent subsegment appears both in the contact region and in the highly curved sides, and for higher bending angle, the J


Article


Figure 9. Six color maps showing the spatial J-factor that gives the concentration of base pair m in a sphere of radius 25 Å around base pair n. Two very different representative configurations are shown for each of the cases with β = 0, 30, and 60 degrees. The protein-induced bent segment, for which 165 ≤ n ≤ 176, is marked in blue in each of them.

In order to investigate the effect of protein-induced bending, a study of a homogeneous DNA molecule was performed. By comparing the results shown in Figure 8 with those of Figure 9, one can see how thermal fluctuations give rise to two competing states, one with the bent site away from the centroid of the configuration and one in which the bent site is in close proximity to it. For high induced curvature the presented bifurcation analysis suggests that the state with the bent subsegment in the region of contact is the stable global minimum energy state; however, the colormap on the lower right-hand side of Figure 9 shows that this position is of very small likelihood. The colormap describing the case with the induced bending angle β = 0 shows that the assumed rigid segment slithers from the contact region up to the top or bottom of each of the two loops but essentially avoids the highly bent center of the side loops. When the induced bending is 30°, the bent segment can slither from the region of contact to each of the looped sides. When the induced bending angle is increased further to 45° or 60°, the resulting colormaps show high tendency of the protein-induced bent segment to be in the looped subsegments away from the contact region. Although the existence of an equilibrium configuration in which the curved segment is positioned along the side loops is possible, the fact that the minimum energy state is so far away from the state of maximum likelihood suggests that the slithering is driven by an entropic force. When the bent subsegment is in the center along the contact region, excluded volume effects reduce significantly the number of available microscopic states, and the system then is driven to the competing state for which the bent subsegment is shifted to one of the side loops. Similar entropic effects that drive a highly supercoiled minicircle from compact states with high writhe to open states with high total twist and low writhe were discussed in the molecular dynamics

bent segment tends to have maximal distance from the contact region. It is suggested that the combination of high bending together with high excess link gives rise to a significant entropic force that shifts by slithering the highly bent region away from the contact sites. Since the simulated DNA is homogeneous, the results must be indifferent to slithering and eversion modes for the protein-free minicircle, and hence the spatial J-factor must give the uniform strip shown on the upper left-hand side of Figure 9. In all the cases shown in Figure 9 couples of base pairs that are sequentially separated by not more than 40 basepair steps are excluded.

5. CONCLUSIONS The effect of protein-induced bending on supercoiled minicircles DNA molecule was studied by the two unique schemes for the calculation and analysis of equilibria and the generation of canonical ensemble of configurations. I first isolated the sequence dependence effects and showed that the sequence itself has a dramatic effect on the probability of contact between sites in a supercoiled DNA minicircle. As shown in Figure 4 the assumed intrinsic geometry of the molecule has high curvature around n = 1 and around n = 200 in nearly opposite directions. For the calculated minimum energy configurations the later subsegment has its highly bent center in the contact region for the two cases with Lk = 29 and Lk = 35. This position explains the lower left patch in each of the colormaps shown in Figure 5 and the three representative configurations shown on the left of the figure. The competing state is the case for which the two intrinsically bent subsegments are located away from the contact region partly composing the two side loops. This is the main character of the three representative configurations shown on the right-hand side of Figure 5. K


Article


(14) Coleman, B. D.; Swigon, D. Theory of supercoiled elastic rings with self-contact and its application to DNA plasmids. J. Elasticity 2000, 60, 173−221. (15) Lavigne, M.; Kolb, A.; Yeramian, E.; Buc, H. CRP fixes the rotational orientation of covalently closed DNA molecules. EMBO J. 1994, 13, 4983−4990. (16) Coleman, B. D.; Olson, W. K.; Swigon, D. Theory of sequencedependent DNA elasticity. J. Chem. Phys. 2003, 118, 7127−7140. (17) Biton, Y. Y.; Coleman, B. D.; Swigon, D. On bifurcation of equilibria of intrinsically curved, electrically charged, rod-like structures that model DNA molecule in a solution. J. Elasticity 2007, 87, 187− 210. (18) Biton, Y. Y.; Coleman, B. D. Theory of the influence of changes in salt concentration on the configuration of intrinsically curved, impenetrable, rod-like structures modeling DNA minicircles. International Journal of Non-Linear Mechanics 2010, 45, 735−755. (19) Biton, Y. Y.; Coleman, B. D. Applications of a theory of sequence-dependent DNA elasticity that accounts for intramolecular electrostatics forces. Continuum Models and Discrete Systems CMDS 11 2008, 359−382. (20) Biton, Y. Y.; Kumar, S.; Dunlap, D.; Swigon, D. Lac Repressor Mediated DNA Looping: Monte Carlo Simulation of Constrained DNA Molecules Complemented with Current Experimental Results. PLoS One 2014, 9, e92475. (21) Harteis, S.; Schneider, S. Making the Bend: DNA Tertiary Structure and Protein-DNA Interactions. Int. J. Mol. Sci. 2014, 15, 12335−12363. (22) Friedman, A.; Fischmann, T.; Steitz, T. Crystal structure of Lac repressor core tetramer and its implications for DNA looping. Science 1995, 268, 1721−1727. (23) Romanuka, J.; Folkers, G. E.; Biris, N.; Tishchenko, E.; Wienk, H.; Bonvin, A. M.; Kaptein, R.; Boelens, R. Specificity and Affinity of Lac Repressor for the Auxiliary Operators O2 and O3 Are Explained by the Structures of Their Protein-DNA Complexes. J. Mol. Biol. 2009, 390, 478−489. (24) Lewis, M.; Chang, G.; Horton, N.; Kercher, M.; Pace, H.; Schumacher, M.; Brennan, R.; Lu, P. Crystal Structure of the Lactose Operon Repressor and Its Complexes with DNA and Inducer. Science 1996, 271, 1247−1254. (25) Swinger, K. K.; Lemberg, K. M.; Zhang, Y.; Rice, P. A. Flexible DNA bending in HU-DNA cocrystal structures. EMBO J. 2003, 22, 3749−3760. (26) Fisher, D. E.; Parent, L. A.; Sharp, P. A. Myc/Max and other helix-loop-helix/leucine zipper proteins bend DNA toward the minor groove. Proc. Natl. Acad. Sci. U. S. A. 1992, 89, 11779−11783. (27) Czapla, L.; Peters, J. P.; Rueter, E. M.; Olson, W. K.; Maher, L. J. Understanding Apparent DNA Flexibility Enhancement by HU and HMGB Architectural Proteins. J. Mol. Biol. 2011, 409, 278−289. (28) Czapla, L.; Grosner, M. A.; Swigon, D.; Olson, W. K. Interplay of Protein and DNA Structure Revealed in Simulations of the lac Operon. PLoS One 2013, 8, e56548. (29) Fogg, J. M.; Kolmakova, N.; Rees, S.; Magonov, I.; Hansma, H.; Perona, J. J.; Zechiedrich, E. L. Exploring writhe in supercoiled minicircle DNA. J. Phys.: Condens. Matter 2006, 18, S145−S159. (30) El Hassan, M. A.; Calladine, C. R. The assessment of the geometry of dinucleotide steps in double-helical DNA: a new local calculation scheme. J. Mol. Biol. 1995, 251, 648−664. (31) Lu, X.-J.; Hassan, M. E.; Hunter, C. Structure and conformation of helical nucleic acids: analysis program (SCHNAaP). J. Mol. Biol. 1997, 273, 668−680. (32) Olson, W. K.; Bansal, M.; Burley, S. K.; Dickerson, R. E.; Gerstein, M.; Harvey, S. C.; Heinemann, U.; Lu, X. J.; Neidle, S.; Shakked, Z.; Wolberger, C.; Berman, H. M. A standard reference frame for the description of nucleic acid base-pair geometry. J. Mol. Biol. 2001, 313, 229−237. (33) Manning, G. S. Limiting laws and counterion condensation in polyelectrolyte solutions: I. Colligative properties. J. Chem. Phys. 1969, 51, 924−933.

simulations performed in ref 65. My results suggest that the statistical tendency of the molecule to slither toward one of these two states can be controlled by the magnitude of the protein-induced bending. This may be part of an adjustable mechanism used by binding proteins to harness the superhelical density and the sequence architecture of DNA in order to bring sequentially far sites to close proximity by bending the DNA at the bound subsegment.

■

AUTHOR INFORMATION

Corresponding Author

*Phone: +972 (0)58 6000595. E-mail: [email protected]. ORCID

Yoav Y. Biton: 0000-0002-5907-5125 Funding

This research was supported by a grant from the Research and Development Department, SCE, Shamoon College of Engineering, Beer Sheva, Israel. Notes

The author declares no competing financial interest.

■

ACKNOWLEDGMENTS The author thanks Professor Bernard D. Coleman for inspiring this research and Professor David Swigon for useful discussions on the theory and possible effects of sequence dependence DNA mechanics.

■

REFERENCES

(1) Drew, H. R.; Travers, A. A. DNA bending and its relation to nucleosome positioning. J. Mol. Biol. 1985, 186, 773−790. (2) Craig, N. L. The Mechanism of Conservative Site-Specific Recombination. Annu. Rev. Genet. 1988, 22, 77−105. PMID: 3071260. (3) Semsey, S.; Virnik, K.; Adhya, S. A gamut of loops: meandering DNA. Trends Biochem. Sci. 2005, 30, 334−341. (4) Han, L.; Garcia, H. G.; Blumberg, S.; Towles, K. B.; Beausang, J. F.; Nelson, P. C.; Phillips, R. Concentration and Length Dependence of DNA Looping in Transcriptional Regulation. PLoS One 2009, 4, e5621. (5) Shapiro, T. A.; Englund, P. T. The structure and replication of kinetoplast DNA. Annu. Rev. Microbiol. 1995, 49, 117−143. (6) Olson, W. K.; Gorin, A. A.; Lu, X. J.; Hock, L. M.; Zhurkin, V. B. DNA sequence-dependent deformability deduced from protein-DNA crystal complexes. Proc. Natl. Acad. Sci. U. S. A. 1998, 95, 11163− 11168. (7) Swigon, D. In Mathematics of DNA Structure, Function and Interactions; Benham, C., Ed.; Springer: Berlin, 2009; pp 293−320, DOI: 10.1007/978-1-4419-0670-0. (8) Tobias, I.; Olson, W. K.; Coleman, B. D. The dependence of DNA tertiary structure on end conditions: Theory and implications for topological transitions. J. Chem. Phys. 1994, 101, 10990−10996. (9) Furrer, P. B.; Manning, R. S.; Maddocks, J. H. DNA Rings with Multiple Energy Minima. Biophys. J. 2000, 79, 116−136. (10) Balaeff, A.; Koudella, C. R.; Mahadevan, L.; Schulten, K. Modelling DNA loops using continuum and statistical mechanics One contribution of 16 to a Theme ’The mechanics of DNA’. Philos. Trans. R. Soc., A 2004, 362, 1355−1371. (11) Antman, S. Nonlinear Problems of Elasticity; Springer-Verlag: New York, 1995; DOI: 10.1007/0-387-27649-1. (12) Tobias, I.; Swigon, D.; Coleman, B. D. Elastic stability of DNA configurations: I. General theory. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 2000, 61, 747−758. (13) Coleman, B. D.; Swigon, D.; Tobias, I. Elastic stability of DNA configurations: II. Supercoiled plasmids with self-contact. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 2000, 61, 759−770. L


Article


structure of minicircles under torsional stress. Nucleic Acids Res. 2017, 45, 7633−7642. (59) Zheng, G.; Colasanti, A. V.; Lu, X. J.; Olson, W. K. 3DNALandscapes: a database for exploring the conformational features of DNA. Nucleic Acids Res. 2010, 38, D267−274. (60) Colasanti, A. V. Conformational states of double helical DNA. Ph.D. Thesis, Rutgers University, 2006. (61) Haran, T. E.; Mohanty, U. The unique structure of A-tracts and intrinsic DNA bending. Q. Rev. Biophys. 2009, 42, 41−81. (62) Drsata, T.; Spackova, N.; Jurecka, P.; Zgarbova, M.; Lankas, F.; Sponer, J. Mechanical properties of symmetric and asymmetric DNA A-tracts: implications for looping and nucleosome positioning. Nucleic Acids Res. 2014, 42, 7383−7394. (63) Levene, S. D.; Crothers, D. M. Ring closure probabilities for DNA fragments by Monte Carlo simulation. J. Mol. Biol. 1986, 189, 61−72. (64) Olson, W. K.; Swigon, D.; Coleman, B. D. Implications of the dependence of the elastic properties of DNA on nucleotide sequence. Philos. Trans. R. Soc., A 2004, 362, 1403−1422. (65) Mitchell, J. S.; Harris, S. A. Thermodynamics of Writhe in DNA Minicircles from Molecular Dynamics Simulations. Phys. Rev. Lett. 2013, 110, 148105.

(34) Gonzalez, O.; Petkeviciute, D.; Maddocks, J. H. A sequencedependent rigid-base model of DNA. J. Chem. Phys. 2013, 138, 055102. (35) Tolstorukov, M. Y.; Colasanti, A. V.; McCandlish, D. M.; Olson, W. K.; Zhurkin, V. B. A Novel Roll-and-Slide Mechanism of DNA Folding in Chromatin: Implications for Nucleosome Positioning. J. Mol. Biol. 2007, 371, 725−738. (36) Fenley, M. O.; Manning, G. S.; Olson, W. K. Approach to the limit of counterion condensation. Biopolymers 1990, 30, 1191−1203. (37) Stigter, D. Interactions of highly charged colloidal cylinders with applications to double-stranded DNA. Biopolymers 1977, 16, 1435− 1448. (38) Rybenkov, V. V.; Cozzarelli, N. R.; Vologodskii, A. V. Probability of DNA knotting and the effective diameter of the DNA double helix. Proc. Natl. Acad. Sci. U. S. A. 1993, 90, 5307−5311. (39) Manousiouthakis, V. I.; Deem, M. W. Strict detailed balance is unnecessary in Monte Carlo simulation. J. Chem. Phys. 1999, 110, 2753−2756. (40) Benham, C.; Mielke, S. DNA mechanics. Annu. Rev. Biomed. Eng. 2005, 7, 21−53. (41) White, J. H. Self-linking and the Gauss integral in higher dimensions. Am. J. Math. 1969, 91, 693−728. (42) White, J. H. An introduction to the geometry and topology of DNA structure. In Mathematical methods for DNA Sequences; Boca Raton, Florida, 1989; pp 225−253. (43) Calugareanu, G. Sur les classes d’isotopie des noeuds tridimensionnels et leurs invariants. Czechoslovak Math. J. 1961, 11, 588−625. (44) Swigon, D.; Coleman, B. D.; Tobias, I. The elastic rod model for DNA and its application to the tertiary structure of DNA minicircles in mononucleosomes. Biophys. J. 1998, 74, 2515−2530. (45) Britton, L.; Olson, W.; Tobias, I. Two perspectives on the twist of DNA. J. Chem. Phys. 2009, 131, 245101. (46) Czapla, L.; Swigon, D.; Olson, W. K. Sequence-Dependent Effects in the Cyclization of Short DNA. J. Chem. Theory Comput. 2006, 2, 685−695. (47) Swigon, D.; Olson, K. W. Mesoscale modeling of multi-proteinDNA assemblies: The role of the catabolic activator protein in Lacrepressor-mediated looping. International Journal of Non-Linear Mechanics 2008, 43, 1082−1093. (48) Nelson, P. C.; Zurla, C.; Brogioli, D.; Beausang, J. F.; Finzi, L.; Dunlap, D. Tethered Particle Motion as a Diagnostic of DNA Tether Length. J. Phys. Chem. B 2006, 110, 17260−17267. (49) Towles, K. B.; Beausang, J. F.; Garcia, H. G.; Phillips, R.; Nelson, P. C. First-principles calculation of DNA looping in tethered particle experiments. Phys. Biol. 2009, 6, 025001. (50) Zoli, M. J-factors of short DNA molecules. J. Chem. Phys. 2016, 144, 214104. (51) Zoli, M. Entropic penalties in circular DNA assembly. J. Chem. Phys. 2014, 141, 174112. (52) Wiggins, P. A.; van der Heijden, T.; Moreno-Herrero, F.; Spakowitz, A.; Phillips, R.; Widom, J.; Dekker, C.; Nelson, P. C. High flexibility of DNA on short length scales probed by atomic force microscopy. Nat. Nanotechnol. 2006, 1, 137−141. (53) Cloutier, T. E.; Widom, J. Spontaneous Sharp Bending of Double-Stranded DNA. Mol. Cell 2004, 14, 355−362. (54) Vologodskii, A.; Du, Q.; Frank-Kamenetskii, M. D. Bending of short DNA helices. Artificial DNA: PNA & XNA 2013, 4, 1−3. (55) Le, T. T.; Kim, H. D. Probing the elastic limit of DNA bending. Nucleic Acids Res. 2014, 42, 10786−10794. (56) Irobalieva, R. N.; Fogg, J. M.; Catanese, D. J., Jr; Sutthibutpong, T.; Chen, M.; Barker, A. K.; Ludtke, S. J.; Harris, S. A.; Schmid, M. F.; Chiu, W.; Zechiedrich, L. Structural diversity of supercoiled DNA. Nat. Commun. 2015, 6, 8440. (57) Perez, P. J.; Olson, W. K. Insights into genome architecture deduced from the properties of short Lac repressor-mediated DNA loops. Biophys. Rev. 2016, 8, 135−144. (58) Wang, Q.; Irobalieva, R. N.; Chiu, W.; Schmid, M. F.; Fogg, J. M.; Zechiedrich, L.; Pettitt, B. M. Influence of DNA sequence on the M


Effects of Protein-Induced Local Bending and Sequence Dependence

Recommend Documents