Global and local mechanical properties control endonuclease reactivity of a DNA origami nanostructure.

We used coarse-grained molecular dynamics simulations to characterize the global and local mechanical properties of a DNA origami triangle nanostructure. The structure presents two metastable conformations separated by a free energy barrier that is lowered upon omission of four specific DNA staples (defect). In contrast, only one stable conformation is present upon removing eight staples. The metastability is explained in terms of the intrinsic conformations of the three trapezoidal substructures. We computationally modeled the local accessibility to endonucleases, to predict the reactivity of twenty sites, and found good agreement with the experimental data. We showed that global fluctuations affect local reactivity: the removal of the DNA staples increased the computed accessibility to a restriction enzyme, at sites as distant as 40 nm, due to an increase in global fluctuation. These results raise the intriguing possibility of the rational engineering of allosterically modulated DNA origami.

Very recently, our group discovered unexpected allosteric behavior in a DNA triangle (19), a DNA origami first introduced by Rothemund, and subsequently used in several studies (20,(28)(29)(30). Specifically, we investigated the action of several sequence-specific, DNA-cutting enzymes (restriction endonucleases or REases) towards their recognition sites present in the M13 scaffold sequence. We showed that REases act in a binary fashion, as certain sites cannot be cut even at extended reaction times, or, for HhaI REase recognition sites in particular, become reactive following the introduction of a distant structural defect in the triangle (i.e. by omitting, in the self-assembly process, four staples at a distance of ∼40 nm from the REase site).
In this study, we computationally examine this phenomenon to obtain insight on the global mechanical behavior of the DNA origami triangle and its influence on local parameters such as site reactivity towards restriction endonucleases. We studied three versions of the triangle: the first contains the full complement of staples; the second and the third are triangles deficient in four and eight staples, respectively, in a localized region on the bottom trapezoid. We determined the global mechanical properties of the three triangles, as well as the detailed dynamics of the three trapezoids substructures composing the triangle, and compared their fluctuations with that of the isolated trape-zoid. We then correlated these properties with the reactivity of individual sites towards cleavage by HinP1I REase (a neoschizomer of HhaI).
Our findings are as follows. First, the absence of selected staples affects the amplitude of the global conformational fluctuations of the triangle. Two metastable states are observed in the triangle, while omission of four staples lowers the free energy barrier between the two states, and omission of eight staples completely abolishes the barrier, with the free energy profile reducing to a single minimum. The metastability originates from an intrinsic property of the isolated trapezoid, which exhibits two alternate conformations. When coupled into the triangle, the three trapezoids fluctuate synchronously between the two states, thereby giving rise to the metastability of the triangle. Second, we devised a computational method that estimates the steric accessibility of REase sites in DNA origami, and obtained good agreement between the experimental and theoretical results. Specifically, the results reveal that site accessibility depends on local conformational fluctuations. Third, we correlated theoretical site accessibility with the global conformational fluctuations of the triangles and of the isolated trapezoid, thereby showing how endonuclease reactivity in a strongly coupled mechanical system can be modulated in an allosteric fashion.

System setup and simulation details
We performed extensive Langevin molecular dynamics simulations to determine the conformational dynamics of three triangles that differ in staple composition, and of the isolated trapezoid substructure too. The fully intact structure, hereafter referred to as 'triangle', has the scaffold organized by the full complement of staples ( Figure 1A). The second and the third structures lack four staples and eight staples, respectively, in a region that bridges the seam of one of the three trapezoidal elements forming the triangle, and thus they are referred to as 'four staple deficient (4sd) triangle' and 'eight staple deficient (8sd) triangle' (Figures 1B and C). The fourth structure is an isolated trapezoid identical to those comprising the triangle ( Figure 3A). The first two structures were previously studied experimentally in (19).
The triangles and the trapezoid are composed of about approximately 14 000 and 4700 nucleotides, respectively, and were modelled using the second version of oxDNA (31), a coarse-grained DNA model that has been shown to accurately reproduce the mechanical properties of doublestranded and single-stranded DNA molecules (32)(33)(34)(35)(36)(37)(38) and that parametrize major and minor grooves and different monovalent salt concentrations.
To obtain the initial configurations, we started from the cadnano (39) constructs of ref (19) (see Supplementary Figures S1-S4), and converted them to the oxDNA representation using the tacoxDNA package (40). In the case of the triangles, the initial structure presented the three trapezoids in sequential order, thus they were roto-translated to obtain the final geometry, and over-stretched bonds were relaxed using the protocol described in ref. (40).
Simulations were carried out with T = 300 K and at 1 M monovalent salt concentration. The time of evolution of the four structures was monitored for about 20 ms for the triangles, and 10 ms for the trapezoid; the mapping between the intrinsic and standard time was previously found to yield LJ ∼ 0.7 ns (41,42) and the configurations were sampled every 7 s (for a total of ∼3000 configurations). Each simulation was run on 140 processors for a maximum of 1.5 months using LAMMPS (43,44). The initial setup files are available as Supplementary files.

Docking of HinP1I to GCGC sites and accessibility measurement
We developed a quantitative method to measure the accessibility of HinP1I REase towards its recognition sequence (GCGC) in the DNA origami. The first step consisted of superimposing an ideal oxDNA B-DNA helix of 10 bp onto the atomistic DNA structure bound to HinP1I in the PDB structure (PBD code 2FL3, (45)), using a rigid body motion (46). The ideal double helix was then superimposed onto each GCGC sequence of the triangle, for each frame of the trajectory. The same transformation was applied to the protein.
To determine whether the protein docked to a restriction site overlaps with an adjacent DNA strand, two distance maps were built and compared. These maps, max i, and min i, are functions of the site base-pairs (i = 1, . . . 10) and of the angle around each of them ∈ [ − , ], see Supplementary Figure S5 for more details. The maximum distance map, max i, , reports the distance between the site and the protein particle that is farthest from it. The minimum distance map, min i, , reports the distance between the site and the closest particle from an adjacent DNA strand. Clearly, if max i, < min i, for every i and , the enzyme satisfies the accessibility requirement and the configuration is classified as 'accessible'.

Global mechanical properties and metastable conformations of the DNA triangle
For a first characterization we computed the average structure and the root mean square fluctuation (RMSF) of each nucleotide around the mean position for the three triangles, see Figures 1D-F. To this end, we applied optimal rotations and translations (so as to minimize the root mean square deviation with respect to the initial configuration) to each frame of the trajectory using the Kabsch algorithm (46). For all triangles, the average structure is globally twisted, due to the unnatural helical pitch of 10.67 bp/turn imposed by the crossover structures (1,47). The fluctuations, highlighted in red, are mostly concentrated at the vertices that exhibit an RMSF >6.5 nm and also near the defects in the 4sd and 8sd triangles.
At first glance, the triangle and the 4sd triangle do not exhibit evident differences that might rationalize the experimental results of (19), except for the area bearing the defect and the edges of the trapezoids, which show slight variations. A deeper analysis of the trajectories, however, show that the triangles fluctuate in different ways over time (see Supplementary Movies S1-S3). To distinguish between the different global fluctuation modes, we performed principal component analysis (48) over the trajectories using GROMACS (49). We constructed the covariance matrix for the three different triangles as where t is the frame index and x n is the nth component of the 3N-dimensional vector of deviations from the average position (n = 1, . . . 3N are the coordinates of the scaffold only, in order to compare fluctuations across triangles using the same number of particles), and we computed eigenvalues and eigenvectors ({ n , p n }, λ 4sd n , p 4sd n and λ 8sd n , p 8sd n ) for the triangle, the 4sd and 8sd triangles, respectively.
We projected the first 10 normalized eigenvectors of the triangle on the first 10 ones of the 4sd and 8sd triangles, p i |p 4sd j and p i |p 8sd j ( Figures 1G and H), and observed a good collinearity for i = j between triangle and 4sd triangle eigenvectors (except for the inversion between eigenvectors 4 and 5), while for the 8sd triangle only the first eigenvector shows good collinearity. Thus, we can quantitatively compare the mean square fluctuations of the triangle and 4sd triangle along the triangle eigenvectors ( Figure 1I), and observe that removal of four staples perturbs mainly the first mode, increasing the fluctuation amplitude of 101 nm 2 . In contrast, removal of eight staples heavily affects the directions of the principal modes, thus a precise comparison could not be made.
To further analyze the first fluctuation mode, we projected the trajectories over p 1 (see Figure 2A). This allowed us to analyze the fluctuations of the first mode as a function of time. What emerges is the existence of two metastable conformations, with transitions observed in the triangle and 4sd triangle. In particular, we observe that the triangle resides in each metastable state much longer than the 4sd triangle, while for the 8sd triangle there are no metastable states. By reconstructing the free energy profiles as a function of x|p 1 ( Figure 2B), we estimate for the triangle an energy barrier between the two minima of about 5k B T and 2k B T going from the left minima to the right one and vice versa, respectively. For the 4sd triangle, these barriers decrease to about 2.5k B T and 0.15k B T. The smaller energy barrier enhances the rate of exchange between the two basins. The 8sd triangle presents only a single minimum.  The two dominant metastable conformations of the triangle that are associated with the two free energy minima are shown in Figure 2B: the triangle in State 1 is concave upward, while in State 2 it is concave downward. Here, the term upward refers to the direction that is normal to the plane containing the triangle, and that points inward relative to Figure 4A, while the term downward refers to the opposite direction (see also Supplementary Movies S1 and S2, in which States 1 and 2 are colored in blue and red).

Mechanical properties of the trapezoid in the DNA triangle
To understand the mechanical origin of the two metastable states, we analyzed separately the time evolution of each of the three trapezoids composing the triangles, designated as left (L), right (R), bottom (B), and compared them with the time evolution of the isolated trapezoid ( Figures 3A and B).
We performed a principal component analysis for each trapezoid (considering only the scaffold), producing a set of eigenvalues and eigenvectors. In particular, the first and most relevant eigenvector of the isolated trapezoid ( p t 1 ) is naturally conserved across all the triangle trapezoids (Supplementary Figure S6), and thus can be used to project all the trajectories to analyze their fluctuation.
The first eigenvector of the isolated trapezoid describes the transition between two alternate conformations (Figure 3C and D), differing in their concavity with respect to the surface normal (with state 1 less probable than state 2). When three trapezoids are embedded in the triangle, the geometrical constraints restrict the fluctuations of the three trapezoids between these two states, as evident from the probability distributions of x|p t 1 ( Figure 3E). Moreover, the fluctuations are highly correlated, and synchronized diffusion between state 1 and state 2 correspond to jumps across the two global metastable states of the triangle. In the 4sd triangle, a similar situation occurs ( Figure 3F), but the distribution is broadened by the presence of a defect in the B-trapezoid, a feature that ultimately results in a lowering of the free energy barrier seen in Figure 2B. In the 8sd triangle the defect is so large that the trapezoids exhibit uncoupled oscillations. Notably the L-trapezoid is preferentially in state 1, while the R-trapezoid is primarily in state 2. This reflects the fact that the L-and R-trapezoids are coupled from two sides, with only the the defective trapezoid side being uncoupled.

Theoretical accessibility of REase recognition sequences
We developed a method to determine the local accessibility of REase recognition sequences in the triangle. We created a structural proxy by docking the crystallographic structure of the REase to the DNA site under study, then evaluated, for each time frame, the degree of volume overlap between the enzyme and the adjacent DNA structures. A site is scored as accessible if no clash occurs between the two volumes, and inaccessible otherwise (see Figure 4B and methods for the exact procedure). Finally, we defined the theoretical site accessibility as the fraction of time frames in which the site is accessible to the enzyme.
Since the HhaI REase used in (19) lacks a PDB structure, we used HinP1I as a proxy enzyme, as it recognizes the same sequence (GCGC), has a similar mass as HhaI (HinP1I, 28.7 kDa; HhaI, 27.8 kDa), and a crystal structure of HinP1I docked to a duplex DNA has been reported (PDB: 2FL3 (45)). Note that HinP1I functions as a monomer, and cuts each strand of the GCGC duplex in separate events (45).
The triangle has twenty GCGC sites (see Figure 4A), with each exhibiting a scaffold side and a staple side. We distinguished three situations, in which HinP1I can access (i) only the scaffold strand, (ii) the staple strand or (iii) both strands simultaneously (see Figure 4C). We excluded sites 10, 11 and 15 (shaded) as they overlap with a DNA crossover junction and therefore are inherently unreactive. The resultant picture is that four sites exhibit high accessibility: 12, 7, 19, 16 (in ascending order).
We can quantitatively correlate theoretical accessibility to experimental reactivity by proposing that the former quantifies the fraction of reactive triangle conformations. The enzymatic reaction for a specific site i can then be described as: where E is the enzyme, T u are the uncut triangles at site i, r is the fraction of accessible T u triangles (0 < r < 1), T c are the cleaved triangles. For a given concentration of T u , each site i will be accessible to the enzyme with a different concentration fraction rT u . This will be reflected in the kinetics of cleavage: the more accessible sites (with a larger r) will be cut more rapidly (50). We define the experimental accessibility as the fraction of triangles cleaved at a single time point, before the reaction at the most reactive site reaches completion. Figure 4D shows the comparison between the theoretical accessibility for HinP1I and the experimental accessibility, obtained from gel electrophoretic data for a one hour reaction by HhaI (19) (see also Supplementary Figure S7). For both sets of data we considered only the accessibility of the scaffold strand, as the experimental data provide information only for that side.
We observed that the theoretical accessibility matches the experimental accessibility within an error margin of 20%, a value comparable to the statistical background 'noise' of the gel electrophoretic data. In particular, site 7 is cut more slowly than sites 16 and 19, as predicted by the model. These results support the assumption that the reactivity of a given GCGC site is related to its accessibility. Only sites 4 and 12 are an exception, for multiple possible reasons. Firstly, the use of monovalent instead of divalent cations (implicit in oxDNA) might locally alter the modeled structure, thereby resulting in subtle changes in helix-helix distances (51) and relative orientation. Secondly, the experiments in Supplementary Figure S7 Figure S7). The Pearson's correlation coefficient is 0.9.
structural differences between HhaI and Hinp1I enzymes may be relevant for specific docking angles. Lastly, there might be a threshold of experimental accessibility, above which each site can be scored as experimentally reactive.

Local mechanical fluctuations affect endonuclease accessibility
To obtain insight on the local factors that influence site accessibility, we projected the center of mass of the adjacent DNA double-helical structures over a plane perpendicular to the GCGC helical axis, using as a reference frame the direction of one of the GCGC nucleotides (see Figure 5A). Figure 5B shows the projection for three representative sites, that also distinguishes between scaffold and staple sides: site 6, non-accessible from both sides; site 7, accessible from both sides; and site 8, accessible only from the staple side. The local accessibility is determined by the degree of overlap between the HinP1I enzyme projection and the adjacent DNA molecules. Complete overlap prevents cleavage, while partial overlap affords accessibility, with a cleavage rate proportional to the percentage of overlap over time. In particular, for site 6, we find that the enzyme always overlaps with the upper adjacent DNA molecule, thus suppressing the site accessibility. For site 7, the enzyme overlaps only 50% of the time with the upper adjacent DNA molecule, due to the higher flexibility of that molecule. For site 8, we find the enzyme overlap with the left adjacent DNA molecule is the most important to suppress site accessibility from the scaffold side, but the same does not apply from the staple side.
We note that the degree of overlap between the enzyme and the adjacent DNA molecules depends on several variables, such as the orientation of the bases in the GCGC sites (from which the orientation of the docking protein depends), and how much the adjacent DNA molecules fluctuate, which in turn is dictated by the distance from a crossover junction.

