Influence of Rigidity and Knot Complexity on the Knotting of Confined Polymers

We employ computer simulations and thermodynamic integration to analyze the effects of bending rigidity and slit confinement on the free energy cost of tying knots, ΔFknotting, on polymer chains under tension. A tension-dependent, nonzero optimal stiffness κmin exists, for which ΔFknotting is minimal. For a polymer chain with several stiffness domains, each containing a large amount of monomers, the domain with stiffness κmin will be preferred by the knot. A local analysis of the bending in the interior of the knot reveals that local stretching of chains at the braid region is responsible for the fact that the tension-dependent optimal stiffness has a nonzero value. The reduction in ΔFknotting for a chain with optimal stiffness relative to the flexible chain can be enhanced by tuning the slit width of the 2D confinement and increasing the knot complexity. The optimal stiffness itself is independent of the knot types we considered, while confinement shifts it toward lower values.


INTRODUCTION
While in the macroscopic world it is clear that the effort needed to tie a knot in wire or string will always increase if it is made more rigid, for equivalent microscopic objects, polymers, the same does not hold. Instead, it is found that the free energy cost of knotting a polymer, ΔF knotting , has a minimum at a nonzero stiffness. 1 This finding is particularly interesting in the context of biological macromolecules, such as DNA or RNA, where on the one hand knotting is known to occur 2,3 and has significant effects on key processes, 4−6 while on the other rigidity may depend sensitively on the base sequence, 7−9 leading to varying flexibility along the polymer. Furthermore, there is evidence of correlations between DNA stiffness and sites preferred by type II topoisomerases, 1 enzymes that regulate knotting. 10 It is expected that the rigidity dependence of ΔF knotting will affect the behavior of knots in DNA with nonuniform flexibility, for example by localizing them in regions with favorable stiffness. However, previous work 1 neglected a key qualitative feature of biological DNA, namely that it is typically highly confined. 11−13 Confinement of a knotted polymer in a good solvent may significantly affect its properties. For example, in contrast to three dimensions where they are weakly localized, knots in polymers adsorbed on a surface are strongly localized. 14,15 Considering the properties of polymers confined in a slit, simulations of DNA found a nonmonotonic dependence of the knotting probability on the slit width, 16 and for flexible polymers evidence was found that the particular topology is important. 17 While previous work on knotting in confinement has focused on polymers that have one specific stiffness, here we apply a simple model for a polymer chain under tension to investigate the dependence of ΔF knotting on rigidity for various widths of the geometrical confinement. We find that a local stretching of the chains at the braiding region of the knots is responsible for the fact that the optimal bending rigidity for knot formation, κ min , differs from zero. Geometric confinement, however, pushes this optimal rigidity toward smaller values. The effect of confinement on κ min as well as the amount by which ΔF knotting is reduced for the optimal rigidity κ min depends sensitively on the tension applied to the polymer chain.
The rest of the paper is organized as follows: We first present our model and details about the simulation in section 2. In section 3 we define and explain the observables that have been measured in our simulations. In particular, we define our notion of bending of the polymer chain and establish its connection to ΔF knotting . Section 4 introduces the analysis of the local bending in the interior of the knot, which is carried out to investigate which part of the knot is responsible for the reduction of ΔF knotting for polymers with nonzero bending stiffness κ min relative to a fully flexible chain. We present our results in section 5, whereas in section 6 we summarize and draw our conclusions.

