Monomer exchange dynamics in ureido–pyrimidinone supramolecular polymers via molecular simulations

Abstract

The use of synthetic supramolecular polymers, built with monomers that self-assemble via non-covalent, reversible interactions, is rapidly growing in many fields, including energy, environmental, and bioengineering applications. Very recently, ureido–pyrimidinone (UPy)-based supramolecular polymers have been used to synthesize biocompatible hydrogels aiming to mimic the dynamic environment of extracellular matrices. Tuning the dynamics, stiffness, and bioactivity of UPy hydrogels effectively influences cellular behaviour and tissue development. However, a complete understanding of UPy-network dynamics over different length and time scales is still lacking, and even the most advanced experimental approaches are unable to capture the dynamics of monomer exchange with atomistic resolution. Here we present a computational study on UPy supramolecular assemblies in water that uncovers the mechanism of monomer exchange between the UPy supramolecular polymers and their surroundings. Our results, based on atomistic molecular dynamics (MD) simulations combined with enhanced sampling and machine-learning (ML) techniques, show that the fine interplay of solute–solvent interactions is the main engine of monomer motion, which makes UPy supramolecular polymer ends more dynamic as compared to the static backbone. This computational work complements the qualitative experimental evidence on supramolecular dynamics with the mechanism of monomer exchange, revealing the most favorable environment for supramolecular polymer damage as well as the underlying principle of self-healing.

---

Living systems possess an intrinsic ability to mitigate various forms of damage, showing self-healing and regenerative properties.^1–3 Inspired by such dynamic features, scientists have advanced the synthesis of functional materials with on-demand reversibility and stimuli-responsiveness.^4–8 One common strategy relies on the design of supramolecular materials,⁹ i.e., assemblies of non-covalently bound monomers, holding great promise in environmental applications,¹⁰ drug delivery,^11–15 regenerative medicine,^16,17 skin-like stretchable electronics,^18,19 and anticorrosive coatings.^20,21 Examples of monomers forming supramolecular structures include benzene 1,3,5-tricarboxamide (BTA), self-assembled in one-dimensional (1D) supramolecular polymers via core–core stacking and three-fold hydrogen bonding,^22,23 benzotrithiophene building-blocks,^24,25 peptide amphiphiles,^26–28 or metal-coordinated porphyrins.^29–31 In this realm, UPy molecules have also been largely used by taking advantage of their ability to self-organize in “hierarchical” fibrillar structures in water. Indeed, UPy monomers dimerize by self-complementary quadruple hydrogen bonding in a donor–donor–acceptor–acceptor (DDAA) fashion (Fig. 1a). These planar dimers can stack on top of each other to form supramolecular polymers, (hereafter also referred to as stacks) which, in turn, can interact and bundle to form more complex fibrillar structures.^32–34

Fig. 1 Molecular dynamics (MD) models of UPy supramolecular polymers. (a) Left: chemical structure of two UPy monomers forming a dimer by self-complementary quadruple hydrogen bonds (HB–dim). R1 indicates a possible functionalization. Right: self-assembly cartoon of UPy dimer building blocks in the axial direction, creating a supramolecular polymer. (b) All-atom MD (AA-MD) simulation of (R1-free) UPy dimer self-assembly. Two MD snapshots at 10 and 300 ns are shown. (c) Chemical structure and AA-MD model of the UPy-C6-u monomer considered in this study. Oxygen, carbon, and hydrogen atoms of the UPy-C6-u core are colored in red, gray, and white, respectively. The side chain is in light gray. (d) UPy-C6-u dimerization free-energy-surface in aqueous solution obtained with well-tempered MetaDynamics (WT-MetaD) simulation. The number of HB–dim and core–core distance are selected as the WT-MetaD collective variables. A, B, and C identify three energy minima corresponding to the stacking, dimerization, and lateral assembly configurations, respectively. (e) AA-MD snapshot of the pre-stacked UPy-C6-u supramolecular polymer made of 20 dimers (20D) considered in this study.

Over the past twenty years, UPy monomers have been functionalized with poly(ethylene glycol) (PEG),^35,36 poly(N-isopropylacrylamide) (PNIPAM),^37,38 glycine (Gly) amino acids,^16,39 among others, in order to improve their biocompatibility and tunability.⁴⁰ The resulting supramolecular polymers have demonstrated intriguing mechanical and dynamic properties in a wide range of multi-component environments, playing a crucial role in the syntheses of biomimetic hydrogels.^41,42 For example, UPy supramolecular polymers have been exploited to mimic the liquid–liquid-phase-separation found in biological fibrils, collagen, and, in general, to reproduce the adaptive behavior of extracellular matrices with tunable stiffness, dynamics, and bioactivity.^35,39,41 However, for inducing cell adhesion in hydrogel networks, tight control over multiscale dynamic processes is required.⁴³ Particularly, a tailored understanding is crucial for both molecular-level dynamics—such as monomer exchange within a single stack—and the bulk dynamics, i.e., supramolecular polymer rearrangements in hydrogels covering larger scales.

Although recent experimental studies have provided estimates for the monomer exchange rate in UPy supramolecular polymers (10% in 1 hour),^35,43 most experimental approaches cannot resolve the molecular and submolecular mechanisms occurring in UPy self-assemblies, where monomeric and oligomeric units continuously exchange, setting a supramolecular “equilibrium dynamics.” This monomer exchange dynamics, on the other hand, can be well detected by multiscale molecular modeling and advanced computational methods.^44–46 In this framework, while BTA supramolecular polymers have been largely investigated with atomistic,^47–49 higher-scale simulations,^44,50,51 and advanced ML tools,^52–54 UPy supramolecular polymers have received much less attention from computational research. Chen et al.⁵⁵ used an umbrella sampling technique to estimate the potential of mean forces between two interacting UPy molecules. While their work successfully highlights the role of the hydrophobic spacer in the UPy dimerization process, it does not capture the range of configurations that UPy building blocks may adopt while self-assembling. Later studies utilized atomistic and coarse-grained simulations to investigate the self-assembly mechanisms occurring in longer UPy supramolecular polymers.^56,57 These approaches, however, are constrained by time and space limitations typical of classical MD, preventing a thorough exploration of the system's configurational space. Additionally, none of the computational studies on UPy self-assemblies conducted so far have investigated the essential phenomenon of monomers' exchange within and outside their self-stacking structure, underlying bioinspired properties such as self-healing and reconfiguration.

