DNA fluctuations reveal the size and dynamics of topological domains

Abstract DNA supercoiling is a key regulatory mechanism that orchestrates DNA readout, recombination, and genome maintenance. DNA-binding proteins often mediate these processes by bringing two distant DNA sites together, thereby inducing (transient) topological domains. In order to understand the dynamics and molecular architecture of protein-induced topological domains in DNA, quantitative and time-resolved approaches are required. Here, we present a methodology to determine the size and dynamics of topological domains in supercoiled DNA in real time and at the single-molecule level. Our approach is based on quantifying the extension fluctuations—in addition to the mean extension—of supercoiled DNA in magnetic tweezers (MT). Using a combination of high-speed MT experiments, Monte Carlo simulations, and analytical theory, we map out the dependence of DNA extension fluctuations as a function of supercoiling density and external force. We find that in the plectonemic regime, the extension variance increases linearly with increasing supercoiling density and show how this enables us to determine the formation and size of topological domains. In addition, we demonstrate how the transient (partial) dissociation of DNA-bridging proteins results in the dynamic sampling of different topological states, which allows us to deduce the torsional stiffness of the plectonemic state and the kinetics of protein-plectoneme interactions. We expect our results to further the understanding and optimization of magnetic tweezer measurements and to enable quantification of the dynamics and reaction pathways of DNA processing enzymes in the context of physiologically relevant forces and supercoiling densities.


Supplementary Methods
A. DNA and restriction enzymes. All measurements used a 7.9-kbp DNA construct described previously (1,2). For specific attachment of the DNA to the magnetic bead and the flow cell surface, 600-bp PCR-generated DNA fragments labeled with multiple biotin and digoxigenin moieties, respectively, were ligated to the DNA. The DNA construct was attached to 1.0 µm-diameter MyOne, streptavidin coated beads (Life Technologies, USA) by incubating 0.5 µL of picomolar DNA stock solution and 2 µL MyOne beads in 150 µL 1× PBS (Sigma-Aldrich, USA; ≈ 150 mM monovalent salt, pH 7.4) for 5 min. Experiments with restriction enzymes were carried out in a Ca 2+ -containing buffer (50 mM potassium acetate, 20 mM Tris acetate, 10 mM calcium acetate, 100 µg/ml BSA, pH 7.0) based on the commercial Cutsmart buffer (NEB) and using commercial preparations of the enzymes BsaXI, NaeI, and SacII (NEB). The enzyme stock solutions were diluted 20, 000 − 50, 000x to ensure that only approximately 30% of all tethers show interaction, thereby minimizing the possibility of tether interactions with more than two enzyme complexes simultaneously.
B. Magnetic tweezers set up. The custom-built MT setup uses a pair of 5 × 5 × 5 mm 3 permanent magnets (W-05-N50-G, Supermagnete, Switzerland), oriented in a vertical configuration and separated by a 1 mm gap (3). The distance between beads and magnets is controlled by a DC-Motor (M-126.PD2, PI, Germany), and magnet rotation is performed by another DC-Motor (C-150.PD, PI, Germany). Beads are observed with a 40× oil immersion objective (UPLFLN 40×, Olympus, Japan) and imaged with a CMOS sensor camera (12M Falcon2, Teledyne Dalsa, Canada). By reducing the field of view to 5% of the original area (to 1792 × 282 pixels, with pixelsize ≈ 110 nm) a frame rate of 1 kHz is achieved. Images are transferred to a frame grabber (PCIe 1433, NI, USA) and analyzed in real-time with a custom-written tracking software (4). A LED (69647, Lumitronix LED Technik GmbH, Germany) is used for illumination and a piezo stage (Pifoc P-726.1CD, PI, Germany) moves the objective to produce the look-up table (LUT) to enable bead tracking during experiment. Forces were calibrated from the transverse bead fluctuations as described (5).
Flow cells were assembled from two microscope coverslips (24 × 60 mm, Carl Roth, Germany). Prior to assembly, the bottom coverslip was coated with (3-glycidoxypropyl)trimethoxysilane (abcr GmbH, Germany), and 50 µL of a 5000× diluted stock solution of polystyrene reference beads (Polysciences, USA) in ethanol (Carl Roth, Germany) was deposited on the silanized slides, slowly dried, and baked at 80°C for 1 min. The top coverslip was processed using a laser cutter, producing openings with 1 mm radius, for liquid exchange. The two coverslips were glued together by a single layer of melted Parafilm (Carl Roth, Germany), precut to form a ∼ 50 µL channel connecting the inlet and outlet opening of the flow cell. Following flow cell assembly, 100 µg/ml anti-digoxigenin (Roche, Switzerland) in 1× PBS was introduced and incubated for 2 h. To minimize nonspecific interactions, the flow cell was flushed with 800 µL of 25 mg/mL bovine serum albumin (Carl Roth, Germany), incubated for 1 h and rinsed with 1 mL of 1× PBS. Subsequently, 50 µL of the bead-coupled DNA constructs were introduced into a flow cell (see above), and allowed to bind for 5 min. Finally, the flowcell was rinsed with 2 mL of 1× PBS to flush out unbound beads, and the magnet was mounted to constrain the supercoiling density of the tethers and to apply an upward force on the beads.

C. Magnetic tweezers measurements.
Prior to each measurement, selected beads were evaluated for the presence of multiple tethers and for torsional constraint. The presence of multiple tethers was assessed by introducing negative supercoiling under high tension (f = 5 pN.) For a single DNA tether, the formation of plectonemes at negative linking differences and f = 5 pN is suppressed due to DNA melting and no height change is observed. In contrast, for the case of multiple tethers, introduction of negative turns braids the DNA tethers, which decreases the z-extension of the bead. Beads bound by multiple tethers are discarded from further analysis. The evaluation of extension fluctuations in the absence of protein is performed in PBS buffer, by recording z-extension at different numbers of applied turns (i.e. changing ∆Lk) for 180 s each. Subsequently, measurements are repeated at different forces.
Prior to flusing enzymes into the flow cell, we introduce 1 mL 10 mM Tris-HCl (pH = 8.0) buffer to replace the phosphate buffer that would otherwise result in the formation of precipitates due to complexation with Ca 2+ in the assay buffer. Next, the assay buffer is introduced in the flow cell, followed by the application of a linking difference ∆Lk +25 in the tethers at 0.5 pN tension. The z-extension of supercoiled DNA tethers in assay buffer is recorded for > 10 min, to verify the absence of anomalous behavior. Next, the z-translation motor that controls the magnet position is moved closer towards the flow cell, to apply a force of 5 pN. Enzyme solution (100 µL) is then flushed in the flow cell at a flow rate ∼ 150 µL min −1 , after which the force is reduced to its original value (0.5 pN).

