Toward a general neural network force field for protein simulations: Refining the intramolecular interaction in protein

Molecular dynamics (MD) simulation is an effective and widely used technique to investigate proteins, such as protein folding and unfolding, protein–ligand binding, and protein conformational and structural analysis.^1–5 One of the most important factors related to the accuracy of results from MD simulations is the force field.^6–8 The traditional molecular mechanical (MM) force fields outperform the quantum mechanical (QM) or QM/MM force fields because they have much simpler energy representations, which are much more efficient on large protein systems. The potential energy described by MM force fields usually consists of bond, angle, dihedral, improper dihedral, and non-bonded energies such as electrostatic energy and van der Waals energy.⁹ The parameters of all energy terms are generally fitted to experimental data, such as macroscopic properties, spectroscopic and crystallographic data, and quantum mechanical (QM) calculations.¹⁰ Four major families of general force fields, AMBER,^11–14 CHARMM,^15–18 GROMOS,^19–22 and OPLS,^23–26 are widely used in MD simulations, and many improvements have been achieved to increase the applicability of these force fields over the years. In addition to the traditional fixed-charge MM force fields, polarizable force fields were developed to account for the polarization effect²⁷ based on fluctuating charges,²⁸ Drude oscillators,²⁹ inducible dipoles,³⁰ and methods including multipole electrostatics.³¹ 

Even though the current protein force fields are applicable to many problems, they are still not accurate enough and could fail in some situations.^32,33 For example, protein folding pathways and thermodynamics cannot be identified very well,^34,35 and the strength of protein–water interactions is not balanced well compared to protein–protein interactions.^36,37 In addition, most force fields have deficiencies in simulating intrinsically disordered proteins (IDPs), which lack well-defined three-dimensional structures but have full functionalities.^38,39 Since an IDP can only be characterized with an ensemble of flexible interconverting conformations, the MD simulation approach would be an ideal approach.^40–42 This suggests the urgency of developing more accurate force fields. The current force fields generally generate a too compact ensemble of IDP,^43–46 and two types of strategies were generally applied to improve force fields for IDP.^47,48 The first one is adjusting backbone dihedral or energy correction map (CMAP) parameters to avoid the high propensity of secondary structures for α-helix or β-sheet, such as ff14IDPSFF,⁴⁹ CHARMM36IDPSFF,^50,51 OPLSIDPSFF,⁵² RSFF2C,⁵³ and ESFF.⁵⁴ The second strategy is refining protein–water interactions by adjusting non-bonded parameters to achieve a better description of protein–water interactions, such as a99SB-UCB,⁵⁵ ff03ws,³⁷ and KBFF.^56–58 Moreover, specific coarse-grained force fields were also developed for IDP,^59,60 such as AWSEM-IDP⁶¹ and MOFF-IDP.⁶² With the continuous development of force fields, the simulation results for IDP were improved, but there are still unsolved issues, such as the limitation of studying the post-translation modifications in IDP and the difficulty of obtaining balanced local and global structural features.^41,48 This is caused by the simplicity and limitations of force field forms in MM representations. According to an investigation on the structural reorganization energy defined by the energy difference between the MM optimal structure and the QM optimal structure, the minimization of QM optimal structures with the major AMBER, CHARMM, GROMOS, and OPLS force fields all resulted in large reorganization energies, which suggests that the MM force field forms need to be improved.⁶³ 

Recently, machine learning techniques were widely and successfully integrated into force field development because of their capability to describe potential energy surfaces (PESs) at a high level of theory at a low computational cost.^64–68 Although machine learning methods can be utilized to adjust MM force field parameters,^69–72 our focus is primarily on the construction of general force fields with machine learning models because of the potential improvement of functional form. The common machine learning models used to build force fields are divided into two categories: kernel-based methods and artificial neural networks (NNs). One example of kernel-based methods is the Gaussian approximation potential, which is approximated by the Gaussian process regression on different descriptors and kernels,⁷³ such as local atomic density and the SOAP kernel.⁷⁴ Another kernel-based method is gradient domain machine learning (GDML), which is trained to predict the force directly rather than the energy.^75–78 For artificial NNs, high-dimensional neural network potential (HDNNP) was the first descriptor-based neural network potential (NNP) that uses atomic-centered symmetry functions as descriptors,^79,80 and subsequent developments included the consideration of long-range interaction and non-local phenomena.^81–83 Some different NNPs were also developed based on different descriptors such as ANI,⁸⁴ TensorMol,⁸⁵ deep potential,⁸⁶ QM/MM-NN,^87,88 water NN,⁸⁹ and embedded atom NN.^90,91 Another type of NN is end-to-end NN, which directly uses the Cartesian coordinates and nuclear charges as inputs instead of descriptors.⁹² End-to-end NN is usually based on graph NN,⁹³ like message-passing NN,⁹⁴ and some prominent examples are deep tensor NN,^95,96 SchNet,^97,98 and PhysNet.^99,100 Machine learning provides a powerful and promising way to develop force fields with QM accuracy and excellent efficiency.

Despite the success of machine learning force fields in a plethora of chemical problems, including electronic effects, thermodynamics, reactions, and spectroscopies,¹⁰¹ the development of general, transferable, and scalable force fields for proteins is still very challenging for the following reasons. First, to obtain plenty of data for the training of a machine learning model, the calculation of QM reference energy is infeasible for proteins since they could have thousands of atoms. Second, the size of proteins for MD simulations could range from short peptides to large proteins, which increases the difficulty of training a general machine learning force field that is applicable to different systems. To address these issues, fragmentation methods can be utilized to systematically express the total energy of large systems with the individual energy of small fragments.¹⁰² Some common fragmentation methods include divide and conquer,¹⁰³ many-body expansion,^104–107 systematic molecular fragmentation (SMF),^108–110 the effective fragment potential approach,¹¹¹ electrostatically embedded generalized molecular fractionation with conjugate caps (EE-GMFCC),^112–115 the X-pol method,^116–119 and the energy-based fragmentation method.^120,121 In 2019, Hao and Yang developed residue-based systematic molecular fragmentation (rSMF) to enable the construction of general fragments for any proteins, and the neural network force field (NNFF) was trained on glycine and alanine dipeptides to demonstrate the transferability to mixed alanine and glycine polypeptides.¹²² More machine learned force fields based on similar fragmentation methods were developed,^123–125 but a comprehensive and transferable force field that can be applied to general protein simulations is still not available.

In this work, we apply the rSMF method to obtain fragments that are trained with NN and develop an NNFF based on the CHARMM force field, named CHARMM-NN. The CHARMM-NN is based on atom types, and the input features are similar to the variables in CHARMM energy terms, such as bonds, angles, dihedrals, and non-bonded features, which is highly compatible with the original MM simulations. The CHARMM-NN is first validated by the geometric data and structural reorganization energies, and more simulations of peptides and proteins are performed to test the quality of this force field.

We approach the construction of the protein force field by partitioning the protein into embedded fragments and developing the NN description of fragments with input data from QM calculations.¹²² Based on the original SMF method,^108,126 residue-based systematic molecular fragmentation (rSMF) was developed by Wang and Yang in 2019,¹²² which generates fragments that can be applied to any protein. The fragments are based on amino acid residues, and different levels of rSMF can be applied. The rSMF at level 1 generates dipeptides, defined by one side chain group and two peptide bonds, as the basic fragments. A higher level of rSMF can be applied to increase the accuracy of fragmentation, but the sampling effort as well as the computation of QM energies and forces will increase significantly. Therefore, the rSMF used in this work is at level 1. Using the rSMF method, a protein with N residues capping with an acetyl (ACE) group on the N terminus and an N-methyl amide (NME) group on the C terminus can be divided into N dipeptides and N − 1 ACE-NME fragments. As an example, a tripeptide P³ is fragmented, as shown in Fig. 1, where A₁ and A₂ are the corresponding dipeptides for each amino acid in the tripeptide, and ACE-NME is the peptide bond between A₁ and A₂.

The test of the rSMF method on homogeneous and heterogeneous alanine and glycine polypeptides suggests that this method can provide stability of errors with increasing size, which is highly desirable and necessary in protein force field development since the protein size could be arbitrarily large. Based on this behavior, we construct NN models for the ACE-NME fragment and 24 dipeptide fragments, including three dipeptides for different protonation states of His (Hse, Hsd, and Hsp), two protonated states of negative charged amino acid dipeptides (AspH and GluH), and 19 other natural amino acid dipeptides (Gly, Ala, Val, Cys, Pro, Leu, Ile, Met, Trp, Phe, Ser, Thr, Tyr, Asn, Gln, Lys, Arg, Asp, and Glu). Even though these fragments are enough for some proteins, many proteins contain disulfide bonds formed by two cysteines, which cannot be described by the fragments. Therefore, to enable the applicability on more proteins, we complement the rSMF method with two more fragments: one is the dipeptide fragment of cysteine with an SMe group connected to the S atom, named Cys-SMe, and another is the dimethyl disulfide (DMDS) fragment that accounts for the overlap between two cysteines that form a disulfide bond. As shown in Fig. 2, the molecule of two cysteines connected by a disulfide bond can be fragmented into two Cys-SMe dipeptides and a DMDS fragment. The CHARMM-NN includes 27 fragments in total, and each fragment has its own NN model.

The NN correction energy is the sum of the energies of all dipeptide fragments minus the energies from ACE-NME and DMDS, and the total energy is the sum of the MM energy and the NN correction energy. This is an approximation because the non-bonded interactions between the capped hydrogen atoms and other atoms in the dipeptide fragment and ACE-NME fragment are not completely the same, even though the additional bond, angle, and dihedral energy terms that are related to the capped hydrogen atoms can be exactly canceled each other. Based on this approximation, we can run normal MD simulations with C36m force fields and apply the additional correction from NNFF, which is very convenient for integration into existing MD programs without the requirement to modify core codes.

All MD simulations to test CHARMM-NN on dipeptides, peptides, and proteins were performed in the modified GROMACS 2021.6,¹⁴⁴ in which the NN corrections were calculated with the TensorFlow C API, and all QM energies and forces were calculated by Gaussian 16.¹⁴⁵ The basic MD parameters, such as cutoff scheme, temperature coupling, and pressure coupling, are the same as they were in the adaptive samplings. In the simulations to validate CHARMM-NN, the bonds involving hydrogens were constrained by the LINCS algorithm¹⁴⁶ with a fourth order and one iteration. The integration time step was, therefore, set at 1 or 2 fs. All systems were first solvated by creating a cubic water box with the length of the longest distance in the protein plus 2 nm buffer and neutralizing it with enough sodium or chloride counter-ions. Next, the solvated systems were minimized with the steepest descent algorithm and equilibrium with a 500 ps NVT simulation using a v-rescale thermostat followed by a 1 ns NPT simulation using a Parrinello–Rahman barostat. Starting from the QM-optimized structures, the conformers of all dipeptides were minimized with the steepest descent algorithm for up to 50 000 steps for C36m and CHARMM-NN. The structural reorganization energies were calculated by subtracting the energy at the last step of minimization, which is the MM-optimized structure, from the energy at the first step of minimization, which is the QM-optimized structure. All production MD simulations were performed in solution with the TIP3P model, and a v-rescale thermostat and Parrinello–Rahman barostat were used. For the simulations of Gly₃, Ala₃, Val₃, and Ala₄, the total simulation time was set as 1 µs, in which the first 100 ns were discarded and the coordinates were saved every 10 ps. The J-coupling constants were calculated with Karplus equations, and the parameters were extracted following the previous work.¹⁴⁷ Slightly different results would be obtained if different sets of Karplus parameters are applied.

All data related to this work can be found at https://figshare.com/projects/CHARMM-NN/163756. The energy, force, and coordinates of the sampled structures, the CHARMM-NN force field parameters, the codes to run MD simulations with CHARMM-NN, training data and scripts, MD simulation files, and analysis scripts are provided.

Next, we examined the detailed performance of CHARMM-NN on dipeptides. One typical system is alanine dipeptide, which has many conformations and has received systematic studies.¹⁴⁸ The graphical structure of Ala dipeptide is shown in Fig. 3. We first compare the geometric information optimized from C36m, CHARMM-NN, and B3LYP. As displayed in Table II, the CHARMM-NN outperforms C36m on all three Ala conformers, C7_eq, C7_ax, and C5. Most angles and dihedral angles from CHARMM-NN are closer to the B3LYP geometries compared to C36m, and almost all bonds optimized from CHARMM-NN are better than C36m. The large corrections occur on C4–C5 and C5–N7 bonds, for which the bond lengths from C36m are generally much smaller than B3LYP structures, and on the C9–C11 bond, which typically has a large bond length from C36m.

Besides the geometric evaluation, we also compared the relative potential energies and reorganization energies between the C36m and CHARMM-NN force fields. As shown in Table III, six Ala conformers, C7_eq, C7_ax, C5, α_L, β₂, and α′, were calculated because the optimization with B3LYP can be obtained for these conformers so that direct comparison can be executed. The optimized backbone dihedrals from QM calculations are shown in Table III for all conformers. According to B3LYP calculations, the C7_eq conformer has the lowest energy, which is selected as the baseline for comparison, and the relative energies are listed in Table III, in which the C5 and α′ conformers have the lowest and highest relative energies, respectively. For the conformers minimized with C36m, the errors of relative energies are between 0.8 and 1.4 kcal/mol for C7_ax, α_L, and α′ conformers compared to B3LYP. For C5 and β₂ conformers, the errors are up to 2.29 and 1.92 kcal/mol, respectively, and the C36m optimized C5 tends to be slightly more stable than C7_eq. Using CHARMM-NN, the relative energies are very consistent with B3LYP results for all conformers, and the largest error is only 0.2 kcal/mol. The accurate relative energies from CHARMM-NN demonstrate that the local minima on the PES of Ala dipeptide are close to QM PES, which can be very important for the simulation samplings being similar to QM results. In addition to the relative potential energies, we also compared the reorganization energies between C36m and CHARMM-NN for these six conformers of Ala dipeptide. In Table III, we can observe that the reorganization energies with respect to B3LYP range from 12.63 to 31.55 kJ/mol for C36m but only between 0.93 and 2.39 kJ/mol for CHARMM-NN. The minimized structures from CHARMM-NN are much closer to B3LYP optimal structures than C36m.

Besides the few tests of reorganization energies for Ala dipeptides with respect to B3LYP, which is the reference QM method that we train the CHARMM-NN on, more thorough validations of the reorganization energies were performed on the YMPJ conformer database,^149,150 which includes the conformers of dipeptides optimized with MP2/cc-PVTZ. As shown in Table IV, the 14 neutral dipeptides have different numbers of conformers at energy minima. For C36m, the reorganization energies for almost all dipeptides are larger than 30 kJ/mol, except for the Gly dipeptide (24.33 kJ/mol), which is the smallest dipeptide. The reorganization energies can be higher than 40 kJ/mol for Gln, Met, Thr, and Tyr dipeptides. For CHARMM-NN, the reorganization energies for most dipeptides are below 10 kJ/mol, except for Gln, Leu, Met, and Phe dipeptides. The Met dipeptide, which includes the most conformers, has the largest reorganization energy (12.07 kJ/mol), but it is still much better than the C36m results. The average reorganization energy for C36m is 38.82 kJ/mol, which is close to the work from König and Riniker.⁶³ Even though the CHARMM-NN is trained with B3LYP, it outperforms C36m on the dataset of MP2/cc-PVTZ as well, and the average reorganization energy is only 9.01 kJ/mol. The low reorganization energy with respect to QM is critical to obtaining accurate thermodynamic properties like free energy that are close to QM since the expected variance is exponentially related to the reorganization energies for the free-energy estimation.⁶³ 

To validate the performance of CHARMM-NN in MD simulations, the J-coupling constants were calculated from the MD simulations for Gly₃, Ala₃, Val₃, and Ala₄, and the results are shown in Table V. For Gly₃, the J-coupling constants are close to the experimental values for residues 2 and 3. The largest error is on ²J(N, C_α), which has a deviation of 2.66 and 1.38 Hz from experiments for residue 2 and residue 3, respectively. However, the experimental values of 10.45 and 9.05 Hz can never be obtained because the maximum value is less than 9 Hz with the parameters for the Karplus equation. For Ala₃, the deviations for ³J(H_N, H_α), ³J(H_α, C′), and ³J(H_N, C_β) are around 0.82 to 1.01 Hz for residue 2, but they are much better for residue 3. Other J-coupling constants have similar errors for residue 2 and residue 3, and the largest deviation is 2.00 Hz on ²J(N, C_α) for residue 2. For Val₃, the ³J(H_N, H_α) is pretty good for residue 2 with a difference of 0.40 Hz, but it is terrible for residue 3, which has an error of 4.32 Hz. In contrast, the ¹J(N, C_α) and ²J(N, C_α) have a better agreement with experiments on residue 3 compared to residue 2. For Ala₄, the results are similar to Ala₃ in that the ³J(H_N, H_α), ³J(H_α, C′), and ³J(H_N, C_β) have less errors from residue 2 to residue 4. Overall, Gly₃ has good results for all J-coupling constants except ²J(N, C_α). The Ala₃ and Ala₄ have similar results, suggesting that the simulations on systems with different sizes provide stable results for CHARMM-NN. The results of Val₃ are good for some J-coupling constants, but they can also have a large error on others. In addition, MM force fields still have better agreement with experimental J-coupling values than the current CHARMM-NN, mainly because it only has corrections for intramolecular interactions of dipeptides, while the protein–protein and protein–solvent non-covalent interactions have not been addressed yet.

We also ran several simulations on several folded proteins with CHARMM-NN. The folded proteins we tested are ubiquitin (PDB: 1UBQ¹⁵¹), crambin (PDB: 1EJG¹⁵²), GB3 (PDB: 1P7E¹⁵³), and lysozyme (PDB: 135L¹⁵⁴). As displayed in Fig. 4, the potential of the mean force plots for the folded proteins is similar to the normal Ramachandran plots. The dominating regions include α region (−160° < ϕ < −20° and −120° < ψ < 50°) in which the right-handed α helix α_R (−100° < ϕ < −30° and −67° < ψ < −7°) are sampled most, β region (−180° < ϕ < −90° and 50° < ψ < 180° plus −180° < ϕ < −90° and −180° < ψ < −120° plus 160° < ϕ < 180° and 110° < ψ < 180°), ppII region (−90° < ϕ < −20° and 50° < ψ < 180° plus −90° < ϕ < −20° and −180° < ψ < −120°), and α_L region (30° < ϕ < 100° and 7° < ψ < 67°).¹⁸ The crambin and GB3 have almost all their points located in these regions, showing the reliability of CHARMM-NN in sampling conformational space. However, results for ubiquitin and lysozyme show several unfavored regions, such as the regions with 100° < ϕ < 120° for ubiquitin and the regions with 0° < ϕ < 30° for lysozyme. The non-ideal results from the Ramachandran plots of some folded proteins and the J-coupling constants of some peptides indicate that the CHARMM-NN is still not ready to be used for general proteins. The difficulty of CHARMM-NN is because of the following reasons. First, it is only trained with all data in the vacuum state, and no background charges are included in QM calculations and NN inputs; therefore, the method can only approach the limit of accuracy at the mechanical embedding QM/MM level without considering the MM charge contribution to the QM region. We do include solvent in all MD simulations of peptides and proteins, while the protein–solvent interactions are based on the MM force fields and not improved with machine learning models. Protein–water interactions are crucial to MD simulations of proteins. However, the protein–water interactions cannot be easily included in the training of NN because the solvent environment of small fragments is largely different from the solvent environment of proteins, and there is yet no general sampling method that can be applied to sample enough small fragments surrounding with solvent to resemble the actual solvent environment in all proteins. Second, the non-bonded interactions between fragments are still calculated at the MM level, which may cause the CHARMM-NN corrections to be imbalanced on proteins.

In summary, we constructed machine learned force fields based on C36m force fields, named CHARMM-NN, by using the rSMF method to generate the elementary fragments that can form any kind of protein. The dataset was generated and enlarged by combining NMS and adaptive sampling based on prediction uncertainty. The CHARMM-NN force fields use the atom-type based NN to calculate energies and forces, and the input features can be obtained or simply derived from the traditional MM variables. The high compatibility between CHARMM-NN and MM force fields enables the convenient implementation of CHARMM-NN in all MD programs without the need to modify core codes. The training error for CHARMM-NN is low for all 27 fragments, and the validations on dipeptides demonstrate that CHARMM-NN can result in good geometric data similar to QM calculations and much lower reorganization energies than traditional MM force fields. For the MD simulations with CHARMM-NN on several peptides and proteins, some results are acceptable, but some results are not sufficiently accurate because the current CHARMM-NN can only achieve corrections at the mechanical embedding QM/MM level, suggesting that the CHARMM-NN needs to be improved further. Future directions include considering the comprehensive solvent effect in data sampling and constructing machine learned force fields to describe non-bonded interactions between fragments beyond the mechanical embedding level.

This work has been supported by the National Institutes of Health (Grant No. R01-GM061870).

The authors have no conflicts to disclose.

Pan Zhang: Data curation (lead); Formal analysis (lead); Investigation (lead); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead). Weitao Yang: Conceptualization (lead); Funding acquisition (lead); Methodology (lead); Project administration (lead); Resources (lead); Supervision (lead); Writing – review & editing (lead).

The data that support the findings of this study are available within the article.