{'fraction_non_alphanumeric': 0.02389294891046818, 'fraction_numerical': 0.011879734374478593, 'mean_word_length': 5.1022195072286705, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Prediction of protein–protein complexes using replica exchange with repulsive scaling

The realistic prediction of protein–protein complex structures is import to ultimately model the interaction of all proteins in a cell and for the design of new protein–protein interactions. In principle, molecular dynamics (MD) simulations allow one to follow the association process under realistic conditions including full partner flexibility and surrounding solvent. However, due to the many local binding energy minima at the surface of protein partners, MD simulations are frequently trapped for long times in transient association states. We have designed a replica‐exchange based scheme employing different levels of a repulsive biasing between partners in each replica simulation. The bias acts only on intermolecular interactions based on an increase in effective pairwise van der Waals radii (repulsive scaling (RS)‐REMD) without affecting interactions within each protein or with the solvent. For a set of five protein test cases (out of six) the RS‐REMD technique allowed the sampling of near‐native complex structures even when starting from the opposide site with respect to the native binding site for one partner. Using the same start structures and same computational demand regular MD simulations sampled near native complex structures only for one case. The method showed also improved results for the refinement of docked structures in the vicinity of the native binding geometry compared to regular MD refinement.


| INTRODUCTION
Biomolecular binding and in particular protein-protein binding processes to form functional complexes are key elements of almost all biological processes. Knowledge of the three-dimensional (3D) structure of protein-protein complexes is a prerequisite for understanding its function. Experimental structure determination as well as prediction of protein-protein complex structures are also of significant interest for the rational design of drug molecules to influence biological processes. Computationally efficient docking algorithms are frequently applied to identify putative protein-protein binding geometries based on surface complementarity or simple pairwise interaction potentials. [1][2][3][4] Molecular docking, however, often largely neglects or only approximately accounts for the flexibility of the binding partners and interactions with the solvent. [3,[5][6][7] It is possible to include a moderate degree of flexibility using for example deformations in soft normal modes at reasonable computational costs. [8][9][10][11][12] In some approaches a refinement stage with side chain flexibility is performed, mainly focusing on interfacial rearrangements. [13,14] In addition, the evaluation of identified binding geometries is largely based on empirical scoring functions applied to single complex conformations neglecting conformational and orientational entropic contributions to binding. Ideally, molecular association should be simulated including full flexibility of both partners and accounting for surrounding water molecules and ions. [15] Molecular dynamics (MD) simulations are in principle well suited for investigating biomolecular association processes including full atomic flexibility. The methodology has already been used to refine potential binding geometries identified in the docking efforts. [16][17][18] In selected cases it is even possible to skip any initial docking but to use ultralong atomistic MD simulations and directly mimic the physical binding process. [19,20]  does not scale well with system size and another method, Hamiltonian REMD (H-REMD) might be more suitable [21] because one can specifically scale force field parameters affecting receptor-ligand interactions. One possibility is to linearly scale the Lennard-Jones (LJ) and electrostatic potential across replicas [22] or reduce the ruggedness of the energy landscape by introducing soft core potentials. [23] The latter method has shown promising results to refine complex geometries close to the native binding mode but do not effectively reduce the problem of trapped binding sites on the receptor surface. [23] Transient binding states in agreement with experiment could be recognized using replica exchange Monte Carlo simulations for three protein-protein complexes using a coarse grained representation of the molecules. [24] It is also possible to use metadynamics methods to reconstruct the free energy surface of association and dissociation of proteinligand systems by gradually adding biasing potentials that destabilize already sampled protein surface regions. [25] In the latter study the choice of only two collective variables (CVs) was enough to identify the binding site of four protein-ligand systems. In general, a higher number of CVs is necessary to completely describe the relative ligand-receptor position and orientation. [26] A larger set of CVs can be chosen using reconnaissance metadynamics that incorporates a self-learning algorithm that gradually pushes linear combinations of the CVs. [27] For a protein-ligand system this method was able to identify multiple binding sites of the protein.
In a recent study by Pan et al., reversible association and dissociation of five protein-protein complexes was observed using tempered binding MD simulations. [28] However, the binding and unbinding events were captured in still expensive computer simulations (simulation times of several hundred microseconds) on the special purpose machine Anton. [29] In tempered binding the interaction strength between the two solutes (but not within each protein partner) is used as tempering coordinate (instead of the total temperature as used in standard simulated tempering). As another alternative it is possible to add an explicit repulsive biasing potential between partners in a series of replicas (BP-REMD) that keeps the ligand and receptor at various distance intervals apart in higher replicas. [30] The higher replicas allow to keep some space between partner molecules and therefore result in fast diffusion. Upon exchange with lower replicas favorable binding sites can be rapidly sampled also in the reference replica. This method showed promising performance to specifically accelerate the search process for identifying ligand binding sites on protein surfaces under realistic conditions. [30] However, the method requires to calculate an ambiguity distance between all pairs of surface atoms of both partners which is computationally demanding and not well suited to run in parallel on many cores such as graphical processing units (GPUs).
In the present study, the possibility of increasing the repulsion between ligand and receptor by specifically increasing the pairwise effective van der Waals (vdW) radii and reducing the vdW attraction along the replicas in an H-REMD simulation is explored. It weakens not only the LJ contribution to binding interactions but also reduces the number of hydrogen bonds and electrostatic interactions due to an increased average distance between ligand and receptor. Hence, the biasing potential in the replicas allows the partners to rapidly dissociate from possible suboptimal binding sites to effectively search the protein surface. The method is promising for its simplicity of implementation and only requires adjusting parameters and can therefore be used with existing simulation software that runs on GPUs. The approach was tested on several protein-protein complexes of different sizes and types. In contrast to regular MD simulations it allowed the identification of near-native complexes even when starting far from the native binding region. In addition, we tested the approach for refinement of complexes starting from geometries in the vicinity of the native binding arrangement. In this case also a slightly better performance than regular MD simulations at the same computational effort was achieved.

| MATERIALS AND METHODS
For all atomistic simulations the Amber16 or Amber18 software packages, [31,32] were used employing the pmemd.cuda module for efficient calculations on GPUs. The ff14SB [33] force field was used together with an implicit water representation using the OBC Generalized Born (GB) model [34] (igb = 8 option in amber) involving an infinite cutoff.
2.1 | Simulations of protein-protein complexes starting far from the binding geometry MD simulations on protein-protein complexes in order to identify the native binding arrangement were started from an initial placement of one partner (termed ligand) at the opposite side of the second (receptor) protein with respect to the native binding site. Six complexes were considered for these simulations (pdb-id of complexes: 2oo9, 2cfh, 7cei, 2sni, 1gcq and 1syx, see also Supporting Information Table S1). In all cases the unbound protein structures were used for the simulations. The six complexes were selected due to the relatively small size from the docking benchmark 3.0. [35] We distinguish between receptor and ligand protein according to the assignment in the benchmark 3.0 [35] (typically the large protein partner is the receptor and the smaller partner is the ligand). The root mean square deviation (Rmsd ligand : Rmsd of the ligand after best superposition of the receptor with respect to the native complex structure) of the initial relative placement of the partner proteins from the native complex geometry was between 30 and 61 Å.
MD simulations were performed using the OBC (Onufriev, Bashford, Case) Generalized Born (GB) implicit solvent model [34] (igb = 8 option) and using an infinite cutoff radius for both the GB radii and nonbonded interactions. A Langevin thermostat with a collision frequency of γ = 5 ps −1 was used to control the temperature. The collision frequency is reduced relative to a more physical value of 50 ps −1 to reduce the apparent viscosity of the solvent and speed up sampling. [36] Equilibration of the start geometry was achieved after energy minimization (50 steps steepest descent followed by 1,500 steps conjugate gradient) and heating in three steps (each 12 ps) to 300 K with positional restraints of 0.05 kcal mol −1 Å −2 applied on the heavy atoms relative to the starting structure. Since we observed in long MD simulations for some proteins a partial unfolding positional restraints on the receptor C α atoms (force constant 0.05 kcal mol −1 Å −2 ) were also included during production simulations. Note, that such weak restraints allow still considerable backbone fluctuations and full side chain flexibility but prevent unfolding or large domain motions in the proteins. To prevent the ligand from diffusing too far away from the receptor, restraints between the center of masses (COM) of the C α atoms of the proteins were employed. The restraining energy was zero for COM distances below a certain threshold and increased quadratically beyond the threshold (force constant 1.0 kcal mol −1 Å −2 ) so that it prevents large receptor-ligand separation but still allows the ligand to dissociate from the receptors up to a certain distance. The COM distance threshold ranged from 27 to 50 Å for the different proteinprotein complexes and was larger than the sum of the largest center to surface distances for the two partner proteins (given for each protein pair in Table 1). It was chosen such that the ligand protein can reach every position on the surface of the receptor protein without violating the COM distance threshold. The mean difference of the applied COM distance restraint threshold and the native COM distances was 11 Å. For avoiding unfolding of the ligand protein additional intramolecular pairwise harmonic distance restraints between the C α atoms of the ligand protein (only distances between 5 and 10 Å) were applied, that prevented the ligand backbone from unfolding (force constant 2.0 kcal mol −1 Å −2 ) but allow full side chain flexibility. Note, that the above harmonic backbone distance restraints do not restrict the orientational and translation freedom of the ligand protein relative to the receptor protein.
In order to perform Hamiltonian replica exchange simulations (H-REMD), 16 replicas for each protein were generated with different Lennard-Jones (LJ) parameters for atom pairs involving atoms from different protein molecules (all intramolecular nonbonded parameters within were preserved). The intermolecular LJ potentials were scaled by a parameter d that adjusts the effective van der Waals radius and a factor e that changes the potential well depths (see next section for a detailed description). The following parameter set for d and e, with a smaller step size between the parameters close to the reference replica that increases in the higher replicas, gave the best results for protein-protein test simulations (see Table 2). For each replica, a short equilibration was performed for 32 ps with no exchange attempts. In the production run every 1,000 T A B L E 1 Simulation setups for each complex indicated by the PDB-id for the repulsive scaling H-REMD (RS-REMD) approach and the regular MD simulations MD steps an exchange between neighboring replicas was attempted, yielding a total simulation time (per replica) ranging from 340 to 845 ns (see Table 1).
Finally, starting from the same equilibration runs, 16 regular MD simulations with no H-REMD but different initial velocities (using the same restraints) were performed for comparable timescales as the H-REMD simulations (see Table 1).

| Refinement of individual protein-protein docking poses in implicit solvent
In addition to simulations starting far away from the native binding geometry, H-REMD and regular MD simulations were also performed for arrangements in the vicinity of the native complex structure obtained by an initial protein-protein docking run using the program ATTRACT. [37,38] The same set of structures and docking procedure as used in a previous study [39] were employed (see Supporting Information Table S2). Since the H-REMD method for refinement of docked complexes is computationally demanding the number of test complexes was limited to 20 complexes from the docking benchmark 3.0. [35] The docking was performed using a standard docking protocol on the unbound partner structures with the program ATTRACT. [37,38] The 300 top-ranked complexes were considered. Out of this set, the 50 models with lowest Rmsd to the native complex structure were used for further refinement using the RS-REMD or regular MD simulations. In order to refine the docking solutions atomistic replica exchange simulations in implicit solvent were performed (OBC model, [34] using the same conditions as described above) starting from the 50 docking poses of 20 protein-protein complexes. Energy minimization consisted of 2,500 minimization steps (400 steps of steepest descent, 2,100 steps of conjugate gradient). The systems were heated gradually in three steps of 15 ps to 300 K using a Langevin thermostat for temperature scaling. For each equilibrated pose eight replicas were generated with increasing bias for higher replica numbers of the intermolecular LJ parameters. As described above for the simulations starting from the opposite side of the receptor protein a parameter d adjusting the effective van der Waals radius and a factor e that changes the potential well depths was varied between the eight replicas (see Table 2).
Each replica was simulated for 0.5 ns with an exchange attempt every 125 steps amounting to 4 ns simulation time per pose. Intramolecular pairwise harmonic distance restraints between the C α atoms of each individual protein were applied (force constant 0. together with a COM distance restraint of interfacial C α atoms between the ligand and receptor (atoms with distances between 10 and 15 Å were considered) with a half parabolic shape (force constant 1.5 kcal mol −1 Å −2 ) that prevents full dissociation in the high replicas and shrinks the possible sampling space for these short simulations. The same simulation conditions and restraints were applied for regular MD simulations of each pose (no replica exchange and bias involved) of the same simulation time (4 ns) following a standard refinement protocol developed previously. [39,40] For evaluating the interaction energy a short MD simulation (30 ps) was applied on the reference replica of the REMD simulations followed by a minimization (500 steps of steepest descent, 2000 steps of conjugate gradient), which was also applied to evaluate the final structures from regular MD simulations. Finally, the minimized structures were scored by subtracting the potential energy of the partners from the energy of the complex. [39] To access the deviation of the refined structures from the native binding site the Rmsd ligand was calculated, the root mean square deviation of the ligand to the native ligand after superpositioning the receptor on the native receptor (only heavy atoms were considered).

| Refinement of a protein-protein docking ensemble in implicit solvent
Multiple docking poses were considered in a single REMD run to perform refinement simulations for each of the 20 protein-protein complexes. Only docking poses with a rmsd ligand above 10 Å (for the complex 7cei 8 Å was chosen, due to a lack of poses with large Rmsd) were considered as starting structures for the subsequent RS-REMD refinement.
Each of these poses was first scored based on the potential energy dif-

| LJ parameter scaling between partner molecules
The LJ interaction consists of an attractive part proportional to 1/r 6 and a repulsive contribution typically modeled by a term proportional to 1/r 12 . The parameters ϵ ij and R ij in the LJ potential determine the magnitude of attractive interaction and the effective (pairwise) van der Waals radius of the interaction between a pair of atoms of type i and j.
Typically, only the parameters between atoms of the same type, ϵ ii and R ii are used and one obtains parameters for pairs of different atom types using the Lorentz-Berthelot rules [41] : By defining new atom types, it is possible to specifically modify the LJ potential for interactions between the ligand and receptor without affecting the LJ interaction within one partner molecule or with the solvent. We used this possibility by adding an adjustable parameter d to an effective pair-wise van der Waals interaction between pairs of ligand and receptor atoms, One might also scale R ij by multiplying with a factor, but this would increase the effective radius of pairs of atoms by different amounts and may strongly distort an interaction interface. A change in the potential well depth ϵ ii is described by a factor e: As this factor enters multiplicative instead of additive the same relative scaling of attractive interactions of different pairs of atom types is possible. One subtle problem with increasing d, however, is that the number of atoms that can interact increases (illustrated in Figure 1).
This increases the total binding strength, even though any individual interaction might be weaker due to an e < 1. ϵ ij has to be decreased further to compensate for that effect. The number of atoms that can exactly fit into the energy minimum around one atom is proportional to the surface area of a sphere with van der Waals radius, which would suggest the following quadratic correction: The LJ potential minimum also gets wider linearly as R ij increases and more atoms can fit into the minimum along the radial direction, leading to a cubic correction: In the cubic case, the binding energy stays approximately constant so that the correction with Equation (6) Table 1). In the smallest test case (2oo9) the native binding site In all but one case sampling of near native arrangements was observed in the reference replica after $20-400 ns (see Figure 2 and Table 3). For the 2oo9 system, it took 55 ns and thus a bit longer than (by $100 ns, see Table 3) and the relative occupancy of the binding site was 78% in the reference replica, higher than the 47% observed when combining all regular MD results (see Supporting Information - Table S1). The RS-REMD simulation on 1syx explored near-native binding arrangements (Rmsd ligand < 10 Å) after 68 ns and thus more  Figure S1).
In all the cases for which the binding site was identified RS-REMD led occasionally to a very close agreement to the native structure with a lowest Rmsd ligand of $3 Å (see Table 3). In particular, the lowest Rmsd ligand using RS-REMD was closer than for regular MD also in those cases were both methods identified the correct binding site.
Thus, in these cases RS-REMD not only performs better than regular MD in the global searching process for the binding site, but also for local rearrangements at the binding site.
In only one case, 1gcq, the near-native binding geometry was not  Table S1). This is not the case for 1syx with a population of the near-native complex of $ 25% (still forming the largest populated cluster; Supporting Information Figure S4) but some alternative binding modes reaching similar population indicating similar binding affinity. The population of ligand placements at the native binding site is highest in the reference replica and decreases for the higher (larger bias) replicas (see Supporting  Information Figure S5 for the example case 7cei) due to the higher repulsive bias. Hence, an advantage of the RS-REMD technique relative to regular MD is that the near-native binding site can be identified by just looking at the population in the different replicas.

| Refinement of individual protein-protein docking poses in implicit solvent
Significant computational demand and simulation times are still necessary to reach near-native binding geometries using RS-REMD from distant initial placements. However, this corresponds to a worst-case scenario.
Instead, it is also possible to first perform a rapid protein-protein docking (not including solvent or partner flexibility) in order to first identify potential binding sites possibly not too far from the native binding geometry. In a second step short RS-REMD simulations are used as a refinement procedure to further improve the docking results. We first performed a docking run on a subset of 20 protein-protein complexes of the proteinprotein benchmark 3.0 [35] using the docking program ATTRACT and obtained 50 top ranked poses with different ligand deviations from the bound complex (Rmsd ligand < 25 Å) around the receptor protein (same as used in a recent study on protein-protein docking scoring [39] ). For each of these poses a short RS-REMD refinement (4 ns) was performed to refine the docking results. In order to limit the computational demand for this procedure an RS-REMD with eight replicas was performed.
Overall, upon RS-REMD refinement of all 50 decoys for the 20 test cases a slightly larger number of models (65%) with higher Rmsd ligand was observed compared to the starting structures (35% of poses had higher Rmsd ligand before refinement) ( Figure 4). Likely, because of the short simulation time, no Rmsd ligand improvements better than 13 Å were sampled. More important than the Rmsd ligand improvement of poses with a high deviation from the native binding mode is the refinement performance of the near-native models. For the docking model closest to the native binding site a smaller Rmsd ligand after RS-REMD refinement was observed for 13 complexes.
Also, the improvement in Rmsd ligand was higher, so an overall improvement of 0.49 Å was found considering the mean difference in Rmsd ligand before and after RS-REMD refinement of the closest to native pose (see Figure 4, left panel).
The RS-REMD scaling results are compared to a well-established atomistic refinement procedure [39,40] The two key poses were shifted by the mean scoring value of all poses and divided by the minimum scoring value in order to obtain comparable results for each protein. In 16 cases RS-REMD was able to identify the near-native binding placement (positive selectivity) and only in four cases the refinement approach resulted in a clearly negative selectivity (1ffw, 2oob, 1ak4, 2i25), identifying an incorrect binding site.
In summary, the RS-REMD refinement was able to improve in many cases the placements of the near-native poses. It overall performed slightly better than an established regular MD refinement procedure (at same computational costs) in terms of the selectivity in identifying the native binding site.  comparison to regular MD refinement (17%) (see Figure 6). For 14 of the 20 structures RS-REMD was able to capture the binding site, in some cases the population was highest at the binding site (see Figure 7). It points out, that even relatively short RS-REMD simulations can capture near-native binding geometries that were not already found in the initial docking search.
Performing only one refinement simulation starting from an ensemble of promising docking solutions and not refining every single pose individually can significantly reduce the computational demand.
Hence, the refinement starting from an ensemble of promising docking solutions can be considered as overall more efficient than refinement of every single pose.
Comparing the Rmsd ligand histograms of the refinement procedure to the histograms of the pure RS-REMD simulations (see Figure 2 and Supporting Information Figure S4)

