How the Choice of Force-Field Affects the Stability and Self-Assembly Process of Supramolecular CTA Fibers.

In recent years, computational methods have become an essential element of studies focusing on the self-assembly process. Although they provide unique insights, they face challenges, from which two are the most often mentioned in the literature: the temporal and spatial scale of the self-assembly. A less often mentioned issue, but not less important, is the choice of the force-field. The repetitive nature of the supramolecular structure results in many similar interactions. Consequently, even a small deviation in these interactions can lead to significant energy differences in the whole structure. However, studies comparing different force-fields for self-assembling systems are scarce. In this article, we compare molecular dynamics simulations for trifold hydrogen-bonded fibers performed with different force-fields, namely GROMOS, CHARMM General Force Field (CGenFF), CHARMM Drude, General Amber Force-Field (GAFF), Martini, and polarized Martini. Briefly, we tested the force-fields by simulating: (i) spontaneous self-assembly (none form a fiber within 500 ns), (ii) stability of the fiber (observed for CHARMM Drude, GAFF, MartiniP), (iii) dimerization (observed for GROMOS, GAFF, and MartiniP), and (iv) oligomerization (observed for CHARMM Drude and MartiniP). This system shows that knowledge of the force-field behavior regarding interactions in oligomer and larger self-assembled structures is crucial for designing efficient simulation protocols for self-assembling systems.

Introduction

Self-assembled fibers have recently drawn attention because of their rich, dynamic behavior, similar to that of many materials occurring in biological systems. Moreover, they often exhibit features that conventional polymers, connected by covalent interactions, do not have.[1] However, there is still little known about the mechanism of the formation of supramolecular fibers. Such knowledge would allow improving control over their structure and function.[2] Therefore, substantial effort is directed toward understanding the self-assembly process. However, self-assembly steps often occur on a temporal and spatial scale beyond the reach of experimental techniques.

Consequently, for processes on short time and length scales, many computational studies are devoted to supramolecular polymers, as shown in recent reviews by Frederix et al.[3] and Bochicchio et al.[4] The main challenge of simulations of supramolecular systems is connecting the small spatial and temporal scales of computational systems to the large spatial and temporal scales of the experiment.[5] Most of the time, the experiment involves large systems that form structures on time scales spanning from nanoseconds to weeks.[3] Such a time scale is often far beyond the capabilities of current computational methods. These spatial and temporal challenges are often mentioned in the literature in the context of molecular simulation[3,5−10] and are the focus of many studies.[3,4,11,12] However, the nature of supramolecular systems (i.e., the molecules are noncovalently connected by the same type of interactions) is such that small errors in a model describing an interaction are amplified by the number of molecules, as has been pointed out in the context of protein modeling.[13] Despite this fact, studies on how different force-fields influence the supramolecular structures are scarce,[3] and mostly limited to aggregation of amyloids[14−18] and surfactants.[19−21]

In this work, we study two currently standard approaches that are employed to give insight into self-assembly: simulation of spontaneous self-assembly starting from randomly distributed molecules[22−25] and simulation starting from a prebuilt model structure of the proposed final assembly.[26−30] As a model example, we use a derivative of 1,3,5-trisamidocyclohexane (CTA; see Figure a), which is known to create long ordered fibers upon self-assembly and for which crystal structure of its analog has been reported.[31] CTA belongs to a large class of supramolecular molecular blocks that form fibers via trifold hydrogen bonding and recently are a subject of intensive studies.[30,32−35] Here, we simulate the self-assembly and stability of CTA fibers using different force-fields (namely: GROMOS,[36] CHARMM general force-field (CGenFF),[37−39] CHARMM Drude,[40−43] General Amber force-field (GAFF),[44] Martini,[45,46] and polarized Martini[47]). In terms of computational performance Martini and polarized Martini are the best (1 ns simulation takes ∼3 min/CPU); CGenFF and GROMOS are approximately 2 orders of magnitude slower (1 ns simulation takes ∼8 h/CPU), and CHARMM Drude is another four times slower (1 ns simulation takes ∼28 h/CPU; see the Supporting Information, SI). This work shows that force-fields capture various aspects of self-assembling systems differently, causing the simulation of the self-assembly processes and the prediction of the stable self-assembled structure(s) to be nontrivial and requiring different strategies.

