Assessing the capability of in silico mutation protocols for predicting the finite temperature conformation of amino acids †

Mutation protocols are a key tool in computational biophysics for modelling unknown side chain conformations. In particular, these protocols are used to generate the starting structures for molecular dynamics simulations. The accuracy of the initial side chain and backbone placement is crucial to obtain a stable and quickly converging simulation. In this work, we assessed the performance of several mutation protocols in predicting the most probable conformer observed in finite temperature molecular dynamics simulations for a set of protein–peptide crystals differing only by single-point mutations in the peptide sequence. Our results show that several programs which predict well the crystal conformations fail to predict the most probable finite temperature configuration. Methods relying on backbone-dependent rotamer libraries have, in general, a better performance, but even the best protocol fails in predicting approximately 30% of the mutations.


Introduction
Mutational protocols have been used successfully for predicting structures of proteins that were not resolved experimentally, 1 for determining the interfaces in protein-protein interactions, 2,3 and for designing de novo peptides. [4][5][6] One widespread application is the computational scanning of alanine or glycine residues, in order to identify hot spots and key amino acids responsible of the protein-protein stabilizing interactions. 7,8 The placement of the mutated residue is crucial to understand the potential effects of the mutation. 9,10 Most mutation protocols require the backbone or C-alpha position of the amino acid, and then generate a side chain conformation. Some protocols, for example Rosetta fixbb 11 and SCWRL4 12 use rotamer libraries from public databases of experimentally-resolved protein structures to predict the side chain configuration. Other protocols use minimization approaches to find the optimal conformation that minimizes an empirical scoring function. 13 Combinatorial approaches have been developed to mutate multiple or single amino acids. 14,15 To improve the accuracy of the side chain prediction, a fundamental step is to sample the conformations of the system. This can be performed using stochastic methods such as Monte Carlo, based on movements constrained by dihedrals, 16 and classical or enhanced molecular dynamics (MD) simulations. 17 However, sampling the rotamer space of the amino acid is time consuming. In protein design applications, where iterative single-point mutations protocols are required for designing novel molecules, running long MD simulations for each mutation is computationally expensive because the number of mutations that have to be explored, even for a small peptide, is extremely large. 18 Protocols able to predict rotamers correctly can diminish the required simulation time to explore the conformational space.
In addition, starting an MD simulation from the wrong rotamer can be problematic because it can lead to conformational changes affecting events such as protein folding 19 and binding. 20 For example, in some cases, the folded conformations of the Trp-cage domain have been simulated with low success using classic and enhanced MD. This has been associated to the erroneous placement of near-native rotamers of the central tryptophan side chain. 21 In the study of proteinprotein interaction, starting the simulation with key residues at the interface in the wrong conformation can be detrimental. In ESI, † we provided an example involving the protein-protein complex Barnase-Barstar, 22 where starting from a single wrong rotamer (compared to the crystal structure) causes a considerable loss of native contacts which are not retrieved along all the MD simulation (see Fig. S1, ESI †). Such effect can, for example, impact the prediction of binding affinities. Therefore, an optimal mutational protocol should predict side chain conformations that are in agreement with the equilibrium distributions from finite temperature simulations.
Here, we assess the performance of different mutation protocols using as a benchmark MD simulations of a set of protein-peptide complexes that differ only by single-point mutations in the peptide sequence. We compare the side chain dihedral angles, the number of contacts and the number of hydrogen bonds resulting from the mutation protocols to the equilibrium distribution of those quantities. The results suggest a rational pipeline to improve the mutational protocols for efficient MD simulations, and protein or peptide design.

Single-point mutation benchmark systems
A set of crystal protein structures of the Oligopeptide Binding Protein (OppA) interacting with a tripeptide of motif K-x-K 23 (where x is one out of 11 amino acids) was used as a reference system (Fig. 1A) to test a set of single-point mutation protocols. The structures were obtained from the Protein Data Bank (PDB) 24 (PDB ids: 1B4Z, 1B51, 1QKA, 1B3G, 1B32, 1B3F, 1B46, 1B5J, 1QKB, 1B5I, 1B40) and formatted to share the same amino acid number and chain identifiers.
We also assessed the protocols with other protein-peptide systems. Two complexes, which had structures differing only by a single amino acid on the peptide sequence, were selected: the HLA class I antigen A-2 alpha chain 25 and the MDM4 protein complexed with a 12-mer peptide. 26 For the HLA class I complex, four crystals were selected from the PDB differing only by a single peptide amino acid at the 5th position (PDB ids: 3GSO, 3GSU, 3GSV, 3GSR). For MDM4-peptide complex, the variable amino acid is at the peptide's 6th position (PDB ids: 3JZP, 3JZO).

Mutation scheme
A combinatorial approach was used to construct the mutants. This consists of mutating the variable amino acid to all the other amino acids available in the crystals. In Fig. 1B, we show the mutation scheme for the OppA system, where the variable amino acid from the tripeptide of each crystal is mutated to all other 10 amino acids. Thus, each structure produces 10 different mutations, resulting in 10 different structures for the same peptide sequence. This strategy allowed us to characterize the possible impact of steric and volume constraints due to the starting structure.

Molecular dynamics simulations
Each protein-peptide complex from the benchmark was submitted to 20 nanoseconds (ns) of MD simulations with previous minimization and NVT/NPT equilibration phases. The system was minimized using the steepest descent algorithm, with 50 000 steps and a maximum force threshold of 10 kJ mol À1 nm À1 . NVT and NPT equilibrations were performed for 100 ps using position restraints on the heavy atoms of the protein to allow for the equilibration of the solvent. GROMACS v5.1 27 was used to perform the MD simulations. The Amber99SB-ILDN protein force-field 28 and TIP3P water model 29 were used. The protein was solvated with a cubic box of water with a distance of 8 Å from the furthest atom of the protein. After solvation, counterions of Na + and Cl À were included in the solvent to make the box neutral. The simulation was run using a modified Berendsen thermostat 30 at 330 K temperature-coupling, and the Parrinello-Rahman barostat. 31 The electrostatic interactions were calculated using the Particle Mesh Ewald (PME) method with 1.0 nm short-range electrostatic and van der Waals cutoffs. 32 The equations of motion were solved with the leap-frog integrator 33 using a timestep of 2 fs.
The convergence of the simulations was monitored by computing several observables: the number of hydrogen bonds between the peptide and protein, the number of heavy atom contacts 34 made by the mutated amino acid, and the all-atom Root-Mean-Square Deviation (RMSD) of the peptide from the reference crystal. In some cases, the simulations were extended to achieve convergence.

Mutation protocols
We analyzed five mutation protocols. Modeller v9.19, 35 which relies on minimization cycles of rotamers derived from homologybased models. SCWRL4 36 and Rosetta fixbb v2017.26, 11 which use backbone-dependent rotamer libraries to predict side chain conformations but with different scoring functions. TLEaP from AmberTools v16, 37 used to generate topology files from protein structures for the Amber program. 37 FoldX v4 38 a protocol implemented in protein folding simulations that depends on energy calculations derived from an empirical force field, used also to study the effects of point mutations.
These programs were configured to generate single-point mutations of the variable amino acids. The mutations were performed starting from the crystal structure and also from the structure obtained in the last frame of the MD trajectories (see ESI †). The prediction of the mutated side chain was made, for all protocols, with the peptide in complex with the protein.