D. Analysis of magnetic tweezers data.
Real-time tracking was performed using the open source software framework developed previously (4). This framework employs the Quadrant Interpolation algorithm to enable accurate and simultaneous tracking of many beads in parallel (6). Further processing of the MT data was carried out using custom-written Matlab and Python routines.
E. Monte Carlo simulations: Model set up. We performed Monte Carlo simulations of a self-avoiding twistable wormlike chain (TWLC) model, similar to the models employed previously (7)(8)(9)(10). The DNA is represented by a chain of coarse-grained beads, each corresponding to a segment of ten base pairs of the DNA molecule (see Fig. S4 (a)). Bend-and twist-deformations are traced by assigning local orthogonal frames of reference (triads),T = { e1, e2, e3}, to each bead, such that the unit vector e3 is oriented along the direction connecting a given bead with the next bead and the remaining vectors, e1 and e2, point in directions perpendicular to the chain contour. The relative rotation that maps consecutive triads into each other is captured by a rotation vector Ω, pointing along the rotation axis and having magnitude equal to the rotation angle. Expressing the three components of this vector in the coordinate system (specified by the respective triad) of the first of these two triads allows for the definition of two bending components, Ω1 and Ω2, (which are the components along the vectors e1 and e2 respectively) and a twist mode Ω3. Further details on the definition of these rotation fields can be found in Refs. (8), and (11).
The TWLC can be expressed in terms of these rotation fields. In its discretized version, the elastic energy of the TWLC is given by where a is the discretization length, and A and C are the bending and torsional stiffnesses as used in Eq. (7) and β = 1/kBT . Stretching forces such as those applied by magnetic tweezers may be included by adding a contribution −f · R to the energy, with R the end-to-end vector and f the force. Furthermore, to appropriately model DNA supercoiling the repulsive electrostatic interactions of the negatively charged DNA backbones are crucial, as they control the coiling density of plectonemic supercoils. Following Rybenkov et al. (12), we model ion screened electrostatic interactions by effective excluded volumes (hard-sphere potentials) attached to the chain-monomers. Saline conditions of the solvent may be taken into account by an appropriate scaling of the bead radii. Combining these contributions, the total energy of a given chain configuration is given by where describes the excluded-volume interaction and where d EV is the effective diameter of the excluded volume beads. We design our simulations to closely resemble the setup of molecules within magnetic tweezers in the fixed linking number ensemble. The molecules are tethered between a surface and a magnetic bead, and a force is applied to the bead in the direction perpendicular to the surface. Fixing the bead orientation prevents twist from diffusing into or out of the system. Furthermore, the magnetic bead is sufficiently large that under the range of applied forces, the DNA chain does not pass around the bead. In the simulations, the bead rotation-constraint is imposed by fixing the orientations of the first and last triads. The bead and surface, respectively, are represented by confining the chain between two parallel impenetrable surfaces attached to the chain termini. Furthermore, since the Monte Carlo cluster moves are discrete and may lead to significant positional changes of individual bead segments, a given move may lead to an internal chain crossing, which would change the linking number by a multiple of 2 linking units. Such linking-violations are avoided by linearly interpolating the trajectory of all moved beads (specified by the excluded volume attached to every segment) between their initial and final position and ascertaining whether any of the beads violate the excluded volume constraint with any other bead along their respective trajectories. Moves leading to a violation of that kind are always rejected. A snapshot of a typical supercoiled configuration containing 2 plectonemic supercoils is shown in Fig. S4 (d).

G. Monte Carlo simulations: Parameterization.
As input, the model requires three parameters: the bending-and torsional stiffnesses A and C as well as the effective bead diameter d EV . We choose these parameters by optimizing the agreement between simulated and experimental rotation curves. After exploring a grid of parameters, we find A = 40 nm, C = 100 nm and a bead diameter of d EV = 4.0 nm, to yield the best agreement. A direct comparison between simulated and experimental rotation curves is shown in Fig. 1(c) of the main text and a comparison highlighting the force-dependent curvatures in the pre-buckling regime as well as the slope in the post-buckling regime of extension and variance is given in Fig. 1(d). We note that the values of A = 40 nm and C = 100 nm fall well into the range of experimentally determined values for the bending and torsional stiffness.

Supplementary Text
A. Phenomenological analysis of extension fluctuations. Here we present a phenomenological discussion of the extension variance, ∆z 2 , of supercoiled DNA under tension. The energy of a given configuration of DNA subject to force f and with end-to-end distance, z, in the direction of the applied force is given by E = E0 − f z. Here, E0 accounts for the elastic contributions of bending and twisting, while the term f z is the contribution of the work due to the applied force. The above form of E implies that averages of z, z 2 . . . can be expressed as suitable derivatives of the system free energy with respect to f , since differentiating the Boltzmann weight d n exp(−E/kBT )/df n introduces multiplicative factors of the form (z/kBT ) n . In particular, for the average extension and variance, these relations hold: where F is the free energy per unit length and L the total length of the molecule. Note that the variance is obtained by differentiating z with respect to f . Let us consider the following expression for the average extension with Q, W , Γ, σs and σp parameters. 0 ≤ σ < σs and σs < σ ≤ σp define the pre-and post-buckling regimes, respectively. Although Eq. (5) can be obtained from specific models (see (14,15) and Sec. B) we use it as a simple phenomenological expression to fit experimental data at low forces, f < 1 pN, where z is well-approximated by a symmetric function of σ.
Eq. (4) implies that z and ∆z 2 (= kBT d z /df ) will have the same functional dependence on σ. The parameter Q(f ) is negative, as experiments show that z in the pre-buckling regime is described by a concave parabola. Moreover, Q(f ) is a monotonically increasing function of f , as the curvature in z decreases as the force is increased, hence ∂ f Q > 0. Experiments in the post-buckling regime show that the DNA extension decreases in the post-buckling regime with increasing supercoiling density σ. Moreover, the absolute value of the associated slopes decreases with increasing force. Therefore, Γ is positive and a decreasing function of f , i.e. ∂ f Γ < 0.
In conclusion, this analysis shows that ∆z 2 is an increasing function of the supercoiling density both in the pre-and post-buckling phases and exhibits quadratic behavior pre-buckling and a linear increase with σ post-buckling.
B. DNA extension fluctuations from the two phase model. In the two phase model of linear DNA supercoiling (15) a DNA molecule of length L stretched by a force f and with a supercoil density σ is assumed to be composed of two phases: The stretched phase with length νL and the plectonemic phase with length (1 − ν)L. Stretched and plectonemic phases can have different supercoiling densities (φ and ψ, respectively) such that the total σ = νφ + (1 − ν)ψ is fixed. Introducing S(φ) and P(ψ), the free energies per unit length of the stretched and plectonemic phases (see a concrete example below), one minimizes the total free energy of the system, which fixes the parameters ν, φ and ψ. One finds that for 0 ≤ σ < σs the DNA molecule is in the pure stretched phase with ν = 1, φ = σ (pre-buckling), while for σs ≤ σ ≤ σp the molecule phase-separates into a stretched phase and a plectonemic phase with 0 < ν < 1 (post-buckling). In the latter case, the average supercoiling densities of the two phases are φ = σs and ψ = σp. For the calculations, we used the following free energies (15, 16) with ω0 = 1.75 nm −1 the intrinsic helical twist. Here S is the free energy of the twistable wormlike chain under stretching force f and fixed supercoiling density φ. Eq. (7) is a large force expansion, valid for kBT /f A 1. This is a good approximation for f > 0.5 pN. The parameters A and C are the bending and torsional stiffnesses of DNA. The free energy of the plectonemic phase Eq. (8) is purely phenomenological. It is characterized by a single parameter P known as the effective torsional stiffness of the plectoneme. The double-tangent construction (minimization) leads to the following free energy at post-buckling (15) The mean extension in either phase is simply the force derivative of the respective free energy Eq. (7) and Eq. (9) while differentiating once more yields the variance Eq. (4) [11] Eq. (10) and Eq. (11) have been used to produce the solid lines (theory) in Fig. 1(c) of the main text.
C. Linking number fluctuations in plectonemic loop. In this section we show how the effective torsional stiffness of the plectonemic state P can be deduced from the linking number distribution within protein constrained loops, assuming that the relative occupancy of the various states solely depends on the free energy of the supercoiling state of the DNA and not on the binding free energy of the protein (i.e. if the protein binding affinity is independent of the torque within the DNA template). Fig. S5 illustrates the looped and unlooped domains formed upon protein binding of lengths ∆L and L * = L − ∆L. The total excess linking number ∆Lk is partitioned between the two domains: ∆Lk loop and ∆Lk − ∆Lk loop are the linking number of the looped and unlooped parts. The total free energy then takes the form: where F l are Fu are the free energy of the looped and unlooped domains, respectively, while F b is the protein-binding free energy. The latter is assumed to be independent of linking number and need not be considered in what follows. The looped domain is in the pure plectonemic state and its free energy is given by Eq. (8) where we used ω0σ = 2π∆Lk L to convert supercoiling density to linking number. The unlooped domain consists of coexisting plectonemic and stretched phases, and its free energy is given by Eq. (9). This free energy is linear in the supercoiling density, thus is a linear function of ∆Lk loop and it will therefore not affect fluctuations of ∆Lk loop which are determined by the quadratic term Eq. (13). Equipartition thus yields ∆∆Lk loop 2 = ∆L 4π 2 P , [14] where we have defined ∆∆Lk loop ≡ ∆Lk loop − ∆Lk loop . Accordingly, fluctuations of the linking number in the loop are fully controlled by P the effective torsional stiffness of the plectoneme, which can be readily deduced from the analysis of the linking number fluctuations.

D. Master equation description of protein dissociation and reassociation.
We interpret the steps observed in the extension time traces for NaeI as linking number exchanges between looped and unlooped domains. These exchanges should happen through (partial) dissociation of the protein, followed by the rebinding in a different conformation. We write the generic Master equation for the process described schematically in Fig. S9, assuming that there is a single dissociated state and several protein-bound states. In this model, there is no direct transition between two protein-bound states. Indicating with w * (t) and wi(t) the probability of being in the dissociated and ith protein bound state, respectively, the Master equation then takes the form [15] dwi(t) dt = P ( * → i)w * (t) − P (i → * )wi(t) (i labels protein bound states). [16] The P (i → * ) and P ( * → i) are the unbinding and binding transition rates, respectively. We use P (i → * ) = κ independent on i, i.e. we assume that dissociation is not influenced by the supercoiled state of DNA. The reverse rate is then fixed by detailed balance P (i → * ) = κe −β∆E i , with β = 1/kBT the inverse temperature and ∆Ei the energy difference between the dissociated state and the ith bound state. This energy accounts for the protein binding energy and for the supercoiling contribution, the latter being different for different linking number partitionings between looped and unlooped domains. Using this definition of the rates and the normalization condition, i wi = 1 − w * , Eq. (15) takes the form with Z = i exp(−βEi) the equilibrium partition function for the bound states. The solution of Eq. (17) for the initial condition w * (0) = 0 is w * (t) = 1 − e −κ(1+Z)t 1 + Z . [18] At long times, κ(1 + Z)t 1, this evolves to the stationary value w * = (1 + Z) −1 . Inserting this limit value in Eq. (16) we obtain . [19] As MT traces cannot detect transition to the dissociated state, we rewrite the process with an effective Master equation, where this state is integrated out. For this purpose, we use the stationary solution for w * and rewrite the normalization condition as follows: The renormalized probabilities are larger than the original probabilities (pi > wi), because they absorb the "unobservable" dissociated state. Multiplying both sides of Eq. (19) by (1 + Z)/Z and using Eq. (20) we get . [21] This effective Master equation is characterized by direct transitions between bound states with associated rates, This effective rate is indicated by a dashed arrow in Fig. S9. The escape rate from a given ith protein-bound state is given by This is smaller than the intrinsic protein dissociation rate κ because, when compared to the full model, the effective transition i → j involves at least two steps i → * → j, but also multiple re-associations to the same state, such as i → * → i → * → j, i → * → i → * → i → * → j and so on . . . As a consequence, the characteristic dwell time, τi ≡ κ −1 i , for the ith protein-bound state measured in a MT experiment is longer than the protein dissociation rate τp ≡ κ −1 . Eq. (23) gives, [24] which is the relation reported in the main paper, where we have used pi = e −β∆E i /Z for the equilibrium probability distribution.

