Knot formation of dsDNA pushed inside a nanochannel

Recent experiments demonstrated that knots in single molecule dsDNA can be formed by compression in a nanochannel. In this manuscript, we further elucidate the underlying molecular mechanisms by carrying out a compression experiment in silico, where an equilibrated coarse-grained double-stranded DNA confined in a square channel is pushed by a piston. The probability of forming knots is a non-monotonic function of the persistence length and can be enhanced significantly by increasing the piston speed. Under compression knots are abundant and delocalized due to a backfolding mechanism from which chain-spanning loops emerge, while knots are less frequent and only weakly localized in equilibrium. Our in silico study thus provides insights into the formation, origin and control of DNA knots in nanopores.

www.nature.com/scientificreports/ simulation studies have been further fueled by a recent experiment that demonstrates that knots indeed occur in compression experiments 16 . In computer simulations studies Orlandini and Micheletti have already investigated equilibrium knot formation of coarse-grained DNA models in nanonchannels 31,32 . Of particular relevance to our work is a recent study 38 in which the non-equilibrium formation of knots and so-called geometrical entanglements as measured by counting crossings under projections are investigated in closed nanochannels exposed to compression and relaxation cycles. It was found that the two types of entanglements evolve with different dynamics and are for the most part uncoupled.
Here, we investigate knot formation when a confined dsDNA is being pushed by a nano-dozer in a nanochannel whose width is much smaller than the contour length of the dsDNA. An important difference from previous studies is that our system is open ended in one direction and that we study the evolution of knots in a constantly moving steady state. We also vary the bond stiffnesses to investigate the influence of changes in salt conditions. The key result is that the confined chain in the nanochannel pushed by a nano-dozer will progressively become highly knotted with delocalized knots. The knotting probability is greatly enhanced compared to corresponding equilibrium simulations, which in addition to compactification can be traced back to a backfolding mechanism for semi-flexible chains. Next, we describe the model, some essential facts about the LD simulation scheme, how our coarse-grained chains can be mapped onto DNA and the method that we use to analyze knots.

Results
Emergence of knots in nanochannels. Figure 1 summarizes the main findings of our study. Applying a pushing force leads to a compactification of the polymers (Fig. 1a), which in turn dramatically increases the occurrence of knots in the steady state in comparison to equilibrium values (Fig. 1b). Likewise, the amount of trefoil knots is reduced for configurations with a higher total knotting probability because the high density induces the formation of multiple or more complex knots. In this compact state, knots are delocalized and span over the whole chain as indicated for the example of trefoil knots in Fig. 1c where for a bending stiffness κ > 20 the average length of the knot is approximately 80% of the contour length or higher, which implicates that knots are formed preferentially near each end (please see Fig. 2e), while knots in equilibrium conformations are significantly smaller. These findings in a sense mirror previous observations, e.g. in Ref. 29 , which demonstrated that a θ-transition from a swollen coil to a globular state is not only accompanied by an increase in knotting but also by a delocalization of the latter. Figure 1d investigates the influence of the piston velocity for the experimentally relevant case of κ = 4 (DNA in a nanochannel, see "Methods" section). Again, compactification with increasing velocity is directly related to an increase of overall knotting. These results suggest that the occurrences of knots can be tuned by the speed www.nature.com/scientificreports/ of the piston and converge towards the equilibrium values for small piston velocities (Fig. 1e). This outcome is similar to results of Michieletto et al. 38 indicating that an increased piston force will lead to an increase in the overall knotting probability and knot complexity as shown by a decrease in the occurrence of simple trefoil knots. As indicated above the decrease of knotting towards the equilibrium state at slow piston velocities is again accompanied by a trend towards a weak localization of trefoil knots (Fig. 1f) 29 . Figure 2 sheds light on these findings from a molecular basis. For κ = 4 the structure is disordered but the position of the monomers is still correlated with their sequence as indicated by the color scheme in Fig. 2a and the bead positions in Fig. 2c. For κ = 20 , the persistence length already exceeds the width of the tube which in conjunction with compactification leads to backfolding (Fig. 2b,d). The backfolding on the other hand creates loops which are a prerequisite for knots and in turn explains the initial rise in knotting with κ as well as their delocalization. In this context it would be interesting to study if backfolding creates a prevalence of torus-type knots as e.g. observed in the DNA located in viral capsids 4,7 . Unfortunately, the statistics of our simulations do not allow for a meaningful comparison. For large persistence lengths, backfolding becomes more difficult resulting in a lower knotting probability (Fig. 1b). In the equilibrium case, the compactification from the piston is no longer present and for κ = 20 , the chain can already spread throughout the channel which leads to a low knotting probability and weakly localized knots (Fig. 2e).

