Chaoyu Quan

Benjamin Stamm

Mathematical Analysis and Calculation of Molecular Surfaces

Keywords: Molecular surface, solvent excluded surface, solvation models, mathematical description, constructive algorithm

In this article we derive a complete characterization of the Solvent Excluded Surface (SES) for molecular systems including a complete characterization of singularities of the surface. The theory is based on an implicit representation of the SES, which, in turn, is based on the signed distance function to the Solvent Accessible Surface (SAS). All proofs are constructive so that the theory allows for efficient algorithms in order to compute the area of the SES and the volume of the SES-cavity, or to visualize the surface. Further, we propose to refine the notion of SAS and SES in order to take inner holes in a solute molecule into account or not.

Introduction

The majority of chemically relevant reactions take place in the liquid phase and the effect of the environment (solvent) is important and should be considered in various chemical computations. In consequence, a continuum solvation model is a model in which the effect of the solvent molecules on the solute are described by a continuous model [START_REF] Tomasi | Quantum mechanical continuum solvation models[END_REF]. In continuum solvation continuum models, the notion of molecular cavity and molecular surface is a fundamental part of the model. The molecular cavity occupies the space of the solute molecule where a solvent molecule cannot be present and the molecular surface, the boundary of the corresponding cavity, builds the interface between the solute and the solvent respectively between the continuum and the atomistic description of the physical model. A precise understanding of the nature of the surface is therefore essential for the coupled model and in consequence for running numerical computations. The Van der Waals (VdW) surface, the Solvent Accessible Surface (SAS) and the Solvent Excluded Surface (SES) are well-established concepts. The VdW surface is more generally used in chemical calculations, such as in the recent developments [START_REF] Cancès | Domain decomposition for implicit solvation models[END_REF][START_REF] Lipparini | Fast domain decomposition algorithm for continuum solvation models: Energy and first derivatives[END_REF] for example, of numerical approximations to the COnductor-like Screening MOdel (COSMO) due to the simplicity of the cavity. Since the VdW surface is the topological boundary of the union of spheres, the geometric features are therefore easier to understand. However, the SES, which is considered to be a more precise description of the cavity, is more complicated and its analytical characterization remains unsatisfying despite a large number of contributions in literature.

Previous Work

In quantum chemistry, atoms of a molecule can be represented by VdW-balls with VdW-radii obtained from experiments [START_REF] Rappé | Uff, a full periodic table force field for molecular mechanics and molecular dynamics simulations[END_REF]. The VdW surface of a solute molecule is consequently defined as the topological boundary of the union of all VdW-balls. For a given solute molecule, its SAS and the corresponding SES were first introduced by Lee & Richards in the 1970s [START_REF] Lee | The interpretation of protein structures: estimation of static accessibility[END_REF][START_REF] Richards | Areas, volumes, packing and protein structure[END_REF], where the solvent molecules surrounding a solute molecule are reduced to spherical probes [START_REF] Tomasi | Quantum mechanical continuum solvation models[END_REF]. The SES is also called "the smooth molecular surface" or "the Connolly surface", due to Connolly's fundamental work [START_REF] Connolly | Analytical molecular surface calculation[END_REF]. He has divided the SES into three types of patches: convex spherical patches, saddle-shaped toroidal patches and concave spherical triangles. But the self-intersection among different patches in this division often causes singularities despite that the whole SES is smooth almost everywhere. This singularity problem has led to difficulty in many associated works with the SES, for example, failure of SES meshing algorithms and imprecise calculation of molecular areas or volumes, or has been circumvented by approximate techniques [START_REF] Chen | Proteinvolume: calculating molecular van der waals and void volumes in proteins[END_REF]. In 1996, Michel Sanner treated some special singularity cases in his MSMS (Michel Sanner's Molecular Surface) software for meshing molecular surfaces [START_REF] Sanner | Reduced surface: an efficient way to compute molecular surfaces[END_REF]. However, to our knowledge, the complete characterization of the singularities of the SES remains unsolved.

Contribution

In this paper, we will characterize the above molecular surfaces with implicit functions, as well as provide explicit formulas to compute analytically the area of molecular surfaces and the volume of molecular cavities. We first propose a method to compute the signed distance function to the SAS, based on three equivalence statements which also induce a new partition of R 3 . In consequence, a computable implicit function of the corresponding SES is given from the relatively simple relationship between the SES and the SAS. Furthermore, we will redefine different types of SES patches mathematically so that the singularities will be characterized explicitly. Besides, by applying the Gauss-Bonnet theorem [START_REF] Carmo | Differential geometry of curves and surfaces[END_REF] and the Gauss-Green theorem [START_REF] Federer | The gauss-green theorem[END_REF], we succeed to calculate analytically all the molecular areas and volumes, in particular for the SES. These quantities are thought to be useful in protein modeling, such as describing the hydration effects [START_REF] Richmond | Solvent accessible surface area and excluded volume in proteins: Analytical equations for overlapping spheres and implications for the hydrophobic effect[END_REF][START_REF] Buša | Arvo: a fortran package for computing the solvent accessible surface area and the excluded volume of overlapping spheres via analytic equations[END_REF].

In addition, we will refine the notion of SAS and SES by considering the possible inner holes in the solute molecule yielding the notions of the complete SAS (cSAS) and the corresponding complete SES (cSES). To distinguish them, we call respectively the previous SAS and the previous SES as the exterior SAS (eSAS) and the exterior SES (eSES). A method with binary tree to construct all these new molecular surfaces will also be proposed in this paper in order to provide a computationally efficient method.

Outline

We first introduce the concepts of implicit surfaces in the second section and the implicit functions of molecular surfaces are given in the third section. In the fourth section, we present two more precise definitions about the SAS, either by taking the inner holes of the solute molecule considered into account or not. Then, based on three equivalence statements that are developed, we propose a computable method to calculate the signed distance function from any point to the SAS analytically. In this process, a new Voronoi-type diagram for the SAS-cavity is given which allows us to calculate analytically the area of the SAS and the volume inside the SAS. In the fifth section, a computable implicit function of the SES is deduced directly from the signed distance function to the SAS and according to the new Voronoitype diagram, all SES-singularities are characterized. Still within this section, the formulas of calculating the area of the SES and the volume inside the SES will be provided. In the sixth section, we explain how to construct the SAS (cSAS and eSAS) and the SES (cSES and eSES) for a given solute molecule considering the possible inner holes. Numerical results are illustrated in the seventh section and finally, we provide some conclusions of this article in the last section.

Introduction to Implicit Surfaces

We start with presenting the definition of implicit surfaces [START_REF] Velho | Implicit objects in computer graphics[END_REF]. In a very general context, a subset O ⊂ R n is called an implicit object if there exists a real-valued function f :

U → R k with O ⊂ U ⊂ R n , and a subset V ⊂ R k , such that O = f -1 (V ). That is, O = {p ∈ U : f (p) ∈ V }.
The above definition of an implicit object is broad enough to include a large family of subsets of the space. In this paper, we consider the simple case where U = R 3 , V = {0} and f : R 3 → R is a real-valued function. In consequence, an implicit object is represented as a zero-level set O = f -1 (0), which is also called an implicit surf ace in R 3 , and the function f is called an implicit f unction of the implicit surface. Notice that there are various implicit functions to represent one surface in the form of a zero-level set.

The signed distance function f S of a closed bounded oriented surface S in R 3 , determines the distance from a given point p ∈ R n to the surface S, with the sign determined by whether p lies inside S or not. That is to say,

f S (p) =      -inf x∈S p -x if p lies inside S, inf x∈S p -x if p lies outside S. (2.1)
This is naturally an implicit function of S.

Implicit Molecular Surfaces

In quantum chemistry, atoms of a molecule are represented by VdW-balls with VdW radii which are experimentally fitted, given the underlying chemical element [START_REF] Rappé | Uff, a full periodic table force field for molecular mechanics and molecular dynamics simulations[END_REF]. In consequence and mathematically speaking, the VdW surface is defined as the topological boundary of the union of all VdW-balls. Besides, the SAS of a solute molecule is defined by the center of an idealized spherical probe rolling over the solute molecule, that is, the surface enclosing the region in which the center of a spherical probe can not enter. Finally, the SES is defined by the same spherical probe rolling over the solute VdW-cavity, that is, the surface enclosing the region in which a spherical probe can access. In other words, the SES is the boundary of the union of all spherical probes that do not intersect the VdW-balls of the solute molecule, see Figure 1 for a graphical illustration.

We denote by M the number of atoms in a solute molecule, by c i ∈ R 3 and r i ∈ R + the center and the radius of the i-th VdW atom. The open ball with center c i and radius r i is called the i-th VdW-ball. The Van der Waals surface can consequently be represented as an implicit surface f -1 vdw (0) with the following implicit function:

f vdw (p) = min i=1,...,M { p -c i 2 -r i }, ∀p ∈ R 3 . (3.2)
Similarly, the open ball with center c i and radius r i +r p is called the i-th SAS-ball denoted by B i , where r p is the radius of the idealized spherical probe. Furthermore, we denote by S i the i-th SAS-sphere corresponding to B i , that is, S i = ∂B i . Similar to the VdW surface, the SAS can be represented as an implicit surface f -1 sas (0) with the following implicit function:

f sas (p) = f vdw (p) -r p = min i=1,...,M { p -c i 2 -r i -r p }, ∀p ∈ R 3 . (3.3)
We notice that the above implicit function of the SAS is simple to compute. It seems nevertheless hopeless to us to further obtain an implicit function of the SES if constructing upon this simple implicit function f sas (p) which is not a distance function. On the other hand, having the signed distance function, see (2.1), at hand would allow the construction of an implicit function for the SES due to the geometrical relationship between the SAS and the SES. Indeed, according to the fact that any point on the SES has signed distance -r p to the SAS, an implicit function of the SES is obtained directly as:

f ses (p) = f sas (p) + r p , (3.4)
which motivates the choice of using the signed distance function to represent the SAS. From the above formula, the SES can be represented by a level set f -1 sas (-r p ), associated with the signed distance function f sas to the SAS. Therefore, the key point becomes how to compute the signed distance f sas (p) from a point p ∈ R 3 to the SAS. Generally speaking, given a general surface S ⊂ R 3 and any arbitrary point p ∈ R 3 , it is difficult to compute the signed distance from p to S. However, considering that the SAS is a special surface formed by the union of SAS-spheres, this computation can be done analytically. We state a remark about another implicit function to characterize the SES, proposed by Pomelli and Tomasi [START_REF] Pomelli | Defpol: New procedure to build molecular surfaces and its use in continuum solvation methods[END_REF]. In [START_REF] Parulek | Implicit representation of molecular surfaces[END_REF], this function can be written as:

f ses (p) = min 1≤i<j<k≤M f ijk (p), ∀p ∈ R 3 , (3.5)
where f ijk represents the signed distance function to the SES of the i-th, j-th and k-th VdW atom. However, this representation might fail sometimes, see two representative 2D examples in Figure 2. Indeed, the formula (3.5) for each molecule in Figure 2 can be rewritten as:

f ses (p) = min 1≤i<j≤3 f i,j (p), ∀p ∈ R 2 , (3.6)
where f i,j represents the signed distance function to the SES of the i-th and the j-th VdW atom. However, each molecular cavity defined by {p ∈ R 2 : f ses (p) ≤ 0} has excluded the region in grey inside the real SES. Further, the region enclosed by the Van der Waals surface is called the VdWcavity, that is, any point p in the VdW-cavity satisfies f vdw (p) ≤ 0. More generally, we call the region enclosed by a molecular surface as the corresponding molecular cavity. In consequence, the region enclosed by the SAS is called the SAS-cavity, and the region enclosed by the SES is called the SES-cavity. Similarly, any point p in the SAS-cavity satisfies f sas (p) ≤ 0, and any point p in the SES-cavity satisfies f ses (p) ≤ 0.

The Solvent Accessible Surface

In the framework of continuum solvation models, using the VdW-cavity as the solvent molecular cavity has the characteristic that its definition does not depend in any way on characteristics of the solvent. In other words, the above-mentioned VdW surface has ignored the size and shape of the solvent molecules, while the definition of the Solvent Accessible Surface includes some of these characteristics. In the following, we first provide two more precise mathematical definitions of the SAS considering possible inner holes of a solute molecule. After that, we will provide a formula for the signed distance function to the SAS, which is indeed based on three equivalence statements, providing explicitly a closest point on the SAS to any point in R 3 . In this process, a new Voronoi-type diagram will be proposed to make a partition of the space R 3 , which, in turn, will also be used to calculate the exact volume of the SAS-cavity.

Mathematical Definitions

In [START_REF] Lee | The interpretation of protein structures: estimation of static accessibility[END_REF], the Solvent Accessible Surface is defined by the set of the centers of the spherical probe when rolling over the VdW surface of the molecule. At first glance, one could think that it can equivalently be seen as the topological boundary of the union of all SAS-balls of the molecule. However, we notice that there might exist some inner holes inside the molecules where a solvent molecule can not be present. In consequence, the SAS may or may not be composed of several separate surfaces, Figure 4: This is a 2D schematic of different molecular surfaces, including the VdW surface, the cSAS, the eSAS, the cSES and the eSES. The yellow discs denote the probes, representing solvent molecules, while the blue discs denote the VdW-balls. SAS 1 is the trace of the probe center when a probe rolls over the exterior of the VdW-balls, while SAS 2 is the trace of the probe center when a probe rolls over the inner holes of the VdW-balls. SES 1 and SES 2 respectively denote the corresponding solvent excluded surfaces to SAS 1 and SAS 2 . The cSAS is the union of SAS 1 and SAS 2 , while the eSAS is just SAS 1 . Similarly, the cSES is the union of SES 1 and SES 2 , while the eSES is just SES 1 . Finally, x i0 is a closest point on the SAS to a given point p.

see Figure 3 for a graphical illustration. This inspires us to propose two more precise surfaces: the complete Solvent Accessible Surface (cSAS) defined simply as the boundary of the union of the SAS-balls, and the exterior Solvent Accessible Surface (eSAS) defined as the outmost surface obtained when a probe rolls over the exterior of the molecule, see Figure 4. In the case where there are no interior holes inside the molecule, the cSAS and the eSAS will coincide. We make a convention that the SAS refers to both the cSAS and the eSAS in a general context.

Since both the cSAS and the eSAS are two closed sets, there exists a closest point on the SAS to any given point p ∈ R 3 , which is denoted by x p sas . Thus, the signed distance function f sas can be written as: In the above formula, x p sas depends on p. When p lies on the SAS, p coincides with x p sas and f sas (p) = 0. It remains therefore to find a closest point x p sas on the SAS to the given point p ∈ R 3 . Note that there might exist more than one closest point on the SAS to the point p and x p sas is chosen as one of them. In the context, this is not an easy task. In the following, for the particular case of the SAS, we propose a way to calculate analytically a closest point on the SAS to a given point p, based on three equivalence statements.

f sas (p) = -p -x p

Equivalence Statements

According to the definitions of the cSAS and the eSAS, these two molecular surfaces are both composed of three types of parts: open spherical patches, open circular arcs and intersection points (formed by the intersection of at least three SAS-spheres), see Figure 3. Note that an SAS intersection point can in theory be formed by the intersection of more than three SAS-spheres. However, these cases can be divided into multiple triplets of SAS-spheres for simplicity as mentioned in [START_REF] Kauker | Rendering molecular surfaces using order-independent transparency[END_REF]. In this spirit, we make an assumption that all SAS intersection points are formed by the intersection of three SAS-spheres. Furthermore, we assume that any SAS-ball is not included by any other (otherwise, the inner SAS-ball can be ignored) In the following analysis.

For an SAS, denote by m 1 the number of the SAS spherical patches, by m 2 the number of the SAS circular arcs, and by m 3 the number of the SAS intersection points. Then, denote by P m the m-th SAS spherical patch on the SAS where m = 1, . . . , m 1 . Denote by l m the m-th SAS circular arc on the SAS where m = 1, . . . , m 2 . Denote by x m the m-th SAS intersection point on the SAS where m = 1, . . . , m 3 . Furthermore, denote by I the set of all SAS intersection points written as:

I = {x m : m = 1, . . . , m 3 } = {x ∈ SAS : ∃ 1 ≤ i < j < k ≤ M, s.t. x ∈ S i ∩S j ∩S k }.
With the above notations, we consider to calculate a closest point on the SAS to a given point p in the case when p lies outside the SAS (the cSAS or the eSAS).

Lemma 4.1. Let the point p lie outside the SAS (the cSAS or the eSAS), i.e. p -

c i ≥ r i + r p , ∀1 ≤ i ≤ M . Further, let i 0 ∈ {1, . . . , M } be such that p -c i 0 -(r i 0 + r p ) = min 1≤i≤M { p -c i -(r i + r p )}. (4.8)
Then, the point

x i 0 = c i 0 + (r i 0 + r p ) p -c i 0 p -c i 0 , (4.9)
is the closest point to p on the SAS and

f sas (p) = p -x i 0 = min 1≤i≤M { p -c i -(r i + r p )}. (4.10)
Proof. The proof involves basis geometric manipulations and is left to the reader.

So far, we have discussed the case when p lies outside the SAS (both the cSAS and the eSAS), which is not too difficult to deal with. Next, we need to consider the case where p lies inside the (complete or exterior) SAS-cavity, to obtain the signed distance function f sas (p) from any point p ∈ R 3 to the SAS. The following analysis can be applied to both the cSAS and the eSAS. This problem is handled inversely, in the sense that we will determine the region in the molecular cavity for an arbitrary given point x sas ∈ SAS, such that x sas is a closest point to any point in this region.

To do this, we first define a mapping R : X → Y , where X is a subset of the SAS, and Y = R(X) is the region in the SAS-cavity, such that there exists a closest point in X on the SAS to any point in Y . That is, Y = {y : y lies in the SAS-cavity and ∃ x y sas ∈ X s.t. x y sas is a closest point to y}.

In the following, we propose three equivalence statements between a point x sas on the SAS and the corresponding region R(x sas ) for three cases where x sas lies respectively on the three different types of the SAS. We recall first, however, a useful inequality between two signed distance functions to two surfaces.

(p) ≤ f S (p), ∀p ∈ R 3 .
With the above proposition, we propose first a result which connects a point on an SAS spherical patch with a closed line segment in the SAS-cavity.

Theorem 4.1. Assume that p ∈ R 3 is a point in the SAS-cavity and x sas is a point on the SAS. If x sas is on an SAS spherical patch P m on the i-th SAS-ball S i , then x sas is a closest point on the SAS to p if and only if p lies on the closed line segment

[c i , x sas ]. That is, R(x sas ) = [c i , x sas ].
Further, the closest point x sas on the SAS is unique if and only if p = c i . If p = c i , then any point on P m is a closest point to p.

Proof. First suppose that x sas is a closest point to p with x sas ∈ P m ⊂ S i . Since P m is an open set, take a small enough neighborhood V of x sas such that V ⊂ P m , see Figure 5, and since x sas is a closest point on the SAS to p, we have

p -x sas ≤ p -x , ∀x ∈ V, (4.11)
which yields that the vector from x sas to p is perpendicular to any vector in the tangent plane of S i at x sas , thus the vector from p to x sas is the normal vector at x sas of S i . In consequence, p must lie on the line passing through x sas and the center c i On the other hand, suppose that x sas ∈ P m ⊂ S i , and p ∈ [c i , x sas ]. In consequence, x sas is obviously a closest point on the sphere S i to p, see Figure 5. The signed distance function f S i (p) is equal to -p -x sas . Notice that the cavity inside S i , i.e. B i , is contained in the SAS-cavity. We can then use Proposition 4.1 by taking S as S i and S as the SAS, to obtain

f sas (p) ≤ f S i (p) = -p -x sas . Therefore, we have -p -x ≤ f sas (p) ≤ -p -x sas , ∀x ∈ SAS. (4.12)
That is, p -x ≥ p -x sas , ∀x ∈ SAS, which means that x sas is a closest point on the SAS to p. Finally, assume that p ∈ [c i , x sas ]. If p = c i , then any point on P m is a closest point to p because the distance is uniformly c i -x sas = r i + r p . If p = c i , then the open ball B r (p) with r = p -x sas < r i + r p is included in B i and ∂B r (p) ∩ S i = {x sas }, which implies that x sas is the unique closest point to p.

Next, we propose another equivalence statement, which connects a point on an SAS circular arc with a closed triangle in the SAS-cavity. Theorem 4.2. Assume that p ∈ R 3 is a point in the SAS-cavity and x sas is a point on the SAS. If x sas is on an SAS circular arc l m associated with S i and S j , then x sas is a closest point on the SAS to p if and only if p lies in the closed triangle x sas c i c j defined by three vertices x sas , c i and c j . That is, R(x sas ) = x sas c i c j . Further, the closest point x sas on the SAS is unique if and only if p does not belong to the edge

[c i , c j ]. If p ∈ [c i , c j ], any point on l m is a closest point to p.
Proof. First suppose that x sas is a closest point to p with x sas ∈ l m ⊂ S i ∩ S j . Since l m is an open circle arc on the circle S i ∩ S j , take a small enough neighboring curve γ 0 of x sas with γ 0 ⊂ l m , see Figure 6. Since x sas is a closest point to p, we have

p -x sas ≤ p -x , ∀x ∈ γ 0 . (4.13)
Denote by p the projection of p onto the plane where l m lies. By substituting p -x sas 2 = p -p 2 + p -x sas 2 and p -x 2 = p -p 2 + p -x 2 into the inequality (4.13), we obtain that

p -x sas ≤ p -x , ∀x ∈ γ 0 , (4.14)
which yields that the vector from x sas to p is perpendicular to the tangent vector of γ 0 at x sas . In consequence, p has to lie on the ray Ox sas starting from O and passing through x sas , where O is the center of the circular arc l m . Thus, p must lie on the closed half plane Π defined by the three points x sas , c i and c j and whose boarder is the line passing through c i and c j . The closed half plane Π contains the ray Ox sas and is perpendicular to the plane where l m lies. We then take another two neighboring curves γ 1 and γ 2 of x sas , with γ 1 ⊂ Π ∩ S i ∩ SAS and γ 2 ⊂ Π ∩ S j ∩ SAS, see Figure 6. In this case, γ 1 has one closest endpoint x sas and is open at the other end, which is the same for γ 2 . From the assumption that x sas is a closest point on the SAS to p, we have the following inequality

p -x sas ≤ p -x , ∀x ∈ γ 1 ∪ γ 2 . (4.15)
This yields that p ∈ x sas c i c j , where x sas c i c j is the closed triangle on Π with three vertices x sas , c i and c j . Otherwise, we can find a point x ∈ (γ 1 ∪ γ 2 )\x sas strictly closer to p than x sas , which contradicts the assumption. On the other hand, suppose that p ∈ x sas c i c j . In consequence, it is not difficult to obtain that x sas is the closest point to p on ∂(B i ∪ B j ), where B i and B j are the corresponding SAS-balls corresponding to S i and S j as mentioned above. Similarly, we know that the signed distance function f ∂(B i ∪B j ) (p) is equal to -p -x sas , and notice that B i ∪ B j is contained in the SAS-cavity. We can use again Proposition 4.1, by taking S as ∂(B i ∪ B j ) and S as the SAS, to obtain

f sas (p) ≤ f ∂(B i ∪B j ) (p) = -p -x sas . Therefore, we have -p -x ≤ f sas (p) ≤ -p -x sas , ∀x ∈ SAS. (4.16)
That is, p -x ≥ p -x sas , ∀x ∈ SAS, which means that x sas is a closest point on the SAS to p. The point x sas on l m is a closest point on the SAS to p. γ 0 is a small neighborhood of x sas on the circular arc l m , whereas γ 1 and γ 2 represent two small neighboring curves of x sas on the SAS, with x sas as the endpoints, γ 1 ⊂ Π ∩ S i ∩ SAS, and

γ 2 ⊂ Π ∩ S j ∩ SAS. If p ∈ [c i , c j ], then any point on l m is a closest point to p since the distance from any point on l m to p is constant. If p ∈ x sas c i c j \[c i , c j ], then the open ball B r (p) with r = p -x sas < r i + r p is included in B i ∪ B j and ∂B r (p) ∩ ∂(B i ∪ B j ) = {x sas }, which implies that x sas is the unique closest point to p.
By mapping with R a whole SAS spherical patch P m to the SAS-cavity, we obtain a spherical sector R(P m ) in the SAS-cavity with cap P m , center c i and radius r i + r p , see Figure 5. Similarly, by mapping with R a whole SAS circular arc l m to the SAScavity, we obtain a double-cone region R(l m ) in the SAS-cavity with the circular sector corresponding to l m as base, c i and c j as vertices, see Figure 6.

Removing now the above-mentioned spherical sectors and double-cone regions from the SAS-cavity, the closed hull of the remaining region, denoted by T , is consequently a collection of closed separate polyhedrons, see an example of three atoms in Figure 7. Each polyhedron is denoted by T n with index n and thus T = T n . Considering now the last case where x sas is an SAS intersection point x m , we have a third equivalence statement as a corollary of the previous two equivalence statements. In other words, for an arbitrary intersection point x m ∈ I, we have the formula of its corresponding region in the SAS-cavity as following:

R(x m ) = T ∩ {p : p -x m ≤ p -x , ∀x = x m , x ∈ I}.
(4.17)

It is not difficult to find that R(x m ) is a closed polyhedron. At the same time, mapping I into the SAS-cavity with R, we obtain R(I) = T . With the three equivalence statements as well as the above-defined map R, we obtain in fact a non-overlapping decomposition of the SAS-cavity, including spherical sectors, double-cone regions and polyhedrons. It should be emphasized that given a point p ∈ R 3 contained in the SAS-cavity and known the region where p lies, we can then calculate a closest point x p sas on the SAS to p according to this decomposition. The signed distance function f sas (p) can therefore be calculated analytically by the formula 4.7, which, in turn, will ultimately provide an implicit function of the SES. In the next subsection, we will investigate further in the partition of the SAS-cavity based on the above equivalence statements, in order to compute efficiently the area of the SAS and the volume of the SAS-cavity.

New Voronoi-type Diagram

The above non-overlapping decomposition of the SAS-cavity can also be seen as a new Voronoi-type diagram, which will be presented comparing with the well-known Voronoi Diagram and the Power Diagram recalled in the following.

a.. Voronoi Diagram

The Voronoi diagram [START_REF] Drysdale | Generalized Voronoi Diagrams and Geometric Searching[END_REF] was initially a partition of a 2D plane into regions based on the distance to points in a specific subset of the plane. In the general case, for an Euclidean subspace X ⊂ R n endowed with a distance function d and a tuple of nonempty subsets (A i ) i∈K in X, the Voronoi region R i associated with A i is the set of all points in X whose distance to A i is not greater than their distance to any other set A j for j = i. In other words, with the distance function between a point x and a set A i defined as

d V (x, A i ) = inf y∈X {d(x, y) | y ∈ A i },
the formula of the Voronoi region R i is then given by:

R i = {x ∈ X | d V (x, A i ) ≤ d V (x, A j ) , ∀ j = i}. (4.18)
Most commonly, each subset A i is taken as a point and its corresponding Voronoi region R i is consequently a polyhedron, see an example of three points in R 2 in Figure 8.

b.. Power Diagram

In computational geometry, the power diagram [START_REF] Aurenhammer | Power diagrams: properties, algorithms and applications[END_REF], also called the Laguerre-Voronoi diagram, is another partition of a 2D plane into polygonal cells defined from a set of circles in R 2 . In the general case, for a set of circles (C i ) i∈K (or spheres) in R n with n ≥ 2, the power region R i associated with C i consists of all points whose power distance to C i is not greater than their power distance to any other circle C j , for j = i. The the power distance from a point x to a circle C i with center c i and radius r i is defined as

d P (x, C k ) = x -c k 2 -r 2
k , the formula of the Power region R i is then given by:

R i = {x ∈ R n | d P (x, C i ) ≤ d P (x, C j ) , ∀ j = i}.
The power diagram is a form of generalized Voronoi diagram, in the sense that you can take the circles C i instead of the centers c i and simply replace the distance function d V in the Voronoi diagram with the power distance function d P , to obtain the power diagram, see an example of three circles in R 2 in Figure 8. Notice that the power distance is not a real distance function. By summing up the volume of each power region inside the SAS-cavity, one can calculate the exact volume of the cSAS-cavity, which is equivalent to calculate the volume of the union of balls, see [START_REF] Avis | Computing the volume of the union of spheres[END_REF][START_REF] Cazals | Computing the volume of a union of balls: a certified algorithm[END_REF].

c.. New Voronoi-Type Diagram

In this paper, we propose a new Voronoi-type diagram for a set of spheres in R 3 (or circles in R 2 ), which is inspired by the non-overlapping decomposition of the SAScavity. We will see that this new diagram allows to calculate exactly not only the volume of the (complete or exterior) SAS-cavity, but also the (complete or exterior) SES-cavity which will be defined in the next section.

We first look at the new Voronoi-type diagram for a set of discs in the case of 2D. Notice that the boundary γ of the union of these discs can be classified into two types: open circular arcs {l 1 , l 2 , . . . , l n } and intersection points {x 1 , x 2 , . . . , x n }, with the number of circular arcs equal to the number of intersection points. Take A 1 = l 1 , . . . , A n = l n , A n+1 = {x 1 , x 2 , . . . , x n } in the Voronoi diagram. In consequence, we obtain n + 1 corresponding new Voronoi-type regions {R 1 , . . . , R n , R n+1 }, where R i is given by (4.18).

From a similar mapping R and similar equivalence theorems, we know that the part of R i inside γ is a circular sector when 1 ≤ i ≤ n, and R n+1 is the remaining region composed of polygons, see an example of three circles in Figure 8. For any point x ∈ R 2 , we have x ∈ R i if and only if there exists a point in A i such that it is a closet point to x on γ.

In the 3D case, the SAS consists of three types of geometrical quantities: open spherical patches, open circular arcs and intersection points. Similarly to the case in 2D, we take

A 1 = P 1 , . . . , A m 1 = P m 1 , A m 1 +1 = l 1 , . . . , A m 1 +m 2 = l m 2 and
A m 1 +m 2 +1 = {x 1 , . . . , x m 3 } with the above-mentioned notations. In consequence, we obtain

m 1 + m 2 + 1 corresponding new Voronoi-type regions {R 1 , . . . , R m 1 , R m 1 +1 , . . . , R m 1 +m 2 , R m 1 +m 2 +1 }.
With the mapping R defined in the last section, we know that the spherical sector R(P m ) corresponds to R m in the SAS-cavity, the double-cone region R(l m ) corresponds to R m 1 +m in the SAS-cavity, and the closed region R(I) coincides with R m 1 +m 2 +1 . For any point x ∈ R 3 , we have x ∈ R i if and only if there exists a point in A i such that it is a closet point on the SAS to x. The new Voronoi-type diagram is a powerful tool for calculating analytically the molecular area and volume, which is a direct consequence of the three equivalence statements. Given the components of the SAS, the new Voronoi-type diagram can be obtained directly according to the three equivalence statements, while the power diagram needs more complicated computations.

SAS-Area and SAS-Volume

In this paper, the area of the SAS is called the SAS-area and the volume of the SAS-cavity is called the SAS-volume. Similarly, the area of the SES is called the SES-area and the volume of the SES-cavity is called the SES-volume. Next, we will use the Gauss-Bonnet theorem of differential geometry to calculate the SAS-area, and then will use the new Voronoi-type diagram to calculate the SAS-volume.

a.. SAS-Area

To calculate analytically the solvent accessible area, the Gauss-Bonnet theorem has already been applied in the 1980s [START_REF] Connolly | Analytical molecular surface calculation[END_REF], which allows one to calculate the area of each SAS spherical patch with the information along its boundary, and then sum them up. We introduce briefly the key formula used in this calculation to keep the article complete. The area of a spherical patch P can be obtained from the Gauss-Bonnet theorem as following (see Figure 5):

i α i + i k e i |e i | + 1 r 2 A P = 2πχ, (4.19)
where α i is the angle at the vertex v i on the boundary between two neighboring circular arcs e i-1 and e i , k e i is the geodesic curvature of the edge e i , |e i | is the length of this edge e i and A P is the area of the patch P . Finally, χ is the Euler characteristic of P , which is equal to 2 minus the number of loops forming the boundary of P .

Through the above formula, we calculate the area of each SAS spherical patch A Pm and sum them up to get the SAS-area A sas :

A sas = m 1 m=1 A Pm . (4.20)

b.. SAS-Volume

In the new Voronoi-type diagram, the SAS-cavity is decomposed into three types of region: spherical sectors, double-cone regions and polyhedrons. For a spherical patch P m on S i with center c i and radius r i + r p , the volume of the corresponding spherical sector R(P m ), denoted by V Pm , can be calculated as following:

V Pm = 1 3 A Pm (r i + r p ).
For a circular arc l m associated with S i and S j (see Figure 6), the volume of the corresponding double-cone region R(l m ), denoted by V lm , can be calculated as following:

V lm = 1 6 r lm |l m | c i -c j .
where r lm is the radius of this circular arc and |l m | is the length of l m . Finally, the volume of the closed region T = R(I) can be calculated according the Gauss-Green theorem [START_REF] Federer | The gauss-green theorem[END_REF]:

V T = t A t (n t • n 0 ),
where t denotes a triangle on the boundary of T , A t is the area of the triangle t, n t is the outward pointing normal vector of t and n 0 is an arbitrary given unit direction vector in R 3 , for example, n 0 = (1, 0, 0). From the above three formulas, we sum up the volume of each spherical sector, each double-cone region and the polyhedron T , to get the solvent accessible volume V sas as below:

V sas = 1 3 m 1 m=1 A Pm (r i + r p ) + 1 6 m 2 m=1 r lm |l m | c i -c j + t A t (n t • n 0 ). (4.21)

The Solvent Excluded Surface

The Solvent Excluded Surface was first proposed by Lee & Richards in the 1970s [START_REF] Richards | Areas, volumes, packing and protein structure[END_REF]. Although the SES is believed to give a more accurate description of the molecular cavity, it is more complicated than the other two molecular surfaces. We first give two more precise definitions of the SES considering inner holes of a molecule. After that, we define mathematically different types of patches on the SES, which helps us to characterize and calculate analytically singularities on the SES. Finally, combining the above-proposed new Voronoi-type diagram and the upcoming singularity analysis, we calculate the exact SES-area and the SES-volume.

Mathematical Definitions

In previous works [START_REF] Lee | The interpretation of protein structures: estimation of static accessibility[END_REF][START_REF] Richards | Areas, volumes, packing and protein structure[END_REF], the SES is defined as the topological boundary of the union of all possible spherical probes that do not intersect any VdW atom of the molecule. In other words, the SES is the boundary of the cavity where a spherical probe can never be present. However, there might exist inner holes inside the solute molecule. Similar to the definitions of the cSAS and the eSAS, the complete Solvent Excluded Surface (cSES) is defined as the set of all points with signed distance -r p to the cSAS, while the exterior Solvent Excluded Surface (eSES) is defined as the set of all points with signed distance -r p to the eSAS. We make a similar convention as before that the SES refers to both the cSES and the eSES. From the geometrical relationship between the SAS and the SES, i.e. SES = f -1 sas (r p ), we propose an implicit function of the SES:

f ses (p) = -p -x p sas + r p if f sas (p) ≤ 0, p -x p sas + r p if f sas (p) ≥ 0, (5.22)
where x p sas is a closest point on the SAS to p, which depends on p and can be obtained directly from the equivalence statements.

In Connolly's work [START_REF] Connolly | Analytical molecular surface calculation[END_REF], the SES is divided into three types of patches: convex spherical patches, saddle-shaped toroidal patches and concave spherical triangles, see Figure 9. The convex spherical patches are the remainders of VdW-spheres, which occur when the probe is rolling over the surface of an atom and touches no other atom. The toroidal patches are formed when the probe is in contact with two atoms at the same time and rotates around the axis connecting the centers of these two atoms. While rolling, the probe traces out small circular arcs on each of the two VdW-spheres, which build the boundaries between the convex spherical patches and the toroidal patches. The concave spherical triangles occur if the probe is simultaneously in contact with more than or equal to three VdW-spheres. Here, the probe is in a fixed position, meaning that it is centered at an SAS intersection point and cannot roll without losing contact to at least one of the atoms.

The intersection between different SES patches might occur which leads to singularities on the SES, referred to as SES-singularities, see Figure 10. Despite some particular cases which have been studied by Sanner [START_REF] Sanner | Reduced surface: an efficient way to compute molecular surfaces[END_REF], a characterization of these singularities is not known. We will provide a complete characterization, using the above equivalence statements as well as the following analysis of the SES-singularities. To start with, we give three new definitions of different SES patches from the mathematical point of view: 1) Convex Spherical Patch: A convex spherical patch on the SES, denoted by P + , is defined as the set of the points on the SES such that there exists According to these new definitions, the three types of patches can be rewritten mathematically as follows:

       P + = {p : f ses (p) = 0, p ∈ R(P m )}, 1 ≤ m ≤ m 1 P t = {p : f ses (p) = 0, p ∈ R(l m )}, 1 ≤ m ≤ m 2 P -= {p : f ses (p) = 0, p ∈ R(x m )}, 1 ≤ m ≤ m 3 (5.23)
where f ses is the implicit function of the SES and R is the mapping defined in the previous section. Actually, a convex spherical patch and a toroidal patch are defined in the same way as Connolly, while a new-defined concave patch might not be triangle-shaped because its new definition takes into account the intersection among different SES patches. The above new definitions of the SES patches ensure that different SES patches will not intersect with each other, for the reason that different patches belong to different new Voronoi-type regions in the new Voronoitype diagram.

SES-Singularities

Before characterizing the singularities on the SES, it is necessary to recall the properties of the signed distance function f S to a surface S in R n (n ≥ 1) as below:

1) f S is differentiable almost everywhere, and it satisfies the Eikonal Equation:

|∇f S | = 1.

2) If f S is differentiable at a point p ∈ R n , then there exists a small neighborhood V of p such that f S is differentiable in V .

3) For any point p ∈ R n , f S is non-differentiable at p if and only if the number of the closest points on S to p is greater than or equal to 2.

We call a point x ses ∈ SES a singularity if the SES is not smooth at x ses , which means that its implicit function f ses (p) = f sas (p) + r p is non-differentiable at x ses in R 3 . If x ses ∈ SES is a singularity, we obtain consequently that f sas is nondifferentiable at x ses . From the last property above, we therefore can characterize the SES-singularities by the following equivalence.

Corollary 5.1. A point x ses ∈ SES is a singularity if and only if the number of the closest points on the SAS to x ses is greater than or equal to 2. We now investigate the different types of singularities that can appear on each of the three patch types. From the definitions of different SES patches, it is relatively easy to calculate the convex spherical patches and the toroidal patches given the components of the SAS. In the following, we illustrate how to calculate P -exactly, which provides therefore a complete characterization of the SES-singularities later. Denote by P 0 the concave spherical triangle (P -⊂ P 0 ) corresponding to an SAS intersection point x m showed in Figure 11 and by K the set of all intersection points in I with distances less than 2r p to x m , that is, K = {x ∈ I : x -x m < 2r p , x = x m }, so that K collects all SAS intersection points "near" x m . According to the definition of P -corresponding to x m and the equation (4.17), we have the following formula:

P -= P 0 ∩ R(x m ) = P 0 ∩ T ∩ {p : p -x m ≤ p -x , ∀x = x m , x ∈ I}. (5.24)
The above formula can be used to calculate P -directly, in which R(x m ) is a polyhedron. However, this formula is not convenient to calculate, since one has to calculate T as well as the union of P 0 and R(x m ). This motivates us to look deep into the relationship between P -and P 0 . In the following theorem, we will present a simpler formula of P -, which allows us to calculate it analytically and more efficiently. The proof of Theorem 5.1 is given in the Appendix. From the definition of P -, it is not surprising to have the inclusion P -⊂ P 0 \ x∈K B rp (x). However, the above theorem gives a stronger result that P -and P 0 \ x∈K B rp (x) are identical. This means that the concave spherical patch can be obtained from P 0 by removing the parts intersecting other "nearby" spherical probes centered at SAS intersection points. This theorem also allows us to characterize the singular circular arcs on the concave spherical patch P -. We propose a theorem of the SES-singularities as following.

Theorem 5.2. The following statements hold [START_REF] Aurenhammer | Power diagrams: properties, algorithms and applications[END_REF] There can not exist any singularity on a convex spherical SES patch P + .

[2] On a toroidal patch P t , two point-singularities occur when the corresponding SAS circular arc has a radius smaller than r p and they can be computed following the sketch in Figure 13.

[3] On a concave spherical SES patch P -, singular arcs occur when P -does not coincide with the corresponding concave spherical triangle P 0 . Further, these singular arcs form the boundary of P -that does not belong to the boundary of P 0 .

Proof. First, if x ses is a point on a convex spherical SES patch, then it has a closest point x sas on an spherical SAS patch. Moreover, x sas is the unique closest point to x ses from Theorem 4.1. This, in turn, implies by Corollary 5.1 that x ses can not be a singularity. In consequence, there can not exist any singularity on a convex spherical SES patch P + . Second, if x ses is a point on a toroidal SES patch, then it has a closest point x sas on an SAS circular arc l m associated with two SAS-spheres S i and S j . Moreover, x sas is not unique if and only if x ses belongs to [c i , c j ] by Theorem 4.2, which happens only when the radius of l m is smaller than r p and x ses is one of two cusps on the toroidal SES patch as showed in Figure 13. By Corollary 5.1, two point-singularities on the toroidal patch P t can only occur when the corresponding SAS circular arc has a radius smaller than r p .

Third, if x ses is a point on a concave spherical SES patch P -corresponding to an SAS intersection point x m , then x sas = x m is a closest point to x ses . If x ses belongs to the boundary of P -but not to the boundary of P 0 , then by the formula characterizing P -in Theorem 5.1 we know that x ses lies on another nearby probe ∂B rp (x) where x ∈ K (and thus x = x m ). In consequence, x is another closest point to x ses implying that x ses is singularity by Corollary 4.1. On the other hand, if x ses ∈ P -but does not belong to x∈K B rp (x), then x ses does not lie on any nearby probe and x sas is the unique closest point among all SAS intersection points. Assume by contradiction that there exist another closest point x to x ses on some SAS spherical patch P m . Since now by Theorem 4.1 the closest point to any point on (c i , x] on the SAS is unique, this implies that x ses = ci (because x m is another closest point). This is however a contradiction since x ses can not coincide with the center of any SAS-sphere as all the VdW-radii are assumed to be positive. Further, assume by contradiction that there exist a closest point x to x ses on some SAS circular arc l m . Consequently, x ses lies on the corresponding toroidal patch. According to Theorem 4.2, x ses belongs to [c i , c j ]. If r p < r lm , then x ses can not belong to the SES which is a contradiction. If r p ≥ r lm , then x ses has to be one of the two cusps. In this case, the two ending points of l m are both closest points to x ses which are also SAS intersection points. This conflicts with the fact that x sas = x m is the unique closest point among all SAS intersection points.

Remark 5.1. In the third case of Theorem 5.2, P -can be calculated according to Theorem 5.1. It is obvious that P -does not coincide with P 0 if and only if

P 0 Ä x∈K B rp (x) ä is nonempty.
Finally, we state a corollary about classifying all points on the SES into four classes according to the number of its closest points on the SAS. Corollary 5.2. For any point x ses ∈ SES, assume that the number of its closest points on the SAS is N (denoting N = ∞ for an infinite number of closest points). Then, there exists four cases:

[1] N = 1: x ses is not a singularity on the SES, and x ses has an unique closest point {x 1 sas } on the SAS.

[2] 2 ≤ N < ∞: x ses is a singularity on a concave spherical SES patch, and its closest points {x 1 sas , . . . , x N sas } are among the SAS intersection points.

[3] N = ∞: x ses is a singularity on a toroidal patch corresponding to an SAS circular arc (or a complete circle) on which any point is a closest point on the SAS.

SES-Area and SES-Volume

With the new Voronoi-type diagram as well as the above singularity analysis we can calculate the SES-area and the SES-volume analytically, which shall be explained in the following.

a.. SES-Area

In Connolly's paper [START_REF] Connolly | Analytical molecular surface calculation[END_REF], the presence of singularities on the SES concave spherical patches made it infeasible to calculate the SES-area exactly. The characterization of the singularities carried out earlier in this paper allows us to calculate the area of each SES patch and sum them up to obtain the exact area of the whole SES. To calculate the area of a convex spherical patch P + , we use a similar Gauss-Bonnet formula as (4.19):

v α v + e k e |e| + 1 r 2 i A P + = 2πχ, (5.25)
where α v denotes the angle at a vertex v between neighboring circular arcs on the boundary of P + , k e the geodesic curvature of an edge e on the boundary of P + , A P + is the area of P + and χ is the Euler characteristic of P + .

To calculate the area of a toroidal patch P t analytically, we consider two cases, and suppose that P t corresponds to an SAS circular arc l m with radius r lm . In the case where r lm > r p , there will be no singularity on P t . With the notations introduced in Figure 12, we therefore have the following formula for calculating the area of P t :

A Pt = r p β m [r lm (θ 1 + θ 2 ) -r p (sin θ 1 + sin θ 2 )] .
(5.26)

In the case where r lm ≤ r p , there are two singular points on P t . With the notations introduced in Figure 13, we have the following formula for calculating the area of P t :

A Pt = r p β m [r lm (θ 1 + θ 2 -2θ 0 ) -r p (sin θ 1 + sin θ 2 -2 sin θ 0 )] . (5.27)
To calculate analytically the area of a concave spherical patch P -, we use the Gauss-Bonnet theorem on the spherical probe again, see Figure 14. In the previous section, we have obtained that P -=P 0 \ x∈K B rp (x), which implies that all information about the boundary of P -is known. It allows us to apply the Gauss-Bonnet theorem on the spherical probe, to obtain:

v α v + e k e |e| + 1 r 2 p A P -= 2πχ, (5.28)
where α v denotes the angle at a vertex v on the boundary of P -, k e the geodesic curvature of an edge e on the boundary of P -, A P + is the area of P -and χ is the Euler characteristic of P -. In summary, the area of each SES patch can be calculated independently and analytically, and we can then sum them up to obtain the exact SES-area A ses .

Figure 12: The yellow patch is a toroidal patch P t on the SES corresponding to an SAS circular arc l m with the radius r lm > r p and the radian β m . θ 1 denotes the angle between the line connecting c i with a point on l m and the disc where l m lies. Similarly, θ 2 denotes the angle between the line connecting c j with the same point on l m and the disc where l m lies.

Figure 13: The two yellow parts form a toroidal patch P t on the SES corresponding to an SAS circular arc l m with the radius r lm < r p and the radian β m . θ 0 denotes the angle between the line connecting a singularity on P t with a point on l m and the disc where l m lies and θ 1 , θ 2 are defined as in Fig. 12.

b.. SES-Volume

According to the new Voronoi-type diagram, we can calculate the exact SESvolume. We propose to subtract the volume of the region between the SAS and the SES from the SAS-volume to obtain the SES-volume. This region that needs to be subtracted and which is denoted by R s can be characterized as below

R s = {p : -r p ≤ f sas (p) ≤ 0}.
From the new Voronoi-type diagram, we decompose R s into small regions, each of which corresponds to a spherical patch P m , a circular arc l m or an intersection point

x m : R s ∩ R(P m ), R s ∩ R(l m ) and R s ∩ R(x m ). The volume of the region R s ∩ R(P m ) is given by V Rs∩R(Pm) = 1 3 A Pm (r i + r p ) Ç 1 - r 3 i (r i + r p ) 3 å , (5.29)
where A Pm is the area of the SAS spherical patch P m , r i + r p is the radius of the corresponding SAS-sphere S i on which P m lies. For R s ∩ R(l m ), denote by r lm the radius of l m and by β m the radian of l m . In the case where r lm > r p , using the notations of Figure 12, the volume of the region R s R(l m ) is given by

V Rs∩R(lm) = β m r 2 p ï r lm 2 (θ 1 + θ 2 ) - r p 3 (sin θ 1 + sin θ 2 ) ò . (5.30)
In the case where r lm ≤ r p , using the notations of Figure 13, it is given by

V Rs∩R(lm) = β m r 2 p ï r lm 2 (θ 1 + θ 2 -2θ 0 ) - r p 3 (sin θ 1 + sin θ 2 -2 sin θ 0 ) ò + 1 3 β m r 2 lm r 2 p -r 2 lm .
(5.31) Consider now R s ∩ R(x m ) corresponding to a concave spherical patch P -. Notice that there might be some flat regions on the boundary of R s ∩ R(x m ), caused by the intersection of the probe B rp (x m ) and its nearby probes, see Figure 14 (right). Denote by D i the i-th flat region with the boundary composed of line segments and circular arcs. Furthermore, denote by d i the distance from x m to the plane where D i lies, and by A D i the area of D i . Then, the volume of the region R s ∩ R(x m ) can be formulated as:

V Rs∩R(xm) = 1 3 A P -r p + i 1 3 A D i d i , (5.32)
where, A P -is the area of the concave spherical SES patch P -. Finally, by summing up the volume of each subtracted region, we obtain the subtracted volume as following:

V Rs = ξ=Pm, lm, xm
V Rs∩R(ξ) .

(5.33) Therefore, the SES-volume denoted by V ses is equal to V sas -V Rs .

Figure 14: On the left, the concave spherical triangle P 0 with vertices (a 1 , a 2 , a 3 ) corresponds to an intersection point x m . In the case where P -coincide with P 0 , there will be no singularity on the concave spherical patch. On the right, the concave spherical patch P -does not coincide with P 0 and there are singular circular arcs as parts of its boundary. The vertices of P -are (a 1 , a 2 , a 3 , a 4 , b 1 , b 2 , a 5 ). The two flat grey regions D 1 and D 2 are formed by the intersection of P 0 with two other nearby spherical probes. These two regions have the boundaries composed of line segments and circular arcs.

Construction of Molecular Surfaces

In the above sections, we have defined and analyzed the cSAS and the eSAS, as well as the cSES and the eSES. However, all work is based on the explicit knowledge of the components of these molecular surfaces. In this section, we will present a method of constructing the cSAS and the eSAS using a binary tree to construct its spherical patches. After that, we will give a brief construction strategy of the cSES and eSES based on the construction of the cSAS and the eSAS. The construction of molecular surfaces ensures that our previous analysis is feasible.

Construction of the cSAS and the eSAS

To start the construction, we need some quantities for representing different components of the SAS. First, an intersection point will be represented by its coordinate in R 3 and an identifier. Second, to represent a circular arc, we use its starting and ending intersection points (resp. the identifiers), its radius, center and radian as well as an identifier. Finally, to represent a spherical patch, we use the SAS-sphere on which it lies on and the identifiers of all circular arcs forming its boundary. Then, we propose to construct the SAS in five basic steps as below: 1) Compute the set of all intersection points {x 1 , • • • , x m 1 }, denoted by I.

2) Calculate each SAS circular arc or circle l m associated with some S i and S j . For each SAS circular arc or circle, we record the information, including the center, the radius, the radian, the corresponding two SAS spheres, the identifiers of the starting and ending intersection points. Notice that each circular arc connects two intersection points. 3) Construct all loops on each SAS sphere S i , which also form the boundaries of SAS spherical patches on S i . Notice that each loop is composed of circular arcs, or a complete circle. We start from a circular arc on S i , find another arc connecting this arc, add it into the loop, and repeat this until we finally obtain a complete loop. The k-th loop on the i-th SAS sphere is denoted by L i k . 4) Construct all spherical patches on each SAS sphere S i . Since the boundary of a spherical patch on S i is composed by one or several loops on S i , we can use the identifiers of these loops to represent a spherical patch. Denote by P i k the k-th spherical patch on S i . The difficulty lies on determining whether two loops on S i belong to the boundary of a common spherical patch or not. This final question will be discussed in the next section. 5) Finally, we should distinguish the cSAS and the eSAS. The cSAS is just the set of all above-constructed SAS spherical patches. To construct the eSAS, we say that two spherical patches P i k and P j l are neighbors if they have a common circular arc or circle on their boundaries. Then, we start our construction of the eSAS, by mapping a faraway point p ∞ onto an SAS spherical patch, which is the initial patch on the eSAS. Then, we add the neighboring spherical patches into the eSAS one by one, to finally obtain the whole eSAS.

Binary Tree to Construct Spherical Patches

In the above construction process, there remains a problem of classifying the loops on an SAS-sphere into several parts, such that each part forms the boundary of a spherical patch. To do this, we need to determine whether two loops on the SAS-sphere belong to the boundary of a same spherical patch or not. Note that two different loops won't cross each other but can have common vertices. We propose to construct a binary tree whose leaves are the different spherical patches.

Given a loop L on a fixed SAS-sphere S i , L divides the sphere into two open parts and it is composed of circular arcs formed by the intersection of S i and other SAS-spheres. We call the open part of S i which is not hidden by those other SASspheres as the interior of L, denoted by L • , while the other open part as the exterior of L, denoted by L c . We say that another loop L on S i is inside L if L ∈ L • , where L • = L ∪ L • is the closed hull of L • on S i . Notice that each loop L being part of the boundary of a spherical patch contains all the other loops that belong to the boundary of the same patch. By testing if a loop L on S i is inside the loop L, we can classify all loops on S i into two parts: the loops inside L (including L itself) and the remaining loops outside L. We do this division repeatedly until each loop is tested which results in a binary tree.

To better understand this process, let's see the example of Figure 15. On an SASsphere S i , there are six loops {L 1 , L 2 , L 3 , L 5 , L 6 , L 6 } forming three spherical patches {P 1 , P 2 , P 3 } with green, red and blue boundaries. The leaves of the binary tree in Figure 15 represent the boundary of the three spherical patches. To construct the tree, consider all loops {L 1 , L 2 , L 3 , L 5 , L 6 , L 6 } and test with L 1 , to divide into the set of loops inside L 1 ({L 1 , L 2 , L 3 , L 6 }) and the set of the remaining loops outside L ({L 4 , L 5 }). Then, we test with L 2 , to find that L 2 is inside itself, while {L 1 , L 3 , L 6 } are outside, which implies that L 2 itself forms the boundary of a spherical patch. Afterwards, we find that {L 1 , L 6 } are both inside L 3 by testing with L 3 and that {L 1 , L 3 } are inside L 6 by testing with L 6 . This implies that {L 1 , L 3 , L 6 } form the boundary of the second spherical patch, as these three loops are inside each other. Finally, we test respectively with L 4 and L 5 to find that they are inside each other and they form the boundary of the last spherical patch.

Interior of a Loop

Let L and L be two loops on S i . It is left to explain how to test whether L is inside L or not. We assume that L is composed by n circular arcs which are formed by the intersection of S i and other SAS-spheres S 1 , . . . , S n . Denote all intersection circles by C 1 , . . . , C n , where C j = S i ∩ S j . The corresponding SAS-balls to S 1 , . . . , S n are denoted by B 1 , . . . , B n . Since L and L do not cross each other, we can determine if the loop L is inside L or not by testing whether a particular point x ∈ L is inside L or not. We denote by x k (1 ≤ k ≤ n) the closest point on C k to x, which can be analytically given. Notice that x k has the smallest Euclidean distance in R 3 from x to any point on C k , which is equivalent to that the shortest path on S i from x to C k which ends at x k . For all circles C k , k = 1, . . . , n, we can then find the circle C k 0 with the minimum path length on S i from x to C k . The closest point on C k 0 to x is thus denoted by x k 0 . With these notations, the following lemma is proposed to test whether x is inside L or not, note that L ∩ n j=1 B j = ∅.

Lemma 6.1. Given an arbitrary point x ∈ S i \ n j=1 B j , x is inside L if and only if

x k 0 ∈ L.
Proof. We first prove the sufficiency. Assume that x is inside L, i.e. x ∈ L • = L • ∪ L. Assume by contradiction that x k 0 ∈ C k 0 does not belong to L. Then, since C k 0 ⊂ L ∪ L c it follows that x k 0 ∈ L c . Then, the shortest path on S i from x ∈ L • to x k 0 ∈ L c must cross L. In consequence, the intersection point between L and the shortest path on S i has a smaller path length to x than x k 0 does, which is a contradiction. Therefore there holds that x k 0 ∈ L.

Second, we prove the necessity. Assume that x k 0 ∈ L. By contradiction again, we assume that x is not inside L, i.e. x ∈ L c , which yields that x ∈ Ω := L c \ n j=1 B j from the lemma's condition. Then, the shortest path on S i between x ∈ Ω and x k 0 / ∈ Ω must cross Γ = ∂Ω. In consequence, the intersection point between Γ and the shortest path on S i has a smaller path length to x than x k 0 does, which is a contradiction. Therefore, x lies inside L.

Construction of the cSES and the eSES

After the construction of the cSAS and the eSAS, patches of the cSES and the eSES can consequently be distinguished since each patch corresponds to one component of the SAS. An SES patch belongs to the cSES if and only if its corresponding SAS patch belongs to the cSAS. In analogy, an SES patch belongs to the eSES if and only if its corresponding SAS patch belongs to the eSAS.

To construct the whole SES, one can construct each SES patch one by one. A convex spherical patch P + on the SES can be obtained directly by shrinking its corresponding SAS-patch P m from the SAS-sphere to the VdW-sphere. A toroidal patch P t corresponding to an SAS circular arc l m can also be computed without difficulty given the radius, center, starting and ending points of l m , considering two cases Table 1: Surface areas and cavity volumes for a two intersecting atomic system obtained respectively from the closed form expression and the proposed method in the paper. The last row records the difference between these two results.

respectively in Figure 12 and Figure 13 (see also [START_REF] Connolly | Analytical molecular surface calculation[END_REF]). To calculate a concave spherical patch corresponding to an SAS intersection point x m , we can use the formula presented in Theorem 5.1:

P -= P 0 \ x∈K B rp (x),
which is one of the main results of this paper. Notice that P -lies on the spherical probe centered at x m and its boundary is constituted by circular arcs, which can be computed analytically. To construct P -, we can consequently use the same treestructure (see also Figure 15 on the right) as was used in the previous subsection for constructing SAS spherical patches to classify loops into different sets on the spherical probe (instead of on an SAS-sphere). Some illustrations of the construction of the SESs are provided in the next section (see for instance Figure 20 -21).

Numerical results

In this section, we illustrate some numerical results with the proposed method based on the previous theoretical results.

A two atomic system

The proposed method gives the analytic expression of surface areas and cavity volumes for both the SAS and the SES if round-off errors are neglected. Here, we test the simplest case of a two atomic system for the comparison of molecular areas and volumes between the closed form expressions based on the properties of a spherical cap (see website [START_REF] Wikipedia | Spherical cap -wikipedia, the free encyclopedia[END_REF]) and the proposed method that relies on the Gauss-Bonnet theorem.

With Matlab, we compare the closed form expressions of the surface areas and cavity volumes with the proposed method. We considered two hydrogen atoms (having radius 1.2 Å) with distance 2 Å and the probe radius is also set to 1.2 Å. Table 1 reports the molecular areas and molecular volumes obtained. We confirm that the result coincides up to a small error due to round-off errors. ∆V ses = V eses -V cses and ∆V sas = V esas -V csas . The probe radius r p is fixed to be 1.5 Å.

Molecular areas and volumes

The complete and exterior surface of a molecule can have different areas and volumes. Table 2 reports the surface areas and cavity volumes of different molecules and also the relative errors between the complete results and the exterior results. We observe that the relative difference between the complete and exterior SES-area and SES-volume is up to 12.5% resp. 2.7% in this test cases.

Comparison with the MSMS-algorithm

We have run the MSMS-algorithm to compare it with the proposed method. It has been observed that if there is no SES-singularity, the exterior molecular areas and volumes from the MSMS coincides with the results of the exterior molecular surface from the proposed method in general. However, in the case where SESsingularities occur on concave spherical patches (which is often the case for complex molecules), the MSMS might fail to give the right SES-area (or SES-volume) due to an incomplete characterization of all SES-singularities.

This can be illustrated on a simple three atomic system where we calculate the molecular areas respectively obtained from the MSMS and the proposed method. Consider that there are three atoms with the same radius 1.2 Å. In the case of a probe radius of r p = 2 Å (see Figure 16), there is no SES-singularity and the molecular areas obtained from the two methods coincide (see Table 3). In the case of a probe radius of r p = 1 Å (see the left of Figure 17), there is a singular circle on the concave spherical patches. The SAS-areas obtained from the two methods still coincide, but The numerical SES-areas reported in Å2 for the three atomic system (r p = 1 Å, see the left of Figure 17) with respect to the number of triangles in the mesh. For each method, the analytical and numerical area is computed.

the SES-areas are different (see Table 3) because the MSMS has ignored the "selfintersection" of the surface by simply treating the concave spherical patches as two triangle-shaped spherical patches (see the right of Figure 17). Furthermore, we have tested the proposed method for sequences of an increasing number of triangles in the mesh for the case of three atoms with r p = 1 Å, see Figure 18. Indeed, based on the provided analysis one can construct meshes at various sizes, see also Section 7.5. Like this, we obtain the analytical SES-area as well as a numerical (approximate) value derived from the mesh. We observe that this numerical value converges to the analytical one if using finer and finer meshes. In the same manner, the MSMS-algorithm provides an analytical and a numerical value for the SES-area. However, the "analytical" value does not account for "selfintersection" and therefore this value is too large. In consequence, the numerical value also converges to this wrong limit.

In addition, we have computed the exterior SAS-areas and the exterior SESareas of molecules at various sizes by both the MSMS and the proposed method for comparison, see Table 4.

Table 4: Exterior molecular areas ( Å2 ) of molecules at different sizes respectively obtained from the MSMS and the proposed method (r p = 1.5 Å), where A * esas and A * eses are results obtained from the MSMS, A esas and A eses are results obtained from the proposed method, ∆A esas = A * esas -A esas and ∆A eses = A * eses -A eses .

Time cost

We have implemented the proposed method of computing the areas and volumes with Matlab for molecules at various sizes. Figure 19 demonstrates the relationship between the total CPU-time and the size of the molecule, where we observe almost a linear relationship. The program was run on a laptop with 2.5GHz quad-core Intel Core i7 processor. This first implementation in Matlab is to be understood as a proof of concept. The important information to retain is that it can be implemented in a linearly scaling fashion. Of course, in order to lower the pre-constant of the linear scaling, one has to refer to a proper and more professional implementation in a better performing language. We therefore expect a better performance of the proposed method in fortran or C++.

Molecular visualization

Molecular visualization is a very important application of all the analysis about the surfaces that we have derived. One of the possibilities is to visualize meshes of the surfaces. In this case, the basic idea is to mesh different molecular patches one by one, given the boundary information for each of them. In this paper, we visualize the meshes of two SESs obtained from the basic advancing-front meshing algorithm [START_REF] George | The advancing-front mesh generation method revisited[END_REF] for each patch as a first implementation, see Figure 20

Conclusion

In this paper, based on three equivalence statements and a new Voronoi-type diagram, an analytical and also computable implicit function of the SES has been provided. Furthermore, all SES patches including SES-singularities can be characterized and computed according to Theorem 5.1 and Theorem 5.2. Consequently, different molecular areas and volumes can be calculated analytically with explicit formulas. In addition, we have also proposed the concept of the cSAS and the eSAS (as well as the cSES and the eSES), considering the possible inner holes inside a solute molecule or not. Finally, a corresponding method to construct (or distinguish) the cSAS and the eSAS has been given, with the use of a binary tree. In summary, the molecular surfaces have been analytically characterized (especially the SES), which allows to do more accurate calculation in solvation models associated with different molecular cavities than the state-of-the-art methods. With the above lemma, we can now prove Theorem 5.1. According to the definition of K, it accounts for all SAS intersection points that might intersect B rp (x m ) and thus P 0 \ x∈K B rp (x) = P 0 \ x∈I B rp (x). We first observe that (A.2)

In the above, the first equality results from Lemma 1.1, the third one from equation (A.1) and the last one from (5.24).

Figure 1 :

 1 Figure 1: This is a 2-dimension (2D) schematics of the Solvent Accessible Surface and the Solvent Excluded Surface, both defined by a spherical probe in orange rolling over the molecule atoms in dark blue.

Figure 2 :

 2 Figure 2: The above figures illustrate two SESs of two artificial molecules respectively containing three atoms. In each of them, f 1,2 denotes the signed distance function to the SES of the 1-th and the 2-th atoms which is depicted with dashed red curves (f 2,3 and f 1,3 are similar).

Figure 3 :

 3 Figure 3: The left figure shows the components of the eSAS of the protein 1B17 (the probe radius r p = 1.2 Å): open spherical patches in blue, open circular arcs in yellow and intersection points in red. The right figure shows the exterior surface (transparent) and the interior surfaces of the eSES of the protein 1B17. The boundary of an exterior spherical patch of the eSAS is composed of circular arcs depicted in yellow, and the boundary of an interior spherical patch of the cSAS is composed of circular arcs depicted in purple.

Proposition 4 . 1 .

 41 Consider two bounded, closed and oriented surfaces S ⊂ R 3 and S ⊂ R 3 with two corresponding signed distance functions f S (p) and f S (p). If the cavity inside S is contained in the cavity inside S , then we have f S

Figure 5 :

 5 Figure 5: Illustration when the point x sas lies on an SAS spherical patch P m which is part of an SAS-sphere S i , and p lies on the segment [c i , x sas ]. V represents a neighborhood of x sas on the spherical patch. α s is the angle variation at the vertex v s between two neighboring circular arcs e s-1 and e s on the boundary of this spherical patch.

Figure 6 :

 6 Figure 6: Illustration when l m is an SAS circular arc associated with two SAS-spheres S i and S j .The point x sas on l m is a closest point on the SAS to p. γ 0 is a small neighborhood of x sas on the circular arc l m , whereas γ 1 and γ 2 represent two small neighboring curves of x sas on the SAS, with x sas as the endpoints, γ 1 ⊂ Π ∩ S i ∩ SAS, and γ 2 ⊂ Π ∩ S j ∩ SAS.

Corollary 4 . 1 .

 41 If x sas is an SAS intersection point x m ∈ I associated with S i , S j and S k , then x sas is a closest point on the SAS to p if and only if p lies in the closed region T and x sas is a closest point in I to p.

Figure 7 :

 7 Figure 7: This figure illustrates the closed region T of three SAS-spheres with the centers (c i , c j , c k ), where the tetrahedron T has five vertices (x sas , c i , c j , c k , x I ). Here, x sas is one SAS intersection point in I, and x I is the other SAS intersection point.

Figure 8 :

 8 Figure 8: The left figure gives the Voronoi diagram of three points {a 1 , a 2 , a 3 } in R 2 . R 1 is the Voronoi region associated with a 1 , while R 2 associated with a 2 and R 3 associated with a 3 . The middle figure gives the power diagram of three circles {C 1 , C 2 , C 3 } in R 2 , with three corresponding power regions R 1 , R 2 and R 3 . The right figure gives the new Voronoi-type diagram of three circles in R 2 , with four regions R 1 , R 2 , R 3 and R 4 respectively corresponding to l 1 , l 2 , l 3 and {x 1 , x 2 , x 3 }.

Figure 9 :

 9 Figure 9: The above figures show both the SAS and the SES of the caffeine molecule (the probe radius r p = 1.2 Å). On the left, the SAS is composed of spherical patches in blue, circular arcs in yellow and intersection points in red. On the right, the patches in red (resp. in yellow or in blue) are the corresponding convex spherical patches (resp. toroidal patches or concave spherical patches) on the SES.

Figure 10 :

 10 Figure10: The above figures show two kinds of singularities respectively on a toroidal patch P t and on a concave spherical patch P -. On the left, each of the two point-singularities (cusps) on the toroidal patch P t has an infinite number of closest points on the corresponding SAS circle in orange. On the right, the SES-singularities on the concave spherical patch P -form a singular circle. The two spherical spheres in green center at two SAS intersection points in red.

Figure 11 :

 11 Figure11: This is the schematic of the concave spherical triangle P 0 corresponding to an intersection point x m , with the boundary composed of three circular arcs in blue ( ā1 a 2 , ā2 a 3 and ā3 a 1 ) on the spherical probe. This spherical triangle touches at the same time three VdW-balls at vertices (a 1 , a 2 , a 3 ).

Theorem 5 . 1 .

 51 P -=P 0 \ x∈K B rp (x) where B rp (x) denotes the open ball (or disc in R 2 ) centered at x with radius r p .

Figure 15 :

 15 Figure15: On the left is is a brief schematic of an SAS-sphere and the loops on it. There are six loops on the SAS-sphere {L 1 , L 2 , L 3 , L 5 , L 6 , L 6 }, and three spherical patches with the boundaries formed by two loops in green {L 4 , L 5 }, three loops in red {L 1 , L 3 , L 6 } and one loop (circle) in blue {L 2 } respectively. The tree on the right illustrates the corresponding binary tree whose leaves identify the boundaries of three different spherical patches.

Figure 16 :

 16 Figure 16: The SES of a three atomic system without SES-singularity (r p = 2 Å).

Figure 17 :

 17 Figure 17: The SES of a three atomic system with a singular circle on the concave spherical patches (left) and the SES with the intersection of two concave spherical patches (right, triangulation provided by the MSMS-algorithm). The probe radius r p = 1 Å.

Figure 18 :

 18 Figure 18: The numerical SES-areas reported in Å2 for the three atomic system (r p = 1 Å, see the left of Figure17) with respect to the number of triangles in the mesh. For each method, the analytical and numerical area is computed.

- 21 .

 21 In future work, we expect to combine the proposed method of generating molecule structures with different meshing algorithms for visualizing different molecular surfaces.

Figure 19 :

 19 Figure 19: Total CPU-times of computing molecular areas and volumes for molecules at different sizes, from tens to tens of thousands (r p = 1 Å).

Figure 20 :

 20 Figure 20: Two eSESs of a three atomic system (left, the MSMS fails to provide a triangulation in this case) with r p = 1.8 Å and the protein 1B17 (right) with r p = 1.5 Å. The green circular arcs constitute the boundaries of all concave spherical patches with singular arcs.

Figure 21 :

 21 Figure 21: Two eSESs of the molecule 101M (left) with r p = 1 Å and the molecule 3WPE (right) with r p = 1 Å. The green circular arcs constitute the boundaries of all concave spherical patches with singular arcs.

Figure A. 22

 22 Figure A.22: A 2D schematics of the two different polygons in the eSAS-cavity of an artificial molecule for the proof of Lemma 1.1. T xm is the triangle with three vertices (x m , c i , c j ) and the polygon in green is the reduced polygon T 0 .

Figure A. 23 :

 23 Figure A.23: 2D schematics of two cases where P 0 ∩ T x m = ∅ with x m being an SAS intersection point different from x m . The green segment [c i , c j ] is an edge of T 0 and a and b are the two intersection points between the circle ∂B rp (x m ) and the segment [c i , c j ]. The grey region is the region where P 0 can possibly appear outside T 0 that may be included in the triangle x m c i c j (left) or not (right).

  c) ∃ x m = x m s.t. P 0 ∩ T x m = ∅: Figure A.23 gives the schematics of the region (in grey) where P 0 can possibly appear outside T 0 . If this region is included in T x m , the proof is again trivial. If parts of this region lies outside T x m , then this region is include by the probe B rp (x m ) by symmetry. In consequence, all such regions can be covered by the corresponding spherical probes, implying that P 0 \T ⊂ x∈K B rp (x).

P 0 ÄP 0

 00 \B rp (x) = P 0 ∩ Ä B rp (x) ä c = P 0 ∩ {p : p -x m ≤ p -x }, ∀x = x m , (A.1)and therefore we obtain thatP 0 \ x∈K B rp (x) = P 0 \ x∈K B rp (x) T = x∈I \B rp (x) ä T = P 0 ∩ T ∩ {p : p -x m ≤ p -x , ∀x = x m , x ∈ I} = P -.

  1A0C from the protein data bank http://www.rcsb.org/pdb/home/ home.do) where we use the notations ∆A ses = A eses -A cses and ∆A sas = A esas -A csas resp.

	Name	Acses	Aeses	∆Ases Aeses	Acsas	Aesas	∆Asas Aesas
	1ETN	956.636	956.636	0	1397.664	1397.664	0
	1B17	3372.577	3372.577	0	4297.098	4297.098	0
	101M	7335.438	7096.387	-0.034	8424.749	8413.275	-0.001
	3WPE	28575.829	27390.001	-0.043	29747.737	29644.146	-0.003
	1A0C	61006.132	54231.268	-0.125	53107.185	51933.798	-0.023
	Name	Vcses	Veses	∆Vses Veses	Vcsas	Vesas	∆Vsas Vesas
	1ETN	1567.728	1567.728	0	3333.077	3333.077	0
	1B17	8630.389	8630.389	0	14435.224	14435.224	0
	101M	25029.270	25195.072	0.007	36942.276	36943.700	0.000
	3WPE	98651.803	99642.385	0.001	143181.258	143201.169	0.000
	1A0C	238282.734	244765.065	0.026	326348.079	326601.929	0.001
	Table 2: Different molecular areas (in	Å2 , top) and volumes (in	Å3 , bottom) for five molecules
	(1ETN, 1B17, 101M, 3WPE,

Table 3 :

 3 Molecular areas ( Å2 ) for a three atomic system with different probe radii (r p = 2 Å, 1 Å), obtained from the MSMS and the proposed method.

		Ases (rp = 2 Å)	Asas (rp = 2 Å)	Ases (rp = 1 Å)	Asas (rp = 1 Å)
	MSMS	63.712	253.122	64.616	153.545
	Proposed method	63.712	253.122	58.857	153.545

Acknowledgments

The authors would like to thank the anonymous referee for the valuable suggestions to improve the paper as well as Pascal Frey, Filippo Lipparini, Yvon Maday and Benedetta Mennucci for the encouraging and fruitful discussions we had in the last months. Ch.Q. acknowledges support from the Fondation Sciences Mathématiques de Paris.

First, we provide a lemma about the geometric relationship between P 0 and T = R(I) as following.

Lemma 1.1. There holds that P 0 \ x∈K B rp (x) ⊂ T , where T = R(I) is the closed region composed of polyhedrons in the SAS-cavity.

Proof. We first recall an equivalence relationship for three sets A, B and C:

In consequence, to prove the lemma is equivalent to prove that P 0 \T ⊂ x∈K B rp (x). For simplicity, we only give the proof for the case of the eSAS in R 2 . There is no essential additional difficulty for the proof in R 3 . We define another polygon called the reduced polygon T 0 contained in T (see Figure A.22), which is obtained by removing from T those vertices which are SAS intersection points, so that only the centers of the SAS-balls are left as vertices. In consequence, we have a simple relationship T = ( xm∈I T xm )∪T 0 , in which T xm is the triangle with the three vertices x m , c i and c j . Each edge of T 0 is associated with an SAS intersection point.

With the above notations, three cases will be discussed as following: a) P 0 ⊂ T xm : Since T xm ⊂ T implies that P 0 \T = ∅, this is a trivial case.

b) P 0 ⊂ T xm ∪ T 0 : Since T xm ∪ T 0 ⊂ T implies that P 0 \T = ∅, this is also a trivial case.


Farida Ouslimani
email: ouslimani.f@gmail.com

Zohra Ameur

Rotation Invariant Features Based on Regional Rank for Texture Classification

Keywords: Texture descriptor, regional rank, rotation invariance, texture classification

In this paper, we present rotation invariant descriptors using regional rank for texture classification. The regional rank presents the rank of the gray level of each pixel in a region whose size and shape depend on the gray level of the treaty pixel and its neighbors. Rotation invariant features are obtained by combining the rank which is found and the treaty pixel code. This latter is calculated by global thresholding. Eight discriminates and rotation invariant features are then obtained. The features size don't increase with scale and kept constant. Tests are performed on the well known Outex database. Compared to LBP method using different schemes, the proposed method achieves good texture classification performance while enjoying a compact feature representation.

INTRODUCTION

Texture analysis is an active research topic in image processing and pattern recognition. However, arbitrary rotation could occur inreal-world textures. It will affect the performance of the texture analysis methods. Thus, extracting texture features that are rotation-invariant is animportant issue to be addressed. Many methods were proposed for rotationinvariance [START_REF] Zhang | Brief review of invariant texture analysis methods[END_REF].Early ones focus on the statistical analysis of texture images. This starts with the co-occurrence matrix method [START_REF] Haralik | Texture features for image classification[END_REF].Later, many model-based methods were proposed for rotation invariance, such as circular autoregressive model [START_REF] Rl | A model-based method for rotation invariant texture classification[END_REF] and markov random fields [START_REF] Cohen | classification of rotated and scaled textured image using Gaussian markov random field models[END_REF][START_REF] Jl | Rotation and gray scale transform invariant texture identification using wavelet decomposition and hidden Markov model[END_REF][START_REF] Deng | Gaussian MRF rotation -invariant features for image classification[END_REF]. Rotation invariant features obtained through filter bank responses are also used [START_REF] Randen | Filtering for texture classification: A comprative study[END_REF][START_REF] Haley | Rotation-invariant texture classification using a complete space -frequency model[END_REF][START_REF] Varma | A statistical approach to texture classification from single images[END_REF]. Recently, Ojala et al. [START_REF] Ojala | Multiresolution Gray-Scale and Rotation Invariant Texture Classification with Local Binary Patterns[END_REF] proposed a local binary pattern (LBP) histogram for rotation invariant texture classification. LBP is a simple yet efficient operator to describe local image pattern, and it has achieved impressive classification results on representative texture databases [START_REF] Ojala | Outex-New framework for empirical evaluation of texture analysis algorithms[END_REF].Many extensions of LBP were presented [START_REF] Lio | Dominant Local Binary Patterns for Texture Classification[END_REF][START_REF] Guo | A completed modeling of local binary pattern operator for texture classification[END_REF][START_REF] Khellah | Texture Classification Using Dominant Neighborhood Structure[END_REF].

In [START_REF] Ouslimani | Texture analysis method using regional rank coding[END_REF], we introduced regional rank coding method. In order to extract more discriminate texture pattern, some modifications were conducted. First, to acquire more precision,regional rank ratio is portioned in four intervals and not two. Second, the gray level information is used for coding after global thresholding. Thus, eight rotation invariant patterns are obtained. Third, different neighbourhood search space dimension are exploited. Also, the features histogram has the same size in different scales. So, they can be easilyfused to improve the texture classification accuracy.

The rest of the paper is organized as follows. Section 2 reviews briefly LBP. Section 3 presents the proposed features. Section 4 reports the experimental results on large public texture database. Finally, Section 5 concludes the correspondence.

BRIEF REVIEW OF LOCAL BINARY PATTERN (LBP)

LBP is initially proposed by Ojala et al. [START_REF] Ojala | A comparative study of texture measures with classification based on features distributions[END_REF] to support the local contrast measure of image.

The spatial structure of a local image texture is characterized by thresholding a 3×3 square neighbourhood with the value of the center pixel. To allow multi-resolution analysis and rotation invariance a more general formulation defined on circular symmetric neighborhood was proposed in [START_REF] Ojala | Multiresolution Gray-Scale and Rotation Invariant Texture Classification with Local Binary Patterns[END_REF] .The LBP pattern is computed by comparing the value of the center pixel 𝑥 𝑐 with those of its P neighbors that are evenly distributed on a circle of radius R centred at center 𝑥 𝑐 , such that the LBP response is calculated as

𝐿𝐵𝑃 𝑅,𝑃 = 𝑠 𝑥 𝑖 -𝑥 𝑐 𝑃-1 𝑖=0 2 𝑖 , 𝑠 𝑥 = 1 𝑥 ≥ 0 0 𝑥 < 0 (1)
The major problem of LBP is the exponential growth of the number of patterns with respect of the neighborhood size. In order to address this problem, Ojala et al. [START_REF] Ojala | Multiresolution Gray-Scale and Rotation Invariant Texture Classification with Local Binary Patterns[END_REF] defined three mapped patterns,which are uniform, rotation invariant and rotation invariant ones, respectively.

The uniform value of an LBP pattern, which is the number of circular spatial transition (bitwise 0/1 or 1/0 changes), can be mathematically computed by

𝑈 𝐿𝐵𝑃 𝑅,𝑃 = 𝑠 𝑥 𝑖 -𝑥 𝑐 -𝑠 𝑥 𝑖+1 -𝑥 𝑐 𝑃-1 𝑖=0 (2)
When 𝑈 𝐿𝐵𝑃 𝑅,𝑃 is less than or equal to 2 i.e U≤2, the pattern is called uniform pattern, denoted as 𝐿𝐵𝑃 𝑅,𝑃 𝑢2 . Uniform patters have P× (P-1) +3 different output values; these patterns are one of the fundamental patterns within image textures [START_REF] Ojala | Multiresolution Gray-Scale and Rotation Invariant Texture Classification with Local Binary Patterns[END_REF].

To achieve rotation invariance, a locally rotation invariance and pattern is defined as follows:

𝐿𝐵𝑃 𝑅,𝑃 𝑟𝑖 = min (𝑅𝑂𝑅 𝐿𝐵𝑃 𝑅,𝑃 , 𝑖 𝑖 = 0,1, … 𝑃 -1) (3)
Where (𝑅𝑂𝑅 𝑥, 𝑖 ) denotes bit-wise right shift on the number 𝑥 𝑖 times. By introducing the definition of rotation invariance, LBP not only has prominent performance for image rotation, but also has fewer patterns [START_REF] Ojala | Multiresolution Gray-Scale and Rotation Invariant Texture Classification with Local Binary Patterns[END_REF].

In order to obtain improved rotation invariance and to further reduce the dimensionality of the LBP histogram feature, building on 𝐿𝐵𝑃 𝑅,𝑃 𝑟𝑖 , Ojala et al [START_REF] Ojala | Multiresolution Gray-Scale and Rotation Invariant Texture Classification with Local Binary Patterns[END_REF] proposed the rotation invariant uniform patterns 𝐿𝐵𝑃 𝑅,𝑃 𝑟𝑖𝑢 2 , the collection of those rotation invariance patterns having a U value of at most 2.

𝐿𝐵𝑃 𝑅,𝑃 𝑟𝑖𝑢 2 = 𝑠 𝑥 𝑖 -𝑥 𝑐 𝑝-1 𝑖=0 , 𝑖𝑓 𝑈(𝐿𝐵𝑃 𝑅,𝑃 ) ≤ 2 𝑃 + 1 𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒 (4)
If P and R are respectively set to 8 and 1 the histogram dimensions of uniform pattern, rotation invariant pattern and rotation invariant uniform pattern are 59, 36 and 10 respectively.

REGIONAL TEXTURE DESCRIPTORS

The regional rank coding is based on the rank of pixel grayscale defined on regional neighborhood. The size and shape of each region depend on the gray level of the treaty pixel.

For each pixel x of the image, we determine the set of pixels forming the tray or plateau containing x noted Pt (x). A plateau is a set of pixels of the same gray level and which are 8 connected (Fig. 1). The regional rank coding assigns the same value to all the pixels of the plateau Pt (x), based on gray levels of its neighbors. The set of neighbors of a plateau Pt (x) be the set V (Pt (x)) defined as follows:

V(Pt(x)) ={y  (Pt(x)) c / B8 (y)  ((Pt(x))  } (5)
Fig. 1 represents a partition of an image, we can see two plateau whose gray level equal to 7

(graycolor) surrounded by all relevant neighbors (blue color). One of the plateau is composed of five pixels surrounded by sixteen neighbors, while the other plateau is composed of three pixels surrounded by nine neighboring pixels.

𝐹 𝑥, 𝑥 𝑣 = 1 𝑖𝑓 𝐼 𝑥 ≥ 𝐼(𝑥 𝑣 ) 𝑂 𝑖𝑓 𝐼 𝑥 < 𝐼(𝑥 𝑣 )

The found rank will not be used directly to encode the texture because it varies depending on the number of neighbors, then we consider the ratio P0

= rank / | | V | |.
Noting that the plateau Pt(x) and the the set of neighbors V(Pt(x) ) remains the same even if image is rotated,so ratio P0 is invariant to rotation.

P0 has a continuous value ranging between 0 and 1, hence a quantization is needed; this was done by adding to gether P0 which belonging to the same intervall :[0; 0.25], [0.25;0.5],

[0.5;0.75]and [0.75; 1]. P0 aquire an important information , for example, if P0 ∈[0; 0.25], the Pt(x) posseds the minima gray level but if P0 ∈ [0.75. 1] , the corresponded plateau Pt(x) gray level is the maxima amoung its neigbors.

The image gray level has discriminate information. So, to calculate the final coding, a value corresponding to global thresholding of the gray level of the treaty pixel iscombined with P0 according to Table I .Thus, eight rotation invariant features are extracted.

M represents the average gray level of the whole image I.

I(x) is the gray level of the treaty pixel.

To acquire more information, other scales can be used: features are extracted by extending the search space around each pixel in square neighbourhood of S×S pixel size . The case explained in this section corresponds to S=3.

Fig. 2 represents a partition of an image using a square neighborhood of 5×5 pixel size (S=5). We can see plateau whose gray level equal to 7 (graycolor) surrounded by all relevant neighbors (blue color).

TABLE I. REGIONAL RANK PATTERN IN FUNCTION OF GRAY LEVEL AND RATIO P0

The features sizes don't increase with scale and kept constant and equaling to 8.RRC S corresponds to the features extracted using S× S neighborhood size.

EXPERIMENTS AND RESULTS

The classification performance of the proposed descriptor is compared with state LBP methods on the well known texture databases Outex [START_REF] Ojala | Outex-New framework for empirical evaluation of texture analysis algorithms[END_REF].The Outex database include two test suites: Outex_TC_00010(TC10) and Outex_TC_00012(TC12).The two suites contain the same 24 classes of textures as shown in Fig. 3.Each texture class was collected under three different illuminants("horizon", "inca", and "t184") and nine different rotation angles (0° ,5° ,10° ,15° ,30° ,45° ,60°, 75° ,90°),there are 20non-overlapping 128*128 texture samples for each class under each situation .The experiment setups are as follows:for TC10, samples of illuminant "inca" and angle 0° in each class were used for classifier training and the other eight rotation angles with the same illuminant were used for testing.Hence, there are 480(24*20) models and 3840(24*8*20) validation samples.For TC12, the classifier was trained with the same training samples as TC10, and it was tested with all samples captured under illuminant "t184" or "horizon".Hence, there are 480(24*20) models and 4320(24*20*9) validation samples for each illuminant.

These databases are challenging due to the large number of texture classes, significant rotation and illumination changes.Each texture image is converted into gray scale and normalized to zero mean and unit standard deviation . The nearest-neighborhood classifier 𝜒 2 distance is adopted for the classification task [START_REF] Khellah | Texture Classification Using Dominant Neighborhood Structure[END_REF].

The proposed method is compared to rotation invariant uniform patterns 𝐿𝐵𝑃 𝑅,𝑃 𝑟𝑖𝑢 2 using different combination. 𝐿𝐵𝑃 𝑅,𝑃 𝑟𝑖𝑢 2 feature size are 10, 18 and 26 for R= 1, 2 and 3 respectively. Second, histograms concatenated over different scales increase the performance of the methods .So, the best classification rats is obtained by concatenated features histograms extracted from the three scales. We can observe that LBP exceeds slightly the proposed method for TC10 and TC12 "t". But our method exceeds LBP for TC12"h" in different cases . Knowing that the fuzzed LBP feature size 54(10+ 18+26) is about two times the feature size of the fuzzed RRC S 24(3×8).

Third,RRC S performance decrease when S increase because the information patterns acquired is more global but less accurate.

CONCLUSION

This paper introduced new rotation invariant texture descriptors based on regional rank coding.

The features are extracted by combining the regional rank order and the pixel gray level information. Eight rotation invariant features are then extracted. Different neighborhood search space dimension are exploited. Also, the features histogram has the same size in different scales.

So, features extracted from each scale are concatenated in a single feature histogram with n*8 feature dimension where n is number of scales.

As demonstrated in the experimental results performed on the Outex database .

RRC 3 outperforms 𝐿𝐵𝑃 1,8 𝑟𝑖𝑢 2 for TC10 and TC12. RRC S features produces the best classification rate for TC12 "h" in different cases .

The concactenated histogram features extracted over different scales performs results. The best rate is obtained by using the tree scales exceeding 94% for TC10 and reaching respectively 88.96% and 89.17% for TC 12"t" and TC 12"h".

We think that RRC 3 features has a large area of perspectives to be exploited. Especially to make features more discriminates and eventually robust to noise.

Fig . 1 .

 1 Fig .1. Representation of plateau and its neighbors

Fig. 2 .

 2 Fig.2. Representation of plateau (green) and its neighbors (blue) using 5× 5 search space .

First, The RRC 3

 3 classification rate exceeds 𝐿𝐵𝑃 1,8 𝑟𝑖𝑢 2 considerably.

Table

  II lists the experimental results by different schemes.Under TC12,"t"represents the test setup of illuminant " t184" and "h" represents "horizon".The best results acquired by the proposed methods and LBP are in bold.From the produced results, we can make the following remarks: