Allosteric interactions in a birod model of DNA

TL;DR

This paper develops a birod model of DNA to quantitatively analyze allosteric interactions, showing that the interaction energy decays exponentially with distance and oscillates with the DNA's helical periodicity, aligning with experimental observations.

Contribution

The study introduces a birod elasticity model to explain allosteric interactions in DNA, linking mechanical properties to interaction decay length and sequence composition.

Findings

Interaction energy decays exponentially with distance

Energy oscillates with DNA's helical periodicity

Decay length depends on DNA's mechanical properties

Abstract

Allosteric interactions between molecules bound to DNA at distant locations have been known for a long time. The phenomenon has been studied via experiments and numerical simulations, but a comprehensive understanding grounded in a theory of DNA elasticity remains a challenge. Here we quantify allosteric interactions between two entities bound to DNA by using the theory of birods. We recognize that molecules bound to DNA cause local deformations that can be captured in a birod model which consists of two elastic strands interacting via an elastic web representing the base-pairs. We show that the displacement field caused by bound entities decays exponentially with distance from the binding site. We compute the interaction energy between two proteins on DNA as a function of distance between them and find that it decays exponentially while oscillating with the periodicity of the double-helix, in excellent agreement with experiments. The decay length of the interaction energy can be determined in terms of the mechanical properties of the strands and the webbing in our birod model, and it varies with the GC content of the DNA. Our model provides a framework for viewing allosteric interactions in DNA within the ambit of configurational forces of continuum elasticity.

Equations

$$\Delta G=\Delta G_A+\Delta G_B+\Delta \Delta G_{AB}(a),$$

$$\begin{array}{rl}& {\mathbf{r}}^{+}=x\text{leavevmode}{\mathbf{e}}_1+\frac{d}{2}\text{leavevmode}{\mathbf{e}}_2+u\text{leavevmode}{\mathbf{e}}_1+w\text{leavevmode}{\mathbf{e}}_2,\\ & {\mathbf{r}}^{-}=x\text{leavevmode}{\mathbf{e}}_1-\frac{d}{2}\text{leavevmode}{\mathbf{e}}_2+u\text{leavevmode}{\mathbf{e}}_1-w\text{leavevmode}{\mathbf{e}}_2,\end{array}$$

$$\begin{array}{rl}& \text{leavevmode}{\mathbf{d}}_1^{\pm }=\cos \theta \text{leavevmode}{\mathbf{e}}_1\pm \sin \theta \text{leavevmode}{\mathbf{e}}_2\approx \text{leavevmode}{\mathbf{e}}_1\pm \theta \text{leavevmode}{\mathbf{e}}_{,}\\ & \text{leavevmode}{\mathbf{d}}_2^{\pm }=\mp \sin \theta \text{leavevmode}{\mathbf{e}}_1+\cos \theta \text{leavevmode}{\mathbf{e}}_2\approx \pm \theta \text{leavevmode}{\mathbf{e}}_1+\text{leavevmode}{\mathbf{e}}_2,\\ & {\mathbf{R}}^{\pm }=\left[\begin{array}{ccc}\cos \theta & \mp \sin \theta & 0\\ \pm \sin \theta & \cos \theta & 0\\ 0& 0& 1\end{array}\right]\approx \left[\begin{array}{ccc}1& \mp \theta & 0\\ \pm \theta & 1& 0\\ 0& 0& 1\end{array}\right].\end{array}$$

$$\begin{array}{r}\mathbf{R}=({\mathbf{R}}^{+}{\mathbf{R}}^{-T}{)}^{1/2}{\mathbf{R}}^{-}={\mathbf{I}}_{3\times 3}.\end{array}$$

$$\begin{array}{rl}& \mathbf{P}=({\mathbf{R}}^{+}{\mathbf{R}}^{-T}{)}^{1/2}=\left[\begin{array}{ccc}\cos \theta & -\sin \theta & 0\\ \sin \theta & \cos \theta & 0\\ 0& 0& 1\end{array}\right]\approx \left[\begin{array}{ccc}1& -\theta & 0\\ \theta & 1& 0\\ 0& 0& 1\end{array}\right].\end{array}$$

$$\begin{array}{rl}& {\mathbf{v}}_0^{\pm }=\frac{\partial {\mathbf{r}}_0^{\pm }}{\partial x}=\text{leavevmode}{\mathbf{e}}_1,\\ & {\mathbf{v}}^{\pm }=\frac{\partial {\mathbf{r}}^{\pm }}{\partial x}=(1+u_x)\text{leavevmode}{\mathbf{e}}_1\pm w_x\text{leavevmode}{\mathbf{e}}_2.\end{array}$$

$$\begin{array}{r}{\mathbf{n}}^{\pm }=EA{u}_x\text{leavevmode}{\mathbf{e}}_1\pm GA(w_x-\theta )\text{leavevmode}{\mathbf{e}}_2.\end{array}$$

$$\begin{array}{rl}& {\mathbf{U}}^{\pm }=\frac{\partial {\mathbf{R}}^{\pm }}{\partial x}{\mathbf{R}}^{\pm T}=\left[\begin{array}{ccc}0& \mp {\theta }_x& 0\\ \pm {\theta }_x& 0& 0\\ 0& 0& 0\end{array}\right],\\ & {\mathbf{\kappa }}^{\pm }=\pm {\theta }_x\text{leavevmode}{\mathbf{e}}_3.\end{array}$$

$$\begin{array}{r}{\mathbf{m}}^{\pm }=\pm EI{\theta }_x\text{leavevmode}{\mathbf{e}}_3.\end{array}$$

$$\begin{array}{rl}2\mathbf{f}=& \mathbf{RH}{\mathbf{R}}^T[\mathbf{w}-\mathbf{R}\check{{\mathbf{w}}_0}]\approx Lw\text{leavevmode}{\mathbf{e}}_2,\end{array}$$

$$\begin{array}{rl}2\mathbf{c}& =\frac{1}{\alpha }\mathbf{RG}{\mathbf{R}}^T(\mathbf{\eta }-\mathbf{R}\check{\text{upeta}})-\mathbf{\eta }\times (\mathbf{w}\times \mathbf{f})\approx K\theta \text{leavevmode}{\mathbf{e}}_3,\end{array}$$

$$\begin{array}{rl}& {\mathbf{n}}_x=0,\\ & {\mathbf{m}}_x+{\mathbf{r}}_x\times \mathbf{n}=0.\end{array}$$

$$\begin{array}{rl}& {\mathbf{n}}_x^c-2\mathbf{f}=0,\\ & {\mathbf{m}}_x^c+{\mathbf{r}}_x\times {\mathbf{n}}^c-\mathbf{c}=0.\end{array}$$

$$\begin{array}{rl}& EA{u}_{xx}=0,\\ & 2GA(w_{xx}-{\theta }_x)-Lw=0,\\ & 2(EI{\theta }_{xx}-d/2EA{u}_{xx})+(1+u_x)\text{leavevmode}{\mathbf{e}}_1\times 2GA(w_x-\theta )\text{leavevmode}{\mathbf{e}}_2-K\theta \text{leavevmode}{\mathbf{e}}_3=0.\end{array}$$

$$EI{w}_{xxxx}-(\frac{EIL}{2GA}+\frac{K}{2})w_{xx}+(\frac{L}{2}+\frac{KL}{4GA})w=0.$$

$$\begin{array}{r}EI{w}_{xxxx}-\frac{K}{2}{w}_{xx}+\frac{L}{2}w=0.\end{array}$$

$$\begin{array}{rl}& w(x)=A{e}^{\lambda x}+B{e}^{\mu x}\text{when}x<0,\\ & w(x)=A{e}^{-\lambda x}+B{e}^{-\mu x}\text{when}x>0,\end{array}$$

$$\begin{array}{r}E[w]=EI{w}_{xx}^2+\frac{1}{2}K{w}_x^2+\frac{1}{2}L{w}^2.\end{array}$$

$$\begin{split}\Delta G=&\frac{L}{2}\Big{(}\frac{e^{-\lambda a}\left(A^{2}\lambda^{2}\mu-A^{2}\mu^{3}+A^{2}\lambda^{3}\mu a-A^{2}\lambda\mu^{3}a-4AB\lambda\mu^{2}\right)}{\lambda\mu(\lambda^{2}-\mu^{2})}+\\ &\quad\quad\frac{e^{-\mu a}\left(4AB\lambda^{2}\mu+B^{2}\lambda^{3}-B^{2}\lambda\mu^{2}+B^{2}\lambda^{3}\mu a-B^{2}\lambda\mu^{3}a\right)}{\lambda\mu(\lambda^{2}-\mu^{2})}\Big{)}+\\ &\frac{K}{2}\Big{(}\frac{e^{-\lambda a}\left(A^{2}\lambda^{3}-A^{2}\lambda\mu^{2}-A^{2}\lambda^{4}a+A^{2}\lambda^{2}\mu^{2}a+4AB\lambda^{2}\mu\right)}{(\lambda^{2}-\mu^{2})}+\\ &\quad\quad\frac{e^{-\mu a}\left(-4AB\lambda\mu^{2}+B^{2}\lambda^{2}\mu-B^{2}\mu^{3}-B^{2}\lambda^{2}\mu^{2}a+B^{2}\mu^{4}a\right)}{(\lambda^{2}-\mu^{2})}\Big{)}+\\ &EI\Big{(}\frac{e^{-\lambda a}\left(A^{2}\lambda^{5}-A^{2}\lambda^{3}\mu^{2}+A^{2}\lambda^{6}a-A^{2}\lambda^{4}\mu^{2}a-4AB\lambda^{2}\mu^{3}\right)}{(\lambda^{2}-\mu^{2})}+\\ &\quad\quad\frac{e^{-\mu a}\left(4AB\lambda^{3}\mu^{2}+B^{2}\lambda^{2}\mu^{3}-B^{2}\mu^{5}+B^{2}\lambda^{2}\mu^{4}a-B^{2}\mu^{6}a\right)}{(\lambda^{2}-\mu^{2})}\Big{)}.\end{split}$$

