The Weighted Gaussian Curvature Derivative of a Space-Filling Diagram

The morphometric approach [HRC13,RHK06] writes the solvation free energy as a linear combination of weighted versions of the volume, area, mean curvature, and Gaussian curvature of the space-filling diagram. We give a formula for the derivative of the weighted Gaussian curvature. Together with the derivatives of the weighted volume in [EdKo03], the weighted area in [BEKL04], and the weighted mean curvature in [AkEd19], this yields the derivative of the morphometric expression of solvation free energy.


Introduction
This is a companion paper to [7] and [4] on the derivatives of the weighted volume and of the weighed area of a space-filling diagram, and particularly of [1] on the derivative of the weighted mean curvature, with which this paper shares the notation. Together with these three papers, we complete the explicit description of the derivatives of the weighted intrinsic volumes in R 3 . All four derivatives are needed to realize the morphometric approach to the solvation free energy; see [10,11]. Getting its inspiration from the theory of intrinsic volumes, this approach writes the solvation free energy as a linear combination of the weighted volume, the weighted area, the weighted mean curvature, and the weighted Gaussian curvature.
Compare this with Hadwiger's characterization theorem, which asserts that any measure of convex bodies in R 3 that is invariant under rigid motion, continuous, and additive is a linear combination of the four (unweighted) intrinsic volumes [9]. Among the four intrinsic volumes, the Gaussian curvature is distinguished by the Gauss-Bonnet theorem, which asserts that it is 2π times the Euler characteristic of the surface [3]. Translated to a set of spheres in motion, this implies that the curvature function is piecewise constant and thus lacks any interesting derivative. The situation changes drastically in the weighted case, in which the contribution of a sphere to the Gaussian curvature is scaled by a real weight, which represents a physical property of the corresponding atom, such as its hydrophobicity. The main result of this paper is an explicit formula for the derivative of the weighted Gaussian curvature of a space-filling diagram, and an analysis of the cases in which this derivative either does not exist or is not continuous.
Outline. Section 2 reviews the background relevant to this paper. Section 3 discusses the weighted Gaussian curvature of a space-filling diagram. Section 4 derives the constituents of the weighted Gaussian curvature. Section 5 states the derivative in terms of the gradient. Section 6 analyzes the configurations at which the gradient is not continuous. Section 7 concludes this paper.