Figure 1

Simulations of eight molecules in a small simulation box for simulations up to 500 ns for different force-fields. (a) Chemical structure of derivative of 1,3,5-cyclohexanetricarboxamide (CTA). Final snapshots of the simulations for: (b) GROMOS, (c) CHARMM Drude, (d) CGenFF, (e) CGenFF mod., (f) GAFF, (g) Martini, and (h) MartiniP. In none of the simulations we have observed formation of long-range ordered structures. (i) Progression of number of amide–amide hydrogen bonds. (j) Progression of solvent accessible surface area (SASA). The dashed line at 100 ns indicates the change of the x-axis scaling.

Results

Spontaneous Self-Assembly Simulations

The most common way of using molecular dynamics (MD) simulations is to simulate a system of interest for as long as possible. In general, this approach should work because a system that starts from a nonequilibrium situation progresses on an energy landscape and explores it, finding local minima, which ideally corresponds to the experimentally stable structures. For many systems, like proteins, it is hypothesized that the energy landscape has the shape of a funnel;[48] therefore, on long enough simulation, exploring it will lead to the most stable structure being visited most often. In practice, however, one might end up in a local minimum and not escape from it. We have attempted simulation of supramolecular self-assembly by performing 500 ns simulations of 8 molecules in a small simulation box (a 4.5 × 4.0 × 4.1 nm3 dodecahedron; see Table S2). Such length of simulations has been successful in self-assembly simulations of amphiphiles in water.[22,49] We have simulated systems using different force-fields at different levels of resolution: the all-atom model including electronic polarizability CHARMM Drude model[40−43] (due to computational limitations only 67 ns simulated), all-atom CGenFF model,[37−39] the CGenFF model with modified charges obtained by mapping effective charges from the CHARMM Drude model (CGenFF mod.), the united-atom GROMOS model,[36] all-atom GAFF model,[44] the coarse-grained (CG) polarized Martini (MartiniP) model,[47] and its parent coarse-grained model Martini.[45,46] The final snapshots of these simulations are shown in Figure b–h. In all simulations, molecules aggregate into a cluster, as observed by visual inspection and from the solvent-accessible surface area (SASA; see Figure j). However, in none of them, we have observed the formation of long-range ordered structures, as characterized by the number of hydrogen bonds between amides of the CTA, which is shown in Figure i. The most ordered structures were observed in MartiniP, where we could observe small, ordered fragments (dimers and trimers). The variations in both the number of hydrogen bonds and the SASA as a function of simulation time informs about the dynamics of the system. It can be seen that, for GROMOS, a compact and stable structure is obtained, which is reflected in a low value of and small variation in SASA. The structure undergoes slow growth, which can be observed by steady growth of the number of hydrogen bonds over 500 ns simulations. A slightly more flexible structure can be observed in GAFF, for which the 500 ns simulation results in the formation of several dimers. The structures are considerably more flexible in the CGenFF, CGenFF Mod., and CHARMM Drude models. Since the volume of the bead for CG models is different, it is difficult to compare the SASA for them with atomistic force-fields, but the analysis shows that the Martini models also form relatively compact structures.

It is worth noting that the self-assembly mechanism depends on the concentration and influence of water content, which has been previously studied for a similar system.[50] We have studied the effect of water content for self-assembly and observed that quadrupling the size of the system in GAFF improved long-ranged order (see the SI). A larger volume decreases the probability of forming a single unordered cluster, which results in the formation of a smaller cluster that can quickly rearrange into ordered structures. Lastly, it is also worth noting that more realistic simulations should be performed under a constant chemical potential ensemble rather than a constant number of molecules. Unfortunately, such simulations are still challenging.[51]

Fiber Stability

The success of conventional MD simulations relies on the assumption that the experimentally observed structure is in the global minimum on the energy landscape obtained from the simulation. For a long enough time scale, the simulation would lead to visiting all possible configurations of the system, including complete fiber. If a force-field models the interactions perfectly, such fiber would be in the free energy minimum and therefore stable[52] (it is worth noting that this might not be a global minimum due to differences between experimental and computational conditions, in particular, the presence of periodic boundary conditions and a different concentration and ensemble). However, self-assembled systems are challenging to model: even small errors in the parametrization of the force-fields are important since they are multiplied by the number of molecular blocks. Therefore, one of the reasons why spontaneous self-assembly simulations might not lead to an experimentally observed final structure might be that the fibers are not in the energetic minimum for the chosen force-field. The question regarding the true minimum of the force-field prompts the second strategy commonly used in supramolecular structure modeling: assess the stability of a proposed architecture.[26−30] In order to check the stability of the supramolecular fiber, we have created a stack of 24 CTA molecules (see Figure a and Table S3 for computational details) from a known crystal structure of its analogue.[31] Then we have simulated the structures for ∼300 ns or until it collapses (and ∼100 ns for the Drude force-field due to computational cost) using the same force-fields as in the previous section. The final snapshots in the NVT ensemble are presented in Figure b–g. To quantitatively measure the stability of the structure, we analyzed the trajectories by calculating the number of hydrogen bonds between CTA amides (Figure h) and solvent accessible surface area (SASA; Figure i).

Figure 2

Simulation of the supramolecular fiber of CTA obtained from the crystal structure. (a) Starting structure. Single simulations (from 50 ns up to 300 ns) of the fiber in different force-fields in the NVT ensemble: (b) GROMOS, (c) CGenFF, (d) CHARMM Drude, (e) CGenFF with charges obtained from a mapping from CHARMM Drude FF (see main text), (f) GAFF, (g) Martini, and (h) MartiniP. Additionally, for panels b and e, the results for the NPT ensemble are shown since they qualitatively differ from the NVT ensemble. The most stable structures were obtained for CHARMM Drude, GAFF and MartiniP, which can also be seen on graphs of the number of hydrogen bonds per molecule (i) and solvent accessible surface area (SASA; j). For stable fibers, the number of hydrogen bonds and SASA are constant. The dashed line at 100 ns indicates the change of the x-axis scaling.

The fiber in CHARMM Drude, GAFF, and polarized Martini force-fields stays in the ordered structure during the simulations, as can be observed by visual inspection and the constant number of hydrogen bonds (Figure h) and SASA (Figure i). The fiber in the GROMOS force-field stays stable for a relatively long time (∼130 ns) when it collapses. However, even in the collapsed structure, it retains partial order (as seen from the number of hydrogen bonds). With the standard CGenFF, the fiber collapses: most of the hydrogen bonds are immediately broken, and the structure rearranges into an unordered, compact agglomerate (see SASA in Figure i). Although CHARMM Drude is able to sustain a stable structure, it is computationally the most expensive force-field. Therefore, we checked if we can use the standard CGenFF with modified charges obtained from the CHARMM Drude force-field that reflect the average polarization of the chemical groups in the stable assembly (see the SI for details of backmapping of the charges). Although the fiber structure using this modified force-field with effective charges (CGenFF mod.) is more stable than in standard CGenFF, it collapses in the course of the simulation. For coarse-grain force-fields, it seems that reproduction of directional interactions of amide groups is necessary to model a stable fiber: for Martini the structure collapses, whereas for MartiniP it stays stable. This contrasts with the results of Bochicchio et al. for a BTA molecule, for which Martini and MartiniP yield similar outcomes.[53]

For GROMOS and CGenFF Mod., simulations in the NPT ensemble (semi-isotropic) are more stable (shown on the bottom of Figure b,e; see the SI). The fiber in GROMOS remains stable over 500 ns, but it does not maintain a straight shape. The increased stability in the NPT ensemble results from suppressing fluctuations of the distances between monomers (in the direction of fiber’s main axis). The fixed length of the simulation box in the NVT ensemble might destabilize fibers because it does not allow correlating fluctuations in total length with fluctuations in orientation and packing of the monomers. For this reason, simulation of fibers should be performed in a semiisotropic NPT ensemble if possible. It is worth noting that simulations of fiber in isotropic NPT ensemble result in the breaking of the fiber in all force-fields except MartiniP (see the SI).

Energy of the Fiber in Vacuum

In order to investigate which interactions are responsible for the stability of the fiber, we analyzed the addition of molecules to a growing small stack. The enthalpy of the creation of dimer, trimer, etc., from monomers in a vacuum is presented in Figure a. All force-fields show a strong cooperative effect. The gain in the energy per molecule upon addition of a further monomer is large for dimers until tetramers and substantially slows down for pentamers and longer stacks. Although these trends hold for all force-fields, the differences between enthalpies for different force-fields also grow with the size of the stack. The differences seem significant, and the addition of a monomer to an 11-mer reaches ∼70 kJ/mol for MartiniP and CGenFF (and ∼50 kJ/mol for all-atomistic force-fields, i.e., between CGenFF and CGenFF Mod.). However, the analysis is done only for a single conformation, and the differences might change upon a proper sampling of the fiber conformations. Then, we analyzed the contributions of van der Waals, electrostatic and bonded interactions to the binding energy; these are presented in Figure b for the CHARMM Drude and GROMOS models. Analysis for the specific interactions shows that a cooperative effect is present in both Coulomb and van der Waals interactions. Interestingly, the most important contribution for the GROMOS force-field comes from the van der Waals interaction (Lennard-Jones, L-J), but for the CHARMM Drude force-field it comes from electrostatic interactions. Other atomistic nonpolarizable force-fields give similar results (see the SI). For coarse-grain force-fields the main contribution comes from L-J interactions. It is important to note that Lennard-Jones interactions are short-ranged (they decay with r–6, where r is the distance between atoms) and therefore can only be weakly directional, but Coulomb interactions are long-range interactions (they decay with r–1) and therefore a combination of different partial charges can make these interactions highly directional. As a result, the driving force for self-assembly in simulation depends on the choice of force-field: for the force-fields studied here, for the coarse-grained ones the strongest interactions are L-J, for nonpolarizable all-atomistic ones mainly L-J, with an important contribution of electrostatic interactions, and for the polarizable all-atomistic one the opposite is the case, i.e., mainly electrostatic with the important contribution of L-J interactions.

Figure 3

(a) Enthalpy of creation of dimer, trimer, etc. per molecule. In all force-fields studied here a cooperative effect upon stacking is observed. (b) Decomposition of enthalpy of creation dimers, trimers, etc. for GROMOS (large circles) and CHARMM Drude (small circles) force-fields (see others in the SI). Although the overall enthalpies are similar, the contributions from coulomb (red) and van der Waals (orange) interactions are reversed. In Drude FF coulomb interactions are the strongest, whereas in GROMOS the van der Waals are strongest. Bonded interactions (blue) do not contribute to the binding.

Energy Landscape of Fiber

We have shown that in some force-fields the fiber structure is stable; however, simulation from randomly distributed molecules did not lead to the formation of small ordered structures. To investigate why simulations do not lead to stable structures, we have performed simulations of creating a dimer from two free molecules and a pentamer (small ordered fiber) from a free monomer and a tetramer. Division of the self-assembly process into these single steps allows studying it on a much shorter time scale, hopefully accessible by conventional MD. We have run multiple independent 10 ns simulations of two free molecules (see Figure a–d and Table S5 for computational details). To simulate pentamerization we have run multiple independent 10 ns simulations of a tetramer kept stable by position restraints (on the cyclohexane core atoms) with one additional free molecule (without restraints; see Figure e–h and Table S5). Here, we show simulations for four different force-fields: GROMOS and CHARMM Drude FF, GAFF, and MartiniP, which were able to maintain the stable fiber. The results for Martini, CGenFF, and CGenFF Mod. are presented in the SI. We have calculated a 2D histogram of the distribution of the position of a single molecule with respect to the other molecule or the tetramer, measured over the final 1 ns. The position is characterized by the collective variables that reflect the distance between the centers of the two entities and the coordinate of the monomer along the stacking direction (defined as the z-direction). The simulations for different force-fields gave different outcomes. For CHARMM Drude, two molecules do not create a stable dimer; however, the addition of one molecule to a tetramer results in preferential attachment of the monomer to the end of the stack. For GROMOS, we observed a stable dimer; however, a free molecule preferentially attaches to the side of a tetramer rather than to its end. GAFF forms neither dimer nor pentamer. In MartiniP, two molecules form a dimer; tetramer and free molecule seem to have a slight preference for the formation of a pentamer.

Figure 4

Simulation of monomer–monomer and monomer–tetramer pairs. (a–h) Histograms of positions in the last 1 ns of simulation of the added molecule to the system. The central molecule/tetramer shows a snapshot of the structure from simulations. Only the core of the molecules are shown, and side chains are shown semitransparent. (a–d) Histogram of distribution of two molecules around each other for different force-fields. In GROMOS and MartiniP molecules prefer to form dimers. For CHARMM Drude and GAFF (see main text), the molecules do not form dimers. (e–h) Histogram of distribution of addition of one molecule to the system with a tetramer (the tetramer is stabilized by position restraints on atoms of the cyclohexane rings). For CHARMM Drude and MartiniP there is a preference to attach to the end of the fiber (with a much stronger preference for CHARMM Drude). For GROMOS, the monomer tends to attach to the side of the fiber. For GAFF there is no preferential attachment. It is important to note that panels a–h show results for short simulation (10 ns) and, therefore, show local minima rather than the global one. (i) Free energy of attachment of monomers to the end of the fiber (orange) and to the side of the fiber (blue). Although the difference between energy levels is similar (∼20 kJ/mol), the energy levels for different force-fields vary substantially (up to 40 kJ/mol). There are no results for CHARMM Drude due to the lack of support of umbrella sampling for this force-field. (j) Mean square displacement of molecule added to tetramer. For all force-fields except GROMOS, the molecules stay mobile.

For long simulations that sample all states, histograms such as shown in Figure a–h would be equivalent to a free energy landscape. However, the limited simulation time (10 ns) could mean that these histograms show local minima rather than global ones. To validate them, we calculated the difference in free energy of adsorption of an unbound monomer to the side of the fiber and to the end of the fiber using umbrella sampling. We can approximate the experimental value of the free energy of elongation (which we interpret as an attachment to the end of the fiber) by the pseudophase approximation,[54] which results in ∼30 kJ/mol (see the SI). Due to the lack of support of umbrella sampling for CHARMM Drude, we do not have results for this force-field. All other force-fields prefer adsorption to the end of the fiber, and the difference between adsorption to the side and to the end of the fiber is similar for all force-fields, being at the level of ∼−20 kJ/mol. However, the difference in free energy between bound and unbound states is different for the different force fields. Surprisingly, GROMOS, GAFF, and MartiniP, which previously showed stable fibers and preferential elongation of tetramers, resulted in the free energies of elongation most different from the approximated experimental value. The largest difference is found for GROMOS, for which adsorption to the side and to the end is ∼−50 and ∼−75 kJ/mol, respectively. During the self-assembly process, the accessible surface area for newly arriving molecules is larger at the side of the tetramer than at its ends. Therefore, we can anticipate that molecules initially adsorb on the side of the fiber and then migrate to its end. However, for GROMOS the strong adsorption to the side of the fiber (50 kJ/mol) might prevent an adsorbed molecule from desorbing or from moving along the side of the fiber. Indeed, the 2D histograms presented before confirm that molecules adsorb preferentially on the side of the fiber and during 10 ns simulations rarely desorb (see Figure f). This can also be seen from the mean square displacement of the monomer over the last 1 ns of the tetramer-monomer simulations (see Figure j). From this calculation, it can be seen that molecules that adsorb to the fiber do not move anymore for the GROMOS model. For other force-fields, monomers still have some mobility.

Discussion and Conclusions

Self-assembly of supramolecular structures is a multistage process where details remain challenging to unravel with experimental and computational techniques. The molecular simulation results presented here show how different behavior can be observed with different force-fields and provide insights into the reasons behind the success or failure of MD simulations of self-assembly systems. In particular, we study the self-assembly and stability of CTA fibers using MD simulations with force-fields with different levels of detail, namely coarse-grained Martini and MartiniP, all-atomistic GROMOS and CGenFF, GAFF, and polarizable all-atomistic CHARMM Drude. In line with other research done in the field, the most challenging issue remains the time scale. Randomly dispersed molecules in solution need to diffuse to each other and explore their mutual orientations. This can lead to the formation of dimers, trimers, larger oligomers, and eventually larger aggregates. Often, the formation of a stable nucleus is a crucial step in the self-assembly: therefore, the propensity of a force-field for oligomers to stick together can strongly influence the required time scale for the formation of structure in a self-assembly simulation. For the CTA molecule, we found that the propensity to form dimers from monomers differs considerably between force-fields (Figure a–d), and we anticipate that these differences greatly influence required time scales for successful spontaneous assembly simulations.

The size of the nucleus may be determined by the range of collective or cooperative interactions. In this study, all force-fields show cooperativity in the stacking of CTA molecules (Figure a), where a convergence of interaction per molecule with size in vacuo, i.e., not accounting for solvent effects, is similar (6–8 molecules), but the strength of the cooperativity spans a range of about 70 kJ/mol per molecule. Nevertheless, the strength of the cooperative interaction in the CTA stack is not necessarily a good predictor for a stack being stable in a simulation starting from a perfectly stacked assembly (see Figure a–h): whereas the CGenFF Mod., Martini, and MartiniP force-fields have the largest value, the CGenFF Mod. and Martini stacks collapse but the MartiniP stack does not. The CGenFF shows only slightly lower cooperativity than the GROMOS force-field, but the former collapses, whereas the latter remains stacked. Interestingly, the driving force for self-assembly depends on the force-field. For Martini and MartiniP it is L-J interactions (although the small electrostatic contribution in MartiniP is essential for stable, ordered structure). For all-atomistic force-fields it is the combination of L-J and electrostatic interactions, with L-J stronger for nonpolarizable force-fields and electrostatic stronger for polarizable force-fields. The awareness of which interactions dominate in the self-assembled structure might be a crucial criterion of the choice of the force-field. However, often this choice is made by personal preference.

After a nucleus is formed, a larger structure evolves by monomers arriving at the nucleus and becoming part of the structure. Also, the initially formed structure may undergo internal reorganization, i.e., the molecules that arrive do not necessarily quickly take up the positions and orientations they have in the final assembly. This opens plenty of opportunities for kinetic traps in the self-assembly process. For the CTA molecules, this work shows that the likelihood of such a trap differs for different force-fields (Figure e–j). The different force-fields show that the free energy of the addition of an arriving monomer to the side of a tetramer stack is less favorable than an addition to the end (directly elongating the stack) by approximately the same amount. However, the binding energy at the side differs considerably, meaning that adsorption–desorption equilibria are strongly different (Figure i). The GROMOS force-field has such strong binding on the side of the fiber that effectively a molecule gets trapped there. The internal reorganization is likely to involve molecules adsorbed at the side of the fiber, given that longer fibers have a larger accessible area at their sides than at their ends. Thus, the ability of a molecule to either desorb from the fiber or to crawl along the fiber may be crucial for the formation of longer-ranged structure. The very low mobility of an adsorbed molecule on the tetramer for the GROMOS force-field (Figure j) may thus impede the success of self-assembly during MD simulation, even when the force-field reproduces the anticipated behavior of the final structure.

The diverse behavior of the systems studied here makes quantifying the quality of the force-field a challenging task. We believe that an efficient simulation of supramolecular fiber formation would require a force-field with two features: the ability to retain the stable fibrous structure and to elongate the existing stack. In this context, only CHARMM Drude and MartiniP were able to fulfill both conditions. Although GAFF can keep a fiber stable, it does not form oligomers. The latter is probably due to a limited time scale, as shown by the free energy of elongating and binding to the side of the stack; however, the slow formation of oligomers make simulations less efficient. The fiber in GROMOS is somewhat less stable, and during oligomerization, molecules interact strongly with the side of the stack, impeding the growth. The fibers in CGenFF and Martini do not retain stability and do not form oligomers. As demonstrated, CGenFF can be slightly improved when the mapping of the Drude charges is applied. We anticipate that further optimization of the force field, e.g., by using tools developed for this purpose, such as FFParam,[55] could further improve parametrization.