Global mechanical fluctuations affect endonuclease accessibility in an allosteric manner
We next compared the local site reactivity between the three triangles and the isolated trapezoid, in order to explore its possible connection with changes in global mechanical fluctuations ( Figure 6). To this aim, we mapped the restriction sites of the triangle onto the equivalent positions on the trapezoid. Qualitatively, we observed that the most remarkable accessibility variation across the studied structures is seen for sites 4 and 8 from the staple side (13% and 18% increase, respectively), and for site 7 from the scaffold side (20% increase).
Sites 7 and 8 happen to be in a region of the triangle about 40nm far from the defective region, and far from the border connected to other trapezoids in the case of the single trapezoid. We therefore examined whether there is a correlate of the increase in accessibility, which is mainly associated with local fluctuations (fast modes), to global fluctuations   Figure 4C for more details). The most consistent increase in accessibility across all sites is for sites 4, 8 from the staple side and site 7 from the scaffold side.
(slow modes). We observed that fluctuations across different eigenvectors are highly correlated to each other (e.g. Supplementary Figure S8). For this reason, we speculate that an increase in fluctuation over the slower modes can propagate over the faster ones. This is corroborated by the fact that the total mean square fluctuation (i.e. the sum of all the eigenvalues) of the trapezoids containing sites 7 and 8 shows the same trend as its accessibility (Figure 7). These results suggest that allosteric effects may be operative in the DNA triangle: while the GCGC sites 7 and 8 are distant from the B-trapezoid, a structural change (staple omission or addition) in the latter can modulate the reactivities of site 7 and 8 by altering its dynamics. In the experiments of ref. (19) it was observed that sites 9, 12 and 13 become reactive in the 4sd triangle. Although only site 12 exhibits an increase in local accessibility in the simulations, multiple other factors mentioned above, and not accounted in the model, can contribute to a change in accessibility along with the observed change in global fluctuations.

DISCUSSION
DNA nanostructures, including DNA origami, constitute a novel class of materials with emergent properties and behaviors yet to be fully understood and applied. Here, we used extensive coarse-grained molecular dynamics simulations to probe the global and local mechanical properties of a triangular DNA origami. We found that two metastable conformations are present, with the triangle capable of assuming either an upward or downward concave shape with respect to the triangle surface normal. Upon omission of four staples, the interconversion rate between the two states increases, due to a decrease in the intervening free energy barrier. By removing eight staples over the same defective region, the metastability is lost. These mechanical fluctuations can be understood in terms of mechanical fluctuations of the trapezoid substructures, which intrinsically possess two distinct alternate states. In the triangle, the three trapezoids exchange between the two states in a synchronous manner, resulting in global transitions between the two metastable conformations of the triangle. With the four staple defect, the fluctuations are synchronous, but more frequent. With the eight staple defect, the synchronicity is lost.
We devised a method to quantify the local accessibility of GCGC sequences to a cognate restriction endonuclease, and showed that site reactivity on the scaffold side correlates with the predicted accessibility. In particular, the reactivity of a site is dependent upon the degree of steric overlap between the docked enzyme and the adjacent double helices, which in turn depends on the local structure and conformational fluctuations of the site.
For the three triangles and the trapezoid structures that were analyzed we found that several sites exhibit a positive trend in increased accessibility. These changes in local fluctuations can be in part correlated to changes in global fluctuations, which in general are due to modifications distant from the sites themselves.
This study provides an initial step towards a deeper understanding of global and local mechanical properties of DNA origami, and how these properties can play a role in the control of the DNA accessibility (and reactivity too) to a protein. Future studies intend to explore the relationship between global fluctuation modes and the design arrangement of staples e.g. using a honeycomb lattice, and the quantitative assessment of the kinetic behavior of enzymes acting on these nanostructures (see, e.g. (20)). These advancements in turn should enable the rational design of allosteric metamaterials, whose shape and conformational fluctuations are controllable by appropriate localized stimuli (52)(53)(54)(55). Such behaviors can be exploited to develop, for example, carrier platforms with programmable reactivities toward nucleases or with dynamic stimuli-responsive properties.