Alexander Sasse

Sjoerd J De Vries

Christina E M Schindler

Isaure Chauvot De Beauchêne

Martin Zacharias
email: zacharias@tum.de

Rapid design of knowledge-based scoring potentials for enrichment of near-native geometries in protein-protein docking

Protein-protein docking protocols aim to predict the structures of protein-protein complexes based on the structure of individual partners. Docking protocols usually include several steps of sampling, clustering, refinement and re-scoring. The scoring step is one of the bottlenecks in the performance of many state-of-the-art protocols. The performance of scoring functions depends on the quality of the generated structures and its coupling to the sampling algorithm. A tool kit, GRADSCOPT (GRid Accelerated Directly SCoring OPTimizing), was designed to allow rapid development and optimization of different knowledge-based scoring potentials for specific objectives in proteinprotein docking. Different atomistic and coarse-grained potentials can be created by a gridaccelerated directly scoring dependent Monte-Carlo annealing or by a linear regression optimization. We demonstrate that the scoring functions generated by our approach are similar to or even outperform state-of-the-art scoring functions for predicting near-native solutions. Of additional importance, we find that potentials specifically trained to identify the native bound complex perform rather poorly on identifying acceptable or medium quality (near-native) solutions. In contrast, atomistic long-range contact potentials can increase the average fraction of near-native poses by up to a factor 2.5 in the best scored 1% decoys (compared to existing scoring), emphasizing the need of specific docking potentials for different steps in the docking protocol.

Introduction

Protein interactions play a key role in almost all biological processes [1] [START_REF] Eisenberg | Protein function in the post-genomic era[END_REF]. While the number of protein-protein interactions discovered by experimental and computational approaches rises rapidly, the number of known complex structures lags behind [START_REF] Stumpf | Estimating the size of the human interactome[END_REF] [START_REF] Rolland | A proteomescale map of the human interactome network[END_REF]. However, experimental structural biology methods such as nuclear magnetic resonance (NMR) spectroscopy and X-ray crystallography have been used successfully to determine many of the unbound constituents.

Protein docking protocols aim to predict the structure of protein complexes from its unbound components. Docking protocols have been developed for single protein-protein, multiple protein, protein-peptide, protein-RNA and protein-DNA interactions [START_REF] Schindler | Fully Blind Peptide-Protein Docking with pepATTRACT[END_REF][6] [START_REF] Setny | A coarse-grained force field for Protein-RNA docking[END_REF] [START_REF] Zacharias | Protein-protein docking with a reduced protein model accounting for side-chain flexibility[END_REF]. State-of-the-art docking programs often achieve satisfactory results for sampling near-native docking geometries, particularly for cases with no or little structural changes in each constituent during complex formation [START_REF] Huang | Exploring the potential of global protein-protein docking: an overview and critical assessment of current programs for automatic ab initio docking[END_REF] [START_REF] Moal | The scoring of poses in protein-protein docking: current capabilities andfuture directions[END_REF].

Docking protocols usually can be divided into two stages: a sampling stage to generate an ensemble of possible complex solutions (decoys) and a scoring stage to select near-native complex structures from the sampled ensemble. However, some docking protocols model conformational adjustments during the binding process by introducing flexibilities during the sampling stage [START_REF] Vries | Flexible docking and refinement with a coarse-grained protein model using ATTRACT[END_REF] [START_REF] Bonvin | Flexible protein-protein docking[END_REF] and most approaches include a flexible refinement step after the initial rigid body search [START_REF] Schindler | iATTRACT: Simultaneous global and local interface optimization for protein-protein docking refinement[END_REF].

The selection of near-native docking solutions usually involves application of one or several combined scoring functions and diverse approaches have been designed to generate successful scoring functions in the last few decades [START_REF] Fiorucci | Binding site prediction and improved scoring during flexible proteinprotein docking with ATTRACT[END_REF] [START_REF] Pierce | ZRANK: reranking protein docking predictions with an optimized energy function[END_REF]. Physical approaches attempt to generate a universally valid scoring function represented by a model for the free energy. A common method is to linearly combine models for each energy contribution. The parameters for these models are fitted to experimentally determined energy values [START_REF] Vangone | Contacts-based prediction of binding affinity in protein-protein complexes[END_REF] [START_REF] Jorgensen | The OPLS [optimized potentials for liquid simulations] potential functions for proteins, energy minimizations for crystals of cyclic peptides and crambin[END_REF] [START_REF] Zhang | Determination of atomic desolvation energies from the structures of crystallized proteins[END_REF]. However, energy models frequently seem to be insufficient, since they often neglect or oversimplify more complicated terms such as entropic contributions or the change in solvation energy upon binding [START_REF] Van Dijk | Solvated docking: introducing water into the modelling of biomolecular complexes[END_REF]. Moreover, energy funnels of protein complexes appear to be narrow towards the native structure [START_REF] Tsai | Folding funnels, binding funnels, and protein function[END_REF], so that already small deviations on the interface cause large differences in energy. Hence, these methods are less accurate when scoring sampled complex geometries with rather low interface quality. Knowledge-based scoring functions for protein-protein docking like ITScore-PP [START_REF] Huang | An iterative knowledge-based scoring function for protein-protein recognition[END_REF], Sipper [START_REF] Pons | Scoring by intermolecular pairwise propensities of exposed residues (SIPPER):a new efficient potential for protein-protein docking[END_REF],

ProBinder [START_REF] Geppert | Protein-protein docking byshape-complementarity and property matching[END_REF], DECK [START_REF] Liu | DECK: Distance and environment-dependent, coarse-grained, knowledgebased potentials for protein-protein docking[END_REF], DARS [START_REF] Comeau | DARS (Decoys As the Reference State) potentials for protein-protein docking[END_REF], Tobi [START_REF] Tobi | Designing coarse grained-and atom based-potentials for protein-protein docking[END_REF] and ATTRACT [START_REF] Zacharias | Protein-protein docking with a reduced protein model accounting for side-chain flexibility[END_REF], are based on parameters extracted from the comparison of the scores between near-native or native structures and incorrect poses. One should note that such scoring functions are affected by the structures in the training set and recognize specific structural properties from the sampling algorithm used to generate the decoys [START_REF] Vajda | Sampling and scoring: A marriage made in heaven[END_REF]. Consequently, docking protocols generally use specific knowledge-based scoring functions for each class of ligands or step in the protocol [START_REF] Setny | Protein-DNA docking with a coarse-grained force field[END_REF][7] [START_REF] Fiorucci | Binding site prediction and improved scoring during flexible proteinprotein docking with ATTRACT[END_REF].

Composite scoring functions such as pyDock [START_REF] Cheng | pyDock: electrostatics and desolvation for effective scoring of rigid-body protein-protein docking[END_REF], RosettaDock [START_REF] Lyskov | The RosettaDock server for local protein-protein docking[END_REF], HADDOCK [START_REF] Dominguez | HADDOCK: A Protein-Protein Docking Approach Based on Biochemical or Biophysical Information[END_REF], Zdock [START_REF] Pierce | ZRANK: reranking protein docking predictions with an optimized energy function[END_REF] [31], FireDock [START_REF] Andrusier | FireDock: fast interaction refinement in molecular docking[END_REF], and FiberDock [START_REF] Mashiach | FiberDock: Flexible induced-fit backbone refinement in molecular docking[END_REF] usually use a linear combination of scoring terms to further improve their scoring by accounting for independent complex features. The linear weights for this purpose are determined by optimization or machine learning techniques on a selected set of decoys. In addition, several efforts were made to compare multiple docking protocols and their scoring functions to show future prospects of possible combinations [START_REF] Huang | Exploring the potential of global protein-protein docking: an overview and critical assessment of current programs for automatic ab initio docking[END_REF].

In general, each of the regarded functions or methods do not seem to work equally well for all cases in the docking benchmark. [START_REF] Moal | The scoring of poses in protein-protein docking: current capabilities andfuture directions[END_REF]. In particular, Vajda and Kozakov pointed out that any scoring function is substantially affected by the properties of the decoys in the training set [START_REF] Vajda | Sampling and scoring: A marriage made in heaven[END_REF]. Each sampling or refinement method creates a different ensemble of solutions with variant interface characteristics such as different dominant contacts or distances between atom-pairs. As a consequence, scoring functions need to be tailored to the targeted molecules or the sampling algorithm to improve their performance to predict near-native solutions.

In this work, we present a tool kit to develop knowledge-based scoring potentials for various decoy sets, sampling methods and problems in protein-docking. The GRid-Accelerated-Directly-SCoring-OPTimizing (GRADSCOPT) tool kit enables a choice between several atomistic and coarse-grained representations for various simple functional forms. The parameters of these potentials can generally be trained on any set of structures of protein complexes by linear regression (LR) models or a directly scoring-dependent Monte-Carlo algorithm (MC). Thanks to pre-calculations of potential-specific feature vectors, re-scoring and re-evaluation of the whole benchmark can be performed very quickly by simple vector multiplications to facilitate the fast testing of several scoring functions and their combinations.

