Metadynamics is a powerful computational tool to obtain the free energy landscape of complex systems. The Monte Carlo algorithm has proven useful to calculate thermodynamic quantities associated with simplified models of proteins, and thus to gain an ever-increasing understanding on the general principles underlying the mechanism of protein folding. We show that it is possible to couple metadynamics and Monte Carlo algorithms to obtain the free energy of model proteins in a way which is computationally very economical.
Metadynamics is an algorithm which coupled to molecular dynamics provides an efficient tool to obtain the energy landscape of systems displaying large energy barriers, and thus whose sampling by standard tools is, at best, problematic. It is based on the knowledge of few slow collective variables s i of the system and on the use of a non-Markovian potential U (s i ) that disfavors the exploration of regions of the phase space already visited by the system (Laio and Parrinello (2002)). This algorithm has been succesfully used to obtain the free energy of molecular systems at atomic detail (Babin et al. (2006)).
In the case of simplified protein models, where the atomic structure of each amino acid is coarse-grained, it is common to sample the conformational space with the help of Monte Carlo algorithms. Such an approach is computationally more economic and more simple to implement than the corresponding molecular dynamics algorithm (see, e.g. Shimada et al. (2001), Kussell et al. (2002), Shimada and Shakhnovich (2002)). It is then natural to try to extend metadynamics so as to make it possible to couple it to a Monte Carlo algorithm.
Of course, other modifications of the straight Monte Carlo sampling have been developed during the last tens of years, including simulated tempering, multicanonical sampling, parallel tempering, etc. All of them are aimed at preventing the system to get trapped in free energy minima. In the following we show that Monte Carlo metadynamics is efficient, accurate and particularly easy to implement.
We apply a scheme to the calculation of the free energy, as a function of the RMSD, of a small domain protein, namely Src-SH3. It is a widely studied domain (Grantcharova et al. (1998), Yi et al. (1998), Riddle et al. (1999)) of the Proto-oncogene tyrosine-protein kinase Src, a 536 residue protein that plays a multitude of roles in cell signalling. Src is involved in the control of many functions, including cell adhesion, growth, movement and differentiation. SH3 is a domain built out of 60 residues, displaying mainly β-strands (see Fig. 1). From calorimetry and fluorescence experiments, it is known to fold according to a two-state mechanism, that is, populating at biological temperature mainly two states (the native and the unfolded state) (Grantcharova and Baker (1997)). Consequently, we expect the free energy landscape to display two minima separated by a barrier.
In Section 2 we present the protein model used in the simulations along with a working description of the algorithm. We devise a formal proof of the correctness of the method in Section 3 and then test it in the specific case of Src-SH3 in Section 4.
The model employed in the simulations describes the protein as a chain of beads centered on the C α of the protein backbone (see Fig. 2). The allowed moves are the flip-move Each Monte Carlo step, we apply a Metropolis algorithm (Metropolis et al. (1953)) where the transition probability is given by
that is, the probability with which the next Monte Carlo move is accepted is calculated on the variation of the energy of the system, plus the variation of the metadynamics potential.
During each fragment of trajectory after the update of the non-Markovian potential at each time T , the collective variable s explores a region A(T ). If one makes the critical assumption that the dynamics has been able to visit this region so extensively that ergodicity holds, then the probability distribution of the collective variable is
.
(2)
After the end of this sampling, the non-Markovian potential is updated, and the new potential reads
where W is the heigth of the energy added to the non-Markovian potential. Further assuming that W is small, that is that the new term does not perturb in an important way the shape of the potential U , then the previous equation can be rewritten as
.
(4)
Once the free energy landscape is completely filled by the non-Markovian energy, then the growth of this non-Markovian energy will be independent on s, that is
where A is the whole interval spanned by the collective variable. Integrating by saddle-point evaluation, leads to
where s 0 is defined by U ′ (s 0 ) = -F ′ (s 0 ). The last equation states that the free energy of the system is, except for an additive constant, equal to the opposite of the non-Markovian potential. A nice property of this algorithm is that the obtained free energy depends logarithmically on any additive error in the determination of U (s, T + τ ) (i.e., if one adds ǫ(s)τ to Eq. ( 4), one obtains an additive term log ǫ(s) in F (s)).
In order to obtain a reference free energy landscape as a function of the RMSD for comparison to the metadynamics reconstructed landscapes, we first carried out a fairly long classical Monte Carlo simulation (90 billions of Monte Carlo Steps (MCS)). The free energy calculated at temperature θ = 0.625 (slightly below the folding temperature, defined as the temperature at which the volume of the na
This content is AI-processed based on open access ArXiv data.