Background
This section reviews the space-filling diagram and its dual alpha shape, both geometric models of molecules.
Space-filling diagrams. Throughout this paper, we let X be a collection of n closed balls in R 3 , and to remind us of the connection to chemistry and biology, we call the union of these balls a space-filling diagram. Writing x i ∈ R 3 for the center and r i > 0 for the radius, the corresponding ball is B i = {a ∈ R 3 | a − x i ≤ r i }, and the spacefilling diagram is X = n−1 i=0 B i . Its boundary consists of (spherical) patches of the form S i \ j =i int B j , in which S i = bd B i . These patches meet along closed (circular) arcs of the form S ij \ k =i,j int B k , in which S ij = S i ∩ S j . Finally, these arcs meet at corners of the form S ijk \ =i,j,k int B , in which S ijk = S i ∩ S j ∩ S k . For ease of reference, we collect and extend the introduced notation in Table 1, which is given in the appendix.
We are interested in the volume of the space-filling diagram and the area, mean curvature, and Gaussian curvature of its boundary. Since this boundary is not necessarily smoothly embedded in R 3 , we use discrete formulas for the curvatures that agree with the smooth formulas in the limit. For the mean curvature, we have non-zero contributions from the patches and the arcs; see Equation (4), and for the Gaussian curvature, we have non-zero contributions from the patches, the arcs, and the corners; see Equation (5). However, according to the Gauss-Bonnet theorem [3], the Gaussian curvature is invariant under deformations that preserve its topology, and it changes by integer multiples of 2π whenever the topology is not preserved. The derivative under motion of the balls is therefore not very informative, namely generically zero and occasionally undefined. This is in sharp contrast of the weighted case.
Voronoi decomposition and dual alpha complex. We need a decomposition of the space-filling diagram and of its boundary to define what we mean by its weighted measures. Recall that B i is the closed ball with center x i and radius r i . The power distance of a . The Voronoi domain of B i is the set V i of points a ∈ R 3 for which π i (a) ≤ π j (a) for 0 ≤ j < n. It is a possibly unbounded and always closed and convex 3-dimensional polyhedron. Generically, is either empty or a point. The Voronoi domains decompose the space-filling diagram into convex sets of the form B i ∩ V i . The boundary of such a set consist of a spherical patch (which also belong to the boundary of the space-filling diagram) and a collection of flat patches (along which the sets meet). The spherical patch meets the flat patches along circular arcs on the boundary of X, and the flat patches meet along straight line segments inside X.
We use the notation in Table 1 for the fractional measures of the pieces in the decomposition. They are conveniently computed using the alpha complex of X, which is the nerve of the collection of B i ∩ V i [8]. We recall that the Delaunay mosaic is the nerve of the collection of V i , which implies that the alpha complex is a subcomplex of the Delaunay mosaic. Just like the Delaunay mosaic, the alpha complex is geometrically realized by mapping its vertices to the corresponding points in R 3 . The alpha shape is the underlying space of the alpha complex, and we stress that this is a set of points while the alpha complex is a set of simplices.
The boundary of the alpha shape consists of all points for which no open neighborhood is contained in the set. It is the underlying space of a subcomplex of the alpha complex, and we refer to its members as boundary simplices. The alpha complex contains a simplex for every non-empty common intersection of Voronoi domains, which makes it a convenient book-keeping device for the intended computations. To appreciate this facility, it is useful to familiarize ourselves with the details of the correspondence between the boundaries of the space-filling diagram and the alpha shape. Whenever S i contributes a non-zero number of patches, its center, x i , is a boundary vertex of the alpha shape. The number and structure of the patches can be extracted from the structure of simplices in the alpha complex that share x i . Whenever S ij contributes a non-zero number of arcs, the edge connecting x i to x j belongs to the boundary of the alpha shape. The triangles and tetrahedra that share the edge are arranged cyclically around the edge, and the number of gaps between contiguous triangles in the cyclic order is equal to the number of contributed arcs. Finally, whenever S ijk contributes a non-zero number of corners -namely one or two -the triangle with vertices x i , x j , x k belongs to the boundary of the alpha shape. The number of contributed corners is equal to the number of sides of the triangle that face the outside of the alpha shape, which is again one or two.