Here we present a computational work where MD, enhanced sampling approaches, and ML techniques are integrated to provide an overview of the structural and dynamics features of UPy supramolecular polymers in water. In this study, the UPy core is functionalized with short carbon spacers terminating with a urea moiety (UPy-C6-u in Fig. 1c). The UPy core forms dimers through quadruple hydrogen bonding, and the hierarchical growth is also driven by bifurcated hydrogen bonds of the urea group flanking the alkyl spacer.^58,59 Starting from pre-assembled UPy dimers—based on literature data and preliminary MD indications—we first investigate the structural properties of these UPy dimer stacks, assessing their size-dependent stability. The dynamics of monomer exchange is then studied via infrequent well-tempered metadynamics (WT-MetaD) simulations,⁶⁰ which accelerate the rupture of the dimerization hydrogen bonds across the supramolecular structure. Afterwards, we employ data-driven analyses^61,62 to detect key dynamic environments, which reveal the mechanisms of monomers' exchange within the longest UPy stack. Finally, we also provide insights into the monomer exchange dynamics occurring between bundled UPy supramolecular polymers, the typical higher-level structure detected in experiments.^35,43,58,63 This study provides solid guidelines to describe the structure and dynamics of UPy supramolecular assemblies in water, allowing us to gather useful indications on how the molecular structure of the system relates to the supramolecular dynamics and, subsequently, to the bioinspired properties of the material.^35,43,59

1. Results

1.1. Structural analysis of UPy-C6-u supramolecular polymers

By using all-atom MD (AA-MD) simulations, we here describe the key physical and chemical mechanisms determining the structural stability and dynamic features of UPy supramolecular polymers. As outlined in the introduction, UPy cores dimerize in water, forming four complementary hydrogen bonds, and grow orthogonally via π–π stacking, building 1D supramolecular polymers (Fig. 1a). Such 1D dimer stacks then interact with each other, aggregating in fibers, which in turn form the bundle network that constitutes the supramolecular material. This hierarchical self-assembly process is shown here by using AA-MD simulations of R1-free (i.e., without side-chain R1) UPy monomers. One single stack having such a stacked-dimer hierarchical structure is spontaneously formed after 300 ns-long MD simulation, starting from a suspension of 42 dispersed UPy monomers (Fig. 1b). This outcome confirms that the combination of hydrogen bonds and π–π interactions governs the mono-directional stacking of UPy cores. However, due to this high aggregation propensity, it is essential to enhance water solubility and prevent solution precipitation. Adding hydrophobic and hydrophilic side chains, i.e., functional R1 groups, to the UPy core is a strategy to promote self-assembly in water and tune the resulting properties.

Building on experimental studies,^58,59 here we select a relatively simple molecular design for the UPy motif, namely including one hydrophobic spacer with 6 carbon atoms (C6-spacer) and one urea terminal (the chemical structure of the UPy-C6-u monomer and the relative all-atom model are in Fig. 1c). This monomeric structure adds the complexity of an interacting side-chain without overly increasing the simulation times required to capture the monomer exchange dynamics across the supramolecular structure. Unlike the self-aggregation of R1-free UPy monomers (Fig. 1a), it is challenging to observe spontaneous and ordered assembly of UPy-C6-u within the timescales accessible through AA-MD. Therefore, our computational protocol includes the setup of pre-stacked UPy-C6-u supramolecular polymers based on the configurational guidelines emerging from the analysis of UPy dimerization free energy surface (FES). To estimate the FES we adopted well-tempered MetaDynamics,⁶⁴ (WT-MetaD), choosing the number of dimerization hydrogen bonds, HB–dim, and the core–core distance as collective variables (see the Methods section). The resulting FES (Fig. 1d) exhibits three minima, A, B, and C, corresponding to the most probable configurations retained by two UPy-C6-u monomers in aqueous solution. The first minimum, A, represents the two-monomer-stacking arrangement in the orthogonal direction, having a core–core distance of ≈ 0.4 nm in conjunction with the formation of ≈ 2 HB–dim. The second minimum, B, corresponds to the dimer assembly driven by self-complementary quadruple hydrogen bonding in a DDAA pattern, i.e., the number of HB–dim = 4, associated with a core–core distance of ≈ 0.6 nm, matching with earlier experimental results.⁶³ Finally, the third minimum, C, lays down in the phase space of zero HB–dim and core–core distances around ≈ 1.2 nm. The MD representative snapshot of C state (Fig. 1d) reveals that the UPy-C6-u monomers are bound head-to-tail, held together by urea–core interactions. This suggests a possible arrangement of the inter-supramolecular polymer aggregation within a larger self-assembly. Based on these favorable dimer configurations and on the experimental evidence of supramolecular polymer structures,^32–34 we arranged UPy-C6-u monomers in stacks formed by pre-assembled dimers, and we tested the stability of such structures via AA-MD (Fig. 1e).

The structural stability and dynamics of the pre-assembled UPy-C6-u supramolecular polymers are investigated by simulating, for 400 ns in aqueous solution, 4 stacks of distinct sizes, i.e. made of 2, 5, 10 and 20 dimers (2D, 5D, 10D and 20D) (Fig. 2a). The resulting MD trajectories are firstly analyzed by computing the radial distribution function (RDF) among all monomers, considering the distance between the geometric centers of UPy cores. Regardless of stacks' length, the RDF profiles show multiple regular peaks obtained at increased core–core distance (Fig. 2b). For each system, the first peak occurs at 0.37 nm, which approximately corresponds to the position of the free energy minimum A (Fig. 1d). This peak therefore identifies the first neighbors along the stacking direction. The next peak, located at ≈ 0.6 nm, corresponds to the monomer–core distance relative to the dimer formation, as confirmed by the free energy landscape in Fig. 1d (B minimum). The subsequent lower and broader peaks include both the higher-order stacking and dimerization configurations (that are degenerate at this order). We observe that longer stacks exhibit higher peaks. Considering that the reported RDF is normalized by the number of monomers, we expect that the purely stacking peaks increase with supramolecular polymer size due to the higher statistics of neighbors. However, the same does not hold for the dimerization peak, and its size dependence is evidence of higher dimer stability in longer supramolecular polymers. We also found that the probability of a UPy-C6-u dimer to be in a perfect planar orientation is higher in a longer stack (20D) as compared to a short one (2D) (Fig. S1 in SI). These elements indicate that an ordered, planar dimer conformation is more likely in longer supramolecular polymers, suggesting cooperativity in the polymerization of UPy-C6-u.

Fig. 2 Structural analyses of UPy-C6-u stacks. (a) AA-MD snapshots of four UPy-C6-u pre-stacked supramolecular polymers of variable size: 20, 10, 5, and 2 dimers (D). (b) Radial distribution function (RDF) computed on each UPy-C6-u monomer forming the four stacks in (a). Each RDF is normalized by the number of stacked monomers. Blue, yellow, orange, and green RDFs identify the 20D, 10D, 5D, and 2D supramolecular polymers. (c) AA-MD snapshot of the 20D pre-assembled stack. The atoms chosen to account for the number of HB–dim are displayed as colored spheres in the zoomed-in view. (d) Number of initial HB–dim (N_HB–dim0) time series, computed on each dimer forming the supramolecular polymers in (a) across the 400 ns equilibrated MD simulations in water. The kernel density estimation (KDE) profile of all HB–dim0 data is plotted on the right-hand side. (e) Number of HB–dim0 averaged over all the dimers in the stack backbones (tips excluded) reported in (a), as a function of t along the 400 ns-long MD trajectories. The 〈N_HB–dim0〉 time series are colored according to the stack size; the time-averaged values of 〈N_HB–dim0〉 are displayed in the boxplot on the right-hand side.