In summary, this work shows crucial aspects that have to be considered when simulating self-assembly systems. A priori knowledge about the final structure might be crucial for tuning the force-field. It is tempting to try to predict the long-range structure based on the preferred binding of two molecules, but it is to be expected that cooperative interaction plays an important role in supramolecular systems, and therefore an exploration of the (free) energy landscape of oligomers is expected to give important insights. For the CTA molecule, the study of dimerization provides little information about the expected success of self-assembly simulation. However, the study of the interactions of a tetramer stack with a monomer gives insight into possible kinetic trapping. Here, four tested force-fields gave four different outcomes of dimerization and elongation: upon formation of a dimer, we observed elongation and its lack (for MartiniP and GROMOS, respectively); upon lack of formation of a dimer, we also observed elongation and its lack (for CHARMM Drude and GAFF, respectively). Moreover, the driving force for self-assembly depends on the force-field. Taken together, the successful exploration of supramolecular assembly in the simulation requires awareness of the different stages of the self-assembly process and the (free) energy landscape relevant for each stage in the force-field of choice. The landscape should guide the choice of simulation technique, which can be as simple as straightforward MD in some stages but could require sophisticated adaptive or enhanced sampling techniques in others. An extensive study of the self-assembly processes of CTA with the CHARMM-Drude force-field will be the subject of a forthcoming paper.

Methods

Simulations were done with GROMACS, for nonpolarizable force-field version 5.1.2[56] and for polarizable CHARMM Drude modified version of GROMACS, allowing simulations using extended Lagrangian dynamics with a dual Nosé–Hoover thermostat.[43] For different simulations, we used different thermostats (Berendsen[57]/v-rescale[58]) and barostats (Berendsen[57]/Parrinello–Rahman[59]).

Details of Simulated Systems

Spontaneous Simulation:

Eight CTA molecules were placed in a 4.5 × 4.0 × 4.1 nm3 dodecahedron simulation box and solvated with water molecules (∼2000 water molecules). Fiber simulations: 24 CTA molecules were placed in elongated 4 × 4 × 12 nm3. Dimer: 2 CTA molecules were placed in a cubic 4.00 × 4.00 × 4.00 nm3 (∼2000 water molecules) or dodecahedron 4.00 × 4.00 × 2.83 nm3 (∼1400 water molecules) simulation box. Pentamer: 4 CTA stack was placed in a 5 × 5 × 5 nm3 simulation box. Additionally, single CTA molecules were randomly placed in the box and solvated with water molecules (2400 water molecules). All of the systems were energy minimized and equilibrated in NVT and NPT ensembles. Production runs were done in the NPT ensemble (except for CHARMM Drude, which was done in NVT). The simulations in the NPT ensemble were performed in isotropic conditions with an exception for simulations of fiber, which were done in semiisotropic conditions. More precise details of the systems are present in the SI, Tables S2–S5.

Parameters for GROMACS

Parameters of CTA molecules were obtained for GROMACS for different force-fields.

GROMOS

The parameters for GROMOS 53A6 force-field[36] were obtained using Automated Topology Builder.[60]

CGenFF

The parameters for CHARMM General Force-field (CGenFF) were obtained using cgenff_charmm2gmx.py script.

CHARMM Drude

The parameters for CHARMM Drude polarizable force-field were obtained on the basis of existing parametrizations of small molecules.[61−64] The final parameters for the CTA are included in the SI.

GAFF

The parameters for General Amber Force-field (GAFF) has been obtained using Antechamber[44] and charges calculated by AM1-BCC.[65]

Martini and MartiniP

The parameters for Martini and MartiniP were obtained according to the official parametrization tutorial available on the Martini FF Web site http://cgmartini.nl. See details in the SI. In MartiniP, the amide bead was treated similarly to the water bead: one bead with L-J potential, connected to two beads carrying charges. The charges were kept on opposite sides of the connected bead by distance constraints.

Hydrogen Bonds

We counted hydrogen bonds between amides groups using VMD and HBonds plugin. For all-atomistic force-fields we counted a hydrogen bond if the distance between hydrogen donor and acceptor was below 0.33 nm, and the angle of donor-hydrogen-acceptor was below 40°. For Martini FF, we counted the hydrogen bond if the donor–acceptor distance was below 0.4 nm and the donor-hydrogen-acceptor angle of 40°. For more details, see the SI.

SASA

Solvent accessible surface area (SASA) has been calculated gmx sasa with probe radius 0.14 nm for all-atomistic force-fields and 0.265 nm for Martini-based force-fields.