As an application, we designed ten different scoring functions with four different potential forms in our newly defined grouped all atom (GAA) representation. We used unbound rigid-body proteinprotein docking decoys obtained from the protein-protein docking program ATTRACT [START_REF] Zacharias | Protein-protein docking with a reduced protein model accounting for side-chain flexibility[END_REF]. The performance was evaluated on three different scoring problems on an independent test set consisting of 77 complexes. Furthermore, the contributions of the most distinctive scoring parameters were examined for some of the potentials to explain differences and similarities in their scoring behaviour. Thereby, we found that scoring functions that find a large fraction of near-native solutions favour an increased amount of hydrophobic groups on the interface, especially aromatic rings and end-groups from other hydrophobic side-chains.

We also demonstrate that atomistic potentials trained on near-native solutions from ATTRACT show very different behaviour from potentials that were trained exclusively on the native structure.

In addition, we show that potentials using short range interactions for their scoring seem to be more adequate to detect the native bound form whereas potentials using long-range interactions are able to predict on average more near-native solutions in the ensemble of generated decoys.

Methods

Set up of scoring benchmark:

In this work, we considered 212 protein complexes from the protein-protein docking benchmark 5.0 to form a benchmark for the training of the parameters of our presented scoring potentials [START_REF] Vreven | Updates to the Integrated Protein-Protein Interaction Benchmarks: Docking Benchmark Version 5 and Affinity Benchmark Version 2[END_REF][START_REF] Hwang | Protein-protein docking benchmark version 4.0[END_REF].

Between 6,000 and 60,000 decoys (depending on the size of the complex) were sampled and scored by ATTRACT's rigid body docking for each complex from its unbound constituents. The scoring benchmark was divided into a training set of 135 complexes and a test set of 77 complexes for evaluation of the genuine performance. We took care that the test and the training sets include approximately the same fractions of hard, medium and rigid body cases to avoid a bias between these two sets in the scoring functions performance (see Supporting Information Table S1A andS1B).

Structural quality assessment:

To determine the quality of each structure in the ensemble of sampled solutions, the fraction of native contacts (Fnat), the interface root mean square deviation (Irmsd), the ligand root mean square deviation (Lrmsd) and the CAPRI-stars [START_REF] Hwang | Protein-protein docking benchmark version 4.0[END_REF] are calculated for all cases in the benchmark. To assess the quality of protein-protein complexes, we employed the CAPRI quality scheme which was defined by the community wide experiment on the Critical Assessment of PRediction of Interactions (CAPRI). It distinguishes acceptable (*: IRMSD < 4 Å and 0.1 < Fnat< 0.3), medium (**: IRMSD < 2 Å and Fnat > 0.3) and high (***: IRMSD < 1 Å and Fnat > 0.5) quality docking solutions. Generally, acceptable, medium and high quality solutions are summarized under the term near-native solutions to define geometries that are close to the native complex.

Calculation of potential-specific feature vectors:

To ensure fast enumeration of the parameters and also rapid re-scoring afterwards, potential-specific feature vectors are generated in advance for each solution in the decoy sets. The content of these feature vectors is generated based on first, the coarse-grained beads or atom types λ assigned to the structures and second, the potential form. All feature vectors contain the sum over potential-specific attributes for each contact-type between two atoms of type A and B.

Our tool kit allows to generate scoring potentials of different functional forms including step potentials, van-der-Waals (vdw) potentials and potentials that are based on atomistic buried surface areas (BSA). For step potentials, the feature vectors contain the number of contacts for each decoy in the desired ranges of the steps, thus single step potentials can also be referred to as contact potentials. For the BSA-potentials, the buried surface area for each coarse-grained or atom type is calculated in advance, by the rolling-probe algorithm using a water radii of 1.4Å [START_REF] Richmond | Solvent accessible surface area and excluded volume in proteins: Analytical equations for overlapping spheres and implications for the hydrophobic effect[END_REF]. The vectors for the vdw-potentials contain the sum over all distances between two atoms for each contact-type (A,B) to the power of -8 for the repulsive and -6 for the attractive part (see equation 1).

E vdw =⃗ α⋅ ( ∑ r 1,1 -8 ⋮ ∑ r λ ,λ -8 ) -⃗ β⋅ ( ∑ r 1,1 -6 ⋮ ∑ r λ ,λ -6 ) ; ⃗ α=⃗ ϵ⋅ ⃗ σ 8 , ⃗ β=⃗ ϵ⋅ ⃗ σ
6 Equation 1The dimension or the number of parameters in the feature vectors is determined by the selected molecular representation. By default, the tool kit offers to use three atomistic representations and the coarse-grained representation used in the ATTRACT docking engine [START_REF] Fiorucci | Binding site prediction and improved scoring during flexible proteinprotein docking with ATTRACT[END_REF]. For the atomistic representations, the atoms of incomplete side-chain residues are rebuilt and atom-types can be assigned according to the optimized potential for liquid simulations (OPLS) [START_REF] Jorgensen | The OPLS [optimized potentials for liquid simulations] potential functions for proteins, energy minimizations for crystals of cyclic peptides and crambin[END_REF], according to Tobi [START_REF] Tobi | Designing coarse grained-and atom based-potentials for protein-protein docking[END_REF] or to our newly defined grouped all-atom model (GAA). The GAA representation defines 27 atom types based on the chemical character of each group of amino acids, non-polar, polar, aromatic, positively charged and negatively charged (see supporting Information, Table S1).

Based on these vectors of each ensemble of solutions, the scoring and re-scoring can be realized by simple and thus quick vector multiplications between the generated potential parameters and the feature vectors of each docking pose (see equation 1 for vdw potentials).

Parameter training:

To train the parameters of the various scoring potentials, directly scoring-dependent Monte-Carlo Annealing as well as linear regression models are applicable. The Monte-Carlo Annealing method generates potential parameters by optimizing one of the various default target functions which are based directly on the comparison between the scoring of near-native and false solutions of each complex. For this purpose, quality-weights are assigned to the training decoys, either to distinguish between near-native and incorrect solutions or, depending on the objective, to allocate more influence to higher quality solutions in the training set. As a quality-weight for instance, one can consider the fraction of native contacts (Fnat), the interface root mean square deviation (Irmsd), the ligand root mean square deviation (Lrmsd), the CAPRI-stars [START_REF] Hwang | Protein-protein docking benchmark version 4.0[END_REF] (see above) or other alternative assessments, such as DockQ, a recently developed continuous docking quality measure that avoids border effects in its classification scheme[todo ref].

To account for diverse quantities of near-native structures among the complexes, the quality weights are normalized for each complex before optimization. The pre-implemented target functions t are calculated by the sum over the number of all complexes Nc and the number of all their training decoys Nd of the product between two weights w r and w q that depend on the score and the quality of each decoy i of complex c. The values of the ranking weights can either be set to rise linearly or quadratically from the last position in the decoy set or be assigned after another functional form, depending on the desired objective for the scoring function, such as enriching a fraction of nearnatives before refinement or a final ranking.

t ( ⃗ E )= ∑ c N c ∑ i N d w q i,c ⋅w r ( E i ) Equation 2
In order to perform simulated annealing in parameter space,in each step s the protocol consists of:

-changing a randomly picked parameter by a constant or pseudo-temperature T dependent value -re-scoring and re-ranking the decoys for each complex -calculating a new value for the target-function ts(E)

-accepting this step by the temperature dependent probability p s accept which is derived from the Metropolis criterion (Equation 3).

p ac c e pt s =mi n(1, e t s -t s-1 T ) Equation 3
After each step, the temperature is decreased by a selected annealing curve down to 0.1% of the starting temperature in a defined number of steps. In this process, the pre-calculated feature vectors make it possible to re-score millions of decoys per second in several steps on a single computer processor.

On the other hand, the linear regression protocol simply fits the parameters on the feature vectors to the negative values or their negative reciprocals of the structural assessments. The resulting potentials assign higher negative values to solutions of higher quality. For that purpose, various regression methods can be considered from the scikit-learn library [START_REF] Pedregosa | Scikit-learn: Machine Learning in {P}ython[END_REF], such as ordinary least squares regression, non-negative least squares regression, support vector regression, robust regressions or Bayesian-Ridge regression. Here, we only show potentials that were generated by ordinary and non-negative least squared regression. Generally, the different linear regressions differ between their cost-functions and some approaches might be more valuable in order to generate simple force-fields on a set of experimentally determined energy values, for instance to predict binding energies or affinities of specific protein complexes [START_REF] Giguère | Learning a peptide-protein binding affinity predictor with kernel ridge regression[END_REF].

Performance Assessment:

For a proper performance evaluation of the generated functions, several simple assessments are performed on an independent test set. We considered the percentage of complexes for which at least one near-native solution can be detected, the average fraction of near-native solutions for each complex and the probability of finding the native structure in the set of generated decoys. The average fraction of near-natives compared to the probability to find one near-native can be used to estimate the specificity of the scoring function towards certain near-native structures in the decoy sets. Finding the native structure within a set of sampled decoys represents an artificial problem and is only used to evaluate the methodology.

Discriminating scoring contribution:

To gain additional insights into the particular scoring performance of BSA-potentials and steppotentials, we regarded the most discriminating parameter contributions dp to the scoring.

Analyzing the contribution of each parameter allowed us to evaluate differences and similarities between their scoring. As the discriminating parameter contribution for an interaction between atoms of type A and B, we defined the product between the parameter σ AB withthe normed absolute difference between the average number of contacts nc (or size of the buried surface area) of nearnative (n-nat) and incorrect (incor) solutions.

d p A B =σ A B |⟨n c A B n-n a t ⟩-⟨n c A B i n c o r ⟩| ⟨n c A B t ot ⟩ Equation 4
By this regard, one examines the average contribution of each parameter by the distinction between near-native and incorrect solutions, taking into account deviating occurrences of each contact-type in total and between incorrect and near-native solutions. Instead of regarding the absolute difference however, we plotted the contribution p of near-native, native and incorrect solutions separately to detect false positive or false negative contributions, for which the parameter is positive/negative but the average number of contacts is higher in near-native/incorrect solutions. More generally spoken, the sign of these parameters correlates with the average distribution of contacts between near-native and incorrect solutions and points out a difference between potentials generated by our methodology and statistical potentials.

p A B nat =σ A B ⟨n c A B nat ⟩ ⟨n c A B t o t ⟩ Equation 5
The discriminating contribution dp of each parameter is simply given by the difference between the contribution p for incorrect and near-native solutions. The contribution to native complexes is shown as well to see possible differences between the scoring of near-native and native complexes.

Detailed training recipe for the generation of example scoring potentials:

We generated ten potentials in four functional forms: five by linear regression (LR) and five by the Monte Carlo Annealing (MC) method. Two long-range step potentials are based on atom-atom contacts between 0-10Å (MC_gaa_10, LR_gaa_10). Further, two short-range step potentials (MC_gaa_4-6_nat, LR_gaa_4-6_nat) consider two steps between 0-4Å and 4-6Å. The four vdwlike potentials ("_vdw") use the distance r to the power of -8 for the repulsive and -6 for the attractive parts (soft Lennard-Jones potential), respectively. Since the decoys were generated in a coarse-grained model from ATTRACT, clashes (too close contacts or structural overlap) may result from the change to an atomistic resolution. Hence, for the atomistic vdw-potentials, the contacts between the receptor and the ligand in the range of 0-2 Å were shifted to 2 Å to reduce their influence on the potentials. The potentials based on the buried surface areas consider only heavy atoms (type("_BSA")).

All potentials created by the MC approach used an adaptive search for parameter optimization, which altered the parameters by a factor d that started at 1. The temperature was decreased in 300,000 and 100,000 steps respectively from 50 down to 0.05 by a "ziczac" annealing scheme in which the temperature fluctuates between two decreasing exponential functions with a sin²function. When a convergence criterion was fulfilled in the last 2,000 steps, the search was stopped automatically. For the vdw-potentials, the parameter space was restricted. Sigma was kept between 1.5-6 Å and epsilon between 0-50 (see Equation 1). As target-functions for the potentials created by MC, a linearly increasing position weight was used for MC_gaa_4-6_nat, MC_gaa_BSA, MC_gaa_vdw and MC_gaa_vdw_nat. For MC_gaa_10 a quadratically increasing weight was used instead.

For the potentials from the linear regression approach, the fraction of BSAs and contacts was taken as features rather than the absolute number in order to prevent over-fitting on complexes with large interfaces. The potentials LR_gaa_10, LR_gaa_4-6_nat and LR_gaa_BSA were created by an ordinary least squares regression. For the vdw-potentials a non-negative least squares fit was used to generate only positive values for the parameters α=σ⁸ε and β=σ ε. ⁶ For all potentials, 5 parameter sets were created by leave-one-out cross-validation. Each parameter set was examined on its validation set to exclude the possibility of potential over-fitting. For step and BSA-potentials, the scaled parameter averages of the five generated sets were taken as our final set. Therefore, each of the five parameter sets was scaled by dividing the parameters through their standard deviation before taking their average. For the vdw-potentials, the set of parameters with the highest performance in the validation set was considered to be a consistent choice.

Results

Generating knowledge-based scoring potentials with the GRADSCOPT tool kit involves the following steps (Figure 1). First, a benchmark is set up with a sampling protocol generating an ensemble of decoys for each complex in it. The benchmark is then divided into a training and test set of complexes. Secondly, atom or coarse-grained residue types are assigned to the 3D structures of the receptor and the ligand. According to this representation and the form of the desired interaction potential, potential-specific feature vectors are calculated for the generated decoys (see Methods subsection Calculate potential-specific feature vectors). Subsequently, the parameters of the potential are trained on a subset of decoys from the training complexes by a directly scoringdependent Monte-Carlo annealing algorithm (MC) or by linear regression (LR). Finally, the whole benchmark is re-scored using the feature vectors of the whole ensemble, and afterwards the scoring performance of the generated potential is evaluated on a training-independent test set. This procedure can be performed in parallel or sequentially to generate several distinct scoring potentials in order to find the best suited variant.

We applied the GRADSCOPT tool kit to design ten different scoring functions based on three different potential types for unbound rigid protein-protein docking with ATTRACT. All ten scoring potentials are based on the GAA representation of the partners (using the 27 GAA atom types, see Methods and Supplementary Information, Table S2). In order to check the dependence of the potentials on the type of optimization method, the parameters were either optimized by the MC or the LR method. (see Table 1, further details on the optimization in Table S3).

The forms of the presented potentials can be considered as general representatives for the types of scoring functions typically used in the protein-protein docking field: we generated two long-range step potentials with core-core distances between 0-10Å (MC_gaa_10, LR_gaa_10), two short-range step potentials (MC_gaa_4-6_nat, LR_gaa_4-6_nat) with two steps between 0-4Å and 4-6Å, four vdw-like potential ("_vdw"), and two potentials based on the buried surface areas of each atomtype("_BSA") (see Table 1 for an overview). To prevent over-fitting, all potentials were created by 5-fold cross-validation and scaled parameter averaging and evaluated on an independent validation set. For the vdw-potentials, the parameter set with the highest performance in the validation set was considered as an appropriate choice. A detailed overview of the parameters used for the generation of the ten scoring potentials is given in Table S3 (see also "Detailed training recipe for the generation of example scoring potentials" in Materials and Methods). The parameter files of the generated potentials are included in the tool kit.

In order to investigate the dependence on the structures, defined as correct during training, we trained our potentials to consider all near-native solutions as successes by using the Capri-stars as their quality-weight, or by exclusively giving weights to artificially inserted native structures. The latter potentials were designed only for the case of step potentials (MC/LR_gaa_4-6_nat) and for continuous vdw-potentials (MC/LR_gaa_vdw_nat). It must be pointed out that in a practical docking algorithm, there is a zero probability of sampling the native structure from its unbound constituents and in addition, only docking geometries with a significant deviation from the native structure are typically included in the set of solutions. Small deviations between the unbound and bound structure in backbone and side-chain atoms prevent the sampling algorithm from finding the native structure or approaching the native complex structure closely without causing clashes (clashes with the unbound structure at the native position are shown in Figure 2b). Instead, structures with different contacts at the interface are usually formed; this includes non-native contacts and different distances between atoms that are in contact in the native complex. The deviation of acceptable, medium and high quality docking solutions relative to the native structure is illustrated in Figure 2a. Thus, as also emphasized by Vajda et al. [START_REF] Vajda | Sampling and scoring: A marriage made in heaven[END_REF], to generate a practically useful scoring function for any molecular docking protocol, training should be performed on complex structures sampled from their unbound constituents.

The performances of the generated protein-protein docking scoring potentials were evaluated on a separate test set (77 complexes) by three different assessments: (i) the probability to detect at least one near-native solution; (ii) the probability to detect the inserted native structure; and (iii) the average fraction of near-native geometries in a subset (e.g. enrichment). The performances in these tests were compared to ATTRACT's coarse-grained scoring [START_REF] Fiorucci | Binding site prediction and improved scoring during flexible proteinprotein docking with ATTRACT[END_REF] and to the scoring of Tobi's atomistic short range step potential [START_REF] Tobi | Designing coarse grained-and atom based-potentials for protein-protein docking[END_REF]. As a further control, the performances for near-native structures were also compared to a random scoring.

Identification of the native structure

