Universal Implementation of a Residue-Specific Force Field Based on

Subscriber access provided by Kaohsiung Medical University

Biomolecular Systems

Universal Implementation of a Residue-Specific Force Field Based on CMAP Potentials and Free Energy Decomposition Wei Kang, Fan Jiang, and Yun-Dong Wu J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b00285 • Publication Date (Web): 15 Jun 2018 Downloaded from http://pubs.acs.org on June 20, 2018

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 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 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.

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 33 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 Theory and Computation

Universal Implementation of a Residue-Specific Force Field Based on CMAP Potentials and Free Energy Decomposition Wei Kang,1,2 Fan Jiang,*,1 and Yun-Dong Wu1,2 1. Laboratory of Computational Chemistry and Drug Design, Laboratory of Chemical Genomics, Peking University Shenzhen Graduate School, Shenzhen, 518055, China. 2. College of Chemistry and Molecular Engineering, Peking University, Beijing, 100871, China.

*

Correspondence and requests for materials should be addressed.

1

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 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

Abstract: The coupling between neighboring backbone ϕ and ψ dihedral-angles (torsions) has been well appreciated in protein force field development, as in the correction-map (CMAP) potentials. However, although preferences of backbone torsions are significantly affected by side-chain conformation, there has been no easy way to optimize this coupling. Herein, we prove that the three-dimensional (3D) free energy hyper-surface of joint (ϕ, ψ, χ1) torsions can be decomposed into three separated 2D surfaces. Thus, each of the 2D torsional surfaces can be efficiently and automatically optimized using a CMAP potential. This strategy is then used to re-parameterize an AMBER force field such that the resulting χ1-dependent backbone conformational preference can agree excellently with the reference protein coil library statistics. In various validation simulations (including the folding of seven peptides/proteins, backbone dynamics of three folded proteins, and two intrinsically disordered peptides), the new RSFF2C (residue-specific force field with CMAP potentials) force field gives similar or better performance compared with RSFF2. This strategy can be used to implement our RSFF force fields into a variety of molecular dynamics packages easily.

2


Page 2 of 33

Page 3 of 33 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


1. Introduction Computational modeling of biomolecular systems are crucial for predicting their structures and understanding their dynamics and functional mechanisms. A vast majority of these studies rely on empirical potential energy functions (force fields). For these potential functions, optimizing the energetics of dihedral-angles (torsions) has attracted extensive attention, from previous developments of either physics-based1-4 or knowledge-based5-8 potentials, to recent modifications of both AMBER9-16 and CHARMM17-20 force fields. Most of them were focused on backbone ϕ, ψ torsions. Indeed, some important conformational behaviors of simulated biomolecules, such as the balance between different peptide/protein secondary structures, can be extremely sensitive to small changes of backbone torsion potentials.9, 12, 15 The use of dihedral-angle potentials can be regarded as a way to fine-tune the description of non-bonded interactions between atoms separated by three covalent bonds (1-4 interactions). Understandably, specific orbital-orbital interactions,21 local steric interactions,22 and electrostatics beyond point-charge approximation can have significant effects on torsional preferences, and these interactions are usually missing or non-optimal in classical force fields. Along this consideration, the coupling between neighboring torsions should also be optimized to correct the inaccurate description of specific 1-5 interactions. Indeed, two-dimensional (2D) grid-based potential energy function CMAP has been developed by MacKerell et al.,17-18 such that the quantum mechanics (QM) (ϕ, ψ) 2D energy surface of dipeptides can be much better reproduced. Furthermore, recently Best et al.23 found that incorporation of CMAP potential can give better cooperativity in α-helix formation. The torsion parameters in these physics-based protein force fields were largely derived from QM calculations of isolated (gas phase) model systems. However, the condensed-phase environment for their applications is quite different. We have tried another source of reference data: backbone and side-chain conformational distributions of each amino acid residue obtained from protein “coil library”. Different from QM results used by other force fields, our reference data are all free energies instead of energies. We have found that, not only backbone ϕ and ψ torsions in each residue are coupled, but there are also strong couplings between side-chain χ1 and backbone torsions. In most situations, different side-chain χ1 rotamers give 3



Ramachandran (ϕ, ψ) plots with difference larger than that between different amino acids under a same side-chain rotamer.24 In our previous development of residue-specific force fields (RSFFs),25-26 we already considered this coupling, through using special Lennard-Jones (L-J) parameters of local non-bonded (1-5 or 1-6) interactions. However, optimizing these parameters requires tedious manual adjustments, and the resulting force fields can only be implemented in MD software supporting these special L-J terms, such as Gromacs.27 Also, we could not rigorously fit the (χ1, ϕ, ψ) 3D free energy surfaces in our previous RSFFs. One important goal of this work is to implement our RSFFs into other MD software packages such that one can take full advantage of their different merits and features. Especially, Amber28 and OpenMM29 provide strong support for MD simulations (in both explicit30 and implicit31 solvent) using graphic processing unit (GPU) with amazing cost-performance. In addition, more implicit solvent models have been implemented in Amber, Charmm,32 OpenMM, etc. Some MD packages, including Amber and Charmm, have a variety of enhanced sampling methods. Furthermore, the openness and flexibility in OpenMM can facilitate the development of new simulation methods. In this work, we introduce an efficient free energy decomposition strategy to decompose these otherwise complex couplings into the couplings of (ϕ, ψ), (χ1, ϕ) and (χ1, ψ), respectively. Following this strategy, we demonstrate that the re-parameterization of coupled torsions can be significantly facilitated by using CMAP potentials, with no manual adjustments required. This optimization was performed on AMBER ff14SB, one of the currently most fine-tuned protein force fields. The resulting force field was termed RSFF2C, with its general non-bonded parameters as same as those in ff14SB and RSFF2 and the “C” indicating CMAP. A series of testing simulations were also performed to evaluate RSFF2C and enable comparisons with ff14SB and RSFF2.

2. Theory and Methods 2.1 Free energy decomposition Among the 20 common amino acids (AAs), 17 of them (except Gly/Ala/Pro) are Cβ-derivatives of alanine (Ala). In peptides and proteins, each of these Ala-derived residues 4


Page 4 of 33

Page 5 of 33 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


has three rotatable bonds around the Cα atom, related to backbone (ϕ, ψ) and side-chain (χ1) dihedral angles. In our theory, the intrinsic ϕ/ψ/χ1 conformational preferences of these AA residues are determined by the interactions among four fragments: two neighboring peptide groups, the Cβ and attached hydrogen atoms, and the substituent(s) on Cβ (R-group). As shown in Scheme 1, the three pair-wise interactions among Cβ and the two peptide groups are functions of backbone ϕ or/and ψ torsions, and they are related to the intrinsic conformational preferences of Gly and Ala. The remaining two interactions are from the amino-acid-specific R-group to the ϕ-related and the ψ-related peptide groups, respectively, which are functions of (χ1, ϕ) and (χ1, ψ), respectively.

