Martin Engelhard

Georg B Üldt

Giuseppe Zaccai

Gr Égory

Durand Phd

Valentin Gordeliy

Sergei Grudinin

The phenotype of every known living organism is determined mainly by the complicated interactions between the proteins produced in this organism. Understanding the orchestration of the organismal responses to the external or internal stimuli is based on the understanding of the interactions of individual proteins and their complexes structures. The prediction of a complex of two or more proteins is the problem of the protein-protein docking field. Docking algorithms usually have two major steps: exhaustive 6D rigid-body search followed by the scoring. In this work we made contribution to both of these steps. We developed a novel algorithm for 6D exhaustive search, HermiteFit. It is based on Hermite decomposition of 3D functions into the Hermite basis. We implemented this algorithm in the program for fitting low-resolution electron density maps. We showed that it outperforms existing algorithms in terms of time-per-point while maintaining the same output model accuracy. We also developed a novel approach to computation of a scoring function, which is based on simple logical arguments and avoids an ambiguous computation of the reference state. We compared it to the existing scoring functions on the widely used protein-protein docking benchmarks. Finally, we developed an approach to include water-protein interactions into the scoring functions and validated our method during the Critical Assessement of Protein Interactions round 47.

Dedication

To my parents.

List of Figures

List of Tables

Role of proteins and their interactions in the cell

The main dogma of biology states: information in a cell flows from DNA to RNA and then to proteins. Although the actual working of the cell is much more comlpex to be contained in any dogma, it captures a large picture of the state-of-art of our understanding of the working of a cell. Proteins in this picture have a crucial place: they are the cogs of the cell machinery. These are the protein expressed in a cell and their interactions that mainly define phenotype of that cell. Therefore gainig new data on the protein interaction in a cell profoundly enhances our capabilities to understand and change living organisms for the needs of the society.

Classification of protein-protein interactions

Commonly used classification of the protein-protein complexes uses the following main criteria: composition of a complex, its lifetime and stability [START_REF] Nooren | Diversity of protein-protein interactions[END_REF].

Protein composition and interface

The complex can be comprised of two identical or highly homologous proteins. An example of such a complex is the D-alanyl carrier protein (pdb code 4BPG), which consists of two identical polypeptide chains (Fig. 1.1). The proteins of this class are called homo-oligomers. On the other hand, if the complex consists of two distinct proteins it is named hetero-oligomer. An example of hetero-dimer (Glycosidase CelD bound to artificial affitin E12 protein, [START_REF] Correa | Potent and Specific Inhibition of Glycosidases by Small Artificial Binding Proteins (Affitins)[END_REF], pdb code 4CJ0) is shown on Fig 1 .2. Homo-oligomers are further divided dependent on the interface of oligomerization. If the two chains bind using the identical binding interace, their complex is called isologous. We can see that D-alanyl carrier protein is an example of isologous homodimer (Fig. 1.2). On the contrary, the homo-oligomeric assemblies where proteins bind through distinct interfaces are called heterologous.

Stability of protomers

The proteins that form a complex can be either stable or not on their own in vivo.

In the first case, when a complex is formed out of stable proteins it is called non- Figure 1.2: Example of a heterodimeric complex, Glycosidase CelD bound to artificial affitin E12 protein, [START_REF] Correa | Potent and Specific Inhibition of Glycosidases by Small Artificial Binding Proteins (Affitins)[END_REF], pdb code 4CJ0

obligate. In many cases these proteins are localized in different compartements of a cell. Usually, interactions between receptor-ligand, enzyme-inhibitor, etc belong to non-obligate interactions. In the other case when a complex of one or both chains do not exist as folded proteins on their own, it is called obligate. An example of such a complex is the Arc repressor dimer, which consists of two peptide chains. Upon dimerization (catalized by DNA [START_REF] Marcovitz | Arc-repressor dimerization on DNA: folding rate enhancement by colocalization[END_REF]) these chains obtain stable secondary structure, that was absent in the monomer.

Complex lifetime

If two proteins form a very stable complex they said to interact permanently. In many cases obligate complexes, such as Arc repressor dimer, are permanent. On the other hand, if a complex is dissociated and formed continuously in vivo, the interaction between its constituting proteins is called transient. Many of non-obligate interactions are also transient. However, there are strong transient interactions that require the presense of some molecule for the complex dissociation. An example of a strong transient interaction is the formation of trimeric G protein, for which guanine triphosphate is the dissociation trigger.

Experimental methods to probe protein-protein interactions

Giving importance and variety of protein-protein interactions in a cell, numerous methods were developed to discover the fact that two or more proteins form a complex by measuring kinetics and free energy of its formation or even by solving the atomic structure of the complex. The openness of the proteomics community also resulted in the construction of many databases of protein-protein interactions.

All experimental methods of identification and characterization of protein-protein interactions can be classified using two parameters: scale(high-throuput, individual) and system( in vivo, in vitro). Usually, a low-throuput approach gives much more information about interaction properties, whereas a high-throuput one only identifies the presense or absense of an interaction.

Yeast two-hybrid method (Y2H)

The Y2H is probably the oldest and most widespread method to identify interaction of proteins in vivo, which was scaled up to proteome level. In the pioneering work by Fields S. and Song O. [START_REF] Fields | A novel genetic system to detect protein-protein interactions[END_REF] it was shown that if the domain of the transcription activator that directs binding to a promoter (BD) is separated from a domain that activates transcription from this promoter (AD), the transcription is deactivated.

Based on this principle, two proteins of interest are fused to BD and AD, respectively, and the transcription of the reporter gene signals if they interact (Fig 1 .3). The are at least two scaling schemes: matrix and library approaches [START_REF] Uetz | A comprehensive analysis of protein-protein interactions in Saccharomyces cerevisiae[END_REF][START_REF] Ito | A comprehensive two-hybrid analysis to explore the yeast protein interactome[END_REF]. In the first one, clones expressing different proteins X i are taken and plated each in the rows of wells on a plate. The same set of clones is plated in colums of wells on a plate. Then, when the two yeast cells in each well mate and make a diploid that expresses both proteins X i and X j . Afterwards, the expression of a reporter gene in a cell i, j means that proteins X i and X j interact. In the library-based approach, clones containig protein X i are screened agains a library of clones with various proteins, which can be expressed from random cDNA or all open reading frames of a certain genome. In this case, the diploids expressing interacting proteins are selected agains specific growing media. The proteins that interact with X i are determined by the DNA sequencing.

UAS

Mass spectrometry and tandem affinity purification

In this method, a certain protein is fused with the DNA sequence that expresses the TAP tag (IgG binding domains of Staphylococcus protein A and calmodulin binding peptide separated by the TEV protease site [START_REF] Puig | The tandem affinity purification (TAP) method: a general procedure of protein complex purification[END_REF]). When the DNA construct is expressed in a host, it produces the protein of interest that forms complexes with other proteins of the host. During the purification, these protein complexes bind to the IgG matrix. Other proteins that did not form complexes with the TAPtagged one, hop through the matrix. Afterwards, using TEV protease, one cleaves of IgG binding domains and purified complexes are eluated from the matrix. The second purification step is the binding to calmodulin-coated beads. The simple representation of the two-step purification is shown on Fig. 1.4. Finally, the eluate is loaded to the SDS-PAGE gel and the resulting bands are cleaved by proteases.

After the purification, one uses mass-spectrometry to identify the fragments of the cleaved protein-protein complexes [START_REF] Di Tullio | Molecular recognition by mass spectrometry[END_REF]. Mass-spectrometry identifies particles based on their charge-to-mass ration. It allows to recognize fingerprints of short peptides and therefore identify the proteins using the solution of peptide fragments. A certain advantage of the TAP-MS method over Y2H is that it can detect not only dimeric protein-protein interaction, but also multimeric complexes. Step 1

Step 2

Gene co-expression

Recently the methods that allowed measuring expression of the genes on the scale of the whole cell were invented. Using this methodology, one can measure the similarity of expression profiles over different conditions of the two or more genes that code for interacting proteins. It turns out that they are significantly more similar than the gene expression profiles of random non-interacting proteins [START_REF] Jansen | Relating whole-genome expression data with protein-protein interactions[END_REF].

Synthetic lethality

The mutations in the genome of an organism affect its phenotype. In many cases this phenotype alternation is caused by the change in the protein-protein interaction network. By introducing random mutations or deletions in the genes of two proteins, one can monitor survival rate of these cells. The lethality of one such mutation can point to presence of interaction between the two target proteins [START_REF] Siew Loon Ooi | Global synthetic-lethality analysis and yeast functional profiling[END_REF].

Fluorescence resonance energy transfer

The fluorescence resonance energy transfer (FRET) occurs between two molecules: one (donor) in an excited state and the orther (aceptor) in the ground state. The energy is transfered through the dipole-dipole interaction between the molecules and does not involve emission [START_REF] Yan | Analysis of protein interactions using fluorescence technologies[END_REF]. The probability of such transfer to occur strongly depends on the distance between the donor and the aceptor. However, it is almost independent of the environmental conditions, which makes this method a great tool to study molecular interactions in vivo. More specifically, the two fluorophores are fused to the two proteins of interest (Fig. 1.5). One is then excited using a laser and if these molecules form a complex, it transfers energy to the second fluorophore that emits light at a certain wavelength. One of examples of FRET application is the investigation of membrane proteins dimerization in vivo, in particular of melatonin receptor types 1A and 1B [START_REF] Mohammed A Ayoub | Preferential formation of MT1/MT2 melatonin receptor heterodimers with distinct ligand interaction properties compared with MT2 homodimers[END_REF].

Isothermal titration calorimetry

This method allows to measure stoichiometry, dissociation constant, enthalpy and entrophy of the binding reaction between two proteins [START_REF] Velázquez | ITC in the post-genomic era. . . ? Priceless[END_REF]. The experimental setup consists of two thermally isolated chambers. One of the chambers is used as a reference and is filled with water. The other contains one of the interacting proteins. Its interaction partner is titrated in a known amount to the chamber. The thermometer measures temperature difference between these two chambers and controls heaters to equilibrate their temperature (while maintaining the temperature of the reference chamber constant). The amount of heat spent to make the two chamber isothermic is measured during the experiment. Figure 1.6 shows a schematic example of data obtained from an ITC experiment. The red line shows the peaks of energy transfer that correspond to titration events. Fitting these peaks with the interaction model (blue dashed line on Fig. 1.6) of a titration substance and a substance in the reservoir one can deduce the parameters of the model. In the simplest case of protein-ligand interaction one obtains enthalpy, dissociation constant and stoichiometry of the reaction.

Energy transfered

Molar ratio

Nuclear magnetic resonance spectroscopy (NMR)

This method is based on the absorbtion and emission of radio-frequency radiation by the nuclei of certain atoms. The emission and absorbtion spectra depend on the environment of nuclei and therefore the measurement allows to reconstruct the distance matrix between certain atoms in a protein in solution. This method also allows to measure the dynamics of proteins and their interaction. The classical approach is based on the Nuclear Overhauser Effect (NOE). One of the examples of the NMR results is the structure of the Twist protein, a transcription factor that plays a key role in the epithelial-mesenchymal transition and the bromodomaincontaining protein 4 [START_REF] Shi | Disrupting the Interaction of BRD4 with diacetylated twist suppresses tumorigenesis in basal-like breast cancer[END_REF] (see Fig. 1.7). Due to the complexity of spectra measured during the experiment, this approach is limited to protein complexes of the size up to 300 amino acids.

Cryo electron microscopy

The broad variety of methods go under the name of Cryo electron microscopy (Cryo-EM). The basic principle underlying this method is the same as in the conventional light microscopy. The difference lies in the radiation source. In Cryo-EM, the electron beam accelerated up to 300kV is used. This allows increasing the resolution that depends on the wavelength of the incident radiation, up to atomic one (wavelength of electrons accelerated by 300kV is about 0.02 Å). However, the increased resolution goes with the increased radiation dammage to the sample. Variety of methods are used to diminish ionizing effect of electrons on the biological samples. More precisely, measurement are conducted at low temperatures [START_REF] Dubochet | Cryo-electron microscopy of vitrified specimens[END_REF], averaging many identical units [START_REF] Kenneth | Electron diffraction of frozen, hydrated protein crystals[END_REF], single-particle microscopy. The last one is probably the most used one in structural biology.

The single-particle Cryo-EM is based on collecting information from many 2D projections of an object to reconstruct its 3D low-resolution model (electron density map, EDM). Afterwards, the supplementary information provided by NMR of X-Ray crystallography is used to construct an atomistic model of the object.

The 3D reconstruction of the EDM is usually based on the central projection theorem. It states that the Fourier images of 2D projections of a 3D object are the central slices of its 3D Fourier images. The relative positioning of two projections can be derived from the common lines in their Fourier images. These algorithms are implemented in the programs like IMAGIC [START_REF] Marin Van Heel | A new generation of the IMAGIC image processing system[END_REF], SPIDER [START_REF] Tanvir R Shaikh | SPIDER image processing for single-particle reconstruction of biological macromolecules from electron micrographs[END_REF], FREALIGN [START_REF] Grigorieff | FREALIGN: high-resolution refinement of single particle structures[END_REF] etc.

The resolution that can be obtained using this technique is highly dependent on the symmetry of the object measured. For example, highly symmetric icosahedral viruses envelopes were measured with up to 3 Å resolution [START_REF] Zhang | 3.3 Å cryo-EM structure of a nonenveloped virus reveals a priming mechanism for cell entry[END_REF]. The other examples of reconstructed protein complexes usually have resolution in the range of 7 Å-15 Å. This method is the main source of information on the large protein-protein assemblies today, like ribosomes [START_REF] Pawel A Penczek | Identifying conformational states of macromolecules by eigen-analysis of resampled cryo-EM images[END_REF] and chperonines [START_REF] Cong | 4.0-Å resolution cryo-EM structure of the mammalian chaperonin TRiC/CCT reveals its unique subunit arrangement[END_REF].

X-ray crystallography

The X-ray crystallography is the most time-and cost-demanding method, but it provides the most detailed atomistic information on the structure of protein-protein complexes. The starting point of this technique is a protein crystal. Obtaining crystal of a certain protein-protein complex is often the major hurdle and sometimes even impossible. Large protein-protein complexes and membrane proteins are especially hard to crystallize. Nontheless, the X-ray crystallography remains a major method to study structures of protein-protein complexes. The method is based on the diffraction of X-rays on the atoms, arranged in a lattice. The crystal of a protein is rotated with respect to the incident beam and the diffraction patterns are measured for each rotation. From these patterns one can reconstruct the absolute values of the Fourier image of a unit cell of the crystal. Afterwards, the phases of the Fourier image have to be reconstructed or measured. One of the most used way to obtain the phasing is the so-called molecular replacement. This means that approximate theoretical molecular structure of the unit cell is fitted into the given dataset. Afterwards, one minimizes discrepancy of the fitted structure with the diffraction pattern and obtains the final one. The quality of the final structure is judged by the R-factor. It shows to what extent the resulting structure explains the diffraction peaks.

Computational methods to probe protein-protein interactions

The experimental methods described in the previous section give rich information on the protein-protein interactions. However, in many cases the interactions being identified using these techniques are incomplete and contradictory owing to certain limitations and biases of the experimental conditions. To validate and cover the unknown spots in the protein-protein interactome maps, a plethora of computational techniques is used. These methods rely on different assumptions and use different information sources to decipher interaction details of proteins. Despite the rich variety one can approximatelly classify them in two major groups: top-down methods that use whole-organism or even evelutionary information and bottom-up, methods that employ the knowledge of single protein structures. They also differ by the amount of information they provide: ranging from a simple fact of interaction down to the details and precise conformation of the interaction interface.

Top-down approaches

This class of methods use the evolutionary and genomic data to predict if two proteins interact and identify domains that contain the interaction interface.

Gene neighbour and gene cluster methods

This pack of methods rely on the assumption that genes encoding for possibly interacting proteins often transcribed as a single operon in procariotes or co-regulated in eukariotes. A simplified example on how co-regulation maintains a certain stechiometry of a complex is shown on Fig. 1.8. For example, it was found that up to 75% of co-regulated genes in bacterial and achaeal genomes interact [START_REF] Dandekar | Conservation of gene order: a fingerprint of proteins that physically interact[END_REF]. The evolution tends to shuffle the order of genes in the distantly related organisms, however co-regulated gene clusters are found to be conserved. A prominent application of this method is the prediction of interaction of exsosome complex, that is capable of degrading viral RNA, and the RNAse P complex by comparing the order of genes in archaeal and eukariotic genomes [START_REF] Eugene V Koonin | Prediction of the archaeal exosome and its connections with the proteasome and the translation and transcription machineries by a comparative-genomic approach[END_REF].

co-regulated A B C A B C co-regulated A A A B A A B C A B C A B A A B A

B possible interaction gene regulation

Figure 1.8: Example of a gene co-regulated cluster and the probable complex with the stechiometry set by the co-regulation.

Phylogenetic profile methods

These methods are based on the hypothesis that the interacting proteins coevolve and therefore have orthologous proteins among sequenced organisms [START_REF] Barker | Predicting functional gene links from phylogenetic-statistical analyses of whole genomes[END_REF]. If the proteins related to the two proteins in question are present in majority of organisms therefore, they probably constitute a pathway or physically interact. The phylogenetic profiles are constructed for the proteins where their presense indicated by 1 or 0. Then, these profiles are clustered and the proteins belonging to one cluster are assumed functionally related or interacting.

Rosetta Stone method

This method relies on the observation that interacting proteins have homologs in other organisms that are fused into one protein. These types of proteins are called Rosetta Stone proteins. This fusion is the limiting case of the co-expression optimization of functionally related proteins. Using this method Marcotte et all [START_REF] Edward M Marcotte | Detecting protein function and protein-protein interactions from genome sequences[END_REF] identified about 7,000 pair of potentially interacting proteins in E.Coli and further analysis of the data revealed that around a half of these pairs are functionally related.

Co-evolution based methods

During the evolution, mutations in one of the proteins of a complex should be compensated by the mutations in its partner in order to maintain the function (see Fig. 1.9). The co-evolution based methods use the similarity measures between phylogenetic trees of two interacting protein families. Studies showed that some implementations of this method can predict up to 50% of real interactions with false positive rate as low as 6.4% [START_REF] Pazos | Similarity of phylogenetic trees as indicator of protein-protein interaction[END_REF].

Protein family A

Protein family B Figure 1.9: Example of a genes co-evolution where evolutionary changes in one protein compensate for the changes in its interaction partner.

Bottom-up approaches

Bottom-up approaches to the prediction of protein-prtoein interactions rely on the known structure of the proteins that constitute the complex. The structures can either be modelled by homology, obtained using evolutionary constraints or taken from such experiments as NMR, X-Ray crystallography or Cryo-EM. These methods are usually applied to infer the missing structural information about the proteinprotein complex rather than to predict the fact of interaction. Because bottom-up approaches use two or more protein structures that they dock into a complex, this class of techniques are usually called docking algorithms. The pioneering work in this field was done by Janin and Wodak [START_REF] Shoshana | Computer analysis of protein-protein interaction[END_REF]. Since then, this field has grown enormously, with its own quality assessment [START_REF] Méndez | Assessment of blind predictions of protein-protein interactions: current status of docking methods[END_REF][START_REF] Méndez | Assessment of CAPRI predictions in rounds 3-5 shows progress in docking procedures[END_REF][START_REF] Janin | Assessing predictions of protein-protein interaction: the CAPRI experiment[END_REF]. The two major parts of all the algorithms in this field are: sampling and scoring [START_REF] David | Recent progress and future directions in protein-protein docking[END_REF][START_REF] Huang | Search strategies and evaluation in protein-protein docking: principles, advances and challenges[END_REF]. As the name suggest, the sampling part generates putative conformations of a complex of proteins and scoring ranks them according to some criteria.

Search strategies

Suppose we have two molecules that are assumed to be rigid. One of the molecules (usually the one of a higher molecular weight) is called receptor and is fixed at the origin of the coordinate frame. The other, called ligand, is moved around the receptor. The global search algorithm explores all possible rotations and translations of the ligand movements. Complexity of this problem is O(N 9 ), where O(N 3 ) comes from rotational, O(N 3 ) from the translational degrees of freedom and O(N 3 ) is the complexity of the integration of the overlap integral that is computed for each rotation and translation of the ligand. A rough estimate gives ≈ 10 10 [55] operations. Several approaches allowed reducing this enormous complexity by at least an order of magitude. These are the fast Fourier transform correlation and clever heuristics in the direct search.

FFT-based rigid body docking

The idea to perform correlations in the Fourier space was first used by Katchalski-Katzir and colleagues [START_REF] Katchalski-Katzir | Molecular surface recognition: determination of geometric fit between proteins and their ligands by correlation techniques[END_REF]. In this section I give a short description of this approach with a comprehesive 2D example. Suppose that proteins receptor (R) and ligand (L) are represented as the numbers on a 3D grid:

R(l, m, n) =      1, (l, m, n) ∈ surface of R ρ, (l, m, n) ∈ inside of R 0, (l, m, n) ∈ outside of R L(l, m, n) =      1, (l, m, n) ∈ surface of L δ, (l, m, n) ∈ inside of L 0, (l, m, n) ∈ outside of L (1.1)
where ρ -1 and 0 < δ < 1, according to the original work by Katchalski-Katzir [START_REF] Katchalski-Katzir | Molecular surface recognition: determination of geometric fit between proteins and their ligands by correlation techniques[END_REF]. The shape complementarity score will therefore read:

C(p, q, r) = N l,m,n=1 R(l, m, n) × L(l + p, m + q, n + r), (1.2)
where the periodic boundary conditions are applied. Function C defined on the same grid gives complementarity score that depends on the shift of the ligand (p, q, r). This function can be computed using the FFT algorithm as follows:

C(p, q, r) = FFT -1 FFT[R] × FFT[L] ,
where FFT -1 stands for the inverse fast Fourier transform and overline for the complex conjugate. This algorithm allows to rapidly sample translational degrees of freedom. The rotations of the protein are separated from the translations in the outer loop of the algorithm. Fig. 1.10 shows a comprehensible example of 2D calculations using this approach.

An example of the docking procedure is shown on Fig. 1.10. The first row shows the picture of receptor and its Fourier image. The first column shows different ligand poses. The border outlined with the dark-blue and the inner space with light-blue. The numbers placed on the grid are described by Eq. 1.1. The second column shows the Fourier images of receptor and individual ligand poses. Third column shows the phases of Fourier image of the convolution FFT[C(p, q, r)]