Discussion
In this manuscript we investigate velocity induced knot "production" in a nanochannel in comparison to those under equilibrium conformations. Both knotting probability and knot sizes depend strongly on piston velocity and resulting compactification as well as chain stiffness which can be, e.g., mitigated by adjusting ionic conditions and screening of charges in DNA. We observe that if the chain's persistence length is greater than the width of the nanochannel, knots form by a backfolding mechanism. Since backfolding becomes harder for larger stiffness, the probability of knot formation decreases which explains the observed non-monotonic characteristic of knot formation in a nanochannel. We also study relative occurrences of complex knots as a function of the piston velocity and the chain stiffness. Our study thus sheds some new light on recent experiments in which DNA knots were created in a flow channel 16 and provides insight on the molecular origin and control of self-entanglements under these conditions. Finally, we would like to point out that the coarse-grained simulation does not include hydrodynamic effects. It is worth considering how the results will change if we had incorporated it in the simulation. Dorfman has argued that for flexible chains hydrodynamic interactions are important. But for the semi-flexible chains with persistence length larger than the pore width, the chain is fully extended and is described by the free-draining limit 50 . Thus, for most channel sizes which result in a significant extension of the DNA compared to its bulk conformation, the hydrodynamic interactions between segments of the chain are almost completely screened. For the parameters used here the chain conformation lies in the transition region between deGennes blobs and Odjik limit and there is no theoretical argument for the effects of the hydrodynamic interaction in this regime. However, from Fig. 2 we observe that the chain conformations are mostly extended and hydrodynamic effects are likely to be small. For large velocities folded conformations are very different from deGennes blobs and therefore, one would expect the hydrodynamic effects to be small also, and the conclusions of this manuscript will essentially remain the same. www.nature.com/scientificreports/

Methods
Coarse-grained polymer model. The coarse-grained polymer model for LD simulations used here is exactly the same as in our previous publication 47 where a bead-spring model polymer chain is confined to an open-ended rectangular channel and pushed from the right with a piston in the negative x-direction (Fig. 3a). The semi-flexible chain (Fig. 3b) is represented by a generalized bead-spring model 51 where the beads (monomers) interact via excluded volume (EV), a Finite Extension Nonlinear Elastic (FENE) spring potential and a bond-bending potential enabling variation of ℓ p as implemented previously 46,47 . The excluded volume interaction between any two monomers i and j of diameter σ is given by a short range truncated and shifted Lennard-Jones (LJ) potential U LJ (Eq. 1) of strength ε with a cutoff distance r c = 2 1/6 σ is given by: where r ij = |� r i − � r j | is the distance between any pair of beads. Successive monomers are connected by a FENE spring potential where k FENE is the spring constant and R 0 is the maximum allowed bond length. The parameters k FENE and R 0 along with ε and σ determine the bond-length. The chain stiffness is controlled by a bond-bending potential is the angle between two successive bond vectors � b i−1 = � r i − � r i−1 and � b i = � r i+1 − � r i , respectively, as shown in Fig. 3b. In three dimensions, for κ = 0 , the persistence length ℓ p of the chain is related to κ via 39,52,53 for the values of κ considered in this work where k B is the Boltzmann constant and T is the temperature.
Langevin dynamics simulation. We use the following Langevin dynamics equation of motion to advance the position of the ith monomer (1) www.nature.com/scientificreports/ where γ is the monomer friction coefficient, and W i is a Gaussian random force with zero mean at temperature T which satisfies the fluctuation-dissipation relation. The numerical integration is implemented using the algorithm introduced by Gunsteren and Berendsen 54 . Our previous experience with LD simulations suggests that appropriate parameter specifications are γ = 0.7 mε/σ 2 , k FENE = 30ε/σ , R 0 = 1.5σ , and a temperature k B T/ε = 1.2 . For a time step �t = 0.01τ (with τ being the standard Lennard-Jones time) these parameter values produce stable trajectories over a very long period of time and do not lead to an unphysical crossing of a bond by a monomer 55,56 . The average bond length stabilizes at b l = 0.970 ± 0.002 with negligible fluctuation regardless of chain size and rigidity 55 . The piston is moved with a constant velocity of v 0 = 0.005 σ τ if not noted otherwise after an initial equilibration of the chain. We ensure that the MD time for the pushing phase is long enough for the chain to attain a steady state shown in Fig. 3c-e that displays the connection between chain extension (Fig. 3c,d) and knot formation (Fig. 3e) in approach to the steady state. Times to reach the latter depend on bond stiffness κ as seen from the behavior of R 2 g (t) in Fig. 3c,d. While for κ = 4 reaching a steady state takes less than 50,000τ , it takes around 160,000τ for κ = 20 . For each κ and v, physical quantities are averaged over at least ten independent runs. Our analysis indicates that knots can form and dissolve both during the initial compression and in the steady state. Even compact structures can still change their knot type in accordance with other results on knot formation under applied force 38 . Of course, the probability of forming a knot, however, is significantly larger in the compressed (steady) state.
Reptation Monte Carlo simulation to study the equilibrium limit. Note that piston speeds in coarse-grained implicit solvent simulations are typically orders of magnitude faster when compared to experiments. Therefore, we have also undertaken reptation Monte Carlo simulations of a slightly simplified model system of a single semi-flexible bead-spring chain confined inside a rectangular nanochannel of fixed size. In the reptation move which resembles the movement of a slithering snake, one segment is deleted from a randomly chosen chain end and attached to the other end 57 . Moves are accepted based on the Metropolis criterion. The repulsive Lennard-Jones and bond-bending potentials were matched with those of the LD simulation as described above. However, contrary to the LD simulation model's confinement potential imposed onto the tube walls we use non-interacting walls and fixed bond length b l = 0.967 . These MC simulations allow for a comparison of our dynamical investigations with equilibrium values (corresponding to piston velocity v → 0).

