Persistent Homology Reveals the Role of Stiffness in Forming Topological Glasses in Dense Solutions of Ring Polymers

MRC Human Genetics Unit, Institute of Genetics and Cancer, University of Edinburgh, Edinburgh, EH4 2XU, UK

Ring polymers exhibit distinctive properties compared to their linear counterparts ^{1, 2, 3, 4}. Despite extensive research, a thorough understanding of topological constraints in ring polymer melts remains a significant challenge in polymer physics ^{5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29}. One key feature thought to define topological constraints in ring polymers is the interpenetrating structure known as “threading”. Threading occurs when one ring polymer penetrates the loop of another ring polymer, with the penetrating ring classified as active and the penetrated ring as passive, illustrating the asymmetric and hierarchical nature of the threading network. For sufficiently long rings, this threading network can evolve, eventually leading to the formation of “topological glasses,” where the relaxation time is expected to increase drastically with respect to the extent of threading ^{30, 31, 32, 33}.

Analyzing threading and clarifying its relationship with glass-like properties is crucial. While several approaches for quantifying threading have been proposed, including methods based on minimal surface ^{34, 35} and geometric analysis ^{36, 37}, Landuzzi et al. introduced a method for quantifying the threading of ring polymers using persistent homology (PH) ³⁸. PH is a mathematical tool that characterize topological features such as “loops” from point cloud ^{39, 40, 41}. Specifically, Landuzzi et al. investigated threading structures using PH from MD simulations with the Kremer–Grest (KG) model ⁴² for ring polymers. Of particular interest was the chain length $N$ dependence of ring polymers up to $N=2048$ at a monomer number density of 0.1, incorporating a bending potential $U_{\mathrm{bend}}(\theta)=\varepsilon_{\theta}(1+\cos\theta)$, where $\theta$ represents the angle formed by consecutive bonds and $\varepsilon_{\theta}=5$ (see Eq. (3) for details). This bending potential effectively models the polymers as worm-like chains, analogous to the Kratky–Porod model ⁴³.

We recently performed MD simulations using the KG with two types of ring polymers: semi-flexible ($\varepsilon_{\theta}=1.5$) and stiff ($\varepsilon_{\theta}=5$) rings with a fixed chain length $N=400$ to investigate the influence of chain stiffness on their dynamic properties ⁴⁴. The rearrangement dynamics of the center of mass (COM) were analyzed, with a focus on dynamic heterogeneity to clarify glassy behavior. Our results demonstrated that stiff ring polymers exhibit pronounced glassy behavior accompanied by dynamic heterogeneity, whereas semi-flexible ring polymers display homogeneous dynamics characterized by a Gaussian distribution of COM displacement. This distinction suggests that the dynamic properties of ring polymers are fundamentally influenced by the chain stiffness, emphasizing the need to examine threading structures across varying degrees of chain stiffness.

The purpose of this study is to elucidate the influence of the chain stiffness and monomer number density of ring polymers on their threading structures. We first analyze the connectivity of COM using PH. Subsequently, we characterize the active and passive threading structures between pairs of ring chains through PH. Through these analyses, we clarify the topological characteristics of ring polymers and their relationship to glassy behavior, informed by insights gained from the rearrangement dynamics of COM.

All MD simulations were performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) ⁴⁵. Length, energy and time are represented in units of $\sigma$, $\varepsilon_{\theta}$ and $(m/\varepsilon_{\mathrm{LJ}})^{1/2}$, respectively. Additionally, the temperature is expressed in units of $\varepsilon_{\mathrm{LJ}}/k_{\mathrm{B}}$, where $k_{\mathrm{B}}$ is Boltzmann constant. We fixed the temperature $T$, chain length $N$, number of chains $M$ as $T=1.0$ and $N=400$, and $M=100$, respectively. Throughout the simulations, temperature was controlled using the Nosé–Hoover thermostat, with a time step of $\Delta t=0.01$. The monomer number density $\rho\sigma^{3}$ ($=NM\sigma^{3}/V$) was varied as 0.1, 0.2, 0.3, 0.4, and $0.5$ for each degree of chain stiffness. Henceforth, $\rho$ will be referred to as density.

Here, we briefly outline PH: A set of coordinates such as beads of chains or COM, denoted as $\{\bm{r}\}=\{\bm{r}_{1},\bm{r}_{2},\ldots,\bm{r}_{\kappa}\}$, is used as input data, where $\kappa$ is the number of coordinates. At each coordinate, assign a virtual sphere with radius $\sqrt{\alpha}$, where $\alpha$ is a parameter. Initially, when $\alpha=0$, all points are treated as disconnected components. As $\alpha$ increases, the spheres begin to overlap, and connected components form, creating edges and facets. During this process, the topological features varies discontinuously with respect to $\alpha$, i.e., the loops will appear and disappear. We record the radii for appearance and disappearance as $b$ (birth) and $d$ (death) respectively for each hole, and introduce persistence diagram (PD) as a collection $(b,~{}d)$ of all holes. In this context, a zero-dimensional hole represents a connected component, while a one-dimensional hole represents a loop. The instance of the PH procedure is illustrated in Fig. 1, which shows the PD of a ring polymer in solution.