1.2. Monomer exchange dynamics in UPy-C6-u supramolecular polymers

The dynamics of supramolecular polymers is directly linked to the stability of their reversible bonds, which regulate the formation and resolution of structural defects—fundamental to self-healing mechanisms. In other words, strong non-covalent bonds among the UPy-C6-u monomers tend to prevent their reshuffling and thus the overall dynamics. In UPy-C6-u, the supramolecular polymeric structure rests both on the HB–dim forming the dimers and on the π–π stacking that fuels polymeric elongation. Therefore, to obtain quantitative insights into the dynamic nature of our UPy-C6-u pre-assembled stacks, we first focused on the reversibility of the number of HB–dim as the source of structural defects. We calculated how many of the 4 initial HB–dim (HB–dim0) are preserved by each monomer along the MD trajectories of the supramolecular polymer chains (Fig. 2a and Fig. 2c). As a result, the number of HB–dim0, i.e. N_HB–dim0, of each monomer is equal to 4 at the beginning, and it decreases if the original HB–dim conformation is disrupted, as shown by the N_HB–dim0 time series (Fig. 2d). As discussed in the Methods section, the estimated N_HB–dim0 intrinsically includes spurious interactions that generate the fluctuations seen in the profiles (Fig. 2d). We employed relative kernel density estimation (KDE) to highlight the main features of the N_HB–dim0 distribution over the MD data. In particular, the KDE evidences the presence of two main peaks: the largest peak contains those dimers for which the N_HB–dim0 fluctuates around 4, corresponding to the full starting number of HB–dim (the self-complementary dimer assembly shape is shown in Fig. 2c); the second peak indicates that a sample of dimers breaks, forming structural defects. Based on previous results on BTA supramolecular polymer dynamics,⁴⁴ we then differentiate how many N_HB–dim0 are preserved between tip and backbone dimers to identify which ones contribute to the formation of defects. Thanks to this classification, we notice that the number of HB–dim0 in the backbone is overall very stable, oscillating within the main peak of the KDE distribution (N_HB–dim0 ∼ 4 in Fig. 2e), thereby demonstrating that dimer rupture events do not occur in the backbone. In contrast, the N_HB–dim0 computed for specific tip dimers fluctuates from 4 to 0, revealing that full dimer breakage takes place here. In the cases of pre-assembled stacks formed by 20, 10, and 5 dimers, the KDE distribution of the tips (Fig. 3a) is also featured by two main peaks, indicating that intermediate states with N_HB–dim0 from 1 to 3 are more ephemeral. In contrast, the 2D stack, where all four monomers are formally part of the tips, displays a constantly fluctuating N_HB–dim0, alternating defect generation and self-healing, with intermediate dimerization states. Therefore, the 2D stack, where distinction between the backbone and tips is absent, shows a substantially higher monomer exchange.

Fig. 3 Dynamics of UPy-C6-u stacks. (a) Number of HB–dim0 time series relative to the dimers at the tips of differently sized supramolecular polymers (the considered tips are colored in the MD snapshots on the left). The KDE highlights the distribution features for the considered tip dimers, showing the breakage of the initial HBs. (b) Correlation between the time-averaged N_HB–dim0 computed on each UPy-C6-u monomer forming the supramolecular polymers and its associated sasa. The data points are colored based on stack size, while a nonlinear fitted curve is shown with a black line. (c) Cumulative distribution functions (CDFs) of the defect creation times (i.e., dimer rupture) resulting from 30 repetitions of infrequent WT-MetaD. Two distributions are shown, either involving a tip dimer (blue fitting curve) or a backbone dimer (red fitting curve). τ is the characteristic time-scale extracted via CDF fitting. (d) AA-MD snapshots of the 10D stack in 1 μs-long MD simulation. Initial tip dimers are highlighted in dark gray and yellow, while the supramolecular polymer backbone is in light gray. The MD snapshots show the defect formation, monomer sliding, and self-healing process involving the system.

The unique behavior of HB–dim depending on whether the dimers are arranged as tip or backbone suggests a crucial role of the physical environment in which the defects occur, specifically the key contribution from the competitive solute–solvent interaction. This is validated by a clear correlation between the time-averaged N_HB–dim0 and the solvent-accessible surface area (sasa) associated with each monomer: the monomers that exhibit a larger contact with water (higher sasa) are more prone to dimer disassembly (lower 〈HB–dim〉) and hence to defect formation (Fig. 3b). The results achieved at this stage definitely confirm that the physical-chemistry source behind defect generation in the UPy-C6-u supramolecular polymers is the hydration; nevertheless, the limited time and space scales considered in these standard MD simulations make it difficult to generalize these findings. For this reason, we used the infrequent WT-MetaD technique to stimulate the rupture of a dimer bond in a 20D stack. This allows estimating the characteristic time associated with defect generation from either backbone or tip dimers. This analysis shows that the formation of a defect from a tip dimer (Fig. 3c, blue curve) is 7 orders of magnitude faster than the formation of a defect from a backbone dimer (Fig. 3c, red curve). Such a difference suggests that tip dynamics also dominates real systems where the tip-to-backbone ratio is lower than in our model supramolecular polymers.

Although the investigation around the number of HB–dim0 allows the detection of initial dimer disassembly and defect formation, in general, such an analysis is not suitable to capture the self-healing events, as N_HB–dim0 cannot detect the dimerization of monomers not coupled in the starting supramolecular polymer (see also discussion in the Methods section). This is, for example, the case reported in Fig. 3d where we compare some configurations of the 10D stack along 1 μs of MD. After 100 ns, one of the tip dimers (in dark gray in Fig. 3d) disassembles, and one of the two unbound monomers slides along the supramolecular polymer, stacking on the opposite end. On approaching 1 μs, this monomer self-assembles with one of the monomers at the tip (in yellow), forming 4 complementary hydrogen bonds (number of HB–dim = 4). This is an example of defect self-healing, which characterizes the properties of this system. A similar dynamics also occurs in the 20D stack (Fig. S2 in SI). To systematically study the mechanisms of defect formation and self-healing, we analyze the dynamics of the system from a different viewpoint. Recently developed descriptors, coupled with ML tools, have been shown to accurately capture diverse structural environments within a self-assembly, including the probability of building blocks to transfer among the detected domains.^52,65 Dynamic environments, on the other hand, have been directly extracted by using descriptors pointing out the time evolution of neighborhood environments.^53,61,62,66 In the latter context, time smooth overlap of atomic position (τSOAP)⁶¹ focuses on the variations in the supramolecular structure of the system along the MD trajectory. More in detail τSOAP provides a scalar quantity (normalized from 0 to 1) that identifies the rate of structural rearrangements that occurred in the surroundings of selected centers. These centers are identified with specific groups, defining the main interactions between building blocks, thus providing a useful classification of the monomer exchange dynamics in supramolecular structures.⁶¹ In our case, we locate the centers of τSOAP calculation at the center of mass of the 4 atoms forming the HB–dim so that, as shown in our previous results, the arrangement of different centers is informative not only about the formation of dimers but also about the stacking. With this definition, τSOAP captures the dynamic arrangement of the UPy-C6-u interacting groups (see the Methods section for further details). Analyzing the statistics resulting from the time evolution of τSOAP for each monomer, we can classify the main dynamic domains featured in the supramolecular polymer.

We thus simulate the 20D stack system for 1.6 μs, computing the τSOAP descriptor for each of the 40 monomers along the MD trajectory. The time evolution of τSOAP, together with its KDE, shows a data region particularly dense around τSOAP = 0.5, identifying the most probable dynamic state, while sparse fluctuations arise both below and above this dense data distribution (Fig. 4a). A proper cluster analysis of these data series detects the presence of three domains, one corresponding to the more dense data region and the other two above and below, respectively. Such cluster analysis also allows us to associate the identified dynamic domains to monomers that are found in a different structural state: (I) the stacked bound monomers, which form the most populated cluster characterized by an intermediate value of τSOAP, i.e., a moderate rate of structural rearrangement of the environment (light gray cluster in Fig. 4a); (II) the stacked unbound monomers, which manifest the highest dynamics in terms of neighborhood reconfiguration (teal area in Fig. 4a); and (III) the monomers dissolved in solution, or sliding along the supramolecular polymer surface, which correspond to the lowest values of τSOAP (cyan area in Fig. 4a). In configuration III, the interacting centers are far from the others, thus, their environment appears static from our definition of τSOAP (cyan region in Fig. 4a). Three MD snapshots of the 20D stack at distinct time steps are reported with the monomers' colors corresponding to the classification in Fig. 4a. In this graphical representation, the ideal assembly configuration (I), characterizing the bound monomers, as well as the defect events, identified as either (II) or (III) states can be clearly visualized. Based on this classification, we can now interpret the τSOAP signal relative to each monomer (center): the defect formation events take place when a monomer transfer occurs from I to II or from I to III—and self-healing events arise when a monomer moves from II or III to domain I (Fig. 4c). Fig. 5 shows some emblematic examples of monomer exchange dynamics. For instance, M1 and M2 are two tip monomers initially bound, as distinctly confirmed by their τSOAP values, averaging around 0.5 for the first 250 ns (green and purple τSOAP signals in Fig. 5b and relative MD snapshots in Fig. 5a). Then, after 250 ns, a defect forms and while M1 remains bound to the tip (with its τSOAP value transitioning to higher values, II), M2 starts sliding along the supramolecular polymer (with τSOAP value shifting towards 0, III). After 1 μs, we observe a self-healing event, in which M2, which in the meanwhile has reconfigured as a stacked unbound monomer (domain II), dimerizes with M40, which travels toward M2 from the other end of the supramolecular polymer (the orange profile in Fig. 5b and final MD snapshot in Fig. 4a).

Fig. 4 Monomer exchange pathway. (a) τ SOAP time series computed on the interacting centers of each monomer of the 20D UPy-C6-u supramolecular polymer in Fig. 2a, along a 1.6 μs MD. The KDE distribution of the data is shown on the right. The cluster analysis carried out on the τSOAP data detects three domains: (I) stacked bound monomers (light gray), (II) stacked, unbound monomers (teal), and (III) travelling monomers (cyan). (b) Three emblematic MD snapshots of the 20D stack, where monomers are colored according to the domains identified. (c) Scheme of possible monomer exchange pathways: starting from an ideal configuration of stacked dimers into a supramolecular polymer (I), a defect formation induces either stacked unbound (II) or travelling monomers (III), which eventually may exchange their configurations or self-heal by restoring a new stacked dimer (I).

Fig. 5 Monomer exchange pathway. (a) MD snapshots of the 20D stack. The backbone is colored in light gray, while the initial tip monomers are evidenced by different colors: M1 in green, M2 in purple, M39 in blue, and M40 in orange. The disassembly of the M1–M2 dimer is shown on the top side, showing defect formation. At 248 ns, the initial dimerization HBs break apart; then M2 starts sliding on the supramolecular polymer, while M1 remains unbound. The self-healing process instead is clarified following the M40 pathway on the downside: the decisive disassembly of the M39–M40 dimer occurs at around 1 μs, then M39 starts travelling, while M40 remains stacked and unbound. At 1455 ns M40 detaches from the supramolecular polymer and reaches M2 monomer, forming a new dimer. (b) τSOAP time series for M1, M2, M39 and M40 as defined in (a). Here, several pathways of defect creation and self-healing can be identified, clarifying the mechanisms of monomer exchange dynamics.

Overall, this second analysis also confirms that monomer exchange dynamics concentrates at the tips of the supramolecular polymers, highlighting that both defect creation and self-healing preferentially occur at the ends (Fig. S2 in the SI). Coherent results are also obtained by combining τSOAP and LENS descriptors (Fig. S3 in the SI).^53,62 This observation is consistent with the relatively slow dynamics of UPy-C6-u supramolecular polymers detected in the experiment.

We finally explored how the monomer exchange dynamics is affected when multiple UPy-C6-u stacks (Fig. 6) interact to form fibers, typical of the UPy supramolecular structure.^35,43,58,63 We therefore carried out an unbiased MD simulation of three 20D pre-assembled stacks initially placed next to each other, following their interaction along the trajectory. As shown by the sequential snapshots reported in Fig. 6a, few events of inter-supramolecular polymer exchange were captured within a 3 μs timeframe.

Fig. 6 Inter-supramolecular polymer monomer exchange dynamics. (a) Three snapshots of AA-MD simulating three parallel 20D stack placed next to each other. Monomers initially located at the tips of each supramolecular polymer, as well as the first neighbors, are colored in green, cyan and blue, depending on the supramolecular polymer they belong to. The other, backbone monomers are colored in gray. The inter-supramolecular polymer dynamics is localized at the tips. (b) Aggregation propensity (AP, see text for the definition) of the three stack along the AA-MD simulation. (c) Average number of inter-supramolecular polymer contacts per monomer (nc) in case of tips/subtips (cyan curve) and backbones (gray curve).

Entering more in depth into the inter-supramolecular polymer dynamics we observe that the aggregation propensity (AP), namely the ratio between the sasa of three ideally isolated stacks over the sasa of the assembly computed at time t, signals the tendency of the supramolecular polymers to aggregate in fibers, progressively reducing their solvent exposure (Fig. 6b). This aggregation is mainly driven by side-chain interactions, as shown by the relatively small number of inter-supramolecular polymer contacts formed by the core UPy motifs, i.e. excluding side-chain contacts (Fig. 6c). Interestingly, monomer exchange within the fiber occurred almost exclusively at the supramolecular polymer tips. Separating the average number of inter-stack contacts “nc” established by each tip monomer (end dimers and their first neighbors) from those established by backbone monomers revealed a striking difference between the two. The value of tips’ nc mostly fluctuates within the 40–60 range, whereas the backbone's nc mostly remains around 5 (Fig. 6c). Besides this quantitative evidence, MD snapshots visually confirm that exchanges are localized at the tips (Fig. 6a). Moreover, by performing long MD simulations (7 μs) of three infinite supramolecular polymers at T = 343 Kelvin, we could observe that backbone exchanges are not strictly forbidden, but they are far less frequent than tip exchanges (Fig. S4). This evidence supports our interpretation of why UPy supramolecular polymers display slower dynamics compared to other supramolecular polymers such as BTA.^35,41 Our analysis also indicates that the monomer sliding observed along the surface of an isolated chain becomes less favorable in fibers, where side-chain interactions between adjacent UPy stacks hinder such motion. These simulations on multiple supramolecular polymers therefore provide valuable context for interpreting single-chain monomer exchange results in the framework of hierarchical structures formed by UPy motifs.

2. Conclusions

We here report for the first time a computational MD study on the monomer exchange dynamics of UPy-C6-u supramolecular polymers in aqueous solution. In particular, we shed light on the processes of defect creation and self-healing pathways, which are key aspects for monomer motion in supramolecular polymers.

First, the dimerization free energy surface between two UPy-C6-u monomers in water has been investigated to assess the supramolecular architecture of the system, based on quadruple hydrogen bonds, i.e. HB–dim, creating dimer units that stack into supramolecular polymers. Structural analysis carried out on such dimers' stacks of distinct sizes has highlighted a cooperative effect, by which the quadruple HB–dim appears to be more stable in longer supramolecular polymers, with subsequent increased stability of the entire structure.

The monomer exchange mechanism in UPy supramolecular polymers was then explored more in detail. By following the initial number of HB–dim0 along MD trajectories, we found that the dimers' disassembly seen as a consistent variation of the number of HB–dim0 involves exclusively tip monomers rather than backbone. To explain this evidence, the solvent accessible surface area was calculated for each monomer along the MD trajectory. The resulting data demonstrate a correlation between hydrogen bond breaking and monomer hydration, thereby suggesting that competitive solute–solvent interactions are the main driving force for dimer rupture, i.e. defect formation. This result is also quantitatively supported via infrequent WT-MetaD simulations, showing that dimer defect formation occurs several order of magnitude faster at the tips than in the backbone.

We then employed τSOAP,⁶¹ a recently developed descriptor of atomic environment dynamics, to obtain further insights into the most probable defect and self-healing mechanisms taking place along these supramolecular polymers. Specifically, τSOAP coupled with ML tools allow classifying the different monomers according to the dynamics of their supramolecular surroundings, thereby unveiling possible pathways of monomer exchange events. We finally presented AA-MD simulations of UPy supramolecular fibers, showing indications on how monomer exchange dynamics takes place when multiple UPy stacks aggregate. Overall, our results show that the origin of UPy supramolecular polymers dynamics always relies on the mechanism of dimers' disassembly, which, while leading to either traveling or stacked monomers (defect creation), offers a suitable local environment for self-healing (defect resolution). In conclusion, beyond the specific case study, our combined computational approaches establish a modeling strategy capable of systematically investigating hierarchical self-assembly in supramolecular systems. Although experimental techniques have made it possible to quantify the supramolecular dynamics (e.g. estimating a 10% of monomer exchange per hour in UPy supramolecular polymers versus 30–40% per hour in BTA supramolecular polymers^35,41), in this study, we complement the experiments by unveiling the most favorable local environment for dimer breakage, the mechanism of monomer exchange, and the underlying principle of self-healing.

3. Methods

3.1. All-atom molecular dynamics (AA-MD) simulations

The atomistic models of UPy monomers, including both UPy without functionalization and UPy-C6-u, were built with Avogadro⁶⁷ following the chemical structure of the molecules. Gaussian⁶⁸ tool, based on the HF/6-31G*, was used to estimate the generated electrostatic potential, and then the RESP⁶⁹ method was applied to obtain the partial charge distribution within the molecule. The complete parameterization was based on the General AMBER Force Field (GAFF),⁷⁰ using Antechamber.⁷¹

3.1.1. Self-assembly of UPy monomers. The self-assembly MD simulation of 42 non-fuctionalized UPy monomers (i.e., without side-chain R1) was carried out in GROMACS 2021.⁷² First, the parameterized UPy monomers were randomly dispersed in a 10 × 10 × 10 nm³ box filled with water molecules described using the TIP3P model⁷³ and periodic boundary conditions were applied in all box directions. The non-bonded interactions among monomers, including van der Waals and short-range electrostatic interactions, were evaluated within a cut-off radius of 1.4 nm, while for the remaining long-range interactions, a particle-mesh Ewald summation was applied to resolve electrostatics in the Fourier space. Two equilibration steps were performed to reach the thermodynamic conditions of 298 K and 1 bar. The self-assembly simulation, lasting 1 μs, was performed by using the v-rescale thermostat⁷⁴ (τ_T = 0.1 ps) coupled with the c-rescale barostat⁷⁵ (τ_p = 0.1 ps).

3.1.2. Pre-assembled UPy-C6-u supramolecular polymers. We promoted the dimerization of two UPy-C6-u monomers and their self-assembly in the axial direction, forming UPy-C6-u stacks of distinct size from 2 to 20 dimers. We studied the stability of a single UPy-C6-u supramolecular polymer in aqueous solution via classical AA-MD simulations carried out with the open-source software GROMACS 2021.⁷² Each single pre-assembled stack was first solvated in a 5 × 5 × 5 nm³ box filled with water molecules described using the TIP3P model⁷³ and periodic boundary conditions were applied in all box directions. Note that for the 20D stack, a 10 × 10 × 10 nm³ box was considered. The non-bonded interactions among monomers, including van der Waals and short-range electrostatic interactions, were evaluated within a cut-off radius of 1.4 nm, while for the remaining long-range interactions, a particle-mesh Ewald summation was applied to resolve electrostatics in the Fourier space. Our MD protocol consisted of a first step of energy minimization and two consequent equilibration steps. Initially, to reach an equilibrium temperature of 298 K, we applied the canonical ensemble (NVT) for 2 ns using a Maxwell Boltzmann speed distribution and the v-rescale thermostat⁷⁴ with τ = 0.1 ps. Subsequently, we set the isothermal–isobaric (NPT) ensemble for 2 ns at an equilibrium pressure of 1 bar and an equilibrium temperature of 298 K. In this step, we used the previous thermostat coupled with the c-rescale barostat⁷⁵ with a time constant of 2 ps. During the equilibration steps, the UPy-C6-u atoms were restrained in their initial positions using a harmonic potential with a force constant of 1000 kJ mol⁻¹ nm⁻². Once the desired thermodynamic conditions were reached, the restraint was removed, and a 400 ns-long MD run (integration step dt = 0.002 ps) was carried out by maintaining the temperature at 298 K with a Noose–Hoover thermostat (τ_T = 0.8 ps) and a pressure at 1 bar (τ_p = 2 ps) by imposing the Parrinello–Rahman barostat⁷⁶ Along the MD simulation, the LINCS algorithm was employed to restrain the covalent bonds involving hydrogen atoms.

