Persistence Homology of Entangled Rings

Here we analyse MD simulations reported in previous works (the interested reader may find more details in Refs. Michieletto and Turner (2016); Michieletto et al. (2017)). Briefly, the polymers are modelled as Kremer-Grest chains made of beads of size $\sigma$. The beads are connected along the backbone by finite-extension-nonlinear-elastic (FENE) bonds and interactions between beads are governed by the Weeks-Chandler-Andersen (WCA) potential (purely repulsive Lennard-Jones). We also introduce polymer stiffness by including a Kratky-Porod term acting on triplets of neighbouring beads $U_{\rm bend}=(k_{B}Tl_{p}/\sigma)(1+\cos{\theta})$, where $\theta$ is the angle formed by consecutive bonds and $l_{p}=5\sigma$.

The motion of the beads is evolved via a Langevin equation that couples the motion of the beads with an implicit solvent. The damping constant is set to $\gamma=\tau_{MD}^{-1}$ where $\tau_{MD}=\sqrt{m\sigma^{2}/\epsilon}$ is the microscopic characteristic time. The thermal noise satisfies the fluctuation dissipation theorem and to integrate the equations of motion we employ a velocity-Verlet algorithm with timestep $dt=0.01\tau_{MD}$.

Using configurations from MD simulations as an input data, the PH analysis generates a diagram called the persistent diagram (PD), from which we can extract hidden topological features Hatcher and Press (2002); Hiraoka et al. (2016); Nakamura et al. (2015). For our present purpose, let us consider the $m$-th ring ($m\in[1,M]$); the input data is $\{{\vec{r}}^{(m)}_{n}\}=({\vec{r}}^{(m)}_{1},{\vec{r}}^{(m)}_{2},\cdots,{\vec{r}}^{(m)}_{N})$, where ${\vec{r}}^{(m)}_{n}$ is the position of $n$-th monomer in the $m$-th ring. Associated to each monomer’s coordinate we assign a ball of radius $r(\alpha)=\sqrt{\alpha}$, which should not be confused with the real radius $\sigma/2$ of the monomers. Initially, we set $\alpha=0$, so the input data is regarded a collection of volumeless points. This method could be generalized to copolymer setting a different initial radius for different types of monomers. We then gradually inflate the balls and for each value of $\alpha$ we connect all the balls that are even partially overlapping. This algorithm generates sets of topological structures made of connected balls that evolve with the balls’ radius $r(\alpha)$.

Persistent homology is a method to compute the topological features of a system at different length scales Hatcher and Press (2002); Hiraoka et al. (2016); Nakamura et al. (2015). The main concepts in this type of analysis are: the persistent homology point, that defines the size at which a certain topology of the system appears/vanishes (birth/death scale), and the persistent diagram (PD), a collection of all the PH points of fixed dimension (see Materials and Methods). Our analysis is focused on entangled solutions of unknotted, unconcatenated ring polymers; in order to capture the constraints related with the ring topology, we consider the homology of dimension 1 that identifies the formation of “loops”. Importantly, here and in what follows, we refer to a “loop” as a one-dimensional (1D) geometric object detected in the PH analysis, while the term “ring” is reserved for the full polymer contour length (for more details see Materials and Methods).

An example of this analysis is shown in Fig. 1. Starting from a single snapshot from MD simulations of ring polymers – made as sterically interacting connected beads (see Materials and Methods) – we associate at each bead position a sphere of size $r(\alpha)=\sqrt{\alpha}$. When $\alpha<\alpha_{min}$ ($\simeq 0.3$ in the example of Fig. 1Bi) there is no loop present because there are no overlapping spheres. At $\alpha=\alpha_{min}$ all the balls are connected with the consecutive along the ring and therefore the ring itself is identified as a loop in the PD analysis (Fig. 1Bii). It should be mentioned that the first loop to be identified by the PH algorithm is not necessarily the full ring contour length and may instead take up only a fraction of the whole ring, in agreement with previous findings using “contact maps” Michieletto (2016); Likhtman and Ponmurugan (2014). For increasing values of $\alpha$, certain existing loops are annihilated while others are formed (Fig. 1C, D) and, eventually – at $\alpha\simeq 96$ for this particular example – the inner space of the ring is filled in and no loop is detected.