MODEL AND SIMULATION DETAILS
For the polymer chain (linear or knotted), we employ a standard, self-avoiding bead−spring model with rigidity κ, confined parallel to the (x,z)-plane and under tension τ. The interaction part of the Hamiltonian thus reads as In eq 1, r i,j = r j − r i is the vector from bead i to bead j, located at position vectors r i and r j , respectively, with the unit vector rî ,j = r i,j /|r i,j |. The first term represents the bending energy of the chain, where κ is the bending rigidity. The second and third terms are the connectivity and steric terms, respectively, whereas Θ(ω) is the Heaviside step function of ω, which renders the Lennard-Jones potential purely repulsive. We choose ε = k B T, k = 30k B T/σ 2 , and R 0 = 1.5σ, preventing the chain from crossing itself and thus conserving its topology. The chain is confined in a slit parallel to the (x, z) plane, which is realized via a harmonic external potential acting on the ycomponent of the coordinate of each monomer, expressed by the fourth term in eq 1. The last term applies a tension τ on the chain along the z-direction of the setup, with L z denoting the extension of the chain along this direction. We used the LAMMPS simulation package 18 to carry out constant-NτT molecular dynamics (MD) simulations. The polymers consist of N = 256 monomers for the chains simulated at tensions τ = 0.8 k B T/σ and τ = 0.4 k B T/σ and of N = 512 monomers for the simulations at tensions τ = 0.2 k B T/σ and τ = 0.1 k B T/σ. The longer polymer chains at the two smaller tensions are necessary due to the larger knot size for these tension values. The chains are placed in a simulation box with volume V = 100σ × 150σ × L z at the two higher tensions and V = 200σ × 300σ × L z for the lower tensions. The tension τ is realized via a barostat coupled to the fluctuating z-length of the simulation box, while the box lengths in the x-and ydirections are fixed to 100σ and 150σ for the higher tensions 200σ and 300σ for the lower tensions, respectively. The polymer is connected across the periodic boundary conditions in the z-direction to guarantee that the knot is preserved. We also use periodic boundary conditions in the x-direction, while for the y-direction confinement prevents the polymer chain from getting outside the simulation box. For both the thermostat and the barostat, we used a Nose−Hoover chain with 3 degrees of freedom. 19 With m denoting the monomer mass and β = (k B T) −1 , t 0 = (mσ 2 β) 1/2 sets the unit of time. We integrated the equations of motion with a time step Δt = 10 −3 t 0 . The equilibration time was 2 × 10 7 time steps, and data were collected during a total of 3 × 10 8 time steps.

DEFINITION AND INTERPRETATION OF OBSERVABLES
Definition and Physical Interpretation of ΔF knotting (κ). Let F lin,knot (κ) be the free energies of an unknotted and a knotted chain, respectively, for given κ, τ, and chain length. We define ΔF knotting (κ) ≡ F knot (κ) − F lin (κ), a quantity that gives a measure for the effort to tie a knot into the chain of N monomers. In our simulations, we do not calculate the absolute value for ΔF knotting (κ) but its value relative to ΔF knotting (κ = 0) of a flexible chain, a procedure that removes the N-dependence for a linear and a knotted chain of the same degree of polymerization N. Thus, we calculate Ψ(κ) ≡ ΔF knotting (κ) − ΔF knotting (0), which has a direct physical interpretation: Let us consider a long polymer chain with various domains, which differ by their respective bending stiffness κ i and the number of monomers they contain N i . Then, the quantity Ψ(κ) allows us to predict the probability for the knot to be found in the ith domain relative to the probability for it being in the jth domain as An implicit assumption entering eq 2 is that the length of chain segments, N i , is much longer than the knot size N K , so that for most configurations the part of the polymer chain that is affected by the knot is localized in only one of the domains. Because of the applied tension τ, the knotted part of the chain will always remain finite. In ref 1 it was shown that at fixed κ the size of the knot does not scale with N, the number of monomers on the chain, but rather as N K ∼ (k B T/τ) α with some exponent α. The reason for this is that the stretched polymer forms a series of tension blobs, 20 whose size scales with the tension τ but not with N. The knot can only be in one of those blobs; therefore, N K will be finite for any nonzero tension τ. Therefore, for sufficiently large N i , eq 2 indeed provides a good prediction for P i . In ref 1 the prediction of eq 2 was tested for a simulation of a polymer chain with two stiffness domains.
Definition and Interpretation of the Chain's Bending B. With r̂j ,j+1 denoting the unit vector between monomers j and j + 1, we define the quantity for any configuration of its monomers and call it the bending of the chain; evidently it holds that B̂≥ 0. One can check that this definition of the bending is sensible for various special cases. For instance, for a configuration with a straight polymer chain this definition gives the minimum value for B, namely B̂= 0. The contribution to V({r i }) defined in (1) due to the bending stiffness is κB. Moreover B T (κ) ≡ ⟨B⟩ is the thermodynamic expectation value of the same for a chain of topology T ∈ {knot, linear}.
With F T (κ) denoting the free energy of the polymer, it holds that ∂F T (κ)/∂κ = B T (κ). This allows us to calculate how the cost of knotting changes with bending stiffness κ: The existence of a minimum of Ψ(κ) for κ ≠ 0 will depend on the sign of the slope of Ψ(κ) at κ = 0, ΔB(0). If ΔB(0) < 0, we expect an optimal knotting rigidity κ min ≠ 0, whereas we anticipate a monotonically increasing function Ψ(κ) in the opposite case, ΔB(0) ≥ 0. This is a reasonable expectation, since we will always obtain ΔB(κ) > 0 for sufficiently stiff chains. Indeed, for βκ ≫ 1, a linear polymer will adopt an almost straight configuration with B̂≈ 0, as configurations with nonzero B̂are penalized by a high bending energy. B̂= 0 is unique for the straight configuration, which is of course not knotted. Therefore, ΔB(κ) > 0 will hold for all knots in the βκ ≫ 1 regime.
To illustrate this further, let us for example consider a stiff polymer chain with a trefoil knot. Its minimal energy is obtained for a straight polymer chain with an approximately circular domain at the location of the knot. The circle is tangent to the braid point and contains N K (κ) ≅ [2π 2 κ/(τl)] 1/2 monomers, where l ≅ σ is the bond length. 21 Accordingly, ΔB(κ) ≅ (2π 2 τl) 1/2 /κ > 0 in this limit.