To compare the stability of the four supramolecular polymers, we first analyzed the production run trajectories by computing the radial distribution function between UPy-C6-u monomers. The HB–dim were estimated in PLUMED 2.6.^77,78 taking into account the coordination (R₀ = 0.12 D₀ = 0.27) among the colored atoms in the zoomed-in view of Fig. 2c. Note that, while performing this estimation, we kept the initial dimer configuration through the complete trajectory analysis. The τSOAP descriptor was instead applied to each monomer of the 20D stack, and specifically on the center of mass of the oxygen and nitrogen atoms involved in possible HB–dim (zoom in Fig. 2c). Thus, for each individual center i, τSOAP_i (t) monitors the i-th local environment changes in terms of neighbor monomers' arrangement along the trajectory, ranging from 0 to 1 for static to highly dynamic neighborhoods, respectively. The instantaneous τSOAP value is defined as:

where p^t_i is the full SOAP feature vector associated with the i-th individual center within a certain cutoff neighborhood (r_cut) at the time step t, as described in detail in ref. 61. Here, r_cut = 0.6 nm was employed. In brief, τSOAP_i (t) tracks the variations of the i-th SOAP vector over time, that is, to what extent the molecular environment related to each center changes at every consecutive time interval Δt in terms of SOAP power spectrum. The unsupervised clustering algorithm of Gaussian Mixture Models⁷⁹ was finally adopted to rationalize the data and to identify the dominant molecular environments in the supramolecular polymer. Similar outcomes are also achieved by applying the recent LEAP analysis,⁶² combining LENS⁵³ and τSOAP descriptors (Fig. S3 in the SI).

3.2. UPy-C6-u dimerization free energy surface

To explore dimerization FES, which shows the thermodynamic phase space of two UPy-C6-u monomers interacting in aqueous solution, we performed extensive 355 ns-long well-tempered MetaDynamics (WT-MetaD) simulations.⁶⁴ We selected as collective variables (CVs) (i) the number of HB–dim (R₀ = 0.12 D₀ = 0.27), and (ii) the core–core distance (see Fig. 1d). The latter CV represents the distance between the centers of the two UPy-C6-u cores. Note that the coordination number was computed as implemented in PLUMED 2.6.^77,78 We chose 10 as a bias factor with an initial Gaussian height of 1.5 kJ mol⁻¹, and a width of 0.5 nm for both the distance and the coordination number, respectively. The Gaussian deposition rate was set to 5000 MD step⁻¹, i.e., every 10 ps. After reaching convergence, we reweighted the FES using the Tiwary–Parrinello estimator⁸⁰ on the same CVs. The WT-MetaD simulations were performed using GROMACS 2021⁷² and PLUMED 2.6.^77,78

3.3. Infrequent WT-MetaD simulations

The formation of defects along the supramolecular polymer is a rare event in the timescales effectively accessible using atomistic models. As demonstrated, the real (unbiased) dynamics of an event is related to the transition time associated with events activated by infrequent WT-MetaD simulations (biased dynamics).^60,80,81 This approach is particularly convenient as it allows one to directly extract information on the kinetics of the activated transition from the biased WT-MetaD simulations. Adapting this approach, we calculated the characteristic timescales, τ, for defect formation both within the backbone and on the tips of the 20D stack. In particular, we run multiple infrequent WT-MetaD simulations where the systems undergo a transition from HB–dim = 4 to HB–dim = 0. The unbiased transition time (t) of each transition can be calculated from each WT-MetaD run as:

t = t_WT-MetaD〈e^βV(s(R,t))〉_WT-MetaD (2)

where, V(s(R), t) is the time-dependent bias, the exponential brackets is averaged over the WT-MetaD run, and β is kT⁻¹. The characteristic time scale, τ, of defect formation is then calculated by fitting the cumulative distribution function (CDF) with a Poisson-like cumulative probability:

Fig. 2d shows the CDF profiles of the defect formation at the supramolecular polymer tip (blue curve), and backbone (red curve).

Author contributions

C. P. and A. C. designed the setup and the computational framework. L. R. and P. D. provided the case study. A. C. and C. C. carried out the simulations and the postprocessing of the data. All the authors worked on the interpretation of the results. C. P. supervised the project. All the authors wrote and approved the final manuscript.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). The SI contains the planar arrangement of the dimers forming 20D and 2D supramolecular polymers, the LENS and τSOAP time series corresponding to the interacting centers of each monomers forming the 20D UPy-C6-u supramolecular polymer, and the AA-MD snapshots of three bundled supramolecular polymers at 343 K. See DOI: https://doi.org/10.1039/d5tb01272d. All data and materials related to the simulation MD trajectories and the analysis conducted herein are available at: https://doi.org/10.5281/zenodo.17236935.

Acknowledgements

The authors acknowledge the computational resources provided by the Swiss National Supercomputing Center (CSCS). P. Y. W. D. and L. R. acknowledge support from the Dutch Ministry of Education, Culture, and Science (gravitation programs 024.003.013 and 024.005.020), an NWO VICI grant from the Netherlands Organization for Scientific Research (NWO, VI.C.222.088), and the European Union's Horizon research and innovation program under grant agreement 101079482 (‘SUPRALIFE’).

References

1. I. I. Cisse, I. Izeddin, S. Z. Causse, L. Boudarene, A. Senecal, L. Muresan, C. Dugast-Darzacq, B. Hajj, M. Dahan and X. Darzacq, Science, 2013, 341, 664–667.
2. G. Askarieh, M. Hedhammar, K. Nordling, A. Saenz, C. Casals, A. Rising, J. Johansson and S. D. Knight, Nature, 2010, 465, 236–238.
3. D. Son and Z. Bao, ACS Nano, 2018, 12, 11731–11739.
4. Y. Li, J. Zhu, Z. Zhang, J. Wei, F. Wang, G. Meisl, T. P. Knowles, E. H. Egelman and F. A. Tezcan, Nat. Chem. Biol., 2025, 1–11.
5. H. Shen, E. M. Lynch, S. Akkineni, J. L. Watson, J. Decarreau, N. P. Bethel, I. Benna, W. Sheffler, D. Farrell and F. DiMaio, et al., Nat. Nanotechnol., 2024, 1–6.
6. W. Zhu, H. Wang, B. Feng, G. Liu, Y. Bian, T. Zhao, Q. Wang and X. Weng, Adv. Sci., 2025, 2406479.
7. M. J. Webber and M. W. Tibbitt, Nat. Rev. Mater., 2022, 7, 541–556.
8. A. Bayón-Fernández, A. Méndez-Ardoy, C. Alvarez-Lorenzo, J. R. Granja and J. Montenegro, J. Mater. Chem. B, 2023, 11, 606–617.
9. T. Aida, E. W. Meijer and S. Stupp, Science, 2012, 335, 813–817.
10. Y. Cheng, E. Hirano, H. Wang, M. Kuwayama, E. W. Meijer, H. Huang and T. Aida, Science, 2024, 386, 875–881.
11. I. Insua, A. Cardellini, S. Daz, J. Bergueiro, R. Capelli, G. M. Pavan and J. Montenegro, Chem. Sci., 2023, 14, 14074–14081.
12. I. Insua and J. Montenegro, J. Am. Chem. Soc., 2019, 142, 300–307.
13. B. Zhu, Y. Cai, L. Zhou, L. Zhao, J. Chen, X. Shan, X. Sun, Q. You, X. Gong and W. Zhang, et al., Nat. Commun., 2025, 16, 687.
14. S. Bernhard and M. W. Tibbitt, Adv. Drug Delivery Rev., 2021, 171, 240–256.
15. A. Cardellini, F. Jiménez-Ángeles, P. Asinari and M. Olvera de la Cruz, ACS Nano, 2021, 15, 16139–16148.
16. J. F. van Sprang, J. G. Aarts, M. G. Rutten, L. Rijns, B. M. Tiemeijer, M. J. Schotman and P. Y. W. Dankers, Adv. Funct. Mater., 2024, 34, 2404786.
17. C. Zhao, X. Li, X. Han, Z. Li, S. Bian, W. Zeng, M. Ding, J. Liang, Q. Jiang and Z. Zhou, et al., Nat. Commun., 2024, 15, 1488.
18. N. T. P. Vo, T. U. Nam, M. W. Jeong, J. S. Kim, K. H. Jung, Y. Lee, G. Ma, X. Gu, J. B.-H. Tok and T. I. Lee, et al., Nat. Commun., 2024, 15, 3433.
19. Z. Lei and P. Wu, Nat. Commun., 2018, 9, 1134.
20. C. Shen, X. Qiu, P. Zhang, J. Liu, Z. Zhang, B. Dong, H. Liu, C. Huang, J. Huang and X. Cui, Chem. Eng. J., 2025, 504, 158709.
21. H. Zhang and Z. Guo, Nano Today, 2023, 51, 101933.
22. S. Cantekin, T. F. de Greef and A. R. Palmans, Chem. Soc. Rev., 2012, 41, 6125–6137.
23. P. Besenius, G. Portale, P. H. Bomans, H. M. Janssen, A. R. Palmans and E. W. Meijer, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 17888–17893.
24. A. Demenev, S. H. Eichhorn, T. Taerum, D. F. Perepichka, S. Patwardhan, F. C. Grozema, L. D. Siebbeles and R. Klenkler, Chem. Mater., 2010, 22, 1420–1428.
25. X. Guo, S. Wang, V. Enkelmann, M. Baumgarten and K. Müllen, Org. Lett., 2011, 13, 6062–6065.
26. R. M. Da Silva, D. Van Der Zwaag, L. Albertazzi, S. S. Lee, E. W. Meijer and S. I. Stupp, Nat. Commun., 2016, 7, 11561.
27. H. Cui, M. J. Webber and S. I. Stupp, Pept. Sci., 2010, 94, 1–18.
28. M. P. Hendricks, K. Sato, L. C. Palmer and S. I. Stupp, Acc. Chem. Res., 2017, 50, 2440–2448.
29. S. A. H. Jansen, E. Weyandt, T. Aoki, T. Akiyama, Y. Itoh, G. Vantomme, T. Aida and E. W. Meijer, J. Am. Chem. Soc., 2023, 145, 4231–4237.
30. S.-P. Wang, W. Lin, X. Wang, T.-Y. Cen, H. Xie, J. Huang, B.-Y. Zhu, Z. Zhang, A. Song and J. Hao, et al., Nat. Commun., 2019, 10, 1399.
31. E. Weyandt, L. Leanza, R. Capelli, G. M. Pavan, G. Vantomme and E. W. Meijer, Nat. Commun., 2022, 13, 248.
32. R. P. Sijbesma, F. H. Beijer, L. Brunsveld, B. J. B. Folmer, J. H. K. K. Hirschberg, R. F. M. Lange, J. K. L. Lowe and E. W. Meijer, Science, 1997, 278, 1601–1604.
33. F. H. Beijer, R. P. Sijbesma, H. Kooijman, A. L. Spek and E. W. Meijer, J. Am. Chem. Soc., 1998, 120, 6761–6769.
34. P. Y. W. Dankers, M. C. Harmsen, L. A. Brouwer, M. J. Van Luyn and E. W. Meijer, Nat. Mater., 2005, 4, 568–574.
35. L. Rijns, J. W. Peeters, S. I. S. Hendrikse, M. E. J. Vleugels, X. Lou, H. M. Janssen, E. W. Meijer and P. Y. W. Dankers, Chem. Mater., 2023, 35, 8203–8217.
36. R. E. Kieltyka, A. C. H. Pape, L. Albertazzi, Y. Nakano, M. M. C. Bastings, I. K. Voets, P. Y. W. Dankers and E. W. Meijer, J. Am. Chem. Soc., 2013, 135, 11159–11164.
37. D. J. Wu, M. G. T. A. Rutten, J. Huang, M. J. G. Schotman, J. F. van Sprang, B. M. Tiemeijer, G. M. ter Huurne, S. P. W. Wijnands, M. Diba and P. Y. W. Dankers, Macromolecules, 2024, 57, 6606–6615.
38. L. Rijns, H. Duijs, R. P. M. Lafleur, R. Cardinaels, A. R. A. Palmans, P. Y. W. Dankers and L. Su, Biomacromolecules, 2024, 25, 4686–4696.
39. H. Fu, J. Huang, J. J. van der Tol, L. Su, Y. Wang, S. Dey, P. Zijlstra, G. Fytas, G. Vantomme and P. Y. W. Dankers, et al., Nature, 2024, 626, 1011–1018.
40. A. F. Vrehen, J. F. van Sprang, M. J. Schotman and P. Y. W. Dankers, Mater. Today Bio, 2024, 26, 101021.
41. L. Rijns, M. B. Baker and P. Y. W. Dankers, J. Am. Chem. Soc., 2024, 146, 17539–17558.
42. C. M. Wallace, M. M. Rovers, R. Bellan, M. G. Rutten, A. Seddon, M. J. Dalby, P. Y. W. Dankers and D. J. Adams, J. Mater. Chem. B, 2024, 12, 9283–9288.
43. L. Rijns, M. G. T. A. Rutten, R. Bellan, H. Yuan, M. L. Mugnai, S. Rocha, E. del Gado, P. H. J. Kouwer and P. Y. W. Dankers, Sci. Adv., 2024, 10, eadr3209.
44. A. L. de Marco, D. Bochicchio, A. Gardin, G. Doni and G. M. Pavan, ACS Nano, 2021, 15, 14229–14241.
45. C. Perego, L. Pesce, R. Capelli, S. J. George and G. M. Pavan, ChemSystemsChem, 2021, 3, e2000038.
46. M. Crippa, C. Perego, A. L. de Marco and G. M. Pavan, Nat. Commun., 2022, 13, 2162.
47. K. K. Bejagam, G. Fiorin, M. L. Klein and S. Balasubramanian, J. Phys. Chem. B, 2014, 118, 5218–5228.
48. M. Garzoni, M. B. Baker, C. M. A. Leenders, I. K. Voets, L. Albertazzi, A. R. A. Palmans, E. W. Meijer and G. M. Pavan, J. Am. Chem. Soc., 2016, 138, 13985–13995.
49. M. B. Baker, L. Albertazzi, I. K. Voets, C. M. Leenders, A. R. Palmans, G. M. Pavan and E. W. Meijer, Nat. Commun., 2015, 6, 6234.
50. D. Bochicchio and G. M. Pavan, ACS Nano, 2017, 11, 1000–1011.
51. D. Bochicchio, M. Salvalaglio and G. M. Pavan, Nat. Commun., 2017, 8, 147.
52. A. Gardin, C. Perego, G. Doni and G. M. Pavan, Commun. Chem., 2022, 5, 82.
53. M. Crippa, A. Cardellini, C. Caruso and G. M. Pavan, Proc. Natl. Acad. Sci. U. S. A., 2023, 120, e2300565120.
54. P. Gasparotto, D. Bochicchio, M. Ceriotti and G. M. Pavan, J. Phys. Chem. B, 2019, 124, 589–599.
55. J. Chen, M. Wu, L. Gong, J. Zhang, B. Yan, J. Liu, H. Zhang, T. Thundat and H. Zeng, J. Phys. Chem. C, 2019, 123, 4540–4548.
56. H. Lv, Y. Song, H. Zhang, Y. He, X. Hou, J. Sun and X. Wang, J. Mol. Liq., 2024, 401, 124750.
57. T. Wu, Y. Wang, R. Zou, H. Tan, Q. Fu, Y. Liu and M. Ding, Polymer, 2024, 303, 127116.
58. P. Y. W. Dankers, T. M. Hermans, T. W. Baughman, Y. Kamikawa, R. E. Kieltyka, M. Bastings, H. M. Janssen, N. Sommerdijk, A. Larsen and M. Van Luyn, et al., Adv. Mater., 2012, 24, 2703–2709.
59. S. I. Hendrikse, S. P. Wijnands, R. P. Lafleur, M. J. Pouderoijen, H. M. Janssen, P. Y. W. Dankers and E. W. Meijer, Chem. Commun., 2017, 53, 2279–2282.
60. P. Tiwary and M. Parrinello, Phys. Rev. Lett., 2013, 111, 230602.
61. C. Caruso, A. Cardellini, M. Crippa, D. Rapetti and G. M. Pavan, J. Chem. Phys., 2023, 158, 214302.
62. C. Caruso, M. Crippa, A. Cardellini, M. Cioni, M. Perrone, M. Delle Piane and G. M. Pavan, PNAS Nexus, 2025, 4, pgaf038.
63. W. P. J. Appel, G. Portale, E. Wisse, P. Y. W. Dankers and E. W. Meijer, Macromolecules, 2011, 44, 6776–6784.
64. A. Barducci, G. Bussi and M. Parrinello, Phys. Rev. Lett., 2008, 100, 020603.
65. A. Cardellini, M. Crippa, C. Lionello, S. P. Afrose, D. Das and G. M. Pavan, J. Phys. Chem. B, 2023, 127, 2595–2608.
66. M. Crippa, A. Cardellini, M. Cioni, G. Csányi and G. M. Pavan, Mach. Learn.: Sci. Technol., 2023, 4, 045044.
67. M. D. Hanwell, D. E. Curtis, D. C. Lonie, T. Vandermeersch, E. Zurek and G. R. Hutchison, J. Cheminf., 2012, 4, 1–17.
68. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, Gaussian 16 Revision C.01, 2016, Gaussian Inc., Wallingford CT.
69. C. I. Bayly, P. Cieplak, W. Cornell and P. A. Kollman, J. Phys. Chem., 1993, 97, 10269–10280.
70. J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman and D. A. Case, J. Comput. Chem., 2004, 25, 1157–1174.
71. J. Wang, W. Wang, P. A. Kollman and D. A. Case, J. Mol. Graph. Modell., 2006, 25, 247–260.
72. M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess and E. Lindahl, SoftwareX, 2015, 1, 19–25.
73. P. Mark and L. Nilsson, J. Phys. Chem. A, 2001, 105, 9954–9960.
74. G. Bussi, D. Donadio and M. Parrinello, J. Chem. Phys., 2007, 126, 014101.
75. M. Bernetti and G. Bussi, J. Chem. Phys., 2020, 153, 114107.
76. M. Parrinello and A. Rahman, J. Appl. Phys., 1981, 52, 7182–7190.
77. M. Bonomi, D. Branduardi, G. Bussi, C. Camilloni, D. Provasi, P. Raiteri, D. Donadio, F. Marinelli, F. Pietrucci, R. A. Broglia and M. Parrinello, Comput. Phys. Commun., 2009, 180, 1961–1972.
78. G. A. Tribello, M. Bonomi, D. Branduardi, C. Camilloni and G. Bussi, Comput. Phys. Commun., 2014, 185, 604–613.
79. D. Reynolds, in Gaussian Mixture Models, ed. S. Z. Li and A. Jain, Springer US, Boston, MA, 2009, pp. 659–663.
80. P. Tiwary and M. Parrinello, J. Phys. Chem. B, 2015, 119, 736–742.
81. M. Salvalaglio, P. Tiwary and M. Parrinello, J. Chem. Theory Comput., 2014, 10, 1420–1425.

---