If a loop of connected balls exists for a given $\alpha$, a small inflation may lead to a change of topological structures yielding disappearance of such a loop and appearance of others. The sequence of changes in these 1D topological features for different $\alpha$ is intimately associated to the hierarchical, multi-scale folding of the rings, and is encoded in the PD generated by the PH algorithm. Below, we explain how to translate the PD into useful information about the folding and topological interactions between rings in solution.

For a given ($k$-th) loop on the $m$-th ring, i.e. $c_{k}^{(m)}$, its birth and death scales are denoted as $b_{k}^{(m)}$ and $d_{k}^{(m)}$ and capture the length-scales at which the loop first appears and then disappears, respectively. It should be highlighted that the birth scale is also related to the maximum spatial distance between two monomers that generate a certain loop; indeed, when $\sqrt{\alpha}$ is larger than the distance of two monomers these are then “in contact” and a loop is created by the connected portion of contour length between these two monomers. The minimum length to achieve this is $\sqrt{b_{k}^{(m)}}$. The death scale is instead related with the geometric 3D size of such a loop. More precisely, $\sqrt{d_{k}^{(m)}}$ corresponds to the radius of the maximum sphere that can be drawn inside the loop $c_{k}^{(m)}$.

The PD of the $m$-th ring is made by the collection of all of the birth and death scales $(b_{k}^{(m)},d_{k}^{(m)})$ of all the loops generated by ranging $\alpha\in[0,\infty)$, i.e. $PD(m)\equiv\left\{(b_{k}^{(m)},d_{k}^{(m)})\in\mathbb{R}^{2}|k=1,2,\dots\right\}$. This set of points is typically represented in the top half triangle of a two-dimensional plane as in Fig. 1C.

For each point ($b_{k}^{(m)},d_{k}^{(m)}$) in the PD (see Materials and Methods), we define its life span $l(c_{k}^{(m)})\equiv d_{k}^{m}-b_{k}^{m}$, which measures the distance from the diagonal line to the point. The points that lie far from the diagonal, i.e. with large life span, are the ones reflecting robust and persistent features of the system, whereas the ones lying close to the diagonal are often associated with noise Kusano et al. (2018). Among other off-diagonal points, we call the loop associated to the point with the largest life span “primal”. More specifically, for open non-crumpled conformations, the primal loop may corresponds to the ring itself and its birth scale is expected to be on the order of $\simeq(\sigma/2)^{2}$ ($\sigma$ is the size of a simulated bead). In general, the primal mirrors the presence of a large opening of the double-folded structure in the ring’s conformation and its life span quantifies the size of the opening.

For example, in Fig. 2 we show two typical conformations and their respective PDs in red (ring 17) and blue (ring 18). The primal loop of the red ring has a large death scale and its life span is much longer than that of the blue ring (notice the logarithmic scale in the plot). This difference reflects the fact that while ring 17 displays a large opening, ring 18 assumes a more crumpled and double folded conformation entailing a much shorter primal life span (see the snapshot in Fig. 2A). Another notable feature of Fig. 2A is that ring 18 displays more points away from the diagonal than ring 17 (compare, e.g., the interval $[0.5,5]\in$ birth scale). This feature of the PD reflects the crumpled nature of the conformation: the more crumpled is a ring, the more balls of a certain radius form bridges between two distal segments of the ring contour in turn generating a moderately robust loops.

Focusing our attention on the primal loop, it is possible to study the average linear size $\sqrt{d_{PL}}$ of the loop at different ring lengths, where $d_{PL}$ is the death scale of the primal loop. In Fig. 2B, we show that this quantity follows a power law dependence similar to the radius of gyration. This result agrees with the existing theories Grosberg (2014); Ge et al. (2016); Sakaue (2018, 2011); Halverson et al. (2014b); Cates (1987), that predict an exponent $1/2\leqslant\nu\leqslant 1/3$, with $\nu=1/2$ for the short ring length, an intermediate regime with $\nu=0.4$ Cates (1987) and a compact long ring regime $\nu=1/3$ Grosberg et al. (1993). In the inset is reported the probability distribution of $\sqrt{d_{PL}}$ for different ring size, to show that not only the average but also the peak value of the distribution is shifted towards larger sizes. For comparison, we show the probability distribution of the radius of gyration $R_{g}$, and its correlation with $\sqrt{d_{PL}}$ in Appendix A.

Having identified features of the PH analysis that characterise the conformation of single rings, we now aim to employ this tool to detect and quantify inter-ring interactions called threadings Michieletto et al. (2014a) for which an unambiguous and unique (i.e. without arbitrary and free parameters) identification algorithm is needed.

Let’s consider two unconcatenated and unlinked rings in 3D space; it is natural to define the arrangement of the two rings depicted in Fig. 3A as one in which the black ring is “threading” the other by crossing twice the surface spanned by the closed blue curve. Bending the blue ring into a baseball-like shape (Fig. 3B), the black ring still seems to pierce twice through the surface. In contrast one may argue that there is no threading in Fig. 3C, where the arrangements of the curves is the same but the spanned surface is forming a gulf that is not pierced by the former curve. Given the symmetry of the baseball seam curve, which side is selected as a spanning surface is a subtle matter. Such a subtlety would be even more enhanced for 3D crumpled rings. Additionally, while the situations in Figs. 3B and D appear different, they can be continuously mapped to the same three dimensional arrangement (as Fig. 3D); furthermore, the physical hindering of the black curve on the dynamics of the blue may be similar. Finally, in both cases, the obstacle can be removed by a continuous deformation, for instance by sliding it to left direction, as threadings are not conserved under isotopy. Note that the figures Figs. 3A-D are a typical example of isotropic change, and there can be other examples with similar ambiguity when we define threading in an intuitive way.

Nevertheless, here we aim to propose a definition for which all arrangements (Fig. 3A-E) are identified as threadings and, because discovered through a persistence homology analysis, we will dub them “homological threadings” (H-threadings).

In order to characterize the topological interaction between two rings, say $i$ and $j$, we construct four PDs from the configuration data of the two rings. Two of them are the PDs of each ring, which we write PD($i$) and PD($j$). The third one is the PD obtained from the union of the configurations of these two rings PD($i\cup j$) and the last one is the union of PD($i$) and PD($j$) which we denote as uPD($i,j$) $\equiv$ PD($i$) $\cup$ PD($j$). From these, we create the following PDs as set difference

where ${\rm PD}(j\rightarrow i)$ represents PD points which are associated to a loop in ring $i$ that is however lost in the PD in which both configurations are included. Such a loop, identified as a lost point of ring $i$, should be interpreted as a loop being threaded by ring $j$. Indeed, if any segment of ring $j$ is passing through an area that was identified as an opening in ring $i$’s conformation, then it must represent an obstacle for the ring $i$ in the sense described in Fig. 3. By symmetry, the same applies to ${\rm PD}(i\rightarrow j)$ with the interchange between $i$ and $j$. These relations constitute our definition of H-threading events. At the same time, ${\rm nPD}(i,j)$ represents new points created by considering both conformations and that are absent in the union of the PDs obtained using isolated conformations; these points represent new loops that are formed when balls of a certain radius bridge segments belonging to different rings and therefore may also be considered as obstacles, restricting the conformational freedom of the individual rings $i$ and $j$.

Our definition of H-threading respects the directionality of this interaction: points in ${\rm PD}(i\rightarrow j)$ identify arrangements in which ring $i$ is actively H-threading $j$ while points in ${\rm PD}(j\rightarrow i)$ find ones in which $i$ is passively H-threaded by $j$. To quantify these events, we can count the number $n^{(A)}$ of rings through which a reference ring is threading (active threading) and the number $n^{(P)}$ of ring threading into a reference ring (passive threading); their statistics are shown in Fig. 4A-B. Due to symmetry, $\langle n\rangle=\langle n^{(A)}\rangle=\langle n^{(P)}\rangle$, and we find that it grows with ring length as $\langle n\rangle\propto N^{1/6}$ (or $\langle n\rangle=\ln N$ see Fig. 4C and discussion below). To motivate such a scaling law, let us count the number of surrounding rings, which have contacts with the reference ring. This quantity, i.e., the coordination number, is estimated as $\sim R^{3}M/V=\rho R^{3}/N\sim N^{\alpha}$ with $\alpha=3\nu-1$, where $R\sim N^{\nu}$ is the ring spatial size. The use of effective exponents $\nu=2/5$, valid in intermediate cross-over $N$ regime, we find $\alpha=1/6$. On the other hand, the compact scaling exponents $\nu=1/3$ is expected for asymptotically long rings, then, we get $\alpha=0$, implying a logarithmic dependence on $N$. Smrek and Grosberg observed a similar chain length dependence for the number of rings penetrated by a single ring in their minimal surface analysis Smrek and Grosberg (2016). The distinction between $n^{(A)}$ and $n^{(P)}$ shows up in their distributions, with larger width for $n^{(P)}$. This may reflect the fact that passive threadings are more sensitive to the opening of loops, i.e., larger $n^{(P)}$ for rings having a primal loop with longer life span, and vice versa.

We finally stress that since $\langle n\rangle$ is the average number of H-threadings per ring, our findings strongly suggest that our systems of rings are percolating, in the sense that all rings are connected to each other through H-threadings.

So far, we have analyzed the H-threading statistics per ring as a function of ring length. However, it is natural to explore the statistics of H-threadings as a function of the loop size into which the H-threading takes place. To this end, we define the Betti number

for a given configuration of the $m$-th ring $\{{\vec{r}}_{n}^{(m)}\}$ with monomer index $n\in[1,N]$.

The Betti number $\beta^{(m)}(\alpha;N)$ counts the number of loops in the $m$-th ring if observed at the spatial resolution $\sqrt{\alpha}$. In Fig. 4D, we compare the average Betti number per ring

with the quantity

where $({\tilde{b}}_{k},{\tilde{d}}_{k})$ are the birth and death scales of the loops associated with all H-threading events, i.e.,

Therefore, while Eq. 1 represents the total number of loops per ring (averaged over rings) if observed at the spatial resolution $\sqrt{\alpha}$, Eq. 2 counts only loops which are associated with a H-threading event.

Mathematically, the set the element of which is $(\tilde{b}_{k},\tilde{d}_{k})$ is no longer a PD in any sense because the right-hand side of Eq. 3 is defined not as a disjoint union but a union. Only for our analysis, the difference between union and the disjoint union does not affect the original PD analysis because all multiplicity value in $PD(i)$ is confirmed to be unity numerically.

By comparing these two quantities (see Fig. 4D, solid and dashed lines, respectively) one can notice that: (i) The number of threaded loops approach the total number of loops from below, and these two quantities converge at large enough $\alpha$; (ii) the scaling function $f(\alpha)\sim e^{-\sqrt{\alpha/\alpha_{0}}}$ well captures the scale-dependence of the Betti number except within an uninteresting small $\alpha\sim$ monomer scale. There is a small deviation from the exponential behavior at $\sqrt{\alpha}\gtrsim 7$. This may be a signal to a self-similar loop statistics as conjectured in loopy globule model Ge et al. (2016) (see also Appendix) ; (iii) the Betti number follows a super-linear scaling with ring length $\beta(\alpha;N)=N^{\kappa}f(\alpha)$ with exponent $\kappa\simeq 1.2\pm 0.1$, indicating that the number of (threaded) loops grow at least extensively with the rings’ contour length (see Fig.4D inset).

We note that this finding is compelling evidence for the argument that even if the number of neighbouring rings H-threaded at any one time is bounded to plateau in the asymptotic limit (Fig. 4C) the number of H-threaded loops (and hence of the induced TCs) increases at least linearly with the rings’ contour length.