Parallel Tempering of Dark Matter from the Ebola ... - ACS Publications

Subscriber access provided by UNIVERSITY OF ADELAIDE LIBRARIES

Article

Parallel Tempering of Dark Matter from the Ebola Virus Proteome: Comparison of CHARMM36m and CHARMM22 Force Fields with Implicit Solvent Mark A. Olson J. Chem. Inf. Model., Just Accepted Manuscript • DOI: 10.1021/acs.jcim.7b00517 • Publication Date (Web): 29 Nov 2017 Downloaded from http://pubs.acs.org on December 1, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Information and Modeling is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 24 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling


Parallel Tempering of Dark Matter from the Ebola Virus Proteome: Comparison of CHARMM36m and CHARMM22 Force Fields with Implicit Solvent Mark A. Olson Molecular and Translational Sciences, USAMRIID, Frederick, Maryland, 21702 ABSTRACT: Intrinsically disordered proteins are characterized by their large manifold of thermally accessible conformations and their related statistical weights making them an interesting target of simulation studies.

To assess the development of a computational

framework for modeling this distinct class of proteins, this work examines temperature-based replica exchange simulations to generate a conformational ensemble of a 28-residue peptide from the Ebola virus protein VP35. Starting from a prefolded helix-β-turn-helix topology observed in a crystallographic assembly, the simulation strategy tested is the recently refined CHARMM36m force field combined with a generalized Born solvent model. A comparison is provided of two replica exchange methods where one is a traditional approach of a fixed set of temperatures and the other is an adaptive scheme of allowing the thermal windows to move in temperature space. The assessment is further extended to include a comparison to equivalent CHARMM22 simulation datasets. The analysis finds CHARMM36m to shift the minimum in the potential of mean force (PMF) to a lower fractional helicity as compared to CHARMM22, while the later showed greater conformational plasticity along the helix-forming reaction coordinate. Among the simulation models, only the adaptive tempering method with CHARMM36m found an ensemble of conformational heterogeneity consisting of transitions between α-helix-β-hairpin folds and unstructured states that produced a PMF of fractional fold propensity in qualitative agreement with circular dichroism experiments of reporting a disordered peptide.

1 ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 24


1. INTRODUCTION Intrinsically disordered proteins (IDPs) that encompass “Dark Matter” proteomes play a fundamental role in the regulation and function of protein association networks.1-4

Their

hallmark is large-scale conformational heterogeneity in free solution while finding a folded topology upon forming a multimeric assembly. An illustrative example is a 28-residue peptide region extracted from the Ebola virus protein VP35.5 The X-ray crystallographic structure of the peptide (designated as NPBP) bound to the Ebola virus protein NP shows a helix-β-turn-helix topology. In free solution NPBP transitions to a disordered ensemble as observed from circular dichroism (CD) spectroscopy.5 What makes the disorder-order transition of NPBP of larger interest is when added to a solution of 50% trifluoroethanol (TFE) the CD spectrum shows the peptide transitions from an unstructured ensemble to structures containing helical folds. The hidden propensity of NPBP and conceivably that of many other IDPs makes atomistic simulation studies challenging to capture a heterogeneous conformational ensemble without being overly biased by the susceptibility to fold. The challenge is amplified by the inherent deficiencies of all-atom force fields and solvent models that tend to over stabilize fold propensities. The topic of this work is to test the recently reported CHARMM36m force-field parameterization6 in a simulation of NPBP using a temperature-based replica-exchange (ReX) method.7 The force field is a refinement of an earlier version to improve the accuracy in polypeptide backbone conformational ensembles for IDPs.

Although the development of

CHARMM36m is primarily intended for explicit solvent simulations, the work here applies the force field with a generalized Born (GB) solvent approximation represented by the GBMV2 model.8 GBMV2 is one of the most accurate GB models9 and reproduces the thermal stability of small α-helical and β-folded proteins observed experimentally.10-12

Given the continued

evaluation of force fields13,14 and optimization studies of GB models,15 it is of general interest to determine if the combined CHARMM36m/GBMV2 simulation strategy provides a framework for modeling IDPs.16-18 While GBMV2 is a computationally efficient model as compared to explicit solvent simulations, questions remain if implicit solvent models can correctly shift the density of states of a prefolded IDP on an energy landscape hindered by multiple kinetic traps to an ensemble of disordered states favored by configurational entropy.


Page 3 of 24 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



To help with the assessment of the CHARMM36m/GBMV2 simulation model, two different strategies of applying ReX methods are compared where one is a traditional approach of a fixed set of sampling temperatures and the other is an adaptive scheme of allowing the windows to dynamically move in temperature space between two temperature extremes. The underlying idea of the adaptive-T method is an enhanced sampling technique to enrich the population of windows and their exchanges near a so-called “phase-transition” temperature to improve basin hopping across conformational states.

The scheme was first proposed by