ANALYSIS OF THE LOCAL BENDING IN THE KNOTTED DOMAIN
Having in mind the goal of localizing which part of the polymer in the vicinity of the knot is contributing to an increased or decreased average bending and therefore to a nonvanishing ΔB(κ), we need to determine the knotted domain on the polymer chain. We define the knotted state of open subsections by introducing a topologically neutral closure scheme, which transforms an open string into a ring polymer that has a mathematically well-defined topological state. 22−24 We used a scheme where the end points of the open polymer are connected to a sphere at infinity in the direction of the vector from the centroid to the respective end point. The knotted part of the polymer chain is then the smallest domain for which the closure yields a ring polymer that has the correct Alexander polynomial.
Once the ends of the knot have been identified, we introduce a new enumeration scheme for the monomers, denoted by the Greek integer index α, which can have positive as well as negative values. The two monomers in the interior of the knot that lie n bonds away from the end points obtain the index α = −n, whereas the two monomers to be found n bonds away at the exterior of the knot are assigned the index α = n; accordingly, α = 0 for the two end points of the knot. There exist, thus, for every value of α two position vectors r α j , j = L, R, where L/R denotes whether the monomer is at separation α from the left/right end point of the knot. Accordingly, we define the local bending contribution from the two monomers carrying the index α as For α < 0, bα measures the local bending of angles in the interior of the knot, while for α > 0, bα measures the local bending outside of the domain that was identified as knotted. The domain limits for α ∈ [α min , α max ] depend on the instantaneous configuration, as the number of monomers on the knot N K determines how negative α can become. The range of α is constant, α max − α min + 1 = N/2. We introduce a characteristic function χα = 1 or 0 depending on whether the index α occurs for a given conformation or not and define the local bending difference between a knotted and a linear chain as  In eq 6, b linear (κ) ≡ 2B linear (κ)/N is the thermodynamic average of the bending of two angles on a linear chain. It follows that where M = N/2 is the smallest value guaranteeing [α min , α max ] ⊂ [−M, M] for all polymer configurations.

