To thread or not to thread? Effective potentials and threading interactions between asymmetric ring polymers

We use computer simulations to study a system of two unlinked ring polymers, whose length and bending stiffness are systematically varied. We derive the effective potentials between the rings, calculate the areas of minimal surfaces of the same, and characterize the threading between them. When the two rings are of the same kind, threading of a one ring through the surface of the other is immanent for small ring–ring separations. Flexible rings pierce the surface of the other ring several times but only shallowly, as compared to the stiff rings which pierce less frequently but deeply. Typically, the ring that is being threaded swells and flattens up into an oblate-like conformation, while the ring that is threading the other takes a shape of an elongated prolate. The roles of the threader and the threaded ring are being dynamically exchanged. If, on the other hand, the rings are of different kinds, the symmetry is broken and the rings tend to take up roles of the threader and the threaded ring with unequal probabilities. We propose a method how to predict these probabilities based on the parameters of the individual rings. Ultimately, our work captures the interactions between ring polymers in a coarse-grained fashion, opening the way to large-scale modelling of materials such as kinetoplasts, catenanes or topological brushes.

From conformations collected in a simulation of this model, we measured Flory's characteristic ratio, C s , [1] where s is the number of monomers in a selected subchain of the polymer, R(s) is the distance between the terminal monomers of the subchain and the average ⟨·⟩ in the numerator is calculated over all subchains of s monomers. Values of C s calculated for s ∈ {1, . . . , 99} are shown in Fig. S1, where we can see that for s ≫ 1, we reach a plateau approximately corresponding to C ∞ = lim s→∞ C s [1]. We note that for the cases βK bend ≥ 20, the characteristic ratio is not properly converged and still drifting at s ∼ 100, nevertheless, we carry out this procedure mainly to give the reader an idea about the length scales, not to calculate C ∞ with high accuracy. In conclusion, the values of C ∞ in Tab. 1 in the main text correspond to the flexiblities of a corresponding ideal chain with the same bending potential as our interacting rings.
To provide an estimate for the persistance length, let us consider a freely rotating chain with C ∞ as measured in Fig. S1. Following Ref. [2] and Ref. [1], it can be shown that and finally arriving at values listed in Tab. 1 in the main text and extended Tab. I in ESI.

B. Minimal surfaces
We compute the minimal surface on each ring separately using Surface Evolver software [3]. The code allows to look for the surface of minimal area given the fixed boundary by moving free vertices of a triangulated trial surface. Here we describe just the essential elements of our implementation of the method, with full details being described in detail in [4].
The surface is initialized as a union of triangles, each defined by two (fixed) vertices coinciding with two successive monomer positions along the ring contour and the third (free) vertex being the center of mass of the ring. The surface is then refined once, by dividing each edge in half by a new vertex, which inherits the property of being fixed if both vertices of the edge were fixed. The new vertices form new edges, subdividing the original triangle into four new ones. Subsequently, the free vertices are progressively propagated proportionally to the total local surface tension force, which is the sum of forces from all the triangles that the vertex belongs to. The force on the vertex v 0 from the triangle formed by three head-to-tail vectors s 0 , s 1 , s 2 labeled in the counterclockwise direction is where the v 0 is the vertex at the tail of s 0 and the surface tension is unity [3]. The total force on a vertex is then multiplied by the step length (i.e. the proportionality constant, typically in the range of 0.1−0.3), which is optimized to reach the minimum area faster. The optimization in essence finds three step lengths that bracket the energy and quadratically interpolates the best one. In some cases this evolution can lead to very thin triangles that then effectively stall the minimization procedure (because the step length is decreasing as a result of the energy bracketing). To avoid that we regularly use vertex averaging (average the position of neighboring vertices), apply equiangulation (a procedure to redefine two adjacent triangles by replacing their shared edge, i.e. the diagonal of the corresponding quadrilateral by the other possible diagonal, to achieve more uniform distribution of the internal angles of the two triangles) and weeding out skinny (if area is smaller than 0.01σ 2 ) triangles (by removing one free vertex of an edge belonging to the skinny triangle and collapsing the two remaining edges to one). The minimization procedure is stopped if the relative area is not changing by more than 0.1% over a course of 240 steps with intermediate equiangulation and vertex averaging (see details in [4]). Note that the topology of our surface is fixed to that of the disc and therefore we do not aim for the true minimal surface, which might have a different topology. We just want to achieve a surface completely contained within the boundary of the ring. We have shown by comparing different minimization methods in [5] that for this polymer model and lengths, the procedure we use here brings us sufficiently close to the true minimal surface if that has a disc topology, for the purposes of defining and quantifying threadings.