| CONCLUSION AND OUTLOOK
A new H-REMD scheme is presented that includes a repulsive scaling potential (RS-REMD) between different protein molecules based on modification of the intermolecular LJ parameters. The bias requires a modification of the simulation parameter file but no changes in the underlying MD program are involved and full GPU support is possible.
F I G U R E 7 Histograms of the Rmsd ligand from the native structure for the reference replica of the RS-REMD refinement (magenta) (for all 20 protein-protein test cases in an individual figure) is compared to the Rmsd ligand histograms of the regular MD simulations (green). Both refinement procedures were initialized from poses that were not at the binding site [Color figure can be viewed at wileyonlinelibrary.com] The replica exchange scheme was applied and tested on three tasks that seek to identify the native binding geometry of protein-protein complexes using an implicit solvent model. First, RS-REMD allowed to sample near-native binding placements in five out of six example complexes, starting from a random placement far away from the native binding site. In contrast to multiple regular MD simulations, which were stuck mostly at locally stable sticky sites, these sticky sites were overcome through several intermediate steps in the RS-REMD. While the higher replicas sampled the whole receptor surface, the reference replica sampled locally favorable sites quickly until the native binding site was captured but depending on the case alternative binding modes were also still sampled. Although much less demanding than regular continuous (c)MD simulations still quite extensive sampling is needed for this approach that may limit its applicability. Our RS-REMD approach is designed to identify native binding geometries given the knowledge that the protein partners form a complex. For two proteins that do not bind in reality the approach will likely also suggest putative binding geometries. It might, however, be possible to predict the likelihood of complex formation for two protein partners by studying the population of associated states in the higher replicas (with larger repulsive bias) or from the lifetime and accumulation of complexed states in the reference replica. This will be subject of future studies.
In addition to starting from a worst-case scenario, we also used the approach for refining pre-docked poses. By applying a short RS-REMD run for each of the 50 poses of a benchmark set of 20 protein-protein complexes, it was possible to decrease the mean deviation from the native binding site. Moreover, the mean selectivity of identifying the native binding site according to a simple scoring function (based on the interaction energy) was increased in comparison to a regular MD refinement.
The simulation effort could be further reduced using a refinement scheme that associates each replica in the RS-REMD run with a different docking pose as starting structure. In contrast to the first refinement procedure, only one refining simulation had to be performed for each protein-protein complex. The population of the ligand protein partners near the binding site was clearly increased with RS-REMD beyond the result achieved by regular MD simulations. The benchmark set also contained difficult test cases (also reported in [39] ), where the identification of the native binding site was not possible in both refinement procedures.
Possible reasons are inaccuracies of the implicit solvent model that may not always favor correct complex structure relative to alternative arrangements. Explicit solvent simulations may help to solve this issue and will be tested in future studies. However, the probably slower diffusion and increase in number of particles will likely demand higher computational efforts. Another limitation of our setup is the inclusion of conformational restraints of the partner molecules with respect to the unbound conformations. This avoids any large-scale conformational change or unfolding of partners but in some cases may prevent conformational adaptations necessary for productive protein-protein complex formation. More global restraining methods like inclusion of backbone Rmsd restraints can help to overcome this issue in future efforts.
In principle, the RS-REMD biasing scheme can also be helpful to study folding/unfolding events or dissociation/association of parts of a protein structure. In such a case only the interactions of the selected part of the protein with other protein segments are scaled in the replica simulations.