{'fraction_non_alphanumeric': 0.028364276180075233, 'fraction_numerical': 0.013317558548719816, 'mean_word_length': 5.039758153169659, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Accelerated Atomistic Simulations of a Supramolecular Polymer in Water

All-atom molecular dynamics has been recently proven a useful tool for the study of supramolecular polymers. While the high resolution offered by the atomistic models may allow for deep comprehension of the assembled structure, obtaining a reliable equilibrated configuration for these soft assemblies in aqueous solution remains a challenging task mainly due to the complexity of the atomistic systems and the long simulation time needed for their equilibration. Here we have tested two well-known advanced simulation methods (Accelerated Molecular Dynamics and Metadynamics) on an atomistic model of water-soluble 1,3,5-benzenetricarboxamide (BTA) supramolecular polymer and have compared the obtained data to classical molecular dynamics. Our results show that both techniques can be very useful to accelerate the equilibration of the supramolecular fiber in solution. Moreover, we demonstrate that combining the two methods in opportune way allows to take advantage of the strong points of both, providing an equilibrated configuration for the BTA supramolecular polymer in a faster and more reproducible way compared to classical simulations. The versatility of this approach suggests that the latter can be adapted to simulate a variety of supramolecular polymers as well as different types of supramolecular assemblies.


INTRODUCTION
Supramolecular polymers, where the monomers are interconnected via non-covalent interactions, have recently received notable attention thanks to their dynamic and adaptive properties. [1][2][3] These selfassembled structures possess self-healing behavior, stimuli responsiveness and dynamic bio-inspired properties that are very promising for the creation of novel advanced materials for bio-and nanoapplications. 4 In this context, gaining deep understanding of the non-covalent interactions between the monomers is a first fundamental step toward the creation of supramolecular polymers with controlled properties. 5,6 However, the experimental characterization of these soft assemblies is extremely difficult, especially in water, because of their reduced size and the limited contrast offered in solution. This produces a general lack of molecular-level insight into their structure, while in most cases it is almost impossible to obtain from the experiments information on the complex interplay between the different types of interactions (H-bonding, electrostatic, van der Waals, hydrophobic, etc.) controlling the assembly in the aqueous environment.
Given these clear limitations, molecular simulations have recently become an important tool to investigate the structure and properties of supramolecular polymers. [6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22] An extensively studied case is constituted by 1,3,5-benzenetricarboxamide (BTA) supramolecular polymers. The BTA monomers directionally self-assemble into supramolecular fibers due to three-fold hydrogen bonding between the amides and stacking of cores(see Figure 1a-c). 23 Previous computational studies based on density functional theory (DFT) focused on the cooperativity of inter-monomer H-bonding in small BTA stacks in the gas-phase. 9,10,19 H-bonding and dipole amplifications into stacks of BTA monomers with short side chains have also been studied by all-atom classical Molecular Dynamics (MD) simulations in organic solvents. 11 A coarse-grained model of such reduced size BTA derivatives has also been recently reported and tested, which allowed to study BTA self-assembly in nonane. 20 While all these computational efforts provided important insights into the BTA assembly in intrinsically ordered conditions (i.e., the apolar solvent and the semi-rigid nature of the monomers favor ordered stacking), complexity increases enormously when moving to water, mainly due to the more complex and dynamic structure of the watersoluble monomers (Figure 1a,b) and to the strong hydrophobic effects that compete with intra-monomer H-bonding. 21 In 2011 the group of Schatz studied a supramolecular polymer composed of peptide amphiphilic monomers in explicit water by means of classical MD simulations. 7 Their MD results provided a detailed picture of the structure of the polymer and H-bonding between the monomers in the assembly. Recently, our group used all-atom MD simulations to study a BTA supramolecular polymer in explicit water ( Figure 1). This model of an infinite (achiral) BTA fiber was compared to a second polymer variant, where the addition of a single methyl group into the side chains of the BTAs (making the monomers chiral) imparted preferential helicity to the supramolecular fiber. 8 These atomistic models provided important insight into the assembly structure, stability, monomer-monomer interactions, competition of hydrophobic effects vs. H-bonding and penetration of water into the fiber structure. 8 The MD results were also found in optimal agreement with the experiments, providing evidence for stable H-bonding between the BTA cores and structural and energetic data useful to rationalize experimentally observed differences in the structure, stability and dynamics of the achiral and chiral BTA supramolecular polymers. In general, due to the lack of a precise experimental structure for these complex statistical materials, the first fundamental simulation challenge is to obtain from the models a reliable structure for the supramolecular polymer in water and in equilibrium conditions. Typically, atomistic MD simulations start from an initial reasonable configuration for the supramolecular assembly 7,8 created based on the available information from experiments or from higher levels of theory (i.e., DFT). 9 The starting configuration is then relaxed by means of MD to obtain an equilibrated structure for the supramolecular polymer. 7,8,26 The latter is a non-trivial task, given the large size and the complexity of these atomistic systems. Although simulation times of hundreds of nanoseconds are nowadays accessible also for the 6 all-atom MD simulation of such large systems (size also exceeding a few 100k atoms), critical points remain (i) the slow dynamics of the supramolecular structure to reach an equilibrated configuration in the solvent and (ii) the intrinsic sampling limitations of classical MD in studying such complex structures, where during the MD run the assembly can remain trapped in local energy minima. For these reasons, reaching the "true equilibrium" configuration for such complex systems by means of all atom models and classical MD still represents a rather challenging goal.
A possible approach to overcome the aforementioned limitations is to coarse-grain the system, simplifying the structural description of the monomers to be able to reach longer time scales. Recently, coarse graining has been successfully used to study a variety of supramolecular systems. 13,[16][17][18]20 In this work we investigate a different approach: the use of advanced sampling techniques to accelerate the equilibration of complex all atom supramolecular systems. Two popular advanced simulation methods, Accelerated Molecular Dynamics (aMD) 27 and Metadynamics (MetaD) 28,29 , have been tested on an atomistic model of water-soluble BTA supramolecular polymer ( Figure 1) that was previously investigated by means of classical MD. 8 27 has proven an effective method to enhance the conformational sampling of proteins and dendrimers in water. [30][31][32][33][34] MetaD has previously proven useful for a vast variety of systems, including for instance the study of complex synthetic branched macromolecules in solution (i.e., dendrimers). 35,36 Importantly, our data demonstrate that these approaches can be useful to (i) accelerate the equilibration dynamics of a supramolecular system and to (ii) increase the confidence in the reliability of the equilibrium configuration obtained with these advanced simulation techniques.

RESULTS AND DISCUSSION
Our previous all-atom MD simulations of the BTA supramolecular polymer shown in Figure 1e started from an initial model for the BTA fiber consisting of 48 initially extended BTA monomers stacked into a simulation box filled with explicit TIP3P water molecules. 37 This system forming an infinite supramolecular polymer along the main fiber axis through periodic boundary conditions (Figure 1e, left) was relaxed in NPT conditions using anisotropic pressure scaling to allow free rearrangement of the supramolecular fiber during the equilibration MD. 8 Herein, we will identify this reference classical MD simulation as MD0. Early during MD0 the side chains of the BTAs collapsed around the hydrophobic BTA cores to minimize the contacts with water (primary folding), while the whole initially extended fiber started to fold (secondary folding) to further reduce the hydrophobic surface exposed to the solvent.
In this case, 400 ns of MD0 were needed to reach a configuration (Figure 1e, right) in which several key structural and energetic parameters used to characterize the fiber (fiber length, monomer-monomer interaction energy, radial distribution functions g(r) of the cores, H-bonds, etc.) showed good convergence in the MD regime. 8 Theoretical small-angle X-ray scattering (SAXS) profiles obtained from the equilibrated phase MD simulation (the last 100ns of MD0) were also found in good agreement with the experimental ones. 8 However, while apparently fiber folding came to an end and the fiber appeared as fully equilibrated during MD0, the structure still presented hydrophobic patches exposed to the solvent and a rather wrinkled and irregular surface. The average fiber cross-sectional radius obtained from MD0 was also found slightly smaller than the SAXS experimental one (3.1±0.2 nm). Considering the limited sampling efficiency of classical MD, especially for very complex systems, these observations suggested that fiber folding could have possibly been still incomplete.
Accelerated MD simulations of a BTA supramolecular polymer in water. The fundamental idea behind the aMD technique is to raise the potential wells of the system to allow faster crossing of the energy barriers dividing the different configurations that the system can access. 27 This is done by adding 8 an energetic bias to the global potential of the whole system, to the sole dihedrals bonded terms or to both (global potential plus an extra bias term dedicated to dihedrals) when these are below certain threshold values. Here we have used the aMD method as implemented in AMBER 14 38 to simulate the BTA supramolecular polymer of Figure 1e and we have compared the results with those obtained from classical MD0. We have applied the bias to both the whole potential and the dihedrals, where the bias 9 parameters have been calculated using the standard approach used in previous aMD studies reported in literature (see Materials and Methods). 27,[30][31][32][33] First we ran an aMD simulation starting from the same initially extended configuration of the BTA fiber used as a starting point for MD0 (Figure 1e, left). We will identify this aMD simulation aMD0.
During aMD0 the BTA fiber undergoes similar primary and secondary folding than in MD0, but the process was sensibly faster. Interestingly, fiber folding did not stop at the MD0 level, but proceeded further for ≈200 ns of aMD0 before coming to an end. The aMD0 run was then prolonged for additional 200 ns to ensure full equilibration in the aMD regime, during which we did not observe any relevant structural or energetic change. The final structure of the fiber produced by aMD0 is shown in Figure 2a.
Compared to MD0, aMD0 produced a more compact/uniform and nearly cylindrical fiber while substantially preserving the order of the BTA cores. In fact, the latters are found as stably stacked into a rather continuous serpentine with a few punctual defects emerging ( Figure 2a).
We investigated to what extent the final equilibrated state for the BTA fiber produced by aMD0 was different from that obtained from MD0. In particular, we focused on structural features as well as on the main factors controlling the BTA assembly in water (hydrophobic effects, H-bonding, etc.). 8 A good indicator of the hydrophobic effects involved in fiber folding can be obtained from the solvent accessible surface area (SASA) of the BTA supramolecular polymer -i.e., the stronger the SASA shrinkage at the equilibrium (after fiber folding), the better the fiber rearranged its structure to minimize the contacts with water. To elucidate the differences between aMD0 and MD0, we plotted the two simulation trajectories in a 2D space (Figure 3a: MD0, black; aMD0, cyan), where the x and y axes respectively report the average radius and the SASA of the BTA fiber during the runs. From the 2D plots of Figure 3a it is clear that in both simulations fiber folding is accompanied by strong SASA reduction, which highlights the dominant character of the hydrophobic effects. Accordingly, augmented folding produces an increase in the average fiber radius. In both cases the systems fall in a minimum where fiber's SASA and radius appear as equilibrated and oscillate around their average values. However, the two MD0 and aMD0 minima are significantly separated in the radius-SASA space -respectively, the equilibrated fiber configuration produced by MD0 lies between ≈220-260 nm 2 of SASA and ≈2.1-2.3 nm of radius, while that from aMD0 between ≈160-190 nm 2 of SASA and radius ≈2.9-3.1 nm. The stronger SASA shrinkage seen in aMD0 suggests that the BTA fiber was able to better rearrange its structure to further reduce the contacts with the solvent compared to MD0. The final fiber radius obtained from aMD0 was also found to fit very well with the SAXS experiments (3.1±0.2 nm, in grey in Figure 3a), 8 whereas that obtained from MD0 was found smaller -a point in favor of aMD.
In principle, proper reweighting 39 of the equilibrated phase aMD0 trajectory would allow to estimate the free-energy surface (FES) of the BTA fiber in the radius-SASA space. While it is known that obtaining precise reweighting of the aMD trajectories can be problematic for such complex systems, 40 it is worth to stress that here we are not interested in obtaining a rigorous estimation of the FES (which, by the way, seems to possess a unique global minimum in the radius-SASA space, as confirmed by the constant average SASA and radius in the last part of aMD0 and by the other simulations performed -see later on), but just to speed-up the equilibration of the system. While the aMD method has been here used with this sole purpose, we delegated more quantitative energetic and structural analyses on the equilibrated BTA fibers to additional 100 ns of unbiased standard MD simulation (MDpost) starting from the end of each biased simulation. During production MDpost the fiber did not change its structure anymore, indicating good stability for the final configuration produced by aMD0. The same approach has been followed also for all other enhanced sampling simulations tested in this work (see below). We calculated the number of contacts between water and the hydrophobic regions of the fiber (BTA core and dodecyl spacers). Coherent with the SASA data, these are found reduced in aMD0 compared to MD0 (Figure 3b). Another parameter that goes in the same direction is the radial distribution function g(r) of the water molecules calculated respect to the BTA cores (Figure 4a). This data represents the probability of finding water molecules within a certain distance from the cores. We observe a reduction of ≈50% of water penetration within 1 nm from the BTA cores in aMD0 compared to MD0. All these data suggest that a key driving force for the augmented folding seen in aMD0 (and also in all other biased simulations, see below) lies in the hydrophobic effects. Another important analysis concerns the level of order in the BTA core stacking, which is a good descriptor of the stability and persistency of the supramolecular assembly. 8,21,26 This can be monitored by calculating the radial distribution function, g(r), of the BTA cores as a function of the inter-core distance, which peaks (at distances multiple of the stacking one, ≈0.35 nm) are the fingerprints of core stacking. Compared in Figure 4b, the core-core g(r) plots are found similar in MD0 (black) and aMD0 (cyan), indicating that the core stacking order has been substantially preserved in aMD0. The first g(r) peak is however slightly reduced, consistent with the presence of some "broken" points emerging along the fiber early during the aMD0 folding.
We have compared the energies of the equilibrated BTA fiber as obtained from MD0 and aMD0. For what concerns the total potential energy of the fiber, only small differences could be detected. Consistent with what said above for the core-core g(r) plots, this indicates that the fiber as seen during the last 100 ns of MD0 was already almost equilibrated from the point of view of the monomer-monomer interactions and overall persistency/stability of the BTA assembly. In Table 1 we report the main features of the BTA supramolecular fiber that are found improved in aMD0 compared to MD0. For example, as anticipated, the equilibrated fiber SASA is found reduced in aMD0 compared to MD0 by ≈27% -∆SASA = SASA(aMD0) -SASA(MD0) in Table 1 -, which confirms that the fiber can undergo better compaction during aMD0 compared to MD0 to optimize its hydrophobic/hydrophilic interactions. The electrostatic term of BTA-BTA interaction (∆Eele) is also found more favorable in aMD0 than in MD0 by ≈28% - Table 1: ∆∆Eele, calculated as ∆Eele(aMD0) -∆Eele(MD0), equal to ≈-5.8 kcal/mol per-BTA monomer (see also SI). This is coherent with an overall stabilization of the H-bonding network in the fiber, where the average number of H-bonds per-BTA is seen to increase by ≈20% in the case of aMD0. These observations suggest that under different aspects the equilibrated state reached by the BTA fiber during aMD0 is energetically more optimized compared to that from classical MD0.  [a] ∆∆Eele is the difference between the electrostatic component of the BTA self-assembly energy (∆Eele) calculated from MDpost (following to each biased simulation) and from MD0 (values per-BTA monomer).
[b] ∆H-bonds parameter is the difference between the average number of H-bonds per-BTA calculated from MDpost and from MD0. [c] ∆SASA is the difference between the SASA of the fiber calculated from MDpost and from MD0.
To increase our confidence in the aMD approach, we performed a second aMD simulation (aMDMD) using the same bias as in aMD0 but starting from a different configuration for the fiber: the preequilibrated one obtained from MD0 (Figure 1e, right). In this way we intended to verify if this method was able to produce a similar equilibrated state for the fiber while starting from two different configurations.
In aMDMD we observed once again that the BTA fiber folded further than in MD0. In this case fiber folding reached completion after ≈150 ns of aMD, while the aMDMD run was prolonged for other 250 ns to verify the convergence of the structure. At the equilibrium, aMDMD provided a final radius and an equilibrated SASA for the fiber substantially identical to those obtained from aMD0, as demonstrated by  Table 1). Aside from subtle differences imputable to normal statistical oscillations, all the extracted parameters indicate that aMD0 and aMDMD produced almost equivalent final states for the fiber. This is demonstrated by the core-core g(r) plots (Figure 4a) showing again nearly identical core stacking in aMDMD (red) and in MD0 (black). Moreover, all data from aMDMD pertaining to fiber hydration are substantially identical to those of aMD0. Similar to aMD0, aMDMD also shows analogous improvements compared to MD0 for what pertains to the average number of H-bonds per-BTA and monomer-monomer electrostatic interaction (Table 1), demonstrating the reliability of the aMD method in producing results that are substantially independent from the starting configuration.

Metadynamics.
We decided to test also a different biased method, Metadynamics, 28,29 for two main purposes: (i) to check if a different biased approach could produce similar results and (ii) to try to accelerate further the fiber folding process. In fact, the latter took considerable simulation time to reach the equilibrium (>200 ns) even in the aMD accelerated regime. Metadynamics (MetaD) is a very popular biased approach which main purpose is to allow crossing the free-energy barriers between different states and to enhance the sampling of the conformational space of the system. For instance, MetaD has been recently proven useful to allow a more exhaustive sampling compared to classical MD of complex hyperbranched polymers (dendrimers) in aqueous solution. 35,36 The MetaD method requires the identification of one (or more) collective variable (CV) able to describe the slow degree of freedom of the system that has to be accelerated. However, the identification of the right CVs is known to be a key and often non-trivial point. 29 While MetaD is typically used to obtain an accurate description of the FES of the simulated system, we stress again that here our main scope in using MetaD was not to obtain a precise sampling of the FES of the BTA fiber, but rather to speed-up in reliable way the slowest macroscopic process (fiber folding) leading to the equilibration of the supramolecular BTA fiber in the solvent. After a few tests, we identified the best suited CV to our aim: the maximum z distance between the centers of mass of the BTA cores, strictly connected to the fiber elongation/folding along its main axis (details on the MetaD parameters used in this work can be found in the Materials and Methods section).
The results of our first tests with MetaD were surprisingly good for what pertained to the folding speed.
Starting from the initially extended configuration for the BA fiber (Figure 1e, left), the folding was strongly accelerated compared to aMD and the fiber z dimension (fiber length) went straight down to values similar to the equilibrium ones of the aMD runs, but in well reduced time. Only ≈5 nanoseconds of MetaD were sufficient to produce comparable fiber folding and average cross-sectional radius to those obtained from aMD0 and aMDMD in hundreds of nanoseconds (150 ns and 200 ns respectively).
Despite the notable folding acceleration produced by MetaD, the core-core stacking (g(r)) was not destabilized, but on the contrary it was found better preserved than in the aMD simulations and nearly identical to MD0. This indicated that the MetaD method allowed to fold the fiber in faster but at the same time more delicate way than aMD, at least in the earlier stages of fiber folding (i.e., very far from the equilibrium). However, one disadvantage of this MetaD approach was that fiber folding was so fast that no sufficient time was left for full reorganization of the BTA lateral chains. From this viewpoint, the final configuration is still reminiscent of the starting one, and this actualized in a non-uniform rough surface and sub-optimal SASA that turned out very difficult to be optimized in efficient way at such accelerated regime and using this sole CV. , of (a) the water molecules calculated respect to the BTA cores (water penetration into the fiber) and of (b) the BTA cores as a function of the core-core distance.
A combinined MetaD-aMD approach. In light of these observations, we decided to take advantage of the strong points of both methods using them in a combined protocol. Namely, (i) we first used MetaD to accelerate the "gross" folding of the initially extended fiber and then (ii) we used aMD to promote fiber rearrangement and to optimize the folded state. As previously done for aMD, we conducted two simulations (MetaD-aMD0 and MetaD-aMDMD) starting from two different configurations for the fiber: the initially extended BTA fiber and the pre-equilibrated configuration from MD0 reported in Figure   1,left and Figure 1,right respectively. As a criterion to end the first phase (MetaD) and chose a starting point for the second one (aMD) we used a very simple one: we took the configuration of the system immediately after the abrupt change in the slope of the MetaD hills seen during the run (see SI) -namely, the point where the fiber reaches the maximum folding level along the z axis. Obtaining these folded yet non-optimized fibers required just ≈5-7 ns of MetaD simulation. Then, these pre-folded configurations underwent 150 ns of aMD, sufficient to optimize the BTA side chains reaching full equilibration of the BTA fibers from the structural and energetic points of view (see SI).
The trajectories obtained from the combined MetaD-aMD0 and MetaD-aMDMD simulations are also reported in Figure 3a (green and yellow respectively). Qualitatively, these show that the fibers follow slightly different paths during these equilibration runs compared to pure aMD simulations, first increasing the fiber average radius (MetaD: fast folding), and then optimizing their SASA (aMD: overall rearrangement). Seen in Figure 3a, the final fiber radii obtained from MetaD-aMD0 and MetaD-aMDMD were found well compatible with the experimental one (grey region) and with the aMD0 and aMDMD results. The final structure of the fiber was very similar too (Figure 2b: MetaD-aMDMD), and detailed analyses of these simulations showed similar improvements compared to MD0 than in aMD0 and aMDMD (Table 1). The ∆SASA (hydrophobic effects) and ∆∆Eele (BTA-BTA electrostatic interactions) were found also slightly improved compared to the pure aMD simulations (by an additional ≈10%). Also, consistently with what said above, the first peak in the g(r) of the BTA cores as calculated from MetaD-aMD0 and MetaD-aMDMD was found in general slightly higher than the corresponding aMD0 and aMDMD simulations (Figure 4b).
Considering statistical fluctuations, we should carefully say that the two approaches produce comparable equilibrated configurations for the BTA fiber. However, it is worth noting that the combined MetaD-aMD approach allowed to reduce the overall simulation time necessary to obtain an equilibrated configuration for the fiber by ≈75% compared to pure aMD. Thus, summarizing, every procedure used to accelerate the equilibration of the BTA supramolecular structure produced similar final configurations in terms of energy, SASA, and fiber radius (see  Thanks to the technical advantages offered by the combined MetaD-aMD approach, we could also create and simulate (3MetaD-aMD0) a BTA fiber three-times longer (144 BTA cores) than that simulated before. 8 This simulation was useful to prove that the overall characteristics of the BTA supramolecular polymer extracted from our 48 BTA polymer model were not influenced to some extent by examining a rather reduced section of the fiber (or, in other words, to exclude spurious PBC effects).
This model was built starting from MetaD folding of the initially extended 48 BTA fiber, which was then triplicated (144 BTA) along the z axis and underwent 150 ns of equilibration aMD. The final equilibrated configuration for this larger system is shown in Figure 5. As previously done for all other systems, 100 ns of unbiased MDpost simulation were run for data analyses.
All structural and energetic parameters demonstrate perfect agreement with the smaller 48 BTA systems (see Table 1 and Figure 4), proving that the combined MetaD-aMD approach can be remarkably useful (i) to obtain a reliable and well-reproducible equilibrium configuration for the BTA fiber in water, (ii) to dramatically speed-up fiber folding avoiding entrapment in local energy minima and as a consequence (iii) to simulate larger systems allowing for a more accurate statistical modeling of these supramolecular polymers.

CONCLUSIONS
In this work we have used two advanced simulation approaches, aMD and MetaD, to obtain the equilibrium configuration for an atomistic model of BTA supramolecular polymer in explicit water in efficient and reliable way. Our data highlight the technical advantages offered by both methods in the simulation of a very complex supramolecular system such as the water-soluble BTA fibers studied herein.
In particular, our results demonstrate that a combined MetaD-aMD approach, where MetaD is first used to accelerate the initial fiber folding and aMD is then used to promote optimization of the supramolecular structure, allows for remarkable speed-up of the fiber equilibration process. Comparison of our simulations starting from different configurations for the BTA fiber and using two different biased techniques proves that this approach can produce an equilibrated configuration of the supramolecular BTA polymer in water that is reliable and well reproducible. This offers notable technical advantages respect to classical MD especially in terms of efficiency (reduced simulation time), possibility of studying larger systems and confidence in providing a high-resolution picture for the structure of these complex supramolecular systems in relevant environments and in equilibrium conditions. Our results also suggest that this simulation approaches can be opportunely adapted for the simulation of soft 20 supramolecular assemblies other than supramolecular polymers, such as, for example, small micelles, 41,42 macrocycles and nanocages 43 to name a few.

MATERIALS AND METHODS
The entire simulation work was conducted using the AMBER14 software (classical MD and aMD) and the (PME) was used to treat long range electrostatics. 52 The bias parameters for the whole potential energy (AlphaP) and for the dihedral term (AlphaD) used in all our aMD simulations have been calculated with the standard empirical formulas previously used in other aMD studies of protein systems in aqueous solution: 27,30-33 where in our case each BTA monomer has been considered as counting for 4 residues (one for the central ring + three side chains, see Figure 1) in the calculation of the extra bias for dihedral angles. This assumption has been made in light of the higher complexity and larger size of the water-soluble BTA monomers studied herein compared to standard amino acids composing proteins, on which aMD has been usually tested. In fact, initial tests counting each BTA monomer in the fiber as one single monomer produced negligible additional acceleration due to the extra bias applied on the dihedrals. The threshold values for the whole potential energy (EthreshP) and for the dihedral term (EthreshD) of the simulated systems have been set to: where Epot and Edih are the reference equilibrated values taken from our previous unbiased classical MD0 simulation.
In the MetaD simulations we have used as a single global variable the maximum z distance between the centers of mass of the BTA cores, calculated taking into account PBC, as defined in PLUMED2 as: where β has been set to 0.1 and the index i runs over all the BTA cores (48). We have set the hills height to 0.5 kcal/mol and the sigma to 0.05 nm. Hills were deposited every 500 simulation time steps (i.e., every picosecond).
The average radius of the fiber has been estimated the following way. We have calculated the z component of the gyration tensor in the final state of the fiber, from which it is possible to derive an underestimate of the radius (approximating the shape of the fiber to a gyration cylinder). We have also calculated an upper limit of the fiber radius looking at the maximum coordinate values of the fiber in the x and y direction. We have taken the mean of these two extreme values for the estimation of our average radius, which has been compared to the experimental one (see Figure 3a).

22
The energetic analysis on the BTA fiber -i.e., the electrostatic term of monomer-monomer interaction, ∆Eele -has been conducted according to the MM-PBSA 53,54 approach on the last 100 ns of unbiased MDpost conducted at the end of each simulation as: where Eele(fiber) is the average electrostatic energy of a BTA monomer into the fiber and Eele(monomer) is the electrostatic energy of a monodisperse BTA monomer. Negative values for ∆Eele thus indicate favorable electrostatic contribution for the BTAs to stay assembled in the fiber rather than disassembled in solution (see also SI). The variation between the ∆Eele energy as measured during the 100 ns of unbiased MDpost at the end of each i accelerated simulation and on the last 100 ns of classical MD0 was calculated as: Negative values of ∆∆Eele identify a better electrostatic state for the assembly reached during the accelerated simulations compared to the classical MD0 (see also SI).
The H-bonding analysis has been conducted using the ptraj module of AMBER 14, while SASA have been extracted using the gmx_sasa utility of Gromacs 5.1.2.