C. Threading analysis
To detect whether a given bond vector intersects the area of a given triangle, we use six-dimensional representation of each line segment using the Plücker coordinates. Let us consider a line segment, a, defined by the endpoints ⃗ p = (p x , p y , p z ) and ⃗ q = (q x , q y , q z ). Each Plücker coordinate of this line segment is one of the determinants of minor 2 × 2 matrices of the matrix: Next, we define the side operator s(a, b) as a permutated dot product of the 6D representations of line segments a and b: Finally, to detect the threading, we invoke the side operator three times obtaining {s(x, y 1 ), s(x, y 2 ), s(x, y 3 )}, where x is the bond vector and {y 1 , y 2 , y 3 } are sides of the triangle. If all of the three values s(x, y) have the same sign, the line intersects the triangle [6].
To locate all of the threadings between the two rings, we just repeat this procedure for all possible combinations of bond vectors of ring 1 (and ring 2) and all triangles of ring 2 (and ring 1).

D. Biased sampling & Effective potentials
The effective potential as defined in Eq. 8 is calculated from the pair correlation function, for which we can write where P 0 and P 0 id. are the probability density distributions of the reaction coordinate, r, in the interacting and ideal system respectively. While the latter can be obtained analytically, the former is calculated using the biased (umbrella) sampling and weighted histogram analysis method [7][8][9].
Following the notation of Ref. [10], let us consider a histogram P 0 j composed of M bins indexed as {1, . . . , i, . . . , M } for each of the studied systems. For each system, we carried out S independent simulations indexed as {1, . . . , j, . . . S}. Each j corresponds to a different choice of the r j ∈ {0, 0.5σ, 1.0σ, . . . 30.0σ} in the equation Eq. 7 in the main text.
In the simulation j, we collect N j independent samples of the reaction coordinate, which we sort into the corresponding histogram bins n ij such that M i=1 n ij = N j . Thus we obtain the biased probability distribution P ij , which is related to the unbiased distribution as: where V j is the biasing potential in the simulation j and f j is a normalization factor It can be shown, that the unbiased probability distribution is given by: Eq. S11 together with Eq. S10 pose a set of M +N non-linear equations, where f i and the sought-for P 0 i have to be determined self-consistently. The data are presented for extra mean ring-ring separations, r/R g , split into three panels for better visibility, for ∼ 1000 configurations each.
FIG. S11: Instantaneous values of prolateness of ring A plotted against the corresponding value for ring B from the same configuration for the asymmetric case. The data are presented for two different mean ring-ring separations, r ≈ 0 (green) and r ≈ 2R g,0 (purple) for ∼ 1000 configurations each. FIG. S13: Probability of ring i threading ring j as a function of ring-ring separation. For the fully symmetric case, we show the probability averaged over both permutations (i → j, j → i) and the distance is normalized by the radius of gyration of a single ring at infinite dilution. For the asymmetric cases, we show the probabilities for both possible permutation respectively and a similar normalization of distance is applied, but using the average of the infinite-dilution radii of the rings from Eq. 9 from the main text. The color identifies the system, whereas line style (solid and dashed) in the asymmetric cases differentiates the permutation. The inset shows the probability of mutual threading at r ≈ 0 ring-ring separation.