High resolution ensemble description of metamorphic and intrinsically disordered proteins using an efficient hybrid parallel tempering scheme

Mapping free energy landscapes of complex multi-funneled metamorphic proteins and weakly-funneled intrinsically disordered proteins (IDPs) remains challenging. While rare-event sampling molecular dynamics simulations can be useful, they often need to either impose restraints or reweigh the generated data to match experiments. Here, we present a parallel-tempering method that takes advantage of accelerated water dynamics and allows efficient and accurate conformational sampling across a wide variety of proteins. We demonstrate the improved sampling efficiency by benchmarking against standard model systems such as alanine di-peptide, TRP-cage and β-hairpin. The method successfully scales to large metamorphic proteins such as RFA-H and to highly disordered IDPs such as Histatin-5. Across the diverse proteins, the calculated ensemble averages match well with the NMR, SAXS and other biophysical experiments without the need to reweigh. By allowing accurate sampling across different landscapes, the method opens doors for sampling free energy landscape of complex uncharted proteins.


Supplementary
* In case of Beta-hairpin the REHT simulation is compared with the REST2 simulation performed earlier (Ref. 1 ) using the same forcefield. # Transit time indicates the time taken for a replica to jump over the complete state space. The complete replica mixing allows the full utilization of the temperature space and provides sufficient energy for a biomolecule to undergo conformational transition across the barrier. However, in complex systems like RFA-H, the REST2 method is shown to be ineffective where the complete mixing of replicas is not achieved.  Figure 6: Convergence of sampling in physiologically relevant base replica of Trpcage with REST2 (red) and REHT (blue) methods. The CPU time including the contributions from all replicas is provided in additional X-axes. The figure clearly indicates that the results obtained from REHT converges faster and better than the REST2 simulations.

Supplementary Note 1: A quantitative comparison of ergodicity with REST2 and REHT simulations
We assessed the convergence of REST2 and REHT at the physiologically relevant base replicas (@Temp 300K) similar to that suggested by Dave Thirumalai's group 2 and Bruce Berne's group 3 . Towards this, the ensemble at the base replica is split into two equal parts, one at the start (part A) and other at the end (part B) of the total simulation (leaving a gap of 50ns between them). The two parts of the simulation represents the two independent trajectories started from different conformations of protein. The ergodicity is then measured by comparing the two simulations as follows: The conformational landscape represented by the two CVs is discretized into m*n uniform bins. The population in each of the (i,j) th bin ( !,# ) is compared between the two parts of the simulations A and B and the overall difference is measured using $ parameter defined as: Figures 6 and 14 depict the $ as the increasing length of simulation time. The figure indicates that the REHT converges faster than the REST2 in both Trp-cage and His-5 simulations. The difference in convergence between the two methods is more accentuated in His-5 simulation.
A more rigorous measure of convergence in replica exchange simulations is to check the convergence of distributions sampled in the independent time continuous replicas. This has been shown earlier with alanine-di-peptide system 4 , where the energy landscape is much smoother. This analysis is extremely challenging on very rugged realistic protein landscapes and demanding from sampling point of view. We performed these analyses on the time-continuous replicas of Trp-Cage (Supplementary Figure 7). The distributions of RMSD to the native structure reveal similar trends for replicas that did not explore the folded state, whereas the replicas that fold, sample distributions that are different than each other. This is observed for both REST2 and REHT methods. We anticipate that because of the complex rugged nature of the energy landscape and the existence of multiple pathways to folding. Achieving convergence of protein folding by sampling all the kinetic pathways in a single replica is a challenging task with the current level of simulation time in replica exchange simulations. 5 We would also like to emphasize that in fact the current state-of-the-art REX methods such as REST2 and gREST explored independent folding in one or two replicas leaving most of the replicas predominantly unfolded. 1 In that respect, the REHT definitely does a far better job by exploring the folding independently in six different replicas. We assume extending the simulation longer or increasing the temperature of water further may help achieving this.
Notably, the converged sampling of a flatter landscape in intrinsically disordered His-5 is easily accessible with the REHT method compared to REST2, as will be discussed in Supplementary

Supplementary Note 2: Dihedral switch in Alanine dipeptide
Being the simplest biomolecular model that possesses multiple conformational basins, alanine dipeptide serves as one of the first choices in enhanced sampling literature to study the barrier crossing events. The peptide exhibits three basins across the phi and psi dihedral space. We employed both the conventional REST2 and the novel REHT methods to sample the switching between these dihedral basins of the peptide. The sampled basins across the phi and psi space at different time intervals are illustrated for the lowest temperature replica in Supplementary Figure  4. While REST2 can sample all the basins within 4 ns timescale, REHT required 6 ns to completely sample this dihedral space. However, the relatively slow transitioning phi angle is frequently sampled in the latter method than the former method (Supplementary Figure 5). Moreover, the transition was sampled by 4 out of 5 replicas in REHT method, but only once by a single replica out of three in conventional REST2 method. This result suggests that optimal heating of both solute and solvent allows effective crossing of the energy barrier and facilitates frequent transition in the conformational space.
Supplementary Estimation of the relative computational cost between REHT and REST2 should be performed in fully converged simulations. However, the results, as plotted on Supplementary Figures 6 and 14 indicate that the REST2 is not well converged as that of REHT in the timescale simulated for both Trp-cage and His-5 simulations.
We are uncertain how the trend flares out upon extending the REST2 simulation for estimating the relative computational cost. However, extrapolation of the existing data, which assumes that the trends will continue, suggests the REST2 would take about 1000ns/replica in order for attaining the convergence level of 0.1 in TRP-cage (Supplementary Figure 16). On the contrary REHT achieves this in 500ns. In terms of computational cost (CPU units), this is only 1.3 folds reduction in REHT in comparison to REST2 considering all the replicas. Surprisingly, in case of the IDP, His-5, REHT is about 12 folds faster than REST2 (2600ns in REST2 vs 150ns in REHT for achieving a $ value of 0.02 and considering a total of 10 vs 15 replicas used).
Supplementary Figure 16: Estimating relative computational cost from the convergence decay for a) Trp-cage and b) His-5 with REST2 (red) and REHT (blue). The best fit lines for each of the curves have been plotted (green and magenta) and their respective functions and parameters are indicated. The CPU time including the contributions from all replicas is provided in additional Xaxes. a

Supplementary Note 4: Quantification of Heterogeneity of ensemble
To quantify the heterogeneity as suggested by Lyle et.al, 6 each of the conformation in an ensemble is represented as a vector of inter-atomic distances of all C-alpha atoms ( ). The distance, , between two conformational vectors ( , -) is then computed with the cosine distance defined as ,-= 1 − , . -/| , || -|. The larger D value indicates more conformational heterogeneity and vice versa. For an ensemble of N conformations, the pairwise distance calculations yield N(N-1)/2 distance values. The distribution of D calculated for His-5 ensemble with REST2 and REHT is plotted in Supplementary Figure 17a. As a control, we also generated the Flory Random Coil (FRC) ensemble for the specified sequence of His-5 using Campari tool and plotted their distance distribution. 5000 conformations were generated with 3*10 7 Monte Carlo steps after discarding first 500000 steps of equilibration. Our results indicate larger heterogeneity for the REST2 ensemble than for the REHT and FRC, which at face value is quite counterintuitive.
To verify this result, we also used a second heterogeneity measure described by Papoian et.al. 7 According to this metric, the heterogeneity (q) between two conformational vector k and l is computed with the differences in pair distance values ( .,/ ) as follows: , where the is a resolution parameter and usually is set to 2Å.
Consistent with the cosine distance distribution (Supplementary Figure 17a), the histogram of the pairwise-q (Supplementary Figure 17b) also reveals larger heterogeneity (low q) for REST2 than the REHT and FRC. However, the free energy landscape of His-5 shown with Rg and SASA (Supplementary Figure 12) suggests a confined sampling of compact states in REST2. These conflicting results needed to be reconciled and we anticipated that one of the possible explanations for this inconsistency could have its origin in the "heterogeneous compact structures" that the REST2 simulations sample for the His-5 ensemble.
To inspect this, we transformed the cosine distances (Dkl) between all pairs of conformations (n*(n-1)/2 values) into a 2-dimensional map of n conformations using multidimensional scaling (MDS). MDS is particularly useful for visualizing the distance matrix in low dimensional space while preserving the between-object distances as much as possible. Supplementary

Supplementary Note 5:
Replica exchange simulation methodology Replica exchange simulations allows for studying equilibrium properties with lesser computational time by its rapid relaxation and improved conformational sampling. Originally developed in Monte-Carlo background (Swendsen and Wang, 1986), 8 the method was later introduced to Molecular dynamics simulations by Sugita and Okomoto in 1999 9 (Temperature replica exchange molecular dynamics (TREM)). The method achieves effective sampling by simulating a series of low and high temperature replicas, while allowing the exchange of configurations at regular intervals. Moreover, due to its stochastic nature of exchanges that ensures the detailed balance, it generates Boltzmann weighted ensemble from which it is straight forward to obtain the thermodynamic averages. The probability of accepting the exchange between replica m and n depends on the difference in the Boltzmann weight factor, that is exponentially related to the difference in energy and temperature (Equation (S1)).
Where ( and ' are reciprocal temperatures (1/( & ( ) and 1/( & ' ) respectively). : ' and : ( are potential energy of replicas m and n respectively. For a larger system the temperature differences between two replicas should be chosen minimal to yield viable exchange acceptance ratio. This demands large number of replicas to cover a sufficient range of temperatures. Replica exchange solute tempering (REST2): Instead of changing the temperature across the replica ladder the advanced replica exchange solute tempering (REST2) method scales the energy function in a particle-wise manner such that the solute is effectively heated up while keeping the water cold. For instance, the potential energy function of replica m is broken down into intramolecular protein interactions ( ;; ), protein-water interactions ( ;< ), and water self-interactions ( << ) whose potentials are scaled as shown in Equation (2).
Where λ = 0 , in which A and ' are reciprocal temperatures of 0 th and m th replicas. While running all the replicas at the same temperature, the method cancels out the energy difference in water self-interaction energy that otherwise hugely contribute for the poor scaling as in TREM. Thus, the acceptance probability of REST2 as shown in Equation (S3) depends only on the energy differences of intramolecular solute energy and intermolecular energy between protein and water. ' , ( and A are reciprocal temperatures of replicas m, n and 0 th replicas respectively. Replica exchange hybrid tempering (REHT): Though the REST2 claims to be efficient for studying larger proteins including that of IDPs as shown recently by Shrestha et.al, it suffers from inability to explore the complex energy landscape with larger energy barriers. We speculated that this could be due to the imbalance between hot solute and cold solvent that causes differential dynamics of central protein and surrounding bulk water. In general, the biomolecular folding and conformational transition is tightly coupled to its surrounding water dynamics. Hence in this work, we introduce a hybrid method that optimally treats the protein as well as the surrounding water. We achieve this by associating the replicas to different bath temperatures in addition to scaling down the potential function of protein in contrary to REST2. Treating the replicas in such a combination doesn't violate the detailed balance condition. More importantly the method yields expedited protein conformational sampling by allowing efficient rewiring of water. For the viable exchanges the temperature gaps between the adjacent replicas are kept minimal. At the same time the protein is allowed to effectively heatedup to a larger extent by additional REST2 scaling factor.