= FFT[R] × FFT[L].
The phases of the inverse Fourier transform of the data from the third column are shown in the fourth one, it was computed according to Eq. 1.3.2.1.1. Finally, in the last column the pose with the best score is shown, here the ligand is shown with the light-blue and the receptor is depicted with the dark-blue color. The images of Fourier transforms were built using the projection of complex values to the HSV colorspace [START_REF] Petrisor | Visualizing complex-valued functions with Matplotlib and Mayavi[END_REF] with the amplitude fixed to V = 100.0. The algorithms that use FFT-based search strategies outnumber the algorithms relying on the other types of search strategies, to name a few: FTDock [START_REF] Henry A Gabb | Modelling protein docking using shape complementarity, electrostatics and biochemical information[END_REF], GRAMM [START_REF] Ilya | Evaluation of GRAMM low-resolution docking methodology on the hemagglutinin-antibody complex[END_REF], ZDOCK [START_REF] Brian G Pierce | Accelerating protein docking in ZDOCK using an advanced 3D convolution library[END_REF], PIPER [START_REF] Kozakov | PIPER: An FFT-based protein docking program with pairwise potentials[END_REF], etc. This technique is also used to accelerate the search not only in the space of translations but also in the rotational space [START_REF] José | FRODOCK: a new approach for fast rotational protein-protein docking[END_REF][START_REF] Kemp | Protein docking using spherical polar Fourier correlations[END_REF][START_REF] David | Recent progress and future directions in protein-protein docking[END_REF] and even in 5D rotation-translational space .

1.3.2.1.2 Sophisticated scoring during exhaustive search In modern algorithms based on the idea described above, a simple shape complementarity is usually coupled to the information about binding site [START_REF] Ben | Weighted geometric docking: Incorporating external information in the rotation-translation scan[END_REF], electrostatic energy [START_REF] Henry A Gabb | Modelling protein docking using shape complementarity, electrostatics and biochemical information[END_REF], atomic desolvation effect [START_REF] Chen | Docking unbound proteins using shape complementarity, desolvation, and electrostatics[END_REF], knowledge-based potentials [START_REF] Mintseris | Integrating statistical pair potentials into protein complex prediction[END_REF], etc. The general idea behind the incorporation of additional scoring approaches can be shown on the example of the ZDOCK program [START_REF] Mintseris | Integrating statistical pair potentials into protein complex prediction[END_REF]. The scoring function that is used in addition to shape complementarity in this algorithm contains the desolvation term, electrostatics and knowledge-based scoring potentials. To estimate the desolvation energy, the authors used atomic contact energies (ACE) [START_REF] Chen | Docking unbound proteins using shape complementarity, desolvation, and electrostatics[END_REF]. It is defined as the free energy change of breaking two protein atom contacts with water and forming a new contact between these two atoms. The pairwise shape complementarity (PSC) and ACE are represented on the grid in the following way: where R stands for receptor and L -for the ligand. The final score is:

R P SC = L P SC =     
S P SC+DE = Re [R P SC × L P SC ] + 1 2 Im [R DE × L DE ]
The electrostatic energy calculation is performed similarily:

R P SC+ELEC = L P SC+ELEC =     
3 surface of a protein 3 2 protein core 0 empty space Im [R P SC+ELEC ] = β × (electric potential of receptor) empty space 0 otherwise Im [L P SC+ELEC ] = -1 × (atom charge) grid point is the nearest grid point of a ligand atom 0 otherwise

Finally, the total score is computed as follows:

S P SC+DE+ELEC = Re [R P SC+ELEC × L P SC+ELEC ] + 1 2 Im [R DE × L DE ]
We see that using not only real part of the grid values but also their imaginary part, one can compute electrostatics, desolvation and shape complementarity using two convolution calculations. [START_REF] Venkatraman | Protein-protein docking using region-based 3D Zernike descriptors[END_REF], GAPDOCK [START_REF] Gardiner | Protein docking using a genetic algorithm[END_REF], PatchDock [START_REF] Duhovny | Efficient unbound docking of rigid molecules[END_REF], etc represent protein as a set of surface patches. Efficient search algorithms match the surface patches on two proteins and recover the rotation and position of a candidate docking pose from the patches.

1.3.2.1.5 Randomized search strategies An example of this type of approaches is the package RosettaDock [START_REF] Jeffrey | Protein-protein docking with simultaneous optimization of rigid-body displacement and side-chain conformations[END_REF]. It generates random starting positions of a ligand around receptor. Afterwards, it minimizes the score while moving the ligand along the line connecting their centers of masses. RosettaDock also uses different representations of the protein at two distinct stages of minimization to reduce the number of degrees of freedom: coarse-grained and full-atom ones. The same methodology is used by the programs like ATTRACT [START_REF] Zacharias | Protein-protein docking with a reduced protein model accounting for side-chain flexibility[END_REF], HADDOCK [START_REF] Dominguez | HADDOCK: a protein-protein docking approach based on biochemical or biophysical information[END_REF], etc. Another example of randomized search approach is the SwarmDock algorithm [START_REF] Li | Detection and refinement of encounter complexes for protein-protein docking: taking account of macromolecular crowding[END_REF]. This program uses population-based search algorithm, called particle swarm optimization. Each copy of the protein-protein complex being optimized is an agent. During optimization of each agent it shares information about its state with the neighbouring agents and they change the parameters of optimization accordingly. The particle swarm optimization algorithm is especially well suited for exploring energy landscapes with large number of local minima.

Scoring candidate conformations

Despite the vast variety of the methods to obtain the scoring functions, we can group them into two major classes -physics-based SFs and statistical SFs. The first class of SFs is constructed as a weighted sum of terms, such as desolvation [START_REF] Wang | Implicit solvent models for flexible protein-protein docking by molecular dynamics simulation[END_REF], electrostatic interactions [START_REF] Sheinerman | On the role of electrostatic interactions in the design of protein-protein interfaces[END_REF], hydrogen bonds [START_REF] Fernández | Insufficiently dehydrated hydrogen bonds as determinants of protein interactions[END_REF], hydrophobic interactions [START_REF] Scarsi | Hydrophobicity at the surface of proteins[END_REF], etc., given as E = i α i E i . Then, the weights α i are usually tuned to match some experiments or to attain a minimum of the SF on a set of known structures of protein complexes. On the other hand, statistical SFs are developed based on the observation that the distances between the atoms in experimentally determined structures follow the Boltzmann distribution [START_REF] Finkelstein | Why do protein architectures have boltzmann-like statistics?[END_REF]. More precisely, using ideas from statistical theory of liquids, effective potentials between atoms are extracted using the inverse Boltzmann relation:

E ij (r) = -k B T ln P ij (r) Z ,
where k B is the Boltzmann constant, P ij (r) denotes the probability to find two atoms of types i and j at a distance r, and Z denotes the probability distribution in the reference state. The latter is the thermodynamic equilibrium state of the protein when all interactions between the atoms are set to zero. The score of a protein conformation is then given as a sum of effective potentials between all pairs of atoms. Although this concept is old ( it originates from the work of Tanaka and Scheraga [START_REF] Tanaka | Medium-and long-range interaction parameters between amino acids for predicting three-dimensional structures of proteins[END_REF], Miyazawa and Jernigan [START_REF] Miyazawa | Estimation of effective interresidue contact energies from protein crystal structures: quasi-chemical approximation[END_REF] and Sippl [START_REF] Sippl | Calculation of conformational ensembles from potentials of mena force: an approach to the knowledge-based prediction of local structures in globular proteins[END_REF] ), it is still under debates [START_REF] Thomas | Statistical potentials extracted from protein structures: how accurate are they[END_REF][START_REF] Sippl | Helmholtz free energies of atom pair interactions in proteins[END_REF][START_REF] Wa Koppensteiner | Knowledge-based potentials-back to the roots[END_REF][START_REF] Ben-Naim | Statistical potentials extracted from protein structures: Are these meaningful potentials?[END_REF].

Particularly, the computation of the reference state is a challenging problem and only recently some attempts to rigorously justify and compute it have been made [START_REF] Hamelryck | Potentials of mean force for protein structure prediction vindicated, formalized and generalized[END_REF]. Some scoring functions from this class were obtained without the computation of the reference state. Among those we should mention SF obtained using linear programming [START_REF] Tobi | Distance-dependent, pair potential for protein folding: Results from linear optimization[END_REF][START_REF] Tobi | Optimal design of protein docking potentials: efficiency and limitations[END_REF][START_REF] Vendruscolo | Can a pairwise contact potential stabilize native protein folds against decoys obtained by threading?[END_REF], quadratic programming, support vector machines [START_REF] Bernauer | A new protein-protein docking scoring function based on interface residue properties[END_REF][START_REF] Hu | Developing optimal non-linear scoring function for protein design[END_REF][START_REF] Martin | Efficient comprehensive scoring of docked protein complexes using probabilistic support vector machines[END_REF], and iterative techniques [START_REF] Huang | An iterative knowledge-based scoring function for protein-protein recognition[END_REF][START_REF] Huang | Inclusion of solvation and entropy in the knowledgebased scoring function for protein-ligand interactions[END_REF].

Here I describe some examples of knowledge-based scoring functions, used in DFIRE [START_REF] Zhou | Distance-scaled, finite ideal-gas reference state improves structure-derived potentials of mean force for structure selection and stability prediction[END_REF], ATTRACT [START_REF] Zacharias | Protein-protein docking with a reduced protein model accounting for side-chain flexibility[END_REF] and ZDOCK [START_REF] Mintseris | Integrating statistical pair potentials into protein complex prediction[END_REF] algorithms. The approach chosen by the authors of DFIRE scoring function uses the ideal-gas reference state that allowed them to unify the folding and docking scoring functions. The pair distribution function has the following dependence on the number of observed pairs:

N obs (i, j, r) = 1 V N i N j g ij (r)4πr 2 δr
Where N obs (i, j, r) is the number of observed atoms of i-th and j-th types at the distance r. The potential of mean-force is connected to the pair distribution function: u(i, j, r) = -RT ln g ij (r). When the interaction u(i, j, r) is set to zero, one obtains the distribution in the reference state N exp (i, j, r) = 1 V N i N j 4πr 2 δr. However, due to the finite size of the protein, the correction coefficient α was introduced:

N exp (i, j, r) = 1 V N i N j 4πr α δr
The potential is assumed to have finite range and the equation for the potential is simplified by employing the pairwise distribution in the reference state at the potential cutoff distance:

N exp (i, j, r cut ) = N obs (i, j, r cut ) = 1 V N i N j 4πr α cut δr cut
Finally, the form of the potential is the following:

u(i, j, r) =    -νRT ln N obs (i,j,r) " r r cut " α " δr δr cut " N obs (i,j,rcut) , r < r cut 0 otherwise
where ν = 0.0157, R is the gas constant, T was set to 300K and α = 1.61. The δr and δr cut are the widths of bins at the distances r and r cut correspondingly. The coefficient ν was tuned to maximize correllation between the experimental and the predicted stability upon the mutations of monomeric proteins. The exponential factor α is determined using the uniformly distributed points in a sphere for each structure. The radii of the spheres were set to cR g , where R g is the gyration radius of a protein and c was tuned to equal the number of pairs within the r cut for the reference state and the experimental structures. Such a choice of the reference state implies that the information about the protein-protein contacts is negligible compared to the information on the interaction withing proteins. This fact let the authors to unify folding and docking scoring functions. D. Kozakov et. al. developed the "decoys as reference state" (DARS) potential for scoring protein-protein interactions. The key idea was to dock the proteins in the training set using only shape complementarity term in the scoring of conformations. These structures simulate the absense of interactions between two proteins and were used as the reference state. However due to computational complexity, the authors chose only 22 protein-protein complexes to derive the reference state.

Zacharias and his team used a different approach: they used only the Leenard-Jones type of interaction and the electrostatic interaction with distance-dependent dielectric constant = 15r. They estimated the parameters of interaction potential using the similar approach as in the work by Miyazawa and Jernigan [START_REF] Miyazawa | Self-consistent estimation of interresidue protein contact energies based on an equilibrium mixture approximation of residues[END_REF].

Knowledge-based potentials in rigid-body search

During the rigid-body search algorithms typically filter out about ≈ 10 4 conformations. Most of them are later refined by the scoring procedure. However, due to a simple scoring function used during the search, the conformations that are close to the solution could be underrepresented in the output of rigid-body docking algorithms. Therefore, two approaches exist to incorporate knowledge-based scoring functions into exhaustive 6D search procedure. One of them, used in the ZDOCK3.0 program [START_REF] Mintseris | Integrating statistical pair potentials into protein complex prediction[END_REF] integrates atomic contact potential into 6D search procedure in the following way. Suppose one has 2N atom types. For the ligand, N functions on a grid are defined in the following way:

Re [L i ] =
1 if grid cell is occupied by a ligand atom type i 0 otherwise

Im [L i ] = 1 if grid cell is occupied by a ligand atom type i+1 0 otherwise

The N functions for the receptor are: where e i,j is the value of the contact potential between atom types i and j and r cut is the contact radius. The sum of contacts looks as follows:

Im [R i ] = e i,
E = 2N i=1 2N j=1 e ij n ij = N k=1 x,y,z L i × R i
Thus, to compute the contact potential energy one has to perform N forward Fourier transforms and one backward (due to additivity of energy).

Despite the usage of both complex and real parts during the computation of contact potentials N can be around 10, which slows drastically the rigid-body docking algorithm. In order to reduce the number of Fourier transforms Kozakov et. al. [START_REF] Kozakov | PIPER: An FFT-based protein docking program with pairwise potentials[END_REF] proposed to decompose the interaction matrix e i,j into the eigenvectors: e i,j = P p λ p u p,i u p,j where P depends on the allowed error rate. Ususally the aproximation of the pairwise potential energy using the grid yelds 10% error in energy, therefore the truncation of the decomposition is well justified. The functions for the receptor and the ligand look like: R p = i u p,i Over all neighbours i within r cut 0 no neighbour atoms

L p = u p,j If atom of type j is in the cell 0 otherwise
This approach allows to compute only four Fourier correlations with the same error as in the previous algorithm.

Modelling of water molecules

An important part of a protein-protein interface constitute the water-mediated interactions. One pronounced example of the protein-protein complex where water molecules play great role is the barnase-barstar complex. In the interface of this complex 18 water molecules are fully burried mediating a considerable amount of sidechain-sidechain interactions [START_REF] Ashley M Buckle | Protein-protein recognition: Crystal structural analysis of a barnase-barstar complex at 2.0[END_REF]. Water molecules also play a role in the interaction between protein and drug-like molecules [START_REF] Ben-Naim | Molecular recognition-viewed through the eyes of the solvent[END_REF][START_REF] David | Systematic placement of structural water molecules for improved scoring of protein-ligand interactions[END_REF].

In spite of the importance of the water molecules prediction for the drug desing, numerous works are devoted to predict positions of water molecules around a known protein structure [START_REF] Forli | A force field with discrete displaceable waters and desolvation entropy for hydrated ligand docking[END_REF][START_REF] Wang | Ligand binding to protein-binding pockets with wet and dry regions[END_REF][START_REF] Gregory A Ross | Rapid and accurate prediction and scoring of water molecules in protein binding sites[END_REF]. However, the amount of papers attempting to predict the explicit solvation of protein-protein interfaces is substantially less [START_REF] Jiang | A "solvated rotamer" approach to modeling water-mediated hydrogen bonds at protein-protein interfaces[END_REF][START_REF] Huynh-Hoa | WATGEN: an algorithm for modeling water networks at protein-protein interfaces[END_REF][START_REF] Panagiotis L Kastritis | Solvated protein-protein docking using Kyte-Doolittlebased water preferences[END_REF][START_REF] Ahmad | Adhesive water networks facilitate binding of protein interfaces[END_REF].

Historically, the first approach to account for the water-mediated interactions during protein-protein docking was the use of the solvated rotamers [START_REF] Jiang | A "solvated rotamer" approach to modeling water-mediated hydrogen bonds at protein-protein interfaces[END_REF]. The key idea of this method is to attach the water molecules to the functional groups of the residues and the backbone and treat different solvation modes as rotamers. The water molecules placement around the protein was derived by the authors from X-Ray crystallographic structures solved at high resolution. To avoid combinatoric explosion of the number of solvated rotamers they restricted the placement of water molecules to the non-adjacent sites in a protein.

N N H W W

The other approach was implemented in the algorithm WATGEN [START_REF] Huynh-Hoa | WATGEN: an algorithm for modeling water networks at protein-protein interfaces[END_REF]. First, the hydrogens are added to the interacting proteins. These hydrogens, that interact with the atoms of a protein are called donors. Afterwards, all possible positions of water molecules are generated and those which have clashes with the atoms of the two proteins are discarded. Then, the water sites are selected based on the number of interactions they are involved into. Additionally, two water sites closer than 0.5 Å are considered equal. Finally, geometry of the water hydrogens is optimized to maximize the number of interactions.

Multimeric docking

Methods and algorithms described previously are applicable mainly to the problem of doking dimeric protein-protein complexes. However, many proteins in a cell form multimeric assemblies that play crucial role in recycling of the proteins in a cell, folding of the proteins, translation, transcription and many other vital cellular processes. There are two major ways in deciphering the structure of a multimeric complex starting from its subunits. The first type of programs do not rely on any other information except for the possible symmetry of a complex and the structures of its subunits. The second one uses low-resolution electron density maps, obtained from cryo-EM or small angle neutron scattering experiments.

The first class of methods deals with complexes of two general types: symmetric and nonsymmetrical. The number of packages that can do symmetry-based multimeric docking is quite small: SymmDock [START_REF] Schneidman-Duhovny | Geometry-based flexible and symmetric protein docking[END_REF], can build complexes with cyclic symmetry; Rosetta program protocol, that uses Monte-Carlo approach and can take into account cyclic, dihedral, helical and icosahedral symmetries [START_REF] André | Prediction of the structure of symmetrical protein assemblies[END_REF]; M-ZDOCK [101] reduces rotational search space assuming cyclic symmetry and relies on FFT-based docking approach, etc.

A few programs can dock proteins into nonsymmetrical complexes: CombDock [START_REF] Inbar | Prediction of multimolecular assemblies by multiple docking[END_REF] uses combinatorial approach, generating all pairwise pairs and finding among them tripples with the optimal score; Kim and Hummer algorithm [START_REF] Young | Coarse-grained models for simulations of multiprotein complexes: application to ubiquitin binding[END_REF], which uses Monte-Carlo approach and coarse-graining the protein models, HADDOCK [START_REF] Dominguez | HADDOCK: a protein-protein docking approach based on biochemical or biophysical information[END_REF]; ATTRACT [START_REF] Zacharias | Protein-protein docking with a reduced protein model accounting for side-chain flexibility[END_REF]; Multi-LZerD [START_REF] Esquivel-Rodríguez | Multi-LZerD: Multiple protein docking for asymmetric complexes[END_REF] and DockTrina [START_REF] Popov | DockTrina: Docking triangular protein trimers[END_REF].

However, despite the variety of methods, they still perform poorly even on a simple benchmark [START_REF] Popov | DockTrina: Docking triangular protein trimers[END_REF]. The most reliable and widely-used technique to obtain atomic structure of a multimeric assembly is the docking of subunits into a lowresolution electron density map, usually obtained from cryo-EM experiments.

Docking proteins into low-resolution electron density map

For this task, a number of software packages have been developed. Most notable of them are Situs [START_REF] Wriggers | Using Situs for the integration of multi-resolution structures[END_REF][START_REF] Chacón | Multi-resolution contour-based fitting of macromolecular structures[END_REF], NORMA [START_REF] Suhre | NORMA: a tool for flexible fitting of high-resolution protein structures into low-resolution electron-microscopyderived density maps[END_REF], EMFit [START_REF] Rossmann | Combining electron microscopic with X-ray crystallographic structures[END_REF], UROX [START_REF] Siebert | UROX 2.0: an interactive tool for fitting atomic models into electron-microscopy reconstructions[END_REF], etc. Despite the differences in the implementation, all algorithms maximize some score that shows the goodness of the fitting using a certain optimization algorithm. An excellent review on different types of the scoring functions used for cryo-EM density fitting is given by Vasishtan and Topf [START_REF] Vasishtan | Scoring functions for cryoEM density fitting[END_REF]. According to them, one of the most popular scoring functions is the cross-correlation function (CCF) between the EDM and the density of the fitted protein.

Given a protein structure that is described by its electron density f (r), and an EDM obtained from e.g. a cryo-EM experiment described by the function g(r), we can minimize the square root discrepancy between them. Precisely, this discrepancy is given by

S = dr T Rf (r) -g(r) 2 , (1.3)
where T and R are the operators of the translation and the rotation respectively, applied to the density f (r). We can rewrite the scoring function S as

S = dr T Rf (r) 2 + dr g 2 (r) -2 dr T Rf (r)g(r) (1.4)
Therefore, the minimization of the score S is equivalent to the maximization of the CCF:

CCF = dr T Rf (r)g(r) (1.5)
with respect to the parameters of the operators T and R. This scoring function has been used in the majority of the algorithms and software packages that perform the fitting into the EDM [START_REF] Wriggers | Using Situs for the integration of multi-resolution structures[END_REF][START_REF] Siebert | UROX 2.0: an interactive tool for fitting atomic models into electron-microscopy reconstructions[END_REF][START_REF] Suhre | NORMA: a tool for flexible fitting of high-resolution protein structures into low-resolution electron-microscopyderived density maps[END_REF]. Another widely used scoring function is the Laplacian-filtered cross-correlation function (LCCF). It originated from the observation that a human performing a manual fitting a structure into an EDM tends to match the isosurfaces of the densities rather than the densities themselves,

LCCF = dr T Rf (r) ( g(r)) (1.6)
This scoring function works better than CCF for low resolution maps (∼ 10 -30 Å) [START_REF] Wriggers | Using Situs for the integration of multi-resolution structures[END_REF] and was used for the first time in the CoAn/CoFi algorithm [START_REF] Volkmann | Quantitative fitting of atomic models into observed densities derived by electron microscopy[END_REF]. Other scoring functions that e.g. penalise symmetry-induced protein-protein contacts, or make use of protein-protein docking potentials, etc., have also been developed [START_REF] Vasishtan | Scoring functions for cryoEM density fitting[END_REF].

Aim of work

This work deals with the protein-protein docking problem. The aim of the work was to propose and validate new algorithms in the area of 6D exhaustive search and in the field of ranking the docking predictions. Particularly, the goal was to invent and test new ways to solve general global rigid-body search approaches that could be more effective than the state-of-art in the field. In order to convincingly demonstrate their applicability, they have to be applied to one of the challenging problems, like fitting of cryo-EM electron density maps or rigid-body protein-protein docking and compared to the existing programs in the field. The scoring method has to be applicable to a wide range of problems: from drug-like molecules docking to crystallographic water prediction. It has to be free of methodological difficulties as in the case of statistical potentials, where the debates about the validity of the reference state computations are still going. This algorithm should be derived from the basic logical statements and have acceptable convergence properties.

Chapter 2

Hermite fitting 2.1 Introduction