Hansmann and coworkers,19 and Troyer and coworkers.20 Here, the implementation of the adaptive-T ReX algorithm is based on earlier work of modeling the unfolding-folding dynamics of a 57-residue protein SH3.21 It was observed from this study that the adaptive-T approach produced a significantly lower melting transition temperature than that of the traditional fixed-T method, leading to excellent agreement with experimental observations. Other noted applications of the adaptive-T ReX method include structure refinement of knowledge-based protein models22-23 and the simulation analysis of trapping a metastable conformational state in a crystallographic protein homodimeric assembly.24 Given the broad applicability of the adaptive-T ReX simulation strategy for modeling protein folding events, it is of general interest to determine whether this enhanced sampling approach is beneficial in generating an ensemble of less-cooperative transitions that are known to govern IDPs. To extend the work reported here, the CHARMM36m/GBMV2 simulation results are compared to a re-examination of CHARMM22/GBMV2 simulation datasets taken from an earlier study.25 The comparison centers on computing potentials of mean force (PMFs) using the parallel tempering weighted histogram analysis method26 (PTWHAM) and the multistate Bennett acceptance ratio (MBAR) method.27

2. METHODS Replica-Exchange Strategies. The temperatures of the replica clients in a traditional approach are fixed at specific T values during the simulations and can be modeled as a geometrical spaced set of N-1 intervals from the minimum temperature Tmin to the maximum Tmax



Page 4 of 24


T1 = Tmin  1 

(1)

 T   N −1 Ti +1 = Ti  max  , T  min 

where i = 1…N and the number of replica clients is given by N. After a specified number of simulation integration time steps, the neighboring replica clients, a and b, swap temperatures with a probability given by the Metropolis energy criteria28

p ( a ↔ b ) = min 1, e( 

βa − βb )( Eb − Ea )

, 

(2)

where βa=1/ kBTa, kB is Boltzmann’s constant, Ta is the temperature of replica client a, and Ea is the potential energy of client a and similarly for client b. In an adaptive scheme Eq. (1) is replaced by an algorithm that attempts to maximize the number of times replica clients progress in round trips from Tmin to Tmax. Depending on the last temperature extreme visited, each replica client is labeled as either cold or hot.21 Histograms ncold (T ) and nhot ( T ) over temperature space are constructed to collect the number of clients

visiting each temperature window. The fraction of cold clients visiting temperature T is modeled as

f (T ) =

ncold (T ) . ncold (T ) + nhot (T )

(3)

To model the movement of clients, a thermal current can be defined from Eq. (3) by j = D ( T )η ( T )

df , dT

(4)

where D (T ) is the diffusivity and η (T ) is the probability that any client will reside at T. The current j can be maximized by altering temperatures such that f ( Ti ) increases linearly as a function of temperature index i. A continuous f ( T ) is constructed from an interpolation of computed values at the current Ti, and then a search is conducted for the new temperatures where f ( T ) = i/(N ˗ 1). To avoid clustering of all of the windows around the same temperature, a

constraint is introduced such that no neighboring temperatures can be more than two geometric spacing units apart,


Page 5 of 24 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



Ti +1  Tmax  ≤  Ti  Tmin 

 2   N −1  

,

(5) 21

with the lower and upper values of Ti set to Tmin and Tmax, respectively.

CHARMM-based Simulations. For conformational sampling of the 28-residue NPBP peptide, the self-guided Langevin dynamics (SGLD) method developed by Wu and Brooks29-30 was combined with ReX using a strategy first reported in modeling the Trp-cage mini-protein.11 The SGLD equation of motion is given by p& = f − γ p + R + λg , i i i i i i

(6)

where p& defines the rate of change of the momentum of particle i, f is the force acting on the i i particle, γ is the friction constant, R defines the random force and g is a memory function, i i i which is scaled by an ad hoc guiding factor λ.

The memory function g is defined by the i

moving average of the momentum seen by the system over an interval of time, L:

g = p i i

,

(7)

L

where K L denotes a local average. The time interval is further defined as L = t / δt , where L t is the local averaging time and δt the time step along the simulation trajectory. L

The simulation methodology and parameter set applied here are similar to that given in an earlier study of the NPBP peptide using the CHARMM22 force field with the GBMV2 solvent model.25 Here, CHARMM36m is applied using the CHARMM simulation program31 (version c41b1) and the utilities and programming libraries of the Multiscale Modeling Tools for Structural Biology (MMTSB).32 An integration time step of 2 fs was used and parameters for SGLD consisted of the friction constant set to 1 ps-1 for all heavy atoms, the guiding factor λ in Eq. (6) was set to a value of 1, and the time average of the momentum ( t ) was set to 1 ps. L Non-bonded interaction cutoff parameters for electrostatics and nonpolar terms were set at a radius of 22 Å with a 2-Å potential switching function. Covalent bonds between the heavy atoms and hydrogen atoms were constrained by the SHAKE algorithm.33 GBMV2 parameters within MMTSB were set to gbmvabeta = –12 and gbmvap3 = 0.65. Modifications of these 5 ACS Paragon Plus Environment


Page 6 of 24


parameters were first explored by Chocholousˇova´ and Feig34 in a study of simulation strategies to enhance energy conservation in the application of GBMV2 in molecular dynamics. The parameter gbmvabeta is defined as a scaling factor to make the molecular volume sharper or smoother at the protein-solvent dielectric boundary. Likewise the gbmvap3 scaling parameter can be adjusted to match Poisson solvation energies. The effect of varying both parameters in molecular dynamics simulations is presented in a study on modeling the heat capacity of SH3 by Yeh et al.10 The last parameter in GBMV2 is the hydrophobic cavitation term which was modeled by applying a phenomenological surface tension coefficient set to a value of 0.015 kcal/mol/Å2. Simulations were carried out using 24 replica clients and the frequency of exchanges was set to every 1 ps of simulation. Temperatures were geometrically spaced between Tmin = 300 K and Tmax = 475 K. The NPBP peptide was modeled for 200 ns of simulation time per client, generating an ensemble of 4.8 µs. Culled structures consisted of 200000 per temperature and were used in the analysis of computing the secondary structure by the DSSP algorithm35 and radius of gyration (Rg). PMFs (denoted as the measure WT at temperature T) were computed using order parameters of fractional helicity (fH), Rg and Z-score of the potential energies (E) from the CHARMM36m and CHARMM22 force fields as input to PTWHAM and MBAR calculations.

Further evaluations were applied using an order parameter of fractional fold

comprised of α-β secondary-structure elements from DSSP. It should be noted that the PMFs generated by the SGLD method are prone to small deviations from a canonical description due to the ad hoc force term in Eq. (6) which is added to achieve enhanced sampling (see, e.g., Refs. 11 and 30). While an algorithmic scheme has been reported to reweight the biases of local averages of forces and momenta,36 the application is cumbersome and exceeds the purpose of this work where the deviation is thought to be small for a peptide.11 Further analysis of the generated conformations was conducted using MMTSB to evaluate the population density of states by clustering methods that included hierarchical clustering (conformational and potential energy values) and k-means clustering (conformational). Additional examination was conducted of the root-mean-square deviation (RMSD) in selected backbone angles Φ and Ψ from the initial folded peptide structure. The RMSD values were


Page 7 of 24 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



computed across the ensemble for two peptide segments and provided input into a PTWHAM calculation.

3. RESULTS and DISCUSSION The initial conformation of the NPBP peptide (numbered from residues 20 to 47) bound to Ebola virus protein NP is a topology of a helix-β-turn-helix fold and shows a fractional helicity of fH = 0.43 from DSSP. The two helices are partitioned as Trp28 to Thr35 and Val40 to Ile43. The initial fold compactness is given by the dimension Rg = 10 Å and appears more collapsed than an estimate for a comparable unfolded 28-residue peptide showing Rg ~ 13 Å.37

Potentials of Mean Force for Helix Formation. Given the reported experimental CD spectrum without TFE showing the lack of α-helix formation in contrast to that found in the bound NPBP and the observation with TFE in stabilizing a hidden helix propensity,5 the simulation models are first assessed by the calculation of PMFs using an order parameter of fractional helicity. While the CD spectrum data is a qualitative description of the ensemble of dynamic conformations, comparison with simulations give valuable insight on the accuracy of the force fields and sampling methods to model the formation of secondary structure and their transitions. Illustrated in Figure 1 are PMFs at 300 K evaluated by using PTWHAM. Comparison between the fixed-T ReX and the adaptive-T method for CHARMM36m shows a helical fold profile with an overall similar formation for the two sampling algorithms. For the fixed-T method, the minimum in WT(fH,Rg) is observed at fH = 0.26 with a Rg = 7.7 Å, while the adaptiveT ReX gives an identical fH = 0.26 with a Rg = 7.8 Å, thus showing convergence in searching for

the minimum. There are, nevertheless, several differences that are detected between the two simulation models. The transition from the minimum in WT(fH,Rg) to a unstructured state where fH = 0 and Rg is set to an identical value as in the minimum yields for the fixed-T ReX a

unfolding free-energy difference of 3.3 kBT and for the adaptive approach, a value of 2.2 kBT. An appearance of this lower free-energy transition is noticeable changes in the potential energy Z-score profiles where the adaptive method provides greater density population at low fractional helicity.



Page 8 of 24


There are significant distinctions in the PMFs calculated from CHARMM36m and CHARMM22, as well as between the fixed-T ReX and the adaptive-T method for CHARMM22. For the fixed-T method, CHARMM22 shows the minimum at fH = 0.53 and Rg = 9.8 Å. A connecting local minimum is detected at WT(fH = 0.37,Rg = 7.9 Å) with a free-energy difference of only 0.05 kBT from the global minimum.

By comparison, the adaptive-T ReX for

CHARMM22 yields a minimum WT(fH = 0.37,Rg = 8.8 Å) with a nearby basin at WT(fH = 0.16,Rg = 8.8 Å) with a 0.13 kBT free-energy difference from the minimum. Conformational transition to an unstructured state WT(fH = 0,Rg) from the minimum for CHARMM22 using a fixed-T approach yields 1.1 kBT and the adaptive-T method yields 1.5 kBT. Both values are noticeably lower than the corresponding values from CHARMM36m, indicating differences in basin funnel depth. From the plots in Figure 1, CHARMM36m lowers the weight of helix propensity as compared to CHARMM22, yet the conformational plasticity along the helix-forming order parameter is favored by CHARMM22. The application of the adaptive-T ReX for CHARMM22 shows a manifold of connecting multiple basins, while sampling from CHARMM36m is more confined to discrete states and shows a strong bias in funneling conformations to a helical basin which can be characterized as a large kinetic trap. The distinction is evident in comparing the potential energy Z-score landscapes, where CHARMM22 simulation shuttles conformations among the major basins through a pathway of favorable exchanges among the probability density distribution. To check the accuracy of PTWHAM and further characterize the energy landscape topology in applying the CHARMM force fields with GBMV2, WT(fH,Rg) was recalculated by using the MBAR method. The analysis shows overall consistency for CHARMM36m with an equivalent minimum WT(fH = 0.26,Rg = 8.0 Å) as calculated from PTWHAM and similar values for CHARMM22 to yield WT(fH = 0.46,Rg = 8.6 Å). One notable difference for both force fields and sampling methods is that MBAR yields deeper helical basins when transitioning to unstructured conformations. The analysis shows ~8 kBT for CHARMM36m and ~4 kBT for CHARMM22. Nevertheless, the rank order of free energies that model transitions among basins from MBAR is similar to that found by PTWHAM.

Finally for each simulation model,

comparing WT(fH,Rg) over the last 50 ns gives an estimate of the statistical variance in 8 ACS Paragon Plus Environment

Page 9 of 24 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



determining contours of the major free-energy basins and shows an approximate upper bound variance of ± 0.8 kBT.

Analysis of Parallel Tempering Strategies. To develop a better understanding of the methodologies in parallel tempering and their application to conformational sampling, Figure 2 displays data extracted from the simulations. The first plot (Figure 2A) illustrates the final replica exchange temperature profiles for the fixed-T ReX and the adaptive-T method. The results show that allowing the windows to move in temperature space during a simulation produces significant deviations from the traditional method and the dissimilarity depends on the applied force field. Application of CHARMM36m yields a broad-spanning cluster of windows at a lower temperature profile than that observed from CHARMM22. Figure 2B shows the average temperature window computed over the 24 clients along the simulation trajectory. Combined the plots offer a unique perspective on the force fields and provide a measure of temperature-induced biases required to achieve conformational diffusion across potential energy barriers in locating the free-energy minimum. Figure 2C illustrates the replica exchange rates and as expected the fixed-T method shows a flat exchange profile for both force fields.

Of the two adaptive-T simulations,

CHARMM36m yields a maximum exchange of ~96% in a critical sampling regime of the potential energy, while for CHARMM22 the maximum is ~93%. In both, the high exchanges are compensated with significantly reduced exchanges at the temperature end points to prevent aggregate clustering of all clients.21 As an estimate of the overall exchanges and the mobility of replica clients, Figure 2D shows the fraction of cold clients f(T) computed from Eq. (3) for the CHARMM36m and CHARMM22 simulations.

The simulations follow an optimal linear

relationship anticipated of IDPs due to a lack of a sharp phase-transition temperature common of two-state folding events.21 The drop at high temperatures, while noticeably less abrupt than that observed for unfolding proteins with well-defined tertiary folds,21 is attributed to limited roundtrip excursions by the clients between the extremes, typifying a universal problem of parallel sampling methods of applying all-atom protein models.

Peptide Fold Analysis. Displayed in Figure 3 is the starting NPBP conformation bound to the Ebola virus protein NP. As noted above, NPBP is composed of two helical regions denoted as α1 of residues Trp28 to Thr35 and α2 of residues Val40 to Ile43. Illustrated in Figure 9 ACS Paragon Plus Environment


Page 10 of 24


3B are conformations extracted from hierarchical clustering of the CHARMM36m/GBMV2 generated ensemble selected at 300 K from the adaptive-T simulation. Shown are population estimates (P) for each cluster together with extracted structures that best represent the median energies from the CHARMM36m/GBMV2 potential energy function. A comparison is made with structures from CHARMM22/GBMV2 adaptive-T sampling (Figure 3C). The CHARMM36m/GBMV2 adaptive-T simulation shows the conformation located at the sampled potential energy minimum (Emin) contains a short helix (fH = 0.25) in the α1 region. For the two largest clusters, short helical conformations are likewise observed in α1; however one of the clusters transitions to conformers exhibiting short β-hairpins in the C-terminal region including α2. Further transitions are observed among the clusters in both extending and shorting the helical length. By comparison, the Emin of the CHARMM22/GBMV2 adaptive-T sampling simulation contains a helix-turn-helix fold of fH = 0.46 and the largest cluster is of similar general fold. Figure 3D reports the Cα RMSD of structures culled from the simulation models at 300 K relative to the prefolded NPBP conformation. The analysis shows the RMSD values are primarily in the range of 5-9 Å and the CHARMM36m/GBMV2 model with adaptive sampling yields the highest net value, while the adaptive method using CHARMM22/GBMV2 accounts for the lowest value. All model simulations show good convergence with statistical deviations no greater than ± 1 Å. The observation of β-hairpins in conformational sampling of NPBP is unique to the CHARMM36m/GBMV2 energy function. Moreover, the population density depends on the ReX method for sampling. Applying the fixed-T simulation approach, structures culled at 300 K are observed to yield a total P = 5% for populating β-hairpins, whereas the adaptive method found a strikingly higher population of P = 42%. Of the latter dataset, sampled conformations are largely 3-stranded β-hairpins.

The significant difference between the two sampling methods can

conceivably be understood from the temperature profiles reported in Figures 2A and B. Maintaining a large pool of thermal windows with high Metropolis exchanges around ~350 K rather than annealing via stepwise hierarchical windows allows for more localized diffusion and thus greater attempts confined at basin-hopping transitions to overcome an entropic-driven barrier from cooperative interactions of forming a short β-hairpin.


Page 11 of 24 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



To further investigate the folding transitions found by CHARMM36m/GBMV2 using the adaptive-T sampling, Figure 4 illustrates a PTWHAM assessment of RMSD in backbone angles Φ and Ψ from the prefolded peptide for 3-residue segments Ser30-Glu31-Gln32 (α1 region) and Val40-Ser41-Asp42 (α2 region). The analysis finds for α1 (Figure 4A) the highest populated basin to be located at small RMSD differences while overall the landscape shows conformational populations spread out among the profile. For α2 (Figure 4B), a strikingly different result is obtained where the landscape shows a “hot” population region of deviations largely removed from the initial structure, containing multiple transitions among short α2 helix folds and βhairpins. Figure 4C shows the conformational probability distribution computed from PTWHAM of applying an order parameter of fractional fold rather than helicity shown in Figure 1B. Changing the order parameter to include β-folds reveals a significant finding of shifting the minimum in WT(fH,Rg) to fH = 0 with a Rg = 10.4 Å. The weighted probability distribution calculated by modifying the order parameter reveals that the potential energies in PTWHAM are strongly correlated with fractional fold propensity. Transition from the minimum to the populated basin at WT(fH = 0.47, Rg = 7.9 Å) in Figure 4C yields a small free-energy barrier of only 0.03 kBT for interconversion between the unfolded state (designated as IV of Figure 4D) and an α1-

helix-β-hairpin fold (designated as III). illustrated in Figure 4D.

Other conformers with low-energy transitions are

The thermally-accessible transition IIII

IV in WT(fH,Rg) is

illustrative of dynamic conformational heterogeneity found of the CHARMM36m/GBMV2 adaptive-T model and is consistent with the CD spectrum reported without TFE.5 Moreover, the detected α-β folded state is surprisingly consistent with unbiased predictions of secondary structure by a consensus approach.36 The lack of β-folded states in the CD with TFE is not unexpected as the additive compound is known to preferentially over-stabilize the formation of helical structures.

Fold Stiffness and Compactness Propensity. Illustrated in Figure 5 are plots of peptide fold formation in terms of helicity and compactness as a function of sampling temperatures. Shown are statistical averages over datasets for CHARMM36m and CHARMM22 generated from both ReX methods.

A separate plot of fold formation is included from the



Page 12 of 24


CHARMM36m/GBMV2 adaptive-T simulation of reporting the contribution of α-β conformations. Also shown in both plots is an assessment of the statistical averages for the CHARMM36m fixed-T ReX simulation by overlaying values of fH and Rg where WT(fH,Rg) = 0 along the temperature profile computed from MBAR. As anticipated from the PMFs, the analysis shows the CHARMM22 models to exhibit an extended thermal stability in helix formation greater than CHARMM36m (Figure 5A). Unlike CHARMM22, CHARMM36m shows only minor difference between the adaptive-T and the fixed-T sampling of helicity. By including β-hairpin formation, the adaptive-T CHARMM36m results show a higher initial contribution followed by decay consistent with the transient nature of cooperative interactions of stabilizing a β-hairpin. This is followed by a profile similar to the fixed-T sampling for CHARMM3m. One of the concerns of applying implicit solvent descriptions to modeling large-scale conformational heterogeneity is that the sampled structures will be overly collapsed (see, e.g., comparison of solvent models in Ref. 25). Figure 5B shows all simulation models are plagued by extended compactness of the generated conformations. As a noted benchmark, the prefolded peptide conformation shows a Rg ~ 10 Å in the multimeric complex. Only the adaptive-T method with CHARMM36m found a free-energy minimum of an unstructured state with Rg ~ 10 Å in reasonable accord with a typical unfolded 28-residue peptide showing Rg ~ 13 Å.37 Conceivably to help avoid generating over-collapsed conformations, incorporating a volumebased nonpolar solvation free-energy term in GBMV2 could offset unfavorable biases.38

4. CONCLUSIONS This study investigated the application of the CHARMM36m force field with the GBMV2 implicit solvent model in a replica exchange simulation of calculating the conformational ensemble of a 28-residue IDP from the Ebola virus protein VP35. A comparison was provided of two parallel tempering methods where one is a standard fixed-T sampling protocol and the other is a strategy where the temperature windows move in temperature space adapting to conformational sampling of the potential energy function.

Comparisons were made with

equivalent CHARMM22/GBMV2 simulation datasets.


Page 13 of 24 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



The central issue of the study is the applicability of potential energy functions applied in parallel tempering algorithms as a computational approach for modeling large-scale conformational heterogeneity. The measure of success was the accuracy to replicate a disordered conformational ensemble of the peptide with a large manifold of thermally accessible states as reported from CD experiments. Starting from a helix-turn-helix topology, the results revealed that CHARMM36m using tempering methods produced a PMF of lower fractional helicity than CHARMM22, yet the latter force field captured more conformational plasticity along the helixforming reaction coordinate. Changing the order parameter to include β-sheet folding in addition to α-helix formation, only the adaptive tempering method with CHARMM36m found highly populated α-helix-β-hairpin transitions that led to a significant increase in conformational heterogeneity derived from multiple discrete states. Using this model, a PMF of fractional fold propensity was calculated and revealed a free-energy minimum containing an unstructured state in qualitative uniformity with circular dichroism experiments showing the lack of folded structures.

■ AUTHOR INFORMATION Corresponding Author *E-mail: (M.A.O.) [email protected]. Notes The author declared no competing financial interest.

■ ACKNOWLEDGMENTS Financial support for this work comes from US Department of Defense Threat Reduction Agency. The opinions or assertions contained herein are the private views of the author and are not to be construed as official or as reflecting the views of the US Army or of the US Department of Defense. This article has been approved for public release with unlimited distribution.



Page 14 of 24


■ REFERENCES (1) Bhowmick, A.; Brookes, D. H.; Yost, S. R.; Dyson, H. J.; Forman-Kay, J. D.; Gunter, D.; Head-Gordon, M.; Hura, G. L.; Pande, V. S.; Wemmer, D. E.; Wright, P. E.; Head-Gordon, T. Finding Our Way in the Dark Proteome. J. Am. Chem. Soc. 2016, 138, 9730–9742. (2) Arai, M.; Sugase, K.; Dyson, H. J.; Wright, P. E. Conformational Propensities of Intrinsically Disordered Proteins Influence the Mechanism of Binding and Folding. Proc. Natl. Acad. Sci. U.S.A. 2015, 112, 9614–9619.

(3) Wright, P. E.; Dyson, H. J. Intrinsically Unstructured Proteins: Re-Assessing the Protein Structure-Function Paradigm. J.Mol.Biol. 1999, 293, 321–331. (4) Wright, P. E.; Dyson, H. J. Intrinsically Unstructured Proteins and their Functions. Nat. Rev. Mol. Cell Biol. 2005, 6, 197–208.

(5) Leung, D.W.; Borek, D.; Luthra, P.; Binning, J. M.; Anantpadma, M.; Liu, G.; Harvey, I. B.; Su, Z.; Endlich-Frazier, A.; Pan, J.; Shabman, R. S.; Chiu, W.; Davey, R. A.; Otwinowski, Z.; Basler, C. F.; Amarasinghe, G. K. An Intrinsically Disordered Peptide from Ebola Virus VP35 Controls Viral RNA Synthesis by Modulating Nucleoprotein-RNA Interactions. Cell Rep.

2015, 11, 376–389. (6) Huang, J.; Rauscher, S.; Nawrocki, G.; Ran, T.; Feig, M.; de Groot, B. L.; Grubmüller, H.; MacKerell, A. D. Jr. CHARMM36m: An Improved Force Field for Folded and Intrinsically Disordered Proteins. Nat. Methods 2017, 14, 71–73. (7) Ishikawa, Y.; Sugita, Y.; Nishikawa, T.; Okamoto, Y. Ab Initio Replica Exchange Monte Carlo Method for Cluster Studies. Chem. Phys. Lett. 2001, 333, 199−206. (8) Lee, M. S.; Feig, M.; Salsbury, F. R., Jr.; Brooks, C. L., 3rd. New Analytic Approximation to the Standard Molecular Volume Definition and its Application to Generalized Born Calculations. J. Comput. Chem. 2003, 24, 1348−1356. (9) Feig, M.; Onufriev, A.; Lee, M. S.; Im, W.; Case, D. A.; Brooks, C. L., III. Performance Comparison of Generalized Born and Poisson Methods in the Calculation of Electrostatic Solvation Energies for Protein Structures. J. Comput. Chem. 2003, 25, 265−284. 14 ACS Paragon Plus Environment

Page 15 of 24 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



(10) Yeh, I. C.; Lee, M. S.; Olson, M. A. Calculation of Protein Heat Capacity from ReplicaExchange Molecular Dynamics Simulations with Different Implicit Solvent Models. J. Phys. Chem. B. 2008, 112, 15064−15073.

(11) Lee, M. S.; Olson, M. A. Protein Folding Simulations Combining Self-Guided Langevin Dynamics and Temperature-Based Replica Exchange. J. Chem. Theory Comput. 2010, 6, 2477−2487. (12) Chaudhury, S.; Olson, M. A.; Tawa, G.; Wallqvist, A.; Lee, M. S. Efficient Conformational Sampling in Explicit Solvent using a Hybrid Replica Exchange Molecular Dynamics Method. J. Chem. Theory Comput. 2012, 8, 677–687. (13) Song, D.; Wang, W.; Ye, W.; Ji, D.; Luo, R.; Chen, H-F. ff14IDPs Force Field Improving the Conformation Sampling of Intrinsically Disordered Proteins. Chem. Biol. Drug Des. 2017, 89, 5–15.

(14) Ye, W.; Ji, D.; Wang, W.; Luo, R.; Chen, H-F. Test and Evaluation of ff99IDPs Force Field for Intrinsically Disordered Proteins. J. Chem. Inf. Model. 2015, 55, 1021–1029. (15) Lee, K. H.; Chen, J. Optimization of the GBMV2 Implicit Solvent Force Field for Accurate Simulation of Protein Conformational Equilibria. J. Comput. Chem. 2017, 38, 1332– 1341. (16) Best, R. B. Computational and Theoretical Advances in Studies of Intrinsically Disordered Proteins. Curr. Opin. Struct. Biol. 2017, 42, 147–154. (17) Chong, S. H.; Chatterjee, P.; Ham, S. Computer Simulations of Intrinsically Disordered Proteins. Annu. Rev. Phys. Chem. 2017, 68, 117–134. (18) Bhowmick, A.; Brookes, D. H.; Yost, S. R.; Dyson, H. J.; Forman-Kay, J. D.; Gunter, D.; Head-Gordon, M.; Hura, G. L.; Pande, V. S.; Wemmer, D. E.; Wright. P. E.; Head-Gordon, T. Finding our Way in the Dark Proteome. J. Am. Chem. Soc. 2016, 138, 9730–9742. (19) Trebst, S.; Troyer, M.; Hansmann, U. H. Optimized Parallel Tempering Simulations of Proteins. J. Chem. Phys. 2006, 124, 174903–174909.



Page 16 of 24


(20) Katzgraber, G.; Trebst, S.; Huse, D. A.; Troyer, M. Feedback-Optimized Parallel Tempering Monte Carlo. J. Stat. Mech. Theory Exp. 2006, P03018. (21) Lee, M. S.; Olson, M. A. Comparison of Two Adaptive Temperature-Based Replica Exchange Methods Applied to a Sharp Phase Transition of Protein Unfolding-Folding. J. Chem. Phys. 2011, 134, 244111–224417.

(22) Olson, M. A.; Lee, M. S. Evaluation of Unrestrained Replica-Exchange Simulations using Dynamic Walkers in Temperature Space for Protein Structure Refinement. PLoS One.

2014, 9, e96638. (23) Olson, M. A.; Zabetakis, D.; Legler, P. M.; Turner, K. B.; Anderson, G. P.; Goldman, E. R. Can Template-Based Protein Models Guide the Design of Sequence Fitness for Enhanced Thermal Stability of Single Domain Antibodies? Protein Eng. Des. Sel. 2015, 28, 395-402. (24) Olson, M. A.; Legler, P. M.; Goldman, E. R.; Comparison of Replica Exchange Simulations of a Kinetically Trapped Protein Conformational State and its Native Form. J. Phys. Chem. B. 2016, 120, 2234-40.

(25) Olson, M. A. On the Helix Propensity in Generalized Born Solvent Descriptions of Modeling the Dark Proteome. Front. Mol. Biosci. 2017, 4, 3. doi: 10.3389/fmolb.2017.00003. (26) Chodera, J. D.; Swope, W. C.; Pitera, J. W.; Seok, C.; Dill, K. A. Use of the Weighted Histogram Analysis Method for the Analysis of Simulated and Parallel Tempering Simulations. J. Chem. Theory Comput. 2007, 3, 26–41.

(27) Shirts, M. R.; Chodera, J. D. Statistically Optimal Analysis of Samples from Multiple Equilibrium States. J. Chem. Phys. 2008, 129, 124105–10. (28) Metropolis, N.; Rosenbluth, A.; Rosenbluth, M.; Teller, A.; Teller, E. Equation of State Calculations by Fast Computing Machines. J. Chem. Phys. 1953, 21, 1087–1092. (29) Wu, X.; Brooks, B. R. Self-Guided Langevin Dynamics Simulation Method. Chem. Phys. Lett. 2003, 381, 512−518.


Page 17 of 24 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



(30) Wu, X.; Damjanovic, A.; Brooks, B. R. Efficient and Unbiased Sampling of Biomolecular Systems in the Canonical Ensemble: A Review of Self-Guided Langevin Dynamics. Adv. Chem. Phys. 2012, 150, 255−326. (31) Brooks, B. R.; Brooks, C. L., III; Mackerell, A. D., Jr.; Nilsson, L.; Petrella, R. J.; Roux, B.; Won, Y.; Archontis, G.; Bartels, C.; Boresch, S.; et al. CHARMM: The Biomolecular Simulation Program. J. Comput. Chem. 2009, 30, 1545−1614. (32) Feig, M.; Karanicolas, J.; Brooks, C. L., 3rd. MMTSB Tool Set: Enhanced Sampling and Multiscale Modeling Methods for Applications in Structural Biology. J. Mol. Graphics Modell.

2004, 22, 377−395. (33) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. C. Numerical Integration of the Cartesian Equations of Motion of a System with Constraints: Molecular Dynamics of n-Alkanes. J. Comput. Phys. 1977, 23, 327−341.

(34) Chocholousˇova´, J.; Feig, M. Balancing an Accurate Representation of the Molecular Surface in Generalized Born Formalisms with Integrator Stability in Molecular Dynamics Simulations. J. Comput. Chem. 2006, 27, 719−29. (35) Kabsch, W.; Sander, C. Dictionary of Protein Secondary Structure: Pattern Recognition of Hydrogen-Bonded and Geometrical Features. Biopolymers 1983, 22, 2577–2637. (36) Wu, X.; Brooks, B. R. Toward Canonical Ensemble Distribution from Self-Guided Langevin Dynamics Simulation. J. Chem. Phys. 2011, 134, 134108. (37) Kohn, J. E.; Millett, I. S.; Jacob, J.; Zagrovic, B.; Dillon, T. M.; Cingel, N., et al. Random-coil Behavior and the Dimensions of Chemically Unfolded Proteins. Proc. Natl. Acad. Sci. U.S.A. 2004, 101, 12491–12496.

(38) Kieslich, C. A., Smadbeck, J., Khoury, G. A., Floudas, C. A. conSSert: Consensus SVM Model for Accurate Prediction of Ordered Secondary Structure. J. Chem. Inf. Model. 2016, 56, 455–461. (39) Lee, M. S.; Olson, M. A. Comparison of Volume and Surface Area Nonpolar Solvation Free Energy Terms for Implicit Solvent Simulations. J. Chem. Phys. 2013, 139, 044119. 17 ACS Paragon Plus Environment


Page 18 of 24



Page 19 of 24 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 1. Conformational and energy landscapes computed from the four CHARMM/GBMV2 replica exchange simulation models where the order parameters are fractional helicity, radius of gyration and Z-score of potential energies. Plots display data extracted at sampling temperature T = 300 K. Simulation results are: (A) CHARMM36m fixed-T ReX; (B) CHARMM36m adaptive-T ReX; (C) CHARMM22 fixed-T ReX; and (D) CHARMM22 adaptive-T ReX. The free-energy scale is displayed where the color blue represents minima in the PMFs and is followed in order by colors green, yellow, and red, where the latter denotes high energy states of low population. 264x367mm (300 x 300 DPI)

ACS Paragon Plus Environment


Figure 2. Data extracted from parallel tempering. (A) Final replica exchange temperature sets for four model simulation methodologies. Blue colored line and symbols represent the adaptive-T method with the CHARMM36m force field; red colored line and symbols depict CHARMM22 with adaptive-T sampling; and black colored line and symbols depict the fixed-T method using either force field. (B) Average temperature over 24 replica clients during the simulation for each force field and simulation sampling method. Colored lines are the same as above. (C) Metropolis exchange rates computed for each simulation. Blue color denotes the CHARMM36m adaptive-T sampling simulation, red color denotes CHARMM22 adaptive method, green color denotes CHARMM36m using the fixed-T method, and purple color represents the CHARMM22 fixed-T method. (D) Fraction cold of replica clients as function temperature index. Blue and red colored lines and symbols denote CHARMM36m and CHARMM22 adaptive-T simulations, respectively, and the black dashed line represents a linear relationship. 166x145mm (300 x 300 DPI)


Page 20 of 24

Page 21 of 24 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 3. Analysis of peptide folds. (A) The initial prefolded conformation observed in the crystallographic assembly where the small peptide represents NPBP and the molecular surface represents Ebola virus NP. (B) Clustering of conformations taken at T = 300 K from the CHARMM36m/GBMV2 adaptive-T ReX simulation. The energy Emin denotes the minimum conformational potential energy sampled from the simulation. Values of P report the population percentage for each illustrated cluster. (C) Conformations extracted from the CHARMM22/GBMV2 adaptive-T ReX simulation. (D) Plots of conformation index versus Cα-RMSD (units of Å) at T = 300 K from the prefolded conformation where the blue colored line and gray data set represent the CHARMM36m/GBMV2 adaptive-T simulation results, black colored line represents CHARMM36m/GBMV2 fixed-T results, red colored line represents the CHARMM22/GBMV2 adaptive-T and the purple colored line represents the CHARMM22/GBMV2 fixed-T simulation. 254x338mm (300 x 300 DPI)




Page 22 of 24

Page 23 of 24 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 4. Population distribution landscapes and extracted conformations from the CHARMM36m/GBMV2 adaptive-T ReX results. (A) Plot of RMSD in Φ and Ψ space computed from the prefolded NPBP conformation for residues Ser30-Glu31-Gln32; (B) Plot of RMSD for residues Val40-Ser41-Asp42. (C) Conformational landscape of fractional fold containing α-helix and β-hairpin structures versus radius of gyration; (D) Conformations taken from the CHARMM36m/GBMV2 adaptive-T simulation. Colors applied in the PMFs and their scales are noted in Figure 1. 219x252mm (300 x 300 DPI)



Figure 5. Fold stiffness and compactness of generated conformations by the simulation models along the temperature profile. (A) Red colored line and circles represents a statistical average of fH from the CHARMM36m/GBMV2 fixed-T simulation, green colored line and symbols represent fH values from the CHARMM36m/GBMV2 adaptive-T simulation, black line and symbols represent the net peptide fold (fH + βhairpins) from the CHARMM36m/GBMV2 adaptive-T simulation, blue colored line and symbols represent values of fH where WT(fH,Rg) = 0 for CHARMM36m/GBMV2 fixed-T dataset computed using MBAR, longdashed line represents CHARMM22/GBMV2 fixed-T, and purple colored line and symbols represents CHARMM22/GBMV2 adaptive-T simulation. (B) Average Rg versus T using the line colors and symbols noted above. The blue colored line and symbols represent values of Rg where WT(fH,Rg) = 0 for CHARMM36m/GBMV2 fixed-T dataset computed using MBAR. Typical standard deviations are given by the CHARMM36m/GBMV2 fixed-T simulation where generating fH at T = 300 K is ~ 0.10 and ~ 0.02 at T = 475 K; for Rgthe standard deviation is ~ 0.38 Å at T = 300 K and for T = 475 K, deviation is ~ 2.64 Å. Error in determining the minima in WT(fH,Rg) from MBAR along the temperature coordinate is ~ 0.2 kBT.


Page 24 of 24

Page 25 of 24 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


232x330mm (300 x 300 DPI)



TOC graphics 34x19mm (300 x 300 DPI)


Page 26 of 24

Parallel Tempering of Dark Matter from the Ebola ... - ACS Publications

Recommend Documents