Weighted Gaussian Curvature
According to the Gauss-Bonnet theorem [3], the (integrated) Gaussian curvature of a closed surface, F, that is smoothly embedded in R 3 is a fixed multiple of the Euler characteristic of the surface: in which κ 1 (x) and κ 2 (x) are the two principal curvatures at x. The boundary of a space-filling diagram is generally not smoothly embedded, but the extension of the Gaussian curvature used in this paper still satisfies (1). This implies that the (unweighted) Gaussian curvature is constant while the spheres move, until a time at which the surface changes topologically and the Gaussian curvature jumps to a possibly different value. In other words, the derivative of the unweighted Gaussian curvature is generally zero and undefined whenever the topology changes. The weighted case is radically different because the contributions of the individual spheres to the weighted Gaussian curvature to not cancel each others changes. In contrast to the unweighted case, the computation of the derivative of the weighted Gaussian curvature of a space-filling diagram is a meaningful as well as challenging undertaking.
Weighted intrinsic volume. Since the boundary of a space-filling diagram is only piecewise smooth, we compute its weighted Gaussian curvature separately for the sphere patches, of the circular arcs, and the corners. To facilitate a comparison, it is convenient to give such formulas for the weighted versions of all four intrinsic volumes in R 3 : the weighted volume, the weighted area, the weighted mean curvature, and the weighted Gaussian curvature. Recall that X is the union of the balls B i , for 0 ≤ i < n. Its state, x ∈ R 3n , is the concatenation of the center vectors. In other words, the (3i + )-th coordinate of x is the -th coordinate of x i . Dividing the space-filling diagram with the Voronoi tessellation, we use the notation introduced in Table 1 to give formulas for the weighted intrinsic volumes.
Proposition 1 (Weighted Intrinsic Volumes). The weighted volume, surface area, mean curvature, and Gaussian curvature of the space-filling diagram, X, are The weighted Gaussian curvature function decomposes into three sums, which account for the sphere patches, the circular arcs, and the corners, respectively. The sums are over ordered index sequences. This implies that the same three indices appear six times, which explains the coefficient of 1 3 in front of this sum. We will see that φ i,jk = α i φ ijk is the area of a spherical quadrangle, so we can simplify the last sum. Formulas for the derivatives of σ i and of σ ij can be found in [1]. Furthermore, λ ij , σ ijk , φ i,jk depend on the state. We write λ ij (τ ) for λ ij at state x + τ t and λ ij for the derivative of λ ij at τ = 0, and similarly for σ ijk and φ i,jk . Derivatives with respect to parameters other than τ are explicitly stated as such.
How to split a corner. Every corner of the space-filling diagram has a non-zero contribution to the Gaussian curvature, and since the spheres have weights, we need to decide how to distribute this contribution among the generically three spheres that meet at the corner. We do this in a way that is consistent with splitting the curvature concentrated along an arc in equal halves, which is motivated by the Appolonius diagram of the spheres; see [1].
Consider spheres S i , S j , S k , and suppose that P is one of the two points in which the spheres generically intersect. Write n i , n j , n k for the unit outward normals of the spheres at P , and recall that φ ijk is the area of the spherical triangle with vertices n i , n j , n k on the unit sphere; see Figure 1. We split the contribution of the corner using α i + α j + α k = 1, Figure 1 A spherical triangle with vertices ni, nj, n k and the cap with center z whose boundary is the unique circle that passes through the three vertices. The spherical triangle is decomposed into three spherical quadrangles by connecting z to the midpoints of the three sides, and into six cones by further connecting z to the three vertices.
in which the three coefficients are the fractions of the areas of the spherical quadrangles obtained by decomposing the spherical triangle with the great-circle arcs that connect the center of the circumcircle, denoted z, with the midpoints of the sides, denoted m i , m j , m k . By connecting z with the three vertices, we further decompose each quadrangle into two spherical triangles, which we refer to as cones of z over the half-sides; see again Figure 1. We remark that z can lie outside the spherical triangle, in which case the areas of two of the cones are considered to be negative. The cones come in pairs of equal area, and each pair glues together to form an isosceles spherical triangle. As remarked, exactly one of the three isosceles triangles has negative area if z is not contained in the original spherical triangle. Since each of the three quadrangles is the union of two cones, we can write its area as half the sum of the area of two isosceles triangles: and similarly for Area(n j m i zm k ) and Area(n k m j zm i ). Equivalently, the area of the quadrangle is half of Area(n i n j n k ) − Area(n j n k z), but we prefer (6) for our computations.

Derivatives
We reuse the expressions for σ i and σ ij in [1], and we develop formulas for λ ij , σ ijk , φ i,jk . The middle of the three, σ ijk , is piecewise zero because σ ijk is piecewise constant. As will be discussed in Section 6, σ ijk changes at non-generic states, where its derivative is not defined. We first consider the derivative of λ ij and second that of φ i,jk .

Derivative of λ ij
Given spheres S i and S j , we recall that λ ij is the combined length of the unit normal vectors at a point P ∈ S i ∩ S j after projection to the line that passes through the centers, x i and x j ; see Figure 2. If the two vectors point in the same direction -after projection -then we take one length negative, so more precisely, λ ij is the distance between the projections of the endpoints of the unit vectors. To compute the derivative of λ ij , it will be useful to recall After projecting the unit normal vectors at an intersection point, λij is the distance between the two endpoints.
that the signed distance of x i from the Voronoi plane of the two spheres is see [1]. The length of the unit normal vector of S i at P , after projection, is ξ i /r i , and since ξ i is signed so is this length, as required. Adding the two signed lengths, we get The derivative with respect to the distance between the two centers is therefore To translate this result to the derivative with respect to the motion vectors t i of x i and t j of The derivative of the combined length of the unit normal vectors at an intersection point of spheres S i and S j -after projection to the line passing through the centers -at state x with momentum t is in which the coefficient is given in (9).