As was already mentioned in the introduction, an important class of algorithms in computer science and structural biology deals with the exhaustive search in the six-dimensional space of translations and rotations of a rigid body. Modern exhaustive search algorithms either implement the fast 3D translational search using the fast Fourier transform (FFT) [START_REF] Chacón | Multi-resolution contour-based fitting of macromolecular structures[END_REF][START_REF] Katchalski-Katzir | Molecular surface recognition: determination of geometric fit between proteins and their ligands by correlation techniques[END_REF][START_REF] Gabb | Modelling protein docking using shape complementarity, electrostatics and biochemical information[END_REF][START_REF] Wriggers | Using Situs for the integration of multi-resolution structures[END_REF][START_REF] Siebert | UROX 2.0: an interactive tool for fitting atomic models into electron-microscopy reconstructions[END_REF] or the fast 3D rotational search by means of the spherical harmonics decomposition and the FFT [START_REF] Kovacs | Fast rotational matching[END_REF] or even the fast 5D rotational search [START_REF] Kovacs | Fast rotational matching of rigid bodies by fast Fourier transform acceleration of five degrees of freedom[END_REF][START_REF] David W Ritchie | Accelerating and focusing protein-protein docking correlations using multi-dimensional rotational FFT generating functions[END_REF]. Exhaustive search is also widely used as a preliminary step preceding the local search or flexible refinement procedures. Thus, the quality and the speed of the exhaustive search algorithms have a great impact on the solution of the vast variety of problems. Therefore, we believe that new directions of research on this topic are very important and highly valuable.

In this section, we present the new HermiteFit algorithm that uses the orthogonal Hermite functions to perform exhaustive search in the six-dimensional space of rigid-body motions. We apply this method to the problem of fitting of a high resolution X-ray structure of a protein subunit into the cryo-electron microscopy (cryo-EM) density map of a protein complex. As a part of the new method, we developed an algorithm for the rotation of the decomposition in the Hermite basis and another algorithm for the conversion of the Hermite expansion coefficients into the Fourier basis. We demonstrate the ability of our algorithm to compete with the well-established approaches by using two examples of different difficulty, the PniB conotoxin peptide and the GroEL complex. The first example illustrates encoding principles and demonstrates the influence of the encoding quality on the goodness of fit. The second example is the gold standard of all electron density map fitting algorithms. Our approach allows to analytically assess the quality of encoding of the Hermite basis using an estimation of the crystallographic R-factor. We then compare this estimation with the one computed numerically for the PniB conotoxin density map. Finally, we compare the speed and the fitting accuracy of our algorithm with the two popular programs, the ADP EM fitting method and the colores program from the Situs package and demonstrate that HermiteFit spends less running time per one search point compared to the two other methods while attaining a similar accuracy.

The HermitFit algorithm can be straightforwardly applied to a broad class of problems in different fields of research. For example, one of the bottlenecks of the algorithms for molecular replacement in crystallography is the computation of the Fourier coefficients (structure factors) of a molecule [START_REF] Navaza | On the fast translation functions for molecular replacement[END_REF]. This operation is to be precise and fast. However, the exact analytical evaluation of the structure factors is too costly [START_REF] Sayre | The calculation of structure factors by Fourier summation[END_REF] when recomputing them for each rotation of the molecule. Therefore, currently one uses the Sayre-Ten Eyck approach to compute the Fourier coefficients [START_REF] Eyck | Efficient structure-factor calculation for large molecules by the fast Fourier transform[END_REF]. Unfortunately, one has to be very careful tuning the parameters of the electron density model and the grid cell size to obtain the desired precision [START_REF] Navaza | On the computation of structure factors by FFT techniques[END_REF][START_REF] Pv Afonine | On a fast calculation of structure factors at a subatomic resolution[END_REF]. Unlike the Sayre-Ten Eyck, our algorithm offers the analytical expression for the structure factors of the Hermite decomposition of a molecule. Finally, our approach allows to analytically estimate the quality of encoding using, e.g., crystallographic R-factors.

Methods

Summary of the standard fitting algorithm

The standard FFT-based 3D fitting algorithm operates according to the workflow shown in Figure 2.2.1 [START_REF] Katchalski-Katzir | Molecular surface recognition: determination of geometric fit between proteins and their ligands by correlation techniques[END_REF][START_REF] Gabb | Modelling protein docking using shape complementarity, electrostatics and biochemical information[END_REF][START_REF] Chacón | Multi-resolution contour-based fitting of macromolecular structures[END_REF]. The input of this algorithm is a protein atomic structure determined experimentally by, e.g., X-ray crystallography or nuclear magnetic resonance (NMR) experiments. Another input is an experimental EDM determined by means of, e.g., cryo-EM. First, the algorithm decomposes the experimental EDM into the Fourier basis using the fast Fourier transform algorithm. Then, it rotates the protein structure to a certain orientation r and decomposes the electron density of the rotated structure into the Fourier basis.

Sort

The electron density is typically computed as a sum of Gaussians centred on non-hydrogen atoms of the protein. Afterwards, the algorithm exhaustively explores translational degrees of freedom of the rotated protein with respect to the EDM. For every translation t, it determines the corresponding score, which is usually given by the correlation between the two densities. This procedure is equivalent to computing the convolution of two functions,

CCF(r, t) = dx f (r, x -t)g(x), (2.1)
where f (r, xt) is the density of the protein rotated by r and translated by t, and g(x) is the experimental electron density map. To speed up this step, the algorithm computes the values of the Fourier transform of the CCF for all translational degrees of freedom at once, using the convolution theorem. Finally, the algorithm computes the inverse Fourier transform (IFT) of the convolution, generates a new rotation of the protein structure, and returns to the second step. This procedure is repeated until all rotational degrees of freedom of the protein with respect to the EDM are explored (see Fig. 2.2.1). The solution of the fitting problem is then given by (r max , t max ) = argmax r,t {CCF(r, t)}.

The bottleneck of the standard algorithm is the re-projection of the protein electron density into the Fourier space after each rotation. To overcome it, we propose to encode the electron density of the protein structure in the orthogonal Hermite basis, prior to performing the rotational search. This allows to speed up the Operation

Complexity

Loop multiplier Decomposition of the step function

O(M 3 logM 3 ) 1 Decomposition of the Gaussian O(N atoms N 3 ) Construction of the rotation matrix O(N rot N 4 ) Rotation O(N 4 ) N rot Evaluation of the Hermite series O(M 3 • N + M 2 • N 2 + M • N 3 ) Multiplication O(M 3 ) Inverse Fourier Transform O(M 3 logM 3 )
Table 2.1: Complexity of the Hermite fitting algorithm. Here, M denotes the order of the Fourier decomposition; N is the order of the Hermite decomposition; N atoms is the number of atoms in the protein; N rot -the number of rotations to be sampled.

projection of the protein density into the Fourier space. Since only the members of the Fourier family of linear transforms can replace O(N 2 ) operations of a convolution in a time domain by O(N ) operations in a frequency domain [START_REF] Stone | Convolution theorems for linear transforms[END_REF], we still need to perform the convolution in the Fourier space.

Hermite functions

Orthogonal Hermite function of order n is defined as:

ψ n (x; λ) = √ λ 2 n n! √ π exp(- λ 2 x 2 2 )H n (λx), (2.2)
where H n (x) is the Hermite polynomial and λ is the scaling parameter. In Fig.

we show several orthogonal

Hermite functions of different orders with different parameters λ. These functions form an orthonormal basis set in L 2 (R). A 1D function f (x) decomposed into the set of 1D Hermite functions up to an order N reads

f (x) = N i=0 fi ψ i (x; λ) (2.3)
Here, fi are the decomposition coefficients, which can be determined from the orthogonality of the basis functions ψ i (x; λ). Decomposition in Eq. 2.3 is called the band-limited decomposition with ψ i (x; λ) basis functions. To decompose the EDM and the protein structures, we employ the 3D Hermite functions:

ψ n,l,m (x, y, z; λ) = ψ n (x; λ)ψ l (y; λ)ψ m (z; λ), (2.4)
which form an orthonormal basis set in L 2 (R 3 ). A function f (x, y, z) represented as a band-limited expansion in this basis reads

f (x, y, z) = N i=0 N -i j=0 N -i-j k=0 fi,j,k ψ i,j,k (x, y, z; λ) (2.5) Ψ n (x, λ) -1 -0.5 0 0.5 1 x, Å -6 -4 -2 0 2 4 6 Ψ10(x,1.0) Ψ4(x,1.0) Ψ n (x, λ) -1 -0.5 0 0.5 1 x, Å -6 -4 -2 0 2 4 6 Ψ6(x,1.0) Ψ6(x,2.0) Ψ6(x,3.0)

Decomposition of electron densities into the orthogonal Hermite basis

One of the advantages of the orthogonal Hermite basis is that we can derive the exact analytical expression for the decomposition coefficients of a molecular structure. This allows to rapidly obtain the exact decompositions without costly numerical integration over the 3D space. In our algorithm, the electron density of the protein (f (x) in Eq. 2.1, upon which rotation and translation operators act) is expanded in the Hermite basis using the Gaussian model. More precisely, we model the electron density of a single atom in the molecular structure as a Gaussian centred at the atomic position r

(i)
0 with the squared variance equal to α 2 /2. Then, the electron density of the whole molecular structure is given by the following sum:

M (r) = Natoms i=1 e -|r-r (i) 0 | 2 /α 2 , (2.6)
where r

(i)
0 is the position of the i-th atom, α/ √ 2 is the variance of the Gaussian distribution, and r = (x, y, z) ∈ R 3 is the sampling volume. Normally, each Gaussian should be weighted with a coefficient corresponding to electron distribution of a particular atom. However, we omit the weights in our approximation. In the section 2.2.4, we provide analytical expressions (Eqs. 2.8 and 2.14) for the decomposition coefficients of M (r) in the 1D and the 3D cases.

Shifted Gaussian expansion

Here we provide the derivation of the expansion coefficients of a shifted Gaussian of the following form:

g(r) = e -|r-r 0 | 2 α 2 (2.7)
into the orthogonal Hermite basis. The well known property of this basis (as well as of any orthogonal basis) is the following:

if f (x, y, z) = f (1) (x)f (2) (y)f (3) (z) and f (k) (t) = N i=0 f (k) i ψ i (t; λ) then fi,j,k = f (1) i f (2) j f (3) k (2.8)
First, we derive the decomposition of a 1D Gaussian into the 1D orthogonal Hermite basis. Then, using property (2.8) we obtain the decomposition of a 3D Gaussian into the 3D orthogonal Hermite basis. More specifically, the 1D Gaussian function reads as:

g(x) = e -(x-ξ) 2 α 2
(2.9)

Its decomposition coefficients are equal to:

ĝn (ξ; λ, α) = g(x)ψ n (x; λ) dx = n! √ λe -ξ 2 α 2 " 1-1 α 2 β 2 " 2 n n! √ π [ n 2 ] m=0 (-1) m m!(n -2m)! e -β 2 " x-ξ α 2 β 2 " 2 (2λ(x - ξ α 2 β 2 ) + 2λξ α 2 β 2 ) n-2m dx, (2.10)
where

β 2 = λ 2 2 + 1 α 2 .
From now on we will, for brevity, write ĝn instead of ĝn (ξ; λ, α). Changing the variables t = x -ξ α 2 β 2 and denoting a = ξ α 2 β 2 , we obtain:

ĝn = n! √ λe -ξ 2 α 2 " 1-1 β 2 " 2 n n! √ π [ n 2 ] m=0 (-1) m (2λ) n-2m m!(n -2m)! e -β 2 t 2 (t + a) n-2m dx (2.11)
Next, we decompose the sum (t + a) k using Newton's formula:

(t + a) k = k i=0 k i t i a k-i (2.12)
Thus, the integral in Eq. 2.11 will read:

e -β 2 t 2 (t + a) n-2m dx = n-2m i=0, i-even (n -2m)! 2 i i 2 !(n -2m -i)! √ πβ -1-i a n-2m-i (2.13)
Substituting it to the formula for ĝn and denoting n-2m i=0, i-even =

[ n-2m 2 ] l=0
(i = 2l), we obtain the following expression for the coefficients:

ĝn (ξ; λ, α) = e -ξ 2 α 2 " 1-1 α 2 β 2 " n! √ πλ 2 n [ n 2 ] m=0 [ n-2m 2 ] l=0 (-1) m 2 n-2m-2l λ n-2m l!(n -2m -2l)!m! β -2n+4m+2l-1 ξ α 2 n-2m-2l
(2.14)

Finally, using Eq. 2.8 we obtain a decomposition of the 3D Gaussian into the 3D Hermite basis. We should note that in order to avoid the rounding error, one should begin the summation with the Gaussians that are located father from the origin.

Expansion of a function defined on a grid

In many docking algorithms the parwise interaction of particles is approximated as the set of functions on the grid. Therefore in many important cases the protein description goes beyond the sum of gaussians as in Eq. 2.6. In this section we provide a way to directly obtain decomposition of a function f (x, y, z) defined on a regular grid. We can represent this function as a sum:

f (x, y, z) = i,j,k f (x i , y j , z k )η ijk (x, y, z)
where η ijk (x, y, z) is a step-function in the position of ijk-th grid cell. The expression for function is the following:

η(x) = 1, 0 ≤ x < h 0, otherwise (2.15)
The decomposition coefficients of the shifted η read: [START_REF] Bernauer | A new protein-protein docking scoring function based on interface residue properties[END_REF])

