Investigating the sequence-dependent mechanical properties of DNA nicks for applications in twisted DNA nanostructure design

Abstract DNA nick can be used as a design motif in programming the shape and reconfigurable deformation of synthetic DNA nanostructures, but its mechanical properties have rarely been systematically characterized at the level of base sequences. Here, we investigated sequence-dependent mechanical properties of DNA nicks through molecular dynamics simulation for a comprehensive set of distinct DNA oligomers constructed using all possible base-pair steps with and without a nick. We found that torsional rigidity was reduced by 28–82% at the nick depending on its sequence and location although bending and stretching rigidities remained similar to those of regular base-pair steps. No significant effect of a nick on mechanically coupled deformation such as the twist-stretch coupling was observed. These results suggest that the primary structural role of nick is the relaxation of torsional constraint by backbones known to be responsible for relatively high torsional rigidity of DNA. Moreover, we experimentally demonstrated the usefulness of quantified nick properties in self-assembling DNA nanostructure design by constructing twisted DNA origami structures to show that sequence design of nicks successfully controls the twist angle of structures. Our study illustrates the importance as well as the opportunities of considering sequence-dependent properties in structural DNA nanotechnology.

where represents the coordinate vector of BP step parameters as = [τ ρ ω Dx Dy Dz] , whose average is denoted as the angle bracket, is displacement vector whose components indicate the displacement from average configuration, and is a positive definite stiffness matrix where the six diagonal and remaining symmetric off-diagonal terms represent mechanical stiffness (bending, torsion, shearing, and stretching) and coupling stiffness, respectively.
The probability distribution of BP step parameters can be approximated to an N-dimensional Gaussian function assuming moderate fluctuations under harmonic potential in contact with a heat bath of a certain temperature. The covariance matrix of the BP step parameters was used to induce the stiffness matrix of a BP step (1-3) as = k B T −1 (2) where k B is the Boltzmann constant, T is the absolute temperature used in MD simulation, and is the covariance matrix. Here, the covariance matrix was obtained using from the sampled MD trajectories of the BP step parameters (4) as The strain energy can be rewritten as the sum of diagonal and off-diagonal components as (4) where i and j represent the BP step parameters as τ, ρ, ω, Dx, Dy , and Dz . Diagonal terms contribute to mechanical rigidities corresponding to the six BP step parameters, and mechanical coupling coefficients are derived from off-diagonal terms between two BP step parameters. Mechanical rigidities of a BP step were then obtained by multiplying the diagonal stiffness components to the axial length of a BP step from linear elastic theory as Likewise, mechanical coupling coefficients of a BP step were calculated by multiplying the offdiagonal stiffness to the axial length of the BP step as g ij = K ij 〈Dz〉 (8) where g ij represents the mechanical coupling coefficients of two different BP step parameters. For example, g ωDz represents the mechanical coupling coefficient of twist and rise, implying the twiststretch coupling.

Supplementary Method S2. Derivation of mechanical rigidities of a base-pair (BP) step.
Mechanical rigidities were derived as assuming DNA structures of a homogeneous linear elastic material. The initial shape was assumed to have an area (A) and length (L) with axial z-direction. Then, the representative mechanical rigidities (B τ , B ρ , C, Y Dx , Y Dy , and S) were derived by elastic theory.
Bending rigidity (EI) is then evaluated as where k bending is the stiffness defined by the ratio of applied bending moment and corresponding bending deformation. Bending rigidities of DNA, therefore, can be calculated in terms of components of stiffness matrix and BP step parameters as where B τ and B ρ represent tilt and roll bending rigidities, and K ττ and K ρρ indicate the diagonal tilt and roll components in the stiffness matrix, respectively.
(B) Torsional rigidity. Torsional moment (M z ) is obtained by integrating the shear stress (σ τ ) multiplied by the moment arm (r) over the cross-section area (A) as where shear stress (σ τ ), shear strain (γ), and the second polar moment of area (J) are defined as Then, the torsional rigidity (GJ) is evaluated as where k torsion is the stiffness defined by the ratio of applied torsional moment and corresponding torsional deformation. Therefore, we obtained the torsional rigidity of DNA corresponding to the twist parameter as where C is generally referred to torsional rigidity of DNA, and K ωω is the diagonal twist component in the stiffness matrix.
(C) Shearing rigidities. Shearing rigidity is derived from the shear modulus (G) defined by the ratio of shear stress (σ τ ) to the shear strain (γ) as where F and ∆L are the applied shear force and corresponding shear displacement. The shearing rigidity (GA) is then obtained as where k shearing is the stiffness defined by the ratio of applied shear force and corresponding displacement. We then obtained the shearing rigidities of DNA as where Y Dx and Y Dy denote shearing rigidities corresponding to the shift and slide, and K DxDx and K DyDy are the diagonal shift and slide components in the stiffness matrix, respectively.
where k stretching is the stiffness defined by the ratio of applied axial force and corresponding displacement. Therefore, the stretching rigidity of DNA is obtained corresponding to the rise mode as where S is generally referred to the axial rigidity of DNA, and K DzDz is the diagonal component corresponding to the rise parameter in the stiffness matrix.

Supplementary Method S3. Derivation of an equivalent isotropic rigidity in bending and shearing.
Assuming a homogeneous linear elastic material, if a transverse load was applied to a beam, whose axis system coincides with the principal axes of bending, three-dimensional beam theory implies the following displacement field as where z is the axial direction of the beam, x and y are the planar directions on the cross-section of the beam, and θ x and θ y are rigid body rotations of the cross-section respectively. The strain field can be evaluated by the differentiation of the displacement field as ε xx = ε yy = γ xy = γ yz = γ xz = 0 where the curvatures are defined as The axial stress is described by Hooke's law as The strain energy (π) stored in the beam with anisotropic bending rigidities is evaluated as where I x and I y are the second moment of area with respect to x-and y-directions on the crosssection.
When an equivalent isotropic beam is assumed, the force equation under a bending moment is found as If the strain energy stored in anisotropic and equivalent beams is the same, the following equation is satisfied π = ∫ We performed MD simulations of 6HB DNA origami structures with BP-insertion (Supplementary Figure   S9). At each end of the simulated structure, 7-BP-long strands were added to eliminate the boundary effect. Since it was difficult to observe the effects of sequence control of nicked BP steps due to the limited system size, MD result of an arbitrary sequence was used to predict an approximately twist angle of 6HB structures. For the fast convergence of equilibrium states, the initial all-atom structures were developed after 100-ps-long the elastic network-guided simulation (5). All-atom explicit solvent simulations were then performed after energy minimization as previously described for DNA oligomers.
To smoothly convert the induced strain energy by BP-insertion into global torsion of structures, we applied weak harmonic constraints to Watson-Crick base-pairing from 0.5 kcal mol -1 Å -2 to zero during 60-ns-long simulation as previous MD studies on DNA origami structures (5,6). Final 20-ns-long equilibrated trajectories were sampled to calculate twist angle (Supplementary Method S5).

structure.
The twist angle of a 6HB structure was calculated by the summation of local twist angles of twisting blocks between sequential hexagonal cross-sections locating at Holliday-junctions. To obtain the local twist angle between hexagonal cross-sections, we defined six vertexes locating at helices of the 6HB structure as the averaged point of local origins in four neighboring BPs next to the hexagonal crosssection. A hexagonal cross-section was determined to satisfy the smallest sum of the distance between each vertex and an arbitrary plane. Auxiliary vectors (A ⃗ ⃗ 1 , A ⃗ ⃗ 2 , A ⃗ ⃗ 3 ) were defined in hexagonal crosssections: A ⃗ ⃗ 1 , A ⃗ ⃗ 2 , and A ⃗ ⃗ 3 were defined as the vectors from vertex 1 to 4, 2 to 5, and 3 to 6, respectively.
Subsequently, the triads in hexagonal cross-sections were determined: x-vector ( x ⃗ ) representing normalized axial direction was defined as a normal vector of each cross-section, y-vector (y ⃗ ) in the hexagonal cross-sections was described as a normalized vector with projection of an averaged vector of auxiliary vectors on the cross-section, and the z-vector (z ) was calculated using the cross product of x-and y-vector. A rotation matrix ( ) between sequential triads was then obtained using the triads as where i and i + 1 represent successive hexagonal cross-sections. The rotation vector ( ⃗ ⃗ ) equivalent to the rotation matrix was obtained as where W ij is a component of the skew-symmetric matrix ( ) defined as = 1 2 ( − T ). The local twist angle between sequential hexagonal cross-sections (Ω) was finally calculated as We obtained the Gaussian mean and standard deviation of the twist angle in the unit twisting block from MD simulation (Supplementary Figure S9). If the twist angle of the unit twisting block (Φ u ) follows Gaussian distribution (N u ) as where Φ u0 and σ u are the Gaussian mean and standard deviation of the twist angle of the twisting block, we could calculate the twist angle (Φ) of a DNA structure, in which same twisting blocks are connected in sequence and the number of twisting blocks is m, as the sum of the independent Gaussian distributions of the twisting blocks based on the probability theory as Since ten twisting blocks were employed in the entire structure, the approximate range of twist angle of the structure was estimated as the mean of 10Φ u0 and the standard deviation of √10σ u using those of unit twisting blocks (Φ u0 and σ u ) by MD simulation.

Supplementary Method S7. Relation of the trans ratio with the twist angle of 6HB DNA structures.
We followed the details of derivation in the previous study (10). Assuming that 6HB origami structures fluctuate in a harmonic potential, the Gaussian distribution of twist angle can be derived from the Boltzmann distribution as where N 0 is the normalization coefficient as 1/(σ√2π) , E is torsional harmonic energy as E = k t (Φ − Φ 0 ) 2 /2, Φ is the twist angle of 6HB structure, Φ 0 is the equilibrium twist angle, σ is the standard deviation of twist angle as k t /(k B T), and k B T is the product of the Boltzmann constant, k B , and the absolute temperature, T. Here, the torsional stiffness of 6HB structure (k t ) has a relation with torsional persistence length (L P ) and the length of 6HB structure (L) as where the torsional persistence length of 6HB DNA structures (L P ) was reported as 530 nm in a previous study (11), and the length of our structure is approximately 300 nm in AFM images (Supplementary The probability to deposit into trans monomer of the 6HB structures is obtained as where k is an integer as a parametric variable (⋯ , −2, −1, 0, 1, 2, ⋯ ). This probability to deposit into trans monomer is a kind of triangular function, which is plotted as The trans ratio, TR(Φ 0 ) in equilibrium is derived by the average probability of twist angle weighted in Φ 0 for all Φ domain as Using the Gaussian integrals, We then obtained the four integral terms in the equation (51) as Hence, by substituting the equation (60a-d) into the equation (51), the relation of trans ratio with twist angle is derived as Here, the arguments in the k-summation rapidly approach to zero when Φ 0 is out of the range of

Supplementary Notes
Supplementary Note S1. The sensitivity of mechanical properties in simulating environments and neighboring sequence variation.
In this study, we investigated the mechanical properties of nicked BP steps, suggesting that BP steps To verify, we simulated all the BP steps for 100 ns employing the parmbsc1 force field (12), another widely used force field in DNA oligomer simulations, to obtain the mechanical properties and the C B ⁄ ratios. As in the previous results obtained using the CHARMM36 force field, the parmbsc1 force field provided the sequence-dependent properties for regular and nicked BP steps (Supplementary Figure S8 and Supplementary Tables S8-S13). It was also confirmed that nicks led to a significant reduction in torsional rigidity regardless of the force field used. On the whole, the geometry, mechanical rigidities, and coupling coefficients predicted by two force fields were comparable for both regular and nicked BP steps. While the parmbsc1 force field predicted relatively smaller torsional and shearing rigidities than the CHARMM36 as found in previous studies (12,13), the C B ⁄ ratios predicted by the parmbsc1 were not significantly deviated from those obtained using the CHARMM36 because the torsional rigidity decreases for both regular and nicked BP steps Figure S8). Therefore, we could confirm that the trend of C B ⁄ ratio (and hence our design approach) is not sensitive to the choice of force field.
While limited, we increased the time duration from 100 to 300 ns for the GC/GnC step as a selected subset. We also performed additional simulations of 300 ns using parmbsc1 force field with 100 mM of NaCl (Supplementary Figure S5). It was apparent that RMSD values were maintained within 2-8 Å , suggesting that oligomers stayed near an equilibrium state. Results indicated that the mechanical properties of BP steps did not considerably alter compared to original simulation either the simulation time was extended or salt concentrations was changed from MgCl2 to NaCl concentration.
Furthermore, to investigate the effect of the next-nearest-neighbor sequence variation, another additional MD simulations performed for nicked GG/CnC and AA/TnT steps as a limited subgroup in order to evaluate the effect of adjacent sequence variation of a BP step (Supplementary Figures S6 and S7). The dinucleotide GG/CnC step analyzed in this study is the same with the tetranucleotide AGGA/TCnCT step when considering neighboring sequences for original DNA oligomer (Supplementary Table S1). Likewise, the AA/TnT step was identical to with the hexanucleotide GAAAAG/CTTnTTC step. We modified the sequence of the original oligomer for GG/CnC step from GACT4-GAC-TCnCT-CAG-TCAG4 to GACT4-GAC-ACnCA-CAG-TCAG4 to have the same dinucleotide GG/CnC step in the center but different tetranucleotide TGGT/ACnCA step. Likewise, the oligomer with AA/TnT step was also modified from GACT4-GA-CTTnTTC-AG-TCAG4 to GACT4-GA-GTTnTTG-AG-

Supplementary Note S2. Mechanical analysis of DNA bundle structures with BP-insertion.
We theoretically investigated mechanical parameters related to the twist angle of a DNA bundle structure with BP-insertion. A two-helix-bundles (2HB) structure was introduced as the smallest twisting structure by BP-insertion (A), where terminal ends were constraint by Holliday-junctions. It was assumed that two identical and straight helices (H1 and H2) have the radius of r, bending rigidity of B, and torsional rigidity of C. Also, there was nothing to induce strain energy since the end of two helices cross exactly at both ends.
When BPs are inserted but one of the Holliday-junctions is not constrained (B), the initial twist angle due to the helicity of DNA can be regarded as where 〈ω〉 is the mean twist angle of each helix and the number of inserted BPs was denoted as N BP .
Also, each helix is elongated to a longitudinal length of L.
In equilibrium with the constraints at both ends of the helices (C), bending and torsional strain energy is induced at each helix to align the mismatch of connectivity and the center of twist. In addition, bending and torsional angle should satisfy a geometric constraint equation (g) as where θ and Φ represent the equilibrium bending and twist angle of a helix. Then, the induced strain energy in the helix can be described using the method of Lagrange multipliers as where the bending (M b ) and torsional moment (M z ) were derived in Supplementary Method S2 as where the zero initial bending and initial twist angle by BP-insertion yield the angle deformation as Then, the strain energy is derived using the above equations as The equilibrium configuration minimizes this strain energy, suggesting that the gradient of the strain energy is zero vector. We, therefore, obtain the system of equations as Here, the bending angle (θ) can be substituted to twist angle (Φ) with the constraint equation (72), providing a relation of twist angle with rigidities of each helix and the number of inserted BPs as This relation suggests that DNA bundle structures could be more twisted when BPs are more inserted or C B ⁄ is increased in helices. We found that nicks have significant effects on the reduction of torsional rigidity (C) rather than the reduction of bending rigidity (B). It is, therefore, speculated that the twist angle of bundle structures could be controlled by adjusting the location of nicks to control C B ⁄ value. Here, we defined that C B ⁄ ratio of nicked BP steps normalized by corresponding regular BP steps as ⁄ ) regular to quantify the effects of nicks from an ideal structure consisting of all regular BP steps. Utilizing the C B ⁄ ratio data from MD simulation, the twist angle of the structures with high C B ⁄ ratio could be deliberately designed greater than that with low C B ⁄ ratio, suggesting that the torsional energy stored by BP-insertion could be controlled by adjusting nicked BP steps. Future studies based on nonlinear analysis is promising to produce a more realistic theoretical prediction on the DNA structures since our approach is limited to only provide the insight of finding parameters to affect bundle structure due to the simple and linear assumption.

Supplementary Note S3. Comparison of the results of the reverse-phase cartridge (Bio-RP) and the polyacrylamide gel electrophoresis (PAGE) methods in the purification of staples.
We investigated whether the purification method used in the synthesis of oligonucleotides affects the trans ratio of assembled DNA origami structures. Here, 60 staples comprising ten twisting blocks of the flexible 1 design with 2-BP-insertion were purified using the polyacrylamide gel electrophoresis (PAGE) method, in which the manufacturer guarantees a higher purity of more than 95%. We obtained the matrix-assisted laser desorption/ionization time-of-flight mass spectrometry (MALDI-TOF) data for representative ten staples (one staple per each twisting block) purified by both methods (Supplementary Table S19). According to the result of MALDI-TOF, purified staples using Bio-RP and PAGE showed almost identical mass spectrum with a single and clear peak, indicating only oligonucleotides with the desired length existed in the stock solution (Supplementary Figure S11). Noting that the maximum error of the MALDI-TOF spectrum is typically around 0.24%, we could confirm that both Bio-RP and PAGE methods were reliable to obtain the full-length and same sequence of staples. In addition, we constructed two versions of the flexible 1 design using the staples purified with Bio-RP and PAGE methods (Supplementary Figure S12). No distinguishable difference in the trans ratio was observed between these two versions of structures (Supplementary Table S16). The mean and standard deviation of the trans ratio by two methods were 0.7245 and 0.0609 for Bio-RP, and 0.7204 and 0.0272 for PAGE, respectively. Accordingly, we concluded that the trans ratio was hardly affected by the purification method for staples because both methods yield sufficiently high purity.

Supplementary Figures
Supplementary Figure S1. Root-mean-square deviation (RMSD) of DNA oligomers. Each RMSD value was calculated from the minimized structure omitting the four terminal BPs on each end. Supplementary Figure S11. The results of MALDI-TOF using Bio-RP and PAGE methods. We obtained the MALDI-TOF results for representative ten staples (one staple per each twisting block). Table   S19).  Supplementary Table S1. List of MD-simulated DNA oligomers. Each oligomer is 42-BP-long DNA double helix and has the same sequence pattern, 5'-CTGA5-MN-AGTC5-3' / 5'-GACT5-PQ-TCAG5-3'.

Supplementary Tables
MN/PQ represents one of ten regular BP steps. Since MN/PQ or PQ/MN indicate identical one of ten BP steps, we selected notations by comparing M and P in order of A > G > T > C. For example, we chose AG/CT notation in AG/CT or CT/AG case. MnN/PQ or MN/PnQ represent one of sixteen nicked BP steps.    Mean rigidity ratio of each structure was calculated by dividing the summation of multiplying each rigidity ratio and the corresponding number of nicked BP steps by sixty as the total number of nicked BP steps.
A figure below the table shows the mean rigidity ratio values for each structure.