E. Effect of linking number exchange on z and ∆z 2 .
In the main text, we have derived the effect of protein binding on the extension average, z , and variance, ∆z 2 . We have shown that z does not change after protein binding and that ∆z 2 drops by an amount proportional to the length of the looped domain ∆L. This result is based on the assumption that the looped domain has a supercoiling density equal to σp. However, the amount of linking number in a given loop prior to protein binding is subject to fluctuations, such that the supercoiling density in the looped domain may assume the slightly different value σp + ∆σp. Here, we reexamine the calculation for the change in mean extension and variance for this more general assumption. We use the same notation as in the main text: L * = L − ∆L is the length of the unlooped domain and σ * its supercoiling density. One has The excess linking number of a domain of length ∆L and supercoiling density σp is given by which implies that a change of ∆∆Lk loop units of linking number in the looped domain amounts to ∆L∆σp = 2π ω0 ∆∆Lk loop . [27] Combining Eqs. (25) and (27) with Eq. (4) of the main text we find for the average extension after protein binding The mean extension does not remain constant but exhibits a small change proportional to ∆∆Lk loop . Using the estimate, Γ = 15, from experimental data for f = 0.5 pN (from Fig. 1 of the main text), we find 2πΓ/ω0 = 53 nm. This value is consistent with the jump size observed for NaeI (Fig. 5 of the main text). A similar calculation can be performed for the variance. For a change in excess linking number of ∆∆Lk loop in the looped domain, we find a variance drop upon protein binding of [29] The first term on the right-hand side of the previous equation is the term discussed in the main paper, showing a drop in the variance proportional to ∆L. For the rotation curve subject to a stretching force of 0.5 pN, we can estimate −kBT ∂Γ/∂f = 32 nm from the slope of the variance data at post-buckling (see Eq. (5) of the main paper). This yields a contribution to the variance of −2πkBT ω −1 0 ∂Γ/∂f = 114 nm 2 per unit of ∆∆Lk loop . While an integer change of linking number amounts to a change of about 5% in the extension, the variance exhibits a significantly lower relative change of only about 1%. Therefore, it is a better strategy to rely on extension changes to record topology-changing events. Table   Table S1. Curvatures and slopes in the pre-and postbuckling regimes, respectively, from both experimental ("MT") and simulated ("MC") rotation curves. These are the same data as in Figure 1d.   (18)(19)(20) and polymerases (21)(22)(23). In such experiments that use extension changes in the plectonemic regime to deduce protein-induced linking number changes ∆Lk, the signal in the experiment is proportional to the slope of z versus ∆Lk, i.e. the slope in the linear regime of extension versus linking number, which decreases with force, approximately as S ∼ f −1/2 (Figure 1d). Conversely, the noise N is given by the fluctuations in z (N = ∆z 2 ), which also decrease with force, with a scaling similar to N ∼ f −3/4 . In a first approximation, we estimate the noise in these experiments as the mean standard deviation of the extension fluctuations over the entire plectonemic regime (Figure 1c, bottom panel, the region indicated by the turquoise lines). For each force, we quantified the signal-to-noise ratio for a minimum of 5 individual DNA molecules, the symbols and error bars are the mean ± standard deviation. The force-dependence of the signal-to-noise ratio was fitted using a power law (red dashed line), which yields a pre-exponential factor 0.76 ± 0.05 and a scaling exponent 0.21 ± 0.06 (errors are 95% confidence intervals). From the scaling considerations, the signal-to-noise (black symbols) is expected to increases with force, albeit slowly, as S/N ∼ f 1/4 , which is in quantitative agreement with the scaling exponent found from our experimental data. A similar reasoning predicts a scaling of the signal-to-noise ratio with DNA contour length L as S/N ∼ L −1/2 . Taken together, the detection of a given change in ∆Lk is facilitated by measuring with short DNA constructs at high forces.    S6. Experimental determination of the supercoiling in the plectonemic regime ∆Lk P and the proportionality factor linking ∆L and ∆ ∆z 2 from rotation curve data. (a) Mean extension z as a function of linking difference ∆Lk. Extrapolation of the data in the plectonemic regime to zero extension gives the supercoiling in the plectonemic regime ∆Lk P . (b) Extension variance ∆z 2 as a function of linking difference ∆Lk. Extrapolation of the data in the plectonemic regime to ∆Lk = ∆Lk P gives the proportionality factor k B T · Γ · (∂σ P /∂f ). (c) Variance of transverse bead fluctuations along the external magnetic field lines ∆x 2 as a function of linking difference ∆Lk. Extrapolation of the data in the plectonemic regime to ∆x 2 = 0 gives the supercoiling in the plectonemic regime ∆Lk P . (d) Extension variance ∆z 2 as a function of linking difference ∆Lk . Extrapolation of the data in the plectonemic regime to ∆Lk = ∆Lk P gives again the proportionality factor k B T · Γ · (∂σ P /∂f ). (e) Comparison of values of ∆Lk P derived for different beads using either the experimental rotation curve data of extension z versus ∆Lk or variance of the transverse fluctuations ∆x 2 versus ∆Lk. Mean and 95% confidence intervals are shown. The data obtained using extrapolation of ∆x 2 are more robust compared to those obtained through extrapolation of mean extension z, but both methods agree withing experimental error. All data depicted here are obtained at force f = 0.5 pN.  Figure 4a) until the first reduction of the variance. There is a broad distribution of wait times, which is expected since both protein binding to DNA and the plectoneme fluctuations to form the specific loop required for binding are stochastic processes. The mean wait time for our SacII experiments is 93 s, which is much longer than for the other enzymes (NaeI and BsaXI) investigated, for which the wait times could not be determined accurately since binding often occurred during the introduction of the enzymes, i.e. before t = 0. A longer wait time for SacII compared to the other enzymes is also expected, since the long loop required for SacII binding is sampled less frequently than the shorter loops expected for NaeI and BsaXI.