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

**Mateusz Kurcinski \*, Sebastian Kmiecik \* , Mateusz Zalewski and Andrzej Kolinski**

Biological and Chemical Research Centre, Faculty of Chemistry, University of Warsaw, 02-089 Warsaw, Poland; mateusz.zalewski.fuw@gmail.com (M.Z.); kolinski@chem.uw.edu.pl (A.K.)

**\*** Correspondence: mkurc@chem.uw.edu.pl (M.K.); sekmi@chem.uw.edu.pl (S.K.)

**Abstract:** 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.

**Keywords:** protein–protein interactions; protein–protein binding; protein–protein complex; coarsegrained modeling; multiscale modeling

#### **1. 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 protein– protein 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–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–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–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–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–22]. Other tools use rotamer libraries to address side-chain flexibility [23] and Elastic Network Models (ENMs) for modeling backbone rearrangements [24–28]. Accounting for backbone flexibility in

**Citation:** Kurcinski, M.; Kmiecik, S.; Zalewski, M.; Kolinski, A. Protein–Protein Docking with Large-Scale Backbone Flexibility Using Coarse-Grained Monte-Carlo Simulations. *Int. J. Mol. Sci.* **2021**, *22*, 7341. https://doi.org/10.3390/ ijms22147341

Academic Editors: Istvan Simon and Csaba Magyar

Received: 19 June 2021 Accepted: 4 July 2021 Published: 8 July 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the search for the binding site significantly increases the docking complexity and makes it practically intractable using conventional all-atom modeling approaches. This enormous computational complexity of flexible docking can be reduced using coarse-grained protein models [29–32]. The best-performing methods that can now include backbone flexibility during the docking calculations use coarse-grained models and/or ENM-driven simulations. These include RosettaDock combining coarse-grained generation of backbone ensembles and all-atom refinement [33–35]; ATTRACT combining coarse-grained docking with ENM and all-atom refinement [36,37]; and SwarmDock using all-atom ENM [25,38]. All these approaches show some advantages in modeling protein flexibility compared to rigid-body docking followed by structure refinements. However, effective modeling flexibility in protein–protein docking remains an unsolved problem, as demonstrated in the recent CAPRI round [25,35,37,39].

In this work, we use a well-established CABS coarse-grained protein model [29] for protein–protein docking. During the CABS docking simulation, one of the docking partners undergoes a long random process of rotations, translations, and extensive backbone conformational rearrangements that significantly modify its fold. Simultaneously, the backbone of the second protein undergoes small fluctuations.

#### **2. 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.

**Table 1.** Summary of the docking simulations. The table characterizes X-ray data used in the docking, average ligand flexibility, and docking results. The table reports the best accuracy models out from all (10,000) and 10 top-scored models. The metrics definitions are provided in the Methods section. The table divides the presented cases on the three categories: low-flexibility, medium-flexibility and highly flexible cases.



**Table 1.** *Cont.*


**Table 1.** *Cont.*


**Table 1.** *Cont.*

\* 4-letter PDB code for the crystal structures used in this study. \*\* The RMSD (in Å) of the interface Cα atoms for input receptor and ligand after superposition onto the co-crystallized complex system.

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

**Figure 1.** 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).

**Figure 2.** Protein–protein docking stages illustrated by barstar/barnase docking case. The figure shows the barnase receptor in magenta and the barstar ligand in rainbow colors. The respective panels show: (**a**) 20 starting structures for each replica of the system; (**b**) 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. Table 1 reports the most accurate models out of the 10 top-scored.

**Figure 3.** Docking trajectory for the selected replica of barnase/barstar system. The presented replica reached the most accurate barnase/barstar complex structure. (**a**) iRMSD (interface RMSD) and LoRMSD (ligand only RMSD) values. Example simulation snapshots illustrate the plot. The ligand is presented in rainbow colors, the receptor in magenta. The lowest iRMSD model (1.9 A from X-ray structure) is presented on the right lower corner superimposed on the X-ray structure (the X-ray structure is shown in thick lines, the predicted model in thin lines). (**b**) Ligand only RMSD (LoRMSD) values for all replicas. The thick red line presents selected replica. (**c**) Exchange of system replicas between different temperatures driven by Replica Exchange Monte Carlo (REMC) system. The thick red line presents selected replica. The replica trajectory is also presented in the Video S1.

**Figure 4.** Characterization of barnase/barstar flexibility in the docking simulation. The figure shows RMSF plots (upper panels) and contact maps (lower panels) for (**a**) the barnase receptor and (**b**) the barstar ligand. The RMSF profile (see Methods) and contact maps showing the frequency of contacts are derived from the entire simulation (derived from 10,000 models).

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

#### **3. 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–28,36–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.

#### **4. Methods**

#### *4.1. 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.

**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:


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


In recent years, the CABS model has been used for modeling the flexibility of globular proteins [44–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–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.

## *4.2. 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).

#### *4.3. 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.

#### **5. 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.

**Author Contributions:** Conceptualization, M.K., S.K. and A.K.; Methodology, M.K. and A.K.; Software, M.K. and M.Z.; Supervision, S.K.; Visualization, M.K. and S.K.; Writing—original draft, S.K.; Writing—review & editing, M.K., S.K., M.Z. and A.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** SK: MZ and AK acknowledge funding by the National Science Centre, Poland [MAE-STRO2014/14/A/ST6/00088]. MK acknowledges funding by the National Science Centre, Poland [501/D112/66 GR-6271]. SK also acknowledges funding in part by the National Science Centre, Poland [UMO-2020/39/B/NZ2/01301].

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