Scheme 1. Pairwise additivity of the five interactions (shown as dashed arrows) determining local backbone (ϕ, ψ) and side-chain χ1 conformational preferences within a dipeptide model (blocked AA residue). Blue ellipsoids, brown sphere, and green sphere represent the peptide groups, Cβ (with attached hydrogen atoms), and the substituent(s) on the Cβ, respectively. For each interaction, the related torsion(s) is also shown.

It is reasonable to assume that the interaction energies of the above five pairs are additive (Scheme 1). Here we take a further step to assume that the related χ1/ϕ/ψ/ free energy surface can also be decomposed into three 2D surfaces: , , , , ,

(1)

where X is any Ala-derived AA type, and A stands for Ala. In the right-hand-side (rhs) of eq 1, the first term describes the intrinsic conformational preference of Ala, whereas the last two 5



Page 6 of 33

terms describe the effect of a certain side-chain substituent at a certain χ1 conformation. Note that the terms in eq 1 are free energies, which integrated over other degrees of freedom, such as bond-angle bending, χi (i>2) torsions, and solvation, including entropic contributions. Thus, its correctness cannot be rigorously proven, but we can test how well observations fit this assumption. Based on the Boltzmann relationship, eq 1 can be re-written in terms of probability:

, , , ∙ , ∙ ,

(2)

Here pX (χ1, ϕ, ψ) is the 3D distribution of residue type X and pA (ϕ, ψ) is the 2D distribution of Ala. Both of them can be obtained from statistical analysis (next section for details) of known protein structures. On the other hand, pX (χ1, ϕ) and pX (χ1, ψ) are hypothetic probabilities describing certain conformational properties of residue type X. We have found that pX (χ1, ϕ) and pX (χ1, ψ) can be solved from the known pX (χ1, ϕ, ψ) and pA (ϕ, ψ) by applying the following two equations iteratively:

, ∑

∑ ,,

, ∑

∑ ,,

, ∙ ,

, ∙ ,

(3a) (3b)

This hypothesis (eq 1 and eq 2) means that the information encoded in a 3D (χ1, ϕ, ψ) distribution has great redundancy. For example, if we use a grid size of 15°, 24 * 24 * 24 = 13824 values are needed to represent such a distribution. On the other hand, only 3 * 24 * 24 = 1728 values are needed to represent three 2D distributions.

2.2 Statistical analysis of protein coil library To validate our hypothesis and parameterize a new force field, 10116 crystal structures (resolution < 2.0 Å, R factor < 0.2, 50% sequence identity cutoff) were selected and downloaded from PDB database (up to 2013/10/22), and all residues having backbone atoms with B factor > 35 were excluded. For Ala and Gly, the same criteria used to construct the Coil-6 (most stringent) in our previous work33 were adopted here. In simple words, all residues preceding Pro, within or adjacent to any secondary structure were excluded, and the 15% most exposed ones and those adjacent to residues with type belonging to SDNVITFYHW (single letter AA code) were also excluded. For other amino acids, Coil-3+ was used instead because 6


Page 7 of 33 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 Coil-6 library does not contain a sufficient number of residues to obtain reliable statistics on 3D (χ1, ϕ, ψ) grids instead of 2D (ϕ, ψ) grids for Gly and Ala. This coil library was constructed from the Coil-333 (only excluding pre-Pro and secondary structure residues) by excluding the residues preceding any residues with -30⁰ < ψ < 30⁰, because the side-chain of Ser, Thr, Asp, and Asn can form H-bond with backbone H(N), causing local conformational bias. For other residue types, the statistics from Coil-3 and Coil-3+ have little difference. To extract the 3D (χ1, ϕ, ψ) distribution, we used a Gaussian kernel density estimator, considering the periodicity of dihedral angles:

∆ min!",# − ", 360° − ",# − ")

(4a)

∆ min|# − |, 360° − |# − |

(4b)

∆ min|# − |, 360° − |# − |

(4c)

+ , , ∑# ,# exp /−Δ1 ∆ 1 ∆ 1 /24 1 5

(4d)

here 6 counts for all residues of a certain residue type in the coil library, and ,# is the 1/m weighting for m identical chains in one PDB structure. χ1,i, ϕi and ψi are the observed dihedral angles of the ith residue in the coil library. 15⁰×15⁰×15⁰ grids of (χ1, ϕ, ψ) were used with σ = 10⁰. The 2D n(ϕ, ψ) distributions were obtained similarly with 15⁰×15⁰ grids. The obtained n(ϕ, ψ) and n(χ1, ϕ, ψ) are unnormalized distributions. To apply eq 3, they were first normalized using: 8, 9:

, ∑ 8,

(5a)

8 ,, 9:

, , ∑ 8

,,

(5b)

here the pseudocount ε = 0.02 was used to avoid infinite free energy.

2.3 Parameterization of RSFF2C Similar to our previous RSFFs, the parameterization is carried out such that the simulated conformational distributions of AA dipeptides (blocked amino acids, Ac-X-NHMe) in TIP3P water match those from statistical analysis of PDB coil library. All parameterizations were based on replica exchange molecular dynamics (REMD) simulation implemented in Amber 16,28 and simulation methods and settings are detailed in SI. The RSFF2C force field was developed by adding CMAP potentials to the parent force field AMBER ff14SB. The fitting 7



Page 8 of 33

procedure is shown in Scheme 2 with backbone CMAP optimization (left) followed by side-chain CMAP optimization (right). One important feature of the procedure is that, unlike the development of RSFF1 or RSFF2, it was free of manual adjustments thus could be fully automatized.

Scheme 2. The overall parameter optimization procedure of the RSFF2C force field. Left: backbone (ϕ, ψ) CMAPs of Ala and Gly were optimized first. Right: the final backbone CMAP of Ala was used for all Ala-derived AA to optimize their side-chain χ1-related CMAPs, including (χ1, ϕ) and (χ1, ψ). Replica-exchange molecular dynamics (REMD) was used for all dipeptide simulations and 120 ns trajectory from 300 K replica was used for subsequent analyses with the first 20 ns discarded.

All non-Gly amino acid residues share the same backbone CMAP of alanine, and Gly has its own CMAP. To minimize the difference between MD (ϕ, ψ) distribution and coil library statistics, the following CMAP term should be added as a compensation:

; , ?@AB , − =C , −DE ln GHIJ

,

KL ,

(6)

This strategy is significantly different from our earlier force field optimization method, 8


Page 9 of 33 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


where matching between MD and reference distributions was achieved by optimizing Fourier expansion coefficients and local 1-4/1-5/1-6 interactions. The actual optimization of ECMAP (ϕ, ψ) was iterative. Initially, ECMAP (ϕ, ψ) ≡ 0 was used with the parent force field ff14SB. During each optimization run, ∆ECMAP(ϕ, ψ), the difference between the (ϕ, ψ) free energy surface in the coil library and that in MD simulations, was calculated (eq 7a) and added to the CMAP from the previous run to obtain current CMAP (eq 7b). This process was repeated until the similarity S between MD and coil-library distribution was > 0.995.

∆;_A9 , ?@AB , − NNOPQ9_A ,

(7a)

;_A9 , ;_A, ∆;_A9 ,

(7b)

The similarity S was defined as same as that in our previous works:25-26

R