η(x -a) = N i α i ψ i (x; λ) (2.
α i = +∞ -∞ η(x -a)ψ i (x; λ) dx = h+a a ψ i (x; λ) dx (2.17)
We can obtain the coefficients α n by simple integration:

α n = 1 2 n n! √ πλ λ(h+a) λa exp(- t 2 2 )H n (t) dt (2.18)
Using the well known result for the finite integral of the Hermite polynomials:

z q-1 e -pz 2 H n (z)dz = n! [ n 2 ] k=0 (-1) k-1 2 n-2k-1 z n-2k+q (pz 2 ) k-n+q
with parameters q = 1 and p = 1/2 we obtain:

e -z 2 2 H n (z)dz = (2.20) n! (Sgn[z]) n+1 [ n 2 ] k=0 (-1) k-1 Γ " n+1 2 -k, z 2 2 " 2 3k-3n 2 + 1 2 k!(n-2k)! + C (2.21)
After deriving the missing coefficients C n we obtain the following expression for the coefficients α n :

α n = 2 n √ n! 2 √ πλ [ n 2 ] k=0 (-1) k-1 2 √ 2 2k k!(n -2k)! (Sgn[h + a]) n+a Γ n + 1 2 -k, (λh + λa) 2 2 -(Sgn[a]) n+a Γ n + 1 2 -k, (λa) 2 2 - (Sgn[h + a]) n+1 -(Sgn[a]) n+1 Γ n + 1 2 -k, 0
Next, we are moving to the 3D case:

η(x, y, z) = 1, 0 ≤ x < h x 0 ≤ y < h y 0 ≤ z < h z 0, otherwise η(x, y, z) = η(x; h x )η(y; h y )η(z; h z ) η(x -a x , y -a y , z -a z ; h x , h y , h z ) = N k N j N i α i α j α k ψ i (x; λ)ψ i (y; λ)ψ k (z; λ)
As we see we have to multiply the coefficients of individual 1D-function with the shifts corresponding to the grid cell position along individual axes. We begin summation staring from the cells that farther from the frame origin because rounding error substantially influences the final result.

Laplacian filter in the Hermite basis

For mid-to low-resolution maps the Laplacian-filtered cross-correlation function gives a better match compared to the CCF [START_REF] Wriggers | Using Situs for the integration of multi-resolution structures[END_REF]. In the Hermite basis, the Laplacian filter has a particularly simple form. Using the well-known recurrence relation for the derivatives of the Hermite functions, we can easily derive the following relation for the second derivative of a 1D basis function:

d 2 dx 2 ψ n (x; λ) = λ 2 2 n(n -1)ψ n-2 (x; λ) + (2n + 1)ψ n (x; λ) + (n + 1)(n + 2)ψ n+2 (x; λ)
(2.22) A similar relationship holds for the coefficients of the decomposition:

ĥ n = λ 2 2 n(n -1) ĥn-2 + (2n + 1) ĥn + (n + 2)(n + 1) ĥn+2 , (2.23)
where ĥn and ĥ n are the n-th order decomposition coefficients of the original basis and its Laplacian representation, respectively. For n < 0 and n > N we let ĥn = 0 and ĥ n = 0. Due to the properties of the Laplace operator and the 3D Hermite decomposition, the contribution of the derivatives along each axis are additive. The derivation of the formula for the 3D decomposition derivative is straightforward and we omit it for brevity.

Rotation of the Hermite decomposition

Recently, Park et al. [START_REF] Park | Accurate image rotation using Hermite expansions[END_REF] presented the method to perform an in-plane rotation of a 2D orthogonal Hermite band-limited decomposition. Here, we extend their method for the 3D case. Let us first consider the decomposition of a 2D function into a 2D orthogonal Hermite function basis:

f (x, y) = N n=0 N -m m=0 fn,m ψ n (x; λ)ψ m (y; λ) (2.24)
The decomposition of a function f θ (x, y) rotated clock-wise by an angle θ reads

f θ (x, y) = N m=0 m k=0 ( m n=0 fn,m-n S m k,n )ψ k (x; λ)ψ m-k (y; λ), (2.25)
where coefficients S m k,n are computed using the following recurrent formulas [START_REF] Park | Accurate image rotation using Hermite expansions[END_REF]:

S m+1 q,n = n m -q + 1 sin(θ)S m q,n-1 + m -n + 1 m -q + 1 cos(θ)S m q,n S m+1 q,0 = m + 1 m -q + 1 cos(θ)S m q,0 S m+1 m+1,n = n m + 1 cos(θ)S m m,n-1 - m -n + 1 m + 1 sin(θ)S m m,n S m+1 m+1,0 = -sin(θ)S m m,0
(2.26)

The key idea that allows to generalize these formulas to a 3D decomposition is that we can factorize a rotation in 3D space into 3 independent in-plane rotations around three different axes, and then rotate each 2D decomposition using Eq. 2.25. Let us consider the following 3D decomposition:

f (x, y, z) = N n=0 ψ n (x; λ) N -n m=0 N -m-n l=0 fn,m,l ψ m (y; λ)ψ l (z; λ) (2.27)
If we rotate this decomposition about x axis, this rotation will be equivalent to N rotations of different 2D decompositions in the yz-plane:

f n (y, z) = N -n m=0 N -m-n l=0 fn,m,l ψ m (y; λ)ψ l (z; λ) (2.28)
This observation means that in order to perform such rotation, we need to recompute rank-3 tensor of coefficients fn,m,l slice by slice N times using Eq. 2.25. Figure 2.2.7 illustrates three subsequent rotations of tensor fn,m,l . Each rotation of the coefficients in one plane corresponds to a multiplication of these coefficients with a rotation matrix. Therefore, a 3D rotation defined with three Euler angles is equivalent to three sequential rotations of coefficients in three planes.

Transition from the Hermite to the Fourier basis

In order to perform a fast convolution as in Eq. 2.1, we convert the decomposition coefficients from the Hermite basis into the Fourier basis. This allows to use the fast convolution algorithm based on the Fourier convolution theorem, which was first introduced in protein-protein docking studies [START_REF] Katchalski-Katzir | Molecular surface recognition: determination of geometric fit between proteins and their ligands by correlation techniques[END_REF][START_REF] Gabb | Modelling protein docking using shape complementarity, electrostatics and biochemical information[END_REF] and then also applied in the EDM fitting [START_REF] Chacón | Multi-resolution contour-based fitting of macromolecular structures[END_REF][START_REF] Wriggers | Using Situs for the integration of multi-resolution structures[END_REF][START_REF] Siebert | UROX 2.0: an interactive tool for fitting atomic models into electron-microscopy reconstructions[END_REF]. The key idea of this algorithm is to compute the Fourier transform of the values of a scoring function on a grid, CCF(r, t) = f (r, x)g(r, xt)dx, using the convolution theorem:

F [f * g] = F [f ]F [g], (2.29)
i.e. to multiply the complex conjugated coefficients of the Fourier transform of the protein electron density with the coefficients of the Fourier transform of the EDM. Then, we obtain CCF(r, t) by taking the inverse Fourier transform of

F [f * g], CCF(r, t) = IFT F [f ]F [g] (2.30)
Now we explain how we convert the decomposition coefficients from the Hermite basis into the Fourier basis. Consider the decomposition of a function f (r) in the 3D Hermite basis with the decomposition coefficients fi,j,k (Eq. 2.5). Orthogonal Hermite functions are the eigenfunctions of the continuous Fourier transform:

ψ n (x; λ)e -2πiωx dx = (-i) n ψ n (ω; 2π λ ) ≡ ψn (ω; λ), (2.31)
where ω is the frequency in the reciprocal space. In order to compute Fourier coefficients of f (r) up to order M , we first compute the Fourier transforms of the basis functions ψ i (x; λ), ψ j (y; λ), and ψ k (z; λ) using Eq. 2.31. After, we substitute these coefficients into Eq. 2.5 and obtain the following expression for fl,m,n , the Fourier coefficients of f (r):

fl,m,n = 1 L x L y L z N i=0 N -i j=0 N -i-j k=0 fi,j,k ψi ( l L x ; λ) ψj ( m L y ; λ) ψk ( n L z ; λ) (2.32)
These values can be computed in

O(M 3 • N + M 2 • N 2 + M • N 3
) steps (see section 2.2.9).

Fast summation

Here we explain the fast summation in Eq. 2.33:

fl,m,n = N i=0 N -j j=0 N -i-j k=0 fi,j,k ψi,l ψj,m ψk,n , (2.33)
with indexes l, m, n ∈ [0, M ]. The summation in this formula can be performed with less operations than a naive estimation O(M 3 N 3 ) suggests. We perform the fast summation by splitting the equation into three consecutive sums:

T 1 i,j,n = N -i-j k=0 fi,j,k ψk,n (2.34)
T 2 i,m,n = N -i j=0 T 1 i,j,n ψj,m (2.35) fl,m,n = N i=0 T 2 i,m,n ψi,l (2.36)
It is easy to see that the construction of T 1 i,j,n matrix takes O(M N 3 ) operations, the construction of T 2 i,m,n matrix takes O(M 2 N 2 ) operations, and the final summation takes O(M 3 N ) operation. In the common use case (N = 15, M N ) the last sum takes much more time than the other two. To optimize it, we used the Gauss method to multiply complex numbers and expressed the whole sum as a generalized matrix product of three real-valued matrices. To implement these operations, we used the ATLAS library.

Implementation details and running time

We chose to demonstrate the potential of the Hermite basis by implementing the rigid-body fitting of an atomistic structure of a protein in an electron density map of low resolution. The HermiteFit algorithm was implemented using the C++ programming language and compiled using g++ with -O3 optimization. The running times of the tested algorithms are measured on a single core of an Intel ® Xeon ® CPU X5650 @ 2.67GHz processor with 24 GB of RAM on a Linux 64-bit operating system.

Our fitting method typically samples some 10 10 rigid-body configurations. Therefore, it is practical to group its fitting solutions into clusters. There are multiple ways to measure the similarity between rigid-body solutions. For example, the pair-wise root-mean-square deviation (RMSD) is a fast and well-accepted similarity measure. Thus, we clustered the fitting solutions using the rigid-body clustering algorithm implemented with the RigidRMSD library [START_REF] Popov | Rapid determination of RMSDs corresponding to macromolecular rigid body motions[END_REF] as follows. First, the fitting solution with the best score (yet unassigned to any cluster) is taken as the seed for the new cluster. Second, the pair-wise RMSDs between the seed and all other predictions are measured and the predictions with the RMSD lower than a certain threshold are put into the cluster. Finally, these two steps are iterated until all fitting predictions are assigned to corresponding clusters.

Analysis

This section provides analytical and numerical analysis of the density encoding in the Hermite basis. More specifically, we provide the choice of optimal model parameters and assess the quality of encoding.

Choice of parameters of the method

Orthogonal Hermite functions (2.2) decay exponentially after a certain distance and thus can encode information only within some interval. We can estimate this interval using the formula for the last root of a Hermite polynomial, ξ 1,N ≈ √ 1+2N λ [START_REF] Ricci | Improving the asymptotics for the greatest zeros of Hermite polynomials[END_REF], which gives an approximation for the half-size of the bounding box that we can successfully encode:

L box /2 √ 1 + 2N λ (2.37)
On the other hand, orthogonal Hermite functions are the eigenfunctions of the continuous Fourier transform (Eq. 2.31). Therefore, Hermite decomposition of order N can encode only a certain interval of frequencies. Using the same approximation as in the case of the real-space interval, we obtain the following equation for the maximum encoding frequency:

ω max = λ 2π √ 2N + 1 (2.38)
In case of the the Fourier series expansion on an interval (0, L box ), we can use the same estimation for the maximum encoding index M max by setting M max = 2L box ω max . Resolution R of an X-Ray electron density map is defined by the size of the reciprocal lattice as R = 1/(2ω max ), or, equivalently, R = L box /M max . Therefore, using resolution of the map R and the order of the Fourier series expansion M , we can estimate the lower bound on the Hermite scaling parameter λ required to encode all the reflexes of the electron density diffraction pattern to be

λ π max(R, L box /M ) √ 2N + 1 (2.39)
Here, we bounded the actual resolution by L box /M , because this will be the limit allowed by the finite Fourier series of order M . The two inequalities (2.37 and 2.39) give approximate bounds on the scaling parameter λ, provided that we know the size of the box L box containing a protein density and the resolution of the map R. Using these inequalities, we obtain the following relationship between parameters λ and N :

π √ 2N + 1 max(R, L box /M ) λ 2 √ 1 + 2N L box , (2.40)
which is valid for sufficiently large values of N . Nonetheless, we can use the following empirical estimation for the optimal value of λ at any N :

λ opt ≈ π 2 max(R, L box /M ) √ 2N + 1 + √ 1 + 2N L box (2.41)
Using dimensionless relative parameters λL box and L box /R we may rewrite the previous expression as

λL box ≈ π min(L box /R, M ) 2 √ 2N + 1 + √ 1 + 2N (2.42)
If at a given expansion order N there is no such parameter λ that satisfies inequality (2.40), then the protein representation might involve information loss. Therefore, we can estimate the minimum order N min of the Hermite expansion that allows this inequality to have solutions to be

N min ≈ π 4 min L box R , M (2.43)
Validity of the provided estimates and the graphical representation of the realspace and the reciprocal-space bounds on parameter λ will be demonstrated in the following sections.

The maximum order of the Fourier expansion M max can be estimated from the resolution and the size of the density map as R = L box /M max . However, when finding the global maximum of the cross-correlation function, we need to sample the space of possible translations of a protein with respect to the EDM with a step several times finer than the EDM resolution R. In protein crystallography, it is the common practice to set the sampling step size to R/3 [START_REF] Pv Afonine | On a fast calculation of structure factors at a subatomic resolution[END_REF]. In principle, we can use the same reasoning in choosing the optimal number of rotations N rot . When using spherical harmonics, the angular search step usually equals to the resolution of the basis, 2π/N [START_REF] José | ADP EM: fast exhaustive multi-resolution docking for high-throughput coverage[END_REF]. In case of the Hermite basis, we propose to use the same criterion.

The transfer matrix

Below we describe an analytical model of encoding by the Hermite basis for the one-dimensional case. Suppose we have a function f (x) that describes an electron density of a non-periodic object. Without loss of generality, we assume that this function is defined on a 1D interval of (-L box /2; +L box /2). This function has the following decomposition into Fourier series:

f exact k = 1 L box +L box /2 -L box /2 f (x)e -2πikx/L box dx (2.44)
We will refer to Fourier coefficients obtained using this expression as exact. The original function is then recovered by the inverse Fourier transform:

f (x) = +∞ k=-∞ f exact k e 2πikx/L box (2.45)
On the other hand, our algorithm computes approximate Fourier coefficients using the Hermite to Fourier transform:

f approx k = 1 L box N n=0 fn ψn ( k L box ; λ) (2.46)
Assuming that function f (x) is zero outside of the bounding interval, Hermite coefficients fn can be written as the finite integral:

fn = +L box /2 -L box /2 f (x)ψ n (x; λ) dx (2.47)
Now, we can express the approximate Fourier coefficients as a linear combination of the exact ones:

f approx k = +∞ l=-∞ T k,l f exact l , (2.48)
where the transfer matrix T k,l reads as:

T k,l = 1 L box N n=0 ψn (k; λ) +L box /2 -L box /2 ψ n (x; λ)e 2πilx/L box dx (2.49)
The transfer matrix acts as a linear filter in the reciprocal space and demonstrates how the input function is distorted by the finite size N of the Hermite basis. We should note that, generally, its values are complex numbers. This matrix can also be seen as a product of two matrices,

T = F (1) F (2) , (2.50)
where the first matrix is a scaled Fourier transform of the basis functions, F

kn = ψn (k; λ)/ L box (2.51) (1)
and the second matrix is a scaled Fourier series of the basis functions,

F (2) nl = +L box /2 -L box /2 ψ n (x; λ)e 2πilx/L box dx/ L box (2.52)
Figure 2.3.2 shows the absolute values of matrices F (1) and F (2) computed with λ = 0.55 and L box = 23 Å. The values of the Fourier series F (2) were computed numerically using adaptive quadrature. The dashed blue line shows the maximum encoding frequency ω max , according to Eq. 2.38, and bounds the encoding region. The solid black line on the right plot demonstrates the maximum order of the Hermite expansion (Eq. 2.37), after which the Fourier series encode mainly the frequencies near ω max . This is because on a finite interval (-L box /2, +L box /2), high-order Hermite basis functions become orthogonal to low-order Fourier basis functions. , the scaled Fourier transform of a 1D Hermite function as given by Eq. 2.51. Right: F (2) , the scaled Fourier series of a 1D Hermite function as given by Eq. 2.52. The dashed blue line highlights the maximum encoded frequency according to Eq. 2.38. The solid black line on the right plot shows the maximum Hermite decomposition order N max , at which the two matrices are still identical (Eq. 2.37). We see, however, that the transfer matrix at small values of λ encodes only loworder reflexes. The index of the last encoded reflex can be estimated from Eq. 2.38 as

k max = √ 2N + 1λL box /(2π).
With the increase in order N and parameter λ, the number of encoded frequencies rises. At the same time, increasing the scaling parameter λ makes the quality of encoding of all the frequencies worse, as we see in the right column. Therefore, it is very important to tune the value of λ according to the class of input functions, such that the quality of encoding becomes optimal. Below we will assess encoding quality by means of the crystallographic R-factor. . Increasing the value of λ, the number of encoded frequencies rises. However, at the same time, the quality of encoding of low frequencies worsens, as can be seen from the values at the diagonal.

Asymptotic behaviour of the transfer matrix

Here we demonstrate that the transfer matrix asymptotically achieves the Kronecker delta function at N → ∞. Recall the Mehler's formula [START_REF] Mehler | Ueber die Entwicklung einer Function von beliebig vielen Variablen nach Laplaceschen Functionen höherer Ordnung[END_REF]:

N n=0 u n ψ n (x)ψ n (y) = 1 π(1 -u 2 ) exp - 1 -u 1 + u (x + y) 2 4 - 1 + u 1 -u (x -y) 2 4 (2.53)
If we rewrite the transfer matrix in the following way:

T k,l = 1 L box +L box /2 -L box /2 dx N n=0 (-i) n ψ n ( k L box ; 2π λ )e 2πilx/L box ψ n (x; λ), (2.54)
and use the fact that

ψ n (x; λ) ≡ √ λψ n (λx), (2.55)
we see that we can use the Mehler's formula to compute the limit lim

N →∞ N n=0 (-i) n ψ n ( k L box ; 2π λ )ψ n ( l L box ; λ) (2.56)
After a simple derivation, we obtain the final result:

lim

N →∞ T k,l = 1 L box +L box /2 -L box /2 e 2πilx/L box e -2πikx/L box dx, (2.57)
which is exactly the Kronecker delta function.

Encoding quality

There are several ways to evaluate the quality of a model encoding with the subsequent reconstruction. For example, in the optimal control theory [START_REF] Stephen P Boyd | Linear controller design: limits of performance[END_REF], the quality of a linear filter is estimated using a certain norm of the transfer matrix. However, in crystallography, the most used quality criterion is the crystallographic R-factor [START_REF] George | X-ray structure determination: a practical guide[END_REF]:

R = l F exact l -F mod l l F exact l , (2.58)
where F exact and F mod are the exact Fourier coefficients of a molecule and the coefficients computed from the Hermite coefficients, respectively. This quantity is a widely used measure of agreement between a crystallographic model and the corresponding experimental X-ray diffraction data. In the case of an ideal electron density encoding, R-factor is equal to zero. In protein crystallography, models with R-factors less than 0.2 are regarded as good when working at a middle resolution. Equations for the transfer matrix allow to estimate the R-factor values for certain classes of electron density distributions. As described above (Eq. 2.6), we use the Gaussian distribution to model the electron density of an atom. Exact Fourier coefficients of a molecule with N atoms atoms at positions r i are then given as:

f exact l,m,n (s) = α 3 π 3 2 Natoms i=1 e -α 2 π 2 s 2 lmn e -2iπr i s lmn , (2.59)
where s l,m,n is the wave vector, s l,m,n = (l/L x , m/L y , n/L z ), with L x , L y , and L z the dimensions of the bounding box along the corresponding axes. Similarly, onedimensional exact Fourier coefficients of the Gaussian function are given as:

f exact l = απ 1 2 Natoms i=1 e -α 2 l 2 /L 2 box e -2iπr i l/L box (2.60)
To see how the Hermite basis encodes Gaussian densities with various level of detail, we built models of electron density map with different parameters α. The width of the Gaussian determines the resolution of the density map according to:

R = πα 2 (2.61)
The derivation of this formula follows the one well known in crystallography, which describes the extinction of diffraction reflexes. For the sake of completeness of the thesis, we provide its derivation in the section 2.3.5.

To estimate R-factor for certain model parameters, we assume that the input electron density is given as a sum of Gaussians with variance of α/ √ 2 equispaced at a distance α. We should stress that due to the properties of the Hermite functions, the whole model is scale-invariant. More precisely, if we keep the product λL box constant, then the relative shape of the Hermite basis functions would not change. Also, if we scale L box and α simultaneously, then the value of R-factor is unchanged. Therefore, it is useful to provide relative resolutions computed as R/L box . Figure 2.3.4 (Left) shows R-factors for the Gaussian parameter α = 0.2 Å, corresponding to the absolute input signal resolution of R = 0.31 Å and the relative resolution of R/L box = 0.014. However, in this case, the actual absolute resolution is cut at L box /M = 0.77 Å, which corresponds to the relative resolution of 0.033. Figure 2.3.4 (Middle) shows R-factors computed using the Gaussian parameter α = 1.0 Å, corresponding to the absolute input signal resolution of R = 1.57 Å and the relative resolution of R/L box = 0.068. Figure 2.3.4 (Right) shows R-factors computed using the Gaussian parameter α = 5.0 Å, corresponding to the absolute input signal resolution of R = 7.85 Å and the relative resolution of R/L box = 0.34. The estimate on the optimal parameter λ (Eq. 2.41) is plotted with the red dashed line. The realspace bound on the optimal parameter λ (Eq. 2.37) is shown with the orange dashed line. The reciprocal-space bound on the optimal parameter λ (Eq. 2.39) is shown with the blue dashed line. We see that lowering the resolution of the input signal, R-factors decrease, as can be expected from general considerations. We can also see that the the lower (Eq. 2.37) and the upper (Eq. 2.39) bounds on the optimal scaling parameter λ follow the isolines of the R-factor map. Therefore, their mean given by Eq. 2.37 provides a reasonable estimation on the optimal value of λ. The input signal is modelled as a sum of Gaussians (Eq. 2.6) with the variance of α/ √ 2 equispaced at a distance α. The number of Fourier coefficients is M = 30, and the input box size is L box = 23.0 Å. These values are chosen to mimic the 1AKG peptide decomposition. The estimate on the optimal parameter λ (Eq. 2.41) is plotted with the red dashed line. The real-space bound on the optimal parameter λ (Eq. 2.37) is shown with the orange dashed line. The reciprocal-space bound on the optimal parameter λ (Eq. 2.39) is shown with the blue dashed line. Left: The Gaussian parameter α = 0.2 Å, corresponding to the absolute input signal resolution of R = 0.31 Å and the relative resolution of R/L box = 0.014. However, in this case, the actual absolute resolution is cut at L box /M = 0.77 Å, which corresponds to the relative resolution of 0.033. Middle: The Gaussian parameter α = 1.0 Å, corresponding to the absolute input signal resolution of R = 1.57 Å and the relative resolution of R/L box = 0.068. Right: The Gaussian parameter α = 5.0 Å, corresponding to the absolute input signal resolution of R = 7.85 Å and the relative resolution of R/L box = 0.34. shape of input electron density and then used Eqs. 2.48 and 2.60 to compute the analytical R-factors. For these plots, we computed the optimal scaling parameter λ using Eq. 2.37. Parameter L box and the size of the transfer matrix M were constant and equal to 23 Å and 30, correspondingly. As in the previous figure, these values are chosen to mimic the alpha-conotoxin PnIB peptide decomposition used in the fitting example below. The scale of the top horizontal axis gives the absolute resolution for L box = 23 Å. The scale of the bottom horizontal axis gives the relative resolution. In order to compute the absolute resolution, its values need to be multiplied by the chosen value of L box . As expected, the values of R-factors diminish as the resolution becomes lower. This is because at low resolutions, lowfrequency columns of the transfer matrix become more important. In the limiting cases of zero and infinite resolutions, R-factor can be computed directly from the transfer matrix as a certain norm of T -I. For the infinite resolution limit, it is given as L 1 norm of the central column of matrix T -I. For the zero resolution limit, R-factor is given by the entry-wise L 1 norm of T -I, R = i,j |T i,j -δ i,j |. Figure 8 also shows an estimation of R-factors for the 3D case. It is based on the assumption that the Hermite decomposition encoding in 3D behaves similar to the 1D case, with the number of coefficients scaled as

N 1D = 3 √ N 3D .

Resolution model

To illustrate the connection between parameter α in the model of electron density (Eq. 2.6) and the resolution of the X-ray diffraction pattern, we use the simplest model. More precisely, we model the electron density as the array of Gaussians in a perfect 1D lattice perpendicular to the incoming radiation beam. Parameter α then plays the role similar to the temperature B-factors. X-ray diffraction intensity depends on the angle between the incoming beam and the direction to the detector θ as:

I ∝ dx f (x) exp 2πix sin θ λ 2 (2.62)
where λ is the wavelength of the incoming radiation. Using the model density (Eq. 2.6), we obtain:

I ∝ α √ πe -(π sin θ λ α) 2 dx ρ(x) exp 2πix sin θ λ 2 , (2.63)
where ρ(x) is the sum of delta functions at the atomic positions. Therefore, the extinction of the diffraction peaks is proportional to e -(π sin θ λ α) 2 2

, where we neglect the quadratic factor before the exponent. According to the definition used in crystallography, resolution is the inter-planar distance in the real space corresponding to the last observable peak in the reciprocal space. Unfortunately, the index of the last peak depends on the detector's noise and strongly depends on the characteristics of the measurement device. Therefore, to give qualitative estimation on the dependence of resolution on the model parameter α, we assume that the last observable peak is the one whose intensity decreases approximately by the factor e 2 . The corresponding angle then reads:

sin θ max = λ πα (2.64)
Therefore, the minimum inter-planar distance, or, the resolution is given by Bragg's law as: R = πα/2 (2.65)

Results and Discussion

We tested and verified our algorithm using two examples of different difficulty. The first example is a small polypeptide alpha-conotoxin PnIB. We generated the EDM for this example from the coordinates of the polypeptide. The second example is the fitting the GroEL domains into the electron density map of the GroEL complex.

Alpha-conotoxin PnIB

First, we explored the relationship between encoding quality and the quality of the fitting. For this purpose, we chose the small 16-residue polypeptide alpha-conotoxin PnIB. We downloaded the X-ray crystal structure of alpha-conotoxin PnIB (PDB code 1AKG) [START_REF] Hu | Crystal Structure at 1.1 A Resolution of alpha-Conotoxin PnIB: Comparison with alpha-Conotoxins PnIA and GI[END_REF] from the PDB database [START_REF] Berman | The protein data bank[END_REF] and simulated the electron density map (2mFo-DFc) using the Uppsala electron density server [START_REF] Kleywegt | The Uppsala electron-density server[END_REF] with the resolution R = 1.1 Å. We computed the protein density according to Eq. 2.6 with the Gaussian width α = 1.0 Å using only the non-hydrogen atoms of the standard amino acids. We rotated the initial 1AKG structure by the arbitrarily chosen Euler angles equal to 76, 234, and 56 degrees, respectively, and used it as the input for the fitting workflow. We used N rot = 500 (corresponding to an angular step of 36°) rotations represented with uniformly distributed Euler angles spanning the space of 2π ×π ×2π. The order of the Hermite expansion was set to N = 15, which is the minimum expansion order allowed at this resolution according to Eq. 2.43. The order of the Fourier expansion was twice the order of the Hermite expansion, M = 30 for each dimension.

To see how the encoding quality influences the fitting algorithm, we studied the dependence of the decomposition on the scaling parameter λ. We chose a range of λ parameters between 0.05 and 2.0. For each λ, we computed the best fitting score along with the average fitting score. Fitting results are shown in Fig. 2.4.1. We see that by choosing λ small, we neglect the details of the protein structure (Fig. 2.4.1 A) and therefore, we can not discriminate between different orientations of the protein (maximum score for λ = 0.05 is very close to the average score). When choosing λ sufficiently large, we obtain satisfactory discriminative power to find the near-native position of the protein (Fig. 2.4.1 C,D). We also see that, e.g., for λ = 0.5, the difference between the maximum and the average score is much larger than in the case of λ = 0.05. Also, when we take λ too large, we can not encode the whole protein (Fig. 2.4.1 E). The red dashed line on Fig. 2.4.1 shows R-factors computed with Eq. 2.58. We see that the choice of parameter λ influences the R-factors and thus determines the quality of the fitting. Notably, the minimum of the R-factor curve corresponds to the maximum of the fitting score.

Due to the strong influence of the scaling parameter λ on the discrimination power of the algorithm, we estimated its optimal value to gain the maximum separation between the score of the correct pose and the average score. Provided that the box that contains all the rotations of the peptide has the size L box = 23 Å and setting the resolution of the EDM R = 1.1 Å, Eq. 2.41 gives an estimate on the optimal value of the scaling parameter λ opt ≈ 0.50. Fig. 2.4.1 shows that this estimation corresponds to the best discrimination between the near-native and all other structures, which can be deduced from the maximum separation between the score of the prediction and the average score. RMSD between the prediction and the solution at this value of λ is 1.03 Å. We should note that the RMSD can be decreased by taking a finer angular search step.

GroEL complex

Here, we demonstrate that our approach obtains essentially the same results as other programs, provided that the scoring function if the same (LCCF in this case). For this purpose, we use a classical test for a fitting algorithm, the GroEL complex map. We downloaded the EDM of the GroEL complex from the Electron Microscopy Data Bank (EMDB), code EMD-2001 with resolution of 8.5 Å. Then, we downloaded the crystal structure of the GroEL subunits from the PDB database. We used the GroEL-GroES complex (PDB code 1AON), from which we extracted the chain A, centered it and arbitrary rotated to exclude any bias. We chose the sampling grid size according to the resolution and the size of the EDM. The EDM was first padded with zeros and then transformed to the Fourier basis using the FFT algorithm. The number of coefficients in the Fourier decomposition M was equal to 105 × 107 × 119. The angular search step was set to 30°. We used the Hermite expansion order of N = 15, which is larger than the minimum expansion order allowed at this resolution, N min ≈ 9 (see Eq. 2.43). We sampled the rotations using the spiral algorithm [START_REF] Edward | Distributing many points on a sphere[END_REF], which generates an equispaced distribution of points on a sphere. Unlike in the previous example, due to the lower resolution of the GroEL EDM, here we fitted Laplacian filtered protein density into the Laplacian filtered EDM.

After the 6D exhaustive search, we clustered the solutions using the clustering threshold of 10 Å and kept the top 14 poses. All the 14 poses corresponded to the individual chains of the complex, which comprises 2 heptameric rings structure. Fig. 2.4.2 shows the result of the fitting. We compared the fitted model with the model provided by the authors of the EDM (PDB entry code 4AAU). The average RMSD between the chains due to flexible deformations measured using Cα atoms was 3.0 Å. More precisely, we super-posed the corresponding chains of both models using rigidbody transformations and then measured RMSD between them. Overall, the average RMSD between Cα atoms was 5.35 Å. This includes both the discrepancy between corresponding chains in the assembly due to flexible deformations and because of the rigid body misfit. The average distance between the centers of mass of the corresponding chains was 2.64 Å (Table 2 ). For each pair of models, RMSD was measured using the C α -atoms and the centers of mass of the corresponding chains and then averaged over all chains comprising the assembly.

Runtime of Hermite-to Fourier-space transition

The use of the fast Fourier transform was the inevitable step in every fitting algorithm up until now. Instead, we introduced the basis from which we can transform a decomposition into the Fourier basis avoiding evaluation of the FFT on a grid. When the grid becomes large, the asymptotic complexity of our algorithm becomes O (M 3 N ) (see Eq. 2.32 ). It is comparable to the complexity of the fast Fourier transform algorithm, O (M 3 log M ). Intuitively, at large orders of the Fourier expansion M , our algorithm should be faster compared to the FFT. However, prefactors preceded the complexities of the two algorithms are different. Thus, we conducted a numerical experiment to compare the actual running times. Fig. 2.4.3 shows the time needed to compute the FFT on a cubic grid of size M and the time needed to transform a Hermite expansion of order N = 15 to the same Fourier grid. We can see that, generally, at large values of M , M 100, the transition from the Hermite into the Fourier space is faster compared to the speed of the FFT. Also, the timing of the transition grows evenly with respect to M in contrast to the timing of the FFT. One has to take into account that we compared our algorithm with the highly optimized FFTW3 library [START_REF] Frigo | The design and implementation of FFTW3[END_REF]. Probably, additional optimization of HermiteFit could improve performance even further. One of the ways to speed up the transition will be to use the Fast Hermite Transform instead of the naive matrix multiplication [START_REF] Leibon | A fast Hermite transform[END_REF]. This implementation will be the subject of our future work.

Comparison with Situs and ADP EM

We compared the HermiteFit algorithm with two popular existing fitting methods, the colores program from the Situs package [START_REF] Chacón | Multi-resolution contour-based fitting of macromolecular structures[END_REF] and the ADP EM fitting tool [START_REF] José | ADP EM: fast exhaustive multi-resolution docking for high-throughput coverage[END_REF]. These two packages represent the two major approaches to the problem of exhaustive search in the six-dimensional space of rigid-body motions. Colores, a widely used CCF-based fitting tool, rapidly scans the translational degrees of freedom using the fast Fourier transform. The rotations, though, are sampled exhaustively by enumerating a list of equispaced distributed rotations on a sphere. ADP EM choses points in real space, places there the atomic structure and then rotationally matches it to the EDM using the Fast Rotational Matching algorithm. The authors of the ADP EM compared their algorithm with the 5D rotational matching and found that the 3D rotational matching works faster in practice [START_REF] José | ADP EM: fast exhaustive multi-resolution docking for high-throughput coverage[END_REF].

For the comparison, we normalized the running time of the fitting algorithms by the sizes of the search space. For colores and HermiteFit, the size of the search space is equal to the number of grid cells (M 3 for a cubic grid in the HermiteFit algorithm) multiplied by the number of sampled angles. The size of the search space of the ADP EM algorithm is the number of points in real space times the number of cells of the angular grid. The latter is built from uniformly sampled Euler angles on a grid of 2π × π × 2π. The size of the angular grid is determined by the order N exp of a spherical harmonics expansion and equals to 4N 3 exp . For colores and HermiteFit, we used the angular search step of 30 . The resolution of the EDM for colores and HermiteFit was set to 8.5 Å. The Fourier grid that was used by colores and the HermiteFit algorithm had dimensions 105 × 107 × 119. For ADP EM, we used the spherical harmonics expansion order of N exp = 16.

Figure 2.13: Running times of the Hermite to Fourier space transition performed using our algorithm and the FFT algorithm on a cubic grid of M × M × M as a function of the Fourier expansion order M . We used the FFTW3 library [START_REF] Frigo | The design and implementation of FFTW3[END_REF] with the double precision real discrete Fourier transform using the flag FFTW ESTIMATE to measure the speed of the FFT. The order of the Hermite expansion was N = 15.

Table 2.3 shows the normalized times of the complete 6D search for the three algorithms in the case of fitting the GroEL subunit into the 8.5 Å GroEL electron density map. Judging by the total running time, ADP EM has a big advantage over the two other algorithms, which exhaustively search all the space of possible translations. However, in terms of running time per one search point, the Hermit-eFit algorithm is more effective than the other two. Interestingly, colores spends about half of the total search time on the computation of the Fourier coefficients of the rotated protein. Therefore, it was very important for us to speed up this step. Nonetheless, all three tested algorithms have their own advantages and drawbacks. For example, ADP EM can use smart heuristics to contract the number of search points in the real space. However, its sample points in the space of rigid body rotations are distributed non-uniformly. In particular, near the poles rotations are sampled more densely, making this sampling scheme less effective [START_REF] Edward | Distributing many points on a sphere[END_REF]. On the other hand, the HermiteFit algorithm along with the colores algorithm sample the rotational space nearly uniformly using the spiral algorithm while the translational space sampling also remains uniform. We would like to stress that the absolute runtimes (shown in Table 2.3) are not very informative. In particular, they dramatically depend on the choice of the FFT library, code optimization, the choice of compiler and compilation options, etc. However, this comparison clearly demonstrates that the new approach paves the way to speed up one of the bottlenecks of fitting methods, the projection of the rotated structure into the Fourier space.

To assess the fitting quality of the tested methods, we measured the RMSDs between the obtained models and the structure obtained by the authors of the electron density map (PDB entry 4AAU). Table 2.2 shows the comparison of the measured RMSDs for ADP EM, colores and HermiteFit. We used two different criteria for the measurements. First, we measured the average RMSD between α-carbons. Second, we measured the average distance between the centers of mass of the corresponding chains. ADP EM produced a model with RMSD of 4.61 Å from the solution, RMSDs for colores and HermiteFit were 5.42 Å and 5.35 Å, respectively. Clearly, Chapter 3

Scoring functions for protein-protein docking

Introduction

As was described in the introduction, scoring method is used to filter out falsepositive predictions of rigid-body docking step and refine those close to the native conformation of a protein-protein complex. Therefore the method of scoring has a decisive role in success of the whole docking workflow.

In this section we propose a new method to derive scoring functions. We base our method on separation of the native structures and computationally generated nonnative conformations of a complex (decoys). Most of previously used algorithms solving this problem separate all the decoys from all the decoys simultaneously.

The key new idea behind our algorithm is that the decoys of a particular complex should be separated from its native structure only. However the form of the scoring potentials should be the same for all the dataset. Based on these prepositions we show that this problem leads to the well-defined convex quadratic optimization problem. We measure the performance of the scoring functions obtained on the commonly used benchmarks. We show that our algorithm has inherent stability against overfitting. Also, due to the properties of the basis we used our scoring potentials show some interesting coarse-graining properties.

Given the widely recognized importance of water molecules at the protein-protein interface we also developed potentials for the prediction of water molecules relying on the general ideology of the scoring potentials we developed for the protein-protein contacts.

Methods

Problem Formulation

Consider N native proteins configurations P nat i , i = 1...N . For each protein complex number i we generate D decoys, P nonnat ij , j = 1...D, where the first index runs over different protein complexes and the second index runs over decoys. Our goal is to find a scoring functional F , defined for all possible protein-protein complex structures (the set P), such that for each native complex i and its nonnative decoy j the following inequality holds:

F (P nat i ) < F (P nonnat ij ) (3.1)
There are many ways to construct the functional F (Eq. 3.1). To outline its form, we are relying on the following assumptions:

1. Functional F depends only on the interface between the proteins. We define the interface as a set of all atom pairs at a distance smaller than a certain cutoff distance r max , such that the first atom in each pair belongs to the first protein and the second atom in each pair belongs to the second protein.

2. The protein is represented as a set of discrete interaction sites that are located at the centers of the atomic nuclei. All interaction sites are divided into M types according to the properties of the corresponding atomic nuclei. In this study we choose M = 20.

3. Functional F depends only on the distribution of the distances between the interaction sites (the number of site pairs at a certain distance), F (P ) = F (n 11 (r), .., n kl (r), .., n mm (r)) = F (n(r)), (3.2)

where n kl (r) is the number density of site-site pairs separated by a distance r, with site k located on the first protein, and site l located on the second protein. For homogeneous systems, such as liquids, functions n kl (r) can be expressed via site-site radial distribution functions g kl (r), which can be obtained experimentally, as n kl (r) = 4πr 2 ρg kl (r)N a , where ρ is the number density and N a is the total number of atoms in the system [START_REF] Hansen | Theory of simple liquids[END_REF]. However, for proteins this is not the case.

4. F is a linear functional, F (αn

1 (r) + βn 2 (r)) = αF (n 1 (r)) + βF (n 2 (r)).
One of the simplest functionals F (n(r)) fulfilling these assumptions can be written as:

F (n(r)) ≡ F (n 11 (r), .., n kl (r), .., n M M (r)) = M k=1 M l=k rmax 0 n kl (r)U kl (r) dr (3.3)
It contains unknown functions U kl (r) that can be determined from the training set of native protein complexes. From now on, we will call these functions scoring potentials. 1 Once the scoring functions are known, to compute the value of F we need to specify site-site number densities n kl (r). In practice, we calculate them as a sum of all k -l distances in a given protein complex using the equation:

n kl (r) = ij 1 √ 2πσ 2 e -(r-r ij ) 2 2σ 2 , (3.4)
where each distance distribution is represented with a Gaussian centered at r ij with the variance of σ2 . The sum is taken over all k -l site pairs i and j separated by the distance r ij smaller than r max , with site k located on the first protein of the complex, and site l located on the second protein. In the limiting case of the variance tending to zero, Eq. 3.4 turns into a sum over Dirac delta functions. In our study we assume the value of σ to be fixed for all site-site distributions. However, if one has an additional information about individual distance distributions, e.g. Debye-Waller factors, molecular dynamics trajectories, etc., it can be used for more precise parametrization of the variance or even instead of the Gaussian approximation in Eq. 3.4. Finally, we compute the score of each conformation using equation 2 :

Score = ij Υ kl (r ij ) (3.5)
where the sum is taken over all pairs of atoms i and j separated by the distance r ij smaller then r max , with atom i of type k located on the first protein of the complex, and atom j of type l located on the second protein. The function Υ kl (r) is the Gauss transform of the scoring potential U kl (x):

Υ kl (r) = 1 √ 2πσ 2 rmax 0 e -(x-r) 2 2σ 2 U kl (x) dx (3.6)

Expansion of U (r) and n(r) in an orthogonal basis

Given a set of functions ξ p (r) orthogonal on the interval [r 1 ; r 2 ] with a nonnegative weight function Ω(r) such that

r 2 r 1 ξ p 1 (r)ξ p 2 (r)Ω(r) dr = δ p 1 p 2 , (3.7)
where δ p 1 p 2 is the Kronecker delta function, scoring potentials U kl (r) and number densities n kl (r) can be expanded on the interval [r 1 ; r 2 ] as:

U kl (r) = p w kl p ξ p (r) Ω(r), r ∈ [r 1 ; r 2 ] (3.8) n kl (r) = p x kl p ξ p (r) Ω(r), r ∈ [r 1 ; r 2 ] (3.9)
Expansion coefficients w kl p and x kl p can be determined from the orthogonality condition (Eq. 3.7) as

w kl p = r 2 r 1
U kl (r)ξ p (r) Ω(r) dr (3.10)

x kl p = r 2 r 1 n kl (r)ξ p (r) Ω(r) dr (3.11)
Using expansions Eq. 3.8 and Eq. 3.9, the functional F (n(r)) can be rewritten as:

F (n(r)) = M k,l rmax 0 p 1 p 2 w kl p 1 x kl p 2 ξ p 1 (r)ξ p 2 (r)Ω(r) dr (3.12)
In this study we use two types of functions ξ p (r) orthogonal on the interval [0; 10] with a unit weight, (i) shifted Legendre polynomials and (ii) traditionally used shifted rectangular functions. These two types of functions are plotted in Fig. 3.1. Other types of orthogonal functions can also be used. If the functions ξ p (r) are chosen to be negligibly small outside the interval [0; r max ] or if their interval of orthogonality [r 1 ; r 2 ] coincides with the interval [0; r max ], as is the case for two sets of our functions, then the scoring functional F (n(r)) can be expanded up to the order P as:

F (n(r)) ≈ M k=1 M l=k P p w kl p x kl p = (w • x), w, x ∈ R P ×M ×(M +1)/2 (3.13)
We will refer to the vector w as to the scoring vector and to the vector x as to the structure vector. Equations 3.4 and 3.11 provide the projection from a protein complex structure into the scoring space R P ×M ×(M +1)/2 . Using these formulas, we can project structural information of each protein complex into a certain structure vector x on R P ×M ×(M +1)/2 . constructed for the case with multiple complexes. In case of two classes of vectors, Vapnik proposed to use a special kind of separator, the so-called optimal separating hyperplane [START_REF] Vapnik | The nature of statistical learning theory[END_REF], which is unique and maximizes the distance to the closest point from either class. We can generalize this idea and formulate the following quadratic programming optimization problem:

Minimize (in w, b j ) 1 2 w • w Subject to y ij [w • x ij -b j ] -1 ≥ 0, i = 1...N, j = 1...D (3.16)
where y ij = -1 when the structure vector x ij is native and y ij = 1 otherwise. Now we are ready to formulate Lemma 1. If exists such a linear scoring functional of form of Eq. 3.13 that correctly discriminates the native structure vectors for all complexes (Eq. 3.15), then the optimal scoring vector is unique and given by the solution of problem (Eq. 3.16).

Remark. The scoring vector is optimal in the sense that it maximizes the separation between native and nonnative structure vectors.

Generally, such a linear scoring functional (with a fixed value of the expansion order P ) may not exist, as demonstrated below. Therefore, we will have to modify the optimization problem (Eq. 3.16).

Case II. No solution exists

In Fig. 3.3B we present an example when no hyperplane can separate the two classes of the structure vectors of a single complex. For this case, Cortes and Vapnik proposed to relax the condition for the optimal separating hyperplane [START_REF] Cortes | Support-vector networks[END_REF], including an additional term. This term minimizes the sum of penalties for misclassified vectors. We again generalize this idea and introduce for each decoy set j = 1...D slack variables ξ ij , which are positive for misclassified structure vectors and zero otherwise. A non-zero value of ξ ij allows the structure vector x ij to overcome the inequality condition in Eq. 3.16 at a cost proportional to the value of ξ ij (see Fig.

3.3B). The new soft-margin quadratic optimization problem reads:

Minimize (in w, b j , ξ ij ):

1 2 w • w + ij C ij ξ ij Subject to: y ij [w • x ij -b j ] -1 + ξ ij ≥ 0, i = 1...N, j = 1...D ξ ij ≥ 0 (3.17)
The solution of this problem provides a trade-off between how large will be the separation between two classes of the structure vectors of each complex and how many misclassified vectors will be in the solution. Parameters C ij can be regarded as regularization parameters. The solution of Eq. 3.17 tends to maximize the structure vector separation for small values of C ij or minimize the number of misclassified structure vectors for large values of C ij . We choose parameters C ij to be different for native and nonnative structure vectors of each complex because fewer native structure vectors should have the larger weight (see for instance [START_REF] Akbani | Applying support vector machines to imbalanced datasets[END_REF]). The following lemma provides the foundation for the numerical scheme used in this work: Lemma 2. The optimal scoring vector is unique and given by the solution of problem (Eq. 3.17).

Remark. Here, the scoring vector is optimal in the sense that it maximizes the separation between native and nonnative structure vectors and minimizes the number of misclassified vectors. Regularization parameters C ij in Eq. 3.17 tune the importance of either factors.

The proof of lemmas (1, 2) can be found, e.g., in [START_REF] Cristopher | Uniqueness of the SVM Solution[END_REF]. Overall, the formulation of the optimization problem (Eq. 3.17) is very similar to the formulation of soft-margin support vector machine (SVM) problem [START_REF] Cortes | Support-vector networks[END_REF]. Therefore, to solve this problem (Eq. 3.17) we will use techniques developed for SVM. A) The case when infinitely many hyperplanes can separate the two classes. B) The case when no optimal separating hyperplane exists. Slack variables ξ i and ξ j for misclassified structure vectors are added, which are the distances to the corresponding margin hyperplanes. The optimal hyperplane, which maximizes the separation between the two classes, is plotted as a dashed line. Two margin hyperplanes are plotted as solid lines.

Algorithm

Dual problem

Properties and solutions of quadratic optimization problems similar to the one stated above (Eq. 3.17) have been extensively studied in the theory of convex programming [START_REF] Boyd | Convex optimization[END_REF][START_REF] Vapnik | The nature of statistical learning theory[END_REF]. For instance, using the Lagrangian formalism , the optimization problem (Eq. 3.17) can be converted into its dual form (Appendix A), and the resulting dual optimization problem is convex :

Maximize L(λ ij ): L(λ ij ) = ij λ ij -1 2 ij kl y ij y kl λ ij λ kl x ij • x kl Subject to: 0 ≤ λ ij ≤ C ij i y ij λ ij = 0, ∀j , (3.18)
where the maximization is performed with respect to the Lagrange multipliers λ ij . This dual problem is similar to the the soft-margin SVM optimization problem [START_REF] Cortes | Support-vector networks[END_REF]. The difference lies in the constraints. For the soft margin SVM, conditions on the parameters written in the same two-indexed form as in Eq. 3.18, are ij y ij λ ij = 0.

Vectors x ij for which λ ij > 0 are called support vectors. Once the dual problem (Eq. 3.18) is solved and the Lagrange multipliers λ ij are found, we can express the solution of the original primal problem (Eq. 3.17) (the scoring vector) as a linear combination of the support vectors:

w = support vectors y ij λ ij x ij (3.19)
The dual representation (Eq. 3.18) of the original primal quadratic problem (Eq. 3.17) allows us to break down the original large quadratic optimization problem into a series of smaller sub-problems. Below we describe an algorithm that solves the dual optimization problem (Eq. 3.18) using a decomposition technique.

Block sequential minimal optimization algorithm

Due to its enourmus size, the quadratic optimization problem can not easily be solved by standard techniques. The quadratic form in Eq. 3.18 involves a matrix with number of elements proportional to the squared number of the training structure vectors. This matrix often exceeds the size of available RAM, for instance, explicit storage of the matrix used in the current study requires about 20GB of memory. Nonetheless, algorithms that deal with large datasets are widely used in machine learning. More precisely, various decomposition techniques have been developed to reduce the requirements of optimization solvers to the size of available RAM [START_REF] Vapnik | Estimation of dependences based on empirical data[END_REF][START_REF] Platt | Sequential minimal optimization: A fast algorithm for training support vector machines[END_REF]. Here, we employ a block-decomposition technique and propose the block sequential minimal optimization (BSMO) algorithm. Briefly, we partition the training set into N blocks, each block containing one native structure vector with its D nonnative structure vectors. Then, for each block i, we iteratively optimize each pair of Lagrange multipliers (λ 1 , λ 2 ), preserving the equality constraint y 1 λ 1 + y 2 λ 2 = const. To do this, we write the Lagrangian (Eq. 3.18) as a function of λ 1 and λ 2 :

L(λ 1 , λ 2 ) = 1 2 ηλ 2 2 -ηλ 2 λ old 2 + λ 2 y 2 (y 2 -y 1 ) + λ 2 y 2 (x i1 -x i2 ) • w old + Const.(3.20) with η = 2x i1 • x i2 -x i1 • x i1 -x i2 • x i2 (3.21)
Then, we analytically maximize this Lagrangian with respect to λ 1 and λ 2 according to the sequential minimal optimization (SMO) algorithm [START_REF] Platt | Sequential minimal optimization: A fast algorithm for training support vector machines[END_REF]. After the minimization, we obtain new values of λ 1 and λ 2 . We provide more details about the SMO algorithm in Appendix B. After each iteration, we recompute the current scoring vector w new (see Eq. 3.19) according to:

w new = w old + ∆λ 1 y 1 x i1 + ∆λ 2 y 2 x i2 (3.22)
For each block i, we continue the iterative optimization of the Lagrangian (Eq. 3.20) until the relative change in its value between two successive inner cycles of iterations is less than the desired tolerance. Each inner cycle consists in the optimization of all pairs of Lagrange multipliers for a given block i. Globally, we terminate the optimization when the relative change in the value of the Lagrangian (Eq. 3.18) between two successive outer cycles is less than the desired tolerance. Each outer cycle consists in the optimization of all blocks of the training set. The flowchart of the BSMO algorithm is presented in Figure 3.2.4.2.

Load block i

Compute:

q ij =1-w⋅x ij C ij h =C ij -α ij C ij l =α ij
Solve the problem using SMO Output: As it is seen from Eq. 3.21, our BSMO algorithm requires only scalar products of the structure vectors within the same block. Therefore, it is sufficient to load each block into RAM sequentially, which results in memory efficiency of our method. Precisely, RAM required for our implementation of the block-decomposition solver is N 2 times less compared to the direct quadratic problem solvers.

δ L> δ L max δ L max =δ L For all i δ L max > ϵ no yes no exit δ L max =0 yes Output: w , α ij d j , δ L w :=w+ ∑ y ij d j x ij α ij :=α ij + d j

Training database for protein-protein interactions

In the present study we used the training database of 851 non-redundant proteinprotein complex structures prepared by Huang and Zou [START_REF] Huang | An iterative knowledge-based scoring function for protein-protein recognition[END_REF]. This database contains protein-protein complexes extracted from the PDB [START_REF] Berman | The protein data bank[END_REF] and includes 655 homodimers and 196 heterodimers. We updated three PDB structures from the original training database: 2Q33 supersedes 1N98, 2ZOY supersedes 1V7B, and 3KKJ supersedes 1YVV. The training database contains only crystal dimeric structures determined by X-ray crystallography at resolution better than 2.5 Å. Each chain of the dimeric structure has at least 10 amino acids, and the number of interacting residue pairs (as defined as having at least 1 heavy atom within 4.5 Å) is at least 30. Each protein-protein interface consists only of 20 standard amino acids. No homologous complexes were included in the training database. Two protein complexes were regarded as homologues if the sequence identity between receptor-receptor pairs and between ligand-ligand pairs was > 70%. Finally, Huang and Zou [START_REF] Huang | An iterative knowledge-based scoring function for protein-protein recognition[END_REF] manually inspected the training database and left only those structures that had no artifacts of crystallization.

Our algorithm requires as input native and nonnative structure vectors (see, e.g., Eq. 3.15). Native structure vectors can be computed from the native protein-protein contacts in the training database using Eq. 3.11. However, for the computation of the nonnative structure vectors for each protein-protein complex from the training database, we need to generate decoys for each complex. Since our optimization algorithm is very general and has no special requirements for nonnative protein-protein contacts, we generated them by "rolling" a smaller protein (ligand) over the surface of a bigger protein (receptor) using HEX protein docking software [START_REF] Ritchie | Protein docking using spherical polar Fourier correlations[END_REF]2]. To do so, we initialized HEX exhaustive search algorithm with the radial search step of 1.5 Å and expansion order of the shape function equal to 31. We used only the shape complementarity energy function from HEX (i.e., electrostatic contribution was omitted). Afterwards, we clustered HEX docking results with a root mean square (RMS) threshold of 8 Å. The top 200 clusters, ranked by HEX surface complementarity function, plus the native protein-protein complex conformation (giving in total 201 structures) were then used to compute the distance distribution functions (Eq. 3.4). Then, we computed the structure vectors using Eq. 3.11 and labeled them as "native" if the root mean square deviation (RMSD) of the corresponding ligand was < 2 Å from its native position. Otherwise, the structure vector was labeled as "nonnative" or "decoy". On average, we obtained about 2.5 native structure vectors (and, correspondingly, about 198.5 nonnative structure vectors) per protein-protein complex. To each structure vector x ij we assigned a regularization parameter C ij according to