Evaluation criteria
To evaluate the performance of the mutation protocols, we compare the dihedral angles and number of contacts predicted by the protocol to the most probable conformations from the MD simulations. The evaluation criteria are exemplified in Fig. 2. The conformational ensemble of protein folds 39 and protein-protein complexes 40,41 obtained from equilibrated finite temperature simulations has proved to depict accurately the interactions between the involved molecules.
Backbone dihedrals (f, c) and side chain dihedrals w 1 and w 2 were calculated for each mutated amino acid. For the side chain dihedrals, the distribution from the MD simulations was used to verify that the predicted rotamers are located within the conformations visited in the MD. The w 1 ([0,360] deg) and w 2 ([0,360] deg) 2D dihedral angle space was binned using a 50 Â 50 grid. Thus, along each dihedral direction a bin size of 7.2 deg was taken, resulting in a 7.2 Â 7.2 deg 2 for the 2D bin size. To indicate if a side chain prediction is visited in the trajectory, we looked at the 2D bin corresponding to the rotamer prediction and compared its population to that of the most populated bin. We count a rotamer as visited during the MD trajectory if the rotamer falls within a bin which has an MD population that is at least 5% of that of the most populated bin. We assume that if the mutated conformation is sampled in the trajectory then it is possible to reach equilibrium in a similar or shorter time. The threshold allows us to compute a success rate for the mutation protocols that can assess which protocols perform better.
Another evaluation criterion was defined using the heavyatom contacts between the mutated amino acid and the protein. For this purpose, we used the information from the bins previously established in the side chain dihedral analysis. Specifically, each mutated complex, which is located in a bin from the w 1 and w 2 grid, was compared to the MD structures belonging to the same bin. The heavy atom contacts were monitored using the contact matrix 42 with different distance thresholds d 0 , including 4 Å, 3.5 Å and 3 Å. The results of the latter are included in the main text, and the others are available in the ESI. † The contact matrix is defined as: where R ij is the distance between atom i and j. The number of conserved contacts (S MR ) between the mutated complex (C M ) and a reference structure from the MD (C R ), was estimated as, where the i sum runs over the peptide atoms and the j sum runs over the protein atoms. The result is a number from 0 to 1, where 1 is the most successful scenario (i.e., all the contacts predicted by the mutation are also present in the MD structure). The average and standard deviation of the S MR for each protocol and each mutated amino acid were obtained by averaging over five structures from the corresponding dihedral bin.

Rosetta mutation-protocol modifications
The Rosetta Commons project (www.rosettacommons.org) has available Monte Carlo approaches to optimize both backbone and side chain dihedrals of protein structures. These methods can be used to refine the system before or after performing a single-point mutation. We implemented the following protocols: relaxing with rigid backbone, 43 prepacking of interface side chains, 44 refinement of the system with FlexPepDock 45 and inclusion of backbone flexibility for both protein and peptide using the BackRub protocol. 9 The differences between these protocols are related to the type of molecular movements, the computational exhaustiveness and the internal constraints used to obtain the lowest energy conformations of the selected amino acids. The performance of these modifications was also evaluated using the previously described criteria.

Results
We tested various protocols to perform single point mutations on peptides bound to protein targets. We used all-atom MD simulations with explicit solvent to sample the conformational space and compare to the results from the mutation protocols. First, we present the results for the OppA-tripeptide complex using the evaluation criteria, and then we show the results for the HLA class I and the MDM4 complexes (see Methods).

Convergence of the molecular dynamics simulations
We used the equilibrium ensemble from MD simulations as a test to evaluate the performance of the mutation protocols. The benchmark system of the 11 OppA crystals presented stable observables, such as the number of hydrogen bonds and all-atom RMSD, during the 20 ns of simulation (see Methods and Fig. S2, ESI †). We used the complete 20 ns trajectory to calculate the equilibrium distributions. We tracked the distribution of the backbone and side chain dihedrals angles of the mutated amino acid during the MD trajectories. We found that the backbone dihedrals remain quite stable during all the MD trajectories (see Fig. S3, ESI †). For the side chain dihedrals (w 1 and w 2 ), we checked that the most frequent conformations are categorized in three on-rotamer groups: gauche(+), trans and gauche(À), centered in 300, 1801 and 601 respectively. 46 In general, for both w 1 and w 2 most of the side chain conformations were classified as gauche(+), which is in fact the most abundant conformation in the PDB, with the gamma side chain pointing in an opposite direction to the main chain nitrogen (see Fig. S4, ESI †). 47

Side chain dihedral prediction
After calculating the distributions from the converged trajectories, we compared the predicted w 1 and w 2 dihedrals of each mutation protocol to the distribution from the MD simulations (see Methods for details). We first evaluated the prediction of w 1 dihedral for all the 11 amino acids. The results for serine and valine are presented in Fig. 3. We also report for each amino acid the percentage of success in predicting the dihedral with a conformation located in a histogram bin with a population larger than 5% of the most populated rotamer bin (Table 1). If only w 1 is considered, Rosetta fixbb and FoldX are the methods with highest success rates.
To better assess the side chain dihedral prediction, we performed a similar analysis taking into account both the w 1 and w 2 dihedrals. In Fig. 4, we show the results for the mutated amino acids isoleucine (I), arginine (R) and asparagine (N) (see Fig. S5 and S6 for other amino acids, ESI †). In Table 2, we report the percentage of mutations that succeeded, for each protocol, in predicting both the w 1 and w 2 dihedrals using a 5% bin-threshold (see Methods). A similar table using a 30% bin-threshold was calculated (Table S1, ESI †), which despite being stricter over the conformational space shows the same tendencies as reported in Table 2. The results show that SCWRL4, Rosetta fixbb and FoldX (see Methods) are the most successful protocols for predicting both the w 1 and w 2 dihedral angles. However, the success rate decreases for all protocols in comparison to the results for only w 1 . The good performance of SCWRL4 and Rosetta fixbb can be related to the fact that both use backbone-dependent rotamer libraries as the basis for rotamer selection.
A similar analysis was also performed for mutations created from the last MD frame (instead of using the crystal structure), and we observed similar trends for the w 1 and w 2 dihedrals prediction (see Table S2, ESI †).
In addition, we analyzed the performance from the perspective of the starting amino acid that was after mutated. This analysis elucidates if the starting amino acid produces an effect, e.g., due to its size or side-chain orientation, over the mutation protocol performance. We find that for the OppA crystals, the results for each protocol are similar for all the amino acids (see Table S3, ESI †). This can be explained by the side chain orientation, which is always pointing to the same direction (Fig. 1A). This implies that the starting amino acid has small impact on the mutation protocols.

Conservation of contacts
For the second evaluation criteria, we calculated the average and standard deviation of the conserved contacts for each selected protocol and each predicted amino acid (Methods). For each variable amino acid, the comparison was made between the 10 mutated structures of the OppA complex and five selected structures from the trajectory with similar dihedral angles (see Methods for details), giving a total of 50 comparisons per amino acid per protocol. The results using a contact threshold of 3 Å are shown in Table 3 (see Tables S4 and S5 for other thresholds, ESI †). Similarly as for the dihedral analysis, we found that methods using backbone-dependent rotamer libraries are better in predicting the contacts observed in the MD simulations, demonstrating the capabilities of the Rosetta fixbb and SCWRL packages to perform single point-mutations in peptides bound to proteins.

Performance of the modifications to the Rosetta mutation protocol
Based on the previous analysis, Rosetta fixbb has one of the best performances. We studied if it can be improved based on some available protocols 9,43-45 to move the backbone and side  chains atoms. We tested five Rosetta protocols. One protocol (pre-Relax) was applied before the mutation, and the other four were made after the mutation. The results of the dihedral analysis and contact conservation using these protocols are described in Table 4. These results indicate that relaxing the complex with a fixed (Post-Relax) or flexible (BackRub) backbone, after the mutation, slightly improves both the dihedral rotamer prediction and the conservation of contacts. The protocols designed for docking (Pre-Pack and FlexPepDock) were not able to improve the rotamer prediction.
Both systems were submitted to single-point mutations using a combinatorial scheme similar to that described in the methods but according to the number of available structures. Using the same assessment as for the OppA system, in Table 5, we describe for each mutation if the mutation protocol was able to predict both the w 1 and w 2 dihedrals correctly. For the case of the MDM4 complex, all protocols, with the exception of TLEaP, were capable to predict the rotamers explored by the MD trajectory (see Fig. S8, ESI †). For the HLA class I the performance was different. The best performing protocols were Modeller, SCWRL4 and Rosetta fixbb after relaxing the complex. However, none of the protocols were able to put the methionine side chains as explored by the MD simulations. This might be influenced by the intrinsic flexibility of the amino acids that are not creating hydrogen bonds with the pockets of the HLA class I a chain. 48