∑ 8STUV , ∙8WX , Y∑ 8STUV, Z ∙ Y∑ 8WX, Z

(8)

For any Ala-derived residues (including Pro), based on our free energy decomposition (eq 1), we can decompose the difference between 3D distributions from coil library and ff14SB simulation into three terms, each of which can be represented by a 2D CMAP term (like eq 6):

,?@AB , , − ,NNOPQ , ,

[,?@AB , ,?@AB , ,?@AB , \ −/,=C , ,=C , ,=C , 5

;, , ;, , ;, ,

(9)

When backbone CMAP of Ala, EA,CMAP (ϕ, ψ), was added to ff14SB, eq 9 becomes:

,?@AB , , − ,NNOPQ9_ , , ;, , ;, , (10) where GX,ff14SB+CMAP_A (χ1, ϕ, ψ) is the dipeptide’s free energy hyper-surface from MD simulation with ff14SB potential combined with backbone CMAP of Ala. Based on Boltzmann relationship, we can rewrite eq 10 into eq 11: ]

, , ∙ ] , , ,?@AB , , /,NNOPQ9_ , , (11)

Similar to eq 3, the δpX,CMAP (χ1, ϕ) and δpX,CMAP (χ1, ψ) were solved from the known pX, coil (χ1, ϕ, ψ) and pX,ff14SB CMAP_A (χ1, ϕ, ψ) by applying the following two equations iteratively: 9



] ]

, ,

∑ ,GHIJ ,,

Page 10 of 33

(12a)

(12b)

∑ ,hhijklmKn_ ,, ∙ o,mKn ,

, , ∑

∑ ,GHIJ ,,

,hhijklmKn_ ,, ∙

o

,mKn ,

Uniform distribution of δpX, CMAP (χ1, ψ) ≡ 1 was used as an initial guess, and the convergence can always be achieved within 10 iterations. Then δpX, CMAP (χ1, ϕ) and δpX, CMAP (χ1, ψ) were converted to 2D CMAP potentials separately:

;, , −DEln ] , ,

(13a)

;, , −DEln ] , ,

(13b)

Similar to the optimization of ECMAP (ϕ, ψ), EX,CMAP (χ1, ϕ) and EX,CMAP (χ1, ψ) were iteratively updated by a few cycles. The final backbone and sidechain CMAPs are shown in Figure S3. We also examined the PMFs of χi (i>1) of each AA residue, and we found that ff14SB gives PMFs in good agreement with coil library statistics for almost all AAs except for Glu, Gln, and Asn. Thus, these residues’ side-chain dihedral parameters were further optimized. Instead of CMAP, traditional Fourier expansions (up to fourth order) (eq 14) were added to the original functions. The PMFs before and after modifications are shown in Figure S4, along with the reference data from the coil library.

∆p# p?@AB # − pNNOPQ # ∑O8s q8 cos 8 #

(14)

Some amino acid residues (Asp, Glu, His, Lys, Arg, Cys) can have alternative protonation states. Currently we do not have sufficient samples from the protein coil library to perform reliable re-parameterization for their “non-standard” protonation states, except for His. For them and other special cases, such as phosphorylated residues or “non-natural” amino acids, previously developed AMBER parameters or generalized AMBER force field (GAFF) could be used compatibly, following the strategy in our previous works.34-35

3. Results and Discussion 3.1 Validation of the free energy decomposition To prove the aforementioned free energy decomposition does apply to the observed (χ1, ϕ, ψ) distribution, we first decompose the 3D distribution from protein coil library into (χ1, ϕ) and (χ1, ψ) by applying eq 3, based on coil library (ϕ, ψ) of Ala, then reconstruct a 3D distribution 10


Page 11 of 33 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


from these three 2D distributions by applying eq 2 (Figure 1A). For all the Ala-derived AA residues, the reconstructed (χ1, ϕ, ψ) distributions are very similar to the corresponding original (χ1, ϕ, ψ) distributions. The similarities (S) of the total (ϕ, ψ) distributions are > 0.99 for AAs except Asn. Also, S > 0.985 for all χ1-rotamer-dependent (ϕ, ψ) distributions were observed except for some cases of Asn, Ile and Val (0.977 for t rotamer of Asn; 0.968 and 0.982 for grotamers of Ile and Val, respectively). In addition, the χ1-rotamer preference can be excellently reproduced (Figure 1C). Interestingly, despite its special chemical structure, the energy decomposition also applies to proline residue excellently, with S = 0.999 and 0.997 for g+ and g- rotamers, respectively.

Figure 1. Our free energy decomposition using the coil library (ϕ, ψ, χ1) distribution. (A) Based on the coil library (ϕ, ψ) distribution of Ala, the residue-specific information encoded in the 3D (ϕ, ψ, χ1) distribution of glutamate residue (Glu, as an example) is extracted as two 2D distributions, and a very similar 3D distribution can be reconstructed from the resulting 2D 11



Page 12 of 33

distributions with the (ϕ, ψ) plot of Ala. The similarity between the original and reconstructed (ϕ, ψ) distribution (summation over all χ1 values) is 0.996. Throughout the paper, the contours are drawn every RT free energy difference (corresponding to the natural logarithm scale in terms of probabilities) (B) Similarities between the reconstructed and the original (ϕ, ψ) distributions, for all Ala-derived AA residues. Similarities for χ1 rotamer-dependent (ϕ, ψ) distributions are also given, with the definition: g+ (χ1 = 300⁰ ± 45⁰), t (χ1 = 180⁰ ± 45⁰) and g(χ1 = 60⁰ ± 45⁰). A few outlier AA residues are labeled with single letter code. (C) Percentages of g+, g-, and t rotamers in the reconstructed distributions against the original coil library statistics.

The obtained pX(χ1, ϕ) and pX(χ1, ψ) indeed show coupling between χ1 and ϕ or ψ dihedral-angles (Figure 1A, Figure S1), indicating side-chain and backbone conformational preferences are inter-dependent and cannot be further separated, in agreement with our previous finding.24 Finally, it should be noted that the obtained ∆G(χ1, ϕ) and ∆G(χ1, ψ) are not independent. They are coupled through χ1 and can be altered with following transformation:

∆ t , ∆ , u

(15a)

∆ t , ∆ , − u

(15b)

where f(χ1) is an arbitrary function of χ1. This affects the interpretation of the physical meaning of these two decomposed terms. It is important to note that our scheme (free energy decomposition) is quite different from the pairwise additivity assumption (potential energy decomposition) in current classical force fields. The ∆G(χ1, ϕ) and ∆G(χ1, ψ) not only describe the direct interactions between side-chain and backbone atoms, but also include important contributions from solvent effect. From our quantum mechanics (QM) calculations (Figure 2), including of solvent effect from water can significantly increase the energy variance at different ϕ or ψ values for a given χ1 side-chain rotamer. Understandably, when non-polar Cγ atom approaches polar backbone atoms, there are not only steric repulsions but also unfavorable desolvation of polar atoms. With the solvent effect, our QM calculations can reproduce observed χ1-dependent (ϕ, ψ) plots and χ1-rotamer preferences quite well,33 which cannot be achieved without solvent effect. These effective atomic interactions in solvent result in potential of mean forces (or free energy functions) 12