Knot analysis.
Knots in a closed chain are typically characterized by the minimum number of crossings observed when projecting a 3D chain onto a plane and can be considered as a fine gauge for the overall structure. Apart from the unknotted ring, the so-called unknot, the simplest knot is the trefoil (3 1 ) knot, which contains three crossings. There is one knot type with four crossings (4 1 ) and two with five crossings, and from there on the number of different knots with the same number of crossings increases exponentially 58,59 . In our setup the polymer chain is open, and therefore, a closure connecting both ends of the chain has to be defined. First, we connect the end-points of each polymer with its center of mass. Along these lines we define a closure which emerges from one terminus follows the first line connects to the second one far away from the polymer and ends at the second terminus 40 . After closure, the Alexander polynomial can be determined as described in detail in 60 (compare Fig. 3b). Knot sizes are determined by successively removing monomers from the ends of a polymer until the knot type changes 29 .
Mapping onto DNA and comparison with experiments. Mapping our semi-flexible chain onto DNA is based on Eq. (4). For simplicity, we assume a solvent-independent persistence length of 50 nm or 150 base pairs. Furthermore, we assume that our beads describe the locus of a double-stranded DNA strand. In high salt conditions (1M NaCl), charges of DNA are completely screened and σ ≈ 2.5 nm. In physiological conditions charges are only partially screened and σ ≈ 5 nm, and for low salt conditions σ increases even further to about 15 nm at 0.01 M NaCl 33,39,61 . With a simulation temperature of T = 1.2 used throughout we obtain (in simulation units) κ = 24 for high salt, κ = 12 for physiological and κ = 4 for low salt conditions. This allows us to put our simulations in the context of recent experiments by Amin et al. 16 undertaken at an estimated ionic strength of 8 mM which corresponds to our low salt scenario. Our chain has a contour length of L = Nσ = 1024σ = 15,360 nm or 46,080 base pairs, while our confining tube has a width of 16σ ≈ 240 nm. This compares to 168,903 base pairs and tube dimensions of 325 × 415 nm used in Ref. 16 . Note that the mapping changes drastically with ionic conditions.
The time scale of the simulated quantities can be translated into an experimental time scale via For low salt conditions we use σ = 15 nm as stated above, a mass m of 618u per base pair, a persistence length of κσ k B T = 3.33 beads and therefore 45 base pairs per bead. We assume room temperature of T = 300 . This results in a time scale of 1 simulation time which equals to approximately 1.6ns. Our simulation time of 400,000τ for κ = 4 is therefore equivalent to 6.4 × 10 −4 s. Hence without explicit solvent one obtains a much faster time scale compared to experiments which often take place at a scale of seconds 16 . A typical experimental piston velocity for these experiments is 0.1-1 µm/s 14,15 . In the simulation we advance the piston with a velocity v 0 = 0.005 σ τ ≈ 0.05 m/s, several order of magnitude faster than the experimental piston speed. Thus dynamics in our coarse-grained (5) mr i = −∇(U LJ + U FENE + U bend + U wall + U piston ) + γ v i + W i ,