Derivative of φ i,jk
Instead of deriving α i and φ ijk separately, we recall that φ i,jk = α i φ ijk is the area of a spherical quadrangle, which we can derive directly. We begin with intrinsic formulas for the area and the radius of a spherical triangle.
Area of spherical triangle. Recall that φ ijk is the area and φ ij , φ jk , φ ki are the lengths of the sides of the spherical triangle with vertices n i , n j , n k . Writing s = 1 2 (φ ij + φ jk + φ ki ) for half the perimeter, we have the following formula for the sine of half the area: see [5, page 230, equation (225)]. To rewrite this formula, we introduce the squared cosines of the half-lengths as parameters. Geometrically, cos 2 x 2 is the distance from the center to the midpoint of the straight edge that connects the endpoints of an arc of length x on the unit sphere. The square of this distance conveniently relates to Pythagoras' theorem for right-angled triangles. In a first step, we rewrite the product of the four sines.
We refer to Appendix A for a proof. We can thus rewrite (11) to give the area of any spherical triangle given the squared cosines of its half-lengths: By (11) and (12), we have φ ijk = S(a, b, c) with parameters as given in Formula 3. We introduce notation for the areas of the isosceles triangles obtained by gluing equal-area cones: A(a, r) = S(a, r, r), B(b, r) = S(r, b, r), C(c, r) = S(r, r, c), with r = cos 2 R ijk 2 , in which R ijk is the radius of the cap bounded by the circle that passes through n i , n j , n k . By symmetry between the arguments of S, all three functions are the same, but it is convenient to distinguish them so we do not lose the connection to the sides of the spherical triangle. Using (6), we get the areas of the spherical quadrangles: in which the signs depend on whether z lies inside or outside the triangle spanned by n i , n j , n k . Specifically, sgm i,jk = +1 if z and n i lie on the same side of the great-circle that passes through n j and n k , and sgm i,jk = −1 if the two points lie on opposite sides of this great-circle.
In the boundary case, when z lies on the great-circle, n i , n j , n k form a Euclidean right-angled triangle. Pythagoras' theorem for its edges is equivalent to sin 2 φij 2 + sin 2 φ ki 2 = sin 2 φ jk 2 . If z and n i lie on the same side of the great-circle, then the left-hand side is strictly larger than the right-hand side. Substituting sin 2 x = 1 − cos 2 x, we get Since S(a, b, c) is symmetric in its three arguments, it suffices to compute the derivative with respect to one of them. To this end, we recall a few basic rules of differentiation, namely d arcsin v(a, b, c) = 4abc, we get du/da = 4bc − 2(a + b + c − 1) and dv/da = 4bc. Furthermore setting w = u v and t = √ w, we have S(a, b, c) = 2 arcsin t and therefore Multiplying the three right-hand sides and simplifying the result, we get the derivative of the area with respect to the squared cosine of one half-length: Permuting the three arguments, we get symmetric expressions for dS/db and for dS/dc.

Radius of spherical triangle.
We need a formula for the radius of the cap, R ijk , in the same parametrization: tan R ijk = tan To turn this into a formula for the squared cosine of the half-radius, we recall two basic trigonometric identities: cos 2 x 2 =  r(a, b, c) r(a, b, c) is symmetric in the three arguments, it suffices to compute the derivative with respect to one of them. Referring to (25), we set which give du/da = 4bc − 2(a + b + c − 1) and dv/da = 2(−a + b + c − 1). Observing that r(a, b, c) = 1 2 + 1 2

u(a, b, c)/v(a, b, c), we see that its derivative with respect to the first argument is (vdu
Symmetric expressions hold for dr/db and for dr/dc.
Using the chain rule. We get the derivatives of the areas of the spherical quadrangles by combining the results with the chain rule. Beyond the relations above, we need the derivative of the length of a side, which we find in [1]: We recall that d cos 2 x 2 /dx = − sin x 2 cos x 2 , and we use (14) to express the area of the spherical quadrangle: φ i,jk = 1 2 [sgm k,ij A(a, r) + sgm j,ki C(c, r)], in which a = cos 2 φij 2 , b = cos 2 φ jk 2 , c = cos 2 φ ki 2 , and r = cos 2 R ijk 2 as given in (25). The derivatives of A with respect to its arguments are dA/da = dS/da and dA/dr = dS/db + dS/dc, both evaluated at b = c = r. The derivatives with respect to the distances between the centers of the spheres are therefore (31) and (32) are symmetric because the same two sides of the isosceles triangle depend on x j − x k and on x k − x i , while (30) is different because all three sides depend on x i − x j . We can now summarize and state the derivative of φ i,jk with respect to the motion vector.
Lemma 4 (Derivative of φ i,jk ). The derivative of the area of the spherical quadrangle with vertices n i , m k , z, m j on the unit sphere is in which the derivatives on the right-hand side are given in or symmetric to (30), (31), (32).

Gradients
For computational purposes, it is convenient to write the derivative of the weighted Gaussian curvature function, gauss : R 3n → R, in terms of the gradient of gauss at x ∈ R 3n , denoted g = ∇ x gauss. Recalling (5), this derivative is gauss = d + e + f + h , with Writing g = [g 1 , g 2 , . . . , g 3n ] T , we recall that g i = [g 3i+1 , g 3i+2 , g 3i+3 ] T is the 3-dimensional gradient that applies to x i . Using boldface letters for the gradients of d, e, f , h, and similar conventions for the 3-dimensional sub-vectors, we have g = d + e + f + h and We get the gradients by redistributing the relevant derivatives given in [1] and stated in Lemmas 2 and 4.
First term. The derivative of σ i is given in Lemma 6 of [1] as in which the first sum is over the boundary edges of the alpha complex incident to x i , and the second sum is over the triangles incident to these edges. We can therefore rewrite (39) as in which the first sum ranges over all (directed) boundary edges of the alpha shape, and second sum ranges over all triangles incident to these edges. We get d i by accumulating the relevant contributions. Using u ji = −u ij and u jik = u ijk , we get in which the first sum is over all boundary edges incident to x i , and the second sum is over all triangles incident to these edges. Compare (47) with the weighted area gradient given in [4] and with the first term of the weighted mean curvature gradient given in [1].
Second term. This case is complicated by the change of motion needed to compute the derivative of σ ij . Following the exposition in [1, Section 4], we introduce some notation before we can rewrite (40).
we introduce three vectors so that the scalar product of the new motion vector with x k − P , can be written as a sum of three scalar products, each involving only one of the original motion vectors: In addition, we recall the derivatives of σ ij , r ij , and α P from [1]: in which α P is the angle parametrizing the motion of the corner, P , along the circle, S ij . With this, we rewrite (40) as in which the sum is over the three ordered versions of all oriented boundary triangles of the alpha shape, and P is the corresponding corner of the space-filling diagram. Redistributing the terms, we get in which the sum is again over the three ordered versions of all oriented boundary triangles of the alpha shape that share x i .
Third term. Using Lemma 2, we rewrite (41) as in which the sum ranges over all (directed) boundary edges of the alpha shape. We get f i by accumulating the relevant contributions, as before: in which the sum is over boundary edges incident to x i .

Fourth term. Using Lemma 4, we rewrite (42) as
in which the coefficients are defined in and after Lemma 4, and the sum is over all (ordered) triangles in the boundary of the alpha shape. We get h i by accumulating the relevant contributions: in which the sum is over all ordered boundary triangles of the alpha shape that are incident to x i . To see this, we observe that the three ordered versions of the same orientation, ijk, jki, kij, all contribute the same term. By symmetry of the spherical quadrangles defined for the opposite orientation, the corresponding three ordered versions, ikj, kji, jik, also contribute the same term.
Summary. We finally get the gradient of the weighted Gaussian curvature function by adding the gradients of the four component functions:

Continuity
Just as the derivatives of the other intrinsic volumes of a space-filling diagram, the derivative of the weighted Gaussian curvature is almost everywhere but not everywhere continuous. Indeed, the situation is sufficiently similar to the mean curvature case that we can follow the analysis in [1] and be brief.
General position. The derivative of gauss : R 3n → R is continuous everywhere except at non-generic states. To define what this means, we recall that x ∈ R 3n encodes a collection, X, of n closed balls in R 3 . We say X is in general position and, equivalently, x is generic, if the following two conditions are satisfied: I. the common intersection of p + 1 Voronoi domains is either empty or a convex polyhedron of dimension 3 − p; II. the common intersection of p + 1 spheres bounding balls in X is either empty or a sphere of dimension 2 − p.
Condition I implies that the Delaunay mosaic of X is a simplicial complex in R 3 . When the balls move, the mosaic undergoes combinatorial changes, and these happen at states that violate Condition I. In addition, Condition II captures states at which simplices in the mosaic are added to or removed from the alpha complex. Each violation of the two conditions corresponds to a submanifold of dimension at most 3n − 1 in R 3n . Write M I and M II for the unions of submanifolds that correspond to violations of Condition I and of Condition II, respectively. Since there are only finitely many submanifolds, and at least some of each type have dimension 3n − 1, M I and M II have dimension 3n − 1 each.
Combinatorial changes. The submanifolds decompose R 3n into cells, and we refer to the connected components of R 3n \(M I ∪M II ) as chambers. As long as the trajectory of states stays inside a chamber, the Delaunay mosaic and the alpha complex are invariant, combinatorially, and so are the formulas for the weighted Gaussian curvature and its derivative. The functions defined by the formulas are continuous, which implies that gauss and ∇gauss restricted to a single chamber are continuous.
When the trajectory passes from one chamber to another, it necessarily crosses a submanifold, and we let x be the corresponding non-generic state. If x ∈ M I , then the associated combinatorial change is a flip in the Delaunay mosaic; see [6] for details. As argued in [1], mean and ∇mean are both continuous at x, and the same argument implies that also gauss and ∇gauss are continuous at x. If x ∈ M II , then the associated combinatorial change is the addition of an interval of simplices in the Delaunay mosaic to the alpha complex or, symmetrically, the removal of such an interval from the alpha complex. Following [1], we distinguish between singular intervals that contain only one simplex and non-singular intervals that contain two or more simplices; see [2] for background on the discrete Morse theory of Delaunay mosaics. A singular interval corresponds to two, three, or four spheres whose common intersection is a single point. If this point happens to belong to the boundary of the space-filling diagram, then the operation changes the topology type of the surface. By the Gauss-Bonnet theorem, this implies that in the unweighted case, gauss changes by a non-zero integer multiple of 2π, and ∇gauss is undefined. Similarly, each non-critical interval corresponds to two, three, or four spheres whose common intersection is a single point, with the difference that now one of the spheres is contained in the union of the balls bounded by the other one, two, or three spheres. There are six basic cases, each with its own implications on the continuity of ∇gauss. The details are not important, so we refer to the analogous analysis in [1].
Theorem 6 (Continuity of Gradient). The gradient of the weighted Gaussian curvature of a space-filling diagram of n closed balls in R 3 is continuous provided the state x ∈ R 3n of the diagram does not belong to M II , which is a (3n − 1)-dimensional subset of R 3n .

Discussion
The main contribution of this paper is the analysis of the derivative of the weighted Gaussian curvature of a space-filling diagram. Specifically, we give an explicit description of the gradient of the weighted Gaussian curvature function, which for n spheres is a map from R 3n to R. In addition, we study the subset of R 3n at which the derivative violates continuity. In summary, this is sufficient information for an efficient and robust implementation of the weighted Gaussian curvature derivative, one that can be added to the inner loop of a molecular dynamics simulation of a physical system.

A Proof of Formula 3
We present a proof of the formula from Section 4 after restating it for convenience.

B Notation
Si = bd Bi sphere bounding ball Sij = bd Bij circle bounding disk S ijk = bd B ijk pair of points bounding line segment xi, xij, x ijk centers of Si, Sij, S ijk ri, rij, r ijk radii of Si, Sij, S ijk uij =