Conformational Entropy as Collective Variable for Proteins

TL;DR

This paper proposes using the NMR $S^2$ order parameter as a collective variable for enhanced sampling in protein folding simulations, enabling efficient folding/unfolding studies and free energy calculations.

Contribution

It introduces the novel application of the $S^2$ order parameter as a collective variable for enhanced sampling in protein simulations, linking entropy and folding dynamics.

Findings

Successfully reversible folding and unfolding of a small protein using $S^2$ with Metadynamics.

Efficient computation of free energy and unfolding rates at various temperatures.

Demonstrates the relation between $S^2$, entropy, and protein folding energetics.

Abstract

Many enhanced sampling methods, such as Umbrella Sampling, Metadynamics or Variationally Enhanced Sampling, rely on the identification of appropriate collective variables. For proteins, even small ones, finding appropriate collective variables has proven challenging. Here we suggest that the NMR $S^2$ order parameter can be used to this effect. We trace the validity of this statement to the suggested relation between $S^2$ and entropy. Using the $S^2$ order parameter and a surrogate for the protein enthalpy in conjunction with Metadynamics or Variationally Enhanced Sampling we are able to reversibly fold and unfold a small protein and draw its free energy at a fraction of the time that is needed in unbiased simulations. From a more conceptual point of view this implies describing folding as a resulting from a trade off between entropy and enthalpy. We also use $S^2$ in combination with the free energy flooding method to compute the unfolding rate of this peptide. We repeat this calculation at different temperatures to obtain the unfolding activation energy.

Equations

$$S^2=\sum _n{S}_n^2$$

$$S_n^2=\tanh \left(0.8\sum _k[\exp (-r_{n-1,k}^O)+\exp (-r_{n,k}^{H_N})]\right)-0.1,$$

Full text

Enhanced sampling method have received great attention since they offer the promise of overcoming the limited time scale that direct simulations can afford. Among the plethora of methods proposed in the literature Voter (1997); Huber et al. (1994); Grubmüller (1995); Torrie and Valleau (1977), Metadynamics (MetaD) Laio and Parrinello (2002); Barducci et al. (2008) and more recently Variationally Enhanced Sampling (VES) Valsson and Parrinello (2014) are finding application in a vast array of problems Palazzesi et al. (2013); Martoňák et al. (2006); Pavan et al. (2013); Branduardi et al. (2007); Tiwary et al. (2015); Bottaro et al. (2016); Shaffer et al. (2016). Like other similar approaches they rely on the identification of appropriate collective variables (CVs) or order parameters. In MetaD or VES the fluctuations of the CVs are enhanced such that transition between different metastable basins are favored Barducci et al. (2011); Valsson et al. (2016).

A vast number of CVs have been suggested and are easily accessible in open source codes Tribello et al. (2014). However the quest for new order parameters continues in the hope of finding CVs that are efficient and yet generic enough such that they can be applied to a larger class of problems, without prejudging the final results. Identifying appropriate CVs is not only a technical issue needed to accelerate sampling but offers a key to the understanding of the underlying physical processes. The need for such CVs is particularly pressing in the field of bio-molecular systems, such as proteins, whose conformational changes are defined by a large number of degrees of freedom like the arrangement of backbone atoms, side-chains and solvent molecules.

At first this appears like a very demanding request. The purpose of this paper is to show that this is not necessarily so, at least for small proteins or selected regions of larger bio-systems. To this effect we introduce a conceptually new CV and illustrate its efficiency in non trivial examples. Our guiding principle is that in the behavior of proteins, and of many other systems Piaggi et al. (2016), entropy plays an important role and if we are able to identify a CV that measures entropy even if in an approximate way, this could go some way towards being able to sample the complex landscape of proteins Chavez et al. (2004).

In this search we shall be helped by the NMR literature in which several attempts have been made at converting the NMR observable into a measure of conformational entropy Akke et al. (1993); Yang et al. (1997); Stone (2001); Sharp et al. (2015); Yang and Kay (1996); Li et al. (1996); Marlow et al. (2010). Without going into the complex detail of the NMR technique, it suffices to say that in nuclear relaxation experiments, it is possible to measure the dynamical behavior of selected bonds, like N-H or C-H, and, within some approximations, extract the so-called order parameter $S^{2}$. This parameter, that can vary between 0 and 1, can provide useful information on the degree of spatial motion of the system Lipari and Szabo (1982a, b). From the knowledge of this parameter several empirical relation between $S^{2}$ and the conformational entropy have been proposed and their validity assessed in comparison with either experiments or molecular dynamics simulations Akke et al. (1993); Li et al. (1996); Yang and Kay (1996); Sharp et al. (2015). Despite being rather empirical Stone (2001), these relationships have been used to study several bio-macromolecular processes, and to calculate protein heat capacity Spyracopoulos and Sykes (2001); Yang et al. (1997); Marlow et al. (2010); Bracken (2001).

It is the existence of these relations and the notion that entropy plays an important role in proteins that have inspired us to use $S^{2}$ as CV. However, to order to proceed one needs to express $S^{2}$ as a function of the atomic coordinates. Luckily such relations are available for both the NH and CH3 groups Bowman (2016); Chatfield et al. (1998); Peter et al. (2001); Maragakis et al. (2008). Since here we shall only bias the N-H bonds related order parameter, we solely report the relevant expression proposed by Zhang and Brüschweiler Zhang and Brüschweiler (2002):

[TABLE]

where $S^{2}_{n}$ is the order parameter for the $n$-th amino acid residue defined as:

[TABLE]

where $k$ runs over all the heavy atoms with the exception of the residues $n$ and $n-1$ and $r^{H_{N}}_{n,k}$ and $r^{O}_{n-1,k}$ are the distances from the heavy atom $k$ from the amide hydrogen in residue $n$ and the carboxyl oxygen in residue $n-1$. In Eq. (2) the distances are to be expressed in Ångstroms.

In this paper we assume that this relation is valid for any atomic configuration and we use it as a CV. From a practical point of view the merit of this choice is in the end to be judged by the results. Once the CV is identified, the fluctuations of $S^{2}$ can be amplified by using an enhanced sampling technique, such as MetaD or VES.

Before starting our calculations we indeed checked that $S^{2}$ is able to distinguish between the folded and unfolded protein conformations. For this reason we calculated the free energy surface (FES) along this CV utilizing the chignolin (CLN025) Honda et al. (2008) trajectory provided to us by the D.E. Shaw Research Group Lindorff-Larsen et al. (2011). As can be seen in Fig. 1, the FES exhibits two well-defined basins corresponding to the folded and unfolded state. This in itself is an interesting result that suggests the attempt at using $S^{2}$ as collective variable is not totally devoid of merits.

We tried to drive the folding transition using only $S^{2}$ (hereafter denoted $s_{S}$) as CV but failed. Thus in the spirit of this work in which we describe folding as a trade off between entropy and enthalpy we introduced a second CV that is meant to be a surrogate for enthalpy. This could have been accomplished by separating from the total energy of the systems those components that describe the protein-protein interactions and use these as CV. However this is somewhat expensive and we preferred to use as surrogate for the enthalpy the native H-bonds contact map (hereafter $s_{H}$). At the end we shall reweight the FES Tiwary and Parrinello (2014) and obtain its projection onto the protein enthalpy defined as the sum of the protein-protein contributions, $E_{pp}$, to the total energy.

We performed the simulation at $T=340$ K as in Ref. 37 and we use the same potential, the Charmm22* Piana et al. (2011) protein force field and the TIP3P Jorgensen et al. (1983) water model. We use GROMACS 5.1.4 MD package Abraham et al. (2015) patched with the PLUMED 2 plug-in Tribello et al. (2014) in which we have implemented the $s_{S}$ CV. In the MetaD calculations we integrate the equation of motion using a time step of 2.0 fs and the temperature is controlled by the stochastic rescaling thermostat Bussi et al. (2007). Gaussians of initial height 2.82 kJ/mol and width 0.05 for $s_{S}$ and for $s_{H}$ were deposited every picosecond. The value of the bias factor $\gamma$ was set equal to 8 Barducci et al. (2008). To speed up the calculation and make use of parallelism we used 4 multiple walkers Raiteri et al. (2006). We evaluate the convergence using the error metric previously used in Ref. Barducci et al. (2008); Branduardi et al. (2012); Valsson and Parrinello (2015) and using the unbiased data of Ref. 37 as reference.

It can be seen from Fig. 2 that the convergence in the FES is reached in about 1.0 $\mu$s, that is a much shorter time relative to $\sim$100.0 $\mu$s reported in Ref. 37. This reflects the fact that in the MetaD run the rate of transition between folded and unfolded states is accelerated about one hundred times.

In this plot we also express the FES as a function of entropy and enthalpy. The entropy is extimated from sS using the relation given by Wand and co-workers Sharp et al. (2015), while as a measure of the protein enthalpy we use Epp. This calculation clearly proves the usefulness of using entropy and enthalpy to drive reversible transitions between folded and unfolded states. We have repeated the calculations using VES obtaining statistically indistinguishable results. The VES simulations are performed using the VES code ves module for PLUMED 2. The results for this calculation are reported in the SM.

An analysis of the structures lying at the bottom of the folded basin in Fig. 1 shows that they deviate from the experimental structure Honda et al. (2008) on average by around 0.8 Å. This suggests that if we are only interested in studying the unfolding process, $s_{S}$ alone could be used to promote unfolding transition and calculate unfolding rates using the approach introduced in Ref. 47. This is the VES analogous of the infrequent metadynamics method of Tiwary and Parrinello Tiwary and Parrinello (2013), that in turn is based on the ideas of Grubmuller Grubmüller (1995) and Voter Voter (1997). In all these approaches one relates the rates of a rare event, as calculated in a biased system, to the physical unbiased rates by a simple relation. The requirement that the bias does not act in the transition region is crucial for this relation to hold. The methods mentioned earlier differ in the way this is achieved. In Ref. 47 the variational approach is used to truncate the bias up to a preassigned free energy level. If this cutoff value is smaller than the free energy barrier this latter region remains free of bias and the condition put forward by Grubmuller Grubmüller (1995) and Voter Voter (1997) applies. Computational details can be found in the SM. The results of our calculation are shown in Fig. 3. The value obtained at $T=340$ K is in good agreement with both the unbiased estimation of Lindorff-Larsen et al. Lindorff-Larsen et al. (2011) and the infrequent MetaD simulations performed by Tung and Pfaendtner on a mutated form of our system, using the root-mean square deviation with respect to the folded structure as biased CV Tung and Pfaendtner (2016).

Contrary to these previous estimations we push our calculation to lower temperatures (320 and 300 K), finding for the unfolding of this simple peptide an Arrhenius behavior, with an activation energy of around 50 kJ/mol. This value is in the right ballpark when compared to the estimation based on a similar $\beta$-peptide Scian et al. (2013).

In conclusion we have shown that the idea of using entropy and enthalpy as collective variables, or rather surrogate expressions for them, is a useful one. One of the secret of the success in using these CVs is not only that it is founded on physical ideas and concepts but also on the fact that is non-local and thus sensitive to the whole structure of the protein. We would like to add that the $S^{2}$ defined in Eq. 2 reflects the backbone structure. If one were interested in the role of the side-chains the use of $S^{2}$ order parameter based on the methyl groups might prove useful.

We believe that here we have introduced a powerful new concept in the simulation of protein. How far it can be pushed when going to larger systems it remains to be seen. Already at this stage, it can be safely said that, without any modification, the folding of small proteins, the study of intrinsically disordered proteins and the conformational flexibility of proteins segments can be handled as described above Kokh et al. (2016). We are confident that with an appropriate adaptation much larger proteins can be similarly handled.