In order to test the performance of our methodology to generate high quality scoring functions, we first tested the capacity of the designed scoring potentials to identify the native complex structure among all other decoy structures (Figure 3, Supplementary Information Table S4b). We expected that a scoring function that was trained entirely on native structures could easily solve this problem, since native interface contacts are typically much better aligned than contacts in sampled nearnative decoys. Indeed, our potentials trained exclusively to distinguish the native structure yielded very impressive results in predicting the native structure in the top-ranked 10 poses for 88% (twostep potentials) and 91% (vdw-potentials) of the cases (Figure 3). The designed potentials closely approached the performance achieved by the scoring function of Tobi [START_REF] Tobi | Designing coarse grained-and atom based-potentials for protein-protein docking[END_REF]. Other potentials not trained to identify the native structure showed much less impressive results (35-62% top 10 success rate, see also Supporting Information Table S4).

Identification of near-native structures

In a second evaluation, we tested the ability to identify near-native (that is at least Capri one-star quality) docking solutions among the decoys in the test set (Figure 4, Supplementary Information Tables S5). Note, that the decoy set did not include the native complex but only docking solutions obtained from systematic docking using ATTRACT. This corresponds to a realistic protein-protein docking experiment. All optimized potentials scored far better than random. The two 10Å step potentials outperformed ATTRACT by detecting a near-native solution in the top 10 for 31 % and near-native solution for 25 % of the cases in the test set. Considering the top 100, other functions such as LR_gaa_vdw and MC_gaa_vdw also showed good performance, identifying structures for 62% and 57% of the complexes, respectively. Interestingly, the two two-step potentials trained on native structures, MC_gaa_4-6_nat and LR_gaa_4-6_nat, performed very well, predicting a nearnative structure in the top 100 for 65 % and 68 % of the cases. In contrast, the vdw-potentials that were trained exclusively on the native structure (and showed great performance in identifying the native complex, Figure 3) lagged far behind these results with only 27 % and 30 % of cases with a near-native solutions scored among the top 100. The comparison between the results for near-native and native structures indicates that vdw-potentials are more sensitive towards certain distances between contacts of the training structures than step potentials. Therefore, the two-step potentials may be able to score near-native solutions well whereas vdw-potentials are unable to show a sufficient scoring. Also here, the results for the training sets were similar to the results in the test set (see Supplementary Information Table S5).

Enrichment of near-native solutions

To improve the quality of docking poses, many groups consider a subsequent refinement of a small fraction of their sampled structures. Refinement protocols allow atoms to adjust to the right position on the interface by introducing more flexibility than in the initial sampling. The chance to generate a high quality structure by refinement increases with the number of near-native structures considered. Therefore, we looked at the average fraction of near-native solutions found in the best scored 0.1%, 1%, 2% and 5% of the decoy sets.

We found that the long-range and the BSA-potentials outperform the standard ATTRACT score by a factor 2.5 and 2 respectively, for the fraction in the best 1% (see Figure 5. and Supplementary Information Tables S6b). On average they predicted 39 % and 33 % of all generated near-natives whereas ATTRACT only scores 17 % in the top 1 %. The vdw-potentials that were trained on nearnatives (LR_gaa_vdw and MC_gaa_vdw) predicted slightly more near-natives as ATTRACT with about 50 % in the top 5 % compared to 39 %(out of 6,000 to 60,000 decoys). Thus, MC_gaa_10 predicts on average as many near-natives in the top 1 % as ATTRACT in the top 5 %. Again, most potentials trained to find only the native structure performed worse: the potential by Tobi and the vdw-potentials identified on average less than 26 % of the near-native structures in the best scoring 5%, while the two-step potentials still placed 48 % and 54 % of the near-natives in the top 5%.

Comparing the long-range step and the BSA-potentials, we observed that potentials created by Monte-Carlo Annealing worked slightly better than their counterpart from ordinary least squares fit.

This might result from the target function in the MC algorithm that aims at improving the score of all near-native structures in the decoy set. The results further indicate that compared to atomistic or coarse-grained vdw-potentials, simple long-range step potentials and BSA-potentials are more likely to identify structurally diverse near-native complexes. Both our generated atomistic potential and the coarse-grained vdw-potential from ATTRACT seem to be more specific towards structures of higher quality, since contacts on the interface have to be within a certain distance in order to achieve a significant score. Hence, vdw-potentials generally do not seem to be a good choice to deal with the diversity of the different near-native solutions of each complex.

Discriminating parameter contributions to near-native enrichment for BSA-and long-range step potentials

After the designed long-range step potentials and the BSA-potentials achieved impressive results

for the fraction of near-native solutions, we wanted to investigate their scoring in more detail. We analysed the 20 most positive and negative discriminating scoring contributions (dp) for the step potentials along with all the scoring contributions of the BSA-potentials (see Discriminating scoring contribution). We looked at the difference between the average contribution p of each parameter of near-native (green) and incorrect solutions (red) to estimate its distinctive power in the scoring. By multiplying the parameters with their normed features, we accounted for the deviating occurrences of each contact-type in general and between incorrect and near-native solutions. First, we considered the discriminating scoring contributions for the 23 parameters of the BSA-potentials, MC_gaa_BSA (Figure 6S1a) and LR_gaa_BSA (Figure S1b). We found both potentials mainly predict near-native structures on increased BSAs of (I) C-rings and CH2 groups from aromatic side chains, (ii) CH3, CH2, CH, and S groups from nonpolar residues and (iii) nitrogen from positively charged amino acids as well as on reduced BSAs of CH2 groups from negatively and positively charged amino acids. The two potentials seem to deviate only slightly in the strength of contributions from less discriminating atom-types, for which the near-native and incorrect contributions are almost equal.

When we looked at the most negative discriminating parameter contributions of the 10Å step potentials, MC_gaa_10 and LR_gaa_10 (Figure 6a, S2a), we detected similarities to the BSApotentials as well as between the two potentials themselves. Eight and seven contributions respectively were between atom-types that were also most significant for the BSA-potentials.

Seventeen and fifteen out of the 20 most favourable contributions involved at least one group from a hydrophobic residue. The two potentials still showed major differences in the 20 most positive parameter contributions (Figure 6b, S2b): The MC_gaa_10 potential included eleven contributions which involved at least one CH2 group of charged amino acids, which represents a penalty for charged residues at the interface. The LR_gaa_10 potential only includes two of that kind.

Additionally, the LR_gaa_10 potential indicates eleven strongly false positive contributions whereas MC_gaa_10 shows only four. False negative/positive contributions were defined as contributions with a negative/positive parameter value but their average number of contacts is higher in incorrect/near-native structures. However, when we changed the signs of parameters with false-positive or false-negative contributions, the overall scoring performance to enrich a subset with near-native structures got worse (data not shown), indicating a difference of our potentials to statistically derived potentials.

Comparing the average scoring contribution of each parameter between native complexes (blue), incorrect solutions and near-native solutions, we found scoring contributions to native and nearnative structures to be very similar. Nevertheless, the difference to the average contribution of incorrect solutions is generally larger, which may be one reason for the better scoring of native structures. One must keep in mind that we looked at average numbers of contacts in a distance of 10 Å or at BSAs. We would expect to see even larger differences in closer ranges, which could explain the great performances of short-range potentials for native structures.

Discussion

The performance of the generated scoring potentials for protein-protein docking showed that both our approaches were able to rapidly create high quality scoring potentials. All our potentials worked significantly better than a random scoring; they even outperformed or competed with two state-ofthe-art functions in all three presented assessments.

The continuous vdw-potentials performed extremely well in scoring the native structure but poorly for enriching near-native docking solutions. Very similar results were found for the popular Tobi score. These results supported our idea that vdw-potentials are extremely dependent on the distances of interface atoms in their training structures and hence are biased towards certain nearnatives or natives in the decoy set.

The scoring performances of our generated potentials for at least one near-native solution are in general comparable to the results obtained from ZDOCK [START_REF] Pierce | ZDOCK server: interactive docking prediction of protein-protein complexes and symmetric multimers[END_REF], ClusPro and SwarmDock [START_REF] Torchala | SwarmDock: a server for flexible protein-protein docking[END_REF]. For example, our generated potential MC_gaa_10 placed a near-native structure for 13 % of the cases at rank 1, for 31 % in the top 10 and for 61 % in the top 100 compared to 10 %, 36% and 65 % for SwarmDock.

Furthermore, we showed that long-range step potentials and simple potentials based on the atomistic buried surface areas were able to detect on average up to 2.5 times more near-native structures in the top 1% than the standard ATTRACT score. Therefore, these potentials seem to be well applicable for selecting a subset of structures for a subsequent refinement in the ATTRACT or other docking protocols.

Interestingly, the analysis of the parameter contributions for these two types of potentials indicated that not only specific contacts between groups from aromatic and non-polar residues were favoured by the step potentials but also contacts between these groups with polar and charged residues or backbone atoms. The preference of these contacts by long-range step potentials may account for the increased presence of these groups in protein-protein interfaces in general, as also detected by the BSA-potentials. This seems to be helpful to select many near-natives, since interfaces of nearnatives with lower quality may not be aligned well enough to select them exclusively on chemically specific contacts.

As a result, we suggest that it would be beneficial to use different scoring functions before and after a refinement. For initial scoring, we suggest that unspecific long-range step or BSA-potentials should be applied for an enrichment of near-native structures. In contrast, short-range or atomisticdistance dependent potentials, such as our vdw-potentials and our multiple-step potentials, performed well to identify native structures. Therefore, the ideal protocol would consist of initial rigid-body sampling, scoring using a long-range potential, flexible refinement to move the structures closer towards the native form, and finally scoring using a short-range potential.

Finally, our results show that the GRADSCOPT training protocol is robust and offers viable strategies to generate problem-adapted solutions, as they are needed in practical docking problems.

Due to the usage of pre-calculated feature vectors for each potential form, the tool kit enables the user to generate and evaluate various scoring functions rapidly to find a problem-adapted solution.

Moreover, it also allows the user to optimize the functional form of the scoring potential (e.g. with respect to inter-atomic distances), as we found that the optimal functional form depends very strongly on the scoring problem. Other methods to create scoring potentials, such as linear programming [42][43], possess the advantage of numerically generating a globally optimal set of parameters with respect to the defined target function. However, the handling of these approaches appears to be more complex and their applicability seems to be narrower since it requires a set of training structures that do not impose impossible conditions on its constraints [START_REF] Tobi | Designing coarse grained-and atom based-potentials for protein-protein docking[END_REF]. Especially for the training on near-native solutions with diverse structures and interface compositions, or larger benchmarks with diverse types of protein-interactions, finding a set of restricting constraints can become difficult or impossible without excluding some important decoys.

The GRADSCOPT tool kit, a detailed manual with example executive scripts, the ATTRACT rigidbody protein-protein docking benchmark and parameter files for the created potentials can be found at http://www.t38.ph.tum.de/-> Downloads. as described in the methods section by the difference between the near-native and incorrect contributions (red and green).

Figure 1 .

 1 Figure 1. Work-flow to generate a knowledge-based scoring function by the GRADSCOP tool kit.Various scoring potentials can be generated and evaluated in parallel or sequentially by changing the decoy set, the protein representation, the potential form or the training method..

Figure 2 .

 2 Figure 2. (A) Comparison of the native complex pdb1KTZ (smaller ligand partner protein as green cartoon; receptor protein as beige vdW surface) and examples of acceptable (*)(red cartoon), medium (**) (purple) and high quality (***) (blue) docking solutions (using the CAPRI criteria, see Methods). (B) Clashes (orange) with the unbound chymotrypsin (pdb1ACB, beige vdW-surface) after superposition of the unbound chymotrypsin inhibtor structure (blue ) onto the bound form (light green ).

Figure 3 .

 3 Figure 3. Fraction of test cases (y-axis) for which the native bound complex structure was identified among the decoy complexes in the top 1 (red), top 10 (orange), top 100 (green), top 500 (light blue) and top 1000 (dark blue) using the scoring function indicated at the x-axis. The performances were sorted after the best performing scoring function in the top 10 from left to right.

Figure 4 .

 4 Figure 4. Fraction of test complexes for which at least one near-native (*,**,*** CAPRI-stars) structure was found in the test set in the top 1 (red), top 10 (yellow), top 100 (light blue), top 1000 (dark blue). The performances were sorted after the best performing scoring function in the top 10 from left to right.

Figure 5 .

 5 Figure 5. Average fraction of near-native structures (*,**,*** CAPRI-stars) in the decoy set whichwere identified in the best 0.1 (red), 1 (yellow), 2 (light blue) and 5 (dark blue) % of all decoys in the test set using the scoring functions indicated at the x-axis. The scoring performances were sorted after the fraction of near-native solutions in the top 5 % from left to right.

Figure 6 .

 6 Figure 6. 20 most negative (A) and 20 most positive (B) scoring contributions of the parameter (xaxis) for the MC_gaa_10 potential shown for the average native (blue), near-native (green) and incorrect (red) solution. The discriminating scoring contribution is defined as described in the methods section by the difference between the near-native and incorrect contributions (red and green).

Figure 1 .Figure 2 .

 12 Figure 1. Work-flow to generate a knowledge-based scoring function by the GRADSCOP tool kit. Various scoring potentials can be generated and evaluated in parallel or sequentially by changing the decoy set, the protein representation, the potential form or the training method.

Figure 3 .

 3 Figure 3. Fraction of test cases (y-axis) for which the native bound complex structure was identified among the decoy complexes in the top 1 (red), top 10 (orange), top 100 (green), top 500 (light blue) and top 1000 (dark blue) using the scoring function indicated at the x-axis. The performances were sorted after the best performing scoring function in the top 10 from left to right.

Figure 4 .

 4 Figure 4. Fraction of test complexes for which at least one near-native (*,**,*** CAPRI-stars) structure was found in the test set under the top 1 (red), top 10 (yellow), top 100 light blue, top 1000 (dark blue). The performances were sorted after the best performing scoring function in the top 10 from left to right.

Figure 5 .Figure 6 .

 56 Figure 5. Average fraction of near-native structures (*,**,*** CAPRI-stars) in the decoy set which were identified in the best 0.1 (red), 1 (yellow), 2 (light blue) and 5 (dark blue) % of all decoys in the test set using the scoring functions indicated at the x-axis. The scoring performances were sorted after the fraction of near-native solutions in the top 5 % from left to right.

Acknowledgements.

We acknowledge Christina Keller who aided our written style and grammar in our manuscript for revision.

Author Contributions

Conceived and designed the experiments: AS CS ICdB SJdV MZ. Performed the experiments: AS.

Analyzed the data: AS CS ICdB SJdV MZ. Contributed reagents/materials/analysis tools: AS CS SJdV. Wrote the paper: AS SJdV MZ. S1A) and for test set of 29 complexes (Table S1B).

Supporting information

S2


Christina Em Schindler

Isaure Chauvot De Beauchêne

Sjoerd De Vries

Martin Zacharias
email: martin.zacharias@ph.tum.de

Protein-protein and peptide-protein docking and refinement using ATTRACT in CAPRI

Keywords: elastic network model, docking minimization, protein-protein complex formation

The ATTRACT coarse-grained docking approach in combination with various types of atomistic, flexible refinement methods has been applied to predict protein-protein and peptide-protein complexes in CAPRI rounds 28-36. For a large fraction of CAPRI targets (12 out of 18), at least one model of acceptable or better quality was generated corresponding to a success rate of 67%. In particular for several peptide-protein complexes excellent predictions were achieved. In several cases, a combination of template-based modeling and extensive molecular dynamics-based refinement yielded medium and even high quality solutions. In one particularly challenging case, the structure of an ubiquitylation enzyme bound to the nucleosome was correctly predicted as a set of acceptable quality solutions. Based on the experience with the CAPRI targets, new interface refinement approaches and methods for abinitio peptide-protein docking have been developed. Failures and possible improvements of the docking method with respect to scoring and protein flexibility will also be discussed.

Introduction

Protein-protein and peptide-protein interactions are abundant in the cell and are involved in virtually all important biological processes. However, so far only a small fraction of complex structures has been characterized experimentally. Especially, structure determination of transient interactions between proteins is experimentally challenging and probably in many cases impossible to obtain. Since atomic structural knowledge is vital for understanding the biological roles of these interactions, protein-protein docking methods that predict the 3D structure of complexes and efficient refinement approaches have become increasingly important in structural biology. The Critical Assessment of PRediction of Interactions (CAPRI) experiment [START_REF] Lensink | Docking and scoring protein complexes: CAPRI 3rd Edition[END_REF][START_REF] Lensink | Docking and scoring protein interactions: CAPRI 2009[END_REF][START_REF] Lensink | Docking, scoring, and affinity prediction in CAPRI[END_REF] has provided a framework for blind testing and comparative assessment of protein-protein docking and refinement methods. Apart from evaluating the current state of the field, CAPRI's diverse, challenging targets have stimulated the development of new and more sophisticated protocols and pushed the limits of what is achievable in protein-protein docking. Our protein-protein docking approach ATTRACT [START_REF] Zacharias | Protein-protein docking with a reduced protein model accounting for side-chain flexibility[END_REF][START_REF] May | Accounting for global protein deformability during protein-protein and protein-ligand docking[END_REF][START_REF] May | Protein-protein docking in CAPRI using ATTRACT to account for global and local flexibility[END_REF][START_REF] Zacharias | ATTRACT: protein-protein docking in CAPRI using a reduced protein model[END_REF][START_REF] May | Energy minimization in low-frequency normal modes to efficiently allow for global flexibility during systematic proteinprotein docking[END_REF][START_REF] Fiorucci | Binding site prediction and improved scoring during flexible protein-protein docking with ATTRACT[END_REF][START_REF] Sj De | ATTRACT-EM: a new method for the computational assembly of large molecular machines using cryo-EM maps[END_REF][START_REF] Christina | iATTRACT: Simultaneous global and local interface optimization for protein-protein docking refinement[END_REF][START_REF] Sjoerd | A web interface for easy flexible protein-protein docking with ATTRACT[END_REF][START_REF] Sjoerd | Cryo-EM data is superior to contact or interface information in integrative modeling[END_REF] can predict proteinprotein and protein-nucleic acid interactions and has already been used successfullly in various rounds of CAPRI [START_REF] Zacharias | ATTRACT: protein-protein docking in CAPRI using a reduced protein model[END_REF][START_REF] May | Protein-protein docking in CAPRI using ATTRACT to account for global and local flexibility[END_REF][START_REF] Fiorucci | Binding site prediction and improved scoring during flexible protein-protein docking with ATTRACT[END_REF][START_REF] De | Flexible docking and refinement with a coarse-grained protein model using ATTRACT[END_REF]]. ATTRACT's main characteristics are its coarse-grained (CG) force field, the ability to incorporate conformational flexibility already during the initial large-scale search [START_REF] May | Accounting for global protein deformability during protein-protein and protein-ligand docking[END_REF][START_REF] Bastard | Accounting for loop flexibility during protein-protein docking[END_REF][START_REF] May | Energy minimization in low-frequency normal modes to efficiently allow for global flexibility during systematic proteinprotein docking[END_REF] and the possibility to dock any number of (protein) partners [START_REF] Sj De | ATTRACT-EM: a new method for the computational assembly of large molecular machines using cryo-EM maps[END_REF]. The CG model was derived from statistical analysis of protein-protein and protein-nucleic acid interfaces and is intermediate between a residue/base-level and full atomistic description. It represents each amino acid by up to four pseudoatoms (two for the backbone and one or two for the side chains) [START_REF] Fiorucci | Binding site prediction and improved scoring during flexible protein-protein docking with ATTRACT[END_REF]. A systematic docking search consists of several potential energy minimizations starting from hundreds of thousands of initial configurations. To speed up the docking calculations, the potential energy can be precalculated on a grid [deVries2016attractgrid ]. Global flexibility (e.g., domain-domain motion) can be included explicitly during docking by energy minimization along the directions of precalculated soft normal modes [START_REF] May | Energy minimization in low-frequency normal modes to efficiently allow for global flexibility during systematic proteinprotein docking[END_REF]. Side chain and loop conformational changes can be accomodated by a multi-copy [START_REF] Bastard | Accounting for loop flexibility during protein-protein docking[END_REF] or an ensemble docking approach. The flexible interface refinement method iATTRACT [START_REF] Christina | iATTRACT: Simultaneous global and local interface optimization for protein-protein docking refinement[END_REF] can be used to further optimize the rigid body docking solutions. We have participated in CAPRI rounds 28-36 and in the following report on our predictions for targets 59-107 and related new methodological developments. This also includes our efforts to design a fully blind peptide-protein protocol and molecular dynamics-based refinement approaches.

Methods

ATTRACT rigid body docking

The protein and peptide structures were converted to the ATTRACT atom type representation [START_REF] Zacharias | Protein-protein docking with a reduced protein model accounting for side-chain flexibility[END_REF] with the ATTRACT tool reduce. The empirical coarse grained force field in ATTRACT represents the amino acid side chains by one or two pseudoatoms and takes all polar backbone atoms without hydrogens into account. Interactions are based on LJ type potentials and can be either attractive or repulsive (using a saddle point instead of an energy minimum in the LJ-potential) [START_REF] Fiorucci | Binding site prediction and improved scoring during flexible protein-protein docking with ATTRACT[END_REF]. In addition, interactions between charged residues are calculated via a Coulomb term with a distance-dependent dielectric constant (ε = 15r). Starting points were generated by choosing random positions and orientations for the association partners with an appropriate center-of-mass distance to prevent steric overlap in the initial configuration. The starting structures are subjected to rigid body optimizations in a potential energy minimization of 1,000 minimization steps with the ATTRACT metric minimizer [START_REF] May | Accounting for global protein deformability during protein-protein and protein-ligand docking[END_REF][START_REF] May | Energy minimization in low-frequency normal modes to efficiently allow for global flexibility during systematic proteinprotein docking[END_REF]. Energy calculation was accelerated using a precalculated grid [deVries2016attractgrid ]. In cases where distance restraints were employed, the optimization of the six rigid body degrees of freedom was preceded by a minimization in which the center of mass positions were fixed and the docking partners could orient towards each other. During this stage, only the restraint potentials were applied ("ghost" mode). If multiple conformations were available for the association partners, these conformations were docked separately (ensemble docking). Finally, the docking models were ranked by their ATTRACT energy evaluated within a squared cutoff of 50 Å2 and highly similar models were removed with the deredundant tool.

iATTRACT flexible interface refinement

The protein and the peptide structures were converted into the OPLS atom type description with the ATTRACT tool aareduce. Missing hydrogens were built with PDB2PQR [START_REF] Dolinsky | PDB2PQR: an automated pipeline for the setup of Poisson-Boltzmann electrostatics calculations[END_REF][START_REF] Dolinsky | PDB2PQR: expanding and upgrading automated preparation of biomolecular structures for molecular simulations[END_REF] and protonation states were determined by PropKa [START_REF] Li | Very fast empirical prediction and rationalization of protein pKa values[END_REF]. Peptide termini were charged (unless the peptides were part of a larger protein), protein termini left uncharged. The atomistic refinement uses a physical force field based on the OPLS parameters to calculate non-bonded and electrostatic interactions between the protein partners. Contacts from the input structure are treated as flexible during a simultaneous potential energy minimization in rigid body degrees of freedom and interface flexibility [START_REF] Christina | iATTRACT: Simultaneous global and local interface optimization for protein-protein docking refinement[END_REF]. A structure-based force field is determined on-the-fly to evaluate intra-protein interactions for the flexible interface residues. Depending on the size of the target, a few hundreds to thousand models from rigid body docking with ATTRACT were selected for iATTRACT refinement. The refinement parameters were chosen as specified in [START_REF] Christina | iATTRACT: Simultaneous global and local interface optimization for protein-protein docking refinement[END_REF]. Structures were in general not rescored after iATTRACT refinement.

Molecular dynamics refinement

Atomistic refinement simulations were performed sing the AMBER14 molecular dynamics package [START_REF] Da Case | AMBER 14[END_REF] in combination with an implicit solvent description. The structures were converted to the AMBER atom type description using the pdb4amber tool. A Generalized-Born (GB) implicit solvent model (igb=8) was used with the newest version of the AMBER force field ff14SB [START_REF] Da Case | AMBER 14[END_REF]. The structures were first minimized with the sander program (500 steepest descend steps) with a short cutoff and a small initial step size to relax possible atom overlap and deformations resulting from the structure-based force field used in iATTRACT refinement. Then two short molecular dynamics (MD) simulations were run with the pmemd.cuda program of the AMBER package for 5000 and 2500 time steps (2 fs) at temperatures T = 400K and T = 350K respectively. During the MD simulations, intra-molecular distances and inter-molecular distances between backbone atoms were restrained to prevent large deformations and dissociation of the binding partners. The intra-molecular distances were restrained with a harmonic potential to the distance found in the unbound structure with force constant of 2 kcal/(mol Å2 ). The inter-molecular distances were allowed to change by up to 10 Å with respect to the distance found in the initial structure. Larger deviations beyond 10 Å were penalized by a harmonic potential with force constant 0.25 kcal/(mol Å2 ). Finally, the structures were minimized for 5000 steps with a large cutoff using the pmemd.cuda program without restraints. The force field energy was evaluated for the complex and the individual protein partners by the sander program. The binding interaction energy score was calculated by subtracting the energy of the free protein partners from the energy of the complex. The final models were ranked by their binding interaction energy score without using any reweighting scheme of the energy contributions. Note, that the pmemd.cuda program employs graphical processing units (GPUs) and results in refinement simulations within seconds to minutes for averaged-sized complexes (100-200 aminos acids per partner protein). This translates to refinement of hundreds to thousands of complexes per day on a single GPU node. For refinement in explicit solvent with GROMACS version 4.6 (www.gromacs.org) [START_REF] Herman | GROMACS: a message-passing parallel molecular dynamics implementation[END_REF][START_REF] Pronk | GROMACS 4.5: a high-throughput and highly parallel open source molecular simulation toolkit[END_REF][START_REF] James | GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers[END_REF], structures were converted into the gro format with the tool pdb2gmx. Simulations were run with explicit solvent using the TIP3P water model and the AMBER99SB-ILDN force field at a temperature of 300 K. The positions of the backbone atoms were restrained with force constant 1000 kJ/mol/ Å2 in x, y and z direction.

Results and Discussion

In CAPRI rounds 28-36, we submitted predictions for targets 59-67, 95-97 and 102-107 (we did not participate in the CASP-CAPRI experiment in round 30 and in CAPRI round 32; round 36 was canceled before the submission deadline). Since the templates for the targets in round 32 (targets 98-101) were all of very low sequence identity we did not attempt to model these targets and to predict possible complex structures. Note, that indeed none of the predictors achieved any successful prediction for round 32. Table 1 shows a summary of the results. The targets vary strongly in terms of interaction type and docking difficulty. Especially, many targets of the most recent CAPRI rounds involved a high degree of flexibility and consequently, results for all predictor groups were rather poor. In the following, we discuss our predictions and the challenges we faced for some of the targets.

Round 28 (Targets 59-64)

Round 28 consisted of one protein-protein complex (target 59) and five peptide-protein complexes (targets 60-64). Target 59 corresponded to the Edc3 LSm domain, an ac- tivator of the mRNA decapping complex, with a motif from Rps28B. NMR ensembles of structures for the LSm domain were available both in the apo form and bound to another motif and we used all NMR models removing the highly flexible regions prior to docking [START_REF] Tritschler | A divergent Sm fold in EDC3 proteins mediates DCP1 binding and P-body targeting[END_REF][START_REF] Simon | The structural basis of Edc3-and Scd6-mediated activation of the Dcp1: Dcp2 mRNA decapping complex[END_REF]. For Rps28B, we used an ensemble of homology models created with MODELLER [] based on available structures [START_REF] Wu | Solution structure of ribosomal protein S28E from Methanobacterium thermoautotrophicum[END_REF][START_REF] James M Aramini | Solution NMR structure of the 30S ribosomal protein S28E from Pyrococcus horikoshii[END_REF]. We performed ATTRACT ab-initio rigid body docking and subsequent molecular dynamics refinement. For this target, we obtained two acceptable structures. The best submitted model deviated from the native structure at the interface by 3.8 Å and retrieved 18 % of the native contacts. To date, a more detailed analysis for this target is not possible, since the experimental structure has not yet been published. Target 60-64 were structures of importin-α binding to different peptides derived from nuclear localization signals (NLSs). NLSs contain one or two clusters of basic residues and are recognized by the import receptor importin-α [START_REF] Chang | Distinctive Conformation of Minor Site-Specific Nuclear Localization Signals Bound to Importin-α[END_REF]. Several structures of importin-α in complex with different peptides were available at the time of round 28 which showed two binding sites: a major and a minor site. The peptides in Round 28 were derived from nuclear localization signals and had been optimized towards binding to the minor site [START_REF] Chang | Distinctive Conformation of Minor Site-Specific Nuclear Localization Signals Bound to Importin-α[END_REF]. In the crystal structure, the major site was also occupied and so evaluation of these targets was carried out for both the major and the minor site. We generated models for the peptide based on the available crystal structures and used this in an ATTRACT rigid body docking with ambiguous distance restraints twoards the major and the minor sites. The best models were then refined by an energy minimization with AMBER. For all but one target, we achieved A Major site.

B Minor site.

Figure 1 Peptide-protein docking model for target 60 (docked peptide indicated as red sticks, importinα receptor shown as grey surface) superimposed on the native structure (PDB 3ZIN, bound peptide in black). Peptide binding to importin-α was modeled using a combination of homology modeling and molecular dynamics refinement. For this target, and related targets 62-64, several three-star quality models were submitted high quality predictions (Figure 1). However, we failed to predict the α-helical turn in the peptide when binding to the minor site. Here, possibly more extensive peptide structure modeling [START_REF] Thévenet | PEP-FOLD: an updated de novo structure prediction server for both linear and disulfide bonded cyclic peptides[END_REF] prior to complex prediction could have improved the results.

Round 29 (Targets 65-67)

In Round 29, three distinct peptide-protein complexes were proposed as targets. Target 65 and 66 were complexes of proteins with a SSB C-terminal peptide. For both targets, we analyzed available structures of other proteins in complex with a SSB Cterminal peptide (PDB 3UF7, 3Q8D, 3C94) [START_REF] Ryzhikov | Mechanism of RecO recruitment to DNA by singlestranded DNA binding protein[END_REF][START_REF] Lu | Structural basis of Escherichia coli single-stranded DNA-binding protein stimulation of exonuclease I[END_REF] and identified a conserved binding mode with an hydrophobic anchor at the C-terminus and a solvent exposed part of the peptide, although in general the peptide appeared relatively flexible. We extracted the peptide conformations from the available complexes, recombined parts of the peptides from different crystal structures to generate additional conformations and used this peptide ensemble in a fully blind coarse-grained search of the entire protein surface combined with two stages of flexible refinement [START_REF] Christina | Fully blind peptide-protein docking with pepATTRACT[END_REF]. Before refinement, the rigid body models were filtered to detect models with a buried phenylalanine at the C-terminus.

Target 65 was a complex with RNase Hi [START_REF] Petzold | Interaction with single-stranded DNA-binding protein stimulates Escherichia coli Ribonuclease HI enzymatic activity[END_REF]. For this target we did not generate any near-native prediction, since we did not sufficiently model the protein's flexibility at the peptide binding site (we only used a single crystal structure for the protein during docking [START_REF] Ishikawa | Stabilization of Escherichia coli ribonuclease HI by cavityfilling mutations within a hydrophobic core[END_REF]) and failed to detect the correct binding pocket. Target 66 corresponded to a PriA helicase in complex with a SSB C-terminal peptide [START_REF] Bhattacharyya | Structural mechanisms of PriA-mediated DNA replication restart[END_REF]. We obtained one model of acceptable quality (top 1) and one model (top 2) which had a very similar orientation of the peptide but an incorrect peptide conformation (α-helical turn) (Figure 2). In our topranked model, the whole peptide was in contact with the protein surface, whereas in the A top 1

B top 2

Figure 2 Predicted structural models (top1 and top2) for target 66 (docked peptide in red, receptor as grey surface) superimposed on the native structure (PDB 4NL8, bound peptide in black).

crystal structure the N-terminal part of the peptide is solvent-exposed. The tendency to maximize the interface between a flexible peptide and a binding region can be attributed to the use of an implicit solvent model during final refinement (see Methods). In future cases, we will consider explicit solvation during refinement of peptide-protein complexes which may help to improve the accuracy of the predictions.

Target 67 corresponded to a WW domain in complex with an extended peptide containing proline (P) residues. Several structures of PPXY motifs bound to WW domains were available in the PDB [START_REF] Bobby | Structure and dynamics of human Nedd4-1 WW3 in complex with the αENaC PY motif[END_REF] and we used these to build initial peptide-protein complexes. Since only a limited number of complexes were generated it was possible in this case to use restrained molecular dynamics simulations in explicit solvent with GRO-MACS for refinement (1 ns simulations at room temperature and normal pressure, see Methods). Finally, the models were energy-minimized with sander of the AMBER package. All our 10 submitted models were at least of acceptable quality and the best model had an IRMSD of 0.8 Å and retrieved 88 % of the native contacts (medium quality).

Round 31 (Targets 95-97)

Round 31 comprised three very challenging protein-protein complexes as targets. For target 95, binding of an ubiquitylation enzyme PRC1 (Bmi1/Ring1b ubiquitin ligase) to the nucleosome had to be predicted [START_REF] Robert K Mcginty | Crystal structure of the PRC1 ubiquitylation module bound to the nucleosome[END_REF]. PRC1 ubiquitylates the histone H2A tail at residue LYS 119 [START_REF] De | Polycomb group proteins Ring1A/B link ubiquitylation of histone H2A to heritable gene silencing and X inactivation[END_REF][START_REF] Wang | Role of histone H2A ubiquitination in Polycomb silencing[END_REF]. Furthermore, several residues on PRC1 (ASP 56, LYS 92, LYS 93, LYS 97 and ARG 98 on Ring1b and LYS 62 and ARG 64 on Bmi1) had been identified as important for binding by mutational experiments [START_REF] Matthew | Recognition of UbcH5c and the nucleosome by the Bmi1/Ring1b ubiquitin ligase complex[END_REF]. There was also experimental evidence for PRC1 binding to DNA (although only for isolated DNA, not for the nucleosome) [START_REF] Matthew | Recognition of UbcH5c and the nucleosome by the Bmi1/Ring1b ubiquitin ligase complex[END_REF] and data pointing to an important role for the acidic patch in PRC1 function [START_REF] Leung | Nucleosome acidic patch promotes RNF168-and RING1B/BMI1dependent H2AX and H2A ubiquitination and DNA damage signaling[END_REF]. We used the unbound protein structures for PRC1 [START_REF] Matthew | Recognition of UbcH5c and the nucleosome by the Bmi1/Ring1b ubiquitin ligase complex[END_REF] and the nucleosome [START_REF] Vasudevan | Crystal structures of nucleosome core particles containing the 601strong positioning sequence[END_REF] and performed a large scale rigid body docking search with a 10 Å upper harmonic distance restraint between the C α atoms of residue LYS 118 on the H2A tail A Top view.

