Protein-protein interactions are involved in nearly all of the biological processes in human health and disease. Understanding the dynamics of binding and the structure of protein complexes at the molecular level can be instrumental in delineating biological mechanisms and developing intervention strategies.
Computational protein-protein docking provides a route to predict the three-dimensional structures of protein assemblies or complexes from known structures of individual monomeric proteins [1] .
Docking methods are tested in the blind prediction challenge known as the Critical Assessment of PRediction of Interactions (CAPRI) [2] , which in recent rounds pushed the field by including a wide array of target types such as transport proteins, higher order assemblies and host-virus interactions [3, 4] . Out of the 28 protein-protein targets evaluated in CAPRI over the past four years [4, 5] , predictors achieved high quality structures for 11 "easy" targets, defined as those with little backbone motion (unbound to Figure 1 : Performance of protein docking approaches on blind targets in CAPRI Rounds 38-46. [4, 5] Distribution of DockQ scores for the best model submitted by each predictor group (points) for each individual target (x-axis). DockQ measures a combination of intermolecular residue-residue contacts, interface RMSD, and ligand RMSD on a scale of 0 (incorrect) to 1 (matching the experimental structure) [5] . Targets are labelled by their CAPRI target number and, when needed, interface number (after the decimal). The targets are classified into rigid (easy) targets (high-homology monomer templates and under 1.2 A unbound-bound backbone motion, and flexible targets (poor template availability and/or over 1.2Å RMSD BU ). DockQ scores are color-coded by CAPRI model quality ranking: blue, high; green, medium; yellow, acceptable; gray, incorrect. Data graciously provided by Marc Lensink [4, 5] .
bound C α root mean square deviation (RMSD BU ) of less than 1.2Å [6, 7] ; Figure 1 ). The remaining 17 targets were categorized as "difficult" (RMSD BU over 2.2Å and/or poor monomer template availability).
For these targets, predictors only achieved acceptable quality in 8 of 17 targets (47%) and high quality in only 2 (12%) [4, 5] . Thus, the intrinsic flexibility of biomolecules still confounds the protein docking community at large.
In this review, we focus on the central docking challenge of capturing larger binding-induced conformational changes. We summarize progress by recent algorithms and frameworks, additionally augmented by growth in databases and computational power (CPU-and GPU-based). These new methods have achieved greater accuracy on more challenging targets and additionally yielded insight into binding mechanisms. We first present progress in binding site identification and then docking methods including molecular dynamics (MD) and Monte Carlo (MC) approaches, normal modes, and machine learning.
Together, these techniques have helped better explore broader regions of conformational space and more thoroughly evaluate the energy landscape to improve protein-protein docking. 
To reduce the complexity of the immense conformation space of flexible proteins, coarse-grained models are frequently used to reduce the degrees of freedom ( Figure 2 ). In the extreme, global docking approaches typically first treat protein partners as rigid bodies by restricting to six degrees of freedom (three rotational and three translational). A prime method to exhaustively sample the global 6D space is enumerating and scoring different rigid-body orientations on a dense grid. Approaches such as Clus-Pro [15, 16] , ZDOCK [17, 18] , PIPER [19] and HexServer [20] [21] . Relative to traditional FFT-based docking, FMFT accelerates calculations ten-fold [13]*. Another shape-based approach is geometric hashing, which indexes point sets or curves to match geometric features under arbitrary transformations like translations, rotations or even scaling [22] . Local 3D Zernike descriptor-based docking (LZerD), one of the top methods in CAPRI, projects 3D surfaces onto spheres to efficiently capture complementarity of protein surfaces [23] . Some rigid-body approaches exploit data from chemical cross-linking experiments [24] or small-angle X-ray scattering (SAXS) [25] to further improve discrimination of generated structures. These approaches provide fast, global exploration of the energy landscape, and in recent CAPRI rounds [4, 5] , many predictors incorporated these approaches as the first step to identify putative binding patches, and they supplement with other refinement tools to capture backbone flexibility.
Molecular dynamics (MD) is one strategy that is often used after grid-search or template-based approaches for refinement ( Figure 3 ) [26, 27] . Unbiased, all-atom MD simulations can provide a highresolution, time-resolved microscopic model of protein-protein interactions. MD calculates Newtonian trajectories using physics-based energy functions to simulate protein association and dissociation events.
MD use for protein docking has been limited because non-native local minima trap proteins, and dissociation is too slow [28] . Over the past decade, two new modifications to capture conformational changes are steered molecular dynamics (SMD) [29] , which utilizes external force constraints, and Markov sampling, which breaks a long MD simulation into multiple short trajectories [30] . To accelerate dissociation of protein partners at sub-optimal binding regions, Ostermeir et al. developed a Hamiltonian replica exchange MD protocol (H-REMD) for protein docking [31] *. In H-REMD, biasing potentials are based on the shortest distance between protein partner atoms (defined as "ambiguity restraints"). As the biasing potential and associated ambiguity restraints vary across replicas, associated protein partners in one replica are forced to dissociate in another. Pan et al. simulated long timescales in a global search space for a benchmark set of five targets on the special purpose machine Anton [32, 33] . Their "tempered binding" protocol updates energy function parameters throughout the simulation: a soft-core van der Waals intermolecular potential is scaled so that long-lived states are dissociated more frequently, improving the 
In contrast to MD approaches that target flexibility with Newtonian dynamics; Monte Carlo (MC) methods sample by random moves often followed by minimization (MCM) [40, 41] . MC allows a wide variety of conformational move types to sample diverse conformations. MC algorithms have emulated the kinetic binding models, namely key-lock, conformer selection (CS) and induced-fit (IF) mechanisms [42, 43, 44] . The CS model chooses protein backbones from a pre-generated ensemble, thus this approach has the advantage of docking one partner's conformations at a time. However, CS docking can fail if the ensemble is devoid of native-like backbone conformations [45] . docked from an ensemble of 10 structures). To diversify backbone conformations, the protocol generates monomer structures by three methods: (1) normal modes [46] (2) backrub motions [47] and (3) all-atom backbone refinement [48] . Further, to discriminate between near-native and non-native structures, they developed a more accurate coarse-grained energy function with 6-dimensional residue-pair data obtained 
Since intrinsic fluctuations in proteins contribute to conformational change, some docking approaches utilize harmonic dynamics to capture protein backbone motions [49, 50, 51] . Normal modes of vibration represent internal motions of a protein based on a Hookean potential between close residues. Normal mode analysis (NMA) is incorporated in docking approaches such as ATTRACT [52] , FiberDock [53] , SwarmDock [54] and EigenHex [55] . To mimic induced-fit, Schindler et al. developed iATTRACT [56] by moving interface residues in Cartesian coordinate space subject to NMA-generated harmonic potentials.
iATTRACT served as a refinement stage and improved the fraction of native contacts predicted by 70%. For targets with unbound to bound interface RMSD over 4Å, iATTRACT can achieve acceptable quality models [56] . Population-based methods such as particle swarm optimization (PSO) have also [57, 4] . Extending the swarm intelligence methods, the LightDock algorithm uses a "glowworm" swarm optimization to sample different backbone conformations in local regions of the protein surface with an anisotropic network model [58] . LightDock additionally uses multiscale modeling to combine all-atom and coarse grained scoring functions.
While normal modes have typically been used on individual protein partners prior to docking, Oliwa and Shen introduced the complex NMA in docking to also sample molecular complex fluctuations [60] . By calculating modes of an encounter complex, this approach focuses on the binding region as it reduces the dimensionality of the search space [61] . One of the problems of NMA is that higher frequency modes often distort protein bonds. To overcome this limitation, Frezza and Lavery developed the internal coordinate NMA (iNMA) approach to move in the torsion angle space, that is, with fixed bond lengths and angles ( Figure 4 ) [62] . With a reduced protein model in an internal coordinate space, they captured larger conformational changes from eigenvectors of low-frequency modes [59] **. iNMA can generate structures within 3Å of the bound state when starting from the unbound for 39% of single-domain and 45% of multi-domain proteins in their benchmark.
Although protein folding has been one prime focus of deep learning methods in biology (e.g., AlphaFold [63, 64] and RaptorX [65] ), in recent years, a few studies have explicitly addressed challenges relevant to protein docking [66] . Protein binding sites can be thought of as an information-rich molecular space that can be mined for elucidating protein interactions [67, 68, 69] .
One approach is to use this information to create score functions for use with traditional docking approaches. For example, Geng et al. used graph representations to train a support vector machine (SVM) on native and non-native protein complex structures to develop a scoring potential (GraphRank) to rank docked poses [70] . And iScore, composed of the GraphRank and HADDOCK [71] scores, achieved top performance in CAPRI scoring rounds (medium or high quality structures for nine out of 13 targets).
Other teams have used deep learning techniques to identify protein interfaces by extrapolating image recognition tools to protein structures. RaptorX-ComplexContact [69] uses a deep residual neural network trained on single-chain proteins to predict contacts between binding partners, achieving the top contact prediction scores in CASP [72] . Another approach is to characterize interaction environments.
Townshend et al. created "voxels," i.e., volumetric pixels with local atomic information for every protein surface residue, and with this 3D representation, they trained a deep 3D convolutional neural network (SASNet) on a curated database of bound protein complex structures [73] . Pittala et al. employed graph convolutions with the nodes representing the amino acid residues and edges connecting residues with a C β − C β distance under 10Å [74] . They placed geometric and chemical features on both nodes and edges and used a graph neural network to predict epitopes and paratopes in antigen-antibody interfaces. In a unique approach by Gainza et al., a geometric deep learning model (MaSIF) used molecular interaction "fingerprints" calculated using geometric and chemical features of protein surfaces [14]** (2). Their deep network was composed from geodesic convolutional layers, and they used it to predict binding sites, evaluate alternate docked interfaces, and assess likelihood of a given protein-protein interaction. Relative to conventional rigid docking methods on protein targets, MaSIF-search can perform ultra-fast scanning to identify true 'binder' with similar accuracy but significantly faster (4 CPU-minutes vs. 45 hours for PatchDock and 93 days for ZDOCK to evaluate a benchmark of 100 bound protein complexes).
In a study to explore how neural networks might be used to generate structures with considerable backbone motion, Degiacomi trained an autoencoder with conformations from MD simulations, compressing the protein motion into a low-dimensional latent space [75] *. By training with simulations of both closed (bound) and apo conformations of a target protein, the autoencoder generated an intermediate closed-apo conformation at 0.8Å RMSD [75] from the native state. However, when the autoencoder was trained only with open conformations, the generator could only create structures far from the closed state (over 4.2Å), limiting the utility of this approach for blind docking. In an approach suitable for blind cases, Cao and Shen developed a Bayesian active learning (BAL) model to quantify uncertainty in protein 8 structure quality, and then they extended their model to flexible protein docking [76] *. The Bayesian framework determines the posterior probability as it samples backbone conformations [60] . Flexibility is captured with low-frequency complex-NMA modes, and in principle it can be extended to higher frequencies that capture loop and hinge motions. Compared to ZDOCK [17] and PSO, BAL improves the interface RMSD of the near-native predictions by 0.5Å.
In conjunction with experimental data, docking has advanced a range of biological and health applications (e.g., Alzheimer's disease [77] , celiac disease [78] , SARS-CoV-2 [79] , influenza [80] , cancer [81] , and heart disease [82] , to name just a few). Over the past few years, docking success rates have improved on "difficult" blind prediction targets, but rates need to be higher for docking to be a reliable stand-alone tool in all cases. Clearly, a diverse and impressive array of tools has steadily advanced toward reliably capturing large conformational changes in protein docking. Docking will be even more impactful when the field finally overcomes this challenge. 