Enhanced sampling of protein conformational states for dynamic cross‐docking within the protein‐protein docking server SwarmDock

Abstract The formation of specific protein‐protein interactions is often a key to a protein's function. During complex formation, each protein component will undergo a change in the conformational state, for some these changes are relatively small and reside primarily at the sidechain level; however, others may display notable backbone adjustments. One of the classic problems in the protein‐docking field is to be able to a priori predict the extent of such conformational changes. In this work, we investigated three protocols to find the most suitable input structure conformations for cross‐docking, including a robust sampling approach in normal mode space. Counterintuitively, knowledge of the theoretically best combination of normal modes for unbound‐bound transitions does not always lead to the best results. We used a novel spatial partitioning library, Aether Engine (see Supplementary Materials), to efficiently search the conformational states of 56 receptor/ligand pairs, including a recent CAPRI target, in a systematic manner and selected diverse conformations as input to our automated docking server, SwarmDock, a server that allows moderate conformational adjustments during the docking process. In essence, here we present a dynamic cross‐docking protocol, which when benchmarked against the simpler approach of just docking the unbound components shows a 10% uplift in the quality of the top docking pose.


| INTRODUCTION
Specific protein-protein interactions are both ubiquitous and essential within biological processes, ranging from immune system surveillance to tissue development and repair. Invariably, changes in the conformational state of the components of protein complexes do occur upon complex formation. 1,2 However, it is difficult to a priori predict the extent of such changes. Nevertheless, much research has provided guiding principles, the most notable for which is the concept of conformational selection and induced fit. 3 Conformational selection can be defined as a set of states, with likely varying degrees of stability, and therefore occupancy, which a protein may sample before complexation with a binding partner. Induced fit, on the other hand, can be defined as the influence that one binding partner may have on the other, often in the final stages of docking, where conformational states may be stabilized that are not possible for the free unbound components to sample. Although recent work has uncovered some interesting trends, 4 just how important the above two mechanisms are for the complex formation of any one particular protein receptor/ ligand pair remains difficult to estimate. In the context of fully automated docking servers, which often cannot capture the full extent of conformational changes displayed by unbound components, it would be desirable to be able to precalculate unbound conformations that are more predisposed toward their bound conformational states, and subsequently take such predictions into account within the automated docking protocol. One relatively successful protocol, in which the above has been partially implemented, is that of cross-docking.
A number of cross-docking methodologies have been reported 5-7 ; however, once a number of diverse receptor and ligand conformations have been selected, the docking process is typically rigid body. Here, we ask the question of how well might flexible cross-docking perform: that is, selection of viable alternative conformations for docking, and then employment of a docking methodology that enables limited flexibility in docking the selected conformations, here termed as "dynamic cross-docking." Clearly, many thousands of preselected conformations could potentially be cross-docked; however, most automated docking algorithms, with some degree of dynamic flexibility allowed during the docking process, would struggle to search such a large space within a turnaround time suitable to run as a publicly available server.
Moreover, sampling of the conformational states space is a tedious process as it depends on the single potential calculation time and the number of samples. With increasing numbers of normal modes, sampling boundaries, and granularity of sampling, the search task rapidly becomes intractable when utilizing standard highperformance computer (HPC) hardware, thereby demanding multicore parallelization and advanced knowledge in distributed computing.
As the importance of the given nontrivial normal mode is generally unknown when the bound structure cannot be compared with unbound, we decided to start with a naïve protocol-docking the conformations at the extreme magnitudes (in positive and negative directions) along the first nontrivial normal mode. This required nine docking runs per complex. After that, we investigated whether knowledge of the theoretically best combination of normal modes improves the docking results. Finally, we employed a method for sophisticated inputs generation, which involves sampling and analyzing the multidimensional normal modes space; here restricted to the first three nontrivial normal modes, but the approach is extendable to any number of dimensions. This enhanced approach requires only four docking runs per complex, but includes a generally expensive sampling step.
Therefore, we employed a novel spatial partitioning library to efficiently search the protein conformational states space to reduce time needed to produce a minimal set of conformations for cross-docking.
All protocols are compared on a set of 55 unbound receptor/ligand pairs, and one recent, difficult CAPRI target, for which the coordinates of the complex are available.

| Overview
The methodology presented here is aimed at identifying a plausible route to identifying the best conformation for each of a receptor and ligand pair to input to a docking server to enhance the quality of resulting docked poses. Our own protein-protein docking server, SwarmDock, 8 a server that allows moderate levels of conformational change, is used to cross-dock selected receptor and ligand conformations. Conformations for docking are selected according to three protocols, as outlined later.

| The data set
The 55 complexes from the docking benchmark version 5, 9 for which there are both experimentally determined unbound receptor and ligand components, as well as the experimentally determined complex, were used in this analysis. In addition, one difficult CAPRI target, T131 (6GBG), was also selected to test the methodology: a CAPRI case for which the SwarmDock server previously failed to return a correct docked pose. Prior to the methods described after the original Benchmark 5 was manually curated to repair missing atoms and residues, these data accompany the DYNACROSS decoy set for scorers and we have generated as a part of this work (see Supporting Information for details).

| Calculation of elastic normal modes
To investigate the potential range of conformational states, a receptor/ ligand pair may sample upon forming a complex we calculate from their unbound states as determined experimentally and their individual set of normal modes. Representing protein conformations as a combination of normal modes is a well-established methodology for both understanding protein-fold conformational state transitions and employing in protein docking pipelines. [10][11][12][13][14] More specifically, a full-atom elastic network normal mode calculation was performed on each unbound receptor/ligand pair with the program elNemo. 15 Throughout this work, references to the lowest modes indicate the lowest nontrivial modes: the first six trivial modes, with zero frequency that corresponds to rigid-body rotation and translation, are not considered.

| Calculation of the solvent-accessible surface area
To investigate how docked pose quality may have improved when normal modes were applied to the inputs, we calculated the solventaccessible surface area (SASA) using the FreeSASA 16 Python package (https://freesasa.github.io/).

| Dynamic cross-docking with SwarmDock
All receptor/ligand pair conformations selected for cross-docking, see later for selection protocols, were submitted to the SwarmDock server. 8 The server and underlying algorithm have been described previously. 8,17,18 Briefly, the server models any missing residues and atoms, reformats nonstandard residues, and minimizes the starting conformation. The core algorithm can be described as a modified version of the particle swarm optimization (PSO) algorithm, 19 which is employed to optimize the binding energy, as determined using the DComplex potential. 20 Each particle search vector consists of the position and orientation of the ligand, as well as normal mode coefficients, corresponding to the five lowest frequency nontrivial modes, to model the conformation of the receptor and the ligand, calculated by the elNemo program. 15 After each PSO iteration, the lowest energy member undergoes a local optimization. 21 Starting swarms are released from approximately 120 evenly spaced positions around the receptor, and swarms are run four times from each of these positions.
It is necessary to mention that the quality of docked poses generated by SwarmDock Server may sometimes be of lower quality than other leading methodologies, for example, ClusPro. 22 However, here we aim to evaluate "dynamic cross-docking," a general methodology that can be applied to any other docking approaches and thus support the wider docking community.
2.6 | Selection protocols 2.6.1 | P1: Naïve approach: extremes of the conformation space of each unbound receptor/ ligand pair Most of the conformational space between unbound and bound conformations can be accounted for within the first few lowest mode combinations of the unbound. 18 First, for each unbound set of coordinates, the maximum extent along each of the three modes separately (maximum mode amplitude) is calculated, in both negative and positive directions. Preliminary runs of this protocol indicated some unwinding of secondary structure elements, and while this is tolerable to an extent, indeed may naturally occur upon protein complexation, a high degree of such unwinding would clearly destabilize the binding partners. Therefore, it was deemed prudent, for both this protocol and protocol P3 described later, to monitor the amount of unwinding and set a threshold for its occurrence.
Hence, each extended conformation is minimized via the program CHARMM 23 and the secondary structure conservation (see later) calculated by comparing the initial and new conformation secondary structure fingerprints. Preliminary runs also suggested initial maximum mode amplitudes of 100 and −100 in the positive and negative directions, respectively, and provided a good balance between the energy of an expanded conformation and extent of secondary structure unwinding. If the structure is minimized within a secondary structure threshold, the maximum mode amplitude is doubled, otherwise, it is halved. The process is complete when the minimum and maximum amplitude ranges match. As discussed earlier, the secondary structure conservation is an additional constraint to ensure that protein folds do not notably unwind; the percentage of secondary structure for the new coordinates is checked to ensure that the percentage of secondary structure relative to the original unbound coordinates does not deviate by more than 10% for any of the three main secondary structure states, helix (H), beta strand (E), and coil (C), as calculated by the program DSSP 24 (all secondary states other than H and E were converted to C) and called via Biopython (https://biopython.org/).
In this protocol, only the lowest frequency mode for each receptor and ligand is searched, and the maximum extent (mode amplitude) is determined given the constraints described earlier.
The first of the nine runs consists of a normal SwarmDock run of docking the unbound coordinates of receptor (Ru) and ligand (Lu).

| P2: Theoretically best normal modes combination as input
The vector of coefficients, β, that represent the closest fit between unbound and bound receptor or ligand conformational states, can be formulated in matrix notation as follows: where M is an m × n matrix of m normal modes and n coordinates (unbound), and T is a vector of the difference in coordinates, after superposition, between bound and unbound conformations. For the full mathematical derivation, please see Moal et al. 18 As a benchmark, the theoretically best combination of the three lowest frequency modes for both receptor and ligand, see Table 1, is used to construct the input conformations to the server. For this protocol, there is no need to monitor unwinding of secondary structure elements as we are simply interested in the parsimonious route between unbound and bound conformations and the combination of mode amplitudes to achieve this. It should be noted that the SwarmDock server is not using any information on known bound structures. The above procedure is for benchmarking purposes only and to investigate whether knowledge of the best combination can be helpful in docking, which is of interest for development of in silico predictors from unbound structures.

| P3: Sophisticated sampling and selection of inputs
Here, we concentrate our efforts on mapping the conformational space explored by the first three lowest mode frequencies of each receptor, and each ligand, for every complex in our set of 56 targets. Once the maximum envelope has been calculated, separately for every normal mode, as described earlier, its volume was sampled at regular intervals along each axis by calculating the relative stability of each volume point, a point that represents the linear combination of the three modes. Sampling was performed with the distance-dependent empirical atom-atom potential program DFIRE2, 25 a pair potential program that can rapidly assign the relative energy of a conformation.
In order to allow for energetic space comparison, the DFIRE2 potential values were normalized and multiplied by 10.0, so the nor- It should be noted that even for just three modes, the sampled conformational envelope for each receptor and ligand pair is large.
We selected a single conformation from each envelope as described  Table 1) indicate, see Figure S1, that although there is a general trend for the smaller eigenvalues to correspond to larger mode amplitudes, a particular eigenvalue can map to a wide range of mode amplitudes, especially at the lower end of the eigenvalue range, thereby indicating that any guidance gleamed from eigenvalue ratios could only be considered weak.  (1); values highlighted in bold if the first amplitude for a receptor or ligand is higher than the following two. Finally, columns P0-P3 show the best obtained docked pose quality (incorrect, acceptable, medium, or high) for basic docking runs (from unbound; (Ru):(Lu) run from P1), naïve approach (P1), theoretical combination of normal modes input approach (P2), and sophisticated sampling approach (P3). When comparing with protocol P0, transitions with improving quality are highlighted in blue, and the ones worsening the quality are highlighted in red. The C++ code was compiled in release mode (contrary to the default cmake debug mode). "sched_setaffinity" was used to pin the processes to specific cores. The results are presented in Table S1.

| Sampling efficiency benchmark
A common property used to measure distributed systems is "speed-up" S N , defined as the time to calculate on 1 CPU/core (T 1 ) vs the time to calculate on N CPUs/cores (T N ) (Equation 2): The speedup of 1.0 means no change in calculation times, whereas the speedup below 1.0 suggests that there is no gain in using parallel processing, for example, due to overheads. Similarly, we can compare the time to compute off-and on-Aether on 1 CPU/core. Traditionally, the theoretical limit for the speedup on N cores was N (eg, 4.0 for 4 cores). However, various special techniques, both hardware-based like Hyper-Threading, and software based, like spatial partitioning, are capable of boosting speedup even higher.
Using Aether Engine on a single core, a maximum 4.37 speed up in sampling has been achieved; this was beneficial for proteins above 4000 atoms, as well as increased acceleration because of the engine's ability to add parallelism to previously unparalleled problems (like DFIRE2 here) and scale to any number of CPUs/cores/machines without any change to the code. The speed-ups on one core were

| RESULTS AND DISCUSION
The 55 complexes from docking benchmark version 5, 9 for which there are both experimentally determined unbound receptor and ligand components as well as the experimentally determined complex, were used in this analysis (see Table 1). Interestingly, as part of reporting this benchmark, four automated protein-protein docking servers were benchmarked as to their success in finding solutions, three ab initio servers SwarmDock, 8  Although we have entirely focused on designing protocols to better generate conformations that enrich the list of docked poses with higher quality models, rather than improve their final ranking, a preliminary analysis to determine whether we can identify the improved poses did not indicate significant improvements (data not shown).
This underlines that development of good ranking schemes continues to remain a major obstacle for server developers. Nevertheless, if a list of solutions does not contain high-quality docking poses, no form of improved ranking will be able to enhance the overall docking pipeline, so we are encouraged by our ability to generate higher quality poses.
3.1 | P1: Cross-docking extreme conformations from the naïve mapping approach Nine-times cross-docking runs (as described in the Materials and Methods section) were performed using inputs generated at the extreme magnitude of the lowest normal mode (in both the positive and negative directions) for both receptor and ligand. In terms of the best FNATs achieved across the board for the 56 targets, we observe approximately a 10% uplift comparing with P0 (see Figure 1). Even more encouraging, we observe that for 21% of targets, a top docking pose can be generated that improves upon simply docking unbound receptor and ligand pairs in a single SwarmDock run (as measured by the CAPRI scoring scheme; see Table 1).
As the SwarmDock server has a number of stochastic elements associated with its central PSO algorithm, such as the generation of random initial orientations for the particles in each swarm, exactly the same docking poses are not obtained for multiple runs with identical starting protein conformations. Therefore, on 10 of the targets selected randomly, ranging from rigid body to the more difficult cases, we tested whether running SwarmDock nine times on each of them, with the standard algorithm (ie, just R(u):L(u)), would also produce improved docking poses as has been noted. Analysis of these results showed that no improved docking poses could be obtained by this simple approach (results not shown). Taking the best FNATs achieved (across the 56 targets) into consideration, we observe a similar (10%) uplift as gained by P1 (when comparing with P0; see Figure 1). We also observed that for 21% of targets, a top docking pose can be generated, which is an improvement upon simply docking unbound receptor and ligand pairs in a single SwarmDock run (as measured by the CAPRI scoring scheme; see Table 1). Considering the CAPRI methodology (which assumes that medium grade structure is also acceptable), we see an uplift in quality

| CAPRI Target 131 (6GBG)
We applied the abovementioned conformational search and crossdocking protocol to one of the CAPRI targets (see Figure 2A): Target 131, a complex between CEACAM1, a cell surface protein receptor, and HopQ, a Helicobacter pylori adhesion protein. 32 This is a target for which SwarmDock could not find an acceptable solution during the CAPRI blind trials, and the best docking pose had 0.18% of native contacts. We found acceptable solutions from (Ru):(L-) run in protocol P1 and from (Ru):(L5) run in protocol P3 (see Figure 2B for P1 run).
The conformational sampling energetic maps for receptor and ligand for the whole sampled three normal modes space and with normalized DFIRE2 ≤ 5.0 are depicted in Figure 3. It should be noted that in normal modes space, the energy distribution is not symmetrical. It must also be stated that with our final ranking scheme, we were unable to rank this improved docking pose within the top 10 runs, indicating that a parallel improvement in our ranking scheme would be needed to take full account of improvements achieved in enhancing the quality of some of our docking poses.

| CONCLUSIONS
The CAPRI blind trials have set many challenges in the macromolecular docking field, ranging from docking of homology models to the construction of multibody complexes. Irrespective of the specific challenge, one clear message always feeds through in the assessment of results, and deterministic modeling of flexibility in the docking process remains an unsolved problem. 33 One of the most exacerbating aspects, particularly for those developing docking servers, is our inability to estimate the likely degree of flexibility for each target.
Counterintuitively, docking from the inputs, when the theoretically best combination of normal modes for the unbound-bound transition is applied, is not shown to be the best docking protocol. Interestingly, the better conformations were obtained for protocols P1 and P3 where a degree of perturbation for the input structures conformations was employed. It is hypothesized that a relatively strained starting conformation could give an entropic impetus for docking.
Here we show, in protocol P1, that by carefully searching the lowest  Table 1 for details). Arrows show the normal modes axes-red, green, blue for the first, second, and third nontrivial normal mode, respectively. The pictures were generated using VMD 26 time for each position in normal modes space) could be mitigated by applying spatial partitioning and multicore/multimachine computation.
Here, we used Aether Engine, which displayed consistent speed-ups as parallelism is added, thereby allowing one to benefit from the full amount of resources available. Therefore, the efficient preprocessing and docking pipeline could be built and run as a publicly available server, with no prior knowledge or access to supercomputing resources needed.
As an additional improvement for the future work, we consider applying Aether Engine to allow the Particle Optimization algorithm to be run on cross-docked structures in one simulation instead of running nine (P1), or four (P3), separate SwarmDock server runs performed here. This will require the communication between swarms with various input conformations (in relation to translational and rotational space) and separately considering the normal mode space.
Future work on scoring strategies would be a good extension to this work as improving the quality of poses in the docking pipeline could only lead to improved ranking results.