$$\mathbf{r}^{+}_{0}=b\Big{(}\cos\omega x\leavevmode\nobreak\ \mathbf{e}_{1}+\sin\omega x\leavevmode\nobreak\ \mathbf{e}_{2}\Big{)}+x\leavevmode\nobreak\ \mathbf{e}_{3},$$

$$\mathbf{r}^{-}_{0}=b\Big{(}\cos(\omega x+\alpha)\leavevmode\nobreak\ \mathbf{e}_{1}+\sin(\omega x+\alpha)\leavevmode\nobreak\ \mathbf{e}_{2}\Big{)}+x\leavevmode\nobreak\ \mathbf{e}_{3}.$$

$$\begin{split}&\mathbf{r}^{+}_{0}=b\Big{(}\cos\omega x\leavevmode\nobreak\ \mathbf{e}_{1}+\sin\omega x\leavevmode\nobreak\ \mathbf{e}_{2}\Big{)}+x\leavevmode\nobreak\ \mathbf{e}_{3},\\ &\mathbf{r}^{+}(x)=(b+r)\Big{(}\cos(\omega x+\beta^{+})\leavevmode\nobreak\ \mathbf{d}_{1}+\sin(\omega x+\beta^{+})\leavevmode\nobreak\ \mathbf{d}_{2}\Big{)}+\int_{0}^{x}dx(1+b\xi)\leavevmode\nobreak\ \mathbf{d}_{3}.\\ \end{split}$$

$$r=r(x),{\beta }^{+}={\beta }^{+}(x),\xi =\xi (x).$$

$$\begin{array}{r}\text{leavevmode}{\mathbf{d}}_i=\mathbf{Z}\text{leavevmode}{\mathbf{e}}_i,\mathbf{Z}=\sum _{i=1}^3\text{leavevmode}{\mathbf{d}}_i\otimes \text{leavevmode}{\mathbf{e}}_i,i=1,2,3\end{array}$$

$$\begin{array}{rl}& \text{leavevmode}{\mathbf{d}}_{ix}=\text{upkappa}\times \text{leavevmode}{\mathbf{d}}_i,\text{where}\text{upkappa}=k_1\text{leavevmode}{\mathbf{d}}_1+k_2\text{leavevmode}{\mathbf{d}}_2+k_3\text{leavevmode}{\mathbf{d}}_3.\\ & \text{leavevmode}{\mathbf{d}}_{1x}=k_3\text{leavevmode}{\mathbf{d}}_2-k_2\text{leavevmode}{\mathbf{d}}_3,\text{leavevmode}{\mathbf{d}}_{2x}=k_1\text{leavevmode}{\mathbf{d}}_3-k_3\text{leavevmode}{\mathbf{d}}_2,\text{leavevmode}{\mathbf{d}}_{3x}=k_2\text{leavevmode}{\mathbf{d}}_1-k_1\text{leavevmode}{\mathbf{d}}_2.\end{array}$$

$$r^{+}(x),k_1(x),k_2(x),k_3(x),\zeta (x),{\beta }^{+}(x)\sim O(\epsilon ).$$

$$\begin{array}{rl}& {\mathbf{R}}_0^{+}=[{\mathbf{n}}_0^{+}{\mathbf{b}}_0^{+}{\mathbf{t}}^{+}],\\ & {\mathbf{n}}_0^{+}=-\cos \omega x\text{leavevmode}{\mathbf{e}}_1-\sin \omega x\text{leavevmode}{\mathbf{e}}_2,\\ & {\mathbf{b}}_0^{+}=-\cos k(-\sin \omega x\text{leavevmode}{\mathbf{e}}_1+\cos \omega x\text{leavevmode}{\mathbf{e}}_2)+\sin k\text{leavevmode}{\mathbf{e}}_3,\\ & {\mathbf{t}}_0^{+}=\sin k(-\sin \omega x\text{leavevmode}{\mathbf{e}}_1+\cos \omega x\text{leavevmode}{\mathbf{e}}_2)+\cos k\text{leavevmode}{\mathbf{e}}_3,\end{array}$$

$$(\cos \omega x\text{leavevmode}{\mathbf{d}}_1+\sin \omega x\text{leavevmode}{\mathbf{d}}_2)=\text{leavevmode}{\mathbf{f}}_1^{+},(-\sin \omega x\text{leavevmode}{\mathbf{d}}_1+\cos \omega x\text{leavevmode}{\mathbf{d}}_2)=\text{leavevmode}{\mathbf{f}}_2^{+},\text{leavevmode}{\mathbf{d}}_3=\text{leavevmode}{\mathbf{f}}_3^{+}.$$

$$\begin{array}{rl}{\mathbf{r}}_x^{+}& =(r_x-b\omega {\beta }^{+})(\cos \omega x\text{leavevmode}{\mathbf{d}}_1+\sin \omega x\text{leavevmode}{\mathbf{d}}_2)+(b\omega +\omega r+b{\beta }_x^{+}+b{k}_3)(-\sin \omega x\text{leavevmode}{\mathbf{d}}_1+\cos \omega x\text{leavevmode}{\mathbf{d}}_2)+\\ & (1+b\xi -b{k}_2\cos \omega x+b{k}_1\sin \omega x)\text{leavevmode}{\mathbf{d}}_3.\end{array}$$

$$|{\mathbf{r}}_x^{+}{|}^2=1+{\omega }^2{b}^2+2b({\omega }^{2}r+b\omega {\beta }_x^{+}+b\omega k_3+\xi -k_2\cos \omega x+k_1\sin \omega x)+O({\epsilon }^2)=|{\mathbf{r}}_{0x}^{+}{|}^2=1+b^2{\omega }^2,$$

$$\xi -k_2\cos \omega x+k_1\sin \omega x=-{\omega }^{2}r-b\omega (k_3+{\beta }_x^{+}).$$

Full text

1 Introduction

Configurational forces that describe the interaction between defects in an elastic solid are those that depend explicitly on the positions of the defects [29, 38]. For example, two parellel screw dislocations at a distance $a$ from each other interact with a configurational force per unit length proportional to $1/a$ or an energy per unit length proportional to $\log a$ [8]. Similarly, the interaction energy of a point defect located at distance $a$ from an edge dislocation varies as $1/a$. Just as defects produce local elastic fields in a solid, proteins binding to DNA also deform it locally. Since DNA behaves like an elastic rod at scales of a few tens of nanometers [25], we expect that if two proteins bind to DNA separated by a distance $a$ then the deformation fields created by them will overlap and lead to an interaction energy which depends on $a$ in a clearly quantifiable way. This problem has not been theoretically addressed so far, but there is experimental evidence of the interaction. Some of this experimental evidence has been extracted by connecting the interaction energy with the kinetics of protein binding/unbinding. In spirit, this is similar to continuum elasticity in which configurational forces often determine defect dynamics through a kinetic law [29, 38]. Kim et al. [4] have exploited this connection of interaction energies to kinetics to show that gene expression, which depends on RNA polymerase binding affinity to DNA in live bacteria, is a function of the proximity of LacR and T7 RNA polymerase bound to DNA. Similarly, the IHF protein affects RNA polymerase activity in E. Coli DNA [23]. Again, the binding of a drug distamycin to calf thymus DNA has been shown to be cooperative i.e., if one drug molecule binds to the DNA, then it becomes energetically favorable for other drug molecules to bind [24]. Similarly, binding of the Hox transcription factor to DNA contributes nearly 1.5 kcal/mol to binding the Exd transcription factor [26]. These effects are called allosteric interactions on DNA. Our goal in this paper is to quantify interaction energies between proteins binding to DNA as a function of the distance $a$ separating them and the boundary conditions imposed by the proteins on the DNA. We will apply our methods to the quantitative experimental results of Kim et al.[4] who measure allosteric effects on gene expression as well as transcription factor affinity to DNA.