PD captures not only the topological features at a specific radius, i.e., a threshold of connection, but also how these features change as the threshold increases. All analysis were performed using the HomCloud ⁴⁶.

Figure 2 demonstrates that $\lambda(\alpha)$ decreases and converges to zero at a specific length scale $\alpha$. This behavior indicates percolation, where clusters are formed by virtually connected COMs. The characteristic length scale is $\alpha=\left\langle R_{\mathrm{g}}^{2}\right\rangle/2$, where a virtual bond is considered to have formed if the distance $r_{ij}$ between the COMs of ring polymer pair $(i,j)$ satisfies $r_{ij}\leq\left\langle R_{\mathrm{g}}^{2}\right\rangle$ (see horizontal lines in Fig. 2). For flexible ring chains with $\varepsilon_{\theta}=0$, $\lambda$ takes finite values for $\alpha\leq\left\langle R_{\mathrm{g}}^{2}\right\rangle$ across all densities, indicating the presence of numerous small clusters. In contrast, stiff ring chains with $\varepsilon_{\theta}=5$ exhibit $\lambda\approx 0$ at $\alpha=\left\langle R_{\mathrm{g}}^{2}\right\rangle/2$, suggesting the formation of percolated networks among COMs of ring chains. Furthermore, the length scale of $\alpha$ exhibiting a plateau of $\lambda\approx 1$ approximately corresponds to the characteristic core length, analogous to the behavior of ring polymers modeled as soft macromolecules. As $\varepsilon_{\theta}$ increases, this core length is reduced, as shown in Fig. 2. These findings indicate that in flexible ring chains, the cores are large and overlap each other, but their COMs are not connected with one another. In contrast, for stiff ring chains, the smaller cores and relatively larger radius of gyration lead to an increase in the number of virtual bonds between COMs. Thus, the density and chain stiffness strongly influence the structural and dynamic behavior of ring polymer systems. See Supporting Information for further discussion on the dependence on chain stiffness and density dependence of the mean square gyration of radius $\left\langle R_{\mathrm{g}}^{2}\right\rangle$, radial distribution function of COMs, $g(r)$, and the number of virtual bonds $Z_{\mathrm{b}}$ defined by the radial distribution function of COMs, as shown in Fig. S1-S3 in Supporting Information, respectively.

The next analysis focuses on the one-dimensional hole, i.e. “loop” structure, characterized by PH. Specifically, PH is performed on each individual ring polymer $i$ by using the monomer coordinates as input, which generates a persistent diagram (PD) denoted as $\mathrm{PD}(i)$. This analysis reveals the birth and death of topological features such as loops within the structure of the polymer. Furthermore, the “life” of the loop is defined as the vertical distance from the diagonal line in the PD, denoted as $l=d-b$, which quantifies how long during the increase of $\alpha$ the loop persists before disappearing. Thus, larger values of $l$ indicates longer-lived loops, reflecting more stable topological features of the system against the change of threshold. Next, we performed PH for $(i,j)$ pairs of ring polymer chains using the set of their coordinates as input, generating PD denoted as $\mathrm{PD}(i\cup j)$. Since threading occurs when the loop of one ring polymer disappears due to penetration by another polymer, $\mathrm{PD}(j\to i)=\mathrm{PD}(i)\backslash\mathrm{PD}(i\cup j)$ allows us to quantify the loops being threaded ³⁸. Here, the set difference operator $\backslash$ represents the subtraction of topological features that vanish when polymer $j$ interacts with polymer $i$. In this context, polymer $i$ is considered “passive” while polymer $j$ is the “active” participant in the threading process. Thus, this approach quantifies the threading structures between pairs of ring polymers. The probability density distributions of $\mathrm{PD}(i)$, $\mathrm{PD}(i\cup j)$, and $\mathrm{PD}(j\to i)$ with $\varepsilon_{\theta}=1.5$ and $5$ at densities $\rho=0.1$ and $0.5$ are shown in Fig. S4-S7 in Supporting Information.

Figure 3(a) shows the density $\rho$ dependence of $\beta_{1}(\alpha)$ and $\tilde{\beta}_{1}(\alpha)$ at the highest density $\rho=0.5$. The stiff ring exhibits a broader peak at larger length scales $\alpha$ compared to that of the flexible ring, indicating the presence of large loops. Peaks emerge and decay over similar length scales for both $\beta_{1}(\alpha)$ and $\tilde{\beta}_{1}(\alpha)$. The peak intensity, however, depends on the chain stiffness $\varepsilon_{\theta}$. For sufficiently large $\alpha$, the convergence of $\beta_{1}\sim\tilde{\beta}_{1}$ indicates that all large loops are involved in threading. The all results of $\beta_{1}(\alpha)$ and $\tilde{\beta}_{1}(\alpha)$ with varying $\varepsilon_{\theta}$ and $\rho$ are shown in Fig. S8 in Supporting Information.

The full width at half maximum (FWHM) of $\beta_{1}(\alpha)$ and $\tilde{\beta}_{1}(\alpha)$ are plotted in Fig. 3(b) and (c). Note that the FWHM provides insights into average behavior of loop sizes. The FWHM of $\beta_{1}(\alpha)$ and $\tilde{\beta}_{1}(\alpha)$ are larger for stiff rings compared to those of flexible rings, indicating that stiff rings have a broader distribution of loop sizes at low density. As the density increases, the FWHM of $\beta_{1}(\alpha)$ and $\tilde{\beta}_{1}(\alpha)$ converge to a common value. This result suggests that the stiffness of ring chains significantly influences the formation of large loops in dilute solutions, while the rings strongly overlap and srink the size of loops as the density increases. Note that flexible ring chains with $\varepsilon_{\theta}=0$ exhibit non-monotnic behavior at low density $\rho\leq 0.3$. This behavior is attributed to the peak of $\beta_{1}(\alpha)$ at $\alpha\approx 0.6$, as shown in Fig. S8(a), which is resulted from short-lived loops along the diagonal line in the PD.

Figure 4 presents the probability density distribution of $N_{\mathrm{a}}$ and $N_{\mathrm{p}}$, respectively. Note that $N_{\mathrm{a}}$ and $N_{\mathrm{p}}$ were calculated by including threading at all length scales, with the threshold $l_{\mathrm{th}}$ set to zero. It is demonstrated that for both $N_{\mathrm{a}}$ and $N_{\mathrm{p}}$, the peak shifts to higher values with increasing chain stiffness $\varepsilon_{\theta}$ and density $\rho$, indicating a greater occurrence of threading. Notably, the density dependence of the distribution becomes more pronounced for stiff rings compared to that of flexible rings. In addition, $N_{\mathrm{p}}$ exhibits a slightly broader distribution than $N_{\mathrm{a}}$ at high density for stiff rings. This asymmetric property between $N_{\mathrm{a}}$ and $N_{\mathrm{p}}$ was found to be pronounced for longer stiff rings, suggesting that the passive threading is significantly influenced by the presence of long-lived loops. In other words, larger loops are likely to be involved in the passive threading.

We further characterize the long-lived active and passive threading structures by introducing the threshold value $l_{\mathrm{th}}$, which has a dimension $\sigma^{2}$. Since points near the diagonal line are considered noisy, we introduce $l_{\mathrm{th}}$ to filter out threading associated with loops of short life, thereby characterizing loops that are mostly correlated with topological constraints. While the results for varying $l_{\mathrm{th}}$ are not displayed, the threshold value $l_{\mathrm{th}}=9$ was determined to capture the most relevant characteristics, and the corresponding results are shown below.

Figure 5 illustrates the density dependence of probability density distribution of active and passive threading numbers, $N_{\mathrm{a}}$ and $N_{\mathrm{b}}$, at $l_{\mathrm{th}}=9$. For flexible ring chains, both $N_{\mathrm{a}}$ and $N_{\mathrm{p}}$ show the tendency of the decrease toward zero as the density $\rho$ increases. This trend is expected to become more pronounced as the threshold value $l_{\mathrm{th}}$ increases. This observation suggests that the number of loops necessary for threading becomes minimal in higher densities, consistent with the overlapping structures between the crumbled globules characteristic of flexible ring chains. In contrast, for stiff ring chains, the distribution of $N_{\mathrm{a}}$ exhibit a peak at $N_{\mathrm{a}}\approx 20$ across all densities, whereas the distribution of $N_{\mathrm{p}}$ shows two distinct peaks, one at $N_{\mathrm{p}}=0$ and another at $N_{\mathrm{p}}\approx 20$. In addition, the latter peak broadens as the density $\rho$ increases. This observation implies that, when focusing on passive threading of stiff ring chains, they can be categorized into two different types: those having large loops facilitate threading and those lacking such structures. The latter rings are regarded as exhibiting more compact characteristic rather than those of the former.

This work was supported by JSPS KAKENHI Grant-in-Aid Grant Nos. JP24H01719, JP22H04542, JP22K03550, JP23K27313, and JP23H02622 We also acknowledge the Fugaku Supercomputing Project (Nos. JPMXP1020230325 and JPMXP1020230327) and the Data-Driven Material Research Project (No. JPMXP1122714694) from the Ministry of Education, Culture, Sports, Science, and Technology and to Maruho Collaborative Project for Theoretical Pharmaceutics. The numerical calculations were performed at Research Center for Computational Science, Okazaki Research Facilities, National Institutes of Natural Sciences (Project: 24-IMS-C051). DM thanks the Royal Society for support through a University Research Fellowship and the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement No 947918, TAP).