C native ij = CD nonnative j /D j C nonnative ij = CD native j /D j (3.23)
We repeated the same procedure for each protein-protein complex from the training database. We used M = 20 atom-centered interaction sites based on the atom types definitions provided by Huang and Zou [START_REF] Huang | An iterative knowledge-based scoring function for protein-protein recognition[END_REF]. These atom types were defined by the classification of all heavy atoms in 20 standard amino acids according to their element symbol, aromaticity, hybridization, and polarity. These 20 atom types result in total of 210 pair potentials.

Our training set has several proteins homologous to the ones from the two widely used docking benchmarks, Rosetta, and Zdock, which we employ below to validate our results. We define two protein complexes to be homologous if for each chain in the first complex there is a chain in the second complex with sequence identity more than 60%. We determined the sequence identity using the FASTA36 program [START_REF] Pearson | Comparison of DNA sequences with protein sequences[END_REF]. Below, we benchmarked our scoring function while both excluding homologs from the training set and leaving it unchanged. The comparison of the benchmark results in these two cases is shown in supplementary Tables S1 andS2.

Training database for water-protein interactions

To obtain potentials for the water molecules we employed the same set of complexes prepared by Huang and Zou [START_REF] Huang | An iterative knowledge-based scoring function for protein-protein recognition[END_REF]. We added a new atom type to the set of 20 atom types for docking corresponding to the water oxygen. We generated a cube with dimensions 130 × 130 × 130 Å of TIP3P water molecules and equlibrated it using molecular dynamics at room temperature with periodic boundary conditions. Afterwards native structures and decoys from the training set were immersed into this water-box. In the case when the size of the protein was more than 130 Å we periodically continued the water box. Then, in case of native structures initial water molecules presented in the crystallographic structures were retained. In the case of decoys we removed all the original water molecules. The placed bulk molecules were then removed if they clashed with either the proteins or retained water molecules. The clashing distance was set to be 2 Å. We also removed each water molecules farther than 12 Å from all protein atoms. Although we aimed to obtain all pairwise potentials we did not use water-water potential for the predictions, because it influences the equilibrium structure of bulk water, thus in principle after obtaining this potential we had to recalculate water box. However this procedure goes beyond our method of obtaining scoring potentials.

Results and Discussion

Overfitting and Convergence

Various methods of derivation of the knowledge-based potentials usually produce results biased towards the training data set. Typically, such algorithms maximize the predictive accuracy of the corresponding potential on a set of training data, which does not mean that the same potential will perform equally well on a new set of data. Indeed, fitting the potential to the training data set also fits the noise in the data. Thus, very often a knowledge-based potential memorizes noisy features of the training data instead of deducing general predictive concepts from it. This phenomenon is usually referred to as overfitting [START_REF] Dietterich | Overfitting and undercomputing in machine learning[END_REF]. A clear indication of an "overfitted" potential can be, for example, the need for post-smoothing techniques applied to the initial knowledge-based potential, as in [START_REF] Huang | An iterative knowledge-based scoring function for protein-protein recognition[END_REF][START_REF] John | BLEEP-potential of mean force describing protein-ligand interactions: I. Generating potential[END_REF]. "Overfitting" is clearly not desirable. In order to avoid it, many optimization techniques (regularization, crossvalidation, etc.) have been successfully proposed to penalize the initial objective function with various additional terms [START_REF] Arlot | A survey of cross-validation procedures for model selection[END_REF][START_REF] Kearns | An experimental and theoretical comparison of model selection methods[END_REF]. These terms serve to achieve a better predictive accuracy on the off-training data based on the predictions of the training data.

To avoid overfitting, we used two techniques -regularization and cross-validation. Regularization penalizes the initial objective function with various additional terms [START_REF] Arlot | A survey of cross-validation procedures for model selection[END_REF][START_REF] Kearns | An experimental and theoretical comparison of model selection methods[END_REF]. We introduced two regularization parameters, σ for the width of the Gaussian distribution of distances in Eq. 3.4, and C for the hinge loss function in Eq. 3.17. To find the best values of these parameters we used the following cross-validation procedure. First, we divided the training set into two parts, consisting of 650 complexes (temporary test set) and 200 complexes (temporary training set). Then, for each value of σ and C, we obtained the scoring potentials using the temporary training set and verified it on the temporary test set. Finally, we chose those values of σ and C that correspond to the maximum number of guessed structures in the temporary test set. We define the structure as guessed if its native complex has the score better than all of its decoys. Figure 3.5 shows the predictive performance of the scoring potential on the two sets as a function of σ and C. Obviously, the maximum predictive performance on the training set is achieved at the highest values of C. However, the validation on the test set highlights the best choice of values of C and σ. These values are C = 10 6 . . . 10 7 and σ = 0.4 Å (Fig. 3.5B).

Figure 3.6 shows the convergence of the success rate on the training set with the number of iterations of the training BSMO algorithm. The success rate was measured as the number of guessed structures divided by the total number of proteinprotein complexes. We can see a fast convergence of the method. In principle, a hundred optimization steps is sufficient to obtain the final result. We have also observed that increasing the regularization parameter C leads to the slower convergence and vice versa. We should note that thanks to the convexity of our optimization problem, its solution is unique and does not depend on the starting point and the optimization method used.

Extracted Potentials

Our method can in principle use any type of orthogonal polynomials to decompose the structural statistics and reconstruct the potentials. However, since rectangular functions are the most widely used to collect statistics, we employed this basis as a reference. Then, we also used the Legendre basis, orthogonal on the interval [0;10]. We chose this basis because of its simplicity, in particular because the weight function for this basis does not depend on the distance.

From now on, we call the obtained scoring potentials the Convex Protein Protein (ConvexPP) potentials. Figure 3.7 shows typical scoring potentials derived using two different orthogonal bases. Obtained potentials are smooth by construction, thanks to the smooth Gaussian kernel in Eq. 3.6. According to the plot, the shape of the potentials does not depend on the basis set that was used to derive it. This is the consequence of the global convergence of the optimization problem (see lemma 2). We can also see that the obtained potentials tend to zero as the interaction distance increases. On the other hand, all the potentials approach zero at short distances. This is due to the absence of statistics for the native structures at short distances and the result of the w • w term in optimization problem 3.17. We discuss this behaviour in more detail below.

Due to the Gaussian smoothing of statistics, it is sufficient to use the expansion order of P = r max /σ. For σ = 0.4 Å and r max = 10 Å, the estimate on the number of basis functions is P = 25. However, due to the adjustment of σ with the cross-validation procedure, in our experiments we used a larger expansion order, P = 40. Figure 3.8 demonstrates how the resulting potentials depend on the expansion order. We should note that the decompositions of orders above 25 are almost indistinguishable and thus are not shown. Indeed, increasing the order of the polynomial P decreases the distance between two consecutive zeros in the polynomial basis. Therefore, at a constant value of σ, the integral in Eq. (3.11) tends to zero with the growth of the value of P due to the oscillatory behaviour of the basis polynomials. Such behaviour of the integral confines all the useful information about the distributions and the scoring potentials in the first few coefficients of the polynomial decomposition. The number of these coefficients depends solely on the value of σ and does not change with the training set or the value of parameter C.

Protein-Protein docking benchmark version 3.0

First we tested the ConvexPP scoring function on the protein-protein docking benchmark version 3.0. It consists of 124 crystallographic structures of protein-protein complexes from PDB database [START_REF] Hwang | Protein-protein docking benchmark version 3.0[END_REF]. They are divided into three groups: rigid, medium and difficult cases. The division criteria is the scale of conformational changes of the proteins upon binding: from minor changes in rigid cases to major in difficult. The nonredundancy of the benchmark was set at the level of family-family pairs. That means that if a complex in Benchmark v3.0 is formed of the protein of family A and another one of family B, then there are no more family A -family B complex in the benchmark. The assignment of a protein to a family was taken according to SCOP database [START_REF] Alexey G Murzin | SCOP: a structural classification of proteins database for the investigation of sequences and structures[END_REF].

The decoys for scoring were generated using ZDOCK3.0 [START_REF] Brian G Pierce | Accelerating protein docking in ZDOCK using an advanced 3D convolution library[END_REF] with the sampling step equal to 6 degrees (we call this set of docking position ZDOCK benchmark below). They were downloaded from the ZLAB website [3]. The docking program ZDOCK3.0 generates the rigid-body protein-protein docking predictions with the corresponding scores. Scoring function used in this program includes shape complementarity, statistical pair potentials and electrostatics. To compare our scoring function to the well established one downloaded the decoy set reranked by ZRANK [START_REF] Pierce | ZRANK: reranking protein docking predictions with an optimized energy function[END_REF]. It is the program for reranking the ZDOCK3.0 predictions. In addition to the factors used in ZDOCK3.0, it computes detailed electrostatics, estimates desolvation and uses additional Van-der-Waals potential to re-score the decoys.

The benchmark 3.0 has several complexes homologous to certain protein complexes in the training set. Therefore, to see the effect of training set and test set similarity we trained our potential both excluding homologs from the training set and leaving it unchanged.

Rosetta Benchmark

Baker, Gray et al generated the Rosetta benchmark using 54 complexes of the ZDOCK Benchmark 0.0 [START_REF] Chen | Docking unbound proteins using shape complementarity, desolvation, and electrostatics[END_REF] using a flexible docking protocol, which is a part of the RosettaDock suite [START_REF] Gray | Protein-protein docking with simultaneous optimization of rigid-body displacement and side-chain conformations[END_REF]. The first step in the protocol is the random translation and rotation of one of the proteins constituting the complex. Afterwards, the side chain is optimized simultaneously with the rigid body displacement. Finally, the full-atom minimization is done to refine the conformation. For each complex, Baker and Gray generated 1000 decoys following the described protocol. The resulting decoys with the corresponding scores assigned by the RosettaDock program can be obtained from [1]. We calculated the success rate of RosettaDock using the same quality criteria as in Critical Assessment of PRediction of Interactions [START_REF] Méndez | Assessment of blind predictions of protein-protein interactions: current status of docking methods[END_REF][START_REF] Méndez | Assessment of CAPRI predictions in rounds 3-5 shows progress in docking procedures[END_REF][START_REF] Janin | Assessing predictions of protein-protein interaction: the CAPRI experiment[END_REF]. The Rosetta benchmark contains 5 complexes homologous to the ones present in the training set. Therefore we trained our scoring function using training sets with and without these homologs (Table S2). Table 3.2 shows that our potentials significantly improve Top1 prediction rate over ITScore-PP and RosettaDock scoring functions while also outperforming them according to the other criteria (Top1 and quality 1 etc.). We computed the percentage of the structures for which the first acceptable prediction was ranked within the top predictions for each complex and plotted it on Fig. 3.10. According to the plot our scoring function outputs the plausible structure (quality ≥3) for more complexes than ITScore-PP and RosettaDock. Unlike the results on the ZDock benchmark, the results on the Rosetta unbound benchmark slightly decrease when we remove homologous complexes from the training set. Among the prediction quality criteria it is the number of predicted high quality structures that changed the most. On the other hand Top1 prediction rate stayed almost the same. This observation signalizes that the number of predicted high-quality structures is amenable to overfitting. Therefore unlike the Top1 criterion, it can not serve as a reliable measure of a scoring function predictive power. Table C.2 shows the per-complex comparison of two scoring functions trained with and without homologs. The percentage of the structures with plausible predictions (quality higher than 3) that are ranked within a certain number of top predictions. The data for ITScore-PP was taken from original publication [START_REF] Huang | An iterative knowledge-based scoring function for protein-protein recognition[END_REF]. Scoring data for RosettaDock was taken from [1].

Water interaction potentials prediction for CAPRI T47 target

We used the scaled Legendre polynomials as the basis to obtain the potentials for water oxygen interactions. Obtained solvation scoring functions for three atom types (aromatic carbon, carboxyl oxygen, and charged nitrogen) are shown in Figure 3.11. There, one can clearly see the difference between hydrophobic and hydrophilic interactions. Another property of our solvation scoring functions is a same peculiarity at short distances as in the other potentials. Precisely, as there is no training data for distances from 0 to about 2 Å, all obtained scoring functions are close to zero in this range. Therefore, if one would like to use these function for minimization, one might need to adjust them at short distances with some additional information. The Critical Assessment of Predicted Interaction (CAPRI) [START_REF] Méndez | Assessment of blind predictions of protein-protein interactions: current status of docking methods[END_REF][START_REF] Méndez | Assessment of CAPRI predictions in rounds 3-5 shows progress in docking procedures[END_REF][START_REF] Janin | Assessing predictions of protein-protein interaction: the CAPRI experiment[END_REF] is the community-wide competition in prediction of protein-protein complexes. To validate our poteintials for water prediction we took part in predicting the structure of the target 47 [START_REF] Marc F Lensink | Blind prediction of interfacial water positions in CAPRI[END_REF]. The participants were asked to predict the interface of the two proteins: DNase domain of colicin E2 and IM2 immunity protein. After docking prediction were performed groups, taking part in this event were invited to predict the positions of the water molecules near the interaction interface of the two proteins. Figure 3.12 shows the crystall structure of this complex that was published after the 23rd round of CAPRI ended. Each contestant was required to submit 10 models that incuded predicted complex and water molecules. We obtained each model as follows. First, we constructed a protein-protein complex by homology using Modeller [START_REF] Eswar | Comparative protein structure modeling using Modeller[END_REF]. Then, we refined the protein-protein interface using our protein-protein docking potentials. Finally, we immersed the complex into the water box and minimized the value of the score, obtained by the summation over all interactions between water pseudo-atoms and other atom types, with respect to the positions of water oxygens.

The criteria of model quality was based on definition of water-mediated contacts between ligand and receptor. A residue of receptor and a residue of ligand were assumed to have water-mediated contact when one of their heavy atoms were closer than 3.5 Å to the same water oxygen. The quantity f wmc (nat) was computed as the fraction of native water-mediated contacts, predicted by the model. Additionally the organizers of the competition computed number of native and model interface water molecules and number of clashes in each model. The two water molecules assumed to clash if the distance between their oxygen atoms is less than 2.5 Å. All models were then checked for the number of clashes and those where the number of clashing molecules exceeded the number of native interface water molecules were rejected. Other models were ranked according to the criteria in the Table 3

.3. 0 Bad 0 ≤ f wmc (nat) < 0.1 + Fair 0.1 ≤ f wmc (nat) < 0.3 2+ Good 0.3 ≤ f wmc (nat) < 0.5 3+ Excelent 0.5 ≤ f wmc (nat) < 0.8 4+ Bad 0.8 ≤ f wmc (nat) ≤ 1.0 Table 3
.3: Classification of model quality according to the fraction of predicted water-mediated contacts.

Among the models that we submitted 9 were of good and one of fair quality according to the criterium described above. Overall our approach was ranked 4th among others competing techniques. The group ranked first (Nakamura et all ) and the group ranked third (Zou et all ) were using the water molecules that were derived from the interfaces of the homologs as the initial positions of their predictions. Afterwards they optimized the positions of water molecules using the AMBER forcefield. The second-ranked group led by Zacharias used ab-initio water prediction technique. They combined energy-like scoring functions with the well-established forcefield (AMBER), energy minimization and short molecular dynamics runs to generate the final predictions. It is worth to note that our method generated almost no false-positive predictions compared to other competing approaches ([79], SI). The key property of a scoring function is the existence of the correlation between the score of a structure and its similarity to the corresponding native structure. Conventionally, the lRMSD is taken as the measure of similarity of the decoys to the native structure. lRMSD is the ligand (the smaller protein in a complex) root mean square deviation of C α atoms of a decoy relative to the native complex structure when receptors (the larger proteins in the complexes) are superposed. To verify that our potentials indeed correlate with the similarity to the native structures, we plotted the ConvexPP score of each decoy versus the lRMSD for all decoys from the ZDOCK and Rosetta benchmarks. Figure 3.13 shows some typical plots for the complexes from the training set and the two benchmarks. Typically, in the training set we see a wide separation between native and non-native structures. This happens because decoys in the training set have only a few near-native structures with lRMSD<10 Å. On the contrary, about 28% of the Rosetta decoys are the near-native structures.

Discussion

The ZDOCK benchmark has few near-native decoys compared to Rosetta, only 1.5% of the decoys have the near-native conformations.

For the decoys from the ZDOCK and Rosetta benchmarks, we computed the Pearson correlation between lRMSD and the score using near-native decoys with lRMSD<10 Å. In this calculation, we considered only the structures that have at least 50 such decoys. All 54 complexes from the Rosetta benchmark fulfilled this criterion whereas only 24 complexes from the ZDock benchmark had more than 50 near-native decoys. The average Pearson correlation for the Rosetta benchmark is 0.31 whereas for the ZDock benchmark it is approximately 0.21. We investigated the reason of this discrepancy by looking at the atom-pairs distance distributions of the decoys from the training set, the ZDock and Rosetta benchmarks. For our analysis, we used decoys of the three common structures from these sets, PDB codes 1PPE, 1CGI, 1ACB. For these decoys, we computed the average total number of atom pairs at a certain distance using Eq. 3.4. Then, we normalized these distributions on a reference number of pairs αr 3 . We tuned coefficient α in such a way to make each plot approach the value 1.0 when the distance approaches 12.0 Å. Figure 3.14 plots normalized atom-pairs distance distributions for the three abovementioned complexes. From this figure we see that the distance distributions for the Rosetta benchmark are much closer to the native distributions compared to the ZDOCK and training set distributions. We can also see that the Hex docking program [START_REF] Kemp | Protein docking using spherical polar Fourier correlations[END_REF], which we used for the generation of the training set, produces fewer short-distance atom contacts compared to ZDOCK. Since Rosetta decoys were additionally minimized using the Rosetta scoring function, they do not have shortdistance atom contacts and generally their distance distributions resemble the native statistics. Native structures neither have statistics at short distances. Therefore, reconstructed potentials in the vicinity of zero are not reliable and can not provide fair scores for e.g. decoys generated with ZDOCK, since these decoys have many short-distance contacts. Ideally, one needs to additionally penalize short-distance contacts using, e.g., empirical scores that cannot be obtained with statistics from the native structures. However, in this study we do not attempt to provide such additional penalization and only focus on the potentials directly obtained from the training set of the protein structures.

Filtering

Some knowledge-based potentials are smoothed with a smoothing filter a posteriori. For example, Mitchel et al. [START_REF] John | BLEEP-potential of mean force describing protein-ligand interactions: I. Generating potential[END_REF] and Huang et al. [START_REF] Huang | An iterative knowledge-based scoring function for protein-protein recognition[END_REF] used a "1:2:4:2:1" filter, DOPE potential is smoothed using cubic polynomials [START_REF] Shen | Statistical potential for assessment and prediction of protein structures[END_REF], etc. On the contrary, our method introduces an assumption about interaction pair distance uncertainty a priory. More specifically, we collect statistics using gaussian events of a variance σ (Eq. 3.4). We determine the value of the variance from the afore-mentioned crossvalidation procedure. Then, according to Eqs. 3.5-3.6, ConvexPP scoring function is smooth by construction. In other words, we do not need to apply a smoothing filter to the obtained potentials, since we introduce the uncertainty when we collect statistics.

Another parameter that indirectly influences the smoothness of the resulting potential is the regularization parameter C (3.23). According to Eq. 3.19, the scoring vector w, from which the scoring potentials are derived, is a weighted sum of the support structure vectors x ij . The more support structure vectors are in the sum, the more regular the scoring vector w will be. On the other hand, this number equals to the number of non-zero Lagrange multipliers λ ij (Eq. 3. [START_REF] Ashley M Buckle | Protein-protein recognition: Crystal structural analysis of a barnase-barstar complex at 2.0[END_REF] , which is uniquely defined by the value of the regularization parameter C [START_REF] Pontil | Properties of support vector machines[END_REF]. Decreasing C results in the increase of the number of non-zero λ ij therefore resulting in smoother scoring potentials. We also determine the value of this parameter by the crossvalidation procedure.

The consistent determination of the two parameters σ and C allows us to obtain smooth potentials (Eq. 3.6) directly as the solution of the optimization problem (Eq. 3.17).

Uniqueness of the solution and the reference state

The concept of the statistical knowledge-based potentials is based on the definition of two states: the observed state and the reference state [START_REF] Tanaka | Medium-and long-range interaction parameters between amino acids for predicting three-dimensional structures of proteins[END_REF][START_REF] Miyazawa | Estimation of effective interresidue contact energies from protein crystal structures: quasi-chemical approximation[END_REF][START_REF] Sippl | Calculation of conformational ensembles from potentials of mena force: an approach to the knowledge-based prediction of local structures in globular proteins[END_REF]. The observed state is usually the state when a single protein or a complex has the native conformation. It can be derived from the crystal structures. Reference state was introduced as an atom pair distance distribution when the interactions between the atom pairs are absent. The knowledge-based potential is then expressed in terms of these two states as:

u ij (r) = -RT ln N obs ij (r)/N obs ij N ref ij (r)/N ref ij ,
where N ref ij (r) and N obs ij (r) are the numbers of atomic pairs i, j at a distance r in the reference and observed states, correspondingly, and numbers N ref ij and N obs ij are the total numbers of pairs i, j in these states. Some widely used approaches to derive the reference state for protein folding are the ideal-gas approximation [START_REF] Zhou | Distance-scaled, finite ideal-gas reference state improves structure-derived potentials of mean force for structure selection and stability prediction[END_REF], the shuffling of atoms [START_REF] Rykunov | Effects of amino acid composition, finite size of proteins, and sparse statistics on distance-dependent statistical pair potentials[END_REF], a random-walk chain [START_REF] Zhang | A novel side-chain orientation dependent potential derived from random-walk reference state for protein fold selection and structure prediction[END_REF], etc. For protein docking Chuang et al. used decoys as the reference state [START_REF] Chuang | DARS (Decoys As the Reference State) potentials for proteinprotein docking[END_REF], Bernard and Samudrala took the average over the atomic pairs and a cumulative distribution function for all pairs as two reference states [START_REF] Bernard | A generalized knowledge-based discriminatory function for biomolecular interactions[END_REF], etc. The very wide variety of approaches to derive the reference state has its roots in the loose definition and the complexity of the problem.

Recently, the new algorithms that avoid the reference state calculation appeared. We should mention the iterative scheme used by Huang et al. [START_REF] Huang | An iterative knowledge-based scoring function for protein-protein recognition[END_REF] and the neural network classifier by Chae et al. [START_REF] Chae | Predicting protein complex geometries with a neural network[END_REF]. These algorithms indeed avoid the definition of the reference state. However, they do not guarantee the uniqueness of their solution. On the contrary, we showed that our algorithm converges to the global minimum of the function (Eq. 3.17). Thus, we avoid dependence on the initial guess of the interaction potential.

Chapter 4 Conclusion and outlook

The aim of the work was the development and validation of new algorithms that could lead to certain advances in the fields of exhaustive rigid-body search and scoring of protein-protein complexes conformations. During this work the Hermite fitting algorithm was proposed that is not only competable with the state-of-art approaches to this problem, but adds a new class of algorithms, that operate in Hermite functions space, to the existing ones, that function if spherical harmonics and grid representations of a function. A new algorithm to obtain a scoring functions was proposed that is derived from basic logical considerations about the nature of the training dataset. It avoids common unsolved problems with the reference state and has a valuable property of global convergence. It was applied with success to the problems of protein-protein conformations scoring as well as to the prediction of positions of crystallographic water molecules at the protein-protein interaction interface. It was validated using well established benchmarks and community-wide critical prediction of protein interaction assessement challenge.

Future developments

With the advent of the post-genomic era, the cost of whole-genome sequencing plummets and the number of sequenced organisms grows rapidly. However, the proteomics field still did not step into the phase of the exponential growth. Therefore I believe, that probing protein-protein interactions on the scale of proteome will revolutionize the fields of interactomics, evolution and systems biology. The methods that currently applied to discover protein interaction network are either large-scale, but give big number of false-positives and negatives or fit for the discovering detailed picture of single protein-protein complex. Bridging the gap between these two classes of methods would be a tremendous leap forward in the field of protein-protein interaction prediction.

Therefore it is worth trying to optimize current rigid-body search algorithms to scale the computations up. The parallel computing paradigm and especially general GPU programming could lead to a considerable HermiteFit algorithm runtime reduction.

In the area of scoring functions I believe that the paradigm of predefined atom types should be overcomed. The newly emerging area of dimentionality reduction and deep learning surely will bring new advances to the scoring field.

Another direction of expanding this work is to integrate the developed algorithms into one user-friendly package for the convenient use by those who do not posess programming skills required to implement them.

Next, SMO computes the constrained maximum by clipping the unconstrained maximum to the ends of the line segment:

λ new,clipped 2 =      L, if λ new 2 ≤ L λ new 2 , if L <

1. 1

 1 Example of a homodimeric complex . . . . . . . . . . . . . . . . . . . 1.2 Example of a heterodimeric complex . . . . . . . . . . . . . . . . . . 1.3 Y2H explanation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 TAP-MS explanation . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 FRET explanation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.6 Isothermal titration calorimetry explanation . . . . . . . . . . . . . . 1.7 The structure of Twist protein with the Brd4 protein . . . . . . . . . 1.8 Co-regulation scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.9 Co-evolution scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.10 FFT-based docking example . . . . . . . . . . . . . . . . . . . . . . . 1.11 Water molecules placement around aromatic nitrogen in histidine residue functional group. . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 Flowchart of the standard fitting algorithm . . . . . . . . . . . . . . . 2.2 Flowchart of HermiteFit . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Hermite functions examples . . . . . . . . . . . . . . . . . . . . . . . 2.4 Hermite functions and sine function . . . . . . . . . . . . . . . . . . . 2.5 Shifted 1D step-function η . . . . . . . . . . . . . . . . . . . . . . . . 2.6 Rotations of the Hermite expansion . . . . . . . . . . . . . . . . . . . 2.7 Matrices F (1) and F (2) . . . . . . . . . . . . . . . . . . . . . . . . . . 2.8 Transfer T -matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.9 Analytical R-factors . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.10 R-factors in one and three dimensions . . . . . . . . . . . . . . . . . . 2.11 Results for the alpha-conotoxin PnIB . . . . . . . . . . . . . . . . . . 2.12 Result of the fitting of the GroEL EDM . . . . . . . . . . . . . . . . 2.13 Running times of the Hermite to Fourier space transition and FFTW 3.1 Two types of orthogonal functions . . . . . . . . . . . . . . . . . . . . 3.2 Structure vectors for three complexes . . . . . . . . . . . . . . . . . . 3.3 Two classes of structure vectors for a single complex . . . . . . . . . . 3.4 The flowchart of block sequential minimal optimization algorithm . . 3.5 Predictive performance of the scoring potential as a function of the smoothing parameter σ and the regularization parameter C . . . . . . 3.6 Performance rate of the scoring potentials on the training set versus the number of iterations of the BSMO algorithm . . . . . . . . . . . . 3.7 The scoring potentials trained in two different polynomial bases . . . 3.8 Dependence of the extracted scoring potentials on the order of the decomposition in the Legendre basis . . . . . . . . . . . . . . . . . . . 3.9 The success rate on ZDOCK benchmark . . . . . . . . . . . . . . . . 3.11 Knowledge-based solvation scoring functions . . . . . . . . . . . . . . 3.10 the success rate on Rosetta benchmark . . . . . . . . . . . . . . . . . 3.12 The crystallographic X-Ray structure of DNase domain of colicin E2 and IM2 immunity protein . . . . . . . . . . . . . . . . . . . . . . . . 3.13 The plots of the ConvexPP versus the lRMSD for the decoy structures 3.14 The normalized atom-pairs distance distributions for three complexes 1PPE, 1CGI, and 1ACB . . . . . . . . . . . . . . . . . . . . . . . . .

4

  Conclusion and outlook 4.0.7 Future developments . . . . . . . . . . . . . . . . . . . . . . .

Figure 1 . 1 :

 11 Figure 1.1: Example of a homodimeric complex, D-alanyl carrier protein, pdb code 4BPG.

Figure 1 . 3 :

 13 Figure 1.3: Schematic representation of the Y2H method. If the two protein X and Y do interact, expression of the reporter gene is high (yellow star), otherwise AD and BD domains of the promoter are separated and the expression is low. UAS stands for the upstream regulating sequence.

Figure 1 . 4 :

 14 Figure 1.4: Schematic representation of the tandem affinity pyrification method. Contaminants are shown with green stars, IgG binding domain is shown with pink pentagon (as well as the IgG coating) and calmodulin (as well as calmodulin coating) is shown with squares.

Figure 1 . 5 :

 15 Figure 1.5: Schematic representation of the fluorescence resonance energy transfer method. If the two proteins interact, the fluorophores (blocks) come close and the energy transfer between them becomes possible. Therefore, exciting one fluorophore, one detects the emission from another. Otherwise, if the proteins do not interact, no emission of the second fluorophore is detected upon excitation of the first one.

Figure 1 . 6 :

 16 Figure 1.6: Schematic representation of the data obtained from an isothermal titration calorimetry device. The red line shows the energy transfered from the reservoir with the protein to the reference reservoir. Blue dashed line is an example of model fitted to the data.

Figure 1 . 7 :

 17 Figure 1.7: The structure of Twist protein(green) with the bromodomain-containing protein 4 (red) measured using NOE-NMR technique [124]. PDB code 2MJV.

Figure 1 . 10 :

 110 Figure 1.10: Example of a FFT-based docking procedure.

  DE ] = Re [L DE ] = sum of ACE and PSC scores for nearby atoms empty space 0 otherwise Im [R DE ] = Im [L DE ] = 1 grid point is the nearest grid point of an atom 0 otherwise ,

  Fig 1.11 shows an examples of water molecules placement around aromatic nitrogen in histidine.

Figure 1 . 11 :

 111 Figure 1.11: Water molecules placement around aromatic nitrogen in histidine residue functional group.

Figure 2 . 1 :

 21 Figure 2.1: Flowchart of the standard fitting algorithm based on the Fourier correlations. Green blocks correspond to the operations in the Fourier space.Density Map

Figure 2

 2 Figure 2.2: Flowchart of HermiteFit, the new fitting algorithm based on the Hermite expansions. Green blocks correspond to the operations in the Fourier space. Blue blocks correspond to the operations in the Hermite space.

Figure 2 . 3 :

 23 Figure 2.3: Left: 1D Hermite functions of order six for three different scaling parameters λ. Right: 1D Hermite functions of two different orders for the scaling parameter λ = 1.

Figure 2

 2 Figure 2.2.2 shows that Hermite functions are very similar to the cosine functions near the coordinate axis origin.

Figure 2 . 4 :

 24 Figure 2.4: Hermite function of the order 50 (solid line) and cosine function with the frequency equal to 50 (dashed line).

Figure 2 . 5 :

 25 Figure 2.5: Shifted 1D step-function η.

Figure 2 . 6 :

 26 Figure 2.6: Sequential rotations of coefficients fn,m,l about different axes. The rotated layer is shown with the solid cubes, other coefficients are shown with the dashed cubes.To perform the complete rotation of the decomposition about one axis, we rotate each layer of coefficients about the corresponding axis in the space of coefficients.

Figure 2 . 7 :

 27 Figure 2.7: Absolute values of two matrices F (1) and F (2) that give the transfer matrix as their product (Eq. 2.50). These matrices are computed with the scaling parameter λ = 0.55 and the input box size of L box = 23.0 Å, which mimics the first fitting example shown below. Left: F (1) , the scaled Fourier transform of a 1D Hermite function as given by Eq. 2.51. Right: F (2) , the scaled Fourier series of a 1D Hermite function as given by Eq. 2.52. The dashed blue line highlights the maximum encoded frequency according to Eq. 2.38. The solid black line on the right plot shows the maximum Hermite decomposition order N max , at which the two matrices are still identical (Eq. 2.37).

Figure 2 .

 2 Figure 2.3.2 shows several examples of the absolute values of the transfer matrix components for three different values of the Hermite scaling parameter λ and three values of the Hermite decomposition order N . The size of the transfer matrix was limited to 60 × 60 and the box size L box was set to 23 Å.The ideal transfer matrix should be identity, which is the case only at N → ∞, as we demonstrate below. We see, however, that the transfer matrix at small values of λ encodes only loworder reflexes. The index of the last encoded reflex can be estimated from Eq. 2.38 as k max = √ 2N + 1λL box /(2π). With the increase in order N and parameter λ, the number of encoded frequencies rises. At the same time, increasing the scaling parameter λ makes the quality of encoding of all the frequencies worse, as we see in the right column. Therefore, it is very important to tune the value of λ according to the class of input functions, such that the quality of encoding becomes optimal. Below we will assess encoding quality by means of the crystallographic R-factor.

Figure 2 . 8 :

 28 Figure 2.8: Nine examples of the absolute values of the transfer T -matrices for three different values of λ and three different values of the Hermite decomposition order N . The number of Fourier coefficients is M = 60, and the input box size is L box = 23.0 Å, which mimics the first fitting example shown below. Hermite decomposition orders are N ∈ {15, 20, 30}, parameter λ takes the values of 0.3 Å-1 , 0.55 Å-1 , and 1.0 Å-1 . The first column corresponds to the relative λL box value of 6.9, the middle column corresponds to the relative λL box value of 12.65, and the right column to the relative λL box value of 23. Notably, at low values of λ the transfer matrix encodes only small order reflexes. The index of the last reflex can be estimated from Eq. 2.38 as k max = √ 2N +1λL box 2π

  Figure 2.3.4 shows analytical R-factors in one dimension computed using Eqs. 2.48 and 2.60 as a function of the Hermite decomposition order N and the scaling parameter λ. We bounded the input and output frequencies by M = 30 Fourier coefficients. The size of the input interval L box is set to 23.0 Å to mimic the alpha-conotoxin PnIB peptide (pdb code 1AKG) decomposition used in the fitting example below.

Figure 2 . 9 :

 29 Figure 2.9: Analytical R-factors in one dimension as a function of Hermite decomposition order N and scaling parameter λ computed at three different resolutions.The input signal is modelled as a sum of Gaussians (Eq. 2.6) with the variance of α/ √ 2 equispaced at a distance α. The number of Fourier coefficients is M = 30, and the input box size is L box = 23.0 Å. These values are chosen to mimic the 1AKG peptide decomposition. The estimate on the optimal parameter λ (Eq. 2.41) is plotted with the red dashed line. The real-space bound on the optimal parameter λ (Eq. 2.37) is shown with the orange dashed line. The reciprocal-space bound on the optimal parameter λ (Eq. 2.39) is shown with the blue dashed line. Left: The Gaussian parameter α = 0.2 Å, corresponding to the absolute input signal resolution of R = 0.31 Å and the relative resolution of R/L box = 0.014. However, in this case, the actual absolute resolution is cut at L box /M = 0.77 Å, which corresponds to the relative resolution of 0.033. Middle: The Gaussian parameter α = 1.0 Å, corresponding to the absolute input signal resolution of R = 1.57 Å and the relative resolution of R/L box = 0.068. Right: The Gaussian parameter α = 5.0 Å, corresponding to the absolute input signal resolution of R = 7.85 Å and the relative resolution of R/L box = 0.34.

Figure 2 . 10 :

 210 Figure 2.10: Analytical R-factors in one and three dimensions as a function of relative resolution R/L box . The absolute resolution at box size L box = 23 Å is shown in the top horizontal axis. Plots for three different Hermite expansions orders are shown, N ∈ {15, 20, 30}. Parameters L box and M were constant and equal to 23 Å and 30, correspondingly. Scaling parameter λ was estimated using Eq. 2.41.

Figure 2 . 11 :

 211 Figure 2.11: Test of the fitting algorithm on artificially generated EDM for the alpha-conotoxin PnIB (PDB code 1AKG). Here, we plotted the dependence of four parameters, the maximum score, the average score, the score of the near-native conformation and the crystallographic R-factor on the scaling parameter λ. Isosurface of the Hermite decomposition at protein model density equal to (ρ max + ρ min )/2 and several values of λ are shown in sub-plots A (λ = 0.05), B (λ = 0.15), C (λ = 0.3), D (λ = 0.55) and E (λ = 2.0).

Figure 2 . 12 :

 212 Figure 2.12: Result of the fitting chain A of the GroEL-GroES X-Ray structure (PDB entry 1AON) to the GroEL complex electron density map (EMD-2001). Two heptameric rings are shown in different colors. The average RMSD measured using the C α -atoms between the two closest chains in the fitted structure and the flexibly refined structure provided by the authors of the EDM (PDB entry 4AAU) is 5.35 Å .

Figure 3

 3 Figure 3.1: Two types of orthogonal functions. Left: shifted Legendre polynomials orthogonal on the interval [0; 10]. Right: shifted rectangular functions.

Figure 3 . 3 :

 33 Figure 3.3: Two classes of structure vectors for a single complex are shown. Native structure vectors are plotted as blue circles. Nonnative structure vectors are plotted as red squares.A) The case when infinitely many hyperplanes can separate the two classes. B) The case when no optimal separating hyperplane exists. Slack variables ξ i and ξ j for misclassified structure vectors are added, which are the distances to the corresponding margin hyperplanes. The optimal hyperplane, which maximizes the separation between the two classes, is plotted as a dashed line. Two margin hyperplanes are plotted as solid lines.

Figure 3 . 4 :

 34 Figure 3.4: The flowchart of block sequential minimal optimization algorithm.

Figure 3 . 5 :Figure 3 . 6 :

 3536 Figure 3.5: Predictive performance of the scoring potential as a function of the smoothing parameter σ and the regularization parameter C. A) Performance obtained if the scoring functions are trained on the whole database and verified on the same database. B) Performance obtained if the scoring functions are trained on 200 protein complexes and verified on the other 650 complexes from the training database. Here the best performance is obtained with σ = 0.4 Å and C = 10 6 . . . 10 7 .

Figure 3 . 7 :

 37 Figure 3.7:The scoring potentials trained in two different polynomial bases. Dashed lines correspond to the scoring potentials that were obtained using the Legendre basis functions. Solid lines correspond to the potentials that were obtained using the rectangular basis functions. Left: Potential between aliphatic carbons bonded to carbons or hydrogens only. Right: Potential between a guanidine nitrogen with two hydrogens and an oxygen in carboxyl groups.

Figure 3 . 8 :

 38 Figure 3.8: Dependence of the extracted scoring potentials on the order of the decomposition in the Legendre basis. After order P = 25, the potentials are indistinguishable from each other and thus not shown for clarity. Left: Potential between aliphatic carbons bonded to carbons or hydrogens only. Right: Potential between a guanidine nitrogen with two hydrogens and an oxygen in carboxyl groups.

Figure 3 . 9 :

 39 Figure3.9: Dependence of the success rate on ZDOCK benchmark on the number of top predictions in consideration for the three methods.

Figure 3 . 11 :

 311 Figure 3.11: Knowledge-based solvation scoring functions U k (r) between the oxygen of a water molecule and protein atoms as a function of the separation distance. Left: water -aromatic carbon scoring function is plotted. Middle: water -carboxyl oxygen (like in aspartic and glutamic acids) scoring function is plotted. Right: water -charged nitrogen (like in lysins) scoring function is plotted.

Figure 3 . 12 :

 312 Figure 3.12: The crystallographic X-Ray structure of DNase domain of colicin E2 and IM2 immunity protein (PDB code 3U43). IM2 imunity protein is show in green and DNAse domain in red. Blue spheres denote the positions of interfacial water oxygens. The interfacial water molecule is the one within 3.5 Å of both receptor and ligand.

Figure 3 . 13 :

 313 Figure 3.13: The plots of the ConvexPP versus the lRMSD for the decoy structures from the training set (1A7N, 1J7L), Rosetta benchmark (1SPB, 2KAI) and ZDOCK benchmark(1N8O, 2CHF). On the left we show the plots that exhibit funnel-like behaviour near the frame origin. On the right side the plots without obvious funnels are shown.

Figure 3 . 14 :

 314 Figure 3.14: The normalized atom-pairs distance distributions for three complexes 1PPE, 1CGI, and 1ACB. For each complex, four plots are shown: average ZDOCK distribution, average Rosetta distribution, average training set distribution and the native distribution. The average is taken over all decoys from the two benchmarks and the training set.

  2.1 Complexity of the Hermite fitting algorithm . . . . . . . . . . . . . . 2.2 Comparison of the models obtained using HermiteFit, Colores and ADP EM algorithms by RMSD . . . . . . . . . . . . . . . . . . . . . 2.3 Comparison of the HermiteFit algorithm with the Colores and ADP EM algorithms by runtime . . . . . . . . . . . . . . . . . . . . . . . . . . Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.1 Alpha-conotoxin PnIB . . . . . . . . . . . . . . . . . . . . . . 2.4.2 GroEL complex . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.3 Runtime of Hermite-to Fourier-space transition . . . . . . . . 2.4.4 Comparison with Situs and ADP EM . . . . . . . . . . . . . . Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Problem Formulation . . . . . . . . . . . . . . . . . . . . . . . 3.2.2 Expansion of U (r) and n(r) in an orthogonal basis . . . . . . 3.2.3 Geometrical interpretation and connection to quadratic programming . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.4 Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.5 Training database for protein-protein interactions . . . . . . . 3.3 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . .

	3 Scoring functions for protein-protein docking
	3.1

3.1 ZDock benchmark 3.0 results . . . . . . . . . . . . . . . . . . . . . . 3.2 Rosetta unbound benchmark results . . . . . . . . . . . . . . . . . . . 3.3 Classification of model quality according to the fraction of predicted water-mediated contacts. . . . . . . . . . . . . . . . . . . . . . . . . . C.1 ZDock benchmark results with homologs . . . . . . . . . . . . . . . . C.2 Rosetta unbound benchmark results with homologs . . . . . . . . . . 2.4 3.3.1 Overfitting and Convergence . . . . . . . . . . . . . . . . . . . 3.3.2 Extracted Potentials . . . . . . . . . . . . . . . . . . . . . . . 3.3.3 Protein-Protein docking benchmark version 3.0 . . . . . . . . 3.3.4 Rosetta Benchmark . . . . . . . . . . . . . . . . . . . . . . . . 3.3.5 Water interaction potentials prediction for CAPRI T47 target 3.3.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Table 2

 2 .2).

	Algorithm RMSD C α , Å RMSD centers of masses,	Å
	ADP EM	4.61	2.29
	Colores	5.42	2.52
	HermiteFit	5.35	2.64

.2: Comparison of the models obtained using HermiteFit, Colores and ADP EM algorithms with the model obtained by the authors of the electron density map (PDB entry 4AAU

Table 2 .

 2 2 demonstrates that the tested algorithms produce equal quality models. However, results of ADP EM are slightly better, presumably because of the finer rotational sampling. Num of rot-space points Num of trans-space points Runtime, s Time per point, ×10 -7 s Table2.3: Comparison of the HermiteFit algorithm with the Colores and ADP EM algorithms. The comparison criterion was chosen to be the total running time and the running time per one point of the search space.

	Algorithm ADP EM	16384	23186	139	3.6
	Colores	4416	1336965	1454	2.5
	HermiteFit	4416	1336965	917	1.5

Table 3

 3 potentials(15 vs 12 hits). Although if we consider top 8 predictions our scoring function outperforms ZRANK (32 vs 26 hits) and gives the same number of hits for top 5 to 8 predictions. Excluding homologs from the training set results in a slight improvement of the results (TableC.1).Figure3.9 shows ROC curves (success rate vs the number of top predictions considered). We see that ConvexPP scoring functions outperform ZRANK and ZDOCK if the number of considered predictions is more than eight.

	.1 shows results of ZDOCK3.0, ZRANK and our
	scoring functions on the ZDOCK benchmark.

  Table3.2 compares the results of RosettaDock[START_REF] Gray | Protein-protein docking with simultaneous optimization of rigid-body displacement and side-chain conformations[END_REF], ITScore-PP[START_REF] Huang | An iterative knowledge-based scoring function for protein-protein recognition[END_REF] and our ConvexPP scoring functions.

		2JEL	RosettaDock 2 1 4.82	2.22	0.52		ITScore-PP 2 1 6.55	1.98	0.65		ConvexPP 2 1 6.42	1.26	0.86
	Complex	Quality 2KAI 2PCC 2PTC	Rank	1 3 1	L rmsd 2 2 4	I rmsd 2.02 9.19 0.82	F nat 0.97 3.39 0.44	Quality 0.67 0.24 0.80	Rank	2 3 2	L rmsd 2 1 4	I rmsd 2.48 9.38 5.86	F nat 1.05 3.84 1.82	0.70 Quality 0.29 0.70	Rank	1 3 2	8 rmsd L 1 1	2.26 I 9.40 rmsd 5.30	0.89 4.31 F nat 1.16	0.88 0.48 0.74
		2SIC		2	1	5.18	1.53	0.85		2	1	5.18	1.53	0.85		1	1	4.84	0.96	0.82
	1ACB	3	3		9.08	3.41	0.13	3	4		9.04	3.39	0.11	2	1		4.31	1.13	0.78
		2SNI		2	1	6.70	2.01	0.60		1	1	2.88	0.99	0.73		2	1	7.38	2.00	0.61
	1A0O	2	1		10.29	4.18	0.57	2	1		10.29	4.18	0.57	3	1		11.33	5.79	0.35
		2TEC		1	1	2.46	0.81	0.74		1	1	2.59	0.85	0.78		1	1	1.97	0.59	0.80
	1AHW	2	1		6.41	2.37	0.51	3	4		6.46	3.00	0.49	2	1		4.82	2.31	0.44
		2VIR		3	1	7.53	4.19	0.26		2	2	5.74	2.17	0.70		2	1	5.78	1.03	0.67
	1ATN	3	1		5.83	2.79	0.49	3	1		9.49	4.70	0.38	1	1		2.34	0.71	0.77
		3HHR		3	50	9.84	3.95	0.26		3	3	8.17	4.03	0.30		3	1	9.41	3.47	0.33
	1AVW	2	1		6.02	2.09	0.67	2	1		5.07	1.81	0.71	2	1		6.01	1.38	0.75
		4HTC		2	1	3.81	1.54	0.61		3	1	5.95	2.25	0.42		2	1	4.04	1.50	0.76
	1AVZ	3	37		8.08	4.05	0.29	3	22		11.47	5.55	0.35	3	20		8.41	4.41	0.12
						Top1: 66.7% (36/54)					Top1: 59.3% (32/54)					Top1: 83.3% (45/54)
	1BQL	1	1		1.57	0.86	0.64	1	5		1.81	0.72	0.65	1	1		1.40	0.96	0.89
		Homologs			Top1 and quality 1: 16.7% (9/54)			Top1 and quality 1: 11.1% (6/54)			Top1 and quality 1: 20.4% (11/54)
	1BRC	2	4		3.77	1.21	0.75	2	1		3.77	1.21	0.75	2	1		7.24	1.62	0.76
		included		Top1 and quality 1 or 2: 48.1% (26/54)		Top1 and quality 1 or 2: 38.9% (21/54)		Top1 and quality 1 or 2: 57.4% (31/54)
	1BRS	2	1		4.78	1.73	0.64	3	1		8.68	4.46	0.33	3	1		8.70	3.79	0.31
			Top10 and quality 1 or 2: 61.1% (33/54)	Top10 and quality 1 or 2: 57.4% (31/54)	Top10 and quality 1 or 2: 63.0% (34/54)
	1BTH	3	4		18.21	5.54	0.30	3	2		5.60	2.35	0.42	3	1		5.59	2.67	0.44
																		Top1: 81.5% (44/54)
	1BVK	3	1		7.91	3.93	0.20	3	1		7.18	3.54	0.20	3	1		7.75	3.45	0.22
		Homologs															Top1 and quality 1: 16.7% (9/54)
	1CGI	2	2		3.79	1.86	0.50	3	8		6.01	2.37	0.42	3	3		6.44	2.11	0.38
		excluded														Top1 and quality 1 or 2: 55.6% (30/54)
	1CHO	3	1		6.31	2.19	0.46	3	1		10.32	3.76	0.18	2	1		6.35	1.81	0.66
															Top10 and quality 1 or 2: 61.1% (33/54)
	1CSE	2	6		10.10	3.12	0.56	2	1		8.81	2.66	0.71	3	1		7.87	2.29	0.40
	1DFJ	2	1		5.69	2.66	0.59	2	1		5.69	2.66	0.59	2	1		5.63	2.55	0.67
	1DQJ	3	1		6.71	3.35	0.31	3	5		5.00	2.12	0.34	3	1		14.01	6.06	0.38
	1EFU	3	44		5.98	3.78	0.16	3	26		7.83	4.21	0.10	3	20		5.90	4.65	0.11
	1EO8	3	1		10.73	5.54	0.31	3	31		6.12	3.36	0.15	3	1		10.90	3.29	0.42
	1FBI	2	1		2.79	1.37	0.54	2	1		3.64	1.86	0.51	3	1		11.03	4.23	0.36
	1FIN	3	364		9.88	5.38	0.12	3	109		8.26	4.19	0.12	3	200		8.27	6.06	0.10
	1FQ1	3	1		11.37	5.43	0.31	3	1		9.33	5.09	0.41	3	1		6.79	4.02	0.43
	1FSS	1	1		3.07	0.97	0.74	2	1		4.46	1.42	0.46	2	1		3.89	1.54	0.63
	1GLA	2	4		5.96	1.85	0.65	3	7		11.96	3.83	0.35	2	1		6.20	1.70	0.84
	1GOT	3	12		7.86	3.95	0.19	3	4		9.71	4.26	0.19	3	2		12.91	3.21	0.17
	1IAI	3	8		6.60	3.42	0.22	2	1		4.18	1.61	0.62	2	1		3.82	1.62	0.33
	1IGC	1	1		2.15	0.63	0.85	1	1		2.15	0.63	0.85	1	1		1.92	0.56	1.00
	1JHL	3	3		8.56	4.51	0.26	2	2		6.09	2.66	0.56	3	2		7.65	3.94	0.31
	1MAH	2	1		3.72	1.19	0.70	3	1		8.66	2.77	0.26	2	1		3.71	1.15	0.78
	1MDA	3	1		8.77	3.59	0.19	3	1		9.84	3.92	0.26	2	2		6.65	2.38	0.52
	1MEL	2	1		8.47	2.62	0.50	2	4		8.47	2.62	0.50	2	1		7.89	2.84	0.52
	1MLC	2	7		4.91	1.37	0.52	3	1		15.36	3.45	0.20	3	1		15.41	3.04	0.24
	1NCA	2	1		3.06	1.53	0.61	1	1		1.24	0.64	0.75	2	1		7.62	2.13	0.66
	1NMB	1	1		0.90	0.44	0.80	1	6		2.66	0.76	0.85	1	1		2.66	0.57	1.00
	1PPE	1	1		1.38	0.52	0.73	3	1		7.42	2.39	0.28	1	1		1.38	0.54	0.89
	1QFU	2	1		3.02	1.10	0.64	1	4		2.93	1.00	0.69	1	1		3.89	0.98	0.67
	1SPB	1	1		1.06	0.62	0.68	1	1		1.47	0.70	0.69	1	1		1.04	0.54	0.82
	1STF	1	1		1.91	0.68	0.89	1	2		1.57	0.54	0.91	1	1		1.57	0.51	0.94
	1TAB	2	1		4.16	1.30	0.74	2	1		4.39	1.37	0.76	2	1		4.11	1.47	0.68
	1TGS	2	1		2.48	1.44	0.59	2	1		2.26	1.38	0.64	3	1		8.41	3.55	0.44
	1UDI	2	1		3.35	1.43	0.63	2	1		2.08	1.01	0.74	2	1		2.08	1.09	0.71
	1UGH	1	1		1.78	0.86	0.67	2	1		4.64	1.90	0.46	2	1		4.61	1.96	0.60
	1WEJ	3	15		9.28	2.92	0.44	2	2		6.90	2.55	0.62	3	1		9.37	4.79	0.23
	1WQ1	3	1		5.53	2.46	0.34	2	1		3.38	1.92	0.39	2	2		3.83	1.59	0.58
	2BTF	3	1		10.03	3.23	0.22	1	1		1.52	0.60	0.75	1	1		1.31	0.60	0.90

Table 3 .

 3 2: Rosetta unbound benchmark results. Proteins homologous to the ones in the training set are shown with the bold font.

  Finally, SMO determines the value of λ 1 from the new clipped value of λ 2 : Table C.1: ZDock benchmark results. Proteins homologous to the ones in the training set are shown with bold font. Absense of hits among first 1000 predictions is shown with hyphens. Table C.2: Rosetta unbound benchmark results. Proteins homologous to the ones in the training set are shown with bold font. Absense of hits among first 1000 predictions is shown with hyphens.

	1STF	1	1EAW 2OT3 1 1.57		-0.51	2.60 -0.94	1.03 -	1	0.70 -	1	1 -	1.57	2.60 -0.51	1.03 -	0.70 -0.94
	1TAB	2	1EWY Homologs 1 4.11		21 1.47 Top1: 10.5% (13/124) 3.16 1.74 0.56 0.68 2	1	49	3.16 Top1: 12.1% (15/124) 1.74 4.11 1.47 0.68 0.56
	1TGS	3	1EZU included 1 8.41		-3.55 Top10: 27.4% (34/124) --0.44 3	-	1		-	-Top10: 29.0% (36/124) -8.41 3.55 0.44	-
	1UDI	2	1	1F34	2.08		38 1.09	3.41 0.71	1.45	2	0.54	1	30	2.08	3.41 1.09	1.45	0.54 0.71
	1HIA 1MAH 1 4.61 1 9.37 Appendix C 1UGH 2 1WEJ 3 1N8O 1WQ1 2 2 3.83		-1.96 4.79 1.59	-3.64 0.60 H, -1.26 0.23 2.94 1.24 0.58	2 3 3	-if λ new -1 2 0.69 2 1 0.74 1 3	λ new 4.61 2 ≥ H -< H 1.96 4.17 1.48 -9.37 4.79 2.94 1.24 8.42 3.29	-0.61 0.60 0.23 0.74 0.35	(B.4)
	2BTF	1	1OPH 1 1.31		-0.60		-	0.90		-	1		-	1	-	1.31	-	0.60	-	0.90	-
	λ new 1 Homologs included = λ old 1 -y 1 y 2 (λ new,clipped 2 4.62 1.52 0.71 1 3.92 -λ old 2 ) 1.43 Homologs excluded 0.74 6.42 1PPE 1 1.26 0.86 2 1 6.42 1.26 0.86 Docking benchmarks results 2JEL 2 1R0R 8.36 2.46 0.40 4 8.36 2.46 0.40 1TMQ 6.11 1.97 0.45 1 6.11 1.97 0.45 1UDI --------1YVB 7.33 1.92 0.71 9 7.33 1.92 0.71 2KAI 1 8 2.26 0.89 0.88 1 7 2.26 0.89 0.88 Complex Quality Rank L rmsd I rmsd F nat Quality Rank L 2PCC 3 1 9.40 4.31 0.48 3 1 9.40 4.31 0.48 rmsd I rmsd 2PTC 2 1 5.30 1.16 0.74 2 1 5.30 1.16 0.74 F nat 2SIC 1 1 4.84 0.96 0.82 1 1 4.01 0.79 0.86	(B.5)
	2SNI 2TEC 2VIR 3HHR 4HTC	1ACB 2 1A0O 1AHW 1 2 Complex 1ATN 3 1AVW 2	With Homologs 2B42 9.44 2 1 4.31 1.13 7.38 2.00 0.61 2MTA 1 125 3.76 2O8V --3 1 11.33 5.79 2 1 4.82 2.31 1 1.97 0.59 0.80 1 5.78 1.03 0.67 Rank L rmsd I rmsd F 2.23 1.88 -nat 2PCC 458 6.29 2.19 1 1 2.34 0.71 1 9.41 3.47 0.33 2SIC 5.92 1.19 2 1 6.01 1.38 1 4.04 1.50 0.76	0.43 0.78 2 Without Homologs 1 9.44 2.23 2 1 4.31 1 7.38 2.00 0.50 197 3.76 1.88 ----0.35 3 1 6.35 0.44 2 1 4.82 1 1 2.10 0.59 2 1 5.78 1.03 Rank L rmsd I rmsd 0.44 622 6.29 2.19 0.77 1 1 2.34 3 1 9.74 3.74 nat F 0.73 2 5.92 1.19 0.75 2 1 6.01 2 1 4.04 1.50	0.43 1.13 0.61 0.50 -3.44 2.31 0.84 0.67 0.44 0.71 0.38 0.73 1.38 0.76	0.78 0.47 0.44 0.77 0.75
		1AVZ	Rigid-Body 8.57 2.45 8.41 4.41 0.12 0.57 Top1: 83.3% (45/54) 2SNI 3 20		2 3 Top1: 81.5% (44/54) 8.57 2.45 16 8.41	0.57 4.41	0.12
	Homologs	1AHW 1BQL Top1 and quality 1: 20.4% (11/54) 547 2.14 0.91 0.79 2UUY 233 13.72 2.38 1 1 1.40 0.96	465 0.66 0.89 Top1 and quality 1: 16.7% (9/54) 2.14 0.91 0.79 346 13.72 2.38 0.66 2 1 3.78 2.35	0.37
	included	1BVK 1BRC Top1 and quality 1 or 2: 57.4% (31/54) ----7CEI 8.41 2.46 2 1 7.24 1.62	-0.65 0.76 Top1 and quality 1 or 2: 55.6% (30/54) ---1 6.39 1.94 0.65 2 1 7.24 1.62	0.76
		1DQJ 1BRS Top10 and quality 1 or 2: 63.0% (34/54) ----1A2K ---3 1 8.70 3.79	-0.31 -Top10 and quality 1 or 2: 61.1% (33/54) -------3 1 8.70 3.79	0.31
		1E6J 1BTH	35 1AK4 3	4.21 -1		2.26 -5.59	0.48 -2.67	63 0.44 -		3.78 -3	1.92 -1	0.68 -5.59	2.67	-	0.44
		1JPS 1BVK	62 1AKJ 3	4.14 395 1	1.11 3.30 7.75 3.45 0.78 1.54	127 0.60 0.22		3.54 388 3	0.98 5.59 1	0.70 1.96 8.72	0.57 3.78	0.22
		1MLC 1CGI	5 1AZS 3		4.70 3		1.12 11.04 6.44 2.11 0.39 1.88	2 0.69 0.38		4.70 5 3	1.12 12.49 5	0.39 1.74 6.46	0.65 2.05	0.47
		1VFB 1CHO	1239 1B6C 2	10.89 1		2.48 3.67 6.35 1.81 0.30 2.23	948 0.90 0.66	10.89 1 2	2.48 3.67 1	0.30 2.23 6.35	0.90 1.81	0.66
		1WEJ 1CSE	1 1BUH 3	4.13 1		1.30 3.42 7.87 2.29 0.75 1.81	1 0.63 0.40		2.20 4 3	1.16 3.42 1	0.75 1.81 7.87	0.63 2.29	0.40
		2FD6 1DFJ	8 1E96 2	15.65 261 1	2.16 6.03 5.63 2.55 0.80 2.18	99 0.58 0.67	15.65 425 2	2.16 5.56 1	0.80 1.99 5.63	0.54 2.55	0.67
		2I25 1DQJ	83 1EFN 3	4.45 -1	1.87 -14.01	0.33 -6.06	195 -0.38		7.92 -3	2.21 -1	0.36 -14.01	6.06	-	0.38
		2VIS 1EFU		617 1F51 3		23.89 16 20		2.37 3.41 5.90 4.65 0.43 1.82	326 0.60 0.11	23.89 11 3	2.37 3.41 14	0.43 1.82 5.90	0.60 4.65	0.11
		1BJ1 1EO8	2 1FC2 3	2.82 -1	0.98 -10.90	0.86 -3.29	8 0.42 -		2.82 -3	0.98 -1	0.86 -10.90	3.29	-	0.42
		1FSK 1FBI	3 1FQJ 3		3.98 -1	1.39 -11.03	0.81 -4.23	2 0.36 -		3.52 -3	1.25 -1	0.74 -11.03	4.23	-	0.36
		1I9R 1FIN	462 1GCQ 3	16.85 118 200	2.30 1.98 8.27 6.06 0.48 1.19	119 0.71 0.10	16.85 141 3	2.30 1.98 240	0.48 1.19 8.27	0.71 6.06	0.10
		1IQD 1FQ1	16 1GHQ 3	4.20 -1		0.97 -6.79	0.67 -4.02	8 0.43 -		6.90 -3	1.95 -1	0.46 -6.79	4.02	-	0.43
		1K4C 1FSS	242 1GLA 2	5.78 57 1		1.31 4.91 3.89 1.54 0.62 2.22	1177 0.37 0.63		5.78 273 2	1.31 4.91 1	0.62 2.22 3.89	0.37 1.54	0.63
		1KXQ 1GLA	1 1GPW 2	3.06 1		1.04 7.10 6.20 1.70 0.88 2.44	1 0.39 0.84		3.06 1 2	1.04 7.10 1	0.88 2.44 6.20	0.39 1.70	0.84
		1NCA 1GOT	12 1HE1 3		4.50 301 2 12.91 1.38 5.93 3.21 0.86 1.74	29 0.38 0.17		4.50 324 3	1.38 5.93 2	0.86 1.74 10.92	0.38 3.22	0.13
		1NSN 1IAI		636 1I4D 2		4.95 -1		2.00 -3.82	0.50 -1.62	409 -0.33		4.95 -2	2.00 -1	0.50 -3.82	1.62	-	0.33
		1QFW 1IGC		1315 1J2J 1	5.12 86 1		1.35 5.59 1.92 0.56 0.73 2.12	1274 0.55 1.00		5.12 97 1	1.35 5.59 2	0.73 2.12 1.92	0.55 0.56	1.00
		1QFW 1JHL		1315 1K74 3	5.12 2		1.35 7.90 7.65 3.94 0.73 2.02	1274 0.48 0.31		5.12 1 3	1.35 6.65 3	0.73 2.30 8.55	0.70 3.90	0.38
		2JEL 1MAH	957 1KAC 2	6.83 287 1	2.30 4.47 3.71 1.15 0.31 2.21	485 0.36 0.78		6.83 105 2	2.30 4.47 1	0.31 2.21 3.55	0.36 1.16	0.75
		1AVX 1MDA	3 1KLU 2	4.85 -2		2.23 -6.65	0.39 -2.38	5 0.52 -	4.85 -2	2.23 -2	0.39 -6.65	2.38	-	0.52
		1AY7 1MEL	185 1KTZ 2	5.73 282 1	1.82 5.39 7.89 2.84 0.45 1.25	126 0.63 0.52	7.85 735 2	2.46 5.39 1	0.45 1.25 7.89	0.63 2.84	0.52
		1BVN 1MLC	3 1KXP 3	4.09 1	1.74 7.43 15.41 3.04 0.50 2.17	2 0.51 0.24	4.09 4 3	1.74 7.43 1	0.50 2.17 15.41	0.51 3.04	0.24
		1CGI 1NCA	61 1ML0 2	3.20 1		2.30 4.47 7.62 2.13 0.49 1.89	37 0.61 0.66	3.20 1 2	2.30 4.47 1	0.49 1.89 7.62	0.61 2.13	0.66
		1D6R 1NMB	-1QA9 1	1	-	-		2.66	-	-	-0.57	-	-1.00 -		-1	-	-1	-	-2.66	-	0.57	-	1.00
		1DFJ 1PPE	1 1RLB 1	5.97 1		2.42 9.11 1.38 0.54 0.50 1.93	1 0.66 0.89		5.97 10 1	2.42 9.11 1	0.50 1.93 1.38	0.66 0.54	0.89
		1E6E 1QFU		9 1S1Q 1		4.01 766 1	2.41 2.35 3.89 0.98 0.42 1.53	8 0.58 0.67		3.87 1187 1	1.59 2.35 1	0.69 1.53 3.89	0.58 0.98	0.67
		1SPB		1		1				1.04			0.54			0.82		1	1	1.04	0.54	0.82

Γ n+q 2 -k, pz 2 k!(n -2k)! (2.19)

Though the scoring function (Eq. 3.3) is similar by the structure to e.g. the excess internal energy[START_REF] Hansen | Theory of simple liquids[END_REF], our scoring potentials U kl (r) are not equal to the potential energy functions between sites k and l.

Generally, if the distance distributions have a non-Gaussian shape, n kl (r) = ij f (r -r ij ), functions Υ kl (r) will be computed as a convolution Υ kl = f * U kl .

1NW9 ------------1GP2 --------

Acknowledgements

I would like to thank my supervisor Valentin Gordely and co-supervisor Sergey Grudinin for their guidance during my thesis. I also thank my parents, brother and all the people I was lucky to work with for their support.

Geometrical interpretation and connection to quadratic programming

Using the expansion of the scoring functional F provided by Eq. 3.13, we can reformulate the scoring problem (Eq. 3.1) as follows -given N native structure vectors x nat i and N × D nonnative structure vectors x nonnat ij , find a scoring vector w ∈ R P ×M ×(M +1)/2 such that:

or, equivalently,

which is a set of N ×D half-space equations in R P ×M ×(M +1)/2 . Each of the half-spaces is defined by a plane in R P ×M ×(M +1)/2 with the common normal w. Thus, finding the scoring vector is equivalent to finding the common normal w to the planes in Eq. 3.15. Geometrical representation of three groups of structure vectors separated by three parallel hyperplanes with the common normal w is given in Fig. 3.2. In the training set, some decoy structures can be very close to the native structures. In practice, we define the native structure as a structure with ligand rootmean-square deviation (lRMSD) smaller than 2 Å. Therefore, for each complex we may have several native structure vectors along with several nonnative structure vectors. Now the question is -how do we determine the set of separating hyperplanes shown in Fig. 3.2 with common normal w? To answer this question we first consider two special cases presented below.

I II III

Case I. Existence of many solutions

In Fig. 3.3A we present an example of a single complex when infinitely many hyperplanes can separate two classes of structure vectors. A similar example can be easily We reranked top 2000 decoys generated by ZDOCK3.0 using our scoring potentials. A hit is a predicted near-native decoy with iRMSD (RMSD of C α atoms of the predicted interaction interface residues after superposition onto the crystallized complex) less than 2.5 Å. The number of hits when only the top one prediction considered (Top1) obtained by ZRANK is higher than that obtained by ConvexPP Appendix A

Dual optimization problem

Optimization problem (Eq. 3.17) can be solved by the classical method of Lagrange multipliers [START_REF] Boyd | Convex optimization[END_REF][START_REF] Cortes | Support-vector networks[END_REF]. If we introduce N × D nonnegative Lagrange multipliers λ ij associated with the first set of inequality constraints from Eq. 3.17 and N × D nonnegative Lagrange multipliers ν ij associated with the second set of inequality constraints from (3.17), the solution of problem in Eq. 3.17 is equivalent to determining the saddle point of the following Lagrangian function: According to the classical Karush-Kuhn-Tucker (KKT) conditions [START_REF] Kuhn | Nonlinear programming[END_REF][START_REF] Boyd | Convex optimization[END_REF], which is a generalization of the method of Lagrange multipliers to inequality constraints, the saddle point of the Lagrangian function (Eq. A.1) satisfies four following conditions:

1. Stationarity conditions:

2. Complementary slackness conditions:

3. Primal feasibility conditions:

4. Dual feasibility conditions:

Using equation A.1 along with the aforementioned KKT conditions (Eq. A.2 -Eq. A.10), we can rewrite the original optimization problem (Eq. 3.17) as:

Sequential minimal optimization algorithm

Here we describe how the SMO algorithm [START_REF] Platt | Sequential minimal optimization: A fast algorithm for training support vector machines[END_REF] solves the problem in Eq. 3.18 for two Lagrange multipliers, λ 1 and λ 2 . All quantities that refer to the first multiplier have a subscript 1 and all quantities that refer to the second multiplier have a subscript 2. SMO first computes the constraints on these multipliers and then solves the problem (Eq. 3.18) for the constrained maximum. The inequality constraints in Eq. 3.18 force the two multipliers to lie within a box [0, C 1 ] × [0, C 2 ], while the equality constraints force the two multipliers to lie on a diagonal line segment:

This equation explains why we need to optimize two Lagrange multipliers simulteniously. Precisely, it is not possible to optimize a single multiplier without breaking the equality constraints in Eq. 3.18 and, subsequently, breaking the constraints (Eq. B.1). Without loss of generality, SMO first computes the second Lagrange multiplier λ 2 and then expresses the ends of the diagonal line segment in terms of λ 2 . The following lower and upper bounds, L 2 and H 2 , apply to λ 2 :

On the next step SMO computes the location of the unconstrained maximum of the Lagrangian with respect to λ 2 :

The corresponding unconstrained λ 2 will be: