Predicting Protein Dimer Structures Using MELD \xc3\x97 MD
It is challenging to predict the docked conformations of two proteins. Current methods are susceptible to errors from treating proteins as rigid bodies and from an inability to compute relative Boltzmann populations of different docked conformations. Here, we show that by using the ClusPro server as a front end to generate possible protein\xe2\x80\x93protein contacts, and using Modeling Employing Limited Data (MELD) accelerated molecular dynamics (MELD \xc3\x97 MD) as a back end for atomistic simulations, we can find 16/20 native dimer structures of small proteins as those having the lowest free energy, starting from good\xe2\x80\x93bound\xe2\x80\x93backbone structures. We show that atomistic MD free energies can be used to identify native protein dimer structures.
INTRODUCTION<!>MELD Accelerates Conformational Exploration by Using External Information.<!>MELD \xc3\x97 MD Setup.<!>MD Simulation Parameters.<!>Dimer Selection and System Flexibility.<!>RESULTS<!>Using Information from CP, MELD \xc3\x97 MD Successfully Docks 16/20 Dimers.<!>MELD \xe2\x80\x9cPicks\xe2\x80\x9d Out Which ClusPro Contacts Are Native.<!>MELD \xc3\x97 MD Often Gives Better Structures than Vanilla CP.<!>Failures.<!>We Can Apply Our Protocol to an Example CASP/CAPRI Case.<!>CONCLUSIONS
<p>We describe a way to improve computer predictions of the structures of protein–protein dimers. Much work has been done in modeling protein docking. The protein docking community has established an event, called CAPRI (Critical Assessment of PRediction of Interactions), for the blind testing of docking algorithms.1–3 This has brought the challenges into clear focus. Proteins have so many degrees of freedom of internal and relative motion that simplifications are needed to tame the computational combinatorial explosion.4–11 One main simplification is to treat each protein as a rigid body. Another simplification is to use quasi-physical "scoring functions" instead of more physically accurate free energies. Both introduce errors into structure prediction.12,13 Recently the ProPOSE approach has been developed to deal with some of these limitations allowing for flexible side chains during the docking prediction.14 Molecular dynamics (MD) simulations with atomistic force fields in principle represent an approach to such limitations, but at a extremely high computational cost. For example, using highly specialized computational resources,15 Shaw helped to refine oligomeric structures using MD16 and observed several undocking–docking events of five protein dimers starting from the bound structure.17 To reduce the computational cost of MD, Hou et al. used simplified coarse grained molecular models to quite successfully compute the binding free energy of different docking poses.18</p><p>A recently developed method called MELD (Modeling Employing Limited Data) considerably accelerates physics based simulations (like MD or Monte Carlo) using generic or vaguely directive information to restrict the search space of the problem.19,20 We call MELD × MD the application of MELD to MD simulations. MELD × MD is useful in protein structure determination,19 computing the poses and affinities of binding a peptide to a protein,21,22 and the folding of small proteins in CASP, the blind native prediction event.20,23–25 The acceleration in MELD × MD comes from a Bayesian integration of external information that can be probabilistic or combinatorial and not specific.</p><p>Here, we first use the ClusPro (CP) rigid-body docking server to estimate sets of protein–protein contacts from its 15 best docked poses.11,26–29 Those contacts are input to MELD × MD, which then explores the internal and relative degrees of freedom of the two proteins with replica-exchange sampling in a physical potential function. We find that MELD × MD adds value in predicting dimer structures by identifying the highest computed conformational population (lowest free energy). Our study here is limited to situations in which both proteins are given to have roughly the correct backbone trace, so that we can learn from our modeling the value of including the flexibility that MD provides and of computing populations (free energies) to pick out native dimer structures among options. During our MELD × MD simulation the conformations of the side chains and the backbone are allowed to fluctuate. Future work will investigate the use of this flexibility to refine the monomer backbone structure during the MELD × MD docking procedure, but this was not done here. Such refinement is particularly important when bound state monomer structures are not available.</p><!><p>Because a protein has so many degrees of freedom, it is computationally expensive to sample them and find its equilibrium states by molecular dynamics simulations in a physical force field. It is even more challenging to model two proteins binding to each other. MELD is a recently developed method that can accelerate MD simulations of proteins when there is some form of target knowledge about the relevant states. Unlike traditional constraining methods, MELD does not require that the constraints be precise, accurate, or deterministic. MELD can speed up the finding and sampling of important states using external information that is even loose, partly wrong, corrupted, probabilistic, or incomplete. In the present situation, we give MELD the directive that "most of the contacts predicted from the ClusPro rigid-body server, in one of top 15 poses, are likely to be essentially right" and that "the monomers should not unfold".</p><p>Here's how MELD uses smart springs to direct MD simulations to focus around promising regions of conformational space. The springs are smart in two ways. (i) Not all springs are active at the same time. MELD uses Bayesian inference to parse springs into different subsets that are active at different times. The active ones are always the ones with the lowest spring energy (i.e., the ones that are least violated). Switching between different sets of springs happens without violating detailed balance. (ii) The springs energy functions are flat-bottomed. This ensures that, when they are satisfied, no energy is added to the force field's Hamiltonian of the system. Relative populations of the different regions of the conformational space that satisfy some of the external information are consistent with the populations of unrestrained simulations and can therefore be used to compute proper free energy differences between them, based only on the force field.</p><p>The smart-springs approach focuses the computational effort around regions of the conformational space that are consistent with the data.19,20 For good sampling, MELD × MD uses the Hamiltonian and temperature replica exchange (HTRE) protocol.30,31 Higher temperatures correspond to moving up the replica ladder, which leads to weakening the springs.</p><!><p>In this paper, we study protein dimers. Our conformational space is the relatively limited number of docking poses of the two monomeric proteins, in addition to some flexibility in the internal structure of each protein. We enforce these two conditions separately. (i) We limit the docking space using the intermolecular contacts proposed by CP. In the lowest energy replicas, our system is forced to explore the conformational space that is compatible with a fraction of the contacts predicted by CP. We do this by computing the intermonomer contacts of the first 15 CP poses. To avoid kinetic biases of our replica exchange protocol that would make simulation convergence harder to reach, we randomly remove contacts from poses until every pose has the same number of contacts as the pose with the least contacts. The contacts of each one of the poses are collected in different groups and our Bayesian approach allows us to have only a fraction (70%) of the springs of one of these groups active at the bottom of the replica ladder at any moment. This allows some flexibility in the final docking contacts of our poses compared to the CP contacts that we use to create the smart springs in the MELD simulations. (ii) Monomers should not unfold during the simulation. We force the monomers to maintain the native secondary structure: we weakly restrain the position of the alpha-carbons of the first monomer around their native monomer positions, and we enforce 80% of the native monomer intramolecular contacts of the second monomer to be satisfied. Our approach allows for more flexibility in the two monomers than is currently used by other docking approaches. In particular, the side chains are totally free to reorient, and the structures of the main chains are free to fluctuate around the restraints. In this work we use the same approach for all 20 pairs of dimers, but in principle it is possible to tune the flexibility of different parts of the systems according to specific needs.</p><p>The restraints that drive the docking (the one at point i) are slowly turned on while descending the replica ladder (between α = 0.83 and α = 0.33); other restraints are on at all time. The monomers are forced at all times to sit within a sphere of radius RgR + 2.5RgL + 0.5 nm, where RgR and RgL are the radius of gyration of the receptor and of the ligand, respectively. The temperature is geometrically increased climbing the replica ladder, from 300 at α = 0 to 550 at α = 0.6, and after the temperature is constant. Due to the way we enforce the MELD restraints at the highest replica, the monomers are not docked and the ligand is free to rotate. This means that, in all our simulations, every system copy experiences several undocking and docking events in different poses while climbing and descending the replica ladder.</p><p>HTRE simulations are run using 30 replicas. Replicas start equally spaced in replica index space, but their positions are adapted during the simulation to optimize exchanges. Exchange between different replicas are attempted every 50, and the simulations are 2.5 long for each replica. The initial pair of structures is the same for all replicas, and they are structures in which both monomers have their native bound structure, but the two monomers are not in contact.</p><!><p>All our MD simulations are carried out using ff14SB-side force field (FF), which uses ff99SB32 parameters for the backbone and the recent ff14SB33 parameters for the side chains. We use the GB-neck2 implicit-solvent model.34 Simulations are run with steps of 4.5 using a Langevin integrator with a friction coefficient of 1.0 ps−1. Hydrogen masses are repartitioned to allow for the long time-steps. A cutoff of 1.8 is used for all interactions. MELD is a plugin to the OpenMM35 simulation package. Non-MELD simulations used to test the dimer stability are run using AMBER.36</p><!><p>Our aim is to show how a physically accurate (but also computationally expensive) approach can accurately identify the native state of dimers when more simple docking protocols fail to do so. Of the 49 small dimers we identified from the PDB37 that are stable in our simulation condition (see SI 1 for an in depth description of the procedure used to screen PDB dimers), vanilla CP (i.e., the CP run using the standard setting of the web server) is able to correctly dock 34 of them (i.e., the CP TOP1 structure is native-like). Our test set was comprised of the 15 dimers that CP failed to correctly predict, and 5 out of the 34 it correctly predicts. For 19 of these 20 systems CP did provide at least a good structure in the first 15 CP guesses, for the twentieth dimer we provided MELD also native contacts. In this way we show how our free-energy-based prediction can add value when it is needed.</p><p>Including protein flexibilities is critical for good estimates of the free energy. We enforce weak restraints on the receptor α carbon positions, and we enforce a fraction of the native contacts within the ligand. Backbone α carbon RMSDs of the receptor and the ligand go easily up to 2.5 Å (see SI Figure 3), and the side chains are totally free to reorient. While it is true that the starting conformation of the monomers are theirs bound structures, this structure evolves quickly to a few angstroms away in the MD protocol. CP predictions are the ones that are mostly influenced for using the bound structure of the monomers, but we use bound-complex docking as a reliable source of intermolecular contacts. In principle, we can use any source of information (experimental or computational) as long as we can define the confidence associated with such information. Also we note that we use a vanilla approach to CP. A more thoughtful use of CP might yield to better results, but this is outside the scope of this paper. Future work will explore how to use less reliable sources of information compared to bound docking (like template based docking or experimental data), and the simultaneous docking and refinement of monomer structures (e.g., starting from protein structure prediction or apo structures instead of holo structures). In this initial work, we are interested in showing how well we can predict native-like dimer structures using a free energy based method. The flexibility of the proteins in the simulations presented here goes exclusively toward the evaluation of accurate free energies of the docking poses.</p><!><p>The quality of our predictions is discussed in the context of the metrics used by the CAPRI community; the reported values in this paper are calculated using the DockQ software.38 For heterodimers, the ligand RMSD (LRMSD) is the ligand (i.e., the shorter chain) backbone RMSD computed after super-imposing the receptor (i.e., the longer chain) backbone. For most of our structures, which are homodimers, we consider the receptor to be the first chain and the ligand to be the second one in the native dimer PDB file. For the computation of the interface RMSD (iRMSD), an interface residue is defined as any residue that has a heavy atom less than 10 Å away from an heavy atom of a residue that belongs to the other monomer chain in the native structure. We also report the percentage of true positive (fTP) and false positive (fFP) contacts, where contacts are defined with a more restrictive 5 Å heavy-atom distance cut off. We report also the CAPRI score and the DockQ score, which are a combination of the RMSD and contact distance parameters.</p><!><p>Our predictions are obtained by clustering the trajectories of the five lowest replicas in the HTRE simulation. The first 250 ns of each trajectory (i.e., the first 5000 exchange attempts) are considered to be equilibration period for the systems and therefore are discarded from this analysis. We cluster the frames using the DBSCAN algorithm39 implemented in sklearn.40 As the distance metric between structures, we use the LRMSD computed on the backbone alpha carbons. For a structure to be included in a cluster, it should have a maximum distance of 3 Å from core structures in the same cluster. A structure is considered a core structure of a cluster when it has at least 200 neighbors. In SI 4 we show the small effect different clustering parameters and methods have on our predictions. The centroid of each cluster is then selected as representative structure, and our predictions are ranked according to cluster population, i.e. MELD × MD TOP1 (abbreviated as MELD TOP1) structure is the centroid of the most populated cluster. No information about the native structures or the CP poses enters the clustering protocol, making this effectively a prediction.</p><p>We assess the quality of our predictions by comparing the structure of the centroid of our most populated cluster (i.e., MELD TOP1 prediction) with the crystal structure of the native protein. We use this seemingly restrictive definition of success (only TOP1) because a correct ranking of the poses is a necessary feature of free energy based methods. In section 3.1.3 we will investigate the origin of the four failed predictions and provide suggestions to fix these failures.</p><p>Figure 1 shows our results. Each box shows the native structure of the dimer in white, with the receptor depicted as surface and the ligand as a semitransparent cartoon. For each system the ligand position and structure of the MELD TOP1 prediction is represented by a blue cartoon. Boxes are colored in green if the CAPRI score is at least acceptable. The five systems on the left of the figure are the ones that CP successfully ranked (i.e., the CP TOP1 is the best pose predicted by CP according to LRMSD metric). As it is possible to notice when CP is successful in predicting a close to native pose we are also successful in predicting the docking pose. When CP fails, we are often able (11/15 cases) to correctly predict a structure close to native, making our method complementary to what is in the field. In Table 1 we report LRMSD, iRMSD, fTP, fTN, CAPRI, and DockQscores for the MELD TOP1 predictions of the 20 dimers. Out of the 16 successes according to the CAPRI score: 2 are acceptable, 13 are medium, and 1 is high quality. This suggests that, for these systems, a physically accurate model that allows for an accurate representation of the entropy and energy fluctuations of the different poses is key to obtain an accurate ranking of the structures. The relatively high LRMSD of 1Z09 is due to a rearrangement of a ligand α helix. This rearrangement happens away from the interface so it does not dramatically affect the quality of the docking prediction.</p><p>Two factors contribute to the quality of our MELD TOP1 prediction. (i) We use the size of structure population as a proxy for a free energy evaluation. Thanks to the quality of the force field and solvent model, we are able to use/identify the most native set of CP intermolecular contacts, and to successfully dock 16/20 dimers. (ii) MELD ensures that our system sits in the vicinity of one of the first 15 CP poses (in terms of intermolecular contacts). The FF is then allowed to explore the conformational space surrounding the sets of contacts. This allows physics to identify better poses than the (vanilla) rigid docking protocol we used to predict contacts in 13/16 cases. Identifying a good set of contacts and letting physics explore the allowed conformational space to provide a good structure is a two step process in all our predictions, but we will look at them separately in sections 3.1.1 and 3.1.2. In section 3.1.3, we look at why we fail to predict the docked conformations of four systems correctly.</p><!><p>MD is grounded in physics. MELD focuses MD in order to sample only certain regions of the conformational space without biasing the relative probability of the allowed conformations. It is therefore possible to use the relative population of different structural clusters in our trajectory to identify the lowest free energy cluster (i.e., the most populated). Since we allow our simulation to explore only dimer conformations that share contacts with one of the first 15 CP poses, we could say that we are able to rerank CP poses. However, in this work we allow the system some freedom to move around CP predicted contacts. So it is not straightforward to use our simulations as a pure CP reranking tool. In principle, by limiting the search space around CP docking poses it is possible to make MELD × MD an almost "pure CP poses reranking method".</p><p>Figure 2 shows an example. For 2A2Y, the TOP1 CP structure (central column) has the dimer on the other side of the protein from the native crystal structure(left column). In the simulation, MELD × MD attempts to dock the two monomers according to the contacts from the 15 different CP poses. The fact that at the end of the simulation the centroid of the most populated cluster is extremely similar to the native structure (right column) shows how MELD was able to identify the best set of contacts to drive the monomer docking. Figure 3 shows how we consistently pick out the best pose between the 15 ones suggested by CP. It is possible to track which pose in the first 15 CP poses is favored by the FF by counting how often a pose is active in the lowest 5 replicas during our simulations. The histogram bars in Figure 3 show the relative populations of the first 15 CP poses in our MELD × MD simulations; the red lines identify the CP TOP15 (by LRMSD) for each system. It is possible to see how MELD × MD is able to consistently identify the best pose. 2LYJ, 1NS1, 4QR0, 2B87, and 1CFP have more than one good pose (LRMSD < 6.5 Å) within the first 15 predicted structures, making the populations split between those. Due to the relative relaxed way we enforce the restraint, there is some uncertainty in selecting which pose is active. First we compute the intermolecular contacts from the CP predictions. Then we force all poses to have the same number of contact as the pose with the fewest contacts by randomly removing contacts from poses that have more contacts. This is necessary to avoid kinetic bias in selecting poses along the RE ladder that could make convergence harder. Finally we enforce 70% of the contacts as active. This allows quite some freedom to explore around the original CP poses, making it relatively easy for similar poses to be selected because the FF can explore the same conformational space starting from different poses. This is not a problem for the clustering protocol we use to select the MELD TOP1 pose because it does not use any CP information, but it makes difficult to use this MELD × MD protocol as a purely CP reranking method.</p><!><p>As noted above, we enforce only 70% of the intermolecular contacts of any given poses, and we do not add a penalty for any additional intermolecular contact formed. This means that we allow some freedom for the monomers to explore the space around the contacts of the CP poses. An interesting example of this is shown in Figure 4. In this case, the CP TOP1 prediction is acceptable according to CAPRI metrics. MELD × MD correctly identifies that the set of contacts belonging to the CP TOP1 pose are closest to native and provides a structure that satisfies 76% of such contacts. The force field is then able to pick up and go the last mile in docking the dimer into a better structure, one that has a medium CAPRI score, a better rank than the CP pose we use to provide this set of contacts. A similar process happens in all the predictions, and only in 3 of the 16 successful predictions MELD × MD predicts slightly worse docking poses than vanilla CP does (see Table 1). This highlights also the quality of the FF for studying dimers.</p><!><p>The four failures to predict correct dimer conformations are linked with failures of the protocol we use, or with limitations of the physical model used to describe the system. The 3T0E prediction suffered a nonconverged RE simulation. Figure 5 show that the different simulations of the RE protocol did not explore a similar set of conformations. 1R7H is also partially affected by a not totally converged RE simulation. 1R7H and 3T0E are the only two systems for which the population of the MELD TOP1 prediction is below 60% of the overall number of clustered frames. As in protein folding, low MELD cluster populations seem to flag systems for which the answer is uncertain. In 2EZX and 1R7H the ligands change conformations, making docking in our simulation time scale impossible to converge. For example, let us say that the monomers look like two hands, and the dimer native structure looks like a handshake. To achieve a correct docking the two hands must be open when the docking is attempted. It is impossible to dock correctly the two monomers if when undocked the ligand closes like a fist. This is what happens to the 2EZX and 1R7H systems. The information that is in place fails to keep the hand open, since the contacts of an open hand are a subset of the contacts of a closed one and the FF/solvent model favors a closed conformation of the ligand. A possible way to recover this is to penalize the formation of too many new intramolecular contacts, or to use a full intramolecular distances matrix instead of only a contact matrix to limit ligand conformation explorations. Interestingly, Figure 2 shows how the most often active sets of CP information for these two systems are in both cases the ones of the CP TOP15 structures. It is reasonable to assume that limiting ligand conformations for these two systems will recover the correct predictions. We briefly note here that properly tuning the flexibility of receptors and ligands has an important effect on MELD × MD ability to successfully dock dimers. Flexibility increases the size of the search space to explore. Rigidity might hinder moves that can in some case be necessary to wrap monomers around each other. In this paper we present a uniform approach for all the 20 test cases that works for most of them, but some docking attempts would benefit from individualized approaches. Finally, 1H0X prediction fails because the FF favors a different docking pose. A small bump in the third panel of the LRMSD probability highlights that the systems visits a set of native-like conformations, but this cluster appears not be the one favored by the FF. In SI 3 we show how the 16 systems that MELD × MD correctly predicts are not affected by the limitation discussed here.</p><p>In Table 2 we look at the same quantities presented in Table 1, computed over the MELD TOP10 structure for the four systems where MELD × MD fails to predict native, and for 1Z09, for which a change in the structure of the ligand affects some of the score. All the systems affected by a change in the ligand structure (1Z09, 1R7H, and 2EZX) have small population clusters of medium quality that improve on the metrics of the TOP1 prediction. This means that when the system samples a correct ligand structure, MELD × MD is able to correctly dock the dimer. These correct ligand structures are simply sampled too rarely compared to incorrect ones, causing native docking poses not to be the most populated ones.</p><p>It is likely that limiting the flexibility of the ligand in our docking protocol would allow these three systems to be predicted correctly. It is also known that biologically relevant homodimers have a high degree of symmetry,42,43 and this implies that the structure of the two monomers needs to be similar. We can in principle use this knowledge to skim our predictions in order to eliminate homodimer predictions with significantly disparate monomer structures. We show in SI 6 that by doing so we can recover the correct predictions of these three dimers without affecting the prediction quality of the other successes.</p><p>1H0X and 3T0E are not dramatically affected by a structural rearrangement of the ligand, and they do not sample structure close to native. Failure in predicting the native state of these monomer is linked to issue of convergence of the REMD protocol and to limitations in the force field and solvent model.</p><!><p>In section 3.1 we showed the success of our protocol in predicting the structure of dimers when we provide good information. Good information comes from the CP contacts prediction based on bound-complex docking and on exploring monomer conformations close to their bound structures (albeit some ligands are able to evolve to structures up to 10 Å away during the simulation). This is consistent with the aim of this paper to use flexibility to evaluate correctly the relative stability of the different binding poses. In principle, the flexibility we introduce in the system can at the same time be used to refine the monomer structures (side chains and backbone). To test to which extent our protocol is able to handle monomer structural refinement, we use the CAPRI T120 target, from the joint CAPRI round 37 - CASP 12 competition.13,44 As starting structure for the monomers we use the model # 1 submission of the Zhang-Server group (TS479),45 based on I-TASSER46–48 for the two monomers (T0921 and T0922). The monomers have a RMSD to their native structure of 3.0 and 3.4 Å, respectively. We feed these to vanilla CP and we then use the predicted contacts to limit the search space of our simulation. As the initial structure of our simulation we use the two predicted monomer structures not in contact. We run our MELD × MD simulation using the same protocol as described above. Table 3 shows a summary of the TOP1, TOP3, and TOP5MELD predictions and the TOP15 CP prediction. Our protocol is not optimized to handle the uncertainty of the monomer structures or in the intermonomer contacts. It is therefore not surprising that we fail to predict a MELD TOP1 structure close to native. However, we did identify some acceptable structures (2 and 3), and the first acceptable structure is ranked higher than the best CP structure (3). This small improvement on the CP prediction is promising. We recognize that, while this result is promising, future work (well outside the scope of this paper) is needed to fully assess the value that MELD × MD provides in more challenging conditions where monomers need to change their shapes to successfully dock. Different classes of problems fit this description: in template based docking monomers need refinement, in docking from monomer apo structure some structural rearrangement is in order, and in unstructured protein docking a folding-like process is necessary upon docking.</p><!><p>We describe a method that gives improved protein–protein structure predictions. It starts with the 15 best structures produced by the ClusPro rigid-docking server. It then uses the interprotein contacts predicted by CP for each structure as restraints that are imposed probabalistically within the MELD Bayesian method. MELD-accelerated MD then gives improved dimer structures and produces populations (free energies) that are predictive of which dimer conformation is the best of the lot. The protocol we presents performs well if we start with the monomers having the bound backbone conformation but also demonstrates some improvement when the backbone conformation come from predictions as demonstrated in CASP/CAPRI example. In this last case, more aggressive sampling of the monomer structure would benefit the prediction. While the MELD × MD component is computationally expensive, nevertheless this approach is a promising way to compute good protein–protein binding structures.</p>