In the experiments of Kim et al. [4] one end of a DNA molecule is attached to the passivated surface of a flow cell and binding sites are provided for two specific proteins to bind. The length of the DNA between these binding sites, $a$, is increased in 1bp increments between 7 base-pairs (bp) and 45bp. First, one type of fluorescently labeled protein (call it A) is flowed into the cell so that it binds to the DNA. Then, the second protein (call it B) is flowed in at a specific concentration. The dissociation times of the fluorescent protein are then monitored as a function of $a$. This dissociation time depends on the free-energy change $\Delta G$ of the DNA + two protein complex from the state when the two proteins are bound to that when protein A is unbound. Now, in general, the free energy $\Delta G$ of the ternary complex formed by the DNA and proteins A and B consists of three parts [4]:

[TABLE]

where $\Delta G_{A}$ and $\Delta G_{B}$ are the free energy changes caused by binding of A and B respectively to the DNA. These are constants. The last term $\Delta\Delta G_{AB}(a)$ is the portion of the free energy change that accounts for the interaction of the two proteins bound to the DNA while being separated by a distance $a$. The off-rate of A, which is affected by this term, is plotted as a function of $a$ in Kim et al. [4] and it is found that it oscillates with a period of 10-11bp with the amplitude of oscillation decreasing as a function of $a$. Similar curves for a free energy as a function of separation between protein binding sites on DNA have been obtained experimentally for the binding of the lac repressor to DNA[12, 13]. It has been shown that these free energy profiles can be reproduced by modeling DNA as an elastic rod which is forced into forming a loop due to stereo-specific binding of the ${\it lac}$ repressor monomers [9] which come together due to thermal fluctuations. However, Kim et al. [4] have ruled out DNA loop formation by careful experimental design and choice of DNA binding proteins. They have also found that the form of the curve is independent of ionic strength (ruling out electrostatic interactions between A and B), but dependent on modifications of the linker DNA. Kim et al. [4] infer that this implies $\Delta\Delta G_{AB}(a)$ largely depends on DNA mechanical properties. However, as yet there is no analytical description of how the interaction energy $\Delta\Delta G_{AB}(a)$ depends on the DNA mechanical properties.

Allosteric effects and their relation to protein DNA interactions have been studied using molecular dynamic (MD) simulations [21, 22]. Gu et al.[21] have studied various kinds of deformations which include shift, roll, rise, twist, slide, and tilt of the DNA bases. They observed a sinusoidal correlation in the major groove widths similar to the one observed by Kim et al.[4]. Furthermore, Gu et al. point out that the presence of GC rich sequences dampens the allosteric effects which is what Kim et al. observe experimentally. Major groove widths have also been implicated in the MD simulations of Hancock et al. [22] who show how bound proteins alter this quantity. In contrast, our approach in this paper is based on elastic energy considerations and could compliment the analysis of major groove widths as an indicator of allostery in DNA.

Our goal in this paper is to quantitatively describe allosteric interactions using the birod model of DNA of Moakher and Maddocks [1] who originally derived it to study DNA melting. This birod model is a double-stranded rod theory in which in addition to the standard variables of a Cosserat rod theory (i.e., center line of the rod cross-section $\mathbf{r}(s,t)$ and a material frame $[\mathbf{d}_{1}(s,t)\quad\mathbf{d}_{2}(s,t)\quad\mathbf{d}_{3}(s,t)]$), there are two micro-structural variables – $\mathbf{w}(s,t)$, a micro-displacement measuring the change in distance between the two strands, and $\mathbf{P}(s,t)$, a micro-rotation measuring the change in orientation of one strand relative to the other. Fortunately, the forces conjugate to these micro-structural variables obey balance laws that look similar to the balance of forces and moments equations of a standard Cosserat rod. They are coupled to the macroscopic balance equations for the center-line of the rod through distributed body forces and moments. Moakher and Maddocks [1] have provided hyper-elastic constitutive laws for these micro-structural variables that are based on quadratic energies.

The theory of birods has been successfully used by Lessinnes et al. [5, 7] to study the growth and evolution of two filaments elastically bound to each other. These authors have demonstrated the utility of the theory in accurately modelling biological structures across multiple length scales –from tissues and arteries, to growth of roots and stems in plants. Manning et al. [3] have used both a discrete base-pair model and a corresponding continuum rod model to study the cyclization of short DNA molecules (150 bp). The results obtained from both these approaches match remarkably well. However, for shorter length scales ($\sim 16$ bp) Lankas et al. [2] have assayed the merit of the assumptions of rigid bases versus rigid base-pairs to estimate the stiffness parameters of a DNA oligomer and found that the simulated data is closely consistent with the assumption of rigid bases, but not rigid basepairs. This inevitably necessitates the inclusion of elasticity of base-pairs via the webbing in an elastic birod [1] to accurately model the local deformations caused by proteins at small length scales.

2 Strategy to compute interaction energy

In this section, we give a concise blueprint of our strategy to approach the problem of calculating the interaction energy for two proteins binding to DNA. We model the DNA as a helical birod [1] which has two elastic components – outer strands and the connecting web. The two sugar-phosphate backbones of DNA correspond to the outer elastic strands which interact by means of complimentary base-pairing represented by the elastic web in our case. We give a stepwise procedure to do the calculation and in the following sections we label each step.

We begin by assuming a form of displacement for each of the outer strands which are assumed to be inextensible and unshearable. 2. 2.

We then use this displacement to calculate the tangent, normal and binormal to the deformed configuration of the outer strands thereby obtaining the rotation matrix attached to the deformed configuration of the outer strands. 3. 3.

Once we get the deformation and rotation of the outer strands, we use these to calculate the extension, shear and rotation of the web. 4. 4.

At this point, we are in a position to substitute these quantities into the balance laws for the birod. We, then, seek non-zero solutions to the resulting system of differential equations. This leads to an eigenvalue problem. 5. 5.

In the next step, we apply the boundary conditions to evaluate the constants. 6. 6.

We carry out this process first when there is a single protein binding onto the DNA, and second when there are two proteins binding. 7. 7.

Finally, we subtract the two energies obtained in the previous step to get our energy of interaction. We find that it takes the form of a decaying exponential oscillating with the periodicity of the underlying DNA helix.

3 Exponential decay of interaction energy in a ‘ladder’

The calculation described above is considerably involved, so we first illustrate the main concepts in a simpler birod model which we call a ‘ladder’ because it is not helical. We mimic the binding of a protein by force pairs that tend to widen the ladder as shown in fig. 1. Our goal in this section is to demonstrate the utility of the apparatus in section 2 by computing the interaction energy for two force pairs separated by a distance $a$ as shown in fig. 1. We work with a planar 2D birod in this section and assume small elastic deformations in the outer strands and web to keep the calculations tractable. We, ultimately, find that the interaction energy between the force pairs decays exponentially with distance $a$.

3.1 Step 1: Kinematic description of the two strands

We use the arclength parameter $x$ to describe the mechanics of the birod. In the reference configuration, both the strands $\pm$ are straight, $\mathbf{r}_{0}^{\pm}=x\leavevmode\nobreak\ \mathbf{e}_{1}\pm\frac{d}{2}\leavevmode\nobreak\ \mathbf{e}_{2}$, separated by distance $d$. Here $\leavevmode\nobreak\ \mathbf{e}_{1}$ is a unit vector along the length of the birod, $\leavevmode\nobreak\ \mathbf{e}_{2}$ is a unit vector perpendicular to each birod bridging the gap between them and $\leavevmode\nobreak\ \mathbf{e}_{3}$ is normal to the plane of the birod as shown in fig. 1. We begin by assuming a general displacement in $\leavevmode\nobreak\ \mathbf{e}_{1}-\leavevmode\nobreak\ \mathbf{e}_{2}$ plane. For the geometry shown in fig. 1 we expect a mirror symmetry for deformation profiles along $\leavevmode\nobreak\ \mathbf{e}_{1}$ such that

[TABLE]

where $u=u(x)$ and $w=w(x)$ are displacements along the $\leavevmode\nobreak\ \mathbf{e}_{1}$ and $\leavevmode\nobreak\ \mathbf{e}_{2}$ directions, respectively.

3.2 Step 2: Rotation of the two strands

At each point $x$ on the $\pm$ strands we attach an orthogonal rotation frame which is simply $\mathbf{R}_{0}^{\pm}=[\leavevmode\nobreak\ \mathbf{e}_{1}\quad\leavevmode\nobreak\ \mathbf{e}_{2}\quad\leavevmode\nobreak\ \mathbf{e}_{3}]=\mathbf{1}_{3\times 3}$ (the identity matrix) in the reference configuration. The vectors $\leavevmode\nobreak\ \mathbf{e}_{1}$ and $\leavevmode\nobreak\ \mathbf{e}_{2}$ map onto $\leavevmode\nobreak\ \mathbf{d}_{1,2}^{+}$ and $\leavevmode\nobreak\ \mathbf{d}_{1,2}^{-}$ in the deformed configuration for the positive and negative strand, respectively. The $\leavevmode\nobreak\ \mathbf{d}_{i}$, $i=1,2,3$ are again unit vectors.

[TABLE]

We assume small $\theta$ to keep the calculations tractable.

3.3 Step 3: Extension and rotation of the web

We decompose the kinematics of the web into a macroscopic deformation and a microscopic deformation [1]. The former describes the rigid displacement and rotation, while the latter is related to the force and moment transferred by the web. The macro- displacement vector $\mathbf{r}$ is defined as $\mathbf{r}=\frac{\mathbf{r}^{+}+\mathbf{r}^{-}}{2}=x\leavevmode\nobreak\ \mathbf{e}_{1}+u\leavevmode\nobreak\ \mathbf{e}_{1}$ [1]. The macro- rotation tensor is $\mathbf{R}$ defined as $\mathbf{R}=(\mathbf{R}^{+}\mathbf{R}^{-T})^{1/2}\mathbf{R}^{-}$ [1], which in our case is

[TABLE]

We define another tensor $\mathbf{P}$ relating $\mathbf{R}^{+}$ and $\mathbf{R}^{-}$ to $\mathbf{R}$. An elastic constitutive relation discussed in further sections connects the micro- rotation tensor $\mathbf{P}=(\mathbf{R}^{+}\mathbf{R}^{-T})^{1/2}$ to the moment transferred by the web.

[TABLE]

We need to calculate the Gibbs rotation vector $\boldsymbol{\eta}=\tan\frac{\lambda}{2}\hat{\mathbf{k}}$, where $\lambda$ is obtained from $1+2\cos\lambda=\text{tr}(\mathbf{P})$ and $\hat{\mathbf{k}}$ is the eigenvector of $\mathbf{P}$ i.e. $\mathbf{P}\hat{\mathbf{k}}=\hat{\mathbf{k}}$. We need $\boldsymbol{\eta}$ in the subsequent section to compute the moment transferred by the web [1]. By direct observations, $\lambda=\theta$ and $\hat{\mathbf{k}}=\leavevmode\nobreak\ \mathbf{e}_{3}$, so that $\boldsymbol{\eta}=\tan\frac{\theta}{2}\leavevmode\nobreak\ \mathbf{e}_{3}$. The Gibbs rotation vector in the reference configuration $\boldsymbol{\eta}_{0}=0$.

The micro- displacement of the web is defined by $\boldsymbol{w}=\frac{\mathbf{r}^{+}-\mathbf{r}^{-}}{2}$, which is $\boldsymbol{w}_{0}=\frac{d}{2}\leavevmode\nobreak\ \mathbf{e}_{2}$ in the reference configuration and $\boldsymbol{w}=(\frac{d}{2}+w)\leavevmode\nobreak\ \mathbf{e}_{2}$ in the current configuration. We need $\boldsymbol{w}$ and $\boldsymbol{w}_{0}$ to compute the force transferred by the web.

3.4 Step 4: Governing differential equations

We calculate various strains and curvatures associated with the deformation and relate them to the contact force and moment, respectively, which go into the governing equations. For detailed discussion on the relations used in this section we refer the reader to Moakher and Maddocks [1]. The governing equations of the birod consist of three kinetic components: the contact forces in the two strands $\boldsymbol{n}^{\pm}$, the contact moments $\boldsymbol{m}^{\pm}$, and the force $\boldsymbol{f}$ and moment $\boldsymbol{c}$ transferred by the $-$ strand onto the $+$ strand. We compute each of these components as follows:

$\boldsymbol{n}^{\pm}$: We need strains in the current configuration $\boldsymbol{v}^{\pm}$ and in the reference configuration $\boldsymbol{v}^{\pm}_{0}$, in the strands to compute $\boldsymbol{n}^{\pm}$. These strains are:

[TABLE]

The contact forces $\boldsymbol{n}^{\pm}=\mathbf{R}^{\pm}\mathbf{C}\mathbf{R}^{\pm T}\boldsymbol{v}^{\pm}$ where $\mathbf{C}$ is a second order tensor such that $\mathbf{C}_{11}=EA$, $\mathbf{C}_{22}=GA$ and $\mathbf{C}_{12}=\mathbf{C}_{21}=0$. Here $E$ is the stretch modulus, $G$ shear modulus and $A$ is the cross-sectional area of the strands. Upon performing the calculation and taking account of the fact that $u,w$ and $\theta$ are small and upon ignoring higher order terms we get,

[TABLE] 2. 2.

$\boldsymbol{m}^{\pm}$: For calculating the contact moments $\boldsymbol{m}^{\pm}$ in the respective strands we need the curvature vector $\boldsymbol{\kappa}^{\pm}$ for the two strands, which can, in turn, be obtained by computing the axial vector of the skew-symmetric matrices $\mathbf{U}^{\pm}=\frac{\partial{\mathbf{R}^{\pm}}}{\partial{x}}\mathbf{R}^{\pm T}$.

[TABLE]

The contact moment $\boldsymbol{m}^{\pm}$ is related to the curvature via a bending rigidity $EI$ such that

[TABLE]

Here, $I$ is the moment of inertia of the cross-section of the outer strands. 3. 3.

$\boldsymbol{f}$ and $\boldsymbol{c}$: The force transferred by the web $\boldsymbol{f}$ is proportional to the change in the dimensions of the web quantified by $\boldsymbol{w}$ and $\boldsymbol{w}_{0}$ in the previous sections such that,

[TABLE]

where $\mathbf{H}$ is a diagonal second order elasticity tensor such that $\mathbf{H}_{22}=L$. Similarly, the moment transferred by the web $\boldsymbol{c}$ is elastically related to $\boldsymbol{\eta}$ and $\boldsymbol{\eta}_{0}$ calculated in the previous sections:

[TABLE]

where $\alpha=\frac{2}{1+||\boldsymbol{\eta}||^{2}}$, $\mathbf{G}$ is a second order diagonal elasticity tensor and $K=\frac{G_{33}}{2}$.

The governing equations from box 4 in [1] are given by,

[TABLE]

In the above equations, $\boldsymbol{n}=\boldsymbol{n}^{+}+\boldsymbol{n}^{-}=2EAu_{x}\leavevmode\nobreak\ \mathbf{e}_{1}$, $\boldsymbol{n}^{c}=\boldsymbol{n}^{+}-\boldsymbol{n}^{-}=2GA(w_{x}-\theta)\leavevmode\nobreak\ \mathbf{e}_{2}$, $\boldsymbol{m}=\boldsymbol{m}^{+}+\boldsymbol{m}^{-}+\boldsymbol{w}\times\boldsymbol{n}^{c}=0$ and $\boldsymbol{m}^{c}=\boldsymbol{m}^{+}-\boldsymbol{m}^{-}+\boldsymbol{w}\times\boldsymbol{n}=2EI\theta_{x}\leavevmode\nobreak\ \mathbf{e}_{3}+(d/2+w)\leavevmode\nobreak\ \mathbf{e}_{2}\times 2EAu_{x}\leavevmode\nobreak\ \mathbf{e}_{1}\approx 2[EI\theta_{x}-\frac{d}{2}EAu_{x}]\leavevmode\nobreak\ \mathbf{e}_{3}$. Upon substituting these values into the governing equations we get,

[TABLE]

We use $\theta_{x}=w_{xx}-\frac{L}{2GA}w$ and $u_{xx}=0$ and get,

[TABLE]

If we further assume that the outer strands are unshearable ($GA\to\infty$ and $\theta=w_{x}$), the above equation reduces to a simpler equation.

[TABLE]

3.5 Step 5,6 and 7: Interaction Energy

We substitute $w=e^{ms}$, and get eigenvalues $m=\pm\lambda,\pm\mu$. For illustration purposes, we assume $\lambda$ and $\mu$ are real numbers (i.e., $K^{2}-32L>0$) and the ladder extends from $-\infty$ in the negative $\leavevmode\nobreak\ \mathbf{e}_{1}$ direction to $+\infty$ in the positive $\leavevmode\nobreak\ \mathbf{e}_{1}$ direction with $w=w_{x}=0$ at $x=\pm\infty$. Hence, for a force pair at $x=0$

[TABLE]

for some constants $A$ and $B$ which could be determined using boundary conditions in step 5. For two force pairs separated by a distance $a$, the displacement profile $w_{2}(x)=w(x)+w(x-a)$. The elastic energy in the deformed configuration is computed in step 6 and is given by,

[TABLE]

Finally, we compute the interaction energy defined by $\Delta G=E[w_{2}]-2E[w]$ in step 7 and find that it decreases exponentially with the distance $a$.

[TABLE]

In the next section we will follow these steps for a helical birod model of DNA.

4 Interaction energy for two DNA binding proteins

4.1 Assumptions

We first outline the various assumptions and set up the underlying framework for our helical birod model of DNA. We assume elastic deformations throughout. When a protein binds to DNA it causes local bending and twisting. We assume that the resulting twist and curvatures are small. These curvatures could possibly add up to produce large displacements and rotations. The two phosphate backbones of DNA constitute the helical outer strands which are out of phase by a phase angle $\alpha=2.1$ radians. We assume these backbones to be inextensible. These outer strands consist of sugar phosphate single bonds. Thus, we assume that they can not support twisting moments. The inextensibility of the outer strands is a strong geometrical constraint which induces a change in the radius and phase angle between the two helices when a protein causes local deformations. We assume that these changes are small and of the same order as the curvatures.

4.2 Step 1: Deformation of the outer strands

DNA consists of two helical strands with radius $b=1$ nm and pitch $p=3.4$ nm, out of phase by $\alpha=2.1$ radians, wrapped around a common axis as shown in fig. 2. We follow the notation used by Moakher and Maddocks [1] and refer to the two stands as $\pm$. The undeformed state of the outer strands denoted by $\mathbf{r}_{0}^{\pm}(x)$ is a helix with a constant radius and pitch. We choose to parametrize both the curves by arclength parameter $x$. Here, $\omega=\frac{2\pi}{p}$ and $k$ is the characteristic angle of the helix such that $\tan k=\frac{2\pi b}{p}=\omega b$.

[TABLE]

[TABLE]

Let us now focus on the two strands separately. The calculations for the $+$ strand are given in this section while the results for the $-$ strand are given in the appendix. We posit a form of displacement wherein the radius of the helix changes and its axis is allowed to take arbitrary shapes within the ambit of the assumptions specified in section 2. Here $[\leavevmode\nobreak\ \mathbf{e}_{1},\leavevmode\nobreak\ \mathbf{e}_{2},\leavevmode\nobreak\ \mathbf{e}_{3}]$ denotes the standard spatial reference frame and $\leavevmode\nobreak\ \mathbf{e}_{3}$ is along the common axis of the two helices $\pm$ in the reference configuration. This common axis in the deformed configuration is defined by the set of orthogonal directors $[\leavevmode\nobreak\ \mathbf{d}_{1}(x),\leavevmode\nobreak\ \mathbf{d}_{2}(x),\leavevmode\nobreak\ \mathbf{d}_{3}(x)]$. The displacement fields which define the undeformed and deformed configuration are,

[TABLE]

where

[TABLE]

Here $r$ is the change in the radius of the helix, $\beta^{+}$ is the change in the phase of the $+$ strand, and $\xi$ can be considered as a stretching of the axis of the helix.

Let $\mathbf{Z}$ be a second order orthogonal tensor which relates the directors of the deformed centerline $\leavevmode\nobreak\ \mathbf{d}_{i}$ to those of the undeformed one $\leavevmode\nobreak\ \mathbf{e}_{i}$, $i=1,2,3$. As stated in section 2, the curvatures $(k_{1},k_{2},k_{3})$ associated with the deformation of the centerline are assumed to be small, nonetheless these could aggregate to potentially produce large rotations. The orthogonal tensor $\mathbf{Z}$ operates as follows.

[TABLE]

and

[TABLE]

In the above equations, we assume that

[TABLE]

Thus, in the treatment henceforth, any product terms such as $r^{2}$ or $\xi k_{3}$ are $O(\varepsilon^{2})$ and are neglected.

4.3 Step 2: Rotation of strands

We proceed in a standard way by attaching a Frenet-Serret director frame consisting of normal, binormal and tangent to each cross-section of the strand as shown in fig. 2. We denote it by $\mathbf{R}^{+}_{0}(x)$ in the reference configuration.

[TABLE]

For the sake of brevity, we use

[TABLE]

As the strand deforms, the frame $\mathbf{R}^{+}_{0}$ evolves into $\mathbf{R}^{+}(x)$ which consists of normal, binormal and tangent to the deformed configuration of the strand. Our next step is to calculate the tangent vector to the deformed configuration. We differentiate eqn. (19) to obtain,

[TABLE]

We assume the strand to be inextensible and unshearable. This means,

[TABLE]

which leads us to the inextensibility condition:

[TABLE]

We will subsequently use this equation to impose boundary conditions. We substitute eqn. (24) into into eqn. (23) to get,

[TABLE]

Now, we need to find the director frame for the strand in the deformed configuration. We start by calculating the tangent vector,

[TABLE]

We differentiate the tangent vector to calculate the normal in the deformed configuration

[TABLE]

We can use the above expression to calculate the curvature $\Omega^{+}$ for the strand. We find that this is equal to the sum of the original curvature ($\omega\sin k$) and the one induced by the process of deformation $\kappa^{+}$. Hence,

[TABLE]

The bending moment $\boldsymbol{m}^{+}$ in the strand is proportional to $\kappa^{+}$.

[TABLE]

Also, the normal is

[TABLE]

where

[TABLE]

Using the above deformed orthogonal frame attached to each cross section

[TABLE]

where $\mathbf{\Theta^{+}}$ is a skew symmetric tensor and $\mathbf{Z}=\sum_{i=1}^{3}\leavevmode\nobreak\ \mathbf{d}_{i}\times\leavevmode\nobreak\ \mathbf{e}_{i}$ as defined in eqn. LABEL:def_Z,

[TABLE]

We can derive all the above quantities $\mathbf{r}^{-},\mathbf{R}^{-}$ and $\kappa^{-}$ etc., for the $-$ strand too. We give the relevant expressions for these quantities in the appendix.

4.4 Step 3: Mechanics of base-pairing

The sugar-phosphate backbones of the DNA molecule are tied together by means of complimentary base-pairing. We model the base-pairing by elastic rods capable of extension, shear, bending and twisting. We attach the orthogonal frame $\mathbf{Q}_{0}=[\leavevmode\nobreak\ \mathbf{f}_{01}\quad\leavevmode\nobreak\ \mathbf{f}_{02}\quad\leavevmode\nobreak\ \mathbf{f}_{03}]$ to the strands such that $\leavevmode\nobreak\ \mathbf{f}_{01}$ is a unit vector pointing from the $-$ strand to the $+$ strand in the reference configuration as shown in fig. 2. Thus,

[TABLE]

We denote the two ends of the rod in the web as $\pm$ such that the $+$ end lies on the $+$ strand and the $-$ end lies on the $-$ strand. The deformation of the web is completely determined by the displacement $(\mathbf{r}^{+}(x),\mathbf{r}^{-}(x))$ and rotation $(\mathbf{R}^{+}(x),\mathbf{R}^{-}(x))$ of its ends. As the outer strands undergo the deformation prescribed by eqn. (19), the strands themselves undergo various kinds of deformation. We describe the rotation of the web via a rigid rotation and a micro-rotation [1]. The micro-rotation encapsulates the information about the difference in rotation of the two ends of the web. We calculate the mechanical quantities associated with the extension and bending of the web in two separate sections below.

4.4.1 Bending and twisting of the web

Our objective in this section is to calculate the micro-rotation tensor $\mathbf{P}$. We attach a copy of $\mathbf{Q}_{0}$ say $\mathbf{Q_{0}^{\pm}}$ on the + and - end of every spoke in the reference configuration. $\mathbf{Q}_{0}^{\pm}$ change to $\mathbf{Q}^{\pm}$ in the current configuration. The ’difference’ between $\mathbf{Q_{0}^{+}}$ and $\mathbf{Q_{0}^{-}}$ gives the bending and torsion of the web while the ’average’ of $\mathbf{Q_{0}^{+}}$ and $\mathbf{Q_{0}^{-}}$ gives the rigid rotation of the web. We relate $\mathbf{Q^{\pm}}$ to the rotations of $\pm$ strands $\mathbf{R}(x)^{\pm}$. The angles between the columns of $\mathbf{Q}^{+}_{0}$ and $\mathbf{R}^{+}_{0}$ should remain same during the deformation which translates into the following condition.

[TABLE]

We are now in a position to calculate the micro-rotation $\mathbf{P}$ responsible for generating elastic moment in the web. Let the micro-rotation tensor in the reference configuration be $\mathbf{P}_{0}$ which changes to $\mathbf{P}$ during deformation. We use an expression for $\mathbf{P}/\mathbf{P}_{0}$ given in Moakher and Maddocks [1].

[TABLE]

This gives

[TABLE]

Note that $\mathbf{\Phi^{c}}$ is a skew symmetric tensor. The next step is to calculate the Gibbs rotation vector of $\mathbf{P}$ [1]. The Gibbs rotation vector $\bar{t}$ of a rotation matrix $\mathbf{T}$ is defined as $\bar{t}=\tan\frac{\theta}{2}\mathbf{k}$ such that $\text{tr}\mathbf{P}=1+2\cos\theta$ and $\mathbf{k}$ is a unit vector such that $\mathbf{T}\mathbf{k}=\mathbf{k}$. Consider $\bar{\mathbf{P}}=\mathbf{1}+\mathbf{\Phi^{c}}$ where $\mathbf{\Phi^{c}}\sim O(\varepsilon)$. The axis of the infinitesimal rotation $\bar{\mathbf{P}}$ is the axial vector of $\mathbf{\Phi^{c}}$. Hence,

[TABLE]

We can not calculate the magnitude of the rotation by taking $\text{tr}\bar{\mathbf{P}}$, since it gives $1+2\cos\theta=3$ which implies $\theta=0$. We consider the following limit.

[TABLE]

Now we take the trace of the RHS and get $\theta=|\boldsymbol{\phi^{c}}|$. Hence, the Gibbs rotation vector of $\bar{\mathbf{P}}$, $\boldsymbol{\bar{\eta}}$ is given as

[TABLE]

The Gibbs rotation vector of $\mathbf{P}$ is simply $\boldsymbol{\eta}=\mathbf{Z}\boldsymbol{\bar{\eta}}$. Note that in the undeformed state $\boldsymbol{\eta}_{0}=\boldsymbol{\bar{\eta}}_{0}=0$. We now proceed to calculate the rigid rotation of the spoke $\mathbf{Q}$.

[TABLE]

Here $\boldsymbol{\eta}\sim O(\varepsilon)$. Now, the micro-moment $\boldsymbol{c}$ is related linearly to the $\boldsymbol{\eta}$ via an elastic tensor $\mathbf{H}$.

[TABLE]

For further reference, let

[TABLE]

4.4.2 Extension of the web

The distance between the two strands is $\boldsymbol{w}=\frac{\mathbf{r}^{+}-\mathbf{r}^{-}}{2}$ and in the undeformed configuration $\mathbf{w}_{0}=\frac{\mathbf{r}^{+}_{0}-\mathbf{r}^{-}_{0}}{2}$. By direct calculation we observe

[TABLE]

where

[TABLE]

The force exerted by the $+$ strand on the $-$ strand $\boldsymbol{f}$ is given by,

[TABLE]

where $\mathbf{\bar{L}}$ is a tensor of mechanical properties of the web. This force $\boldsymbol{f}$ causes the web to extend and shear. For further reference let,

[TABLE]

4.4.3 Stacking energy

DNA consists of consecutive base-pairs stacked on top of each other in a regular fashion. The resistance to external forces and moments not only comes from the elastic deformation of the strands and the webbing but also from the change in alignment of the base-pairs. We call the energy associated with this change in bases’ position and spatial orientation ‘stacking energy’. Stacking energy plays a critical role in various phenomena such as melting of DNA [32, 33]. We prescribe a form of free energy which is quadratic in the twist $k_{3}$ and stretch $\xi$.

[TABLE]

There are other sophisticated expressions for the stacking energy [33], but we use the quadratic form for two reasons: one, the non-quadratic terms in the energy of [33] account for effects such as base-pair severing which are crucial to DNA melting which does not occur in our problem, two, a quadratic energy keeps our problem linear. This interaction energy results in a distributed body force $\boldsymbol{l}$ and distributed body moment $\boldsymbol{h}$ on the strands.

[TABLE]

4.5 Step 4: Governing equations

We are now in a position to solve the governing equations for the mechanics of our helical birod. These equations consist of balance of linear momentum and angular momentum for both the strands. In the balance equations eqn. (48b) and eqn. (49b):

• •

$\boldsymbol{m}^{\pm}=EI\kappa^{\pm}$ (eqn. 29) denotes the elastic moment in the $\pm$ strand. $\boldsymbol{n}^{\pm}$ are the contact forces for which there is no constitutive relation since the outer strands are assumed to be inextensible and unshearable.

• •

$\boldsymbol{f}$ and $\boldsymbol{c}$ are the distributed force and moment, respectively, exerted by the $+$ strand on the $-$ strand.

• •

$\boldsymbol{l}$ and $\boldsymbol{h}$ are the distributed force and moment exerted by base-pairs on the $+$ and $-$ strand.

The balance equations are:

[TABLE]

[TABLE]

Let $[\leavevmode\nobreak\ \mathbf{f}_{1}\quad\leavevmode\nobreak\ \mathbf{f}_{2}\quad\leavevmode\nobreak\ \mathbf{f}_{3}]=\mathbf{Z}\mathbf{Q}_{0}$. This gives

[TABLE]

We decompose the forces, $\boldsymbol{n}^{+}=(\boldsymbol{n}+\boldsymbol{n}^{c})\sim O(\varepsilon)$ and $\boldsymbol{n}^{-}=(\boldsymbol{n}-\boldsymbol{n}^{c})\sim O(\varepsilon)$. $\boldsymbol{n}=n_{1}\leavevmode\nobreak\ \mathbf{f}_{1}+n_{2}\leavevmode\nobreak\ \mathbf{f}_{2}+n_{3}\leavevmode\nobreak\ \mathbf{f}_{3}$ and $\boldsymbol{n}^{c}=n_{1}^{c}\leavevmode\nobreak\ \mathbf{f}_{1}+n_{2}^{c}\leavevmode\nobreak\ \mathbf{f}_{2}+n_{3}^{c}\leavevmode\nobreak\ \mathbf{f}_{3}$. Now, $\boldsymbol{n}_{x}=(n_{1x}-\omega n_{2})\leavevmode\nobreak\ \mathbf{f}_{1}+(n_{2x}+\omega n_{1})\leavevmode\nobreak\ \mathbf{f}_{2}+n_{3x}\leavevmode\nobreak\ \mathbf{f}_{3}+O(\varepsilon^{2})$. Similarly for $\boldsymbol{n}^{c}_{x}$. We use $\boldsymbol{c}=c_{1}\leavevmode\nobreak\ \mathbf{f}_{1}+c_{2}\leavevmode\nobreak\ \mathbf{f}_{2}+c_{3}\leavevmode\nobreak\ \mathbf{f}_{3}$ and $\boldsymbol{f}=f_{1}\leavevmode\nobreak\ \mathbf{f}_{1}+f_{2}\leavevmode\nobreak\ \mathbf{f}_{2}+f_{3}\leavevmode\nobreak\ \mathbf{f}_{3}$ from eqn. (44) and eqn. (41). Then, the balance equations become:

[TABLE]

We have 12 differential equations in the 12 unknowns $(r,f,\xi,k_{3},\beta^{+},\beta^{-},n^{c}_{1},n^{c}_{2},n^{c}_{3},n_{1},n_{2},n_{3})$. We substitute the following ansatz into the equations.

[TABLE]

This results in an eigenvalue problem. We find $23$ eigenvalues, but retain only $6$ for reasons explained in the appendix. Let those $6$ eigenvalues be $\pm\lambda,\pm\mu,\pm\delta$ and the corresponding eigenvectors $\mathbf{v}_{\pm\lambda}$ and $\mathbf{v}_{\pm\mu}$. Let

[TABLE]

Hence,

[TABLE]

Here, $p_{1}$, $p_{2}$, $p_{3}$, $p_{4}$, $p_{5}$ and $p_{6}$ are the constants which are determined using boundary conditions.

4.6 Step 5: Boundary conditions

We assume that the impact of a protein binding to DNA is two fold: a) the protein fixes the curvatures at the binding site as in [10, 18, 19], and b) the protein causes a change in the radius of the DNA helix [4] as shown in the inset of fig. 6 (b). Thus, we apply boundary conditions on the curvatures $k_{1}$, $k_{2}$ and the change in radius $r$ of the DNA helix. We discuss two cases, first, when one protein binds to the DNA, and second, when two proteins bind to it.

One protein: Let us assume that the protein binds at $x=0$. The boundary conditions for this case are:

[TABLE]

The second boundary condition says that the DNA is straight far away from the protein and that the perturbation in DNA radius occurs only in the vicinity of the bound protein. 2. 2.

Two proteins: Let us assume that the two proteins bind at $x=0$ and $x=a$, respectively. We divide our domain into three parts $-\infty<x<0$, $0<x<a$ and $a<x<\infty$ each of which has different boundary conditions attached to it.

Region 1: $x\in(-\infty,0)$

[TABLE]

Region 2: $x\in(0,a)$

[TABLE]

Region 3: $x\in(a,\infty)$

[TABLE]

4.7 Step 6: Energy of the birod

We assume small elastic deformations throughout, hence the resulting energy is quadratic in the strain variables. The elastic energy has contributions from the bending of the outer strands eqn. (29), the extension, bending and twisting of the web eqn. (42), (43) and the stacking energy eqn. (47).

[TABLE]

We are especially interested in the interaction energy $\Delta G$ which is the elastic energy of interactions between the two proteins.

[TABLE]

where $E^{2}_{a}$ is the energy of two proteins bound to DNA, one at $x=0$ and other at $x=a$, and $E^{1}_{a}$ and $E^{1}_{0}$ are the elastic energies corresponding to a single protein binding at $x=a$ and $x=0$, respectively.

5 Elastic constants

Our model has $9$ elastic constants $L_{1},L_{2},L_{3},H_{1},H_{2},H_{3},K_{c},K_{e},EI$. The experimental values for these constants are not known. In order to get some idea about the magnitude of the elastic constants we calculate the extensional modulus, torsional modulus and twist-stretch coupling modulus for a double-stranded DNA within our birod model. The explicit calculation is presented in the appendix. We choose

[TABLE]

This choice of elastic constants gives the extensional modulus $S\approx 1245$ pN, torsional modulus $C\approx 490$ pNnm2 and twist-stretch coupling modulus $g\approx-90$ pNnm which are close to actual values for ds-DNA [20] measured in experiments. We point out that this choice of elastic constants is not unique, nonetheless we use them to make further calculations.

When we substitute these constants into the governing equations (eqn. (51)) and solve the eigenvalue problem involving $\lambda$, we get the following eigenvalues

[TABLE]

Other eigenvalues are either very large ($\to\pm\infty$), very small ($\sim 0$) or purely imaginary. Purely imaginary and zero eigenvalues when substituted in $e^{\lambda x}$ give a sinusoidal and a constant function, respectively, which do not decay to zero as $x\to\pm\infty$. As mentioned in section 4.6, the curvatures $k_{1},k_{2}$ and change in radius $r$ must go to zero at $\pm\infty$. Thus, zero or purely imaginary eigenvalues cannot satisfy our boundary conditions, and are, therefore, not useful. We refer the reader to the appendix for further discussion on the choice of eigenvalues.

Consider a situation in which two proteins bind DNA, one at $x=0$ and other at $x=a$. In the region $a<x<\infty$ the solution eqn. (53) consists of only negative eigenvalues. There are three negative eigenvalues $\lambda_{2,3,4}$ and consequently three unknown constants. We have three boundary conditions on $k_{1}$, $k_{2}$ and $r$ at $x=a$ to determine those constants. Similarly in the region $-\infty<x<0$, the solution consists of only positive eigenvalues $\lambda_{7,8,9}$, so the constants can again be evaluated from three boundary conditions. We use this scheme to evaluate the strain parameters which we substitute into the expression for the elastic energy functional eqn. (55). Notice that the dominant eigenvalue $\pm 0.36$ nm*-1* corresponds to a decay length of $2.8$ nm ($\approx 10$ bp) which is what Kim et al.[4] report in their experiments.

6 Results

The experimental evidence for allosteric interactions when two proteins bind to DNA is documented in Kim et al. [4]. Many earlier papers have also described allostery in DNA, but Kim et al. present exquisite quantitative details which call for a quantitative explanation.

To unravel the physics behind these allosteric interactions, we begin by examining the case when one protein binds to DNA. As discussed in section 4.5, the strain variables ($r,\zeta,\beta^{\pm},k_{1,2,3}$) are linear combinations of decaying exponentials. For instance, consider $k_{3}(x)$ for a protein binding at $x=0$:

[TABLE]

where $\lambda=0.36\text{ nm}^{-1}$, $\mu=0.42\text{ nm}^{-1}$, and $\delta=0.68\text{ nm}^{-1}$. $\mathbf{v}_{\pm\lambda}$, $\mathbf{v}_{\pm\mu}$, and $\mathbf{v}_{\pm\delta}$ are the eigenvectors associated with eigenvalues $\pm\lambda$, $\pm\mu$, and $\pm\delta$, respectively. The constants $p_{i}$ and $q_{i}$ ($i=1,2,3$) are evaluated using the boundary conditions at $x=0$. It is not difficult to see that the strain variables decay to zero as $x\to\pm\infty$. We can replace $k_{3}$ in the above equation by other strain variables ($r,\xi,\beta^{\pm}$) and recover similar behavior. We discuss a few characteristics of the variation of the strain parameters as functions of position. The results are plotted in fig. 3 and fig. 4. The strain parameters ($r,k_{3},\beta^{\pm}$) decay exponentially with distance from the site of protein binding. The curvatures exhibit an exponentially decaying sinusoidal character with a period of $11$ bp. This periodic decay of the curvatures manifests itself as sinusoidal variations in the interaction energy. We find that these plots are slightly asymmetric about $x=0$. We attribute this to the structural asymmetry in the right-handed double-helix with phase angle $\alpha=2.1$ radians. If we choose the phase angle $\alpha=\pi$ radians instead, we find that the plots are exactly symmetric about the site of protein binding as shown in the appendix.

We now consider the case when two proteins bind to DNA, one at $x=0$ and the other at $x=a$. We proceed in a similar manner as above and express the strain profiles as linear combinations of exponentials:

[TABLE]

The constants $p_{i}$ and $q_{i}$ ($i=1,2,3$) are determined by three boundary conditions (on $k_{1},k_{2}$ and $r$) at $x=0$ and $x=a$, respectively. The constants $m_{j}$, ($j=1,2,3,4,5,6$) are determined by six boundary conditions at $x=0$ and $x=a$. The behavior of the strain variables for two proteins is similar to that for one protein as shown in fig. 5. When two proteins are separated by a large distance $a>10\times 3.4$ nm (i.e., more than 10 helical turns of DNA), the strain profile looks like a concatenation of the profiles of two proteins binding separately. Their strain fields do not interact at such distances, thus there is little interaction energy. When the distance decreases, the strain fields of the two proteins overlap, and this is responsible for the interaction energy.

As discussed in section 3, two defects on a straight ladder fig. 1 interact via an interaction energy that decays exponentially with the distance between them. Now, we focus on the double-helical birod and examine the behavior of different boundary conditions on the interaction energy $\Delta G$ in fig. 6. We assume for simplicity that both proteins apply the same boundary conditions on the DNA, the exact numerical values are given in the figure. If we choose the change in radius $r_{0}=0$ and apply the boundary conditions only on the two curvatures $k_{1},k_{2}$, the interaction energy decays exponentially while varying sinusoidally with a period of $5.5\approx 11/2$ bp. This case corresponds to proteins that bend DNA as shown in the inset of fig. LABEL:deltaG(b). On the other hand, if the curvatures $k_{1},k_{2}$ are zero while the change the change in radius $r_{0}$ is non-zero, we get an exponentially decaying profile devoid of any oscillatory character, which is similar to the results for the ladder in section 3. The exponentially decaying component originates from the elasticity of the web, and the sinusoidal behavior comes from the double-helical structure of DNA. From this exercise we conclude that in order to get a sinusoidally varying interaction energy a protein must change the local curvature in the DNA, a mere change in radius of the DNA is not sufficient to give rise to the interaction energy profiles observed in experiments.

In our model the magnitude of the interaction energy increases monotonically with increase in the magnitude of the changes in curvatures or radius caused by the two proteins. Thus, by systematically varying the boundary conditions imposed by the proteins we can establish agreement of our theoretical results for $\Delta G$ with the experimental values documented by Kim et al.[4]. This is done in fig. 6(b). The values of the curvatures that give the best fit to the experimental data are $k_{11}=k_{21}=0.02$ $nm^{-1}$, $k_{12}=k_{22}=0.05$ $nm^{-1}$ and $r_{1}=-r_{2}=0.02\text{ nm}$. This choice is, however, not unique and it is coupled with the choice of stiffnesses of the webbing in our birod model. Be that as it may, our exercise above demonstrates that a birod model can capture the dependence of interaction enery on the distance between proteins bound to DNA. Calibration of the model and faithfully connecting it to experiment will require deeper analysis, and perhaps also, computation.

The period the interaction energy in fig. 6(a) is approximately 5.5 bp while that in fig. 6(b) is 11 bp as in the experiment. Why? Note that the strain variables in a two protein complex shown in fig. 7 (b) are a function of both the parameter $x$ and the distance between the two proteins $a$. We fix $x$ ($=2$ nm from protein P1) and focus on the dependence on $a$. We assume that both the proteins apply identical boundary conditions. If the proteins do not cause any change in the radius such that $r_{0}=0$, then the strain parameters involved in the elastic energy (eqn. (55)) $\propto e^{-\Gamma x}\psi(\omega a)$, where $\psi(\omega a)$ is a sinusoidal function oscillating with a period 11 bp, and the elastic energy of the two protein complex $\propto(e^{-\Gamma x}\psi(\omega a))^{2}$ oscillates with a period 5.5 bp. On the other hand, when the protein causes both a change in radius $r_{0}$ and a change in curvature $k_{20}$, the strain variables are $\propto(e^{-\Gamma_{1}a}\psi(\omega a)+e^{-\Gamma_{2}x})$ and the elastic energy of the two protein complex $\propto(e^{-\Gamma_{1}a}\psi(\omega a)+e^{-\Gamma_{2}a})^{2}$ oscillates with a period of 11 bp due to the cross term $e^{-(\Gamma_{1}+\Gamma_{2})a}\psi(\omega a)$. We plot the interaction energy $\Delta G(a)$ between the two proteins constituting the protein complex in fig. 7(b) and verify the periods for respective boundary conditions which resolves the apparent discrepancy in the periods in fig.6 (a) and (b). As a final application of our birod model we examine the sequence dependence of the allosteric interaction energy $\Delta G$. While there is overwhelming qualitative evidence, both experimental [4] and numerical [21], showing that AT-rich sequences exhibit stronger allosteric interactions compared to GC-rich ones, a theoretical explanation is still lacking. Stronger interactions are associated with longer decay lengths. Using our theory we can find the dependence of the decay length on the elastic constants of the web. Since, AT base-pairs consist of two hydrogen bonds, the corresponding elastic constants for the web are expected to be lower than GC base-pairs which comprise of three hydrogen bonds. In an attempt to simulate such a scenario we replace the elastic constants for the web ($K_{c},K_{e},L_{i},H_{i}$ $i=1,2,3$) in eqn. (57) with ($\chi K_{c},\chi K_{e},\chi L_{i},\chi H_{i}$ $i=1,2,3$) while keeping $EI$ fixed, and vary the parameter $\chi$ in the range $0.5\leq\chi\leq 1$. We define a measure of the decay length $l_{d}$ to be the inverse of the eigenvalue having the least non-zero magnitude, obtained in eqn. (58). For instance, if $\chi=1$, decay length $l_{d}=1/0.34$ nm$\approx 10$ bp. We plot the variation of $l_{d}$ with $\chi$ in fig. 8. We find that the decay length increases with the decrease in elastic constants of the web. We plot $\log l_{d}$ versus $\log\chi$ and deduce that $l_{d}\sim\frac{1}{\chi^{2/3}}$.

