{'fraction_non_alphanumeric': 0.03502916678759033, 'fraction_numerical': 0.01391589517369475, 'mean_word_length': 4.564392410173597, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 34, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Coarse-grained modelling of strong DNA bending II: Cyclization

DNA cyclization is a powerful technique to gain insight into the nature of DNA bending. The worm-like chain model provides a good description of small to moderate bending fluctuations, but some experiments on strongly-bent shorter molecules suggest enhanced flexibility over and above that expected from the worm-like chain. Here, we use a coarse-grained model of DNA to investigate the thermodynamics of DNA cyclization for molecules with less than 210 base pairs. As the molecules get shorter we find increasing deviations between our computed equilibrium j-factor and the worm-like chain predictions of Shimada and Yamakawa. These deviations are due to sharp kinking, first at nicks, and only subsequently in the body of the duplex. At the shortest lengths, substantial fraying at the ends of duplex domains is the dominant method of relaxation. We also estimate the dynamic j-factor measured in recent FRET experiments. We find that the dynamic j-factor is systematically larger than its equilibrium counterpart, with the deviation larger for shorter molecules, because not all the stress present in the fully cyclized state is present in the transition state. These observations are important for the interpretation of recent experiments, as only kinking within the body of the duplex is genuinely indicative of non-worm-like chain behaviour.


INTRODUCTION
As the mechanical properties of DNA play an important role in its biological capacities, there has been much activity to accurately characterize these properties, not only in the elastic regime of small to moderate fluctuations but also for more strongly stressed systems. For example, DNA is found to overstretch beyond a salt-dependent critical force [1]. Similarly, in response to twist DNA forms plectonemes beyond a critical buckling superhelical density [2]. Here, we are interested in the response of DNA to strong bending.
The worm-like chain (WLC) model provides a good description of small to moderate bending fluctuations in DNA [3][4][5][6]. However, although there is a consensus that for sufficiently strong bending the stress will be localized within small regions, often termed "kinks", much about this crossover to non-WLC behaviour remains controversial. For example, a recent review by Vologodskii et al. [7] highlighted a number of open questions, including what is the free energy cost of kink formation, how does the free-energy of a kink depend on bend angle, what is the critical curvature that causes the double helix to kink?
One way to address these questions is with molecular simulations of a coarse-grained DNA model, as such simulations are able to directly probe the relevant factors. In this and an accompanying paper [8] we have begun this task for the oxDNA model [9][10][11], which is able to describe well the thermodynamics of DNA hybridization and basic mechanical properties such as the persistence length and the torsional modulus [9] and which has been productively applied to study a wide variety of biophysical and nanotechnological systems [12]. In the first paper, we addressed the thermodynamics of kink formation and one particular experimental system, a "molecular vice" [13], that probed DNA strong bending. Here, we consider DNA cyclization, in particular looking at the length of DNA molecules at which duplex kinking begins to play a role in this process.

DNA Cyclization
DNA cyclization is a convenient model system used to probe DNA bending. Cyclization experiments were first reported in 1966, albeit on 48 500 base pair (bp) -DNA [14,15]. In 1981, Shore et al. developed a method to probe the bending of shorter 126-4361 bp fragments [16], later noting periodicity in the cyclization efficiency of 237-254 bp fragments [3].
A typical cyclization experiment, as depicted in Figure 1(a), uses a cyclization substrate with complementary sticky ends, N s bases in length, on both ends of a N d base-pair duplex. Cyclization leads to the formation of a N bp -base-pair duplex, where N bp = N s + N d . The resultant structure is not a closed minicircle -two back- bone 'nicks' are present at either end of the sticky ends. Either the forward rate or equilibrium constant of the cyclization reaction is reported. Experiments differ in how exactly they probe cyclization: methods based on ligation [3,5,16,17], FRET [19][20][21] and multimerization [36] have been reported.
In ligase-based experiments, cyclized molecules are fixed in the cyclized state by ligation of the two backbone nicks. The open and ligated cyclized molecules can then be resolved by gel electrophoresis and the concentration of different products measured. FRET-based experiments can be performed in equilibrium, with the molecules allowed to cyclize and uncyclize indefinitely. Fluorophores are attached to both ends of the molecule as FRET reporters: a high FRET signal will be reported when the duplex ends are in close proximity (cyclized), low FRET when apart (open). Although non-WLC behaviour has been suggested by the ligase-based experiments of C&W [17] and the FRET-based experiments of Vafabakhsh & Ha (V&H) [19], these results have been contested [5,7,37]. There is not yet a consensus on their interpretation.
Cyclization efficiency is traditionally reported in terms of a j -factor, first introduced by Jacobson & Stockmayer [38], which is a measure of the effective local concentration of duplex ends at zero end-to-end separation. The j -factor enables the ring closure probability to be calculated, and importantly, may be related to a ratio of equilibrium constants: where K cyc eq and K dim eq are the equilibrium constants for cyclization and dimerization, respectively.
Multimerization of the cyclization substrate yields a mixture of linear and circular products [5,17]. To avoid this complication, a separate dimerization substrate may be prepared with the same sequence as the cyclization substrate, but with only one N s -base complementary sticky end per molecule [19] (Figure 1 (b)). Following hybridization, the total length of the system is then 2N d + N s base pairs, with blunt as opposed to sticky ends. The consequences of this choice for K dim eq and hence j eq are discussed in Supplementary Section S1 C.
By assuming the contribution from base pairing to K cyc eq and K dim eq is the same, the WLC model can be used to estimate j eq . A common assumption is that the cyclized state is fully stacked, with coaxial stacking across the two nicks. For this situation, the analytic expression derived by Shimada & Yamakawa (SY) [39], which includes both the bending energy cost of bringing the two ends together and the twist energy cost of bringing the two helix ends into register, is the most appropriate.
Some of the confusion surrounding claims of cyclization efficiency greater than that predicted by the WLC model revolves around the use and interpretation of j -factors. While the j -factor relation using the ratio of cyclization to dimerization equilibrium constants is well established, experimental studies usually report the ratio of forward rate constants. In the case of ligase-based assays (reviewed in reference [40]): where k ligase cyc and k ligase dim are the forward rate constants for the formation of the ligated circle and dimer, respectively.
In the experimental limit where the ligation rate is very slow compared to the rate constants for uncyclization (k uncyc ) and undimerization (k undim ), the concentrations of unligated circles and dimers will reach an equilibrium with reactants. If this condition is met, j ligase dyn should be equivalent to j eq . In practice, this limit is valid for low ligase concentrations ([Ligase] ≪ [DNA]) and short sticky ends [5,26]. The importance of this condition is illustrated by Du et al. [5], who suggested that the apparent non-WLC behaviour in the C&W experiments was rather due to an insufficiently low ligase concentration (reviewed in reference [7]).
An additional complication with ligase experiments relates to the structure of the substrate, specifically the DNA surrounding the nick. It is unclear whether the ligase will act equally on all nicked duplexes, or only the subset that happen to adopt the right configuration at the nick, be that coaxially stacked or kinked, to allow the ligase to bind. As we show in this work, to alleviate stress, cyclized systems are more likely than dimers to break coaxial stacking at a nick. The ligation rate of hybridized complementary sticky ends may therefore vary substantially depending on the system.
The FRET experiments of V&H [19] have the advantage of directly monitoring the transition between cyclized and open states. While it is possible to report thermodynamics from FRET experiments, V&H report dynamics: where k cyc and k dim are the forward rate constants for the formation of the unligated circle and dimer, respectively. Dynamic j -factors extracted from FRET-based experiments must also be interpreted with care. Making a comparison between j FRET dyn and the SY prediction for j WLC eq , V&H make a claim of much greater than WLC flexibility at N bp 100 bp. However, j FRET dyn ≈ j eq only in the limit where k uncyc ≈ k undim , a condition that V&H, as well as another more recent FRET experiment [21], have shown not to be met. Given that k uncyc ̸ = k undim , one should not expect j FRET dyn ≈ j eq at short N bp . Thus, the observed deviation of j FRET dyn from j WLC eq is not necessarily an indication of non-WLC behaviour.
The interpretation of cyclization experiments is not straightforward; in particular, the microscopic states, responsible for the putative non-WLC flexibility, cannot be directly observed. Additionally, non-WLC behaviour has been reported in a number of other systems, including DNA minicircles [41,42], a "molecular vice" [13], and a "strained duplex" [43][44][45]. As we emphasised in the accompanying paper [8], establishing whether these observations are mutually consistent is a challenging task due to the distinct interplay of mechanics, geometry and topology inherent in each experimental system.
Fortunately, simulation can help bridge this gap. Models that reproduce the basic mechanical, structural and thermodynamic properties of DNA can be used to simulate the relevant systems. The results can be used to establish: (i) whether experimental results are truly indicative of non-WLC behaviour; (ii) whether the behaviour can be attributed to certain types of structure, such as 'kinks' [46][47][48]; (iii) whether different systems present consistent evidence for non-WLC behaviour; and (iv) whether the inferred occurrence of any disruptions to duplex structure is consistent with our current understanding of DNA thermodynamics.
Here we use oxDNA to probe the cyclization equilibrium as a function of duplex length. This approach provides direct access to microscopic states, enabling a structural interpretation of experimental observations. OxDNA is particularly well-suited to this task, as it provides a good description of both single-and double-stranded DNA, including hybridization thermodynamics, persistence length, torsional modulus and basic structure, and has previously been shown to reproduce a number of stress-induced transitions in DNA [49][50][51][52].

MATERIALS & METHODS oxDNA model
OxDNA [9][10][11] is a nucleotide-level coarse-grained model of DNA that has been employed successfully for a wide variety of systems [12], beginning with the thermodynamic and structural characterization of DNA nanotweezers [53]. Briefly, the model consists of rigid nucleotides with three interaction sites per nucleotide, interacting via Watson-Crick base pairing, base stacking, excluded volume, and a potential to represent backbone connectivity. The model is parameterized to reproduce the thermodynamics of duplex melting at high-salt ([Na + ]− −500 mm), where backbone-backbone electrostatic repulsion is shortranged due to counter-ion screening. Two parameterizations of the model are available. We use the the averagebase parameterization [9,10], in which the strength of the base-pairing and stacking interactions are independent of the identity of the bases, to highlight the basic thermodynamics of DNA cyclization. By contrast, the sequence-dependent parameterization [11] has stacking and base-pairing interactions that depend on the base identity and is used to compare more directly to the V&H experiments.

Simulations
Simulations of the cyclization equilibrium were performed with a virtual-move Monte Carlo (VMMC) algorithm [54] at 298 K. As the free-energy barrier between typical open and cyclized states is large, the transition between the two macrostates constitutes a rare-event. Umbrella sampling, a technique that allows for the biased sampling of states with respect to an order parameter [55], was employed to sample the barrier crossing in reasonable computational time.
We use a two-dimensional order parameter Q = (Q ee , Q bp ) to characterize the transition. Q ee is a discretized measure of the distance of closest approach between the complementary sticky ends. Q bp is the number of base pairs formed between complementary sticky ends, where 0 ≤ Q bp ≤ N s .
To further improve computational efficiency, umbrella sampling was windowed to separately sample the open and cyclized states of each molecule. For the window associated with the open state, the system was restricted to Q bp = 0; for the window associated with the cyclized state, the system was restricted to Q ee = Q min ee (the value corresponding to the shortest distances between sticky ends). Simulations were run until convergence to within ±5 % for each window.
The sampling windows overlap at Q = (Q min ee , Q 0 bp ); results were combined by normalizing each window so that the free energies were equal for this value of the order parameter. As there is only one well-defined overlap between the values of the order parameters for both windows, more complex approaches, such as the weighted histogram analysis method [56,57], were unnecessary. To further simplify sampling, we forbade the formation of base pairs that are not intended in the design of the system (non-native base pairs). Further details of the cyclization simulations can be found in Supplementary Section S1 A. The simulation of dimerization equilibrium is roughly analogous to cyclization, and is elaborated in Supplementary Section S1 B. To compute the equilibrium constants, we deemed all states with Q bp ≥ 1 to contribute to the cyclized and dimerized states (details in Section S1 C).
A complete list of sequences is available in Table S1. Error bars represent the standard error of the mean from 5 independent simulations.

Structural analysis
For sufficiently strong bending stress, localized structural disruptions to the DNA double-helical structure are expected. We define three such disruptions, namely fraying, bubble formation and kinking, which are elaborated in detail in the accompanying paper [8]. Briefly, both fraying and bubble formation involve the breaking of base pairs; the difference is in the location along the duplex. Fraying involves disruption of base pairing at the duplex ends, while bubbles occur in otherwise fully base-paired contiguous stretches away from the duplex ends.
Conceptually, a kink is an area of strong bending localized to a small segment of DNA, and can occur both at a nick and within the duplex. When kinking occurs in a duplex region, it is nearly always accompanied by bubble formation. Similarly, kinking at a nick can be accompanied by fraying. We do not attempt to distinguish different types of kinks, as has been done when analysing atomistic molecular dynamics simulations [47,58]. Instead, we simply define a kink as present within a duplex region when there is a change in orientation of consecutive bases of greater than 90 ∘ on either strand. For nicked regions, only consecutive bases on the intact (unnicked) strand are considered. Although this cutoff is somewhat arbitrary, since in the current system, kinks arise to localize bending stress, they are usually very sharply bent and the criterion works well [8]. Sometimes, however, it gives rise to false negatives, particularly in the case of kinks at nick sites, and when fraying is present. Nonetheless, it is satisfactory as an indicator of behaviour for our purposes. More details, along with subtleties related to kinking at a nick, are discussed in Supplementary Section S1 D.

RESULTS
We simulate a large range of system sizes; some illustrative configurations are shown in Figure 2. Quantitatively, we first consider the behaviour of j eq as a function of length. For the dimerization system we only computed the equilibrium constant K dim eq for a few lengths (monomers N bp = 30, 67, 73, 101). As expected, we found K dim eq to be length-independent to within numerical error (Supplementary Section S2 A, Table S2); therefore, we use an average value of K dim eq = 0.92 ± 0.20 × 10 12 m −1 in our j -factor calculation. In contrast, we found K cyc eq , and thereby j eq (Eq. (1)), to vary substantially with length.
In Figure 3, we show j oxDNA eq values, calculated from our measured values of K cyc eq and K dim eq using Eq. (1), for 81 different lengths in the range N bp = 30 − 207 for fixed N s = 10 using the average-base parameterization of oxDNA. These results are compared to j WLC eq predictions based on the Shimada & Yamakawa (SY) expression [39] using previously calculated values for the relevant structural and mechanical properties of oxDNA [9]. The SY expression is appropriate for the formation of a fully stacked "circle" configuration with coaxial stacking at both nicks (Figure 2 (a)). Note that the comparison in Figure 3 is fit-free. The effect of varying N s as well as the role of nicks and mismatches are discussed in Supplementary Section S2 C.
The behaviour of the SY expression is well understood. In the regime of interest, shortening N bp tends to make cyclization less favourable as a result of increased bending stress within the duplex. This effect becomes particularly acute for DNA lengths below the persistence length (41.82 nm/126 bp for the curve plotted in Figure 3). On top of this systematic behaviour, a periodic oscillation is associated with the need to over-or under-twist the duplex when the natural twist is not commensurate with that required to form a closed circle. The magnitude of this oscillation increases at shorter N bp because a stronger twist per base pair is required.
When considering the behaviour of j oxDNA eq in light of the SY expression, three length-scale dependent regimes become apparent: long (N bp 80), intermediate (N bp ≈ 45 − 80 bp) and short (N bp 45).
In the long length regime (N bp 80), oxDNA reproduces the periodic oscillations predicted by the SY expression, and values of j eq coincide at the maxima of these oscillations. However, even at the longest lengths we consider, the magnitude of the oscillation is smaller for oxDNA than predicted by the SY expression. At shorter N bp , the magnitude of the oscillation decreases for oxDNA, at odds with the SY expression. The difference stems from the possibility of adopting alternative "teardrop" configurations ( Figure 2 (b)), in which most of the twisting stress, and some of the bending stress, can be relieved by kinking at one of the nicks, thus breaking coaxial stacking. In these configurations there is still the constraint of DNA continuity at the nick and this can now be more easily satisfied, rather than through overor undertwist, by out-of-plane bending (Figure 2 (d)).
The possibility of adopting this alternative teardrop configuration reduces the free-energy penalty for incommensurate values of N bp (i.e. (n + 1 /2) × pitch length), thus suppressing the oscillations. Further, at shorter N bp , the bending stress increases more slowly for these kinked teardrop configurations than for the coaxially stacked circles. Thus, the difference in j eq between the "on-register" (coaxially stacked) and "off-register" (kinked at a nick) molecules decreases at shorter N bp , rather than increasing as predicted by the SY model. Although not taken into account for most analyses of cyclization, the possibility of the cyclized molecule exhibiting a "teardrop" configuration has previously been suggested by Vologodskii et Configurations with a kink in the duplex and either (e) a kink at one of the nicks or (f ) kinks at both nicks. For short duplexes where N d is not that much larger than Ns, the sticky ends can associate either by (g) relatively minor bending of the duplex or (h) fraying a few base pairs.
Further evidence in support of this analysis is given in Figure 4 (a), which shows the probability of kinking as a function of N bp . At long lengths (N bp 80), kinking does not occur in the duplex regions, but can occur at the nick sites. There is clear periodicity in kinking at a nick. For example, at N bp = 145 ≈ 14 × pitch length, the probability of kinking at either of the nicks is negligible and the system virtually always adopts a coaxially stacked circle configuration; however, at the longer N bp = 201 ≈ 19.5 × pitch length, the probability is ∼ 40 %. As N bp is shortened, the probability of kinking at a nick gradually increases for these "off-register" lengths; it is not until N bp = 114 ≈ 11 × pitch length that the bending stress along the duplex is sufficient to cause ∼ 10 % of molecules to kink at a nick for an "on-register" length.
In the intermediate-length regime (N bp ≈ 45 − 80 bp), we observe enhanced cyclization efficiency compared to the SY prediction. Although j oxDNA eq continues to decrease with N bp , it does so more gradually than the SY expression would predict; consequently, j oxDNA eq is in excess of the peak envelope of the SY expression. In this regime, bending stress in the circular coaxially stacked configuration is sufficiently large that kinking at one of the nicks occurs for even the on-register systems (Figure 4 (a)).
Oscillations due to on-register effects also seem to contribute to j oxDNA eq in the upper end of the intermediatelength regime, with shallow maxima occurring at N bp = 63 and N bp = 74, approximately 6 and 7 times the pitch length. However, these sizes no longer correspond to minima in the probability of kinking at a nick ( Figure 4 (a)) so the structural underpinnings of these variations is less clear. For shorter lengths, although there are sizedependent variations in j oxDNA eq (e.g. maxima at N bp = 43 and 49), there is no longer a simple relationship to the pitch length, instead reflecting more complex geometric compatibilities that allow cyclized states at these lengths to be particularly stable compared to nearby lengths.
In the intermediate regime, most cyclized molecules are kinked at one of the two nicks. At shorter N bp , the bending stress in the duplex region of the teardrop configurations increases, with the highest curvature being localized opposite the nick that is kinked. Consequently, it becomes   [39] (grey solid line) is plotted using values of torsional stiffness (4.75 × 10 −28 J m −1 ), persistence length (41.82 nm) and pitch length (10.36 bp/turn) appropriate to oxDNA [9,51]. The dashed line gives the maxima envelope for the SY prediction (grey dashed line). Sequence-variation in j oxDNA eq , computed using oxDNA's sequence-dependent parametrization, is illustrated for sequences considered by V&H (brown triangles), including six at N bp = 73 bp, as well as poly-AT and poly-GC.
increasingly favourable to localize bending stress into a kink in the duplex (unnicked) region, with the probability of this duplex kinking (Figure 2 (e)) increasing from near zero at N bp = 81 to near one at N bp = 43 ( Figure 4 (a)). A kink in the duplex will generally be located opposite a kink at a nick because this arrangement minimizes the residual bending stress in the unkinked portions of the duplex by equalizing the lengths of the double-helical segments between the two kinks. In addition to a loss of stacking, kinking in the duplex typically involves breaking 1-2 base pairs [8]. A typical configuration for such a kinked duplex state is shown in Figure 2 (e). As with kinking at the nick, the availability of this configuration lowers the free-energy cost of cyclization, raising j oxDNA eq further above the SY prediction.
At the lower end of the intermediate-length regime, the cost of bending without kinking in the duplex region is so high that molecules with duplex kinking dominate the cyclized state. States with two kinks localize the majority of the bending stress at the kinking sites, as shown in Figure 2 (e), with long duplex sections relatively relaxed.
Consequently, the free-energy cost of looping is largely independent of N bp in this regime, causing j oxDNA eq to level off and reach an approximately constant value. In contrast, the SY j WLC eq prediction decreases very rapidly. For example, at N bp = 43, j oxDNA eq is 10 9 times greater than the SY prediction.
There exists a rich landscape of structures informed by slight differences in local geometries as a function of N bp . In addition to the canonical two-kink structures, containing one kink in the duplex and one kink at a nick ( Figure 2 (e)), we observe several non-trivial arrangements, albeit with relatively low probability. For example, the configuration in Figure 2 (f) contains a kink at both nicks as well as in the duplex.
In the short-length regime (N bp ≤ 42), we observe an increase in j oxDNA eq , in stark contrast to the rapidly decreasing j WLC eq predicted by the SY expression. Given the many WLC assumptions that are violated at this length scale, a deviation is unsurprising; however, we did not anticipate an increase in j oxDNA eq . As N d is now not that much larger than N s , kinking at both nicks allows the single-stranded sticky ends to hybridize without duplex kinking. The system now generally adopts a conformation of two parallel duplexes, with the stress now borne by a mixture of continuous bending (Figure 2 (g)) and fraying of a few base pairs at the ends of the duplexes (Figure 2 (h)). Figure 4 (a) shows a very clear crossover to cyclized states with kinks at both nicks and no kinks present in the duplex, occurring abruptly at N bp = 40 − 42 . At the same time, fraying at the duplex ends also increases (Figure 4 (b)). In this regime, as N bp is shortened, the difference in length of the duplexes decreases and the stress in the system tends to drop, leading to higher j eq (Figure 3). Overlaid on this overall trend are non-trivial geometric effects associated with whether the lengths of the two duplexes relative to the pitch length is convenient for connecting them (the thickness of the double helix is now significant compared to the duplex lengths), which leads to non-monotonic behaviour of j eq and kinking. Note that the value of N bp at the crossover between the short-and intermediate-length regimes is expected to be very dependent on N s , occurring at smaller N bp for smaller N s . So far we have only reported results using the oxDNA average-base parameterization, in which the strength of base pairing and stacking interactions are independent of base identity. As the free-energy cost of disrupting a duplex to form a kink is sequence-dependent, it is important to consider how sequence might perturb the general trends we have elaborated thus far. We therefore studied a variety of sequences used in the V&H [19] experiments including six at N bp = 73 using the parameterization of oxDNA that includes sequence-dependent thermodynamics ( Figure 3).
As duplex kinking for oxDNA typically involves the breaking of base pairs, we observe that this kinking preferentially occurs at A-T base pairs. Base pairs that are weaker than average introduce preferred locations for kinking; thus, when duplex kinking is relevant, j oxDNA eq is expected to be larger for the sequence-dependent than for the average-base parameterization. This is indeed the case for N bp = 73. The sequence-induced variation for the V&H N bp = 73 sequences is a factor of ∼ 4 (GC-content 11-52 %). This compares to a factor of ∼ 8 between the extrema in GC-content, poly(AT) and poly(GC).
That we find sequence heterogeneity generally makes duplex kinking easier is consistent with our explicit investigation of the free energy of duplex kinking in the accompanying paper [8]. Our results imply that the crossover to cyclized configurations with duplex kinking occurs at slightly longer N bp when sequence-dependence is included than for the average-base parameterization.
As well as equilibrium constants, free-energy profiles as a function of the number of bases pairs were computed for each system we considered, as these can give more insight into the pathway for association. Example profiles for cyclization and dimerization at N bp = 101 are illustrated in Figure 5. One interesting feature of the profiles is that the free-energy gain from hybridizing the complementary sticky ends once an initial base pair has formed is less for cyclization than for dimerization: ΔG ‡ uncyc < ΔG ‡ undim . Physically, this indicates that substantial additional bending stress develops as subsequent base pairs form for the cyclization system, reducing the free-energy gain upon zippering of the sticky ends relative to dimerization. For example, the Q bp = 1 configuration in Figure 2 (c) is clearly less bent than the Q bp = 10 configuration in Figure 2 The activation energies derivable from these profiles will be particularly useful in the next section when we consider the dynamic j -factor.

Comparison with experiment
As noted in the Introduction, the results of ligase experiments performed in the low ligase concentration limit should, in principle, be comparable to equilibrium j -factors (although there are subtleties related to the ensemble of states actually detected by the ligation enzymes). Indeed, for long DNA molecules (N bp much longer than Free energy profiles of cyclization (solid) and dimerization (dashed) for a N bp = 101 system (Ns = 10, N d = 91). The activation free-energy barriers, ΔG ‡ for the forward (cyclization, dimerization) and reverse (uncyclization, undimerization) reactions are labelled. ΔG ‡ cyc reflects the free-energy cost of bending to form the first base pair in a cyclization system, whereas ΔG ‡ dim reflects the entropic cost of bringing two monomers together within the simulation volume. The dimerization simulations are for a cubic simulation box of dimension 85.18 nm, corresponding to a duplex concentration of 2.69 µm. Note that for clarity, we have depicted ΔG ‡ uncyc and ΔG ‡ undim as the free energy difference between the fully base-paired closed state Q bp = 10 and the transition state Q bp = 1; however, in practice, frayed states Q bp = [2,9] do contribute to the closed state. While the distinction does not significantly impact our results (< 1 /4 k B T ), we do include the contribution of frayed states in both ΔG ‡ uncyc and ΔG ‡ undim .
the persistence length) and low ligase concentrations, there has been consistent agreement between experiment and the WLC model [3,5,16,17,19,21]. However, Cloutier & Widom (C&W) [17] reported results for N bp shorter than the persistence length (N bp = 93, 94, 95, 105, 116), showing an apparent deviation from WLC behaviour, with j ligase dyn (Eq. (2)) enhanced over the SY WLC prediction j WLC eq [39] by a factor of 10 2 − 10 4 . The differences between the maxima in j oxDNA eq and j WLC eq in this size range are much smaller ( Figure 3) and so cannot account for this discrepancy.
In contrast, Du et al. [5] found no deviation from WLC behaviour for N bp = 105−130 . Furthermore, they presented evidence suggesting that the C&W experiments used too high a ligase concentration to enable j ligase dyn to be compared with j WLC eq . The results of Du et al. are in good agreement with the SY WLC expression, albeit with somewhat different materials parameters (torsional stiffness, persistence length and pitch length) than for oxDNA, owing in part to different buffer conditions (Supplementary Section S2 D, Figure S4 (a)).
FRET measurements on DNA cyclization, as pioneered by Vafabakhsh & Ha (V&H) [19], provide a more direct measure of cyclization because k cyc , k uncyc and K cyc eq are obtainable, although V&H mostly report k cyc . V&H claim enhanced cyclization at N bp 100, based on a comparison between their j FRET dyn (Eq. (3)) and the SY WLC expression j WLC eq . However, as noted earlier, this is only a fair comparison if k uncyc = k undim , which, as k undim is expected to be length independent, also implies that k uncyc should be independent of N bp . However, since both V&H [19] and more recent FRET measurements [21] suggest that k uncyc increases with N bp , this condition is not met.
Although we do not directly simulate the dynamics of cyclization, we can estimate the relative rates of processes using free-energy profiles such as those in Figure 5 and unimolecular rate theory for activated processes [59]. In agreement with experimental investigations, previous work on oxDNA has shown that duplex formation has an effective "transition state" involving a very small number of base pairs [60]. We therefore make the assumption that uncyclization and undimerization rates are given by: where A is a constant for DNA melting, and ΔG ‡ uncyc and ΔG ‡ undim are defined in Figure 5. The physical content of this assumption is that an increased favourability of base-pair formation is manifested in slower unbinding rates. This is important because bending stress in cyclized systems reduces ΔG ‡ uncyc ( Figure 5). The rate constants will be equal (k uncyc = k undim ) in the very long length limit (N bp much longer than the persistence length).
ΔG ‡ cyc and ΔG ‡ uncyc are plotted as a function of N bp in Figure 6. In particular there is a general decrease in ΔG ‡ uncyc at shorter lengths, suggesting k uncyc increases with shorter N bp , in agreement with experimental results [19,21]. Additionally, any torsional stress in the cyclized state is relieved as base pairs are disrupted, leading to the oscillations in ΔG ‡ uncyc at long lengths, with the minima occurring at the more torsionally stressed offregister lengths.
To compare with V&H's results for j FRET dyn , we note that Therefore, using Eq.
at each N bp .
We are now in a position to estimate j oxDNA dyn from oxDNA's equilibrium constants and activation free-energy barriers. We expect K dim eq and ΔG ‡ undim to be lengthindependent as the excluded volume of the duplex far from the complementary single-stranded sticky ends is likely to have little effect on the dimerization process. Indeed, this appears to be the case to within less than 0.5 k B T in ΔG ( Figure 6(b), Supplementary Table S2). Our resulting j oxDNA dyn values are plotted in Figure 7, where they are compared to our j oxDNA eq values, and V&H's j FRET dyn . We observe that oxDNA's j oxDNA dyn values lie substantially above j oxDNA eq , except in the limit of long N bp for onregister lengths. Only if all the stress-induced destabilization of the fully cyclized state is exhibited in the forward rate would j dyn = j eq . Although the activation free-energy for cyclization increases as N bp decreases (Figure 6(a)), the full bending and twisting stress present in the cyclized state is not yet present in the transition state where the two complementary sticky ends have formed their first base pair, e.g. Figure 2(c), and so j oxDNA dyn > j oxDNA eq . Instead, stress which is not present in the transition state leads to a decrease in the free-energy barrier for uncyclization ( Figure 6(b)) and accelerates uncyclization relative to undimerization. Alternatively, one can consider that the stress in the cyclized state subjects the duplex formed between the sticky ends to a shear force [21] which is well known to lead to a more rapid rupture of a duplex [52,61].
The observed behaviour is similar to that seen by V&H; indeed, j oxDNA dyn and j FRET dyn agree remarkably well in the range N bp = 90 − 105 , with reasonable agreement extending to N bp = 70. Our results suggest that variation of uncyclization rates with N bp may be an important contribution to apparent non-WLC behaviour in j FRET dyn . To begin to understand the difference between V&H's dynamic j -factor and j WLC eq a number of authors have tried to account for the role of the single-stranded tails in cyclization by incorporating a "capture radius" into the WLC j -factor calculation, as an approximation for how close the duplex ends must be in order for the sticky ends to hybridize [7,19,21,37]. Indeed, this approach leads to a significantly enhanced j -factor when a value that is taken to be roughly appropriate for the 10-base sticky ends of the V&H experiment, namely 5 nm, is used [19,37].
Assuming this approach is an attempt to capture the weaker constraints at the transition state, and hence to estimate the activation free-energy barrier relevant to j dyn (note j dyn and j eq are often not clearly differentiated in discussions of V&H's results), we can use our oxDNA results to test the reasonableness of this approximation by measuring the separation of the duplex ends in the transition state ensemble. At long lengths our measured capture radius is approximately constant with a value just under 4 nm, but increases at shorter lengths because some of the stress at the transition state is partitioned into stretching the single-stranded tails, reaching a maximum of about 7 nm at about N bp ≈ 40 ( Supplementary Figure S5). Thus, although the 5 nm value used previously [19,37] is not unreasonable, unsurprisingly this approach does not capture the full complexity of the transition state to cyclization, and nor does it account for the strain present in the sticky ends for them to achieve contact.
Therefore, although explaining a dynamic j -factor in terms of a capture radius is physically well-motivated, it does not provide a full understanding. Moreover, as its value is not tightly constrained, and certain contributions to stability are neglected, it is difficult to judge whether WLC behaviour is violated or not by comparing a calculated j -factor curve to an experimental j dyn .
We note that j oxDNA dyn , in contrast to j oxDNA eq , varies comparatively smoothly with N bp , with only a weak periodicity on the length scale of the pitch length. The shallow maxima at larger lengths occur at N d = (n + 1 /2) × pitch length, because the two sticky ends are then on the same side of a torsionally unstressed duplex. Consistent with this smooth variation, ΔG ‡ cyc also varies relatively smoothly with N bp . The strong periodicity in  (red squares) [19]. Results for the oxDNA sequence-dependent parameterization using V&H 14 variable N bp sequences, in addition to 6 sequences at N bp = 73 to highlight the role of sequence-variation. For reference, j oxDNA eq and j WLC eq are also shown (grey). j oxDNA eq comes from the free-energy gain when zippering up the complementary sticky ends (i.e. ΔG ‡ uncyc ), which is greater when N bp is an integer multiple of the pitch length, allowing the formation of a relatively relaxed coaxially stacked circle. V&H suggest that their j FRET dyn data at N bp = 93 − 106 displays a strong oscillation with a period of about one pitch length (Figure 7). In agreement with Vologodskii et al., [7] we find no physical mechanism for such a strong oscillation and would suggest that the experimental evidence for this oscillation is not compelling.
At the shortest lengths investigated, j FRET dyn still lies above j oxDNA dyn . Some of this difference might be due to our use of the oxDNA average-base parameterization. Sequence does play a role in DNA flexibility, as noted in the accompanying paper [8], but the impact of sequence variation is non-trivial. In particular, because kinks tend to localize to AT base pairs, kinking in the duplex is easier for a sequence with 50 % GC-content than reported for our average-base parametrization.
Using the oxDNA parameterization with sequencedependent thermodynamics we found j dyn to increase compared to the values for the average-base parameterization. Consistent with our results for j oxDNA eq , we find a sequence-induced variation in j oxDNA dyn of a factor of ∼ 4 for the six V&H N bp = 73 sequences (GC-content 11-52 %) [19] compared to a factor of ∼ 8 between the extrema in GC-content, poly(AT) and poly(GC). This compares to a factor of ∼ 60 for experimental looping rates, a discrepancy which may be explained by V&H's use of poly(A) tracts, a sequence-motif well-known to introduce intrinsic curvature in duplex DNA [62]. As oxDNA's sequence-dependent parameterization is based on the nearest-neighbour thermodynamics of SantaLucia et al. [63,64], alternative structural motifs such as poly(A) tracts are outside the scope of the model. OxDNA also does not reproduce sequence-dependent structural (e.g. the difference in size between purine and pyrimidine) or mechanical (e.g. flexibility) properties.
The origin of the remaining discrepancy between j FRET dyn and j oxDNA dyn at the shortest N bp is not yet clear, but may indicate enhanced flexibility for V&H versus oxDNA due to WLC (i.e. lower persistence length) or non-WLC behaviour (i.e. kinking within the duplex at slightly longer N bp ). The oxDNA persistence length of 41.82 nm is within the range of experimental observations at [Na + ]− −500 mm, but smaller values at [Na + ]− −750 mm are not implausible [65,66]. It is also possible that oxDNA slightly underestimates the prevalence of kinking within duplex regions [67]; an onset of kinking at slightly lower stress (longer N bp ) would make cyclization at shorter lengths more favourable. Finally, it is worth noting that the presence of fluorophores may cause perturbations in the V&H experi-ments.
V&H do report K cyc eq for some systems; oxDNA results are in good agreement for those lengths (Supplementary Figure S2). CONCLUSION Cyclization is a system dependent manifestation of the general thermodynamics of strong DNA bending, elaborated in the accompanying paper [8]. The remarkable range of behaviour in cyclized systems is explicable by the interplay between three specific deformation modes of stressed duplexes: continuous bending, kinking and fraying.
OxDNA reveals that each of these modes is present at a characteristic length-scale with respect to cyclization: continuous bending at long lengths (N bp 80), duplex kinking at intermediate lengths (N bp ≈ 45 − 80 bp) and fraying at short lengths (N bp 45 bp). In addition, as N bp is shortened, there is an increase in kinking at the two nicks that remain after the hybridization of the sticky ends. At longer lengths, kinking at a nick is only observed for "off-register" molecules that cannot form torsionally relaxed coaxially stacked circles. The ability of said nicks to relax bending, as well as torsional stress, means that they become increasingly prevalent for shorter N bp . At the shortest lengths, kinking at both nicks is dominant.
We use oxDNA to probe the reported observation of non-WLC behaviour in FRET-based cyclization experiments [19]. In agreement with experiment, we observe that for shorter values of N bp , the apparent j -factor lies substantially above the predictions of the Shimada & Yamakawa (SY) [39] WLC model. We also observe that the periodic oscillations predicted by the SY model are suppressed. This behaviour arises from two conceptually distinct phenomena.
Firstly, highly stressed cyclized systems can adopt configurations that relax stress more effectively than through continuous bending, thereby reducing the overall freeenergy cost of cyclization relative to a direct estimate based on a simple WLC-based model. At various values of N bp , oxDNA identifies kinking at nicks, kinking within the duplex region and fraying of base pairs as key relaxation modes.
Secondly, oxDNA suggests that not all of the reduction in K cyc eq relative to K dim eq is due to stress manifest in the cyclization rate; uncyclization rates are also substantially increased relative to undimerization rates. The result is that dynamic j -factors based on the ratio of cyclization and dimerization rates lie even further above the SY prediction than their equilibrium j -factor counterpart.
Of the above effects, only kinking within the duplex can reasonably be described as truly non-WLC behaviour. WLC and related statistical models do not predict absolute rates directly. Kinking at nicks and fraying can only occur when the DNA backbone is discontinuous, and any resultant effects are unrelated to whether WLC mod-els accurately describes the body of the DNA duplex. In oxDNA (with the average-base parameterization), kinking in the duplex only has a substantial effect for N bp 70 bp. Note, if sequence dependence is taken into account, then kinking may occur for slightly longer lengths.
Our results suggest that much of the apparent "extreme bendability" reported by Vafabakhsh and Ha [19] can be attributed to factors that are not strictly speaking evidence of non-WLC behaviour. We cannot account for deviations at their very shortest lengths (N bp 70). In oxDNA, kinking within the duplex region is present, but not completely dominant. It is possible that oxDNA slightly overestimates the difficulty of kinking within a duplex -if this is the case, the data for the very smallest values of N bp studied by V&H may be indicative of duplex flexibility over and above that predicted by the WLC, although we note that to give a substantial effect on j eq , kinking must not only be present, but must dominate the ensemble. Kinking dominating the ensemble at N bp ≈ 70 is inconsistent with oxDNA predictions, and available experimental evidence for a "molecular vice" [8] and DNA minicircles [67].
Various authors have incorporated a "capture radius" into the WLC j -factor calculation to capture phenomenologically some of the effects listed above [7,19,21]. Indeed, this approach leads to an enhanced j -factor, but the choice of capture radius is somewhat arbitrary and imprecise, making it difficult to assess whether non-WLC behaviour is present. Nonetheless, we find that the 5 nm capture radius that has been used in the interpretation of V&H's experiments is not unreasonable.
OxDNA is only a model, and good correspondence with experimental results should not be over-interpreted. Nonetheless, the relaxation mechanisms identified are clearly physically plausible. For example, enhanced uncyclization rates have previously been noted in the literature [19,21] (and also for transcription-factor mediated looping [68]). It is clear that much of the apparent discrepancy between the data of V&H and the predictions of WLC-based models are due to effects that are not true violations of the WLC model of duplex DNA flexibility. Our study also helps to reconcile the results of V&H with previous ligasebased assays which saw no evidence of enhanced flexibility at N bp ≈ 100 [5], and previous studies of minicircles which detected no evidence of duplex disruption at these length scales [41].
To explore whether the shortest lengths studied by V&H do show evidence of kinking and enhanced flexibility, we would propose experiments of shorter sequences and systematic collection of both dynamic and equilibrium data. The latter is extremely important; statistical WLC models make equilibrium predictions, and so the breakdown of a WLC description can only be confirmed with equilibrium data. Indeed, elucidating the subtleties of dynamic and quasi-dynamic (C&W) j -factors is one of the key issues addressed in this work. Supplementary material for "Coarse-grained modelling of strong DNA bending II: Cyclization"

A. Cyclization simulations
Cyclization simulations were performed in three phases: exploratory, equilibration and production. In all cases, we use a virtual-move Monte Carlo (VMMC) algorithm [S1] in combination with umbrella sampling [S2].
In the exploratory phase, we iteratively adjusted the umbrella sampling bias for the windows associated with the open and cyclized states, yielding a flat population of states for both windows. We use discrete potentials for both dimensions of the order parameter (see main text). For Q ee , we use variable width distance increments: 0-1.7036 nm, 3.4072, 5.1108, 8.518, 12.777, 17.036, 21.295, 25.554, 34.072, 42.59, 51.108 and > 51.108 nm. For Q bp , the increment is fixed, and is simply the number of base pairs formed (0 ≤ Q bp ≤ N s ). As the weighting iterations were semi-automated, simulation times varied in the range 10 6 -10 7 VMMC steps per particle. Simulations were performed with one molecule in a cubic box of dimension 170.36 nm, which corresponds to a unimolecular concentration of 336 nm. To further simplify sampling, we forbid the formation of base pairs that are not intended in the design of the system (non-native base pairs).
In the equilibration phase, we equilibrated the system for 10 7 VMMC steps per particle using the aforementioned umbrella weights. Adequate sampling (number of transitions in Q ee and Q bp ) and decorrelation (via a block average decorrelation method) of potential energy, bubble size, fraying and structural kinking were checked. For the production phase, each simulation (five independent simulations per measurement) was initialized with a statistically independent starting configuration and a unique random seed. This is accomplished by randomly drawing starting configurations from the equilibration phase, one for each production trial, ignoring the first 2 × 10 6 VMMC steps per particle of equilibration. The production trial run time is 10 7 VMMC steps per particle. For reference, the characteristic decorrelation time for the potential energy is ∼ 10 4 VMMC steps per particle, while the decorrelation time for kinking is ∼ 10 5 VMMC steps per particle.
The "seed moves" used to build clusters in the VMMC algorithm were: • Rotation of a nucleotide about its backbone site, with an axis chosen uniformly on the unit sphere, and with an angle drawn from a normal distribution with a mean of zero and a standard deviation of 0.10 radians.
• Translation of a nucleotide, where the displacement along each Cartesian axis is drawn from a normal distribution with a mean of zero and a standard deviation of 0.085 18 nm.

B. Dimerization simulations
Dimerization simulations follow a very similar procedure as cyclization simulations (Section S1 A), but adapted for a bimolecular system. With the exception of the starting configuration (bimolecular), exploratory, equilibration and production procedures are identical. As discussed in the main text, each molecule has one complementary sticky end and one blunt end. This precludes the formation of multimers and circular dimers, allowing only the formation of linear dimers. Simulations were performed with these two complementary, non-palindromic, stickyended duplexes, in a cubic box of dimension 170.36 nm, which corresponds to a dimer concentration of 336 nm. To estimate the impact of excluded volume effects, dimerization systems were also performed in a box of half the size and therefore 8 times the concentration (85.18 nm per side, which corresponds to a dimer concentration of 2.69 µm). As for the cyclization simulations, non-native base pairing was disallowed.
Analogously to the cyclization simulations, the dimerization simulations are windowed using the order parameters Q ee and Q bp , respectively the distance of closest approach and the number of base pairs formed between complementary sticky ends. The window associated with the dimerized state is given by Q bp ≥ 1 bp (Q ee = Q min ee , the small separation between base sites of the base paired nucleotides), while the window associated with the undimerized state is given by Q bp = 0 bp.

C. Computation of equilibrium constants
The cyclization reaction is unimolecular (Figure 1 (a)): where A and B are the open and cyclized states respectively, k cyc is the forward rate constant (cyclization) and k uncyc is the reverse rate constant (uncyclization). The equilibrium constant for cyclization, K cyc eq , can be estimated directly from simulations of a unimolecular, isolated system as Here P A and P B are the probabilities with which uncyclized and cyclized systems are observed in simulation (P A + P B = 1).
The dimerization reaction is a bimolecular association of distinct molecules: In this case, the bimolecular equilibrium constant can be inferred from a simulation of a single pair of dimerizing monomers via in which [A 0 ] is the total concentration of strand type A in simulation, and P A/B and P AB are the probabilities with which monomers and dimers are observed in simulation. This expression follows from Eq. 7 of Ref. [S4]. The resultant K dim eq obeys the standard relation for bulk systems in equilibrium, Recall from the main text that we use a definition for the equilibrium j -factor j eq that is consistent with Vafabakhsh and Ha (V&H) [S3]: This is a measure of the effective concentration of one sticky end in the vicinity of the other when the dimerization reaction involves distinct monomers, forming a heterodimer, with only one sticky end per monomer. Researchers often estimate K dim eq using the dimerization of identical monomers, forming a homodimer, with palindromic sticky ends. In this case, given the same underlying interaction strength between sticky ends (and hence the same K cyc eq ), K dim eq would be increased by a factor of two due to combinatorial (4 times the number of possible dimers) and symmetry effects (a factor of 1/2 for a homodimer rather than a heterodimer) [S4, S5]. As a result, with this approach the equilibrium j -factor is estimated as In the case of identical monomers with two nonpalindromic sticky ends, there is no pre-factor in the j -factor estimate, as the combinatorial (twice the number of possible dimers) and symmetry effects (a factor of 1/2 for a homodimer) cancel. Overall, under the approximation of ideal behaviour of separate complexes, these three approaches are equivalent.

D. Structural criterion for kink detection
While it is often visually straightforward to identify kinks, automating their detection is not without difficulty.
We have developed both energetic and structural criteria for identifying kinks, which are based on disruption of stacking interactions and changes in base orientation, respectively. In this study, we rely on the structural criterion; we have discussed the differences between the two criteria in detail in the supplementary material of the accompanying paper [S6].
Structurally, we define a kink using the relative orientation of consecutive nucleotides. In oxDNA, the orientation of each nucleotide is unequivocally determined by two orthogonal unit vectors, the base-backbone vector and the base-normal vector. The base-backbone vector connects the backbone and base interaction sites of each nucleotide. The model is designed such that, in a relaxed duplex, the base-normal vector at index i ,â i , and at the consecutive index i + 1,â i+1 , are approximately parallel; that is,â i ·â i+1 ≈ 1. A kink is defined to be present ifâ i ·â i+1 < 0, a condition implying a more than 90 ∘ change in orientation of consecutive nucleotides along one strand. Naturally, for duplex regions, we consider kinking along either strand; although, usually if a kink is present in the duplex, then both strands are kinked. To limit false positives due to fraying, we do not include the first and last 3 pairs of nucleotides in the duplex region in our analysis.
In the fully cyclized configuration there are two "nicks" where the two ends of each strand meet. Kinks in nicked regions are treated slightly differently, as our criterion is only well defined along an intact strand. A kink at a nick is detected on the opposite strand. As kinks at a nick may diffuse slightly, we define a region of 3 base pairs on either side of the nick on the intact strand. If the intact strand is kinked in this 6 base pair region, then the molecule is considered kinked at that nick.
Since distinguishing fraying from kinking at low Q bp is problematic, we compute kinking only for the most probable values, e.g. Q bp = 8 − 10 for the N s = 10 molecules depicted in Figure 4. The lower Q bp configurations appear in our simulations only because of biased sampling and have a negligible contribution to the equilibrium kinking probability.

A. Dimerization equilibrium
The equilibrium constant for dimerization K dim eq is length-independent to within two times the standard error of the mean for the oxDNA average-base parameterization (Table S2). Varying concentration suggests that the role of excluded volume effects is small between 336 nm and 2.69 µm.
For the sequence-dependent parameterization, there is some variation in K dim eq , but this does not represent length-dependence per-se, rather it is most likely due to sequence variation in the bases at the interface between the duplexes and the sticky ends.  CCG TGC TAC ACC TCC ACC GTT TCA   67   CAG AAT CCG TGC TAG TAC CTC AAT ATA GAC TCC CTT TGA CCT GAC TAT CCT CAC CTC CAC CGT TTC A   CGA TCA TGG AGT TAT ATC TGA GGG AAA CTG GAC TGA TAG GAG TGG AGG TGG CAA AGT GTC TTA GGC A   69   CAG AAT CCG TGC TAG TAC CTC AAT ATA GAC TCC CTT TGA CCC ATG ACT ATC CTC ACC TCC ACC

E8A17
CAG AAT CCG TTT TTA TTT ATC GCC TCC ACG GTG CTT TTT TTT TTT TTT TTT GGC CGT GTT ATC TCG AGT TAG T   73 E8A26   CAG AAT CCG TTT TTA TTT ATC GCC TCC TTT TTT TTT TTT TTT TTT TTT TTT TTC CGT GTT ATC TCG AGT TAG T   73 E8A38   CAG AAT CCG TTT TTA TTT ATC GTT TTT TTT TTT TTT TTT TTT TTT TTT TTT TTT TTT TTT ATC TCG AGT TAG TGC TAG TAC CTC AAT ATA GAC TCC CTA | TGT TAA TAC TTC TCC TAT GAC TTC TAA TTG ACC CAT GAC TAT CCT CAC CTC CAC CGT TTC A   CGA TCA TGG AGT TAT ATC TGA GGG ATA CAA TTA TGA AGA GGA TAC TGA AGA TTA ACT GGG TAC TGA TAG GAG TGG AGG TGG CAA AGT GTC TTA GGC A   97 n Nick 2   CAG AAT CCG TGC TAG TAC CTC AAT ATA GAC TCC CTA TGT TAA TAC TTC TCC TAT GAC TTC TAA TTG ACC CAT GAC TAT CCT CAC CTC CAC CGT TTC A   CGA TCA TGG AGT TAT ATC TGA GGG ATA CAA TTA TGA AGA GGA TAC TGA AGA TTA ACT GGG TAC | TGA TAG GAG TGG AGG TGG CAA AGT GTC TTA GGC A   97 n Double Nick CAG AAT CCG TGC TAG TAC CTC AAT ATA GAC TCC CTA | TGT TAA TAC TTC TCC TAT GAC TTC TAA TTG ACC CAT GAC TAT CCT CAC CTC CAC CGT TTC A   CGA TCA TGG AGT TAT ATC TGA GGG ATA CAA TTA TGA AGA GGA TAC TGA AGA TTA ACT GGG TAC | TGA TAG GAG TGG AGG TGG CAA AGT GTC TTA GGC A   Table S1. Unless otherwise stated, ‡ cyclization sequences are from V&H [S3], with sticky end (Ns) highlighted in red, and the remainder of the sequence (N d ) in black. Sequences are written 5 ′ to 3 ′ . For conciseness, other than the initial example of N bp = 67, a second strand is omitted where it is the reverse complement of the first strand. ‡ Non-V&H sequence. m C:C Mismatch highlighted (blue). n Nicks highlighted with bar | (blue).
For the computation of j eq and j dyn , we use average values for 336 nm. For sequence-dependent results at N bp = 73, a length where the dimerization equilibrium constant was computed, we use the appropriate lengthspecific value instead of the average.

B. Initial stress in cyclized system
There is a general trend towards a larger activation free-energy barrier to cyclization (ΔG ‡ cyc ) at shorter N bp . This is due to the bending stress imposed upon the system by the formation of the initial base pair (Q bp = 1). At the very shortest lengths N bp ≈ 30−45 , this trend is reversed because the complementary single-stranded sticky ends (of length N s = 10) are sufficiently long relative to N d that an initial base pair can form with less bending stress in the duplex region. States at the top of the free-energy barrier (i.e. Q bp = 1) for cyclization and dimerization at N bp = 30, 101 are compared in Figure S1, clearly showing bending in the cyclized molecules, but not the dimerized molecules. Importantly, at N bp = 30, the relatively long single-stranded region (N s = 1 /2×N d = 10) reduces the requirement to bend the duplex ( Figure S1(c)), compared to the strong bending necessary at N bp = 101 ( Figure S1(a)).

Length of complementary sticky ends
We compare oxDNA to available V&H equilibrium data [S3] for N d = 91, N s = 8 − 10 . Specifically, V&H report the fractional occupancy of the cyclized state f cyc . Akin to the nomenclature in Section S1 C, Fractional occupancy can be a somewhat misleading measure by which to compare systems, because the entire transition region f cyc ≈ 0.25 − 0.75 is covered by a free-energy difference of ∼ 2 k B T , while a similar freeenergy difference at the extrema (f cyc → 0 or 1) would be indistinguishable. We therefore make our equilibrium comparison in terms of ΔG cyc , where and OxDNA appears to be in good agreement with experiment; although, this may be partly coincidental as V&H do not report either the temperature (presumed to be 298 K) or a salt concentration (V&H give their imaging buffer cation concentration as [Na + ] = 500 − 1000 mm or [Mg 2+ ] = 10−30 mm; recall that oxDNA is parameterized to [Na + ] = 500 mm). At N s = 9 and N s = 10, oxDNA and experiment differ by ∼ 1 k B T , while at N s = 8, oxDNA under-reports experiment by ∼ 3 k B T (cyclized state is less likely in oxDNA than experiment).
As expected from the SantaLucia model [S7, S8] of nearest-neighbour thermodynamics, the linear relationship between oxDNA results in Figure S2 (a) is reflective of the ∼ 3 k B T stabilization of each additional base pair. A possible explanation for the discrepancy between oxDNA and V&H [S3] is that hairpin formation in the N s region could decrease the fraction of cyclized molecules; but given the sequence, hairpin formation should be negligible. In particular, we cannot explain the discrepancy between oxDNA and the V&H data for N s = 8, although we note that identifying a very low yield in experiment can be challenging, for instance due to the presence of impurities inducing large changes in ΔG cyc .

Nicks and mismatches
We examine the role of structural defects, namely nicks and mismatches, on cyclization using the average-base oxDNA parameterization. Unsurprisingly, we find that nicked and mismatches molecules cyclized more readily than their intact counterparts ( Figure S3). V&H briefly examine the impact of nicks and mismatches on the rate of cyclization, recording f cyc versus time ( Figure S3 in Ref. [S3]). It does not appear that all molecules have reached an equilibrium plateau in f cyc , therefore a concrete equilibrium comparison to oxDNA is difficult. We do, however, observe good agreement between the oxDNA and the experimental observations of V&H regarding the stabilizing effect of various motifs (Table S3). However, we do not reproduce the unexpected similarity between the experimental values for intact sequences at N bp = 69 and N bp = 97. This is consistent with our inability to reproduce V&H's j -factors for their shortest sequences.   Table S3. K cyc eq for oxDNA and experiment (V&H [S3]). The experimental assays were not intended to probe the equilibrium behaviour of the systems; as such, the last point of kinetic cyclization versus time assay may not represent equilibrium values. To emphasize the qualitative difference between intact molecules and those with structural defects (nicks, mismatches), oxDNA results are for the average-base parameterization. * Converted from yield fcyc to K cyc eq using Eq. (S9). † ΔΔG is computed as the difference in the ΔG of cyclization between intact states and those with structural defects (nick or mismatch).  Figure S3. Free-energy profiles for cyclization for (a) an intact duplex compared to a duplex with a single mismatch at N bp = 69, and (b) intact versus nicked duplexes at N bp = 97, using the oxDNA average-base parameterization and Ns = 10 for all sequences (Table S1).

D. Comparison to ligase experiments
Comparison to the ligase experiments of Du et al. [S9] and Cloutier & Widom (C&W) [S10] highlight some interesting features of the oxDNA model ( Figure S4). For example, while both the Du et al. experimental and oxDNA model results are similar to the SY WLC expression for N bp 100, there is a significant absolute deviation in j -factor due to differences in the relevant mechanical properties for oxDNA and those seemingly appropriate for the conditions of the experiment of Du et al . Specifically, they use NEB T4-ligase buffer, which has a salt concentration of [Mg 2+ ] = 10 mm, while oxDNA is parametrized to [Na + ] = 500 mm. Du et al.'s fit to their data gives a weaker torsional stiffness (2.4 × 10 −28 versus 4.75 × 10 −28 J m), longer persistence length (47 versus 41.82 nm) and longer pitch length (10.54 versus 10.36 bp/turn) than for oxDNA. The comparison also highlights how relatively small deviations in DNA mechanical properties (in particular a ∼ 10 % change in the persistence length) may shift the apparent j eq by an order of magnitude. The differences in SY WLC expression parameters are reasonable given the different solution conditions (the persistence length decreases with increasing salt concentration [S11, S12] and temperatures (e.g. j eq is enhanced by a factor of ∼ 10 between 5 ∘ C and 42 ∘ C [S13]). We note that the oxDNA persistence length is reasonable for the high salt conditions used by V&H [S12].
Du et al. only investigated lengths near the maxima in j WLC eq . As a consequence, their experiments are unable to confirm our prediction of a decrease in the magnitude of the oscillations in j oxDNA eq , which we attribute to the affect of teardrop configurations. Whether the experiments could potentially see the effects of teardrop configurations also depends on how the ligase activity depends on the nature of the DNA confirmation. Indeed, Vologodskii et al. mention the possibility of poor ligation efficiency for teardrop configurations [S14]. We also note that because of the potential role of teardrop configurations, the torsional modulus extracted from fits to the SY expression should be interpreted with caution, and it is interesting to note that the values typically obtained from these fits [S5, S9, S15] are significantly lower than those from singlemolecule experiments with magnetic tweezers [S16, S17].
Since it is likely that C&W conducted their experiment at too high a ligase concentration, it is difficult to characterize their j -factor in terms of either j dyn or j eq . To isolate high ligase concentration as the culprit behind an apparent deviation from j WLC eq , Du et al.

E. Capture radius
Comparing results for j dyn with a WLC model requires that the finite dimensions of the sticky ends be taken into account. Defining a "capture radius", a measure of the distance within which the sticky ends can hybridize, is a common way of accounting for the finite size of DNA sticky ends. As oxDNA naturally accounts for the finitesize of the sticky ends, we have the ability to test the capture radius assumption.
In oxDNA, we define the capture radius R capture as the distance between the first and last bases of the duplex region (N d ) at Q bp = 1 ( Figure S1 (a,c)). Nucleotide positions are used to report distance, instead of the midpoint along the helical axis, because fraying may occur at the ends of the duplex. The choice of nucleotide strand does not impact our results.
Our ensemble average capture radius, ⟨R capture ⟩ is shown in Figure S5. ⟨R capture ⟩ is roughly constant at N bp ≈ 100 − 200 , increases slowly between N bp ≈ 50 − 100 , and then more rapidly for N bp ≈ 40 − 50 , before decreasing for N bp 37. The saturation value (⟨R capture ⟩ ≈ 4 nm) is quite close to the capture radius of 5 nm that has been used to interpret V&H's experiments with 10-base sticky ends [S3, S18].
The free-energy cost associated with forming the initial contact is partitioned between that for bending the duplex and that for stretching the single-stranded stickyends. At shorter N bp , the amount of bending required for initial binding increases, as does the stretching of the single-stranded sticky-ends. At short lengths, this trend is reversed when the length of the unperturbed duplex approaches the capture radius at Q bp = 1. Interestingly, the distribution of R capture is roughly Gaussian, with similar standard deviation for all lengths.
While we explicitly disallow misbonding, allowing misbonding may slightly increase the capture radius, an effect that should be more prominent for longer sticky-ends. Overall, the 5 nm value for the capture radius used in the interpretation of V&H's experiments appears not unreasonable.