Noël Malod-Dognin

Achin Bansal

Frédéric Cazals

Characterizing the Morphology of Protein Binding Patches

Keywords: Protein complex, binding patch, interface morphology, structural comparisons; Voronoi models, shelling tree, tree edit distance, dynamic programming based comparisons Complexes protéiques, patchs de liaisons, morphologie des interfaces, comparaisons de structures, modèles de Voronoi, arbre de couches, distance d

Let the patch of a partner in a protein complex be the collection of atoms accounting for the interaction. To improve our understanding of the structure-function relationship, we present a patch model decoupling the topological and geometric properties. While the geometry is classically encoded by the atomic positions, the topology is recorded in a graph encoding the relative position of concentric shells partitioning the interface atoms. The topological -geometric duality provides the basis of a generic dynamic programming based algorithm comparing patches at the shell level, which may favor topological or geometric features. On the biological side, we address four questions, using 249 co-crystallized heterodimers organized in biological families. First, we dissect the morphology of binding patches, and show that Nature enjoyed the topological and geometric degrees of freedom independently while retaining a nite set of qualitatively distinct topological signatures. Second, we argue that our shell-based comparison is eective to perform atomic-level comparisons, and show that topological similarity is a less stringent than geometric similarity. We also use the topological versus geometric duality to exhibit topo-rigid patches, whose topology (but not geometry) remains stable upon docking. Third, we use our comparison algorithms to infer specicity related information amidst a database of complexes. Finally, we exhibit a descriptor outperforming its contenders to predict the binding anities of the anity benchmark. The softwares developed with this paper are available from http://team.inria.fr/abs/vorpatch_compatch/.

Biology rests on macro-molecular complexes, so that a central question consists of understanding the determinants of the stability and the specicity of binding. These questions have been approached from two complementary perspectives, namely experiments and modeling. On the experimental side, structures resolved by X ray crystallography and NMR are of fundamental importance, as they lay the ground for modeling studies [START_REF] Janin | Protein-protein interaction and quaternary structure[END_REF], but also pave the way to protein engineering [START_REF] Fleishman1 | Computational design of proteins targeting the conserved stem region of inuenza hemagglutinin[END_REF]. Structural information is also complemented by directed mutagenesis and binding anity measurements, which convey information of biological and thermodynamical nature [START_REF] Fersht | Structure and Mechanism in Protein Science: A Guide to Enzyme Catalysis and Protein Folding[END_REF], and by evolutionary information [START_REF] Lichtarge | Evolution: a guide to perturb protein function and networks[END_REF]. Structural modeling work on the other hand, aims at developing explanatory and predictive models, and may be classied into two veins. To describe them, the following terminology is used to describe a binary complex: a binding patch (patch for short in the sequel) refers to a collection of atoms on one partner, responsible for the interaction; the union of two such patches denes the interface of the complex.

Dissecting the morphology of interfaces and patches. The design of explanatory and predictive models rests on the identication of structural parameters (geometric and topological) which best describe the biological and biophysical properties of interfaces [START_REF] Janin | Protein-protein interaction and quaternary structure[END_REF]. The rst task when studying an interface is to identify the atoms contributing to the two binding patches, as the buried surface area of these atoms in the complex often reliably hints at the stability of the interaction [START_REF] Lo Conte | The atomic structure of protein-protein recognition sites[END_REF]. To study the amino-acid composition and more generally the biochemical properties of patches, the core-rim model was introduced based on the accessibility of atoms in the complex [START_REF] Chakrabarti | Dissecting protein-protein recognition sites[END_REF]. This model was also used to show that conserved residues tend to locate in the core [START_REF] Guharoy | Conservation and relative importance of residues across protein-protein interfaces[END_REF]. In a more biological perspective, double mutant cycles were used to evidence the modular structure of binding patches [START_REF] Reichmann | From The Cover: The modular architecture of protein-protein binding interfaces[END_REF]. On the prediction side, algorithms computing putative patches on a molecular suface have been developed. Their strategy consists of generating patches according to specic model, and the putative patches are assessed against those observed in cocrystallized complexes. In [START_REF] Jones | Analysis of protein-protein interaction sites using surface patches[END_REF] and [START_REF] Albou | Dening and characterizing protein surface using alpha shapes[END_REF], the patch model consists of picking the k-nearest exposed neighbors of a central atom, resulting in disk-shaped (i.e. isotropic) patches.

Comparing patches. The design of structural comparisons and alignment tools has been carried out with two privileged applications. The rst one is the analysis of bound complexes and the classication of their interfaces. In [START_REF] Keskin | A new, structurally nonredundant, diverse data set of proteinprotein interfaces and its implications[END_REF], a classication of all interfaces of the PDB results into 103 classes split into three groups, based on the structural similarities inferred using geometric hashing, and the folds of the subunits. In a nearby vein, the SCOPPI database classies patches using successive criteria related to the SCOP domains, structural similarities, and sequence identity [START_REF] Winter | SCOPPI: a structural classication of protein-protein interfaces[END_REF]. These analyses, which are concerned with the specicity of biological interactions, ultimately aim at inferring whether two proteins interact, using information from interacting homologous proteins. The second one is the detection of similar binding patches on two unbound partners, a problem reminiscent from docking. While similarity detection is at the very heart of docking [START_REF] Halperin | Principles of docking: An overview of search algorithms and a guide to scoring functions[END_REF], or particular interest are the methods inferring similar patches from the proteins' exposed surface. In particular, in [START_REF] Konc | ProBiS algorithm for detection of structurally similar protein binding sites by local structural alignment[END_REF], similar patches are inferred by merging small graphs (graphlets) whose RM SD d is upper-bounded, and which exhibit comparable bio-chemical properties.

Predicting binding anities. The prediction of binding anities is also a critical endeavor, as reliable predictions would have a major impact on applications to protein engineering and the Inria study of interactomes [START_REF] Kastritis | Are scoring functions in protein-protein docking ready to predict interactomes? clues from a novel binding anity benchmark[END_REF]. The size of an interface, measured by the buried surface area ∆ASA upon complex formationor the number of interface atoms since both quantities are known to be proportional [START_REF] Chakrabarti | Dissecting protein-protein recognition sites[END_REF][START_REF] Bahadur | A dissection of specic and nonspecic protein-protein interfaces[END_REF], is known to correlate to the binding anity, at least for rigid docking cases [START_REF] Kastritis | A structure-based benchmark for protein-protein binding anity[END_REF]. Intuitively, the size hints at the complex stability [START_REF] Lo Conte | The atomic structure of protein-protein recognition sites[END_REF], as small interfaces may favor a fast turn-over of the interactions (e.g. in the enzyme -substrate case), while larger interfaces may favor more stable interactions (e.g. in the enzyme -inhibitor case). In addition to predictions based on such simple structural descriptors, a variety of scoring functions have been developed, their power to single out native-like complexes being assessed within the CAPRI experiment [START_REF] Lensink | Docking and scoring protein interactions: Capri 2009[END_REF][START_REF] Feliu | How dierent from random are docking predictions when ranked by scoring functions?[END_REF]. Their ability to estimate binding anity was also recently evaluated in the context of the binding anity benchmark [START_REF] Kastritis | Are scoring functions in protein-protein docking ready to predict interactomes? clues from a novel binding anity benchmark[END_REF][START_REF] Kastritis | A structure-based benchmark for protein-protein binding anity[END_REF].

Contributions

Questions addressed. As just recalled, a variety of geometric models have been proposed to investigate specic biological -biophysical properties of binding patches, and to compare binding patches. Following our analysis of protein interfaces, where enhanced geometric resulted in rened insights on biophysical properties [START_REF] Cazals | Revisiting the voronoi description of protein-protein interfaces[END_REF][START_REF] Bouvier | Shelling the voronoi interface of proteinprotein complexes reveals patterns of residue conservation, dynamics and composition[END_REF][START_REF] Cazals | Revisiting the Voronoi description of protein-protein interfaces: Algorithms[END_REF], we undertook the task of designing a hierarchical binding patch model amenable to a variety of structural studies. More specically, the patch model introduced herein was designed to address the following questions.

Dissecting the morphology of patches. As evidenced by previous work, the morphology of molecular shapes has essentially been described in geometric terms, using the atomic positions. On the other hand, in mathematics, topology is prosaically dened as the geometry of rubber made objects. In other words, topological features can be conserved, while geometric features are not. Our goal has therefore been to develop a patch model accommodating both the geometry and the topology of patches. Conceptually, enumerating the binding morphologies used by Nature is appealing. Application-wise, such morphologies are especially interesting wherever (isotropic) patch models are used, in particular to seek putative patches on orphan proteins.

Comparing patches. Taking for granted a patch model encoding both geometric and topological features, our second goal has been to design comparison algorithms operating at the atomic level, and privileging either features. Geometrically, the ability to compare patches at the atomic level is of interest to search similar patches, and in the context of structure-specicity studies, to assess if similar patches are involved in similar biological functions. Also, the possibility to assess geometric versus topological deformations undergone during docking is of interest to assess the diculty of docking problems. One indeed expects a rising diculty, from geometrically rigid partners to fully exible partners, with topologically rigid (but geometrically exible) cases as a third tier.

Identifying patches and specicity analysis. The specicity of protein interactions relates to the ability of a protein to bind selected partners, and a key problem in post-genomic studies consists of pulling back functional annotations from a given protein to homologous proteins. Conservation information can be used to do so, but such approaches are unable to handle the cases where similar patches are accounted for by dierent amino-acids. One of our goals has therefore been to investigate the ability of our patch model to detect a specicity related signal amidst a database of biological complexes organized into biological families.

Predicting binding anities. Finally, a key strength of a patch model lies in its ability to support binding anities estimations, at least in selected cases. We investigated this ability using the binding anity benchmark.

RR n°7743

Methodological positioning. Our patch model supports the questions just discussed thanks to its multi-scale structure. In a nutshell, we use the topological incidences of the spherical caps making up the solvent accessible surface of the interface atoms to dene an integer-valued depth of these atoms, which goes beyond the classical core-rim model. Next, we assign atoms of identical depth into so-called shells, and use the relative position of these shells to dene a graph coding the topology of the patchthe geometry of a shell pertaining to its associated atoms. This encoding is used to solve a relaxed structural alignment problem through a maximum clique calculationa NP-hard problem, as initially proposed in [START_REF] Brint | Algorithms for the identication of three-dimensional maximal common substructures[END_REF]. Assuming that two patches are given as lists of atoms BP 1 and BP 2 , we seek the largest subsets A 1 = {a i } i=1,...,n and A 2 = {b i } i=1,...,n of these lists, together with a one-to-one correspondence a i ↔ b i between them, such that the RMSD between internal distances is bounded by a user dened threshold ε, that is

RM SD d (BP 1 , BP 2 ) = 2 n × (n -1) i<j |d ai,aj -d bi,bj | 2 ≤ ε. (1)
This problem has long been known to be equivalent to the maximum clique problem in the so-called product graphs [START_REF] Brint | Algorithms for the identication of three-dimensional maximal common substructures[END_REF][START_REF] Malod-Dognin | Maximum clique in protein structure comparison[END_REF]. The maximum clique is a well known NP-Hard problem [START_REF] Karp | Reducibility among combinatorial problems[END_REF] that can be tackled with enumeration algorithms [START_REF] Bron | Algorithm 457: nding all cliques of an undirected graph[END_REF][START_REF] Cazals | A note on the problem of reporting maximal cliques[END_REF] or optimization algorithms [START_REF] Bomze | The maximum clique problem[END_REF][START_REF] Östergård | A fast algorithm for the maximum clique problem[END_REF][START_REF] Malod-Dognin | Maximum clique in protein structure comparison[END_REF].

Given that the size of the product graphs is quadratic in the number of atoms, and that typical patches involve from 100 to 500 atoms, the product graphs have a number of vertices beyond tens of thousands (up to 143385 vertices with our benchmark), which is intractable for exact algorithms. To fudge around this diculty, we use the aforementioned graph to localize the application of the maximal clique algorithm to shells. As we shall see, the problems solved still challenge the state-of-the-art maximum clique algorithms.

Methods

We rst present our encoding of patches as graphs, the associated comparison algorithms, and the application to database analysis.

A Hierarchical Encoding of Patches

Outline. Consider a binary complex, and assume that the interface atoms have been identied on both partners. Focusing on a given partner in the Solvent Accessible model, where the radii have been expanded by r w = 1.4, we dene its binding patch, called patch for the sake of conciseness, as the solvent accessible surface (SAS) of its interface atoms. A patch therefore consists of spherical polygons called faces; a face is bounded by circle arcs, and two faces are called incident if they share a circle arc; a circle arc is itself bounded by the points found at the intersection of three spheres. The peripheral atoms of such a patch make up its rim, and to measure the distance of a face to the patch boundary, each face is assigned an integer called its shelling order (SO). Finally, the SO are used to decompose the patch into concentric shells, whose relative positions are encoded in a graph called the face shelling tree, from which another graph called the atom shelling tree is derived. These notions are respectively illustrated and explained on Fig. 1 and Fig. SI-1. We now sketch the main steps of the atom shelling tree construction, and refer the reader to the supplemental section 6.1 for the details.

Inria

Dening patches. We identify the interface atoms with our Voronoi interface model [START_REF] Cazals | Revisiting the voronoi description of protein-protein interfaces[END_REF][START_REF] Bouvier | Shelling the voronoi interface of proteinprotein complexes reveals patterns of residue conservation, dynamics and composition[END_REF],

whose denition and software are presented in [START_REF] Loriot | Modeling macro-molecular interfaces with intervor[END_REF] and [START_REF] Cazals | Revisiting the Voronoi description of protein-protein interfaces: Algorithms[END_REF]. Let A and B be the two species of the complex, also called partners or subunits, and denote W the water molecules squeezed inbetween the partners. Let the restriction of a ball B i be the 3D region dened by the intersection between B i and its Voronoi region. (To be precise, the Voronoi region refers to the region of the ball in the power diagram of the balls of the SAS model.) Two atoms are called neighbors provided that their restrictions intersect. A water molecule is called interfacial provided that its has neighbors on both partners. An interface atom is an atom which is neighbor to the other partner's atoms, or to interfacial water molecules.

Having identied the interface atoms, we process the two subunits separately. The patch of a subunit is dened as the SAS of this partner restricted to its interface atoms. The patch is encoded in a so-called Half-edge Data Structure (HDS) [START_REF] Berg | Computational Geometry: Algorithms and Applications[END_REF], that gives access to the faces, the circle-arcs, and their incidences.

Face Shelling Tree. Shelling consists of assigning an integer value called shelling order or SO to each face, and is best presented in terms of graph distance. Term a face a boundary face if one of its bounding circle arcs is incident to one interface atom and one non-interface atom. Such faces are assigned a SO of zero. Consider now the dual graph of the patch: the nodes of this graph are the faces; two nodes are connected by an edge provided that the associated faces share a circle-arc. The shelling order or SO of a face is its shortest distance to a boundary face in the dual graph. Note that the contribution of an atom to the patch may consist of several faces with dierent SO. The SO and the dual graph are used to compute the following topological encoding. First, we dene a shell as a maximal connected component of the dual graph involving faces with the same shelling order, so that the patch is partitioned into shells. Second, we encode the relative position of shells within the so-called face shelling tree. This tree contains one node N SG (s) for each shell s, the node size being the number of faces. To see how the edges of the shelling tree are dened, consider two incident faces whose SO dier of one unit. Let s and t be the shells containing these faces, and assume that SO(s) + 1 = SO(t) : the face shelling tree contains one arc from N SG (s) to N SG (t). The number of outgoing arcs of a node is called its arity.

Atom Shelling Tree. In order to base the comparison of patches on atoms rather than faces, we edit the face shelling tree into an atom shelling tree. The process consists of substituting atoms to faces, with the following special cases: if an atom is present several times in the same shell, it is counted once; if an atom belongs to several shells in a branch of the face shelling tree, it is assigned to the shell closest to the root of the tree. Finally, the sons of a node are sorted by increasing size i.e. number of atoms, resulting in an ordered atom shelling tree, called shelling tree for short in the sequel.

Note that the atom shelling tree encodes topological information namely the relative position of the shells, while the 3D coordinates of the atoms within the shells encode the geometry.

Comparing Patches: a Generic Dynamic Programming based Approach

Encoding a patch as an ordered tree whose nodes contain shells paves the way to patch comparison using dynamic programming [START_REF] Bille | A survey on tree edit distance and related problems[END_REF]. More precisely, to compare two trees, we edit one into the other, computing the so-called Tree Edit Distance (TED). The TED, whose details are recalled in the supplemental section 6.2.1, is based on three operations, namely node deletion, node insertion, and node morphing. The TED calculation delivers an Ordered Edit Distance Mapping, namely

RR n°7743

a set M ⊂ V ertices(T 1 ) × V ertices(T 2 ) such that for any pair (v 1 , v 2 ) ∈ M and (w 1 , w 2 ) ∈ M , one has: (i) v 1 = w 1 i v 2 = w 2 , or (ii) v 1 is an ancestor of w 1 i v 2 is an ancestor of w 2 , or (iii) or v 1 is to the left of w 1 i v 2 is to the left of w 2 . (Recall that trees are ordered.)

In our case, given that a node corresponds to a shell of atoms, adjusting the substitution cost yields two strategies to compare the topology and the geometry of patches, respectively. Topological comparison. To identify common patterns of nested shells encoded within the shelling tree, we compute the TED with the following costs. Adding or deleting a node associated with a shell s has a cost of |s|, namely the number of atoms in the shell. Morphing a shell s 1 into a shell s 2 corresponds to matching min(| s 1 |, | s 2 |) atoms in-between the two shells, or equivalently, to a cost max(|

s 1 |, | s 2 |) -min(| s 1 |, | s 2 |).
At the patch level, the atoms matched by the TED calculation are called isotopologic since they belong to nodes of the shelling trees satisfying the constraints (i,ii,iii) of the edit distance mapping. Denoting SIM t (T 1 , T 2 ) the number of isotopologic atoms between the patches, the TED cost is the following symmetric dierence

TED t (T 1 , T 2 ) =| T 1 | + | T 2 | -2 SIM t (T 1 , T 2 ). (2)
This number being upper-bounded by |T 1 | + |T 2 |, it yields the dissimilarity score ∈ [0, 1]:

DIS t (T 1 , T 2 ) = TED t (T 1 , T 2 )/(| T 1 | + | T 2 |), (3)
which can be interpreted as the percentage of non-common atoms.

Geometric comparison. Consider now the problem of comparing the geometry of two patches, as specied by Eq. ( 1). Because a brute-force attempt to solve this problem for the whole patch is intractable, we restrict the identication of quasi-isometric subsets to pairs of shells. That is, we dene a second TED calculation as follows.

As previously, the cost of inserting/deleting a shell s is | s |. For the morphing cost between shells s 1 and s 2 , assume that | s 1 s 2 | quasi-isometric atoms have been identied by a maximum clique calculation, as specied by Eq. (1) (see details in the supplemental section 6.2.2). The morphing cost is equal to the size of their symmetric dierence, namely

| s 1 | + | s 2 | -2 | s 1 s 2 |. Denote SIM g (T 1 , T 2 )
the number of atoms matched across all pairs of nodes in the edit distance mapping. The corresponding tree edit distance counts the number of un-matched atoms, namely:

TED g (T 1 , T 2 ) =| T 1 | + | T 2 | -2 SIM g (T 1 , T 2 ). (4)
Mimicking Eq. (3), we dene:

DIS g (T 1 , T 2 ) = TED g (T 1 , T 2 )/(| T 1 | + | T 2 |) (5)
Topology versus geometry. The previous two criteria both rely on a TED calculation, and report isotopologic atoms. Yet, the geometric comparison is more stringent, and SIM g (T 1 , T 2 ) ≤ SIM t (T 1 , T 2 ). Equivalently, the topological dissimilarity is a lower bound of the geometric dissimilarity, that is DIS g (T 1 , T 2 ) ≥ DIS t (T 1 , T 2 ).

Identifying Patches and Specicity Analysis

A dataset of n co-crystallized protein complexes yields a database P of 2n patches. We assume that the database is organized into biological families corresponding to biological functions. We further split each family by distinguishing the ligand and the receptor of each complex. Thus, Inria the database of patches is partitioned into typed families. This decomposition scheme aims at performing structural comparisons in conjunction with the analysis of biological functions. Prosaically, we wish to investigate whether it makes sense to speak of the patches of say immunoglobulin -peptides complexes. Denoting P the patches of a typed family, let P be the set of patches that are the partners of the ones in P , and P c be the set of patches neither in P nor in P . Note that P = P ∪ P ∪ P c . For a given patch p, dene P \p such that P = {p} ∪ P \p . We shall use the following partition of the database induced by any patch to test hypothesis on the similarity of patches:

P = p ∪ P \p ∪ P ∪ P c . ( 6
)
Given a patch p and a dissimilarity score s(p, q), practically that of Eq. (3) or Eq. ( 5), we denote by p the nearest neighbor i.e. the patch of the database with lowest dissimilarity:

p = arg min q∈P\{p} s(p, q). (7)
We distinguish the following cases: case I: p ∈ P \p , case II: p ∈ P , and case III: p ∈ P c . Case I directly measures the compatibility between the dissimilarity score and the typed family classication, and the possibility of using the dissimilarity score for generating automatic classications of patches. It is also related to the morphology -specicity relationship, since the proteins in this case meet two criteria, namely they bind the same family of ligand, and their binding patches possess similar morphologies. Case II is related to the symmetry (or lack of) between partner patches across an interface, a feature especially interesting for heterodimers homodimers are symmetrical, albeit not always exactly. This information is of particular interest for docking, where the steric complementarity between the partners is often used as a matching criterion. Case III highlights contradictions between the typed family classication and the dissimilarity score values.

3 Results

Database

The results presented in this section were obtained on a database of 498 patches generated from 249 complexes. These complexes are heterodimers and one of our goals is to study the symmetry of patches. The set of complexes was assembled from two sources. The rst one is the IMGT 3D structure database [START_REF] Ehrenmann | Imgt/3dstructure-db and imgt/domaingapalign: a database and a tool for immunoglobulins or antibodies, t cell receptors, mhc, igsf and mhcsf[END_REF], from which we extracted 116 immunoglobulin -ligand structure, with resolution ≤ 2.0. The second one is the recently assembled binding anity benchmark [START_REF] Kastritis | A structure-based benchmark for protein-protein binding anity[END_REF], a manually curated dataset involving 133 complexes with experimentally measured binding anity, and resolutions in the range 1.1 -3.3, the median being 2.4. (In fact, the anity benchmark contains 144 complexes, but 11 were redundant with the ones that we extracted from the IMGT 3D database.) The reader is referred to the Tables SI-1 and SI-2 for the presentation of the typed families and for the exhaustive list of PDB ids.

The rationale in assembling this database has been the following. First, all the complexes are known to be biological complexesas opposed to crystallization artifacts. Second, the database strikes a balance in terms of diversity and homogeneity: on the one hand, it involves a large variety of complexes and biological functions, in particular the so-called O set of the anity benchmark; on the other hand, the subset from IMGT 3D gives the opportunity to compare a signicant number of proteins involved in the same biological function. Third, the anity RR n°7743 benchmark allows assessing the geometrical and topological changes undergone by the patches upon binding, since the bound and unbound forms of each partner are known.

We used this set of complexes for two purposes. On the one hand, the morphological study of section 3.2 was carried out on the whole database. On the other hand, the clustering and identication study of section 3.4 was performed on the subset consisting of complexes involving enzymes or immunoglobulins. (The variety of complexes involved in the aforementioned O set precludes the identication of patches with signicant similarity.)

All the patches but two (496 over 498) possess between 14 to 294 atoms and 1 to 20 shells. The two exceptions are the patches of the signaling complex 2oza, since chain B contributes 395 atoms and 28 shells (Fig. SI-4), while chain A contributes 365 atoms and 47 shells (Fig. SI-5). Also, the 498 patches yield a total of 498 2 + 498 = 124251 pairwise comparisons.

Dissecting the Morphology of Patches

Canonical Morphologies. To analyze the morphology of patches beyond the core-rim model, we rst assess the repartition of atoms in a patch by plotting the number of atoms against the number of shells (Fig. 2). Since this plot exhibits a continuous variation, we extract typical morphologies by examining extreme cases for a xed number of shells and atoms, respectively.

Minimizing and maximizing the number of atoms for a xed number of shells yields tubular and pyramidal shapes (Fig. 3). Similarly, minimizing and maximizing the number of shells for a xed number of atoms yields anisotropic and isotropic shapes, respectively.

To understand the specicity of these shapes, we plot for each patch the variation of the number of atoms as a function of the SO. Inspection of all curves (data not shown) allows us to single out ve cases (Fig. 4). A curve contained in a narrow horizontal slab corresponds to a tubular shape. A curve involving an increasing section followed by a decreasing section corresponds to a pear-like shape. In the remaining cases, the curve is decreasing. The maximum SO may be used to distinguish isotropic (large max SO) from anisotropic (small max SO) patches. As illustrated on Fig. 3, in each case, the geometry can be used to dene at versus non-at (pyramidal) patches.

Asymmetry Between Partner Patches. Having singled out these shapes, we compare the two patches of a complex resorting to the average shelling order or SO of a patch, dened as the sum of the shelling order of each atom of the patch divided by the number of atoms. This quantity is maximized for a linear atom shelling tree implying that most of the atoms are deep inside the patch, and is minimized for at trees, implying that most of the atoms are located around the patch periphery. For a given complex, the asymmetry of its patches are witnessed by dierent SO, and the two families AA_Pept and the AA_Prot appear as very asymmetric (Figs. [START_REF] Lo Conte | The atomic structure of protein-protein recognition sites[END_REF] and Fig. 6).

Comparing Patches

Ecacy of the Atom Shelling Tree Encoding

Given the hardness of the geometric matching problem, as specied by Eq. ( 1), we rst report the running times of the topological and geometric comparisons, whose punchline consists of localizing calculations at the shell level.

On a 100 nodes cluster equipped with Intel Xeon processors at 2.66Ghz, the 124251 pairwise comparisons of our database took 622 seconds with TED t , 41901 seconds with TED g ( = 1Å), Inria and 1076390 seconds with TED g ( = 2Å). For the latter, a time limit of 2 hours per instances was used, and 13 instances remained unsolved.

To compare the topological and geometric matchings, we plot DIS g against DIS t . While it has been observed in section 2.2 that DIS g ≥ DIS t , Fig. 7 illustrates the variation DIS g as a function of DIS t . Of particular interest are pairs of patches with a similar topology (low DIS t ) but dierent geometries (high DIS g ). Two such patches, having respectively 82% and 32% of common atoms from the topological and geometric standpoints (Fig. SI-6).

Finally, we assess the accuracy of our geometric matchings, an important issue since algorithm TED g provides quasi-isometric matchings at the shell level, but does not provide any guarantee on the RM SD d at the patch level. The RM SD d globally increases with the geometric dissimilarity, and so does the variance of the RM SD d for a xed DIS g (Fig. SI-7). Interestingly, the RM SD d remains moderate (less than 4Å) for a geometric dissimilarity smaller or equal to 0.25 (i.e. between patches having at least 75% of geometrically common atoms), showing that seeking quasi-isometric subsets at the shell level is eective to compare whole patches. Note that from now on, we assume that two patches are geometrically (resp. topologically) similar if their DIS g (resp. DIS t ) is less or equal to 0.25.

Patches and their Preimages

Computing preimages. The topological similarity being less stringent than the geometric one, we use it to scale the conformational changes undergone by the partners of the binding anity benchmark. To this end, consider a partner of a complex, in its bound and unbound forms: the patch being dened from the bound partner, we dene the corresponding atoms on the unbound partner as the patch preimage.

Preimages were generated for the complexes coming from the full Anity benchmark (144 complexes) as follows. First, a one-to-one chain mapping between the chains of the bound and unbound partners was sought, and retained in case of sequence identity higher than 90% for any two chains put in correspondence. (Alignments were generated exactly i.e. without any substitution matrix.) These criteria dismissed 43 complexes. Second, for the 101 valid complexes, we produced atom mappings between the residues matched by the sequence alignments. Then, we generated the preimage, requiring the presence on the preimage of at least 90% of the patch's atoms. Over the 202 patches, only 126 met this requirement.

On docking diculties. For each pair (patch, preimage), we plotted the topological dissimilarity against the geometrical dissimilarity (Fig. 8). While this plot exhibits a continuous distribution and a monotonic correlation (Spearman coecient of 0.87), three cases can be singled out (Fig. 9). Two extreme cases correspond to rigid and exible docking, where both the topology and the geometry hardly change and signicantly change, respectively. The intermediate situation arises when the patch and its preimage have a similar topology but a dierent geometry: such a patch is called topo-rigid. In contrast to our analysis, with an interface RMSD (I-RMSD) of 0.17Å, 0.48Å, and 0.35Å respectively, the rigid, topo-rigid and exible patches appear as rigid from the alpha carbon point of view [START_REF] Kastritis | A structure-based benchmark for protein-protein binding anity[END_REF].

On morphologies and docking predictions. In a second experiment, we tested whether similar preimages lead to similar patches upon binding (Fig. 10(A)). For any two receptor and any two ligand patches, we plotted the topological dissimilarity between preimages against the topological dissimilarity of the resulting patches (Fig. 10(B)), and the geometrical dissimilarity between preimages against the geometrical dissimilarity of the resulting patches (Fig. 10(C)). The corresponding Pearson correlation coecients are respectively of 0.69 and 0.91, with a p-value RR n°7743 smaller than 1e -99 in both cases. Interestingly, topologically similar preimages (DIS t ≤ 0.25) yield topologically similar patches in 88% of the cases (over the 1238 pairs of similar preimages, 1093 lead to similar patches). On the other hand, we could not verify whether geometrically similar preimages produce geometrically similar patches, due to the paucity of geometrically similar preimages in our dataset (see Computing preimages, above).

Identifying Patches and Specicity Analysis

The identication results, as specied in section 2.3, are summarized in Table 1 for TED t and TED g , with dierent thresholds for the latter parameter. For all methods, a large discrepancy is observed between cases I and II, further stressing the lack of symmetry of patches. Moreover, immunoglobulin and enzyme typed-families exhibit extreme behaviors (Table 1 and Table SI- [START_REF] Fersht | Structure and Mechanism in Protein Science: A Guide to Enzyme Catalysis and Protein Folding[END_REF]. About 70% of the immunoglobulin patches have their nearest neighbor coming from their own biological family, showing that proteins binding the same type of ligand possess a similar morphology.

To further understand the dierent identication rates of typed families, we proceed as follows. First, we compute the Wilcoxon-Mann-Whitney rank-sum test of s(P, P ) against s(P, P c ), also known as the U test, the null hypothesis being that the two series of scores come from the same distribution. Second, we compute the Spearman's correlation coecient between the negative log of the U-test and the identication rate, as this test aims at detecting monotonic correlations. The correlation coecient obtained, equal to -0.893, shows the strong coupling between the family consistency (measured by the U test) and the identication rate. In other words, on consistent families, our comparison achieves good identication rate (Fig. SI-8).

Predicting Binding Anities

Binding anities versus structural parameters. The binding anity of a protein complex can be described by its equilibrium dissociation constant K d , which is related to the dissociation Gibbs free energy by ∆G = -RT ln K d /c • (in the c • = 1 M standard state). We thus investigate the correlations betweenln K d and our descriptors. In addition to the number of shells and the depth of a binding patch, we dene its internal path length as IPL = #atoms × SO. Equivalently, the IPL is the sum for all the atoms of their depth in the atom shelling tree, the root being at depth zero. From a statistical standpoint, we computed three correlations coecients: Pearson's coecient C Pea which aims at detecting linear correlations, Spearman's coecient C Spe which is better suited to detect non ane monotonic correlations, and the recently designed Maximal Information Coecient C MIC [34]. The latter remarkably targets general functional associations, and also behaves consistently across functional associations plagued with the same noise level. P-values based on permutation tests were also computed. (Most of the p-values of the MIC coecient are missing, though, since the tables delivered with the Supporting Information of [START_REF] Reshef | Detecting novel associations in large data sets[END_REF] do not feature p-values for small correlations: e.g., no value is provided for C MIC ≤ 0.27 in the table of size n = 140.)

We computed these correlations at the database level (Table SI-5), and also for the three classes of low, moderate and high exibility introduced in [START_REF] Kastritis | A structure-based benchmark for protein-protein binding anity[END_REF] (Tables SI-6 and SI-7). The IPL consistently yields the best correlation, both in terms of correlation coecient and p-value, the contenders being the number of atoms and the interface depth, and to a lesser extent the interface area ∆ASA. For example, C Spe (-ln K d , IPL) = -0.43 and C Spe (-ln K d , #atoms) = -0.37, while C MIC (-ln K d , IPL) = 0.35 and C MIC (-ln K d , #atoms) = 0.24. If the correlations coef- cients are moderate, their signicance is strong. Correlations computed on classes of varying exibility yield better results, in particular for the class I-RMSD ≤ 1. With C Spe (-ln K d , IPL) = Inria -0.59 versus C Spe (-ln K d , #atoms) = -0.58, and C MIC (-ln K d , IPL) = -0.48 versus C MIC (-ln K d , #atoms) = -0.43, the IPL still provides the best correlation.

Binding anities versus structural changes. A question often faced consists of understanding why related complexes exhibit very dierent anities. This issue is discussed in detail for nine pairs of complexes in [START_REF] Kastritis | A structure-based benchmark for protein-protein binding anity[END_REF], where it is noted that the two complexes of a pair have a similar geometry. Using our encoding of patches, we computed the geometric dissimilarity of the receptor and ligand patches (Table SI-8).

Except for complexes 2ptc__E-I and 2tgp__Z-I, the fact that DIS g ≥ 0.35 shows that the receptors are actually geometrically dierent, and likewise for the ligands, despite seemingly modest changes in amino-acid sequence. Indeed, in the case of 2vir__AB-C vs 2vis__AB-C, only one residue substitution in the receptor side occurs, leading to geometric dissimilarities of 0.35 and 0.23 between the receptors' and ligands' patches, respectively; for 1efn__B-A vs 1avz__B-C, a one residue substitution in the ligand side results in geometric dissimilarities of 0.45 and 0.48 between the receptors' and ligands' patches, respectively.

Discussion

Dissecting the morphology of patches

The question of understanding which features of patches account for the specicity of biological interactions has been examined in two veins. On the one hand, a number of works performed correlation studies between structural parameters (interface size, planarity, modularity, organization into a core and a rim), and biophysical properties (composition, conservation, solvation, ∆∆G) [START_REF] Lo Conte | The atomic structure of protein-protein recognition sites[END_REF][START_REF] Chakrabarti | Dissecting protein-protein recognition sites[END_REF][START_REF] Bahadur | Dissecting subunit interfaces in homodimeric proteins[END_REF][START_REF] Bahadur | A dissection of specic and nonspecic protein-protein interfaces[END_REF][START_REF] Guharoy | Conservation and relative importance of residues across protein-protein interfaces[END_REF]. Selected such parameters have also been traced along molecular dynamics [START_REF] Mihalek | On itinerant water molecules and detectability of proteinprotein int erfaces through comparative analysis of homologues[END_REF]. While trends have emerged for collections of complexes [START_REF] Mihalek | On itinerant water molecules and detectability of proteinprotein int erfaces through comparative analysis of homologues[END_REF][START_REF] Guharoy | Conservation and relative importance of residues across protein-protein interfaces[END_REF][START_REF] Janin | Protein-protein interaction and quaternary structure[END_REF], conclusions from meta-analysis, when rened under the lens of sharper structural parameters, may not apply to isolated complexes [START_REF] Bouvier | Shelling the voronoi interface of proteinprotein complexes reveals patterns of residue conservation, dynamics and composition[END_REF]. On the other hand, the search of binding sites on orphan proteins motivated the development of patch models, which when instantiated on a protein surface, allow the comparison between these instantiations and the patches observed in a co-crystallized complex [START_REF] Jones | Principles of protein-protein interactions[END_REF][START_REF] Jones | Analysis of protein-protein interaction sites using surface patches[END_REF][START_REF] Albou | Dening and characterizing protein surface using alpha shapes[END_REF]. So far, the patch models proposed are isotropic ones, as they consist of tracing a geodesic disks on the molecular surfaces.

In this context, this work elaborates on our Voronoi interface model [START_REF] Cazals | Revisiting the voronoi description of protein-protein interfaces[END_REF][START_REF] Bouvier | Shelling the voronoi interface of proteinprotein complexes reveals patterns of residue conservation, dynamics and composition[END_REF][START_REF] Loriot | Modeling macro-molecular interfaces with intervor[END_REF], which oers a unied way to rene classical interface parameters [START_REF] Janin | Protein-protein interaction and quaternary structure[END_REF], and proved instrumental to transpose conclusions from the database level to the single complex level, regarding the biochemical properties of interfaces [START_REF] Chakrabarti | Dissecting protein-protein recognition sites[END_REF], the geometry of conservation [START_REF] Guharoy | Conservation and relative importance of residues across protein-protein interfaces[END_REF], as well as the solvation of residues along molecular dynamics trajectories [START_REF] Mihalek | On itinerant water molecules and detectability of proteinprotein int erfaces through comparative analysis of homologues[END_REF].

The atom shelling tree construction extends the shelling of Voronoi interfaces, and brings novelties in the two directions. With respect to correlation studies, the atom shelling tree is a hierarchical encoding of the patch, replacing a binary attribute (location of an atom in the rim or the core) by an integer-valued one (the atom shelling order). The shelling tree makes it possible to study the topology of a patcha dimension ignored so far, independently from its geometry. While a continuous distribution of patches is observed with respect to topological and geometric features, we have shown that typical patch morphologies, namely tubular, pear-like, pyramidal, isotropic and anisotropic, could be singled out. These typical morphologies show that the isotropic patch models used so far do not account for the variety of morphologies encountered. Phrased dierently, using isotropic patch models is a hindrance to identify potential binding patches. Our encoding also allows assessing the symmetry of the two patches of a complex,

RR n°7743

either directly from the dissimilarity score, or indirectly based on statistics derived from the atom shelling treethe average shelling order SO. In particular, we believe that the ability to evidence the lack of symmetry for selected biological families will shed new light on the problem of determining surface complementarity for rigid body docking.

Three developments should prove particularly useful. The rst one is related to the inverse problem. The atom shelling trees are generated from the observed patches of complexes, but the number biological complexes known is very limited compared to the number of known unbound protein structures and crystal contacts. The inverse problem would consist in nding the optimal mapping of a given atom shelling tree onto a given protein surface. Because, as illustrated by our canonical morphologies, a patch cannot be reduced to an isotropic shape, a solution to this problem would enhance the search of putative patches on orphan proteins, in the spirit of [START_REF] Jones | Analysis of protein-protein interaction sites using surface patches[END_REF][START_REF] Albou | Dening and characterizing protein surface using alpha shapes[END_REF]. The second one is related to the asymmetry detection and to partner retrieval. The asymmetry reported in this work is related to the non-atness of interfaces. This hinders the possibility to retrieve the possible partners of a given patch, in particular for docking applications. Technically, the asymmetry detection comes from the fact that our generic dynamic programming based matching algorithms produces and Ordered Edit Distance Mapping, which is a one-to-one atom mapping. Developing a more general, say k-to-k atom mapping would allow adjusting the level of non-symmetry tolerated as a function of the parameter k. This extension poses challenging graphtheoretical problems, but would prove useful for docking, in combination with lters avoiding steric clashes and forcing the bio-chemical compatibility between the atoms matched.

Comparing patches: geometry versus topology

The problem of comparing and clustering interfaces and patches motivated work in two directions. On the one hand, approaches have been developed for co-crystallized complexes. In [START_REF] Keskin | A new, structurally nonredundant, diverse data set of proteinprotein interfaces and its implications[END_REF][START_REF] Keskin | Similar binding sites and dierent partners: Implications to shared proteins in cellular pathways[END_REF], interfaces are clustered both from the patch and the whole structure point of view, using geometric hashing techniques applied to the C α carbons. Identical interfaces are assigned to so-called class I or class II clusters depending on whether they involve chains with similar or dierent fold, while class III clusters regroup similar patchesthe patches are similar but the interfaces are not. In a nearby vein, the SCOPPI database [START_REF] Winter | SCOPPI: a structural classication of protein-protein interfaces[END_REF] classies patches using a two step approach. First, domain sequences from the same SCOP family are aligned, and then all patches are mapped over their aligned sequence. This 0-1 vector called Interface Tag (IFT) is used as the signature of the patch0 codes a non-interface residue while 1 codes an interface residue. Patches are then clustered into the same family if the cosine angle distance between their two IFT is larger than 0.8. (We note in passing that the IFT comparison does not take into account the gaps induced by the multiple sequence alignments, and thus does not convey any information on the coverage of interface atoms, as opposed to our dissimilarity score.) Second, the obtained families are clustered using selected geometric criteria. On the other hand, tools have been developed to compare solvent accessible patches. In Probis [START_REF] Konc | ProBiS algorithm for detection of structurally similar protein binding sites by local structural alignment[END_REF], graphlets encoding the proximity of functional groups are rst dened. Selected graphlets are compared using a maximum clique approach, and the global match between two patches is obtained by merging graphlets.

Our matching algorithms depart from these works in two major ways. First, we accommodate independently topological and geometric comparisons, based on the atom shelling tree encoding. We have seen that the former comparison if more lenient than the latter, and that both benets from the atom shelling tree encoding to perform the comparison at the shell rather than whole patch level. Moreover, algorithm TED g is accurate to perform atomic scale comparisons, since a geometric dissimilarity DIS g < 0.25 yields RM SD_d < 4.

Inria

Second, the ability to perform the matching at the atomic level, as opposed to the residue and functional group levels in [START_REF] Keskin | A new, structurally nonredundant, diverse data set of proteinprotein interfaces and its implications[END_REF][START_REF] Winter | SCOPPI: a structural classication of protein-protein interfaces[END_REF], motivated the study of patches and their preimages on the binding anity benchmark. We have seen in particular that cases which exhibit comparable diculty level in terms of interface RMSD actually correspond to dierent cases, namely rigid, exible, and topo-rigid. This latter case is actually somewhat expected, as in the conformer selection model [START_REF] Grünberg | Flexibility and conformational entropy in proteinprotein binding[END_REF], the presence of a topologically pre-formed patch may ease the formation of the encounter complex, and thus of the nal complex. Knowing that a patch is topo-rigid is also of interest for sampling the conformational space, as the extra constraints imposed by the stable topology may enable the sampling of a sub-space of the whole conformational space of that partner. In the same spirit, our topological comparison should also prove useful to assess docking results, as done in the CAPRI experiment, to complement the purely geometric criteria currently used [START_REF] Lensink | Docking and scoring protein interactions: Capri 2009[END_REF]. Another striking point outlined by our analysis is the stability of the topological signature of the patches upon docking, since we observed that 88% is similar preimages are resulting in similar bound patches. This information is clearly of high interest for template-based docking, to pull-back information on a putative bound partner, from the bound structure of a similar patch.

Identifying patches and specicity analysis

The coherence analysis of typed families, as revealed by the identication rates and the clustering properties, shows that the atomic level geometric and topological comparisons provide stringent criteria. In particular, algorithm TED t being more lenient than algorithm TED g , the topological consistency is better than the geometric one. However, diverse typed families such as the so-called O set of the anity benchmark currently appear as very diverse. The identication rates also shed light on the morphology -specicity relationship, since they highlight the ability of patches of a family of receptors to bind a specic ligand. The fact that 70% of the immunoglobulin patches have their nearest neighbor coming from their own biological family (Table 1) is interesting in two respects. First, this shows that immunoglobulin patches that binds the same kind of ligand possess similar patch morphologies (i.e. low dissimilarity scores), while immunoglobulin patches that binds dierent kinds of ligand (e.g. a peptide versus a small chemical) possess dierent morphologies (i.e high dissimilarity scores). Second, the characteristics of immunoglobulin patches are dierent from those of the 122 enzyme-related patches coming from the Anity benchmark. These specicities call for further developments. In the clustering of interfaces and patches of [START_REF] Keskin | A new, structurally nonredundant, diverse data set of proteinprotein interfaces and its implications[END_REF][START_REF] Keskin | Similar binding sites and dierent partners: Implications to shared proteins in cellular pathways[END_REF], class III clusters gather geometrically similar patches involved in dissimilar interfaces, i.e. patches having more than one binding function. In a related vein, the specicity of interactions has been studied in [START_REF] Panjkovich | Predicting proteinprotein interaction specicity through the integration of three-dimensional structural information and the evolutionary record of protein domains[END_REF] using 3D templates pre-processed into a database of interaction types based on PFAM domains, together with information on the conservation of residues. As patches can be similar at the residue level but dissimilar at the atomic level, our comparison tools should be useful to rene such analysis. More generally, our ability to handle coherently topological and geometric criteria calls for the development of hierarchical classication of patches, based on a combination of topological, geometric and biological (sequence) information, in a manner similar the classication of quaternary structures performed in [START_REF] Levy | 3D complex: a structural classication of protein complexes[END_REF]. A solution to the aforementioned inverse problem would allow performing such studies not only on the biological complexes of the PDB, but also on isolated proteins of known structure. In a related spirit, we envision applications of such analysis to patches from the immune system, in the context of the IMGT_3D database, see http://www.imgt.org/. Such studies would be particularly meaningful in the context of the collier de perles annotations [START_REF] Ruiz | Imgt gene identication and colliers de perles of human immunoglobulins with known 3d structures[END_REF], which assigns a unique numbering to the amino-acids of the complementarity determining regions (CDR). The CDR are responsible for the specicity of the RR n°7743 immune response in general and for anity maturation in particular. Structural modeling of this latter phenomenon is especially challenging since conservation information is irrelevant, and a unique mutation may result in a signicant change of the anity [START_REF] Lippow | Computational design of antibody-anity improvement beyond in vivo maturation[END_REF]. We foresee that our patch model will prove useful both to assess the geometric changes associated to mutationsas reported herein for seemingly geometrically similar complexes of the anity benchmark, and also to understand with canonical morphologies account for specic interactions for a given family of antigens.

Predicting binding anities

Dissociation constants of biological complexes are known to span the range 10 -5 < K d < 10 -14 , and are also known to depend on factors such as the temperature, the ionic strength, molecular crowding, and in particular the pH [START_REF] Kastritis | A structure-based benchmark for protein-protein binding anity[END_REF]. Because the free energy change ∆G has enthalpic and entropic terms, estimating the binding anity from features of the bound complex is doomed to fail in the most general setting, as the entropic changes cannot be estimatedletting alone the entropic changes of the solvent, and complexes formed by (semi-)rigid docking seem more accessible [START_REF] Kastritis | A structure-based benchmark for protein-protein binding anity[END_REF].

Recently, the Pearson correlation coecients betweenln K d and various scores were reported in [START_REF] Kastritis | Are scoring functions in protein-protein docking ready to predict interactomes? clues from a novel binding anity benchmark[END_REF], on a database of 81 protein complexes. The best correlation obtained is that of the FIREDOCK score [START_REF] Mashiach | Firedock: a web server for fast interaction renement in molecular docking[END_REF], with a value of C Pea = -0.32 and a p-value in the range 0.01 ≤ p-val ≤ 0.05. On the whole anity benchmark, we observe that the same correlation drops down to C Pea = -0.17, a value much worse than C Pea = 0.31 observed with our internal path length (Table SI-5).

While the latter correlation remains mediocre, the ability of our purely geometric and topological descriptor to outperform the complex top-scoring function of [START_REF] Kastritis | Are scoring functions in protein-protein docking ready to predict interactomes? clues from a novel binding anity benchmark[END_REF] is interesting in several respects. On the one hand, the IPL can be seen as a weighted version of the number of interface atoms, the weight of an atom being is distance from the rim of the binding patch. Given the expected correlation between the interface size (or the number of atoms) and the anity, the gain provided by the encoding of the depth should not come as a surprise. The corresponding plots (Fig. SI-9) also show that minimum values of the parameters are required to reach a given high anity. For example, for K d ≤ 10 -10 or equivalentlyln(K d ) ≥ 23.02, one observes that the minimum average number of atoms is 110, that the minimum average depth is 3 (i.e 4 concentric shells), that the minimum IPL is 160.5, and that ∆ASA ≥ 1684 2 . In other words, the typical values of the parameters required to match the prescription K d ≤ 10 -10 correspond to the typical average values observed for protein-protein complexes [1, Table 2].

On more speculative grounds, the depth of atoms within a binding patch is likely related to dynamic properties, both of the partners and of the solvent. For the latter, it has also been observed in [START_REF] Bouvier | Shelling the voronoi interface of proteinprotein complexes reveals patterns of residue conservation, dynamics and composition[END_REF] that the Voronoi shelling order of an interface has a strong correlation with the dynamics of water molecules squeezed in-between the partners. For the partners, it is appealing to believe that the increase of internal entropy observed upon binding is related to the shape of the interface coded by the two shelling trees and their matching, since this matching precisely describes the coupling between the two sub-units. Developing quantitative models exploiting our encoding of depth might therefore prove fruitful for the prediction of dynamic properties and the improvement of binding anity predictions, at least for selected classes of complexes. Table 1 Identication rates via the typed family of the most similar patch. For each method, columns 2 to 4 present, in percentage, the number times p comes from the family of p (column 2), from the partner family of p (column 3), or from an unrelated family (column 4).

The values between parentheses present the identication rates obtained when p belong to one of the immunoglobulin typed-families (rst value) or to one of the enzyme typed-families (second value).

Method case I: p ∈ P -{p} case II: p ∈ P case III: p ∈ P c TEDt 41.6% (50.4%, 23.8%) 9.8% (6.5%, 13.4%) 48.6% (43.1%, 62.8%) TEDg: = 1Å 57.6% (69.5%, 33.6%) 3.3% (1.6%, 6.6%) 39.1% (28.9%, 59.8%) TEDg: = 2Å 58.1% (70.3%, 33.6%) 5.2% (2.8%, 9.8%) 36.7% (26.9%, 56.6%) Inria Figure 2 Morphology of the patches: number of atoms versus number of shells. The bold symbols identify the canonical morphologies presented on Fig. 3. For the sake of readability, the displayed range for the number of shells is [0, 20], so that the gure does not show the patches of the signaling complex 2oza (chain A: 26 shells and 395; chain B: 47 shells and 365 atoms).

RR n°7743

Figure 9 Rigid, topo-rigid, and exible patches: the three cases singled out on Fig. 8. On each row, the left and right columns respectively display the patch preimage and the patch. Top: for complex 2jel_HL-P, the chains HL dene a rigid patch, whose topology and geometry are preserved between the unbound and bound form (DIS t = 0.026, and DIS g = 0.058, associated with a RM SD d of 0.90Å). Middle: for complex 2i25_N-L, the chain L denes a topologically-rigid patch whose topology is preserved, but whose geometry is not (DIS t = 0.081, and DIS g = 0.505, associated with a RM SD d of 4.42Å). Bottom: for complex 1iqd_AB-C, the chain C denes a exible patch whose topology and geometry undergo signicant changes (DIS t = 0.464, and DIS g = 0.608, associated with a RM SD d of 2.46Å).

Inria

Figure 10 On morphologies and docking predictions. (A) Two complexes yield two comparisons between patches and their preimages, namely DIS x (R

(b) 1 , R (b) 2 ) versus DIS x (R (u) 1 , R (u)
2 ) and DIS x (L

(b) 1 , L (b) 2 ) versus DIS x (L (u) 1 , L (u)
2 ), where DIS x refers to the topological or geometric dissimilarity. (B) Topological dissimilarity (C) Geometric dissimilarity.

R (b) 1 R (b) 2 R (u) 1 R (u) 2 L (b) 1 L (b) 2 L (u) 1 L (u) 2

Receptors Ligands

Dening Patches

The identication of interface atoms is carried out by seeking edges present in the α-complex of the expanded balls, and whose endpoints belong to the two partners or involve an interfacial water molecule, see [START_REF] Cazals | Revisiting the voronoi description of protein-protein interfaces[END_REF][START_REF] Bouvier | Shelling the voronoi interface of proteinprotein complexes reveals patterns of residue conservation, dynamics and composition[END_REF][START_REF] Loriot | Modeling macro-molecular interfaces with intervor[END_REF]. Practically, the α-shape is computed using the Alpha_shape_3 package of the Computational Geometry Algorithms Library, see www.cgal.org. The underlying algorithm has randomized complexity O(n log n + k), with n the number of input balls and k the size of the outputthe number of simplices of the regular triangulation underlying the α-shape. This done, the complex is dissociated, and another 0-complex is computed for each partner, in the Solvent Accessible Surface (SAS) model. On a per-partner basis, this 0-complex is used to compute a combinatorial representation of the boundary of the expanded atomic balls [START_REF] Akkiraju | Triangulating the surface of a molecule[END_REF], stored in a half-edge data structure (HDS) [START_REF] Berg | Computational Geometry: Algorithms and Applications[END_REF]. The HDS consists of faces, half-edges, and vertices, and the connectivity information between these items allows in particular to nd the connected component of the patch, and to nd the cycles bounding a given connected component. For example, a patch consisting of a surface patch with a hole in the middle is bounded by two cycles, each consisting of a consecutive circle-arcs. In the sequel, such cycles are called Connected Components of the Boundary, or CCB.

From a geometric standpoint, a robust 3D embedding of the HDS is obtained using exact degree two algebraic numbers to represent the coordinates of the point lying at the intersection of three spheres [START_REF] De Castro | Design of the cgal spherical kernel and application to arrangements of circles on a sphere[END_REF].

Selected Properties of Patches

Before presenting the details of the topological encoding sketched in section 2.1, we discuss selected properties of patches.

Patches with multiple rims. In all generality, a connected component of the patch is not simply connected i.e. it may contain holes. This is illustrated on Fig. SI-1(c,d), where the packing defect in the middle of the interface is such that the patch is topologically equivalent to an annulus. To understand the implications of this fact on the shelling process, consider a connected patch with several CCB, and assume that the faces incident on the half-edges of these CCB have been initialized to one. Computing the shells as described in section 2.1 would result in a directed acyclic graph (DAG) rather than a tree, with a number of roots equal to the number of CCB.

To fudge around this diculty, the outer cycle only is used to initialize the SO calculation by tagging selected faces with a SO of zero. In doing so, the resulting DAG is a tree.

Patches with multiple connected components. Since a patch is dened as the SAS of the interface atoms of a sub-unit, it can be disconnected. This happens if the contact between partners has two distinct regions, which typically occurs for large protein interfaces [START_REF] Bahadur | A dissection of specic and nonspecic protein-protein interfaces[END_REF]. But as illustrated on Fig. SI-1(a,b), this also happens due to packing defects within one subunit. In this case, the shelling graph may contain several small connected components (cc).

Inria

Supporting Information Figure 1 However, one expects one component to contain signicantly more atoms than the remaining ones. Processing all patches and plotting the histogram of the size of the atom shelling trees yields the Fig. SI-2, which has an empty gap between the range [START_REF] Janin | Protein-protein interaction and quaternary structure[END_REF][START_REF] Halperin | Principles of docking: An overview of search algorithms and a guide to scoring functions[END_REF] and [START_REF] Kastritis | Are scoring functions in protein-protein docking ready to predict interactomes? clues from a novel binding anity benchmark[END_REF]438].

Practically, having computed the connected component of the patch, we remove all components containing less than 15 faces results in the removal of at most 15 atoms. In all cases processed, we are left with a unique cc, which is the largest one. This fact accounts for the name face shelling tree as opposed to face shelling forest.

RR n°7743

Step 5: Computing the Shelling Order. Using the connectivity of faces encoded in the HDS, a priority queue is used to assign the SO to all the faces. The queue is initialized with the boundary faces identied at step 4.

Step 6: Computing the shells. A shell being a connected component of faces having the same S0, a Union-Find algorithm is also called to create the shells.

Step 7: Computing the Face Shelling Graph. A parent-child relationship between two shells is witnessed by a half-edge incident on two faces having a SO which diers by one unit. Collecting all such pairs requires a linear pass over all half-edges. Constructing the Face Shelling Graph from the parent-child list is then straightforward.

Step 8: Selecting the Face Shelling Tree. So far, one tree has been computed for each connected component of the patch. We select the tree selected corresponds to the largest component in the Face Shelling Graph. As discussed above, this settles the case of patch with several connected components.

Step 9: Computing the Atom Shelling Tree from the Face Shelling Tree. Editing the atom shelling tree from the face shelling tree just requires handling atoms contributing several faces to the patch, as discussed in section 2.1.

Step 10: Ordering Atom Shelling Tree. This step requires sorting the sons of a node by increasing size.

6.2 Comparing Patches: a Generic Dynamic Programming Based Approach

The Tree Edit Distance

The generic TED. Given two ordered trees T 1 and T 2 , i.e. trees such that the children of each node are ordered, the Tree Edit Distance calculation aims at editing or morphing one tree into the other [START_REF] Bille | A survey on tree edit distance and related problems[END_REF]. The TED computation is actually based on three operations, namely deleting a node, inserting a node, and morphing a node of the rst tree into a node of the second tree. The output of the TED consists of an ordered edit distance mapping, namely a set

M ⊂ V ertices(T 1 ) × V ertices(T 2 ) such that for any pair (v 1 , v 2 ) ∈ M and (w 1 , w 2 ) ∈ M , one has: (i) v 1 = w 1 i v 2 = w 2 , or (ii) v 1 is an ancestor of w 1 i v 2 is an ancestor of w 2 , or (iii) or v 1 is to the left of w 1 i v 2 is to the left of w 2 .
(Recall that trees are ordered.) Call a node of a tree a paired node provided that it is involved in a morphing operation, and let N 1 (resp.

N 2 ) the nodes of T 1 (resp. T 2 ) which are not paired, and let λ be the empty node. If γ() refers to the cost of an insert/delete/morph operation, the cost of the edit distance mapping M is the following:

γ(M ) = (v,w)∈M γ(v → w) + v∈N1 γ(v → λ) + w∈N2 γ(λ → w) (8)
From which one denes the TED as:

T ED = min M :Edit Distance Mapping γ(M ). (9)
It can be shown that the TED calculation is amenable to a dynamic programming approach, which we sketch the sake of completeness.

The recursive structure of the TED. Computing the TED actually requires handling ordered forests rather than trees. Indeed, removing the root of an ordered tree leaves a forest of ordered treesone tree for each son of the root. From now on, we consider two forest F 1 and F 2 . Denoting v and w the rightmost (if any) roots of F 1 and F 2 , and F 1 (v) the sub-tree rooted at vand likewise for F 2 . It can be shown that the TED calculation has the following recursive structure:

T ED =                    δ(∅, ∅) = 0 δ(F 1 , ∅) = δ(F 1 -v, ∅) + γ(v → λ) δ(∅, F 2 ) = δ(∅, F 2 -w) + γ(λ → w) δ(F 1 , F 2 ) = min      δ(F 1 -v, F 2 ) + γ(v → λ) δ(F 1 , F 2 -w) + γ(λ → w) δ(F 1 (v), F 2 (w)) + δ(F 1 -F 1 (v), F 2 -F 2 (w)) + γ(v → w) (10
) These equations show that: the value of δ(F 1 , F 2 ) depends on a constant number of problems of smaller size; each sub-problem can be computed in constant time.

The optimal algorithm developed in [START_REF] Demaine | An optimal decomposition algorithm for tree edit distance[END_REF] has cubic time complexity, and quadratic memory requirements. Note that the TED problem for unordered trees is in general NP-hard [START_REF] Bille | A survey on tree edit distance and related problems[END_REF].

Instantiation in the context of Atom Shelling Trees. The three operations insert/delete/morph are generic in the sense that they depend on the semantics associated to the nodes. In our case, a node of the shelling tree corresponds to a set of atoms, and dierent interpretations can be used, as we have seen in section 2.2: while focusing solely on the number of common atoms yields a topological comparison, namely algorithm TED t , focusing on quasi-isometric subsets of atoms yields a geometric comparison, namely algorithm TED g . In the next section, we explain how the geometric comparison of two shells is carried out, which is underlying the morph operation of TED g .

Comparing Shells

Our strategy to compare shells is a modication of the one proposed in [START_REF] Malod-Dognin | Maximum clique in protein structure comparison[END_REF] in the context of protein's alpha-carbon backbone comparison, and reduces to a maximum clique calculation. To present it, we shall need the following denitions and notations.

Denition. 1 A m × n 2D graph G = (V, E) is a graph in which the vertex set V is depicted by a (m-rows) × (n-columns) array T , where each cell T [i][k] contains at most one vertex i.k from V (note that for both arrays and vertices, the rst index stands for the row number, and the second for the column number). Two vertices i.k and j.l can be connected by an edge (i.k, j.l) ∈ E only if i = j and k = l.

Denition. 2 A clique of a graph G = (V, E) is a subset of its vertex set V such that any two vertices in it are adjacent (i.e. connected by an edge in E).

Inria

Denition. 3 The maximum clique problem (also called maximum cardinality clique problem) is to nd a largest, in terms of vertices, clique of an arbitrary undirected graph G.

In matching two shells s 1 ⊂ BP 1 and s 2 ⊂ BP 2 , the goal is to nd a one-to-one correspondance between two sets of atoms m 1 ⊆ s 1 and m 2 ⊆ m 2 . Since we aim, following Eq. ( 1) at nding quasi-isometric subsets, we shall use constraints. Assume that we wish to match atom i ∈ s 1 with atom k ∈ s 2 , and similarly atom j ∈ s 1 with atom l ∈ s 2 , and let d 1 i.j (d 2 k,l ) be the distance between atoms i and j (resp. k and l). The compatibility constraints between the two pairs go as follows:

1. i = j and k = l; this constraints ensures that an atom of s 1 can be matched with at most one atom of s 2 , and vice versa;

2. for a distance threshold , we impose |d 1 i.j -d 2 k.l | ≤ ; this constraints ensure that the RM SD d of internal distances is upper-bounded by .

A feasible matching is thus a sequence of matching pairs i 1 ↔ k 1 ,i 2 ↔ k 2 ,. . . , i n ↔ k n such that any two pairs are compatible. Searching the largest feasible matching can be rephrased in a |s 1 | × |s 2 | 2D graph G = (V, E) in the following way. Each row i of V represents an atom i ∈ s 1 , and each column k represents an atom k ∈ s 2 . For all possible matching pairs i ↔ k, we create a vertex i.k ∈ V , on row i, column k. For all compatible couples of matching pairs i ↔ k and j ↔ l, we create an edge (i.k, j.l) ∈ E. A feasible matching corresponds to a clique in G, and the longest feasible matching to a maximum clique in G (Fig. SI-3).

Comparing two shells is modeled as nding a maximum clique in a graph. The maximum clique problem is one of the rst problem shown to be NP-Complete [START_REF] Karp | Reducibility among combinatorial problems[END_REF], and it has been studied extensively in literature. Interested readers can refer to [START_REF] Bomze | The maximum clique problem[END_REF] for a detailed state of the art about the maximum clique problem. In our current implementation, the maximum clique in the 2D graph is computed using the Cliquer library [START_REF] Östergård | A fast algorithm for the maximum clique problem[END_REF]. As specied by Eq. ( 1), the algorithm TED g involves a distance threshold . More precisely, recall that matching atom i ∈ s 1 with atom k ∈ s 2 , and atom k ∈ s 1 with atom l ∈ s 2 requires |d i.jd k.l | ≤ . The parameter aects both the quality of the comparisons and their computation times. Quality-wise, the larger , the larger the internal distance discrepancies allowed, whence the larger and the less similar the common subsets of atoms returned. Computationally, the larger , the more dicult the identication of quasi-isometric subsets. On computers with Intel Xeon processors at 2.66Ghz, computing the 124251 pairwise comparisons of our database was done in about 622 seconds by TED t , in about 41901 seconds by TED g ( = 1Å), and 1076390 seconds by TED g ( = 2Å). For the later, a time limit of 2 hours per instances was used, and 13 instances remained unsolved. In this study, we compared TED t against TED g with = 2Å. As shown while discussing the identication ratesTable 1, = 2 strikes a balance between the quality of the comparison and the hardness of the computations.

Supporting Information

Dataset of Biological Complexes

Our approach is validated on a dataset of 498 patches generated from 249 complexes. The set of complexes was assembled from two sources. The rst one is the IMGT 3D structure database [START_REF] Ehrenmann | Imgt/3dstructure-db and imgt/domaingapalign: a database and a tool for immunoglobulins or antibodies, t cell receptors, mhc, igsf and mhcsf[END_REF], from which we extracted 116 immunoglobulin -ligand structure, with resolution ≤ 2.0.

The second one is the recently assembled anity benchmark [START_REF] Kastritis | A structure-based benchmark for protein-protein binding anity[END_REF], a manually curated dataset involving 133 complexes with experimentally measured binding anity, and resolutions in the range 1.1 -3.3, the median being 2.4. (In fact, the anity benchmark contains 144 complexes, but 11 were redundant with the ones that we extracted from the IMGT 3D database.)

Note that while most of the patches possess between 14 to 294 atoms and 1 to 20 shells, the two patches from the signaling complex 2oza are larger: 2oza chain B possesses 395 atoms and 28 shell (Fig. SI-4), and 2oza chain A possesses 365 atoms and 47 shells (Fig. SI-5).

By distinguishing the type (receptor, ligand) of each partner, and using biological information on each complex, these patches are classied into 19 so-called typed families. Supporting Information Figure 8 The quality of the identication is related to the consistency of the families. For each typed family P , the correct identication rate is plotted as a function of -log(P value ), P value being the p-value of the Wilcoxon-Mann-Whitney ranksum test between s(P, P ) and s(P, P c ). Very low P values (consistent families) relate to high identication rates.

RR n°7743

Supporting Information Figure 9 Scatter plots of binding anity against various structural parameters (A) Anity versus Internal Path Length (C Pea = 0.31 and C Spe = 0.43). (B) Anity versus number of atoms (C Pea = 0.27 and C Spe = 0.37) (C) Anity versus depth (C Pea = 0.29 and C Spe = 0.35) (D) Anity versus ∆ASA (C Pea = 0.22 and C Spe = 0.33).

Note that A, B and C use the average value of the considered parameter across the two partner patches of each complex. 6.9 Software: application and le formats

In this section, we describe the tools that we designed for generating and comparing atom shelling trees. The two softwares, VORPATCH and COMPATCH, are available from http://cgal.inria.fr/abs/vorpatchcompatch/.

RR n°7743

Figure 1

 1 Figure 1 Shelling a patch: illustration. (a) Side view of the protein complex 1vfb, an immunoglobulin -antigen complex, with interface atoms displayed in red for partner A (chain A and B) and in blue for partner B (chain C). Gray atoms correspond to the interfacial water molecules. (b,c) Rotated view of the patch of partner A, shelled into concentric shells of atoms, and corresponding atom shelling tree. The colors of the atoms match those of the nodes of the shelling tree, the non interface atoms being represented in gray-blue, and the atoms of a shell are represented with a given color.

Figure 3 Figure 4 Figure 5

 345 Figure 3 Illustration of the ve canonical morphologies. The ligands are represented as cartoons when they do not clutter the picture. (a) tubular (pdbid 3eys, chain Q) (b) pear-like (pdbid 2dqu, chain 1) (c) isotropic-at (pdbid 2ih3, chains AB) (d) isotropic-pyramidal (pdbid 3a6c, chain Y) (e) anisotropic-at (pdbid 3h0t, chains AB) 1 5 2 8 3 12

Figure 6

 6 Figure 6 An example complex whose asymmetry is revealed by the average Shelling Order. Complex 3i, from A_Pept. Left: the receptor (chains HL). Non interface atoms are displayed in grey, and ligand is displayed with a cartoon representation. Right: the ligand (chain P) alone. The SO of the receptor vary in the range [0, 4], with an SO of 1.17, while for the ligand the SO vary in range [0, 7] with an SO of 3.47.

Figure 7

 7 Figure 7 Topological (DIS t ) versus geometrical (DIS g ) dissimilarities. With one point per pairwise comparison, the plot illustrates the fact that the topological dissimilarity is a lowerbound of the geometrical dissimilarity. The black rectangle singles out instances having similar topologies (low DIS t ) but dierent geometries (high DIS g ).

Figure 8

 8 Figure 8 Topological versus geometric dissimilarities for pairs (patch, patch preimage on the unbound partner) of the binding anity benchmark. The square and circle correspond to rigid and fully exible docking, respectively. The triangle corresponds to a patch whose geometry but not topology changes upon binding: such a patch, called topo-rigid, is preformed on the unbound partner. The three highlighted examples are presented in gure 9

1 A

 1 Hierarchical Encoding of PatchesIn this section, we present the details of the atom shelling tree construction, sketched in section 2.1, and illustrated on Fig.SI-1.

  Dening and shelling a patch: 2D illustration. (a) The endpoint of the dashed purple edges, which are dual of the solid purple edges of the Voronoi diagram, identify the interface atoms. (b) The patch of the blue subunit is the Solvent Accessible Surface of its interface atoms, represented as solid circle arcs. Note that the packing defect in-between the atoms centered at b 1 , b 2 , b 4 , b 5 is such that the patch has two connected components cc 1 and cc 2 . (c) A packing defect in-between the partners dismisses atom a 3 as interface atom. (d) On this 2D example, the patch has two c.c.; when the same occurs in 3D, the patch is connected but is topologically equivalent to an annulus i.e. has two rims called connected component of the boundary or CCB. The largest one, namely the outer one, only is used to initiate the shelling order calculation.

Figure 3 4 .

 34 Classical versus 2D graph representation of feasible matching computation. Left: The red arrows correspond to the feasible matching 1 ↔ 1, 2 ↔ 2, 3 ↔ 4, which implies that both d 1 Right: the same matching is represented in a 2D-graph.

Table 1

 1 Table SI-1 displays the PDB ids and the partner specications of our dataset, while Table SI-2 presents the typed families. The 249 protein complexes used in this study. Geometric dissimilarity (DIS g ) versus RMSD of internal distances (RMSD d ) computed at = 2Å. The RM SD d globally increases with the geometric dissimilarity, and so does the variance of the RM SD d for a xed DIS g . The black box singles out instances with low geometric dissimilarity but high RM SD d values.

	Inria

Table 3

 3 Correct identication rates per typed family. The typed family of a patch is identied by the typed family of its nearest neighbors p, according to the geometric dissimilarity (with a distance threshold of 2Å). Columns 2 (resp. 3) presents for each typed family the number (resp percentage) of correctly identied binding patches. The low identication rates can be explained by inconsistencies in the original classication. Such an inconsistency is presented on Fig.SI-6, where both 1jps chains HL and 1dqj chains BA RR n°7743 come from the AA_Prot_R typed family, but only have about 32% of geometrically common atoms.Supporting Information Table4p-values for the Wilcoxon-Mann-Whitney rank-sum test for similarity scores obtained by TED g with =2Å. The null hypothesis (identical distributions) is rejected for all the immunoglobulin families, while in the case of the enzymes, it is only rejected for the E_Inhi family. Note that these probabilities closely follow the identication rates observed in table SI-3.Family (=P) s(P, P ) vs s(P, P ) s(P, P ) vs s(P, P C )

	Typed family Correctly identied patches percentage
	A_Carb_R	4/9	44.4%
	A_Carb_L	3/9	33.3%
	A_Chem_R	37/40	92.5%
	A_Chem_L	34/40	85.0%
	A_Pept_R	7/21	33.3%
	A_Pept_L	8/21	38.1%
	A_Prot_R	39/53	73.6%
	A_Prot_L	41/53	77.4%
	E_Inhi_R	20/40	50%
	E_Inhi_L	21/40	52.5%
	E_Regu_R	0/11	0%
	E_Regu_L	0/11	0%
	E_Subs_R	2/10	20%
	E_Subs_L	0/10	0%

Table 5

 5 The Internal Path Length yields the best correlation against the binding anity (-lnK d ). The values were computed for the 144 complexes of the anity benchmark. Note that MIC p-values are available only for coecient larger than 0.275, see http://www.exploredata.net/Downloads/P-Value-Tables.

		Pearson	Spearman	Maximal Information
	Parameter IPL	C Pea p-value 0.31 1.3e-4	C Spe p-value C MIC 0.43 7.6e -8 0.35	p-value 7.6e -4
	#Atoms	0.27	1.2e-3	0.37	4.7e -6	0.24
	Depth	0.29	4.8e-4	0.35	1.5e -5	0.26
	∆ASA	0.22	8.9e-3	0.33	6.6e -5	0.25
	Firedock score	-0.17 4.2e-2	0.20	1.8e -2	0.23
	I_RMSD	-0.11 2.0e-1	0.17	4.3e -2	0.24
	#Shells	0.092 2.7e-1	-0.16	5.4e -2	0.16
	DIS g	0.16	5.8e-2	-0.14	8.5e -2	0.24
	Assymetry	0.045 5.9e-1	-0.094 2.6e -1	0.19
	DIS t	0.029 7.2e-1	-0.089 2.9e -1	0.20
	Supporting Information

Table 6

 6 Spearman's correlation coecient as a function of the exibility observed upon binding. Spearman coecient between the binding anity ln K d and ∆ASA, #Atoms, Depth and IPL, for the three exibility classes introduced in [17]: I-RMSD < 1 (rigid interfaces, 71 cases), I-RMSD ∈ [1 , 1.5 [ (semi-rigid interfaces, 38 cases), and I-RMSD ≥ 1.5 (exible interfaces, 35 cases). Spe p-value C Spe p-value C Spe p-value C Spe p-value Supporting Information Table 7 Maximal Information Coecient as a function of the exibility observed upon binding. MIC between the binding anityln K d and ∆ASA, #Atoms, Depth and IPL, for the three exibility classes introduced in [17]: I-RMSD < 1 (rigid interfaces, 71 cases), I-RMSD ∈ [1 , 1.5 [ (semi-rigid interfaces, 38 cases), and I-RMSD ≥ 1.5 (exible interfaces, 35 cases).

		∆ASA	#Atoms	Depth		IPL
	I-RMSD (Å) C < 1 Å 0.52	3.5e -6	0.58	1.4e -7	0.54	9.0e -7	0.59	5.9e -8
	in [1Å,1.5Å[	0.18	2.7e -1	0.11	5.0e -1	0.054 7.5e -1	0.23	1.7e -1
	≥ 1.5Å	0.26	1.2e -1	0.34	4.7e -2	0.34	4.2e -2	0.41	1.5e -2
		∆ASA	#Atoms		Depth		IPL
	I-RMSD (Å) C MIC p-value C MIC p-value C MIC p-value C MIC p-value < 1 Å 0.39 1.5e -2 0.43 3.6e -3 0.43 3.2e -3 0.48 4.4e -4
	in [1Å,1.5Å[	0.60	3.6e -4	0.26	_	0.28	_	0.34	_
	≥ 1.5Å	0.22	_	0.36	_	0.33	_	0.31	_

RR n°7743

route des Lucioles -BP 93 06902 Sophia Antipolis Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Supporting Information Figure 2 Frequencies of connected component of given size.

Computed over the 498 shelling graphs generated from our dataset, the histogram presents an empty gap between the range [START_REF] Janin | Protein-protein interaction and quaternary structure[END_REF][START_REF] Halperin | Principles of docking: An overview of search algorithms and a guide to scoring functions[END_REF] and [START_REF] Kastritis | Are scoring functions in protein-protein docking ready to predict interactomes? clues from a novel binding anity benchmark[END_REF]438]. On patches and atom selected. A nal comment is in order to qualify the atoms which are involved in our patch construction, in particular with respect to interface model based on the loss of solvent accessibility ∆ASA > 0. As observed in [START_REF] Cazals | Revisiting the voronoi description of protein-protein interfaces[END_REF] and explained in [START_REF] Cazals | Revisiting the Voronoi description of protein-protein interfaces: Algorithms[END_REF], some interface atoms selected by the Voronoi model may no lose solvent accessibility. This happens in particular for atoms which are buried in their sub-unit: such atoms are interface atoms in the Voronoi model, but are excluded from the patch model since they do not contribute to the SAS of the sub-unit. Such cases represent less than 10% of all interface atoms [START_REF] Cazals | Revisiting the voronoi description of protein-protein interfaces[END_REF].

Algorithm

We are now ready to detail the algorithm sketched in section 2.1.

Step 1: Computing the HDS. The half-edge data structure encodes the boundary of the union of balls, as computed in [START_REF] Akkiraju | Triangulating the surface of a molecule[END_REF]. A certied embedding in 3D is obtained thanks to the robust geometric operations described in [START_REF] De Castro | Design of the cgal spherical kernel and application to arrangements of circles on a sphere[END_REF].

Step 2: Computing the Connected Component of the Boundary (CCB). The CCB are the cycles bounding the patches. Given the HDS, nding all CCB of a patch requires running a Union-Find algorithm [START_REF] Tarjan | Data Structures and Network Algorithms[END_REF], which has (almost) linear complexity.

Step 3: Computing the Connected Component of Half-edges (CC). To identify the connected components of a patch, we run a Union-Find algorithm on all the half-edges of the patch. Note that each c.c. will yield a shelling graph/tree.

Step 4: Initializing the Shelling Order. From steps 2 and 3, the largest CCB of each connected component is selected, and the corresponding faces are assigned a SO of zero. This step settles the case of connected component with several rims.

Inria

Supporting Information Figure 4 The patch having the largest number of atoms.

2oza chain B possesses 395 atoms and 28 shells.

Supporting Information Figure 5 The patch having largest number of shells. 2oza chain A possesses 365 atoms and 47 shells.

Inria

Supporting Information Table 2 The typed family classication of the 498 patches from our database. Stars indicate the typed families used for the identication and clustering experiments. The OG, OR and OX classes are from [START_REF] Kastritis | A structure-based benchmark for protein-protein binding anity[END_REF]. Supplemental Figure 1 Overview of the patch generation with COMPATCH.

Family of complex

Supplemental Figure 2 The atom shelling tree (.ast) le format. An atom shelling tree is described by the list of its nodes (or shells). Each node is rst described by a header line containing three integers: the node ID (starting from 1), the number of atom of the corresponding shell, and the node's father ID (0 if a node is a root one). A header line is then followed by lines describing the shell's atoms: one line per atoms, each containing the pdb ID of the atom, its x, y and z coordinates and its expanded radius. # header of the first node (node ID, #atoms, father ID) Given two patches (.ast le format), COMPATCH use the tree-edit-distance based methods presented in section 6.2.1 to measure their dissimilarity, and also record the optimum tree-edit-script (the sequence of tree-edit operations) in the custom le format described in the supplemental Fig. 4. The numerical values (dissimilarity scores, size of the two trees, tree-edit-distance values and running times) are printed into the console or the log le.

Supplemental Figure 4 The tree-edit-script (.tes) le format. The rst line recall the lename of the two input atom shelling trees. The consecutive lines present the optimum treeedit script (sequence of tree-edit operations) for transforming the rst tree into the second one, and for each operation the associated cost is given. If the comparison was done with TED g , then the mapping operations are followed by the corresponding lists of atom matchings.

.