RESULTS
Results for the local bending of unconfined polymers are shown in Figures 1a,b for two different values of the applied tension. As can be seen, the quantity Δb α (κ) vanishes within less than 10 beads outside the knot; outside this region, a knotted chain hardly differs in its bending from an unknotted one. According to eqs 4 and 7, the area under the curves in Figures 1a,b determines the slope or Ψ(κ) at any κ value. The major negative contribution to ΔB(κ) for κ = 0 comes from a small number of monomers close to the knot ends, i.e., in the braiding region of the knot, in which monomers start getting into the knotted domain. The bending suppression is therefore related to the interaction of the different strands at the crossings, which effectively confines the strands and hence reduces their random bending. Thus, the negative slope of ΔF knotting (κ) at κ = 0 arises from this additional straightening of the knotted chain with respect to its flexible, unknotted counterpart. As ΔB(κ = 0), which is equal to the area under the respective curves in Figures 1a,b, is negative, a flexible knotted chain is, on average, less bent than its unknotted counterpart, contrary to the intuitive expectation that knotting inevitably increases the total bending of a polymer. In Figure 2 we show a simulation snapshot of the knotted domain of a flexible chain, in which the parts of the molecule that get "straightened out" due to the knot are highlighted. As κ grows, we eventually reach the intuitively expected regime in which knotting increases bending; see, e.g., the curve for κ = 20 k B T in Figure 1a, for which ΔB(κ) > 0.
Comparing the curves of Figures 1a,b, one sees that ΔB(κ) for κ = 0 is more negative for τ = 0.8 k B T/σ than for τ = 0.1 k B T/σ. This is due to the fact that at smaller τ the braid region is looser and the bending suppression due to the different strands at the crossings is reduced. On the contrary, if we increase κ, we arrive at a regime where ΔB(κ) is more negative for smaller tensions τ. At βκ = 3.0, ΔB is already positive for τ = 0.8 k B T/σ, due to the positive contribution of Δb α (κ) in the interior of the knot, which arises as the knot enforces the polymer to form a loop. The bending throughout a loop is larger if the loop is smaller. For τ = 0.1 k B T/σ the knot size is increased, which is the reason why the positive Δb α (κ) contribution in the interior of the knot is significantly smaller than for τ = 0.8 k B T/σ. Accordingly, for lower tensions the negative net result for ∂Ψ(κ)/∂κ persists for higher κ-values than for higher tensions. As can be seen in Figure 1b, we see a reversal of ∂Ψ(κ)/∂κ from negative to positive values only at a rigidity as high as βκ ≅ 9.
We now turn our attention to the confined case. It turns out that the effects of confinement are most transparent at a small, but nonzero, value of the bending rigidity; thus we show in Figures 1c,d results for Δb α (κ) for βκ = 0.5, which is characteristic for all values of κ ≤ κ min , the latter being the value of the rigidity for which Ψ(κ) attains its minimum value Ψ min . Here, striking differences between the tensions τ = 0.8 k B T/σ and τ = 0.1 k B T/σ show up. While for τ = 0.8 k B T/σ there is hardly any difference between the confined and unconfined polymers for slit widths as small as d = 5σ, for τ = 0.1 k B T/σ, confinement enhances Δb α at the beginning of the knot by almost a factor 2. For both tensions, the suppression of bending through knotting becomes even stronger as a result of the geometric constraints, and thus ΔB(κ) is more negative for the polymer in the slit than it is for the free polymer. This is consistent with the interpretation given above for the case of the unconfined polymer, since the slit confinement forces the strands at the braid region to come closer which further reduces the random bending in the braid region. This effect, however, is more pronounced for smaller tensions, where the braiding region without confinement is looser than for polymer chains under higher tensions.
There exists a correlation between the slopes of Ψ(κ) and N K (κ) shown in Figure 3. In the high-κ domain,Ψ(κ) ∼ N K (κ) ∼ √κ, as is evident from the discussion following eq 4. For κ = 0, the knot is swollen due to the presence of steric interactions, maximizing in this way its entropy. However, for nonzero κ ∼ k B T, the fluctuations of the monomers are restricted in the first place, enabling thus a tighter braided region and a concomitant reduction of knot size. Thus, also for small κ the slopes of Ψ(κ) and N K (κ) are expected to have the same sign. Note, however, that the value κ̃that minimizes N K (κ) does not coincide with κ min , e.g., for very strong confinements κ min ≠ 0, whereas κ̃is, within simulation resolution, vanishingly small.
The dependence of the rigidity κ min for which Ψ(κ) has its minimum on the degree of confinement for different applied tensions is summarized in Figure 4a. We find that for a tension of τσ = 0.8k B T κ min is only affected by confinement for slit widths lying at the monomer scale. In this regime of ultrastrong confinement, the energy cost that chain segments would have to pay to go one above the other in a gradual fashion at the braiding regions are too high. This is caused by the external potential, which assigns an increasingly high energetic cost for every monomer that deviates strongly from the y = 0-plane. Accordingly, it is preferable for the system to form localized "kinks" of one or two monomers in the braiding region, which expose a minimal number of monomers to the regions of high external potential, while at the same time creating strong bending there.
The situation is quite different, however, for lower tensions. In this case, also a moderate confinement of the order of 10 bond lengths can significantly affect the value of κ min . To better quantify the effects of confinement, we employ a simple, roughand-ready separation of the data points shown in Figure 4a into two groups: for high values of d, the points form plateaus at the bulk values of κ min , which we connect by horizontal lines.
Through the other groups of points straight lines are drawn by hand, which intersect the horizontal ones at tension-dependent crossover confinement widths d × (τ). These values denote, by construction, the crossover of the behavior of κ min from bulkdominated, for d > d × (τ), to confinement-affected, for d <  d × (τ). The results are summarized in Figure 4b, where it can be seen that d × (τ) is significantly increased for lower tensions. As we discuss below, the reason the situation is strikingly different for lower tensions seems to be related to the fact that the knot size is then significantly increased with respect to higher tensions.
The increased effects of confinement as τ decreases are also manifested on the value of Ψ min as well as on the knot size N K (κ min ). The former quantity is shown in Figure 5a and the latter in Figure 5b. As can be seen in Figure 5a, and in contrast to the values κ min itself, even in the case of higher tension the corresponding depth of the minimum Ψ min is influenced by a confining slit width of the order of 10 bond lengths. However, for lower tensions, the effect on Ψ min is felt at even larger slit widths. Furthermore, the difference between the value Ψ min in the bulk (d/σ ≫ 1) and the one for the optimal slit width is enhanced. The fact that the effect of the confinement on the knot size is more pronounced for smaller tensions, as shown in Figure 5b, correlates well with the finding that confinement shifts κ min to lower values for sufficiently small tensions. As we have seen in Figures 1a,b one of the contributions that eventually render ∂Ψ(κ)/∂κ positive is the bending in the interior of the knot, which arises as the knot enforces the polymer to form a loop. This contribution is larger for smaller knot sizes, and it is therefore consistent that κ min will be shifted by confinement if the latter is able to significantly reduce the knot size.
Up to now, all results have been derived for the simplest, trefoil knot; real polymers can, however, display a large variety of increasingly complex knots. 25,26 Considering other knots allows us on the one hand to put the general character of our results to the test and also to corroborate our assertion that the crossings at the braiding region are responsible for the reduction of Ψ(κ) at finite κ-values. Indeed, more complex knots have more crossing points, where the strands of the polymer chain interact with each other. Thus, according to our analysis above, one should expect that a more complex knot will lead to lower values of ΔB(κ) and Ψ min for small but finite κ. Our findings for different knot topologies (denoted in the Alexander−Briggs notation 27 ) are summarized in Figures 6a,b. The data in Figure 6a confirm that the main effect of the increased knot complexity is the addition of crossing points which all result in a similar bending suppression as the crossing points of the trefoil knot. Accordingly, for the knot topologies investigated, Ψ min is approximately proportional to the number of minimal crossings of the respective knot diagram, at least for confinements with d ≥ σ. It is also striking that for the unconfined case Ψ min is, within error bars, identical for the 5 1 and 5 2 topologies. A test of whether Ψ min is in good approximation proportional to the minimal number of  crossings of arbitrarily complex knots is beyond the scope of this work. However, the number of strand crossings will increase with knot complexity. Accordingly, the free energy penalty for putting a knot on a stiff polymer (κ ≠ 0) can be much lower than the one for putting it on a fully flexible polymer (κ = 0), by amounts that grow with the knot complexity. Whereas Ψ min is sensitive to the knot type, κ min is not, as can be ascertained from the results shown in Figure 6b. All data fall within a narrow band of width Δκ min ≅ 0.2k B T irrespective of the knot topology.
The tensions considered in our work are of the order of k B T/ σ. At room temperature and a monomer length scale of 1 nm this corresponds to the pN scale. These tensions result in values of κ min ≈ 5 k B T. As it was found in previous work, 1 without confinement κ min scales approximately as ∼τ −1/2 . A reduction of the tension down to the fN scale, which is typical of doublestranded DNA molecules, 28 will bring κ min at the order of 100 k B T. Our results imply that at these lower tensions the effect of confinement on κ min can be expected to be even more pronounced.

CONCLUSIONS
In summary, we have demonstrated that the local stretching at the braiding region and close to the crossing points is the physical mechanism responsible for the minimization of the free energy penalty of knotting of a linear polymer for nonvanishing values of the bending rigidity. Confinement can affect the location of the optimal rigidity for sufficiently low tensions, when it is at the same time significantly affecting the knot size. We therefore expect the geometrical reduction of dimensionality to become relevant for the location of the knots of chains with variable rigidity if the latter are under sufficiently small tensions. For tensions at the fN scale, which are typical of double-stranded DNA molecules, 28 we therefore expect that confinement to strongly influence the value of the optimal rigidity. On the other hand, the amount of reduction of the knotting of free energy by rigidity strongly depends on the topology of the knot, and it increases with the knot complexity, scaling roughly with the number of minimal crossings of the knot. Accordingly, we anticipate that more complex knots will localize more strongly in the optimal regions of a chain than simpler ones. Recent advances in tying knots on polymers by optical tweezers 29,30 and adsorbing them on mica surfaces 15 should allow for experimental testing of our predictions.

Notes
The authors declare no competing financial interest.

■ ACKNOWLEDGMENTS
This work has been supported by the Austrian Science Fund (FWF), Grant 23400-N16.