B Side view. (LYS 119 was not resolved in the structure) and residue CYS 85 in the active site of the ubiquitin ligase. Subsequently, we refined the best clusters with iATTRACT and short molecular dynamics refinement with AMBER. We found two types of solutions among the top-ranked models: one with PRCC1 contacting the acidic patch and one with PRCC1 binding to nucleosomal DNA similar to an earlier model [START_REF] Matthew | Recognition of UbcH5c and the nucleosome by the Bmi1/Ring1b ubiquitin ligase complex[END_REF]. We hence submitted 5 models for each type of solution and achieved 5 acceptable model with the best model being very close to medium quality (LRMSD 5.03 Å, fnat 0.5 Å, Figure 3). A posterior comparison of our models with the experimental crystal structure of the complex indicated a conformational adjustment of a histone tail upon PRCC1 binding not accounted for in our docking efforts. Hence, the sterical clash of the docked PRCC1 with the histone tail in the unbound nucleosome structure limited the accuracy of our models to the level of acceptable CAPRI solutions. Similar to other CAPRI cases it appears to be necessary in the future to carefully check for highly flexible segments (termini or loops) in partner structures prior to docking and possibly eliminate these elements during docking followed by rebuilding after docking.

Targets 96 and 97 were complexes of designed α-repeat proteins binding to GFP (PDB 4XL5 and 4XVP) [START_REF] Chevrel | Specific GFP-binding artificial proteins (αRep): a new tool for in vitro to live cell applications[END_REF]. We used our standard ab-initio docking protocol in combination with additional scoring with the Rosetta force field which has been used in the past for protein design [START_REF] Wang | Protein-protein docking with backbone flexibility[END_REF][START_REF] Kaufmann | Practically useful: what the Rosetta protein modeling suite can do for you[END_REF]. We achieved acceptable quality for target 96. However, we failed to accurately predict the smaller target 97. In this case docking sampling gave solutions where the ligand protein (GFP) in our models was rotated by 180°with respect to the native structure. This failure can be attributed to deficiencies of the scoring of the docking models selected for further refinement.

Round 33 and 34 (Targets 102-105)

In Round 33, prediction of two protein-protein complexes was offered (targets 102 and 103). Results for target 102 were not available at the time of writing. For target 103, we did not achieve any successful predictions although for one submitted model an overlap with the correct binding interface was achieved for model 10 (fnat=0.27, Table 1). In Round 34, the challenge was to model the placement of interfacial waters. The targets were pyocin DNase domains in complex with immunity proteins which were structurally similar to previously solved structures in colicin. We built an initial model of the complex by homology modeling using a colicin complex [START_REF] Aleksandra | Structure of the ultra-high-affinity colicin E2 DNase-Im2 complex[END_REF] and refined this with backbone restrained molecular dynamics simulations in explicit solvent. We achieved a high quality model for target 104 and medium quality models for target 105. For the prediction of the interfacial water placements we followed our molecular dynamics based protocol described previously [START_REF] De | Flexible docking and refinement with a coarse-grained protein model using ATTRACT[END_REF][START_REF] Marc F Lensink | Blind prediction of interfacial water positions in CAPRI[END_REF]. Briefly, in this protocol water molecules are placed using the AMBER leap module allowing for partial overlap of waters with solute atoms at the interface (resulting in a slight over-hydration of the interface). The hydration structure is then allowed to relax during short MD simulation including positional restraints on the protein backbone. Finally, all waters outside the interface are removed followed by energy minimization to optimize the position and orientation of each water molecule at the interface (see details in reference [START_REF] Marc F Lensink | Blind prediction of interfacial water positions in CAPRI[END_REF]). Interestingly, the prediction of interfacial waters was more accurate for target 105 (medium quality model), probably because there were less water molecules to predict for this target. Interestingly, in the overall CAPRI evaluation it was found that high accuracy water predictions were not necessarily generated only for high quality models of the protein-protein complex (Lensinck, unpublished data). This is in contrast to previous results for interfacial water predictions in CAPRI [START_REF] Marc F Lensink | Blind prediction of interfacial water positions in CAPRI[END_REF] and indicates that probably only certain key groups at the interface have to be in near-native position to allow correct placement of waters.

Round 35 (Targets 102-103 and 106-107)

In these two rounds, the targets were a complex of haemopexin with the haemopexin utilization protein (huxA). The complex was solved in the apo state and with haem bound. It was previously speculated that the C-terminal domain of huxA interacts with haemopexin [START_REF] Fournier | Haem release from haemopexin by HxuA allows Haemophilus influenzae to escape host nutritional immunity[END_REF]. However, at the time of submission for Round 33, structural data was only available for the N-terminal secretion domain [START_REF] Baelen | Structure of the secretion domain of HxuA from Haemophilus influenzae[END_REF]. Hence, the target was proposed again as a target in Round 35 offering an unbound structure of full-length huxA. Still, we were not able to generate any near-native predictions, since binding of huxA to haemopexin involves a conformational change in a loop on HuxA which is partly disordered in the unbound form. Unfortunately, this loop yields steric clashes with haemopexin when superimposing the unbound huxA structure on the complex. A strategy involving detailed flexibility analysis prior to docking, removal of flexible loops and a-posteriori loop rebuilding would have certainly yielded improved results for this target.

Conclusion and Outlook

In CAPRI Rounds 28-36, we submitted predictions for 18 targets (16 distinct targets) and achieved at least acceptable predictions in 12 cases which can be considered as highly successful given the difficulty of many of the targets. The coarse-grained ATTRACT docking approach in combination with different refinement schemes proved to be versatile in dealing with a variety of different targets that ranged from peptide-protein interactions to docking of large proteins to the nucleosome. The CAPRI challenge has also triggered the development of several extensions of the original ATTRACT approach in the area of peptide-protein docking (pepATTRACT) [START_REF] Christina | Fully blind peptide-protein docking with pepATTRACT[END_REF] and of refinement (iATTRACT) [START_REF] Christina | iATTRACT: Simultaneous global and local interface optimization for protein-protein docking refinement[END_REF] within the last rounds. For several CAPRI targets, we found that homology modeling in combination with restrained molecular dynamics refinement yielded high quality predictions. In some cases, significant conformational change or inaccuracies in the homology modeling have contributed to the failure of the docking search. More detailed flexibility analysis prior to docking is needed to select appropriate conformational ensembles and to identify flexible loops. It might be beneficial in many cases to eliminate such highly flexible loops prior to docking if the chances to correctly model the bound forms are small and to tackle the generation of a bound loop structure after docking. Such approaches may also be useful for improving protein-protein complexes predicted based on sequence similarity to a known template complex. Furthermore, improvements in scoring towards distinguishing structures of medium or higher quality from just acceptable solutions is highly desirable and at the focus of our future research.

Figure 3

 3 Figure 3 Best model (top 3) for target 95 superimposed on the native structure (PDB 4R8P, black). PRC1 binds to the nucleosome at the histone acidic patch.

Table 1

 1 Results for ATTRACT predictions in CAPRI rounds28-36.

	Target	Type	Best model	fnat	IRMSD [ Å]		Classification
	59	protein-protein	5	0.18	3.8		acceptable (2 )
	60	peptide-protein	5/6	1.0/0.94	0.44	high (2	,4 )
	61	peptide-protein	1	0.76	0.5		medium (5 )
	62	peptide-protein	1	0.92	0.39	high (2	, 3 )
	63	peptide-protein	1	0.81	0.49		high (2	, 3 )
	64	peptide-protein	2	0.87	0.42	high (3	, 2 )
	65	peptide-protein	2	0.27	3.7		incorrect
	66	peptide-protein	1/1	0.5/0.75	1.3/2.1		acceptable (1 )
	67	peptide-protein	5	0.88	0.8	medium (2 , 8 )
	68-94	protein-protein	-	-	-		no submission
	95	protein-protein	3	0.5	3.3		acceptable (5 )
	96	protein-protein	3	0.15	2.96		acceptable (1 )
	97	protein-protein	8	0.05	12.2		incorrect
	98-101	protein-protein	-	-	-		no submission
	102	protein-protein	-	-	-	not yet available
	103	protein-protein	10	0.27	12.2		incorrect
	104	interfacial water	4	0.68	0.92	high (1	,9 ), water (5++)
	105	interfacial water	3	0.66	1.2	medium (10 ), water (7+++)
	106	protein-protein	-	-	-	not yet available
	107	protein-protein	7	0.05	17.42		incorrect
	108-109		-	-	-		canceled

Acknowledgements

The authors thank the organizers of the CAPRI challenge for this opportunity and the assessors for the hard work to evaluate the predictions. They also thank all structural biologists who made their structures available as targets for the CAPRI experiment. CS acknowledges funding by the Center for Integrated Protein Science Munich.