Enhanced sampling method with coarse graining of conformational space

The sampling of conformations in the molecular simulations for systems with complicated free energy landscapes is always difficult. Here, we report a method for enhanced sampling based on the coarse-graining of conformational space. In this method, the locally converged region of the conformational space is coarse-grained with its population characterized by the related average residence time and visiting number, and at the same time, the direct simulations inside it are eliminated. The detailed balance is satisfied by updating the visiting number and generating outgoing trajectories of this region. This kind of coarse-graining operation can be further carried out by merging all the neighboring regions which are already converged together. The global equilibrium is achieved when the local equilibrated regions cover all the interested areas of the landscape. We tested the method by applying it to two model potentials and one protein system with multiple-basin energy landscapes. The sampling efficiency is found to be enhanced by more than three orders of magnitude compared to conventional molecular simulations, and are comparable with other widely used enhanced sampling methods. In addition, the kinetic information can also be well captured. All these results demonstrate that our method can help to solve the sampling problems efficiently and precisely without applying high temperatures or biasing potentials.

Figure 1

Schematic diagram of the CGCS simulations. (1) The conformational space is divided into small regions by gridding along the predefined reaction coordinates. The locally equilibrated region (LER) by regular dynamic simulations is designated as a coarse-grained unit, in which further molecular simulations are eliminated. Consequently, the regular dynamic simulations only focus on the nonequilibrated regions (NERs) and the transitions in and out of the LER (blue arrow lines). (2) The adjacent LERs in conformational space can be further merged into a higher level LER (hLER). If the hLER is large enough, a barrier crossing event can be sampled, which leads to the emergence of a new LER (3). (4) The global equilibrium is achieved when all the interested areas of the landscape are covered by LERs, and sufficient number of barrier crossing events are sampled. A typical trajectory segment divided by the boundary of the region is also shown here, and the sampling in each region can be seen as the summation of the trajectory segments.

Figure 2

The free energy landscapes of the double-basin system based on (a) the reference direct simulation and (b) the simulation with the CGCS. Here in the CGCS simulation, the threshold $n_{\text{th}}$ is set as 50. The unit of the free energies is $k_{B}T$. In the direct simulation, $2\times {10}^{11}$ steps were involved in order to sample 300 barrier crossing events. Whereas in the corresponding simulation with the CGCS, only $2\times {10}^6$ steps were needed to sample the same number of barrier crossing events.

Figure 3

The snapshots of the LERs during the simulation with the CGCS, at the simulation time of (a) $4.8\times {10}^5$, (b) $1.1\times {10}^6$, (c) $1.3\times {10}^6$, and (d) $2.2\times {10}^6$ steps.

Figure 4

The correlation plot between the residence time of the growing hLER calculated based on the CGCS simulation (${\tau }_{\mathrm{CGCS}}$) and the reference direct simulation (${\tau }_{\mathrm{Ref}}$). The perfect correlation suggests that the CGCS simulation can well reproduce the residence time of the hLERs with different sizes. Here the simulation starts from the left basin. When the hLER grows to the one labeled by the arrow line, the hLERs initiated from the two energy basins of the model potential get merged, and the global convergence is achieved.

Figure 5

The computing cost for achieving a certain number of barrier crossing events $N_{\mathrm{BC}}$ based on the reference direct simulations (Ref) and the simulations with the CGCS. The efficiency enhancement ($R_E$) is also shown. In our simulations, the global convergence (the hLERs from two energy basins get merged) is achieved when the $N_{\mathrm{BC}}$ exceeds 318.

Figure 6

The free energy landscapes of the four-basin model potential based on (a) the reference direct simulation and (b) the simulation with the CGCS. All the parameters are the same as those used in the double-basin model. The unit of the free energies is $k_{B}T$. In the direct simulation, $3.1\times {10}^{11}$ steps were involved in order to sample 50 barrier crossing events over the highest energy barrier. Whereas in the corresponding simulation with the CGCS, only $1.4\times {10}^6$ steps were needed to sample the same number of barrier crossing events.

Figure 7

(a),(b) The cartoon representations of the three-dimensional structures of the closed (PDB: 1cfd) and open (PDB: 1cll) state of the C-terminal domain of calmodulin [52, 53]. (c),(d) The 2D free energy landscapes along the reaction coordinates $RMS{D}_{\mathrm{open}}$ and $RMS{D}_{\mathrm{closed}}$ based on (c) the reference direct simulations and (d) the simulations with the CGCS. The reaction coordinates $RMS{D}_{\mathrm{open}}$ and $RMS{D}_{\mathrm{closed}}$ represent the root mean square deviations of a given structure from the native structures at the open and closed state, respectively. O, open state; C, closed state; PD, partially disordered state. Here in the CGCS simulations, the threshold $n_{\text{th}}$ is set as 20. The unit of the free energies is $k_{B}T$. The results were calculated based on 400 independent simulations. In each direct simulation, $5\times {10}^9$ MD steps were involved in order to sample four barrier crossing events. Whereas in each corresponding simulation with the CGCS, only $2\times {10}^6$ MD steps were needed to sample the same number of barrier crossing events.

Figure 8

(a) The standard deviation of the free energy differences from 100 independent metadynamics simulations as a function of the parameter $w$. The standard deviation is represented by the averaged value ($\overline{{\sigma }_{\mathrm{\Delta }F}}$) among all the seven free energy differences listed in Table S2 [46]. The dashed line represents the corresponding value of the $\overline{{\sigma }_{\mathrm{\Delta }F}}$ for the CGCS simulations with $n_{\mathrm{th}}=50$ (0.46). The metadynamics simulations with $w=0.0054$ give comparable accuracy as the CGCS simulations with $n_{\mathrm{th}}=50$. (b) The $\overline{{\sigma }_{\mathrm{\Delta }F}}$ (red) and averaged free energy barrier height for the transition from state C to state D ($\mathrm{\Delta }{F}_{C\to D}^{†}$) (blue) from 100 independent metadynamics simulations as a function of simulation time with the parameter $w=0.0054$. The dashed lines represent the converged values of the standard deviation (red) and free energy barrier height (blue). When the simulation time is larger than $1.2\times {10}^6$, the standard deviation and free energy barrier height get converged.