7 Conclusion

Kim et al. [4] have presented compelling quantitative evidence for allosteric interactions between two proteins bound to DNA at distant locations. They showed that the interaction energy for two proteins separated by distance $a$ on DNA is a decaying exponential oscillating with period of 11 bp. Various attempts to numerically simulate the allosteric interactions have been made [21, 22] and have associated the oscillating interaction energy to the major groove width in the double-helical structure of DNA. We approach the problem from a purely mechanical standpoint. We conjecture that the local deformation field in DNA caused by a bound protein is similar to that produced by a defect in an elastic solid. We begin by computing the interaction energy for two defects on a ladder and find that it decays exponentially with the distance between them. We, then, proceed to replicate the same calculation for DNA by modelling it as a double-helical birod [1]. We assume that the outer phosphate backbones represented by $\pm$ strands to be inextensible and unshearable while the base-pairs are capable of elastic extension, shear, bending and, twisting. We assume a general form of displacement for these strands (eqn. 19) which we use to calculate the micro-displacement and micro-rotation for the base-pairs. We, then, use these expressions to solve the governing equations for our birod. A crucial factor in our treatment is the boundary conditions. We follow Kwiecinski et al. [19], Kim et al. [4] and Liang and Purohit [11] and impose boundary conditions on the curvatures and the radius of the DNA double-helix. The question, “what kind of boundary conditions a protein could possibly apply”, is not yet comprehensively addressed in the literature and is not the central issue of this study either. Rather our message is that after solving the governing equations and plugging in boundary conditions, we recover the exponentially decaying profile that oscillates with a period of 11 bp. We end by examining the sequence dependence of allosteric interactions and show that AT-rich sequences exhibit stronger interactions than GC-rich sequences.

Even though our birod model does surprisingly well by capturing the dependence of interaction energy on distance there are many important caveats that we must point out. First, we do not expect our birod model to be accurate near the site of protein binding. The deformations near the binding site could be large enough that a linear elastic theory may not be applicable. Our assumptions that the outer strands are inextensible and the web is elastic could also break down in the vicinity of the binding site. Second, we have little knowledge of the elastic constants of the web. We have assumed some stiffness parameters for the web that gave the right experimentally verified moduli for the DNA, but there could have been another set of parameters that would have given similar results. One may have to appeal to molecular simulations [2, 34, 35, 36, 37] to get these parameters. Third, the boundary conditions applied by the proteins on the DNA are not clear. One may have to look for guidance from molecular simulations or protein-DNA co-crystal structures to get a clearer picture. Finally, we have not accounted for fluctuations or entropic interactions in our model. This is partly justifiable because the length of DNA between two protein binding sites for which significant allosteric interactions are observed is often much smaller than the persistence length of the DNA. However, a rigorous calculation should be done to verify this assumption. In spite of these shortcomings, our model could provide a starting point for analyzing allosteric interactions in DNA within the broad framework of configurational forces in elastic solids.

*We acknowledge insightful discussion with Yujie Sun who is one of the authors in Kim et al. [4].

Appendix A Appendix

A.1 Kinematics of the $-$ strand

In the main text we gave detailed derivations for the strains, curvatures, etc., for the $+$ strand in our birod. We now shift our attention to the complimentary $-$ strand. The reference configuration of this strand is denoted by position vector $\mathbf{r}^{-}_{0}$.

[TABLE]

Along the same lines as the $+$ strand, we conceive the deformed configuration to be a helix wrapped around a curved axis defined by curvatures $k_{1},k_{2}$ and $k_{3}$ along the directors $\leavevmode\nobreak\ \mathbf{d}_{1},\leavevmode\nobreak\ \mathbf{d}_{2}$ and $\leavevmode\nobreak\ \mathbf{d}_{3}$, respectively.

[TABLE]

We use the same apparatus mutatis mutandis described for the $+$ strand to calculate various quantities of interest. The results are:

[TABLE]

where $\mathbf{\Theta^{-}}$ is a skew symmetric tensor.

[TABLE]

We compute curvature $\kappa^{-}$ as follows,

[TABLE]

We obtain the moment $\boldsymbol{m}^{-}$ as follows,

[TABLE]

where $\leavevmode\nobreak\ \mathbf{f}_{1},\leavevmode\nobreak\ \mathbf{f}_{2},\leavevmode\nobreak\ \mathbf{f}_{3}$ are given by eqn. 40.

A.2 Evaluation of material properties of the web

In this section, we consider a deformation of the double-helical structure induced by a stretching force $F$ and torque $T$ on one end. We assume that the helix retains its helical configuration, but with changed geometrical parameters. Thus, $r$, $\beta$ and $e$ are independent of $x$. Our goal is to compute the strains and curvatures, then evaluate the energy, and then identify the stretch modulus, twist modulus and twist-stretch coupling modulus of the double-helical structure from this energy expression. The computation of strains, curvatures, etc., of the helix proceeds as in the main text.

[TABLE]

We assume $r,\beta,e\sim O(\varepsilon)$, hence

[TABLE]

The inextensibility condition gives,

[TABLE]

$\boldsymbol{t}^{+}_{0}$, $\boldsymbol{n}^{+}_{0}$ and $\boldsymbol{b}^{+}_{0}$ are the tangent, normal and binormal to the $+$ strand in the reference configuration. We calculate tangent $\boldsymbol{t}^{+}$ to the deformed configuration.

[TABLE]

Next, we calculate the curvature $\kappa^{+}$.

[TABLE]

We go on to calculate the normal in the deformed configuration $\boldsymbol{n}^{+}$.

[TABLE]

We are now in a position to calculate the deformed Frenet-Serret frame $\mathbf{R}^{+}$.

[TABLE]

where $\mathbf{\Theta^{+}}$ is a skew symmetric tensor.

[TABLE]

For the negative strand we follow the same procedure.

[TABLE]

After performing all the calculations

[TABLE]

We substitute $r=-\frac{e}{\omega^{2}a}-a\beta$ from eqn. (69) and compute the elastic constants as follows.

[TABLE]

Then, by trial and error we pick values of $L_{1},L_{2},L_{3},H_{1},H_{2},H_{3},K_{c},K_{e},EI$ to match the $S,g,C$ known from experiments. Our choice of the material parameters $L_{1},H_{2},K_{c}$, etc., is not unique.

A.3 Choice of eigenvalues obtained in section 5

In section 5, we solve the governing differential equation eqn. 51 by substituting $y(x)=y_{0}e^{-\lambda x}$ where $y=(r,f,\xi,k_{3},\beta^{\pm},n^{c}_{i},n_{i})\quad$ $i=1,2,3$. We look for the values of $\lambda$ corresponding to a non-trivial solution of the governing equations. For this we need to solve the eigenvalue problem $\mathcal{A}(\lambda)\mathbf{v}_{0}=0$, where $\mathcal{A}$ is a function of $\lambda$ and elastic constants (eqn. 57) and $\mathbf{v}_{0}=[r_{0},f_{0},\xi_{0},k_{30},\beta^{+}_{0},\beta^{-}_{0},n^{c}_{i0},n_{i0}]^{T}\quad$ $i=1,2,3$. We set $\det\mathcal{A}(\lambda)=0$ and get following solutions for $\lambda$.

[TABLE]

Among these 23 eigenvalues we neglect the eigenvalues $x_{1,2,3,4,5,20,21,22,23}$ whose magnitude is $>10^{3}$ because the corresponding decay length is tiny which leads to large numerical errors given that we need to compute third derivatives. Then, there are small eigenvalues $x_{9,16}$ whose magnitude is close to zero $(<10^{-3})$ and purely imaginary eigenvalues $x_{10,11,12,13,14,15}$ which when substituted in $e^{-\lambda x}$ result in a constant or a sinusoidal function, respectively, that do not decay to [math] as $x\to\pm\infty$. Hence, we must neglect these too. This leaves us with $x_{6,7,8,17,18,19}$, which are used in section 5.

A.4 Results for $\alpha=\pi$ radians