Page 13 of 33 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


instead of potential energies. Similar free energy additivity of distance-based or dihedral-angle-based PMFs is also the implicit assumption underlying knowledge-based potentials.

Figure 2. Local side-chain effects on ϕ and ψ torsions at three χ1 rotamers, calculated using simplified models (left) at MP2/6-31+G** level in gas-phase (middle) and with CPCM solvent model of water (right). Plus sign (+) and minus sign (-) are for g+ and g- rotamers respectively, and open diamond with solid line is for t rotamer. It should also be noted that our force field parameterization is not a direct application of the free energy decomposition in eq 1 ~ eq 3. The use of CMAP is to compensate errors of parent force field in describing certain torsion-torsion couplings by using normal local non-bonded interactions alone. The actual procedure used in our force field parameterization is described by eq 7 (for Ala and Gly) and eq 9 ~ eq 13 (for other AA residues). But we can still validate our parameterization strategy by applying a similar decomposition-and-reconstruction procedure to the difference between coil library and ff14SB free energy surfaces (Figure S2). Indeed, using both backbone and side-chain CMAPs to compensate this difference, we obtained S > 0.97 for all χ1-dependent (ϕ, ψ) distributions. On the other hand, when using only the backbone CMAP of Ala with one-dimensional ϕ/ψ/χ1 torsion potentials, the lowest S value is 0.82 and more than half of the obtained χ1-dependent (ϕ, ψ) distributions have S < 0.95 to the reference data. This clearly indicates the advantageous effect of using side-chain-related (χ1, ϕ) and (χ1, ψ) CMAPs. Without using any CMAP, even larger deviations from the reference data 13



were obtained.

3.2 CMAPs significantly improve similarities to coil library distributions Figure 3 shows the simulated (ϕ, ψ) distributions of the Gly and Ala dipeptides from AMBER ff14SB and RSFF2C. AMBER ff14SB gives (ϕ, ψ) distributions of Gly and Ala in good agreement with the coil library statistics (similarity S > 0.9), with S of Ala significantly higher than that from previous AMBER version ff99SB (S = 0.923 and 0.855 for Gly and Ala, respectively). However, closer examination reveals two discrepancies between ff14SB and coil library (ϕ, ψ) plots: 1) the diagonal shape of αR and PII (for Gly) basins in the coil library plots are not well reproduced by ff14SB; 2) for Ala with ϕ at around -160°, coil library plot gives probability at ψ = 160° (extended β) much higher than ψ = 0° (α′), which cannot be accurately described by ff14SB. These problems are related to ϕ/ψ coupling, and cannot be solved by optimizing dihedral-angle parameters alone. After one cycle of CMAP optimization, similarities to coil library plots can be improved to S = 0.99 (both Gly and Ala) using CMAP corrections. S > 0.995 can be achieved after the second cycle, and no further optimization is needed. Also, RSFF2C achieves significantly higher S values than RSFF2 (S = 0.97), demonstrating the advantages of using CMAP potentials over special 1-5/1-6 Lennard-Jones potentials.

14


Page 14 of 33

Page 15 of 33 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. Backbone (ϕ, ψ) distributions during the optimization of RSFF2C, compared with those from the coil library statistics. The (ϕ, ψ) distributions from AMBER ff14SB (the parent force field) and the RSFF2 (which is also parameterized using coil library statistics as the target) are also shown. The similarity to the corresponding coil library plot is given for each simulated (ϕ, ψ) plot. (A) (ϕ, ψ) distributions of Gly and Ala during the optimization of backbone CMAP. (B) χ1-rotamer-dependent (ϕ, ψ) distributions of aspartate residue (Asp) during the optimization of side-chain CMAPs.

Using of side-chain CMAPs improves the overall (three χ1-rotamers combined) and χ1-dependent (ϕ, ψ) plots (Figure 4A, Table S2 for details). In RSFF2C, 14 out of the 18 AA residues give overall similarities S > 0.990 and the lowest S = 0.985, while the resulting S 15



values from ff14SB are around 0.894 with large variations between 0.782 (Asn) and 0.979 (Pro). Results from dipeptide simulation with RSFF2 has an average S of 0.975, significantly lower than that of RSFF2C (0.993). Especially, the t rotamer of Asp and g- rotamer of Leu from RSFF2 simulations give S values around 0.80, which have been improved to be over 0.98 in RSFF2C. RSFF2C can also achieve good agreement between simulated and experimental 3

JHNHα coupling constants of dipeptides (Figure S6). When each rotamer is considered separately, ff14SB simulations give all g+ rotamers with

S > 0.80, while S values of t rotamers of Asp and Asn are as low as < 0.6 and g- rotamers of 7 AAs are < 0.80. This supports our notion that a single set of backbone torsion parameters may be insufficient to describe varied (ϕ, ψ) distributions of different rotamers simultaneously. On the other hand, RSFF2C maintains higher S values for all rotamers, better than RSFF2. Take aspartate residue (Asp) as an example, the rotamer-dependent (ϕ, ψ) distributions during optimization are shown in Figure 3B along with coil library reference and the corresponding RSFF2 distributions. The most striking difference among various plots with different force fields was observed in the case of t rotamer. Under t rotamer, both ff14SB and RSFF2 favor PPII conformation but coil library statistics shows a clear preference shifted to the C7eq conformer. The addition of residue-specific side-chain-related CMAPs in RSFF2C resolves this discrepancy in RSFF2. From RSFF2C simulations, only the g- rotamers of Val and Ile have S < 0.97 (Figure 4A). This indicates side-chain CMAPs in RSFF2C accurately captured the rotamer-dependent (ϕ, ψ) free energy profiles, despite the striking differences among different rotamers. For Val and Ile, each has two Cγ groups on the Cβ atom, and these β-branched AAs may not fit into our free energy decomposition scheme (Scheme 1) as perfect as others (Figure 1B). Nevertheless, RSFF2C still gives much better agreement with coil library distributions compared with RSFF2.

16


Page 16 of 33

Page 17 of 33 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. Similarities between backbone ϕ/ψ (A) and side-chain χ1 (B) distributions from simulations using three AMBER-related force fields (left: ff14SB, middle: RSFF2, right: RSFF2C) and those from coil library statistics. (A) Overall and χ1-dependent (ϕ, ψ) distributions. Outlier AA residues are labeled with single letter code. (B) Percentages of g+, g-, and t χ1-rotamers obtained from dipeptide simulations plotted against corresponding coil library statistics.

As expected, both RSFF2C and RSFF2 give intrinsic preferences of three side-chain χ1 rotamers (g+, t, and g+) very similar to corresponding distributions from coil library (Figure 4B). On the other hand, ff14SB gives significantly different results compared to coil library statistics, and, on average, t rotamers tend to have a larger proportion.

3.3 ab initio folding of various systems with RSFF2C To validate the performance of RSFF2C by ab initio folding simulations, we selected a series of peptides/proteins with a variety of sizes and secondary structures. Two Trp-cage variants (TC5b, TC10b),36 the chicken villin headpiece (villin HP35),37 and beta sheets of various sizes (CLN025,38 Trpzip-2,39 GB1 hairpin40 and a WW domain GTT variant41) were 17