Discussion
For comparing the performance of the mutational protocols it was necessary to use a robust benchmark system with various crystals differing only by a single mutation. We found the OppA system ideal for these purposes due to the wide availability of crystal structures. Moreover, working with tripeptides allowed the system to reach equilibrium in a short computational time (B20 ns), which is not guaranteed for other systems where proteins are usually longer. 49 We analyzed two additional systems to assess the performance with longer and structurally different peptide chains. We found that the mutation protocol performances are consistent for all test sets.
By comparing the predicted side chain dihedral angles and contacts to the equilibrium distributions from finite temperature simulations, we found that the protocols based on rotamer libraries derived from protein structures 50 perform better. In this context, SCWRL4 and Rosetta fixbb are suitable methods for performing single-point mutations. Both mutation protocols use backbone-dependent rotamer libraries but differ in how each rotamer is scored. SCWRL4 uses a scoring scheme based on single and pairwise rotamer energies computed from attractive and repulsive hydrogen bond and van der Waals terms. 51 Rosetta, in addition, includes weighted terms related to statistical energies derived from distance-dependent pair potentials and solvation energies. 52 The results for other mutation protocols were in some cases comparable or better than Rosetta and SCWRL4, but on average their performance was worst. One of the closest in performance to SCWRL4 and Rosetta fixbb is Modeller, which relies on cycles of minimization with the inclusion of homology-derived dihedral angle restraints, but without using rotamer libraries. 35 This may be a factor that hinders its performance. TLEaP, the worst performing program, is used by AMBER to generate topology files for MD simulations, with the option to add missing atoms based on force field information. 53 The main issue with TLEaP is that most of the dihedrals are predicted as trans/trans conformation without optimizating the dihedrals based on environment interactions. For FoldX, 54 we obtained variable performances, showing successful results in cases where other protocols were unable to predict correctly the conformation, for example for aspartic acid. However, it is important to remark that neither the protocols that use rotamer-libraries, nor those that do not, are able to consistently predict the side chain conformations sampled in the simulations.
Amino acids such as aspartic acid and asparagine were more complex to predict, possibly because the hydrogen bonding capacity of their side chains, which impacts the w 2 dihedral. 55 Another aspect to take into account is the chosen force field. The Amberff99SB-ILDN force field has been parameterized to explore the conformational space of proteins in long MD simulations, using quantum mechanics and experimental parameter data. 28 Interestingly, the equilibrium distributions for aspartic acid and asparagine for this force field differed the most from expectations based on the Protein Data Bank statistics. 24 Computational time is also a relevant consideration for selecting a mutation protocol. Both the time required to predict the rotamer and the computational time of the MD simulation to reach a stable conformation are important variables. All five mutation protocols tested are able to predict the mutated side chain in just a few seconds. When refinements are added to the Rosetta fixbb prediction, the computational time increases to a few minutes, with the exception of the BackRub method that is approximately 10 times longer than the standard Rosetta fixbb protocol. For the case of the MD simulation time, the goal of this work is to select the protocols able to predict the most probable rotamers explored by MD simulations. Consequently, the best performing mutation protocols will reduce the computational sampling time, and consequently the time required to reach a stable conformation. Importantly, we note that very wrong rotamer predictions could destabilize so much the structure that it might be very difficult to obtain a converged MD simulation. Finally, we found that refining the protein-peptide complex after the mutation can improve the side chain prediction. 56 The modularity of the Rosetta protocols to perform the refinement was useful at this scope. We found that relaxing the side chain without having to move the backbone improves the performance whilst maintaining a reasonable computational time.

Conclusions
The assessment of the single-point mutation protocols to predict side chains from equilibrium distributions has shown that, although some protocols are able to predict the most probable rotamer explored in MD, there is still large room for improvements. This is of key importance to the MD community, which highly relies on homology modelling and rotamer prediction. In addition, these protocols are also essential for peptide design, where filtering mutations in a random or guided way can contribute dramatically to the design protocol efficiency. 57 Our work sets a basis to assess, and to further improve and optimize, the mutation protocols to be in accordance with finite temperature simulations. Previously, all strategies have been optimized to predict crystal-structure rotamers. We propose to use also MD simulations as training sets to reparameterize and optimize the mutation algorithms.

Conflicts of interest
There are no conflicts to declare.