Frédéric Cazals

Sébastien Loriot

Computing the Arrangement of Circles on a Sphere, with Applications in Structural Biology *

Balls and spheres are the simplest modeling primitives after affine ones, which accounts for their ubiquitousness in Computer Science and Applied Mathematics. Amongst the many applications, we may cite their prevalence when it comes to modeling our ambient 3D space, or to handle molecular shapes using Van der Waals models. If most of the applications developed so far are based upon simple geometric tests between balls, in particular the intersection test, a number of applications would obviously benefit from finer pieces of information.

Consider a sphere S0 and a list of circles on it, each such circle stemming from the intersection between S0 and another sphere, say Si. Also assume that Si has an accompanying ball Bi. This paper develops an integrated framework, based on the generalization of the Bentley-Ottmann algorithm to the spherical setting, to (i)compute the exact arrangement of circles on S0 (ii)construct in a single pass the half-edge data structure encoding the arrangement induced by the circles (iii)report the covering list of each face of this arrangement, i.e. the list of balls containing it. As an illustration, the covering lists are used as the building block of a geometric optimization algorithm aiming at selecting diverse conformational ensembles for flexible protein-protein docking.

Introduction 1.Modeling with Balls and Spheres

Balls in Geometric Modeling. Balls and spheres are the simplest modeling primitives after affine ones. This accounts for their ubiquitousness in computer science and applied mathematics, especially when it comes to modeling our ambient 3D space. Given a known object, balls can be assembled to provide hierarchical approximations. In computer graphics, such approximations have been used to perform multi-scale visualization [START_REF] Rusinkiewicz | QSplat: a multiresolution point rendering system for large meshes[END_REF], while in robotics, they have been instrumental to perform efficient collision detection [START_REF] Agarwal | Collision detection for deforming necklaces[END_REF]. Balls are also key in inferring geometric and topological informations of an object known from sample points. For example, under mild sampling assumptions, mimicking the medial axis transform of a sampled surface allows one to reconstruct the domain enclosed within this surface [START_REF] Amenta | The power crust, unions of balls, and the medial axis transform[END_REF]. From a more geometric perspective, a number of quantities can be inferred from balls centered at the sample points.

As an illustration, one may cite the so-called boundary measure of a point cloud, from which singular points and sharp features of the sampled model can be estimated [START_REF] Chazal | Stability of boundary measure[END_REF]. Additional scenarios involving balls can be found across the computer science literature in general, and we refer for example the reader to chapter 11.12 of the Visionbib bibliography project at http://www.visionbib.com/bibliography/ describe482.html for numerous pointers. Balls in Structural Biology. Balls are also central in structural biology, especially when it comes to manipulating Van der Waals (VdW) models. In a VdW model, each atom is represented by a ball whose radius depends on the atom type and its environment. Such models are instrumental to investigate implicit solvent models for electrostatics [START_REF] Koehl | Electrostatics calculations: latest methodological advances[END_REF], to compute energies and statistical potentials [START_REF] Summa | An Atomic Environment Potential for use in Protein Structure Prediction[END_REF], to define and model geometric shapes for docking [START_REF] Shoichet | Molecular docking using shape descriptors[END_REF], and to define molecular surfaces [START_REF] Connolly | Molecular surfaces: a review[END_REF]. This latter aspect has been of special interest in computational geometry, since two of the most widely used molecular surfaces read directly from union of balls [START_REF] Akkiraju | Triangulating the surface of a molecule[END_REF]. Given a collection of atomic balls, its VdW surface is the boundary of its VdW balls, while its Solvent Accessible Surface (SAS) is the boundary of the union of balls after expanding their radii by r w = 1.4 Å to account for a solvation layer. Arrangements of Circles on a Sphere. Most of the applications aforementioned are based upon simple geometric tests between balls, in particular the intersection test. A number of them, though, benefit from finer pieces of information. To see which, consider a sphere S 0 and a list of circles on it, each such circle stemming from the intersection between S 0 and a sphere S i . The arrangement of circles on S 0 is the partition of S 0 into regions whose interior is connected. Assuming that each sphere S i is associated a ball B i , the covering list of a face of the arrangement is the list of balls that contain it, and the multiplicity of the face is the size of the covering list. Since our interest comes from structural biology, let us mention briefly two problems directly concerned with arrangements and covering lists.

The first one is the problem of encoding multi-body contacts between an atom and its neighbors. In representing a molecule as a collection of balls, we recalled above the importance of molecular surfaces. For a given atomic sphere within a molecule, consider the arrangement induced by the intersection circles with neighboring atoms. The contribution of this atom to the molecular surface, say the VdW or the SAS surface, consists of the cells of null multiplicity-the surface features the boundary of the union of balls i.e. the exposed spherical caps. Fig. 1 features the cumulative area (over all atoms of a complex) as a function of the multiplicity. Notice in particular that the surface area of the SAS corresponds to a mere 4.5% of the total computed surface area. Given that all previous studies overlooked faces of multiplicity ≥ 1, arrangements of circles hold great promises to refine our understanding of inter-atomic multi-body contacts.

The second problem is related to flexible docking. Recall that docking is concerned with the prediction of the geometry of a complex, say a protein-protein complex, from those of its constitutive partners. For partners undergoing significant conformational changes upon binding, a number of docking algorithms resort to pre-generated conformations of the partners called conformers. To maximize the chances of a docking experiment, one naturally wishes to deal with conformers providing a sampling of the conformational space as good as possible. Given a large number of conformers, we shall present in section 5 an algorithm based on arrangements of circles on a sphere to select diverse conformational ensembles. One such ensemble is presented on Fig.

Contributions and Paper Overview

Arrangement of Circles on a Sphere and Related Problems. We consider a collection of n balls B i,i=1,...,n intersecting a given ball B 0 . The sphere associated to ball B i is denoted S i , and the intersection circle, if any, between S 0 and S i is denoted C i . For a collection of n circles on sphere S 0 , recall that the arrangement induced by the circles is the decomposition of S 0 into faces, circular arcs and vertices. In the following, we shall be concerned by the following two problems: reporting the arrangement on S 0 , and reporting for each cell of this arrangement its covering list-the list of balls that contain it. Contributions and Paper Overview. The punchline of the paper is to present an integrated framework allowing one to construct the arrangement of circles on a sphere and to compute a compact representation of the covering lists called the covering tree. Thus, we make three contributions.

First, we present the first effective algorithm to compute the exact arrangement of circles on a spheresection 2. The algorithm, which generalizes the Bentley-Ottmann [START_REF] Bentley | Algorithms for reporting and counting geometric intersections[END_REF] sweep to the spherical setting, uses a segregation of events specific to circles, and gathers events into a data structure called the event site, so as to handle degeneracies. In passing, we mention the fact the complexity of the algorithm involves faces bounded by two arcs (the so-called lenses and lunes). Second, we present the construction of the half-edge data structure (HDS) storing the arrangement, a construction performed on the fly during the sweep process, and using Union-Find-section 3. As an application of this construction, we present a calculation of the covering tree of the arrangement-section 4. Finally, section 5 presents an application of covering lists to flexible docking. Position with Respect to Previous Work. For the calculation of the arrangement, the only effective alternative to our Bentley-Ottmann variation consists, to the best of our knowledge, of computing a trapezoidal map, based upon an explicit perturbation of the spheres so as to get rid of degeneracies [START_REF] Halperin | A perturbation scheme for spherical arrangements with application to molecular modeling[END_REF][START_REF] Eyal | Dynamic maintenance of molecular surfaces under conformational changes[END_REF]. Our strategy cumulates several advantages. First, the computation of the arrangement is exact, exactness being an obvious route to warrant robustness. (We note that while exactness is not a prerequisite in computational structural biology, it might be a must in other applications.) Second and most importantly, running Bentley-Ottmann on the sphere provides a direct route for the construction of the arrangement and of the covering lists. While the classical way to obtain the arrangement consists of merging the cells of a trapezoidal decomposition [START_REF] Boissonnat | Algorithmic Geometry[END_REF], we perform it on the fly using union-find-an original alternative to the best of our knowledge.

As a general comment, notice that using the stereographic projection, the arrangement of circles on a sphere could be obtained from an arrangement of circles and lines in the plane [START_REF] Wein | Advanced programming techniques applied to CGAL's arrangement package[END_REF]. However, the spherical setting alleviates the calculations within predicates. For example, the z coordinates of the socalled critical points have rational coordinates on the sphere [START_REF] De Castro | Design of the CGAL spherical kernel and application to arrangements of circles on a sphere[END_REF], while the x and y coordinates of their counterparts in the plane are algebraic numbers. In passing, we refer the reader to our companion paper [START_REF] De Castro | Design of the CGAL spherical kernel and application to arrangements of circles on a sphere[END_REF] for the predicates and constructions required by the algorithms presented herein.

2 Bentley-Ottmann on a Sphere

Algorithm

Sweeping θ-monotone Circular Arcs on a Sphere. To report line-segments intersections in the plane [START_REF] Berg | Computational Geometry: Algorithms and Applications[END_REF], the Bentley-Ottmann (BO) algorithm [START_REF] Bentley | Algorithms for reporting and counting geometric intersections[END_REF] requires two data structures which are an event queue E featuring the line segment endpoints and intersection points, and another sorted data structure V providing the ordering along the vertical sweep line. The data structure V features segments intersected by the sweep-line, which are pairwise checked for intersection when they become adjacent in V.

Consider a sphere S 0 . For circles on that sphere, the extension of the BO algorithm consists of sweeping the sphere with a meridian M θ , using cylindrical coordinates. This meridian is anchored at the poles and revolves around the sphere while θ spans the interval (0, 2π]. The sweep process handles θ-monotone circular arcs, which are stored in the vertical ordering V. To specify these arcs from the circles, we first introduce the following classification of circles: Definition 1 A circle C i is called polar if it goes through a single pole of S 0 ; bipolar if it goes through the two poles of S 0 ; threaded if the intersection point between the plane of C i and the z-axis belongs to the open disk bounded by C i ; normal otherwise.

The tangency point(s) between M θ and a circle generates one or two θ-monotone circular arcs. More precisely, a normal circle yields two θ-monotone circular arcs; a polar circle defines one such arc, the second one reducing to a point (its pole). From now on, these θ-monotone circular arcs are just called arcs. During the sweep process, an arc is active while it is in V. For homogeneity, we consider that a threaded circle defines an arc which is always active. Event Sites. To accommodate degeneracies, the event queue E stores event sites, which are defined from two types of event points (critical points and intersection points) as follows.

The normal critical points (θ c , z c ) associated to a normal circle are the two points where the meridian intersects it in a single point. Amongst these points, the start point (respectively end point) is the one where the intersection between the circle and the meridian starts (respectively ends). Consider a polar circle. Denoting z p the z coordinate of the corresponding pole, the polar critical points are defined as the pairs (θ S , z p ) and (θ E , z p ), with θ S and θ E the values such that the meridian is tangent to the circle. The start point is distinguished from the end point as for normal circles. This extension carries to bipolar circles, yielding bipolar critical points, the θ values being those where the circle is included in the plane containing the meridian-the z coordinate is irrelevant. For such a circle, the start (respectively end) critical point is the one corresponding to the smallest (respectively largest) value of θ. An intersection point between two arcs is either a tangency point, if the arcs are tangent in their interior, or a crossing point, if the arcs intersect transversely in their interior.

The event points just described are encapsulated within events, each featuring either a circle and a start/end tag, or a pair of arcs and a first/second tag to define the corresponding event point. Notice that the qualifiers introduced above for event points are inherited by events. Equipped with these notions, events are gathered into an event site as follows: Definition 2 A normal event site is a collection of events with the same event point, partitioned into the following three data structures: -One list F for normal end (finish!) events, sorted by increasing circle radius.

-One list S for normal start events, sorted by decreasing circle radius.

-One list CT which contains the crossing events and tangency events. The restriction of CT to crossing (tangency) events is denoted C(T ).

Note that a (bi)polar event site associated to one (bi)-polar event is the event itself. A collection of event sites corresponding to distinct event points are sorted using the lexicographic order on cylindrical coordinates. Difficulties arise due to collisions between events at a given θ value and/or at the poles: The fact that this total order is compatible with the construction of the arrangement and the calculation of the covering lists will become apparent later.

Complexity Analysis

To state the algorithm complexity, one needs considerations on the structure of the arrangement induced on S 0 by the intersection circles. A vertex of this arrangement is either an intersection point, or a critical point. Denoting n the number of circles, and k (respectively v) the number of intersection points (respectively vertices), we have v = O(n + k). An edge is an arc in-between two vertices. A face is a region of S 0 whose interior is connected. If the interior of a face is not simply connected, then the face contains hole(s). Using Euler's relationship on simply connected regions and holes of the decomposition of S 0 induced by the intersection circles, one can prove the following lemma [START_REF] Cazals | Computing the exact arrangement of circles on a sphere, with applications in structural biology[END_REF]:

Lemma 1 Denote v and e the number of vertices and edges of an arrangement of n circles on the sphere, and let l stand for the number of faces bounded by exactly two edges. One has e ≤ 3(v -1) + l.

Notice this lemma involves the number of faces bounded by exactly two edges, also known as lenses (convex faces) and lunes (concave faces). For n circles of arbitrary radii in the plane, this number is known to be O(n 3/2+t ), for any t > 0, where the constant of proportionality depends on t [START_REF] Alon | On the Complexity of Arrangements of Circles in the Plane[END_REF]. Using this lemma, the following theorem is proved in [START_REF] Cazals | Computing the exact arrangement of circles on a sphere, with applications in structural biology[END_REF]:

Theorem 1 Reporting the k intersection points of a family of n circles on a sphere requires O(n) storage and has O((n + k + l) log n) complexity.

Qualifying Arcs and Half-edges

Having sketched the BO algorithm, we present a classification of arcs and half-edges, to be used in sections 3 and 4.

Qualifying an Arc w.r.t. Circles: Upper or Lower. A threaded circle is called north threaded if its center has a z coordinate larger than or equal to that of the center of S 0 , and south threaded otherwise. In the same way, a polar circle containing the north pole is called a north polar circle, and south polar circle otherwise. Considering the two arcs of a normal circle, we call them the upper and the lower arcs w.r.t. the z axis. A north (respectively south) polar circle has a single non trivial arc, which is lower (respectively upper ). Similarly, we consider that a north (respectively south) threaded circle has a single arc which is lower (respectively upper ). Qualifying a Half-edge w.r.t. an Arc: Upper or Lower. For all but bipolar circles, to each arc A i , we associate two active half-edges qualified w.r.t. the unique intersection between the meridian and the arc: The upper (respectively lower ) half-edge associated to A i is the half-edge leaving to its left the portion of S 0 lying above (respectively below) A i . We now address particular cases: (i)for a threaded circle which is not intersected by any other circle, at the end of the sweep process, the source and the target of the half-edges associated to its arcs are fixed to a common null vertex (ii)for a bipolar circle, no arc is defined. Half-edges are created on the fly at a bipolar event, so as to handle, if any, intersections with active arcs.

Qualifying a Half-edge w.r.t. a Circle: Inner or Outer. For a circle which is not a great circle, consider the spherical cap of S 0 of smallest area induced by this circle. If a half-edge on this circle induces this cap, it is called inner, and outer otherwise. For a great circle, we break the tie and call the cap of smallest area that induced by inner half-edges, that is: (i) if the great circle is bipolar, an half-edge is inner if it induces the spherical cap swept by M θ between its θ S and θ E associated values, and outer otherwise; (ii) if the great circle is threaded, an half-edge is inner if it induces the spherical cap containing the north pole, and outer otherwise.

The relationship between the qualifiers upper/lower and inner/outer of half-edges is illustrated on Fig. 3, and we have:

Observation 1
The spherical cap of smallest (respectively largest) area bounded by a circle is described by inner (respectively outer) half-edges of the circle. The lower/upper half-edge of an upper (respectively lower) arc of a all but bipolar circle is always inner/outer (respectively outer/inner).

Cn An An Ct 1 At 2 At 1 Ct 2 (a) (b) (c) (d) Cbp Cp 1 Ap 2 Ap 1 Cp 2 x y z

Constructing the Arrangement During the Sweep Process

In this section, we develop the topological operations required to construct the arrangement induced by the circles and stored it into an extended half-edge data structure [START_REF] Kettner | Halfedge data structures[END_REF] (HDS) possibly featuring holes in faces.

Describing the Arrangement

Arcs and Half-edges. The vertices of the HDS correspond to event points, its edges are arcs delimited by such vertices, while the faces are the regions of S 0 bounded by vertices and edges. Recall that a face is a two-dimensional region whose interior is connected, and that half-edges are oriented so as to leave the interior of the face to its left. Over the course of the sweep algorithm, an half-edge is created when its first vertex is fixed, and remains active until its second vertex is also fixed. Following classical terminology, notice the first vertex may be the source or the target vertex of the half-edge. In the following, stitching two half-edges should be understood as fixing the next pointer of one half-edge to the second. Faces and Holes. To describe the arrangement, we use sequences of connected half-edges (SCH), such a sequence being called closed if its topology is that of a circle. Two active half-edges associated to arcs in V, with respect to the ordering along V, are termed adjacent if their arcs are adjacent in V, and if they bound the same segment along the meridian M θ . (Out of the four pairs of half-edges associated to two consecutive arcs along V, a single pair corresponds to adjacent half-edges.) A face of the arrangement is represented by a collection of closed SCH. Each such sequence is called a Connected Component of the Boundary or CCB. A CCB is oriented and always induces a contractible region on the sphere S 0 . The set of all CCB describing a face can be split into two categories: one CCB called principal defines the spherical cap containing the face; the remaining ones define holes in the face. See Fig. 4.

(c) (b) S 0 CCB3 CCB2 CCB1 (a) CCB2 CCB1 S 0 S 0 CCB2 CCB1 CCB3

Handling Half-edges

The initialization of V consists of finding the ordering of arcs along M θ for θ = 0 + . To be able to define the faces cut by M 0 + , to each arc in V after initialization, we attribute a pair of opposite half-edges with one null vertex implicitly representing the intersection between M 0 + and the corresponding arc. This collection of half-edges is stored in a list H 0 . At the end of the sweep process, we have a one-to-one correspondence between half-edges of arcs in V and the sequence of half-edges in H 0 , so that a merge can be performed.

To describe the operations underwent by half-edges at a normal event site, consider the arcs involved in the three lists of a normal event site, in the neighborhood of the corresponding event point, say p. Two types of operations need to be performed. The first type, to the left (respectively to the right) of p, consists of stitching each pair of adjacent half-edges of the arcs involved, before (respectively after) updating V. The second type consists of stitching the half-edges of the highest and lowest arcs along V. These operations are straightforward, and the details are left to the reader.

Operations involving polar and bipolar circles are also straightforward: Following Def. 3, half-edges incident to a pole are stitched while rotating around that pole; When handling a bipolar event, from the north to the south pole, all active arcs in V are intersected and their half-edges are updated accordingly. If the bipolar circle goes through a point corresponding to an event site, then the event site is handled at the same time taking into account the half-edges of the bipolar circle.

Building the Faces

Building the faces of the arrangement requires two steps, namely creating CCB and creating faces by joining the CCB. To do so, we resort to two independent union-find algorithms.

Creating a CCB. A SCH becomes a CCB whenever stitching two half-edges creates a topological circle. As the intersection between a CCB and the meridian may feature several connected components (see Fig. 5), we merge SCH using union-find, which requires endowing each half-edge with a pointer to a master half-edge-called the CCB master. Stitching two half-edges is accompanied by a union of their CCB masters (if different), and a CCB appears when two half-edges being stitched already have the same CCB master. Creating a Face. A face consists of a principal CCB and of CCB defining holes. We construct faces using union-find on SCH, which requires a union-find pointer called the face master. While closing a SCH, the resulting CCB is principal if it is its own face master; if not, the CCB becomes a hole of the face referred to by the face master. Moreover while stitching two half-edges, if the face masters are different, then a union operation is done (see Fig. 5). We complete the description by the initialization of the face master.

i 2 i 1 e 1 e 2 s 2 s 1 i 4 i 3 s 4 s 3 i 2 i 1 e 1 e 2 s 2 s 1 i 4 i 3 s 4 s 3 i 2 i 1 e 1 e 2 s 2 s 1 i 4 i 3 s 4 s 3 i 2 i 1 e 1 e 2 s 2 s 1 i 4 i 3 s 4 s 3 M θ M θ M θ
Whenever two arcs are adjacent in V, a pair of active adjacent half-edges bounds the same face. When creating a new SCH, the SCH either contributes to a face in progress (i.e. an adjacent half-edge already handled can be found), or starts a new face. Difficulties arise if one of the half-edge creating the new SCH is adjacent to the north pole. This particular case is carried out by recording the face containing the north pole, denoted F N (M θ )-the subscript θ indicates that this face may evolve during the sweep. If no circle is bipolar or north polar, F N (M θ ) is set once and does not change over the course of the algorithm. Otherwise, the face F N (M θ ) is defined as follows. If the meridian M θ is not tangent to any circle at the pole, F N (M θ ) is defined as the face having the north pole on its boundary and containing an infinitesimal portion of M θ anchored at the north pole. If meridian M θ is tangent to a (bi)polar circle, F N (M θ ) is undefined. But we denote F N (M θ -) (respectively F N (M θ + )) the face containing the north pole for an infinitesimally smaller (respectively larger) value of θ. See Table 1 for initialization and update of F N (M θ ).

Complexity Analysis

To conclude, we give the cost of constructing the HDS storing the arrangement. The analysis consists of counting the number of find and union operations used to maintain the topological data structures i.e. SCH/CCB and faces. Denoting α the inverse of Ackermann's function, recall that the complexity of performing M union-find operations on a N elements set, with M ≥ N , is Θ(M α(M, N )) [START_REF] Tarjan | Data Structures and Network Algorithms[END_REF]. M ≤ U , then α(M, N ) ≤ α(U, U ).

Upon completion of the sweep, the union-find data structure features f connected components corresponding to the f faces of the arrangement. Denoting h the total number of holes, e the number of Position of M θ Operation(s) at north pole Initializing V FN (M θ ) is set to be the face started by the upper hedge of the top arc of V at θ = 0 + . Bipolar or a north FN (M θ + ) is the face started by the inner hedge. polar start event FN (M θ -) is the face described by the outer hedge. North polar end event FN (M θ + ) is the face described by the outer hedge. FN (M - θ ) is the face described by the inner hedge. Bipolar end event FN (M θ + ) is the face started by the outer hedge. FN (M - θ ) is the face described by the inner hedge. Normal start event or FN (M θ ) is set to be the face started by the outer south polar start event hedge(s) (if it is not already done).

Table 1: Initializing and updating the face containing the north pole F N (M θ ), and its relatives

F N (M - θ ), F N (M + θ )
. See text for details. Notice that hedge stands for half-edge.

edges in the arrangement, n the number of circles and n 0 the number of arcs found in V at θ = 0 + , the following theorems are proved in [START_REF] Cazals | Computing the exact arrangement of circles on a sphere, with applications in structural biology[END_REF]:

Theorem 2 Constructing CCB has complexity O((e + n 0 )α(6e + 8n 0 , 6e + 8n 0 )).

Theorem 3 Grouping CCB into faces has complexity O((f

+ h)α(f + h + 20n, f + h + 20n)).

Reporting Inclusion into Balls

Let B i be the ball associated to the sphere S i which generates the intersection circle C i = S 0 ∩ S i . In this section, we describe an algorithm reporting the covering list of each face of the arrangement on S 0 , that is the list of balls containing it.

Enclosing Balls and Unique Circles

In reporting the covering lists, two difficulties are faced. First, a number of degeneracies must be accommodated. Notice that (i) if a sphere is tangent to S 0 the reduces to a point, and if S 0 is covered by the associated ball, so are all the faces of the arrangement (ii) if two spheres intersect S 0 along the same circle, and their balls cover (respectively do not cover) the same part of S 0 , a face covered by one is covered by the other (respectively is not covered by the other). Second, the BO algorithm operates under the assumption that all circles are different. To meet this requirement, we sort the m intersecting spheres using a total ordering returning equality when two spheres intersect S 0 along the same circle. While sorting the spheres, each unique intersection circle is endowed with two special balls: the primary ball, which is associated to the sphere first defining such a circle; and the opposite ball, which is the first to intersect S 0 along the same circle, but is opposite to the primary w.r.t. the plane of the circle. Notice that the opposite ball may not exist. Each primary and opposite ball is attached a list of balls yielding the same intersection circle and covering the same part of S 0 . We call these the lists of friends in the sequel.

Let us get back to the operator used to sort spheres. Two non-great circles are identical iff they have the same center: We first resort to a lexicographic sorting using intersection circle centers. Two great circles (their center being the center of S 0 ) are identical iff the centers of the two spheres defining these circles are aligned with the center of S 0 -the spheres generating these circles belong to the same pencil. Thus, to distinguish great circles, we use the lexicographic order on a canonical vector describing the pencil of spheres yielding a given circle.

Inclusion into Balls

Outline. While running the sweep process, we construct an implicit encoding of the covering lists. In the following, the description focuses on the primary and opposite balls, as the remaining balls are accessible thanks to the lists of friends.

More precisely, the covering lists are represented by a tree, the covering tree. Each node in this tree corresponds to one face of the arrangement, and the edge connecting two nodes corresponds to an arc found in V, stating which primary/opposite ball is added to/removed from the covering list of the father node. Since a face consists of one or several CCB and a CCB is incrementally built from SCH, the principle used to set the covering lists consists of setting one such list for each face master upon creation of the corresponding SCH. This strategy entails two things. First, upon extension at an intersection point of SCH describing a given face, the covering list does not change since the contribution (covering balls) of the corresponding circle has already been taken into account at the creation of the face using a father node. Second, consider the case of the union of two SCH having different face masters: The two face masters must be merged. The consequences are twofold: (i) one of the two covering lists can be discarded; (ii) in the covering tree, the sons of the node suppressed are attached to the node that remains. Initialization. Let LN B θ be the list of balls covering the face F N (M θ ). This list is initialized at θ = 0 + , as indicated in the first section of Table 4.2. The root of the covering tree is precisely LN B θ at θ = 0 + . Updates. First, upon creation of a new face started by a new SCH, its covering list is created, using Observation 1, by updating that of its ancestor in the covering tree. This ancestor corresponds to the face pointed by the face master of the upper half-edge of the highest arc along V involved in the creation of the new SCH. Second, consider the updates of LN B θ . For a circle C, let CLP(C) (respectively CSP(C)) be the function returning, if any, among the primary and opposite balls, the one Covering the Largest (respectively Smallest) Part of S 0 bounded by circle C. Operator + (respectively -) means that the ball is added to (respectively removed from) LN B θ . The cases to be accommodated are listed in the second section of Table 4.2 and illustrated on Fig. 6.

Step At the end of the sweep process, the number of nodes is exactly the number of faces in the arrangement. Moreover, if n stands for the number of intersection circles, the maximum distance between the root and a node is 2n as V never contains more than 2n arcs.

An example of such a tree is presented on Fig. 7. Consider the case of face 6. This face got started at the start event of circle C 4 using its two inner half-edges. The highest one of the pair is the lower half-edge of the upper arc of C 4 . The upper half-edge of that arc describes face 4. Therefore, the node of face 6 is a son of the node of face 4, and since the arc considered is an upper one, the edge between the two nodes indicates that the ball, among primary and opposite, covering the smallest (respectively the largest) part of S 0 according to C 4 must be added (respectively removed).

Complexity Analysis

To analyze the complexity of the implicit encoding, m denotes the number of input balls. The following is proved in [START_REF] Cazals | Computing the exact arrangement of circles on a sphere, with applications in structural biology[END_REF]:

F N (M θ + )/F N (M θ -) is not covered by a ball covering the smallest/largest part of S 0 . N C 1 C 2 a 0 S 0 C 3 C 4 1 3 2 4 7 6 C 5 5 M 0 x y CSP(C 3 ) CLP(C 3 ) + - CSP(C 1 ) CLP(C 1 ) + - CSP(C 5 ) CLP(C 5 ) + - CSP(C 4 ) CLP(C 4 ) + - CSP(C 2 ) CLP(C 2 ) + - CSP(C 1 ) CLP(C 1 ) + - 2 7 4 5 6 1 3 { { LNB at θ = 0 CLP(C 2 ) CSP(C 1 ) CLP(C 3 ) CLP(C 4 ) CLP(C 5 )

Application to Flexible Docking

In this section, we discuss one structural biology problem involving spherical arrangements, namely that of selecting diverse conformational ensembles. The molecular models used are Van der Waals (VdW) models.

Selecting Diverse Conformational Ensembles

Flexible Protein Docking with Conformer Ensembles. Protein-protein complexes are paramount to all biological processes, and predicting the conformation of a complex from the unbound partners is known as the docking problem. Unfortunately, docking is especially challenging, as evidenced in the CAPRI experiment, by the low number of medium and high predictions-as opposed to incorrect and acceptable ones [START_REF] Janin | The Third CAPRI Assessment Meeting Toronto[END_REF]. (The CAPRI experiment is a community wide experiment during which participants are asked to predict the conformation of an unpublished co-crystallized complex from the geometries of the isolated partners. Since the structure of the complex has been solved experimentally, the predictions can be compared to the crystallographic reality.) For such challenging cases, which feature large amplitude motions beyond the time ranges accessible to all-atom molecular dynamics (MD) simulations, discrete ensembles of conformations known as conformer ensembles can be pre-generated and considered simultaneously. (A MD simulation consists of simulating the atomic motions by integrating Newton's equations of motion. Since the integration time-step is of the order of the femtosecond, the simulation times are of the order of tens of nanoseconds for classical systems.) Such ensembles are particularly appropriate when dealing with macromolecular docking, since one wishes to explore the relative position and orientation of the partners, but also their conformational space. In the Monod-Wyman-Changeux interpretation [START_REF] Monod | On the nature of allosteric transitions: a plausible model[END_REF], the unbound proteins are considered as two collections of conformers in thermodynamic equilibrium. When the partners bind, the equilibrium is shifted towards the structure observed in the complex. Implementing this strategy poses two problems, namely generating conformational ensembles, and handling them for docking. We examine this latter problem in the sequel. Diversity from Geometric Optimization. Given a collection of conformers for each partner, flexible docking algorithms manipulating conformational ensembles aim at identifying the two conformations most likely to form a complex. (For a co-crystallized complex, ideally, these conformations should contain at least one conformer as close as possible to the one observed in the crystal structure.) But since the total number of conformations of a given molecule is in general infinite, the question of selecting an ensemble of reasonable size arises. To further specify the problem, assume we are given a large pool C of stable candidate conformations. By stable, we mean that the conformations are energetically favorable, as trying to dock conformations with serious flaws (e.g. steric clashes) would be meaningless. Given such a pre-computed collection C of conformers, we introduce the following two selection problems: ⊲ Problem #1: Find out a subset S of s ≪ n conformers, called the selection, such that the volume occupied by the union of the conformers in the selection is maximized. ⊲ Problem #2: Find out a subset S of s ≪ n conformers, called the selection, such that the area of the VdW surface of the union of the conformers in the selection is maximized.

Notice that the rationale underlying these optimization problems is rather simple: maximizing the volume or the surface exposed by the selection are two ways to ascertain that the conformers are non redundant. For a fixed budget of conformers, this non redundancy warrants the conformational diversity called for by the docking algorithms.

The first problem is intimately related to max-k cover [START_REF] Garey | Computers and Intractability: A Guide to the Theory of NP-Completeness[END_REF][START_REF] Feige | A threshold of ln n for approximating set cover[END_REF], a well-known NP-Complete combinatorial optimization problem. In fact, the problem generalizes max-k cover, since the weight of each cell in the 3D arrangement decomposing the union of balls is not a unit weight but the Euclidean threedimensional volume. The second problem is more tricky, since the surface area of the selection is not, as opposed to the volume, a monotonic quantity: add one conformer, and the surface area may decrease. Because of these difficulties, algorithms solving exactly these problems in time polynomial in both n and s cannot be expected, so that approximation strategies must be sought. An obvious such strategy is the greedy one, which consists of incrementally selecting the s conformers: The conformer selected at the i-th stage is that yielding the best increment. These strategies are examined in [START_REF] Loriot | On the characterization and selection of diverse conformational ensembles[END_REF], where the following theorems are proved: Theorem 5 For Problem #1, the greedy strategy has an approximation guarantee of 1 -(1 -1/s) s > 1 -1/e, which is optimal. (e is the base of the natural logarithm.) Theorem 6 For Problem #2, the greedy approach may have an approximation guarantee as bad as 1/s 2 .

Implementing the Greedy Strategy using Arrangements of Circles on a Sphere. In the following, we focus on Problem #2 for a practical reason: Problem #2 can be solved using arrangements of circles on a sphere, which can be done efficiently as sketched in section 2, using the spherical kernel developed in [START_REF] De Castro | Design of the CGAL spherical kernel and application to arrangements of circles on a sphere[END_REF]; on the other hand, Problem #1 requires decomposing the 3D volume occupied by a collection of balls, and we are not aware of any effective algorithm to perform this task.

To solve Problem #2, assume we have computed for each atom (i) the decomposition of its surface induced by the intersection circles with all atom from all conformers in C, as explained in section 3; (ii) the covering lists for all faces of this arrangement, as explained in section 4.

Implementing the greedy strategy from these pieces of information is elementary: one maintains a priority queue based on the increments each conformer would bring to the selection; upon selecting the top conformer C t , using the covering lists involving balls of C t , one updates the weight of conformers remaining in the queue which feature spherical caps covered by balls of C t . To state the algorithm complexity, define

τ = Ci∈C Sj ∈Ci ( F k ∈Sj s T (F k ) + 1), ( 1
)
where C is the set of all conformers, S j a sphere of a conformer C i , F k a face on the arrangement on sphere S j and s T (F k ) the total number of balls covering F k , as introduced in Thm. 4.

Implementing the priority queue with Fibonacci yields the following theorem, as proved in [START_REF] Loriot | On the characterization and selection of diverse conformational ensembles[END_REF]:

Theorem 7
The greedy selection based on the spherical arrangements of balls of the conformers in C has amortized complexity O(τ + s log n).

Alternatively, one may skip the calculation of the covering lists, in which case the greedy selection can be carried out by recomputing for each candidate conformer C k the surface area of S j-1 ∪ C k , with S j-1 the selection obtained so far. In that case, one has [START_REF] Loriot | On the characterization and selection of diverse conformational ensembles[END_REF]:

Theorem 8
The naive greedy selection of conformers has complexity O(ns 3 ).

From a practical standpoint, and since one has s ≪ n, the best strategy depends on the value of τ given by Eq. ( 1): for sparse arrangements, the algorithm coming along with Thm. 7 should be used; for dense and cluttered arrangements, as τ takes over the cost of running the priority based selection, the algorithm associated to Thm. 8 is the alternative of choice.

Validation

Docking Protocols. To validate the conformer selection strategy based upon surface area maximization, we ran docking simulations between one rigid protein called the ligand (L), and one flexible called the receptor. The receptor itself decomposes into a rigid template (R) and a flexible loop (F). While performing flexible protein docking with conformer ensembles, the strategy consists of using a conformer ensemble for the flexible loop F, this ensemble being selected from a larger pool. Thus, specifying a docking protocol requires specifying the triple R/L/F. To see how, recall that a binary complex used for docking validation features two molecules which have been crystallized under two forms: on their own, i.e. the unbound forms, and in complex i.e. the bound forms. Thus, to specify the rigid parts (R and L), we provide a tag indicating the origin of the partner, namely U for Unbound and B for Bound. To specify the ensemble associated to F, we provide three pieces of informations: (i) the bound/unbound tag which indicates the loop geometry used to generate the pool of conformers (ii) the selection size, and (iii) the algorithm used to select the conformers from this pool. Two selection algorithms were used, namely a standard one based on a hierarchical clustering of the conformers (called HClust) [START_REF] Gordon | Classification[END_REF], and the greedy strategy presented above (called Greedy). For example, F=B-Greedy-10 refers to 10 conformers selected by algorithm Greedy, out of a pool of conformers generated from the Bound structure of the receptor.

To summarize, we report on the following four docking protocols: two using the Bound form of the receptor, namely B/B/B-HClust-10, B/B/B-Greedy-10; and two using the Unbound form of the receptor, namely U/B/B-HClust-10, U/B/B-Greedy-10. In passing, notice that the incentive for using the Bound conformation of the flexible region to generate the conformers is the following: for very flexible systems, the reconstruction of the unbound conformation of the flexible loop from the crystallographic data may not be possible. (If the conformations of the loop changes across the crystallographic units, the signal is not strong enough for the reconstruction to be carried out.) Initial Conditions for a Protocol. For a given protocol, we ran N t docking tests using algorithm ATTRACT [START_REF] Zacharias | Protein-protein docking with a reduced protein model accounting for side-chain flexibility[END_REF]. Each such test corresponds to a specific position and orientation of the ligand with respect to the receptor. Given these initial conditions, ATTRACT performs a sequence of minimizations so as to explore the six degrees of freedom of the ligand. At each stage, the energy of each conformation of the complex is computed, from which a fitness score (between 0 and 1) is attributed to each loop. Upon termination, the loop selected is that having the highest fitness score. An assessment of the quality of the proposed complex is then based upon two figures: (i) the interaction potential energy E of the complex (ii) the i-rmsd of the atoms of the ligand. (The i-rmsd is the standard deviation on the atomic positions of the alpha carbons of the extended interface of two conformations: that output by the docking program, and that found in the co-crystallized complex. The extended interface is defined as the set of residues of the ligand having an atom of the receptor within 5 Å in the X-ray structure.) For the N t tests associated to a given protocol, the plot of the pairs (E,i-rmsd) defines the energy landscape of the docking experiment. Thus, a conformer ensemble is satisfactory if the landscape features at least one conformer yielding a large number of points (E,i-rmsd) next to the bottom left corner of the energy landscape. In particular, following the CAPRI assessment rules, a medium (respectively acceptable) prediction corresponds to a i-rmsd smaller than 5 Å (respectively 10 Å) [START_REF] Lensink | Docking and scoring protein complexes: CAPRI 3rd Edition[END_REF]. Practically, we represent an energy landscape using buckets. For a given bucket B i and conformer C j , let s i,j be the number of times conformer C j yields a complex whose energy and i-rmsd fall in bucket B i . (Notice that i,j s i,j = N t .) Finally, for a given bucket B i , denote l i the index of the conformer that yields the largest value of s i,j , and let r i = j=1,...,n;j =li s i,j . In bucket B i we display the score s i,li , together with r i when r i = 0. The color used to do so is that associated to conformer l i -one color per conformer.

System Tested and Assessment. As an illustration, we report docking energy landscapes for complex 1BTH, whose receptor is a thrombin mutant, the ligand being the pancreatic trypsin inhibitor. The flexible region is a loop of 10 amino acids on the receptor. A pool of n = 500 conformers was generated using Loopy [START_REF] Xiang | Evaluating conformational free energies: The colony energy and its application to the problem of loop prediction[END_REF], from which s = 10 were selected using either the HClust or Greedy algorithms aforementioned. For each selection method, a total of N t = 30, 000 docking tests were run using the same 10 selected conformers. Out of the N t tests, Fig. 8 and Fig. 9 present plots of the results corresponding to a i-rmsd ≤ 15 Å and an energy in the range [-30, 0] units.

Fig. 8 shows that the conformer selection provided by Greedy is of higher quality than that provided by HClust. To see why, consider conformations of the complex with i-rmsd less than 3 Å and energy below -20 units: the number of hits for the top scoring conformer are respectively of about 40 against 160 for HClust and Greedy respectively. On Fig. 9, which corresponds to the unbound protocol, the results are more mitigated. Consider conformers of energy below -15 units, and i-rmsd in the range [START_REF] Boissonnat | Algorithmic Geometry[END_REF][START_REF] Chazal | Stability of boundary measure[END_REF]. Two of the ten conformers selected by HClust are well represented, while one conformer selected by Greedy features a high concentration in-between -17 and -18 energy units. In such cases, a finer inspection of the structures is in order to select the best complex(es) [START_REF] London | Funnel Hunting in a Rough Terrain: Learning and Discriminating Native Energy Funnels[END_REF].

As a conclusion, on this system, using diverse conformational ensembles increases the chance of ending up with a putative complex with low i-rmsd and low energy with respect to the co-crystallized complex. This trend is general, and a more detailed discussion comparing protein-protein complex will appear elsewhere.

Conclusion

This paper presents an integrated framework based upon the generalization of the Bentley-Ottmann algorithm to the spherical setting. This framework accommodates the calculation of the exact arrangement of circles on a sphere, as well as the construction of the corresponding half-edge data structure on the fly. Moreover, assuming that each circle comes from the intersection between the central sphere and neighboring spheres, each coming with an accompanying ball, we explain how to efficiently compute the covering lists of faces of the arrangement, i.e. the lists of balls that contain it. In passing, we notice that while exactness of the arrangement is not a goal per se, it is one way to achieve robustness, at a modest computational overhead as testified in a companion paper. From an application perspective, we use the framework in structural biology, namely for the problem of flexible docking. Application-wise, a number of modeling situations should benefit from the pieces of information encoded in the covering lists. From an algorithmic perspective, as our algorithm is dedicated to circles on a sphere, the question of coming up with more general algorithms featuring the same performances deserves investigations.

Figure 1 :

 1 Figure 1: Modeling atomic multi-body contacts-PDB code 1ACB. Left: A protein-protein complex with two partners; Right: Cumulative surface area, over all atoms of the complex, as a function of the multiplicity of the cells. Multiplicity 0 corresponds to the SAS surface, and accounts for a mere 4.5% of the cumulative area. See text for details.

Figure 2 :

 2 Figure 2: Conformers and flexible protein -protein docking -PDB code 1BTH. Left: Side view of rigid template (cartoons), together with the backbone of a collection of 20 conformers (polylines); Right: Top view of five shallowly intersecting conformers (VdW representation).

Definition 3

 3 Event sites of different types with same θ value are ordered as follows: e Pep < e B < e N < e Psp where a < b means a occurs before b; e B and e N stand for respectively bipolar and normal event site; e Psp and e Pep stand for respectively start and end polar event sites. Moreover, among polar start (respectively end) event sites at the same pole, the one whose associated circle is of largest (respectively smallest) radius occurs first.

Figure 3 :

 3 Figure 3: Qualifying arcs and half-edges: upper and lower arcs are respectively denoted as A and A; inner and outer half-edges are respectively represented using dotted and dashed lines; (a) normal circle C n and its upper and lower arcs; (b) north threaded circle C t1 defining a lower arc, and south threaded circle C t2 defining an upper arc; (c) north polar circle C p1 defining a lower arc, and south polar circle C p2 defining an upper arc; (d) bipolar circle C bp inducing half-edges (no arc is defined). The θ-coordinate of the start point of C bp lies in (0; π].

Figure 4 :

 4 Figure 4: Defining a face with a set of CCB. Each arrow represents a CCB. In the four cases, only one face is defined by the CCB represented.

Figure 5 :

 5 Figure5: Illustration of Union-Find: one color per face, one line-style per CCB. Two SCH get created at events corresponding to i 1 and i 2 , together with two faces. When the meridian reaches start point s 1 , (and s 2 ) a new SCH is created and it contributes to the definition of an existing face. The handling of the event corresponding to the end point e 1 results on the one hand, and on the other hand to the union of two CCB masters and to the union of two face masters: The three remaining CCB describe a unique face. Thanks to this operation, the CCB started with points s 3 and s 4 also describe that face. Finally, the unions of CCB masters at events associated to i 3 and i 4 allows to detect the closure of the CCB at point e 2 .

Table 2 :

 2 Type of event Action(s) on LNB θ Filtering 1a S ∩ S0 = a point and S0 ⊂ B + B spheres 1b B covers the largest part of S0 + B Classifying 2a North polar circle with θE < θS + CSP(C) and -CLP(C) circles 2b North threaded circle + CSP(C) and -CLP(C) Handling 3a North (bi)polar circle Starting + CSP(C) and -CLP(C) events 3b North (bi)polar circle Ending + CLP(C) and -CSP(C) Updating the balls in LN B θ : initialization and updates. The ball/sphere/circle processed is denoted B/S/C. See text for the definition of functions CLP and CSP, and Fig. 6 for the details.

Theorem 4

 4 Computing the covering tree has complexity O(f + m) and O(f + m) space. Denote s T (F ) the total number of balls covering face F . Constructing the covering lists for all faces of the arrangement requires O( F ∈f aces s T (F ) + f ) time.

Figure 6 :

 6 Figure 6: Updating the balls of LN B θ : the six cases of Table4.2. 1(a) S 0 is covered by a ball whose sphere is tangent to S 0 ; 1(b) A sphere intersects S 0 , and the associated ball covers the largest part of S 0 ; 2(a) Meridian M 0 intersects a polar circle. The ball covering the smallest part of S 0 covers F N (M 0 ); 2(b) A north threaded circle. The north pole is always (respectively never) covered by a ball covering the smallest (respectively largest) part of S 0 ; 3(a) The dashed meridian M θ stopped at a (bi)polar start event:F N (M θ + )/F N (M θ -)is covered by a ball covering the smallest/largest part of S 0 ; 3(b) The dashed meridian M θ stopped at a (bi)polar end event: F N (M θ + )/F N (M θ -) is not covered by a ball covering the smallest/largest part of S 0 .

Figure 7 :

 7 Figure 7: Implicit encoding of covering lists. Circle C 1 is north threaded, C 2 is bipolar, while C 3 , C 4 and C 5 are normal circles. Faces of the arrangement depicted on the left are identified using numbers from 1 to 7. Notice that node associated to face 2 reflects the modification of LN B θ induced by bipolar circle C 2 .

Figure 8 :

 8 Figure 8: 1BTH: binning the docking tests using the Bound form of the receptor. See text for details. 15

Figure 9 :

 9 Figure 9: 1BTH: binning the docking tests using the Unbound form of the receptor. See text for details. 16

  .

This work is partially supported by the IST Programme of the 6th Framework Programme of the EU as a STREP (FET Open Scheme) Project under Contract No IST-006413 -ACS (Algorithms for Complex Shapes with certified topology and numerics)http://acs.cs.rug.nl/