chosen. As shown in Figure 5, RSFF2C can successfully fold all the seven peptides/proteins to their experimental structures with low backbone root-mean-square deviation (RMSD). The performance was similar to that of RSFF2. This is expected since that both force fields have used nearly the same reference data (coil library statistics) in their parameterizations.

Figure 5. Representative structures (rainbow) from RSFF2C simulations starting from fully unfolded conformations, compared with corresponding experimental structures (magenta) from PDB. The clustering analysis was performed with the replica at the lowest temperature. The centroid structure of the most populated cluster for each system is shown with its backbone RMSD given below, the percentage of this cluster is also given in parentheses.

As a comparison, similar simulations were performed with ff14SB for all testing systems except TC10b and villin HP35. Representative structures from clustering analysis were shown in Figure S7. ff14SB can fold all model peptides/proteins except WW GTT to their experimental structures with backbone RMSD values comparable to those from RSFF2C simulations, although the fraction of folded in GB1 hairpin is lower compared to the experimental value. For WW GTT, our long-time simulation of 24 µs (1 µs * 24 replicas), failed to sample the native state. Representative structures from major clusters are very different from the native topology (Figure S7).

18


Page 18 of 33

Page 19 of 33 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 6. Temperature dependence of helicity per residue in A14 (Ac-(Ala)14-NHMe, A) and AQ15 (Ac-(AAQAA)3-NH2, B). Melting curves of Trp-cage (TC5b, C), Trpzip-2 (D), and GB1 hairpin (E) from REMD simulations compared with experimental data.36, 39-40 Folding free energy surfaces (F) of GB1 hairpin at 303 K as the function of the Cα-RMSD to the reference structure (extracted from PDB structure of 1GB1) and the fraction of native contacts (Q). Representative structures for the top left basin in the free energy surface with ff14SB are shown on its right.

Figure 6 shows the calculated melting curves of five model peptides: two poly-alanine-based peptides (A14 and AQ15), TC5b, Trpzip-2, and GB1 hairpin. For A14 at 300 K, the helicity given by RSFF2C is approximately equal to the result (35%) predicted by AGADIR algorithm, slightly higher than the helix contents given by RSFF2 and ff14SB. For AQ15, RSFF2C and RSFF2 give quite similar results, and they both show enhanced cooperativity of helix formation over ff14SB.

Noticeably, ff14SB overestimated the stability

of TC5b, with the melting temperature (Tm) ~73 K higher than the experimental value. This is quite surprising, considering the low helicities of A14 and AQ15 given by the same force field. As for β-structures, ff14SB cannot well stabilize the folded state of GB1 hairpin (only 23% at 19



280 K), although it can reproduce the melting curve of Trpzip-2 quite well. In general, RSFF2C and RSFF2 give similar melting curves, especially when considering that they can be very sensitive9 to small changes in conformational energies. However, RSFF2C give somehow higher Tm for Trpzip-2 β-hairpin and slightly lower for Trp-cage (TC5b) and GB1 β-hairpin. It should be noted that these two force fields have same non-bonded parameters except for all the special 1-5/1-6 L-J interactions in RSFF2, and both were optimized against the conformational distributions from protein coil library. Two factors can contribute to this difference: 1) an updated protein coil library was used for RSFF2C, 2) a more rigorous fitting to the reference distributions from protein coil library was achieved with CMAPs in RSFF2C. Regarding the folding free energy surface of GB1 hairpin (Figure 6F), RSFF2 has a deep native-state basin and a very shallow basin on the top left, while in RSFF2C’s energy surface, this shallow basin gets slightly wider and deeper. ff14SB biased the sampling toward the top left energy basin in this system, which mainly contains α-helical conformations (Figure 6F). Overall, ff14SB over-stabilizes Trp-cage, gives the best performance in Trpzip-2, and underestimates the stability of GB1 hairpin. It may suffer from a mild secondary structure bias toward helical conformation. RSFF2 over-stabilizes all the three systems. RSFF2C can simultaneously stabilize the tested systems, although the two β hairpins are over-stabilized to some extent. From our recent study,42 the over-stabilization of β-structures may come from too strong intra-protein interactions43-46 due to the insufficient dispersion47 in common water models. The well-described stability of α-helical systems may come from the error cancellation with the lack of electronic polarization in helix formation. As Best et al.23 pointed out, the relation between force field parameters and enthalpy and entropy is not straightforward and torsions or CMAPs is just only one of many factors affecting the thermodynamics. Besides, the solvent model is another important factor because it plays a key part in striking a delicate balance between protein-protein and protein-solvent interactions. With a revised water model TIP4P-D47, our earlier work42 showed that RSFF2 gave melting curves of Trpzip-2 and the GB1 hairpin in good agreement with experiments. The inclusion of a small H-bond correction to the RSFF2 was needed to stabilize α-helix systems, to partially compensate the missing many-body effects. A similar strategy can be applied to RSFF2C. 20


Page 20 of 33

Page 21 of 33 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


However, an extensive investigation on that is beyond the scope of this paper. Also, in RSFF2C, all van der Waals parameters and atomic charges were adopted from ff14SB without modification. These aspects may need to be re-optimized to achieve better performance. Having considered the multiple factors affecting protein/peptide conformational behaviors, we should point out that more rigorous fitting of dihedral parameters alone does not guarantee improved performance in all systems if other parameters remain unchanged. Future optimization efforts in non-bonded interactions are needed to further improve current force fields.

3.4 Backbone dynamics behavior To investigate the backbone dynamics behavior with the new force field, 1-µs-long simulations for three folded proteins (GB3, Ubiquitin, and Lysozyme) were performed (see Figure S8 for their stabilities during simulations). Lipari−Szabo S2 order parameters of backbone NH were calculated for the folded proteins using AmberTools. First, rotdif analysis was performed to obtain the tumbling correlation time (τC). Then the 1-µs-long trajectory was evenly divided into five sections, and for each section (200 ns), S2 order parameters were obtained by averaging iRED calculation results with windows of length 5 * τC, following the procedure used by Maier et al.15 The calculated S2 order parameters for the three folded proteins are shown in Figure 7. RSFF2C and RSFF2 give similar results for GB3 (0.055 RMSD against NMR versus 0.051 RMSD), showing improved agreement with NMR data around residue K10 (located in the 1st turn) compared to ff14SB. With ubiquitin, the calculated order parameters with RSFF2C are comparable to the results with ff14SB (0.045 RMSD versus 0.046 RMSD for ff14SB), while RSFF2 gives considerably low values around residue G10 (located in the 1st turn) with overall RMSD as 0.058. In lysozyme, both ff14SB and RSFF2 show low values for residues 100‒104 (located in the 4th loop), but the same region is better reproduced by RSFF2C. However, compared to ff14SB and RSFF2, the results are slightly worsened in the 1st loop (around residue Y20) with RSFF2C. Overall, we conclude that RSFF2C shows similar or sometimes better order parameter reproduction compared to ff14SB and seems to achieve better agreement with NMR data than RSFF2 in loop regions. 21



Figure 7. Calculated NH S2 order parameters of GB3, ubiquitin, and lysozyme simulated with ff14SB, RSFF2, and RSFF2C, in comparison with NMR experimental data. RMSD against NMR data is shown for each plot.

In addition, two intrinsically disordered peptides (Histatin 5 and HIV-1 Rev) were simulated with TIP4P-D water. 3JHNHα coupling constants of ubiquitin, Histatin5, and HIV-1 Rev peptide were calculated via the Karplus relationship with the parameter set from Vögeli et al.48 3JHNHα coupling constants of ubiquitin are calculated from another simulation trajectory at 303 K, consistent with the experimental condition. As shown in Figure 8, with ubiquitin, RSFF2C reproduces the NMR data with the lowest RMSD value (0.91 Hz) among the three force fields, considerably better than ff14SB (1.33 Hz). RSFF2 gives results comparable to RSFF2C with 1.01 Hz RMSD. For Histatin 5 and HIV-1 rev, similar results were obtained with RSFF2 and RSFF2C with the former being slightly better (1.13 Hz and 0.80 Hz with RSFF2 versus 1.19 Hz and 0.86 Hz with RSFF2C). Interestingly, ff14SB/TIP4P-D gives results considerably different from NMR data (1.75 Hz and 1.02 Hz RMSD), systematically underestimating 3JHNHα coupling constants in the two disordered peptides.

22


Page 22 of 33

Page 23 of 33 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 8. Calculated 3JHNHα coupling constants of ubiquitin, Histatin5, and HIV-1 Rev peptide plotted against the corresponding experimental data. Ubiquitin simulations were performed with TIP3P water, while TIP4P-D water was used for the two IDP peptides. RMSD (in Hz) against NMR experimental data is shown for each line. A same color scheme was used for each plot, and the RMSD values are from ff14SB, RSFF2, RSFF2C (left to right) simulations.

From the testing results, we conclude that the backbone dynamics behaviors are well captured by RSFF2C/TIP3P. Also, without further re-parameterization, RSFF2C, optimized in TIP3P water, gives improved performance compared to ff14SB when used with TIP4P-D water model. Actually, conformational distributions of dipeptides from RSFF2C/TIP3P and RSFF2C/TIP4P-D simulations are very similar (Table S3), demonstrating certain compatibility with different water models.

3.5 Using RSFF2C with Gromacs and OpenMM Theoretically, due to the use of CMAPs instead of special L-J terms, the resulting force field can be used in any software supporting CMAP, including Amber,49 Charmm, Gromacs,50 23



OpenMM, NAMD,51 LAMMPS,52 etc. Here we choose Gromacs, and OpenMM, largely because the former and the latter are highly efficient in CPU and GPU calculations, respectively. The implementation of the RSFF2C force field in Gromacs is also based on the AMBER ff14SB force field. New atom types were created with their non-bonded and bonded (bond, angle, torsion) parameters copied from those of the original atom types, to allow the use of different CMAP potentials for different amino acid residues. Following this strategy, RSFF2C can be implemented in Gromacs as a stand-alone force field package (RSFF2C.ff). For our previous RSFF1/RSFF2/RSFF2+ force fields, in-house Python scripts are needed to convert the generated topology file, to add special pair (local non-bonded) interactions. Now, this is no longer needed, making the usage of our force field more convenient. Furthermore, the topology file generated by RSFF2C.ff can be directly parsed by the OpenMM software. Then, we simulated two model dipeptides (Ac-X-Nme, X = Ala or Gln) and α-helical peptide AQ15 using RSFF2C with different MD simulation packages. For AQ15 (Figure 6B), nearly the same melting curves (temperature-dependent helicities) were obtained from Amber simulations on both CPU and GPU, and Gromacs simulation on CPU. For small dipeptide systems (Figure 9), simulations using Amber, Gromacs and OpenMM all shows (ϕ, ψ) distributions highly consistent with coil library statistics, with similarities ≥ 0.99. Thus, it is feasible to use our residue-specific force fields cost-effectively on current GPU (here GTX1080 cards were used), achieving high performance.

24


Page 24 of 33

Page 25 of 33 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 9. (ϕ, ψ) distributions at 300 K from simulations with RSFF2C in different MD packages, compared with those from the coil library statistics. The similarity to the corresponding coil library plot is given for each simulated (ϕ, ψ) plot. (A) (ϕ, ψ) distributions of Ala. (B) χ1-rotamer-dependent (ϕ, ψ) distributions of Glutamine residue (Gln).

4. Conclusion The coupling of three neighboring torsions around the α-carbon is important in determining conformational behaviors of peptides and proteins.24 In this work, we first proposed that the 3D free energy hyper-surface of (ϕ, ψ, χ1) torsions within a residue can be decomposed into three 2D free energy surfaces. This assumption was then validated using observed conformational distributions in the coil library of protein crystal structures. This strategy can greatly reduce the numeric data needed to represent such hyper-surfaces. This dimensionality reduction cannot be entirely expected from the assumption of pairwise 25



additivity of interaction energies, because our reference data are free energies including the contributions from solvation. Following this novel free energy decomposition method, 2D CMAP potentials were derived from fitting each 2D free energy surface, and we obtained a new residue-specific force field, RSFF2C, as an improvement of AMBER ff14SB. Nearly perfect agreement with reference distributions from protein coil library can be achieved. Different from the early developed RSFF1 and RSFF2, the parameterization of RSFF2C is straightforward and can be fully automatized. The new force field RSFF2C can be easily implemented into different MD software packages supporting CMAP potentials. In the future, all our efforts in further development and applications of residue-specific force fields (RSFFs) will be based on the force field implementation described here. Finally, the couplings among three torsions around one tertiary carbon should also exist in many other organic and biological molecules. We hope this strategy can also be useful in developing new force fields for molecules other than proteins.

ASSOCIATED CONTENT Supporting information The detailed methods of our parametrization and validation of RSFF2C, and additional tables and figures: simulations settings (Table S1), similarities between coil-library and simulated distributions (Table S2) and between distributions from different water models (Table S3), results from free energy decomposition of reference 3D distributions (Figure S1), similarities between reference and reconstructed distributions using different decomposition strategies (Figure S2), backbone and side-chain CMAPs (Figure S3), side-chain PMFs for Glu, Gln and Asn (Figure S4), helicities of AQ15 (Figure S5), experimental and simulated 3JHNHα couplings of dipeptides (Figure S6), folding performance of AMBER ff14SB (Figure S7), the stabilities of three proteins with different force fields (Figure S8). The Supporting Information is available free of charge via the Internet at http://pubs.acs.org.

26


Page 26 of 33

Page 27 of 33 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


AUTHOR INFORMATION Corresponding Authors *E-mail: [email protected] (F.J.) Notes The authors declare no competing financial interest.

Acknowledgements We thank the financial supports from the Shenzhen Science and Technology Innovation Committee (JCYJ20170412150507046 for Y.-D.W.), and the National Natural Science Foundation of China (21573009 for F.J.).

Reference (1) Scheraga, H. A. Theoretical and experimental studies of conformations of polypeptides. Chem. Rev. 1971, 71, 195-217. (2) Gelin, B. R.; Karplus, M. Side-chain torsional potentials: effect of dipeptide, protein, and solvent environment. Biochemistry 1979, 18, 1256-1268. (3) Weiner, S. J.; Kollman, P. A.; Case, D. A.; Singh, U. C.; Ghio, C.; Alagona, G.; Profeta, S.; Weiner, P. A new force field for molecular mechanical simulation of nucleic acids and proteins. J. Am. Chem. Soc. 1984, 106, 765-784. (4) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. J. Am. Chem. Soc. 1996, 118, 11225-11236. (5) Sun, S. Reduced representation approach to protein tertiary structure prediction: statistical potential and simulated annealing. Journal of Theoretical Biology 1995, 172, 13-32. (6) Jha, A. K.; Colubri, A.; Freed, K. F.; Sosnick, T. R. Statistical coil model of the unfolded state: Resolving the reconciliation problem. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 27



13099-13104. (7) Fujitsuka, Y.; Chikenji, G.; Takada, S. SimFold energy function for de novo protein structure prediction: consensus with Rosetta. Proteins 2006, 62, 381-398. (8) Rohl, C. A.; Strauss, C. E. M.; Misura, K. M. S.; Baker, D. Protein Structure Prediction Using Rosetta. Methods Enzymol. 2004, 383, 66-93. (9) Best, R. B.; Hummer, G. Optimized Molecular Dynamics Force Fields Applied to the Helix−Coil Transition of Polypeptides. J. Phys. Chem. B 2009, 113, 9004-9015. (10) Lindorff-Larsen, K.; Piana, S.; Palmo, K.; Maragakis, P.; Klepeis, J. L.; Dror, R. O.; Shaw, D. E. Improved side-chain torsion potentials for the Amber ff99SB protein force field. Proteins: Struct., Funct., Bioinf. 2010, 78, 1950-1958. (11) Li, D.-W.; Brüschweiler, R. NMR-Based Protein Potentials. Angew. Chem. Int. Ed. 2010, 49, 6778-6780. (12) Nerenberg, P. S.; Head-Gordon, T. Optimizing Protein-Solvent Force Fields to Reproduce Intrinsic Conformational Preferences of Model Peptides. J. Chem. Theory Comput. 2011, 7, 1220-1230. (13) Best, R. B.; Mittal, J. Protein Simulations with an Optimized Water Model: Cooperative Helix Formation and Temperature-Induced Unfolded State Collapse. J. Phys. Chem. B 2010, 114, 14916-14923. (14) Li, Y.; Gao, Y.; Zhang, X.; Wang, X.; Mou, L.; Duan, L.; He, X.; Mei, Y.; Zhang, J. Z. A coupled two-dimensional main chain torsional potential for protein dynamics: generation and implementation. J. Mol. Model. 2013, 19, 3647-3657. (15) Maier, J. A.; Martinez, C.; Kasavajhala, K.; Wickstrom, L.; Hauser, K. E.; Simmerling, C. ff14SB: Improving the Accuracy of Protein Side Chain and Backbone Parameters from ff99SB. J. Chem. Theory Comput. 2015, 11, 3696-3713. (16) Song, D.; Luo, R.; Chen, H. F. The IDP-Specific Force Field ff14IDPSFF Improves the Conformer Sampling of Intrinsically Disordered Proteins. Journal of chemical information and modeling 2017, 57, 1166-1178. (17) MacKerell, A. D.; Feig, M.; Brooks, C. L. Improved Treatment of the Protein Backbone in Empirical Force Fields. J. Am. Chem. Soc. 2004, 126, 698-699. (18) Mackerell, A. D., Jr.; Feig, M.; Brooks, C. L., 3rd Extending the treatment of backbone 28


Page 28 of 33

Page 29 of 33 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


energetics in protein force fields: limitations of gas-phase quantum mechanics in reproducing protein conformational distributions in molecular dynamics simulations. J. Comput. Chem. 2004, 25, 1400-1415. (19) Best, R. B.; Zhu, X.; Shim, J.; Lopes, P. E. M.; Mittal, J.; Feig, M.; MacKerell, A. D. Optimization of the Additive CHARMM All-Atom Protein Force Field Targeting Improved Sampling of the Backbone ϕ, ψ and Side-Chain χ1 and χ2 Dihedral Angles. J. Chem. Theory Comput. 2012, 8, 3257-3273. (20) Huang, J.; Rauscher, S.; Nawrocki, G.; Ran, T.; Feig, M.; de Groot, B. L.; Grubmuller, H.; MacKerell, A. D., Jr. CHARMM36m: an improved force field for folded and intrinsically disordered proteins. Nat Methods 2017, 14, 71-73. (21) Pophristic, V.; Goodman, L. Hyperconjugation not steric repulsion leads to the staggered structure of ethane. Nature 2001, 411, 565-568. (22) Mo, Y. Rotational barriers in alkanes. Wiley Interdisciplinary Reviews: Computational Molecular Science 2011, 1, 164-171. (23) Best, R.; Mittal, J.; Feig, M.; D Mackerell, A. Inclusion of Many-Body Effects in the Additive CHARMM Protein CMAP Potential Results in Enhanced Cooperativity of α-Helix and β-Hairpin Formation. Biophys. J. 2012, 103, 1045-1051. (24) Jiang, F.; Han, W.; Wu, Y.-D. Influence of Side Chain Conformations on Local Conformational Features of Amino Acids and Implication for Force Field Development. J. Phys. Chem. B 2010, 114, 5840-5850. (25) Jiang, F.; Zhou, C. Y.; Wu, Y. D. Residue-specific force field based on the protein coil library. RSFF1: modification of OPLS-AA/L. J. Phys. Chem. B 2014, 118, 6983-6998. (26) Zhou, C. Y.; Jiang, F.; Wu, Y. D. Residue-specific force field based on protein coil library. RSFF2: modification of AMBER ff99SB. J. Phys. Chem. B 2015, 119, 1035-1047. (27) Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E. GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX 2015, 1-2, 19-25. (28) Case, D. A.; Cerutti, D. S.; Cheatham, T. E., III; Darden, T. A.; Duke, R. E.; Giese, T. J.; Gohlke, H.; Goetz, A. W.; Greene, D.; Homeyer, N.; Izadi, S.; Kovalenko, A.; Lee, T. S.; LeGrand, S.; Li, P.; Lin, C.; Liu, J.; Luchko, T.; Luo, R.; Mermelstein, D.; Merz, K. M.; 29



Monard, G.; Nguyen, H.; Omelyan, I.; Onufriev, A.; Pan, F.; Qi, R.; Roe, D. R.; Roitberg, A.; Sagui, C.; Simmerling, C. L.; Botello-Smith, W. M.; Swails, J.; Walker, R. C.; Wang, J.; Wolf, R. M.; Wu, X.; Xiao, L.; York, D. M.; Kollman, P. A. AMBER 2017, University of California: San Francisco, 2017. (29) Eastman, P.; Swails, J.; Chodera, J.; T McGibbon, R.; Zhao, Y.; A Beauchamp, K.; Wang, L.-P.; Simmonett, A.; Harrigan, M.; D Stern, C.; P Wiewiora, R.; R Brooks, B.; Pande, V. OpenMM 7: Rapid development of high performance algorithms for molecular dynamics. PLoS Comput. Biol. 2017, 13, e1005659. (30) Salomon-Ferrer, R.; Götz, A. W.; Poole, D.; Le Grand, S.; Walker, R. C. Routine Microsecond Molecular Dynamics Simulations with AMBER on GPUs. 2. Explicit Solvent Particle Mesh Ewald. J. Chem. Theory Comput. 2013, 9, 3878-3888. (31) Götz, A. W.; Williamson, M. J.; Xu, D.; Poole, D.; Le Grand, S.; Walker, R. C. Routine Microsecond Molecular Dynamics Simulations with AMBER on GPUs. 1. Generalized Born. J. Chem. Theory Comput. 2012, 8, 1542-1555. (32) Brooks, B. R.; Brooks, C. L.; MacKerell, A. D.; Nilsson, L.; Petrella, R. J.; Roux, B.; Won, Y.; Archontis, G.; Bartels, C.; Boresch, S.; Caflisch, A.; Caves, L.; Cui, Q.; Dinner, A. R.; Feig, M.; Fischer, S.; Gao, J.; Hodoscek, M.; Im, W.; Kuczera, K.; Lazaridis, T.; Ma, J.; Ovchinnikov, V.; Paci, E.; Pastor, R. W.; Post, C. B.; Pu, J. Z.; Schaefer, M.; Tidor, B.; Venable, R. M.; Woodcock, H. L.; Wu, X.; Yang, W.; York, D. M.; Karplus, M. CHARMM: The Biomolecular Simulation Program. J. Comput. Chem. 2009, 30, 1545-1614. (33) Jiang, F.; Han, W.; Wu, Y. D. The intrinsic conformational features of amino acids from a protein coil library and their applications in force field development. Phys. Chem. Chem. Phys. 2013, 15, 3413-3428. (34) Hu, K.; Geng, H.; Zhang, Q.; Liu, Q.; Xie, M.; Sun, C.; Li, W.; Lin, H.; Jiang, F.; Wang, T.; Wu, Y. D.; Li, Z. An In‐tether Chiral Center Modulates the Helicity, Cell Permeability, and Target Binding Affinity of a Peptide. Angew. Chem. Int. Ed. 2016, 55, 8013-8017. (35) Zeng, J.; Jiang, F.; Wu, Y.-D. Mechanism of Phosphorylation-Induced Folding of 4E-BP2 Revealed by Molecular Dynamics Simulations. J. Chem. Theory Comput. 2017, 13, 320-328. (36) Barua, B.; Lin, J.; D Williams, V.; Kummler, P.; Neidigh, J.; Andersen, N. The Trp-cage: Optimizing the stability of a globular miniprotein. Protein Eng. Des. Sel. 2008, 21, 171-185. 30


Page 30 of 33

Page 31 of 33 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


(37) Chiu, T. K.; Kubelka, J.; Herbst-Irmer, R.; Eaton, W. A.; Hofrichter, J.; Davies, D. R. High-resolution x-ray crystal structures of the villin headpiece subdomain, an ultrafast folding protein. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 7517-7522. (38) Honda, S.; Akiba, T.; Kato, Y. S.; Sawada, Y.; Sekijima, M.; Ishimura, M.; Ooishi, A.; Watanabe, H.; Odahara, T.; Harata, K. Crystal Structure of a Ten-Amino Acid Protein. J. Am. Chem. Soc. 2008, 130, 15327-15331. (39) G. Cochran, A.; J. Skelton, N.; A. Starovasnik, M. Tryptophan Zippers: Stable, Monomeric Beta -hairpins. Proc. Natl. Acad. Sci. U.S.A. 2001, 98, 5578-5583. (40) Munoz, V.; Thompson, P. A.; Hofrichter, J.; Eaton, W. A. Folding dynamics and mechanism of β-hairpin formation. Nature 1997, 390, 196-199. (41) Piana, S.; Sarkar, K.; Lindorff-Larsen, K.; Guo, M.; Gruebele, M.; Shaw, D. E. Computational design and experimental testing of the fastest-folding beta-sheet protein. J. Mol. Biol. 2011, 405, 43-48. (42) Wu, H.-N.; Jiang, F.; Wu, Y.-D. Significantly Improved Protein Folding Thermodynamics Using a Dispersion-Corrected Water Model and a New Residue-Specific Force Field. Journal of Physical Chemistry Letters 2017, 8, 3199-3205. (43) Piana, S.; Klepeis, J. L.; Shaw, D. E. Assessing the accuracy of physical models used in protein-folding simulations: quantitative evidence from long molecular dynamics simulations. Current opinion in structural biology 2014, 24, 98-105. (44) Best, R. B.; Zheng, W.; Mittal, J. Balanced Protein-Water Interactions Improve Properties of Disordered Proteins and Non-Specific Protein Association. J Chem Theory Comput 2014, 10, 5113-5124. (45) Rauscher, S.; Gapsys, V.; Gajda, M. J.; Zweckstetter, M.; de Groot, B. L.; Grubmuller, H. Structural Ensembles of Intrinsically Disordered Proteins Depend Strongly on Force Field: A Comparison to Experiment. J Chem Theory Comput 2015, 11, 5513-5524. (46) Henriques, J.; Skepo, M. Molecular Dynamics Simulations of Intrinsically Disordered Proteins: On the Accuracy of the TIP4P-D Water Model and the Representativeness of Protein Disorder Models. J Chem Theory Comput 2016, 12, 3407-3415. (47) Piana, S.; Donchev, A. G.; Robustelli, P.; Shaw, D. E. Water Dispersion Interactions Strongly Influence Simulated Structural Properties of Disordered Protein States. J. Phys. Chem. 31



B 2015, 119, 5113-5123. (48) Vögeli, B.; Ying, J.; Grishaev, A.; Bax, A. Limits on Variations in Protein Backbone Dynamics from Precise Measurements of Scalar Couplings. J. Am. Chem. Soc. 2007, 129, 9377-9385. (49) Crowley, M. F.; Williamson, M. J.; Walker, R. C. CHAMBER: Comprehensive support for CHARMM force fields within the AMBER software. Int. J. Quantum Chem 2009, 109, 3767-3772. (50) Bjelkmar, P.; Larsson, P.; Cuendet, M. A.; Hess, B.; Lindahl, E. Implementation of the CHARMM Force Field in GROMACS: Analysis of Protein Stability Effects from Correction Maps, Virtual Interaction Sites, and Water Models. J. Chem. Theory Comput. 2010, 6, 459-466. (51) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kalé, L.; Schulten, K. Scalable Molecular Dynamics with NAMD. J. Comput. Chem. 2005, 26, 1781-1802. (52) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1-19.

32


Page 32 of 33

Page 33 of 33 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


TOC graphic

33


Universal Implementation of a Residue-Specific Force Field Based on

Recommend Documents