Protein–Protein Docking with Large-Scale Backbone Flexibility Using Coarse-Grained Monte-Carlo Simulations

Most of the protein–protein docking methods treat proteins as almost rigid objects. Only the side-chains flexibility is usually taken into account. The few approaches enabling docking with a flexible backbone typically work in two steps, in which the search for protein–protein orientations and structure flexibility are simulated separately. In this work, we propose a new straightforward approach for docking sampling. It consists of a single simulation step during which a protein undergoes large-scale backbone rearrangements, rotations, and translations. Simultaneously, the other protein exhibits small backbone fluctuations. Such extensive sampling was possible using the CABS coarse-grained protein model and Replica Exchange Monte Carlo dynamics at a reasonable computational cost. In our proof-of-concept simulations of 62 protein–protein complexes, we obtained acceptable quality models for a significant number of cases.


Introduction
Protein-protein interactions are fundamental in many biological processes. Their structural characterization is one of the biggest challenges of computational biology. A variety of docking methods are currently available for structure prediction of proteinprotein complexes [1,2]. They can be divided into free (global) and template-based docking. Free (global) docking methods are designed to generate many distinct binding configurations. Template-based methods restrict docking to a binding mode found in a structural template. As demonstrated in the blind docking challenge, Critical Assessment of Prediction of Interactions (CAPRI), template-based methods generate more accurate results but only if a good quality template exists [1][2][3][4][5]. In some cases lacking useful templates, free global docking can yield acceptable results. According to recent estimates, the best free docking methods find adequate models among the top 10 predictions for around 40% of the targets [1]. The CAPRI analysis also indicates that protein backbone flexibility is a big challenge; protein complexes that undergo substantial conformational changes upon docking get no successful predictions from any method [3][4][5].
Presently, most of the free docking methods treat the backbone of input protein structures as rigid. This approximation reduces the protein-protein docking problem to a 6D (three rotational and three translational degrees of freedom) search space. Rigid-body search for the binding site most often rely on the Fast Fourier Transform [6][7][8]. Other successful approaches include 3D Zernike descriptor-based docking [9,10] or geometric hashing [11]. These rigid-body methods are often used as a first docking step, followed by scoring [12][13][14][15][16], using experimental data [17] and/or structural refinement to capture backbone flexibility [5,18]. Molecular Dynamics is perhaps the most common refinement strategy, either in classic or enhanced sampling versions [17,[19][20][21][22]. Other tools use rotamer libraries to address side-chain flexibility [23] and Elastic Network Models (ENMs) for modeling backbone rearrangements [24][25][26][27][28]. Accounting for backbone flexibility in

Results
The most accurate models (out of the sets of 10,000 generated models and 10 topscored) are characterized in Table 1. The table presents different metrics of similarity to the experimental structures for the set of 62 protein-peptide complexes (divided into three categories: low, medium and high flexibility cases). To assess the sampling performance, below we will use the iRMSD values for the best models out of all models. According to the iRMSD values the CABS-based docking algorithm produced a significant number of near-native protein-protein arrangements of acceptable quality (iRMSD < 4 Å, according to CAPRI criteria) for most protein-protein cases in the categories of low and medium flexibility cases. However, in the high-flexibility category, the best iRMSD values were noticeably higher (in the range of 4-12 Angstroms). This resulted from the adopted distance restraint scheme (see the Methods section), which was uniform for whole proteins and introduced a penalty for deviations of more than 1 Å from the input structures (unbound experimental structures). This penalty was very small for the protein ligands. Thus, the distance-restraints scheme allowed for the large-scale conformational changes, however, they might have prevented binding-induced conformational changes in the high-flexibility category. Therefore, there is the need to modify such a scheme for the most challenging targets. The results analysis below focuses on the sampling performance for the selected low-flexibility barnase/barstar case. Figure 1 characterizes iRMSD versus CABS model energy values for the barnase/barstar (1BRS) and another low-flexibility case with clearly the lowest iRMSD value 1.09 Angstroms (2SNI). The results analysis below focuses on the sampling performance for the selected low-flexibility barnase/barstar case. Figure 1 characterizes iRMSD versus CABS model energy values for the barnase/barstar (1BRS) and another low-flexibility case with clearly the lowest iRMSD value 1.09 Angstroms (2SNI). Characterization of docking results using RMSD to the X-ray structure and system energy. The left panels show the interface-RMSD versus CABS energy values. Point color represents the temperature-from yellow (high) to pink (low). The molecular visualizations show X-ray structures and ensembles of predicted models corresponding to selected energy minima (numbered in the picture from 1 to 3). As presented in the picture, the minima numbered as 1st corresponds to near-native protein-protein arrangements, others to non-native ensembles, as presented in the picture. Characterization of docking results using RMSD to the X-ray structure and system energy. The left panels show the interface-RMSD versus CABS energy values. Point color represents the temperature-from yellow (high) to pink (low). The molecular visualizations show X-ray structures and ensembles of predicted models corresponding to selected energy minima (numbered in the picture from 1 to 3). As presented in the picture, the minima numbered as 1st corresponds to near-native protein-protein arrangements, others to non-native ensembles, as presented in the picture. The presented ensembles are the sets of similar models found in the structural clustering of contact maps (see Methods). The figure shows two modeling cases: 1BRS and 2SNI. Figure 2 shows example ensembles of barnase/barstar models and the most accurate model (iRMSD 1.9 Å). A single system replica could explore an ample conformational space that involved significantly different binding configurations and protein-ligand conformations, as demonstrated in Figure 2c and Movie S1. Figure 3a further characterizes this single replica's using iRMSD and LoRMSD (RMSD for ligand only) values. As presented in the figure, the ligand structure fluctuated around 5 Å (the same fluctuations in the context of all replicas are shown in Figure 3b). The ligand became significantly more closer to the X-ray structure after binding to the native binding site as reflected by iRMSD values. Namely, after correct binding, LoRMSD values got noticeably lower to around 2 Å (see Figure 3a). In the following sections, we do not discuss this aspect of our method; however, it is worth mentioning that the proposed method enabled a detailed analysis of plausible docking trajectories. The described docking procedure uses REMC protocol enhanced by simulated annealing of all 20 replicas. Figure 3c shows their evolution through different temperatures. Figure 4 provide more detailed pictures of structural flexibility for protein "receptor" and "ligand". Protein-protein contacts defining the complex assembly are characterized in Figure 5. In the presented example, the most persistent protein-protein contacts occurred in about 15% of snapshots. Therefore, they were significantly less stable compared to intramolecular protein contacts ( Figure 4). space that involved significantly different binding configurations and protein-ligand conformations, as demonstrated in Figure 2c and Movie S1. Figure 3a further characterizes this single replica's using iRMSD and LoRMSD (RMSD for ligand only) values. As presented in the figure, the ligand structure fluctuated around 5 Å (the same fluctuations in the context of all replicas are shown in Figure 3b). The ligand became significantly more closer to the X-ray structure after binding to the native binding site as reflected by iRMSD values. Namely, after correct binding, LoRMSD values got noticeably lower to around 2 Å (see Figure 3a). In the following sections, we do not discuss this aspect of our method; however, it is worth mentioning that the proposed method enabled a detailed analysis of plausible docking trajectories. The described docking procedure uses REMC protocol enhanced by simulated annealing of all 20 replicas. Figure 3c shows their evolution through different temperatures. Figure 4 provide more detailed pictures of structural flexibility for protein "receptor" and "ligand". Protein-protein contacts defining the complex assembly are characterized in Figure 5. In the presented example, the most persistent protein-protein contacts occurred in about 15% of snapshots. Therefore, they were significantly less stable compared to intramolecular protein contacts ( Figure 4).
An essential and unique feature of the presented docking simulations is the level of backbone flexibility during docking. In the example above, the ligand backbone fluctuations (LoRMSD) were in the range of 2-7 Å (Figure 3b), with the average LoRMSD value of 3.3 Å from the entire docking simulations. In other cases, the ligand fluctuations were at a similar level or higher (see LoRMSD values in Table 1).  10,000 models combined from 20 replicas (500 models per replica) in which the highly flexible ligand is covering the entire surface of the flexible receptor; (c) 500 models from one replica only, (d) the best model obtained for barnase/barstar system (the X-ray structure of the ligand is shown in thick ribbon, the modeled in thin ribbon).  10,000 models combined from 20 replicas (500 models per replica) in which the highly flexible ligand is covering the entire surface of the flexible receptor; (c) 500 models from one replica only, (d) the best model obtained for barnase/barstar system (the X-ray structure of the ligand is shown in thick ribbon, the modeled in thin ribbon).
An essential and unique feature of the presented docking simulations is the level of backbone flexibility during docking. In the example above, the ligand backbone fluctuations (LoRMSD) were in the range of 2-7 Å (Figure 3b), with the average LoRMSD value of 3.3 Å from the entire docking simulations. In other cases, the ligand fluctuations were at a similar level or higher (see LoRMSD values in Table 1).
Finally, using structural clustering of contact maps (see Methods), we attempted to select the set of 10 top-scored models for each case.  Example simulation snapshots illustrate th ligand is presented in rainbow colors, the receptor in magenta. The lowest iRMSD mo from X-ray structure) is presented on the right lower corner superimposed on the X-ra (the X-ray structure is shown in thick lines, the predicted model in thin lines). (b) Li RMSD (LoRMSD) values for all replicas. The thick red line presents selected replica. (c) E system replicas between different temperatures driven by Replica Exchange Monte Car system. The thick red line presents selected replica. The replica trajectory is also presen Video S1.   Finally, using structural clustering of contact maps (see Methods), we attempted to select the set of 10 top-scored models for each case. Table 1 reports the most accurate models out of the 10 top-scored.

Discussion
This work demonstrates a significant improvement in the sampling of large-scale conformational transitions during global protein-protein docking compared to other Figure 5. Characterization of barnase/barstar contacts. Panels show barnase/barstar models and contact maps for entire simulation (all models, 10,000 models) and single selected replica (replica 6, 500 models) that reached a near-native arrangement. In the maps, green circles mark the native contacts.

Discussion
This work demonstrates a significant improvement in the sampling of large-scale conformational transitions during global protein-protein docking compared to other stateof-the-art approaches. We show that modeling the large conformational changes is possible at a relatively low computational cost. The presented simulations took between 10 and 80 h (depending on the system size) using a single standard CPU. The proposed modeling protocol can be used as the docking engine in template-based and integrative docking protocols using experimental structural data and additional information from various sources [2,40]. We focused on the free docking of protein ligands with a highly flexible backbone in the present test simulations. Using unbound structures as the input, we produced acceptable accuracy models (iRMSD around 4 Å or lower) in low-flexibility and medium-flexibility cases. However, the selection procedure of the most accurate models needs further improvement. Namely, selecting the best-ranked models led to acceptable models in about half of the tested cases.
Presently, the most common approach to account for conformational changes in protein docking is using ENM [24][25][26][27][28][36][37][38]. The applicability of ENM to modeling protein flexibility is limited to specific systems and depends on how collective the protein motions are. Our method presents a conceptually different approach that seems to be more realistic (see review discussing coarse-grained CABS dynamics in the context of ENM approaches [24]). We demonstrated that it is possible to simulate effectively free docking of highly flexible protein ligands to quite elastic protein receptor structures. Such a significant degree of flexibility was achieved using a highly efficient simulation engine based on the coarse-grained representation of protein structures, Monte Carlo dynamics, and knowledge-based force field. CABS coarse-graining, enhanced by the discretized protein model and interaction patterns, significantly reduces the search space. Monte Carlo dynamics, enhanced by Replica Exchange annealing, leads to huge speed-up of the search procedures. Additionally, a significant (although acceptable for many problems) flattening of energy surfaces by statistical potentials of CABS model simplifies simulations. As a result the flexible docking using CABS-dock is orders of magnitude faster than equivalent simulations based on classical modeling methods. Obviously, the new method also has several limitations that have to be considered when designing new computational experiments. First, since the "ligand" protein is treated as a very elastic object (what is necessary to guarantee efficient search of the binding sites and poses) the cost of computations rapidly grows with the protein size. Thus, completely free global docking of protein ligands larger than 150 residues (see Table 1) may be impractical. Second, the coarse-graining of the sampling space and simplifying interaction patterns (so important for the huge acceleration of the simulations) makes the docking energetics less sensitive. For these reasons, the clustering procedures, refinement of the resulting structures, and final model selection become challenging and need further development. Additionally, speeding-up the entire protocol can be useful. We estimate that the simulations could be easily speeded-up at least 10 times or more through algorithm parallelization. The speed-up would enable making the protocol available as the publicly accessible and automated web service.

Docking Simulation Protocol
In this work, we present the protein-protein docking simulation protocol that relies on the CABS coarse-grained model. The CABS design and applications have been recently described in the reviews on protein coarse-grained [29] and protein flexibility [24,41] modeling. Here we outline only its main features. The CABS model uses a coarse-grained representation of protein chains (see Figure 6), Replica Exchange Monte Carlo (REMC) dynamics, and knowledge-based statistical potentials. Representation of protein chains is based on C-alpha traces, restricted to an underlying high-resolution lattice. The lattice spacing allows slight fluctuations of the C-alpha-C-alpha distances and many pseudobonds orientations. Virtual pseudo-atoms are placed in the centers of these C-alpha-C-alpha bonds and are used to locate the main-chain hydrogen bonds. Additionally, the positions of the two pseudo-atoms representing side chains are defined by the geometry of C-alpha traces and amino-acid identities. Such fixed positions of side chains (taken from the statistics of protein databases) reduce the model's resolution. However, this limitation is less serious than it may appear since even small movements of the main chain (allowed due to the soft nature of the assumed geometrical restrictions) leads to large moves of the side chains. This way, the packing of side chains can be quite accurate. The interaction scheme of CABS consists of statistical potentials mimicking effects of main chain rotational preferences, main-chain hydrogen bonds, and side-chain contacts. All statistical potentials, derived from structural regularities observed in PDB structures, have relatively broad minima compensating the low-resolution effects and allowing a fast search for global energy minima. The solvent is treated implicitly, and its averaged effects are encoded within the above-mentioned contact potentials. Energy computation for protein chain models is very fast since many interactions could be pre-computed (and coded in large tables) due to the discretized patterns of main chains geometry. The Monte Carlo sampling of CABS uses a set of local movers. The resulting model dynamics is quite realistic for large-scale distances, allowing coarse-grained modeling of protein structures, dynamics, and protein-protein interactions. contacts. All statistical potentials, derived from structural regularities observed in PDB structures, have relatively broad minima compensating the low-resolution effects and allowing a fast search for global energy minima. The solvent is treated implicitly, and its averaged effects are encoded within the above-mentioned contact potentials. Energy computation for protein chain models is very fast since many interactions could be pre-computed (and coded in large tables) due to the discretized patterns of main chains geometry. The Monte Carlo sampling of CABS uses a set of local movers. The resulting model dynamics is quite realistic for large-scale distances, allowing coarse-grained modeling of protein structures, dynamics, and protein-protein interactions. Figure 6. Comparison of the all-atom (left) and the CABS coarse-grained model representation (right) for an example tripeptide. In the CABS model, protein residues are represented using C-alpha, C-beta, united side-chain atom, and the peptide bond center [29].
The modeling protocol consists of the following steps: 1. Preparing input structures of a protein-ligand and a protein-receptor. The protocol requires the input of two protein structures (single-or multi-chain) in the PDB format. One of them has to be indicated as a ligand and the second as a receptor. The ligand undergoes large conformational fluctuations, translations, and rotations around the receptor within the proposed protocol. The "ligand" should be a smaller protein because the computational cost of searching its conformational space rapidly grows with the chain length. That is because the motion of the entire structure (including fold relaxation, rotation, and translation of the entire molecule) is simulated by a random sequence of local moves. The accuracy of such sampling is acceptable for not too-large proteins. On the other hand, treating the "ligand" as a fully flexible object allows approximate studies of entire docking trajectories. In some cases, it would be perhaps worth treating a larger protein (but not too large) as a flexible "ligand", although this was out of range of the present studies.
2. Generating starting structures. Starting conformations are built using C-alpha coordinates only (in the CABS model C-alpha traces define the position of other united pseudo-atoms, see details [29]). The algorithm places the protein-ligand center at 20 random positions around the protein receptor at the approximate distance of 20 Å from the protein receptor's surface. Next, these protein-ligand Figure 6. Comparison of the all-atom (left) and the CABS coarse-grained model representation (right) for an example tripeptide. In the CABS model, protein residues are represented using C-alpha, C-beta, united side-chain atom, and the peptide bond center [29].
The modeling protocol consists of the following steps:

1.
Preparing input structures of a protein-ligand and a protein-receptor. The protocol requires the input of two protein structures (single-or multi-chain) in the PDB format.
One of them has to be indicated as a ligand and the second as a receptor. The ligand undergoes large conformational fluctuations, translations, and rotations around the receptor within the proposed protocol. The "ligand" should be a smaller protein because the computational cost of searching its conformational space rapidly grows with the chain length. That is because the motion of the entire structure (including fold relaxation, rotation, and translation of the entire molecule) is simulated by a random sequence of local moves. The accuracy of such sampling is acceptable for not too-large proteins. On the other hand, treating the "ligand" as a fully flexible object allows approximate studies of entire docking trajectories. In some cases, it would be perhaps worth treating a larger protein (but not too large) as a flexible "ligand", although this was out of range of the present studies.

2.
Generating starting structures. Starting conformations are built using C-alpha coordinates only (in the CABS model C-alpha traces define the position of other united pseudo-atoms, see details [29]). The algorithm places the protein-ligand center at 20 random positions around the protein receptor at the approximate distance of 20 Å from the protein receptor's surface. Next, these protein-ligand systems are used as starting conformations for the 20 replicas in the REMC CABS sampling scheme (each replica starts from a different ligand-receptor arrangement).

3.
Docking simulations using CABS coarse-grained model and REMC dynamics. During simulations, the protein receptor structure is kept close to the starting structure using distance restraints. Distance restraints are generated using the input coordinates of the C-alpha atoms. Two residues are automatically restrained if two conditions are met. First, their separation along the sequence has to be at least five residues. Second, the distance between their C-alpha atoms must be within the range of 5-15 Å. During simulations, the receptor restraints imply small-scale fluctuations of the protein receptor backbone in the range of 1 Å and, accordingly, more significant fluctuations of the side-chain atoms. A similar restraints scheme is applied to the protein-ligand but with tenfold weaker weights. During simulations, the ligand moves freely within the vicinity of the receptor and internal restraint allows for large-scale fluctuations of its structure. Usually, the ligand fluctuations are within the range of 2 and 12 Å to the input structure although folding-unfolding events are possible at highest temperatures. The docking simulation is conducted using CABS REMC pseudo-dynamics with simulated annealing. In this work, 20 replicas and 20 annealing steps have been used. All the REMC scheme parameters have been adjusted to allow for largescale conformational transitions, rotations, and translations of the protein-ligand in a reasonable computational time. The modeling protocol collects trajectories from all 20 replicas. The protocol saves only a small fraction (2%) of the generated models for further analysis i.e., 500 models from each replica, thus 10,000 models in total.

4.
Reconstructing to CABS coarse-grained representation. The set of 10,000 models in C-alpha traces are reconstructed to complete CABS model representation using CABS algorithm [29]. In CABS, positions of C-beta and Side-Chain united atoms are defined by the positions of the three consecutive C-alpha atoms and the amino acid identity (the most probable positions from the PDB database are used).

5.
Clustering of contact maps. First, for all of the 10,000 models the contact maps between the receptor and the ligand proteins are calculated. Two residues are considered to form a contact if their Side Chain pseudoatoms are at most 6 Å apart (for Alanines the C-beta atoms are considered as the Side Chain; for Glycines-it's the C-alpha atoms). Next, the algorithm sorts the models according to the number of the receptor-ligand contacts, and the set of top 1000 is kept for further processing. This way the transient and weakly bound complexes are removed from the solutions pool.
In the next step, the 1000 contact maps are clustered together to identify the most frequently occurring contact patterns. The complete link hierarchical clustering was used with the Jaccard index as the distance metric between contact maps. Finally, the identified clusters are ranked according to their density, defined as the number of the cluster members divided by the average metric between them.

6.
Reconstructing to all-atom representation. Representative models from the ten most dense clusters are reconstructed to all-atom representation using Modeller-based rebuilding procedure [42] (or can be reconstructed using other rebuilding strategies, see review [43]).
In recent years, the CABS model has been used for modeling the flexibility of globular proteins [44][45][46][47] and various processes leading to large-scale conformational transitions. These included: ab initio simulations of protein folding mechanisms [48,49], folding and binding mechanisms [49,50], and free protein-peptide docking within the CABS-dock tool [51][52][53][54][55][56][57]. The CABS-dock is a well-established peptide docking tool that has been made available as a web server [51,52] and, most recently, as a standalone application [54]. Its distinctive feature among other tools is the possibility of fast simulation of the large backbone rearrangements of both peptide and protein receptors during binding (see the review on protein-peptide docking tools [58]). In addition, the CABS-dock has been used in multiple applications (recently reviewed [56]), including docking to receptors with disordered fragments [41,59], GPCRs [60], and modeling proteolysis mechanisms [61].
The presented protocol for protein-protein docking utilizes the CABS-dock standalone package [54] developed primarily for protein-peptide docking. In order to tackle the protein-protein docking problem, key changes have been made to the docking algorithm that aimed mainly at the improvement of the conformational sampling. First of all, the temperature distribution between replicas in the REMC scheme was adjusted. Instead of constant temperature increment between consecutive replicas, as in the original CABSdock, here we've implemented progressive geometric raise of the temperature increment. Furthermore, the number of simulation replicas was increased to twenty versus ten in the original CABS-dock. Besides the sampling improvement, a new clustering protocol was introduced. The original CABS-dock used RMSD-based clustering, which worked well for peptides. For the protein-protein complexes, however, purely geometrical similarity condition such as the RMSD is too severe. Namely, for two binding poses, where the mobile protein was docked in the exact same pocket but is slightly tilted in one of them, the RMSD difference would be considerable. Despite representing similar binding poses, the two structures would end up in different clusters. To overcome this, the current protocol uses clustering based on the similarity between receptor-ligand contact maps.

Results Analysis and Quality Metrics
The docking simulation analysis was performed using Python and NumPy (Python library). Structural differences between experimentally determined structures and generated models were evaluated using Root Mean Square Deviations (RMSDs). Interface RMSD (iRMSD) is an RMSD calculated for interface residues of the receptor and the ligand separated by no more than 6 Angstroms. Ligand RMSD (LRMSD) is an RMSD computed for the ligands after the superimposition of the receptors. Ligand only RMSD (LoRMSD) is an RMSD computed for the ligand structure only. Root Mean Square Fluctuation (RMSF) is a measure of the amino acid's flexibility. It is calculated for every residue as the square root of this residue's variance around the reference residue position. The fraction of native contacts (fNAT) was calculated as a number of experimental structure contacts found in the generated structure divided by the total number of contacts found in the experimental structure. Rather restrictive contact criterion, distance up to 6 Å between side-chain centers, was used. All figures presented in this work were generated using PyMOL, UCSF Chimera, and Matplotlib (Python library).

Dataset
In this docking study, we used protein-protein cases from the ZDOCK benchmark set [62] (cases in which a smaller size protein-a protein-ligand-contained more than one protein chain, or chain gaps, were discarded from our set). The set comprises the three flexibility-based subsets: low-flexible (almost rigid), medium-flexible, and highly flexible with available unbound X-ray structures of both the protein-receptor and the protein-ligand. The unbound structures were used as the docking input. As the reference for calculating various similarity measures, we used the X-ray structures of the protein-ligand complexes. Table 1 lists all the PDB IDs of X-ray structures used in the study.

Conclusions
In summary, the described docking procedure accounts for large-scale protein structure fluctuations during unrestrained protein-protein docking search for the binding site. The exploration of such vast conformational space has not been demonstrated before to the best of our knowledge. The approach shows unprecedented sampling possibilities; however, the accuracy of the obtained complexes is still lower than observed for state-ofthe-art docking tools. Definitely, the balancing of the structural restraints scheme needs further developments and tests. Therefore, this work is the first step towards a mature protein-protein docking tool. The next development steps would involve modifications of the distance restraints scheme, which allow for different degrees of flexibility for appropriate protein fragments (now the presented algorithm treats the entire protein-ligand as very flexible) and force-field improvements. The proposed approach is also very promising in the refinement applications when searching for the binding site is not needed, and only the protein-protein interface needs to be optimized.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/ijms22147341/s1, Video S1. The trajectory of a single replica from the protein-protein docking simulation of barnase/barstar system. The movie shows the barnase receptor in surface representation and the barstar ligand in ribbon. The presented replica reached the model with interface RMSD value 1.9 Angstrom from the complex X-ray structure, shown as transparent ribbon.