The Weighted Mean Curvature Derivative of a Space-Filling Diagram

TL;DR

This paper derives a formula for the derivative of the weighted mean curvature of a space-filling diagram, enabling more accurate calculations of solvation free energy in molecular simulations.

Contribution

It provides a new formula for the derivative of weighted mean curvature, complementing existing derivatives for volume, area, and Gaussian curvature, to improve morphometric energy calculations.

Findings

Derived the formula for the derivative of weighted mean curvature.

Combined with existing derivatives for other geometric measures.

Facilitates more precise computation of solvation free energy.

Abstract

Representing an atom by a solid sphere in $3$-dimensional Euclidean space, we get the space-filling diagram of a molecule by taking the union. Molecular dynamics simulates its motion subject to bonds and other forces, including the solvation free energy. The morphometric approach [HRC13,RHK06] writes the latter 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 mean curvature. Together with the derivatives of the weighted volume in [EdKo03], the weighted area in [BEKL04], and the weighted Gaussian curvature [AkEd19], this yields the derivative of the morphometric expression of the solvation free energy.

Tables

Table 1. Table 1: The order of change of the weighted mean curvature and of its gradient at states that violate Condition II.

	${𝐂}_{\mathrm{𝟏}}$	${𝐂}_{\mathrm{𝟐}}$	${𝐂}_{\mathrm{𝟑}}$	${𝐍}_{\mathrm{𝟎𝟏}}$	${𝐍}_{\mathrm{𝟎𝟐}}$	${𝐍}_{\mathrm{𝟎𝟑}}$	${𝐍}_{\mathrm{𝟏𝟐}}$	${𝐍}_{\mathrm{𝟏𝟑}}$	${𝐍}_{\mathrm{𝟐𝟑}}$
$\Delta \mathrm{𝑚𝑒𝑎𝑛}$	$\sqrt{\epsilon }$	$\sqrt{\epsilon }$	$\epsilon$	$\epsilon$	$\sqrt{\epsilon }$	$\epsilon$	$\sqrt{\epsilon }$	$\epsilon$	$\epsilon$
$\Delta \Vert \nabla \mathrm{𝑚𝑒𝑎𝑛}\Vert$	$\infty$	$\infty$	$1$	$1$	$\infty$	$1$	$\infty$	$1$	$1$

Table 2. Table 2: Notation for concepts, sets, functions, vectors, variables.

$S_i=\mathrm{bd}{B}_i$	sphere bounding ball
$S_{ij}=\mathrm{bd}{B}_{ij}$	circle bounding disk
$S_{ijk}=\mathrm{bd}{B}_{ijk}$	pair of points bounding line segment
$x_i,x_{ij},x_{ijk}$	centers of $S_i$, $S_{ij}$, $S_{ijk}$
$r_i,r_{ij},r_{ijk}$	radii of $S_i$, $S_{ij}$, $S_{ijk}$
${𝐮}_{ij}=\frac{x_i-x_j}{\Vert x_i-x_j\Vert }$	unit vector between centers
${𝐮}_{ijk}=\frac{{𝐮}_{ik}-⟨{𝐮}_{ik},{𝐮}_{ij}⟩{𝐮}_{ij}}{\Vert {𝐮}_{ik}-⟨{𝐮}_{ik},{𝐮}_{ij}⟩{𝐮}_{ij}\Vert }$	unit normal to ${𝐮}_{ij}$ with positive component in direction ${𝐮}_{ik}$
${\nu }_i=\frac{\mathrm{𝑉𝑜𝑙𝑢𝑚𝑒}\left(B_i\cap V_i\right)}{\mathrm{𝑉𝑜𝑙𝑢𝑚𝑒}\left(B_i\right)}$	volume fraction of ball
${\nu }_{ij}=\frac{\mathrm{𝐴𝑟𝑒𝑎}\left(B_{ij}\cap V_{ij}\right)}{\mathrm{𝐴𝑟𝑒𝑎}\left(B_{ij}\right)}$	area fraction of disk
${\nu }_{ijk}=\frac{\mathrm{𝐿𝑒𝑛𝑔𝑡ℎ}\left(B_{ijk}\cap V_{ijk}\right)}{\mathrm{𝐿𝑒𝑛𝑔𝑡ℎ}\left(B_{ijk}\right)}$	length fraction of line segment
${\nu }_{ijk\ell }=\frac{\#\left(B_{ijk\ell }\cap V_{ijk\ell }\right)}{\#\left(B_{ijk\ell }\right)}$	$0$ or $1$
${\sigma }_i=\frac{\mathrm{𝐴𝑟𝑒𝑎}\left(S_i\cap V_i\right)}{\mathrm{𝐴𝑟𝑒𝑎}\left(B_i\right)}$	area fraction of sphere
${\sigma }_{ij}=\frac{\mathrm{𝐿𝑒𝑛𝑔𝑡ℎ}\left(S_{ij}\cap V_{ij}\right)}{\mathrm{𝐿𝑒𝑛𝑔𝑡ℎ}\left(S_{ij}\right)}$	length fraction of circle
${\sigma }_{ijk}=\frac{\#\left(S_{ijk}\cap V_{ijk}\right)}{\#\left(S_{ijk}\right)}$	$0$, $\frac{1}{2}$, or $1$
$W=W_{\mathrm{elec}}+W_{\mathrm{np}}$	effective solvation potential
$K,K_r=K\oplus r{𝔹}^3$	convex body, parallel body
$H\in {𝒢}_i^3,\chi \left[K\cap H\right]$	$i$-plane, Euler characteristic of intersection
$X=\left\{B_i\right\},\bigcup X$	set of balls, space-filling diagram
${𝕊}^2,{𝔹}^3,{ℝ}^3;f:{𝕊}^2\to ℝ$	unit sphere, unit ball, Euclidean space; height function
$𝐱,𝐭;𝐯,𝐚$	state, momentum; gradients of volume, area
$F:{ℝ}^{3n}\to ℝ$	intrinsic volume function
$v_{ij},{\overline{v}}_{ij},a_{ij},a_{ijk}$	constants
$x,a,b,c,\gamma$	abstract variables
${\xi }_i,{\xi }_j$	signed distances
${𝐭}_k,{𝐓}_{ijk}$	velocity vectors
${𝝎}_j;{𝐭}_j^{\mathrm{str}},{𝐭}_j^{\mathrm{rot}}$	angular momentum; components of motion
$A;P;v;{𝐮}_{ij}^P;{\alpha }^P$	area; point; speed; vector; angle
${𝐕}_{ij},{𝐕}_{ijk},v_k,h_k$	auxiliary variables
${\mathrm{𝑚𝑒𝑎𝑛}}^{\prime }=p^{\prime }+q^{\prime }+s^{\prime }$	derivatives
$𝐩,𝐪,𝐬;{𝐩}_i,{𝐪}_i,{𝐬}_i$	gradients; as they apply to $x_i$
${𝐚}_{ijk},{𝐛}_{ijk},{𝐜}_{ijk},{𝐝}_{ijk}$	auxiliary vectors
$p_{ij},p_{ijk},q_{ij},{\overline{q}}_{ij},s_{ijk},{\overline{s}}_{ijk}$	constants
$\angle x_i{x}_j,\angle x_i{x}_j{x}_k;\epsilon$	angle parameters; motion parameter
${𝕄}_I,{𝕄}_{\mathrm{II}}$	subsets of ${ℝ}^{3n}$

Equations

$$\displaystyle{\mbox{$\it Volume$}}(K_{r})$$

$$a_i$$

$$W_{\mathrm{np}}$$

$${\sigma }_i$$

$$\mathbf{x})$$

$$\displaystyle{\mbox{$\it area$}}({\bf x})$$

$$\displaystyle{\mbox{$\it mean$}}({\bf x})$$

$$\displaystyle{\mbox{$\it gauss$}}({\bf x})$$

$$\displaystyle{\mbox{$\it mean$}}({\bf x})$$

$$\displaystyle{\mbox{$\boldsymbol{\omega}$}}_{j}$$

$${\xi }_i$$

$${\sigma }_i$$

$$\frac{\mathrm{d}{\sigma }_i}{\mathrm{d}\Vert x_i-x_j\Vert }$$

$$\frac{\mathrm{d}{\sigma }_i}{\mathrm{d}\angle x_i{x}_j}$$

$${\sigma }_i^{\prime }$$

$$A$$

$$\frac{\mathrm{d}A}{\mathrm{d}\Vert x_i-x_j\Vert }$$

$$=\frac{1}{8A}\Vert x_i-x_j\Vert \left(r_i^2+r_j^2-{\Vert x_i-x_j\Vert }^2\right).$$

$$\frac{\mathrm{d}{r}_{ij}}{\mathrm{d}\Vert x_i-x_j\Vert }$$

$$=\frac{2{\Vert x_i-x_j\Vert }^2\left(r_i^2+r_j^2-{\Vert x_i-x_j\Vert }^2\right)-16A^2}{8A{\Vert x_i-x_j\Vert }^2}=\frac{{\left(r_i^2-r_j^2\right)}^2-{\Vert x_i-x_j\Vert }^4}{8A{\Vert x_i-x_j\Vert }^2}$$

$$=\frac{{\left(r_i^2-r_j^2\right)}^2-{\Vert x_i-x_j\Vert }^4}{2{\Vert x_i-x_j\Vert }^2\sqrt{2\left(r_i^2{r}_j^2+(r_i^2+r_j^2){\Vert x_i-x_j\Vert }^2\right)-\left(r_i^4+r_j^4+{\Vert x_i-x_j\Vert }^4\right)}}.$$

$${\varphi }_{ij}$$

$$\frac{\mathrm{d}{\varphi }_{ij}}{\mathrm{d}\Vert x_i-x_j\Vert }$$

$$=\frac{2\Vert x_i-x_j\Vert }{\sqrt{2(r_i^2+r_j^2){\Vert x_i-x_j\Vert }^2-(r_i^2-r_j^2{)}^2-{\Vert x_i-x_j\Vert }^4}}.$$

$$r_{ij}^{\prime }$$

$$x_{ij}$$

$$=\frac{\left({\Vert x_i-x_j\Vert }^2-r_i^2+r_j^2\right)x_i+\left({\Vert x_i-x_j\Vert }^2+r_i^2-r_j^2\right)x_j}{2{\Vert x_i-x_j\Vert }^2}.$$

$$\frac{\mathrm{d}{x}_{ij}}{\mathrm{d}\Vert x_i-x_j\Vert }$$

$$\displaystyle{\mbox{$\bf u$}}_{ij}^{P}$$

$$\frac{\mathrm{d}{\sigma }_{ij}}{\mathrm{d}\angle x_i{x}_j{x}_k}$$

$$\frac{\mathrm{d}{\alpha }^P}{\mathrm{d}{r}_{ij}}$$

$${\sigma }_{ij}^{\prime }$$

$$p^{\prime }$$

$$q^{\prime }$$

$$s^{\prime }$$

Full text

††volume: (

1 Introduction

This paper makes a significant step toward turning the morphometric approach to modeling the solvation free energy in molecular dynamics into practice. We recall that molecular dynamics uses Newton’s second law of motion to simulate the dynamic behavior of molecules. We focus on the case in which the motion is computed within an implicit solvent model, in which the effect of water is captured by an effective solvation potential, $W={W_{\!{\rm elec}}}+{W_{\!{\rm np}}}$. The first term on the right-hand side accounts for the electrostatic polarization, and the second term for the van der Waals interactions and the formation of a void in the solvent. Referring to a large body of work on ${W_{\!{\rm elec}}}$ [19], we focus on the non-polar or hydrophobic effect of water, namely ${W_{\!{\rm np}}}$. In an effort to quantify this effect, Lee and Richards introduced the solvent-accessible surface of a molecule [15]. It is the boundary of the space-filling diagram in which each atom is represented by a ball whose radius is its van der Waals radius plus $1.4$ Angstrom for the approximate radius of a water molecule. Based on this concept, Eisenberg and McLachlan defined the solvation free energy as a weighted sum of the surface area, ${W_{\!{\rm np}}}=\sum_{i}{{w_{i}}}A_{i}$, in which ${{w_{i}}}$ is the atomic solvation parameter and $A_{i}$ is the area of the sphere modeling the $i$-th atom that is accessible to the solvent [10]. There is, however, disagreement about the precise formulation fueled in part by evidence that for a small solute, ${W_{\!{\rm np}}}$ is more closely related to the solvent-excluded volume rather than the solvent-accessible area [16, 20]. To resolve this controversy, the morphometric approach of Mecke and Roth [12, 13, 14, 17, 18] suggests that ${W_{\!{\rm np}}}$ be a linear combination of the intrinsic volumes, which in ${{\mathbb{R}}}^{3}$ are the volume, area, mean curvature, and Gaussian curvature. This approach is inspired by the mathematical theory of intrinsic volumes, which we now summarize.

To begin, we let $K\subseteq{{\mathbb{R}}}^{3}$ be a compact convex body and write $K_{r}=K\oplus r{{\mathbb{B}}}^{3}$ for the parallel body obtained by thickening $K$ by $r$ all around. As proved by Jakob Steiner, the volume of the parallel body is described by a degree-$3$ polynomial,

[TABLE]

in which $a_{0}$, $a_{1}$, $a_{2}$, and $a_{3}$ are the volume, surface area, mean curvature, and one third of the Gaussian curvature of $K$ [21]. Scaled and re-indexed versions of the coefficients are referred to as the intrinsic volumes of $K$. According to Hugo Hadwiger’s characterization theorem, every measure of $K$ that is invariant under rigid motions, continuous, and additive is a linear combination of the intrinsic volumes [11]. Since the solvent-accessible model is obtained by thickening the van der Waals model, the analogy seems compelling, except that molecules are generally not convex. To bridge this gap, we mention a result by Morgan Crofton, which states that the $i$-th coefficient of the Steiner polynomial satisfies

[TABLE]

in which $c_{0}=1$, $c_{1}=\tfrac{2}{\pi}$, $c_{2}=1$, $c_{3}=\tfrac{4\pi}{3}$ [5]. To explain (2), we note that ${{\mathcal{G}}_{i}^{3}}$ is the Grassmannian that consists of all $i$-planes passing through the origin in ${{\mathbb{R}}}^{3}$, $H$ is such an $i$-plane, $x$ is a point on the orthogonal $(3-i)$-plane, $H+x$ is the $i$-plane parallel to $H$ that passes through $x$, and $\chi$ is the Euler characteristic. Importantly, Crofton’s integral formula (2) is not restricted to convex bodies and can thus be used to extend the definition of intrinsic volume to non-convex bodies. Returning to the modeling of molecules, we embrace the morphometric approach. In other words, we write

[TABLE]

in which the $a_{i}$ are the volume, area, mean curvature, and one third of the Gaussian curvature of the solvent-accessible model, and physical meanings of the $\mu_{i}$ can be found in [13, 18]. To model physical properties, such as hydrophobicity, we consider weighted versions of the four measures. To embed this approach into molecular dynamics software, we need the derivative of ${W_{\!{\rm np}}}$ with respect to the motion or, equivalently, the derivatives of the weighted volume, the weighted area, the weighted mean curvature, and the weighted Gaussian curvature. The first two of those derivatives have been studied in [8] and [4], where we find formulas based on the Voronoi decomposition of the space-filling diagram as well as characterizations of the configurations at which the derivative is not continuous. The main result of this paper is a similar analysis of the weighted mean curvature derivative, and we refer to [2] for the last piece of the puzzle, which is the weighted Gaussian curvature derivative.

Outline. Section 2 reviews the background relevant to this paper. Section 3 derives the constituents of the weighted mean curvature of a space-filling diagram. Section 4 states the derivative in terms of the gradient. Section 5 analyzes the configurations at which the gradient is not continuous. Section 6 concludes this paper.

2 Background

The background needed for the technical results reported in this paper include two classic theorems in geometry attributed to Archimedes and to Heron of Alexandria, the concept of a space-filling diagram to represent a molecule, the corresponding alpha complex to do the book-keeping, and the analysis of the weighted volume and weighted area derivatives, which are heavily based on these concepts.

Two area formulas. Write ${{\mathbb{S}}}^{2}$ for the unit sphere in ${{\mathbb{R}}}^{3}$ and recall that its area is $4\pi$. Accordingly, the area of a sphere with radius $r$ is $4\pi r^{2}$. Spheres in ${{\mathbb{R}}}^{3}$ enjoy a property discovered thousands of years ago by Archimedes, which is unique to three dimensions. To describe it, we use Cartesian coordinates and let $f\colon{{\mathbb{S}}}^{2}\to{{\mathbb{R}}}$ map every point to its third coordinate, which we call its height. By Archimedes, the area of the cap consisting of all points at height at most $-1\leq a\leq 1$ depends linearly on $a$ and is therefore $4\pi$ times $1+a$ divided by the diameter of the unit sphere, which is $2$. {formula}[Archimedes’ Theorem]

For every $-1\leq a\leq 1$, the area of the cap $f^{-1}[-1,a]\subseteq{{\mathbb{S}}}^{2}$ is $2\pi(1+a)$.

A slice of the sphere can be written as the preimage of an interval, $[a,b]\subseteq[-1,1]$. The area of $f^{-1}[a,b]$ is the area of $f^{-1}[-1,b]$ minus the area of $f^{-1}[-1,a]$, which is $2\pi(b-a)$. We note that the area of the slice thus depends on the width, $b-a$, but is independent of the location of the interval within $[-1,1]$.

A formula that gives the area of a triangle in terms of the lengths of the three edges is attributed to Heron of Alexandria, although it has been suggested that Archimedes knew about the formula two centuries earlier. {formula}[Heron’s Theorem]

The area of a triangle with edges of length $a,b,c$ is $A=\sqrt{s(s-a)(s-b)(s-c)}$, in which $s=\tfrac{1}{2}(a+b+c)$.

There are various proofs of this formula and a number of equivalent statements. We will use $A=\tfrac{1}{4}\sqrt{2(a^{2}b^{2}+b^{2}c^{2}+a^{2}c^{2})-(a^{4}+b^{4}+c^{4})}$.

Space-filling diagrams. Mathematically, such a diagram is merely the union of finitely many closed balls in ${{\mathbb{R}}}^{3}$. Its significance derives from the use as a geometric model of biomolecules [15], such as proteins, DNA, and the like. There are competing models, such as the level sets of a sum of isotropic Gaussian distributions, but we focus on the geometrically more concrete space-filling diagrams. This section introduces the concept along with most of the notation needed for its discussion in the technical sections.

Throughout this paper, we assume a set of $n$ closed balls in ${{\mathbb{R}}}^{3}$, $X=\{B_{i}\mid 0\leq i<n\}$, in which $x_{i}\in{{\mathbb{R}}}^{3}$ is the center and $r_{i}\in{{\mathbb{R}}}$ is the radius of $B_{i}$. The space-filling diagram is the union of these balls, $\bigcup X=\bigcup_{i=0}^{n-1}B_{i}$. We are interested in the volume, area, mean curvature, Gaussian curvature of this diagram. The volume should be clear. To describe the other three measures, we note that the boundary of $\bigcup X$ consists of patches of spheres that remain after we remove open caps, $S_{i}\setminus\bigcup_{j\neq i}{\rm int\,}{B_{j}}$, in which $S_{i}={\rm bd\,}{B_{i}}$. These patches meet along portions of circles that remain after we remove open arcs, $S_{ij}\setminus\bigcup_{k\neq i,j}{\rm int\,}{B_{k}}$, in which $S_{ij}=S_{i}\cap S_{j}$. These portions of circles consist of closed arcs that meet at corners of the form $S_{ijk}\setminus\bigcup_{\ell\neq i,j,k}{\rm int\,}{B_{\ell}}$, in which $S_{ijk}=S_{i}\cap S_{j}\cap S_{k}$. For ease of reference, we summarize and extend the introduced notation in Table 2, which is given in the appendix. The area is simply the sum of areas of the sphere patches. The mean curvature is more subtle because the boundary of the space-filling diagram is not a smooth surface. Instead of defining it as an integral of the point-wise mean curvature, we use a discrete formula according to which the contribution of an arc in the boundary of $\bigcup X$ is its length times the dihedral angle between the outward normals of the two spheres that share the arc; see (7) for the formula in the weighted case. Similarly, we write the Gaussian curvature as a sum of contributions of sphere patches, circular arcs, and corners; see (8) for the formula in the weighted case. Alternatively, we can make use of the Gauss–Bonnet theorem according to which the Gaussian curvature is $2\pi$ times the Euler characteristic of the surface, but no such shortcut is available for the weighted case, which we discuss next.

Voronoi decomposition and dual alpha complex. To generalize the measures to the weighted case, we introduce the Voronoi decomposition of the space-filling diagram and its dual, the alpha complex, which serves as a book-keeping device.

Recall that $X$ is a set of $n$ closed balls $B_{i}$ with centers $x_{i}$ and radii $r_{i}$. The power distance of a point $a\in{{\mathbb{R}}}^{3}$ from $B_{i}$ is ${\pi_{i}}(a)={\|{a}-{x_{i}}\|}^{2}-r_{i}^{2}$, and we note that $B_{i}={\pi_{i}}^{-1}(-\infty,0]$. The Voronoi domain of $B_{i}$ is ${{V}_{i}}=\{a\in{{\mathbb{R}}}^{3}\mid{\pi_{i}}(a)\leq{\pi_{j}}(a)\mbox{\rm\leavevmode\nobreak\ for\leavevmode\nobreak\ }0\leq j<n\}$, and we write ${{V}_{ij}}$, ${{V}_{ijk}}$, and ${{V}_{ijk\ell}}$ for the pair-wise, triple-wise, and quadruple-wise intersections. The Voronoi domains decompose the space-filling diagram into convex sets of the type $B_{i}\cap{{V}_{i}}$. We will need notation for the fraction of the ball represented by this set, and similarly for the fractions of common intersections, which we collect in Table 2, which is given in the appendix. Letting ${{w_{i}}}\in{{\mathbb{R}}}$ be the weight of the $i$-th closed ball, the formulas for the weighted measures of a space-filling diagram are given in (5), (6), (7), and (8) below.

The fractional measures listed in Table 2 are conveniently computed using the alpha complex of $X$, which is the nerve of the sets $B_{i}\cap{{V}_{i}}$ [9]. Recall that the Delaunay mosaic is the nerve of the sets ${{V}_{i}}$, which implies that the alpha complex is a subcomplex of the Delaunay mosaic, and generically, both are simplicial complexes realized in ${{\mathbb{R}}}^{3}$. It will be convenient to denote a simplex by the set of indices of its vertices, or by the sequence if we need an ordered version of the simplex. By alpha shape we mean the underlying space of the alpha complex. To explain the connection between the simplices and the measures, we note that a tetrahedron $ijk\ell$ in the Delaunay mosaic belongs to the alpha complex iff $\nu_{ijk\ell}>0$, and the same rule applies to the triangles, edges, and vertices of the mosaic. Furthermore, a triangle $ijk$ in the alpha complex belongs to the boundary of the alpha shape iff $\sigma_{ijk}>0$, and the same rule applies to the edges and vertices of the complex. To compute these fractions, we use inclusion-exclusion over subsets of simplices in the alpha complex. For example,

[TABLE]

in which the sums are over all edges, triangles, tetrahedra in the alpha complex that share $x_{i}$ as a vertex. Proofs and additional formulas can be found in [6].

Weighted intrinsic volume. Recall that $\bigcup X$ is the union of the balls $B_{i}$, for $0\leq i<n$. Its state, ${\bf x}\in{{\mathbb{R}}}^{3n}$, is the concatenation of the center vectors. In other words, the $(3i+\ell)$-th coordinate of ${\bf x}$ is the $\ell$-th coordinate of $x_{i}$. We use the Voronoi tessellation to divide the space-filling diagram into individual contributions and get the weighted intrinsic volume by multiplication with the corresponding weight.

Proposition 2.1** (Weighted Intrinsic Volumes).**

The weighted volume, surface area, mean curvature, and Gaussian curvature of the space-filling diagram, $\bigcup X$, are

[TABLE]

To get the above formula for the mean curvature, we smoothen the crevices of the space-filling diagram by rolling a ball of radius ${\varepsilon}>0$ about the surface. For sufficiently small ${\varepsilon}$, the resulting surface is everywhere differentiable, and we take the limit of its total mean curvature as ${\varepsilon}$ goes to zero. Along circular arcs, we partition the mean curvature in equal parts to the two intersecting spheres. The rationale for this rule is that the surface swept out by the centers of the rolling balls of different radii partitions the exterior dihedral angle at the arc in equal parts. Using the same geometric intuition, we partition the contribution of a corner of the space-filling diagram to the Gaussian curvature in proportions $\alpha_{i}+\alpha_{j}+\alpha_{k}=1$; see [2] for details. All variables have been defined above, except for $\phi_{ij}$ (the angle between the unit normals of the spheres at a point of $S_{ij}=S_{i}\cap S_{j}$), $\lambda_{ij}$ (the combined length of the two normals after projection to the line passing through $x_{i}$ and $x_{j}$), and $\phi_{ijk}$ (the solid angle spanned by the unit normals of $S_{i}$, $S_{j}$, $S_{k}$ at a point of common intersection, $P\in S_{ijk}$).

Weighted volume and area derivatives. To state the derivative of the weighted volume, we introduce the momentum, ${\bf t}\in{{\mathbb{R}}}^{3n}$, which is the concatenation of the velocity vectors. In other words, the vector ${\bf t}_{i}=[t_{3i+1},t_{3i+2},t_{3i+3}]^{T}$ is the velocity vector of the $i$-th ball. The derivative of a function $F\colon{{\mathbb{R}}}^{3n}\to{{\mathbb{R}}}$ can be given in terms of the gradient: ${\rm D}F_{\bf x}({\bf t})={\langle\nabla_{\!{\bf x}}f,{\bf t}\rangle}$. Writing {\bf v}=\nabla_{\!{\bf x}}and$$\bf vi = [v_3i+1, v_3i+2, v_3i+3]^T,\cite[cite]{[\@@bibref{}{EdKo03}{}{}]}givesthederivativeoftheweightedvolume:\begin{proposition}[Gradient of Weighted Volume]Thederivativeoftheweightedvolumeofthespace-fillingdiagramatstate{\bf x}$withmomentum${\bf t}$is${\rm D}{\bf x}({\bf t})={\langle{\bf v},{\bf t}\rangle}$with\begin{aligned} {\bf v}_{i}&=\sum_{j}\left[(v_{ij}+v_{ji}){\bf u}_{ij}+(\bar{v}_{ij}-\bar{v}_{ji}){\bf U}_{ij}\right],\\ v_{ij}&=\tfrac{1}{2}\pi{{w_{i}}}r_{ij}^{2}\nu_{ij}\left(1+\frac{r_{i}^{2}-r_{j}^{2}}{{\|{x_{i}}-{x_{j}}\|}^{2}}\right),\\ \bar{v}_{ij}&=\tfrac{2}{3}\pi{{w_{i}}}r_{ij}^{2}\nu_{ij}\frac{1}{{\|{x_{i}}-{x_{j}}\|}},\end{aligned}inwhich${\bf U}{ij}$isthevectorfrom$x_{ij}$tothecentroidof$B_{ij}\cap{{V}{ij}}$,andthesumisoveralledgesinthealphacomplexthatareincidentto$x{i}.\end{proposition}Writing$$\bf a= ∇​$\bf x$ \it area$$and$$\bf ai = [a_3i+1, a_3i+2, a_3i+3]^T,\cite[cite]{[\@@bibref{}{BEKL04}{}{}]}givesthederivativeoftheweightedarea:\begin{proposition}[Gradient of Weighted Area]Thederivativeoftheweightedareaofthespace-fillingdiagramatstate{\bf x}$withmomentum${\bf t}$is${\rm D}{\it area}{\bf x}({\bf t})={\langle{\bf a},{\bf t}\rangle}$with\begin{aligned} {\bf a}_{i}&=\sum_{j}(a_{ij}+a_{ji}){\bf u}_{ij}+\sum_{j,k}(a_{ijk}-a_{jik}){\bf u}_{ijk},\\ a_{ij}&=\pi{{w_{i}}}\sigma_{ij}r_{i}\left(1-\frac{r_{i}^{2}-r_{j}^{2}}{{\|{x_{i}}-{x_{j}}\|}^{2}}\right),\\ a_{ijk}&=2{{w_{i}}}r_{ijk}\nu_{ijk}r_{i}\frac{1}{{\|{x_{i}}-{x_{j}}\|}},\end{aligned}inwhichthefirstsumisovertheboundaryedgesinthealphacomplexincidentto$x{i},andthesecondsumisoverthetrianglesincidenttotheseedges.\end{proposition}\par\par\par

3 Derivatives

We recall (7) and decompose the weighted mean curvature function into a first term, which accounts for the curvature within the sphere patches, and a second term, which accounts for the curvature concentrated along the circular arcs:

[TABLE]

While not reflected in the notation, $\sigma_{i}$, $r_{ij}$, $\phi_{ij}$, and $\sigma_{ij}$ depend on the state. It is convenient to write $\sigma_{i}(\tau)$ for $\sigma_{i}$ at state ${\bf x}+\tau{\bf t}$, and similarly for the other functions. Furthermore, we write $\sigma_{i}^{\prime}$ for the derivative of $\sigma_{i}$ at $\tau=0$, and similarly we write $r_{ij}^{\prime}$, $\phi_{ij}^{\prime}$, $\sigma_{ij}^{\prime}$. Derivatives with respect to parameters other than $\tau$ are explicitly stated as such.

3.1 Derivative of $\sigma_{i}$

To derive $\sigma_{i}$, we follow [4] and decompose the motion into a direction preserving component and a distance preserving component. The direction preserving motion stretches the distance between two centers, while the distance preserving motion rotates one center about the other. We therefore write ${\bf t}_{j}-{\bf t}_{i}={\bf t}^{\rm str}_{j}+{\bf t}^{\rm rot}_{j}$, in which ${\bf t}^{\rm str}_{j}={\langle{\bf t}_{j}-{\bf t}_{i},{\bf u}_{ij}\rangle}{\bf u}_{ij}$. Observe that ${\bf t}^{\rm rot}_{j}=({\bf t}_{j}-{\bf t}_{i})-{\bf t}^{\rm str}_{j}$ is the velocity vector of the rotation with angular momentum

[TABLE]

in which $\times$ is the cross or outer product that maps vectors $a$ and $b$ to the vector $c=a\times b$ of length ${\|{a}\|}{\|{b}\|}\sin\angle ab$ normal to both in such a way that the sequence $a,b,c$ forms a right-handed system in space; see Figure 1.

The analysis of the direction preserving motion is similar to the analysis of the area derivative in [4], and we repeat the main steps for completeness. Since little in terms of a proof of the distance preserving motion is offered in [4], we give a complete argument, which can be translated back almost verbatim to a proof in the area case.

Direction preserving motions. Consider two balls, $B_{i}$ and $B_{j}$, and the change of surface area under a motion of $x_{i}$ that preserves the direction, ${\bf u}_{ij}$. Let ${\|{x_{i}}-{x_{j}}\|}=\xi_{i}+\xi_{j}$, in which the terms on the right-hand side of the equation are the signed distances from the centers to the plane bisector. As noted in [4], we have

[TABLE]

and similarly for $\xi_{j}$. Plugging (11) into Formula 2, we get

[TABLE]

Differentiating with respect to the distance between the two centers, we get the rate of area change, which happens right next to the circle of intersection. In the general case, the cap $B_{j}\cap S_{i}$ may overlap with other caps so that only a fraction of the bounding circle belongs to the boundary of the space-filling diagram. Accordingly, the derivative is the same fraction of the derivative without any such overlap:

[TABLE]

Multiplying the right-hand side of (13) with ${\langle{\bf t}_{j}-{\bf t}_{i},{\bf u}_{ji}\rangle}$, we get the derivative with respect to the direction preserving motion of $x_{j}$, and taking the sum over all $j\neq i$, we get the derivative of $\sigma_{i}$ with respect to the direction preserving motions of all $x_{j}$; compare with the first sum on the right-hand side of (15) in Lemma 3.1.

Distance preserving motions. Consider first the case of a single cap, $B_{j}\cap S_{i}$. The rotation defined by ${\boldsymbol{\omega}}_{j}$ keeps the cap as well as $\sigma_{i}$ at constant size. This is because $S_{i}$ loses area along the front of the circle that bounds the cap, it gains area along the back of this circle, and the loss equals the gain. In the more general case depicted in Figure 2, it is possible that the loss and the gain no longer cancel each other.

To calculate the net loss or gain, we let $N$ and $S$ be antipodal points on $S_{i}$, which we refer to as north- and south-pole, such that the cap center lies on the equator — the great-circle halfway between $N$ and $S$ — and the motion follows the equator. We recall that ${\bf t}^{\rm rot}_{j}$ is the velocity vector of $x_{j}$, and $r_{i}/{\|{x_{i}}-{x_{j}}\|}$ times ${\bf t}^{\rm rot}_{j}$ is the velocity vector of the cap center. In Figure 2, we draw $N\!S$ vertical, we place $x_{j}$ and the cap center on top of $x_{i}$, and we let the velocity vector go horizontally from right to left. By Formula 2, the area loss along an arc of the front does not depend on its position but only on the length of its projection onto $N\!S$. We are interested in arcs of $S_{ij}$ that lie on the boundary of the space-filling diagram, and we compute the net loss or gain indirectly, by measuring the projections of the edges in the Voronoi tessellation defined by $B_{i}$, $B_{j}$, and a third ball, $B_{k}$. Recall that $B_{ijk}$ is the line segment that connects the two points in which $S_{i},S_{j},S_{k}$ meet, and that $\nu_{ijk}$ is the fraction of this line segment that belongs to the Voronoi decomposition of the space-filling diagram. If this edge belongs to the front of the clipped cap, then its contribution is positive because it subtracts from the loss along the front. Symmetrically, if the edge belongs to the back of the clipped cap, then its contribution is negative. To quantify, we recall the definition of ${\bf u}_{ijk}$ given in Table 2. The signed fraction of the projection of the relevant portion of $B_{ijk}$ onto $N\!S$ is its length, $2r_{ijk}\nu_{ijk}$, times the cosine of the projection angle, ${\langle{\bf t}^{\rm rot}_{j},{\bf u}_{ijk}\rangle}/{\|{{\bf t}^{\rm rot}_{j}}\|}$, divided by the diameter of the sphere, $2r_{i}$. By Formula 2, this is also the area fraction of the corresponding strip around the sphere. Dividing by the length of the circular orbit, $2\pi{\|{x_{i}}-{x_{j}}\|}$, we get the derivative:

[TABLE]

in which we write $\angle x_{i}x_{j}$ for the angle that parametrizes the rotation of $x_{j}$. We get the contribution to $\sigma_{i}^{\prime}$ by multiplying with ${\|{{\bf t}^{\rm rot}_{j}}\|}$ and observing that ${\langle{\bf t}^{\rm rot}_{j},{\bf u}_{ijk}\rangle}={\langle{\bf t}_{j}-{\bf t}_{i},{\bf u}_{ijk}\rangle}$ for all $k$; compare with the second sum on the right-hand side of (15) in Lemma 3.1.

Summary. A clipped cap whose entire boundary consists of straight line segments does not touch the boundary of the space-filling diagram and therefore has no contribution to the derivative. These caps correspond to interior edges of the alpha complex. Omitting them from the sum, we get the derivative of $\sigma_{i}$.

Lemma 3.1** (Derivative of $\sigma_{i}$).**

The derivative of the fraction of $S_{i}$ on the boundary of the space-filling diagram at state ${\bf x}$ with momentum ${\bf t}$ is

[TABLE]

in which the first sum is over the boundary edges of the alpha complex incident to $x_{i}$, with coefficient given in (13), and the second sum is over the triangles incident to these edges.

3.2 Derivatives of $r_{ij}$ and of $\phi_{ij}$

Recall that $r_{ij}$ is the radius of the circle $S_{ij}=S_{i}\cap S_{j}$, assuming the two spheres have a non-empty intersection. As illustrated in Figure 3, it is also the height of the triangle with base of length ${\|{x_{i}}-{x_{j}}\|}$ and edges of length $r_{i}$ and $r_{j}$.

The area of the triangle is $A=\tfrac{1}{2}r_{ij}{\|{x_{i}}-{x_{j}}\|}$. Alternatively, we can compute the area using the version of the Heron formula given right after Formula 2:

[TABLE]

The derivative of the area with respect to the distance between the two centers is

[TABLE]

From this, we get the derivative of the height of the triangle:

[TABLE]

Switching our attention, we observe that the angle between the outward normals, $\phi_{ij}$, is also the angle opposite the base inside the triangle in Figure 3. Recall that the Law of Cosines generalizes Pythagoras’ Theorem beyond right-angled triangles: $c^{2}=a^{2}+b^{2}-2ab\cos\gamma$, in which $a,b,c$ are the lengths of the sides and $\gamma$ is the angle opposite to the side of length $c$. In our application, we have $a=r_{i}$, $b=r_{j}$, $c={\|{x_{i}}-{x_{j}}\|}$, and $\gamma=\phi_{ij}$. Hence,

[TABLE]

Recall that the derivative of $\arccos x$ is $-1/\sqrt{1-x^{2}}$. Together with the chain rule for differentiation, this gives the derivative of the angle with respect to the length of the base:

[TABLE]

We summarize the results of this subsection for later reference.

Lemma 3.2** (Derivatives of $r_{ij}$ and of $\phi_{ij}$).**

The derivatives of the radius of the intersection circle of two spheres and of the angle between the outward normals at the intersection at state ${\bf x}$ with momentum ${\bf t}$ are

[TABLE]

with the coefficients in the two equations given in (21) and (24).

3.3 Derivative of $\sigma_{ij}$

Recall that $\sigma_{ij}$ is the fraction of the circle $S_{ij}=S_{i}\cap S_{j}$ that belongs to the space-filling diagram. We compute its derivative under the motion ${\bf t}$ in several steps, the first of which modifies the motion. Without altering the derivative, we do this such that the center and the plane of the circle are fixed.

Modifying the motion. We begin by fixing $x_{i}$ in space, which we do by changing the velocity vector of $x_{j}$ to ${\bf t}_{j}-{\bf t}_{i}$ for every $j$, as before. Next, we fix the normal direction of the plane that contains $S_{ij}$ by removing the angular momentum. Specifically, we change the velocity vector of $x_{j}$ to ${\bf V}\!_{ij}={\bf t}_{j}-{\bf t}_{i}-{\boldsymbol{\omega}}_{j}\times(x_{j}-x_{i})$, in which ${\boldsymbol{\omega}}_{j}$ is given in (10); see Figure 1. After this modification, the point $x_{j}$ moves with speed $v$ along ${\bf u}_{ji}$, in which $v$ is the length of the new velocity vector. Accordingly, we change the velocity vector of $x_{k}$ to ${\bf V}\!_{ijk}={\bf t}_{k}-{\bf t}_{i}-{\boldsymbol{\omega}}_{j}\times(x_{k}-x_{i})$ for every $k$, and note that ${\bf V}\!_{iji}=0$ and ${\bf V}\!_{ijj}={\bf V}\!_{ij}$. We finally fix the plane of the circle. To this end, we recall that $\xi_{i}$ and $\xi_{j}$ given in (11) are the signed distances of $x_{i}$ and $x_{j}$ from the plane of $S_{ij}$. We have $\xi_{i}+\xi_{j}={\|{x_{i}}-{x_{j}}\|}$ and use them to write the center of the circle as an affine combination of the centers of the spheres:

[TABLE]

Suppose now that $x_{i}$ is fixed and $x_{j}$ moves with speed $v$ in the direction of ${\bf u}_{ji}$. To compute the speed of $x_{ij}$ moving in the same direction, we set $x_{i}=0$, write $x_{j}={\|{x_{i}}-{x_{j}}\|}$, and simplify (27) to $x_{ij}=\tfrac{1}{2}({\|{x_{i}}-{x_{j}}\|}+(r_{i}^{2}-r_{j}^{2})/{\|{x_{i}}-{x_{j}}\|})$. Its derivative with respect to the distance is

[TABLE]

we write $D$ for the derivative in this special, $1$-dimensional scenario, and note that $x_{ij}$ moves with speed $Dv$. Subtracting the corresponding multiple of ${\bf u}_{ji}$ from the velocity vector of every point $x_{k}$, we get the final collection of vectors.

Lemma 3.3** (Change of Motion).**

Replacing the velocity vector ${\bf t}_{k}$ by ${\bf T}_{ijk}={\bf V}\!_{ijk}-D{\bf V}\!_{ij}$, for every $k$, fixes the center and the plane of the circle $S_{ij}=S_{i}\cap S_{j}$ while preserving the derivative of $\sigma_{ij}$. The coefficient of ${\bf V}\!_{ij}$ is given in (28).

Movement of $x_{k}$. The circle $S_{ij}$ alternates between the boundary and the interior of the space-filling diagram, and we are interested in the arcs that belong to the boundary. The endpoints of these arcs belong to [math]-spheres of the form $S_{ijk}=S_{i}\cap S_{j}\cap S_{k}$. Recall that $\sigma_{ijk}$ is [math], $\tfrac{1}{2}$, or $1$. If $\sigma_{ijk}=1$, then both points of $S_{ijk}$ lie on the boundary of the space-filling diagram and therefore serve as endpoints of the relevant arcs of $S_{ij}$, if $\sigma_{ijk}=\tfrac{1}{2}$, then only one of the two points serves in this capacity, and if $\sigma_{ijk}=0$, then neither point serves in this capacity. Consider a sphere, $S_{k}$, such that $\sigma_{ijk}=\tfrac{1}{2}$ or $1$, and let $P\in S_{ijk}$ be an endpoint of an arc; see Figure 4.

We are interested in the speed in which the point $P$ moves along $S_{ij}$. To compute this speed, we let ${\bf u}_{ij}^{P}$ be the unit vector that is tangent to $S_{ij}$ at the point $P$ such that $x_{k}$ and $P+{\bf u}_{ij}^{P}$ lie on the same side of the tangent plane, that is:

[TABLE]

and the correct sign is the one for which ${\langle x_{k}-P,{\bf u}_{ij}^{P}\rangle}>0$; see Figure 4. When $x_{k}$ moves along ${\bf T}_{ijk}$ with speed ${\|{{\bf T}_{ijk}}\|}$, then the tangent plane moves along $x_{k}-P$ with speed $v_{k}={\langle x_{k}-P,{\bf T}_{ijk}\rangle}/{\|{x_{k}}-{P}\|}$. At the same time, the tangent plane moves along ${\bf u}_{ij}^{P}$ with a speed $v$ such that substituting $v{\bf u}_{ij}^{P}$ for ${\bf T}_{ijk}$ in this equation gives the same speed, namely $v_{k}$. In other words, $v{\langle x_{k}-P,{\bf u}_{ij}^{P}\rangle}={\langle x_{k}-P,{\bf T}_{ijk}\rangle}$ and therefore $v={\langle x_{k}-P,{\bf T}_{ijk}\rangle}/{\langle x_{k}-P,{\bf u}_{ij}^{P}\rangle}$. This implies

[TABLE]

in which $\angle x_{i}x_{j}x_{k}$ is the angle parametrizing the circular motion of $x_{k}$ about the line passing through $x_{i}$ and $x_{j}$, with $k$ is such that $P\in S_{ijk}$; see the first sum on the right-hand side of (25) in Lemma 3.2.

Change of $r_{ij}$. Recall that Lemma 3.2 gives the rate of change of the radius. We ask how it affects the derivative of $\sigma_{ij}$. To begin, we observe that $h_{k}={\langle x_{ij}-P,x_{k}-P\rangle}/{\|{x_{k}}-{P}\|}$ is the signed distance of $x_{ij}$ from the tangent plane. It is positive if $x_{ij}$ and $x_{k}$ lie on the same side of the plane, and negative if they lie on opposite sides. Letting ${\alpha}^{P}$ be the angle between $x_{ij}-P$ and the line in which the tangent plane intersects the plane of $S_{ij}$, we get $h_{k}=r_{ij}\sin{\alpha}^{P}$. Hence, ${\alpha}^{P}=\arcsin(h_{k}/r_{ij})$, and we are interested in the derivative with respect to the radius. Remembering that the derivative of $\arcsin x$ is $1/\sqrt{1-x^{2}}$, we get

[TABLE]

Using the chain rule, we get the derivative with respect to ${\|{x_{i}}-{x_{j}}\|}$ by multiplying (31) with (21) and then taking the sum over all points $P$; see the second sum on the right-hand side of (25) in Lemma 3.2.

Summary. Adding the contribution of the changing radius to those of the movements of the $x_{k}$ and normalizing by the length of the circle, we get the desired derivative.

Lemma 3.4** (Derivative of $\sigma_{ij}$).**

The derivative of the fraction of $S_{ij}$ on the boundary of the space-filling diagram at state ${\bf x}$ with momentum ${\bf t}$ is

[TABLE]

in which the coefficients in the second sum are given in (31) and (21). Both sums are over all oriented boundary triangles of the alpha shape that share the edge from $x_{i}$ to $x_{j}$, and $P\in S_{ijk}$ is the corresponding corner of the space-filling diagram.

4 Gradients

We write the derivative of the weighted mean curvature function, ${\it mean}\colon{{\mathbb{R}}}^{3n}\to{{\mathbb{R}}}$, in terms of the gradient of ${\it mean}$ at ${\bf x}\in{{\mathbb{R}}}^{3n}$, denoted ${\bf m}=\nabla_{\!{\bf x}}{\it mean}$. Recalling (9), this derivative is ${\it mean}^{\prime}=p^{\prime}+q^{\prime}+s^{\prime}$, with

[TABLE]

Writing ${\bf m}=[m_{1},m_{2},\ldots,m_{3n}]^{T}$, we recall that ${\bf m}_{i}=[m_{3i+1},m_{3i+2},m_{3i+3}]^{T}$ is the $3$-dimensional gradient that applies to $x_{i}$. Using boldface letters for the gradients of $p,q,s$, and similar conventions for the $3$-dimensional sub-vectors, we have ${\bf m}={\bf p}+{\bf q}+{\bf s}$ and ${\bf m}_{i}={\bf p}_{i}+{\bf q}_{i}+{\bf s}_{i}$ for $0\leq i<n$. We get the gradients by redistributing the derivatives stated in Lemmas 3.1 to 3.4.

First term. To begin, we use Lemma 3.1 to rewrite (33) as

[TABLE]

in which the first sum is over all directed boundary edges of the alpha shape, and the second sum is over all triangles incident to these edges. Observe that for fixed $i$, we get possibly non-zero contributions to all ${\bf p}_{j}$. Symmetrically, we get ${\bf p}_{i}$ by accumulating contributions from all $j$. Using ${\bf u}_{ji}=-{\bf u}_{ij}$ and ${\bf u}_{jik}={\bf u}_{ijk}$, we get

[TABLE]

in which the first sum is over all boundary edges of the alpha shape incident to $x_{i}$, and the second sum is over all triangles incident to these edges. Not surprisingly, the result is similar to the weighted area gradient given in Proposition 2.3. Specifically, we get $a_{ij}=p_{ij}r_{i}$ and $a_{ijk}=p_{ijk}r_{i}$.

Second term. We use Lemma 3.2 to rewrite (34) as

[TABLE]

in which the sum is over all directed boundary edges of the alpha shape. Redistributing the terms, we get

[TABLE]

in which the sum is over all boundary edges of the alpha shape incident to $x_{i}$.

Third term. The redistribution of the pieces on the right-hand side of (32) to get the gradient is complicated by the change of the motion. We therefore begin by rewriting ${\bf T}_{ijk}$. To this end, we recall that the cross product is distributive: $a\times b+a\times c=a\times(b+c)$, and that the Lagrange formula turns a triple cross product into two scalar products: $(a\times b)\times c=b{\langle a,c\rangle}-a{\langle b,c\rangle}$. Recalling the definition of $D$ from (28), we can now rewrite the new motion vector:

[TABLE]

in which the third line is obtained by applying the Lagrange formula and writing ${\bf d}_{ijk}=(x_{k}-Dx_{j}+(D-1)x_{i})/{\|{x_{i}}-{x_{j}}\|}$. We need the scalar product of ${\bf T}_{ijk}$ with $x_{k}-P$:

[TABLE]

in which

[TABLE]

We are now ready to rewrite (35) using Lemma 3.4 as

[TABLE]

in which the sum is over the three ordered versions of all oriented boundary triangles of the alpha shape. Letting $i,j,k$ be the vertices of such a triangle, the ordered versions are given by the index triplets $ijk,kij,jki$, and we write $P$ for the corresponding corner of the space-filling diagram, noting that $P\in S_{ijk}$. Observe that an unordered triangle that belongs to $0,1,2$ tetrahedra in the alpha complex occurs $2,1,0$ times in this sum. Redistributing the terms, we get

[TABLE]

in which the sum is over all oriented boundary triangles of the alpha shape that share $x_{i}$. As before, $P\in S_{ijk}$ is the corresponding corner of the space-filling diagram.

Summary. We finally get the gradient of the weighted mean curvature function by adding the gradients of the three component functions:

Theorem 4.1** (Gradient of Weighted Mean Curvature).**

The derivative of the weighted mean curvature of the space-filling diagram at state ${\bf x}$ with momentum ${\bf t}$ is ${\rm D}{\it mean}_{\bf x}({\bf t})={\langle{\bf m},{\bf t}\rangle}$, in which ${\bf m}_{i}={\bf p}_{i}+{\bf q}_{i}+{\bf s}_{i}$ as given in (39), (44), and (58), for all $0\leq i<n$.

5 Violations of Continuity

To embed our formulas in the inner loop of a molecular dynamics application, the implementation must be efficient and robust. We address the latter requirement by identifying the subset of ${{\mathbb{R}}}^{3n}$ where the gradient of the weighted mean curvature function is not continuous. This subset is contained in the subset of non-generic states, which we describe first.

General position. We distinguish two types of degenerate states: where the Delaunay mosaic is ambiguous, and where the Delaunay mosaic is unambiguous but the alpha complex is ambiguous. Recall that $X$ is a collection of $n$ closed balls in ${{\mathbb{R}}}^{3}$, and that these balls can move individually but their radii are fixed. We say $X$ is in general position and, equivalently, that its state is generic if the following three conditions hold:

1. I.

the common intersection of $p+1$ Voronoi domains is either empty or a convex polyhedron of dimension $3-p$; 2. 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 any five Voronoi domains have an empty common intersection, and Condition II implies that any four spheres have an empty common intersection. Each violation of Condition I corresponds to a $(3n-1)$-dimensional submanifold of ${{\mathbb{R}}}^{3n}$ and so does every violation of Condition II, except when two radii are the same, in which case we get a submanifold of dimension $3n-3$. Let ${{\mathbb{M}}_{\rm{I}}}$ and ${{\mathbb{M}}_{\rm{II}}}$ be the union of submanifolds that correspond to violations of Conditions I and II, respectively. Since there are only finitely many such submanifolds, ${{\mathbb{M}}_{\rm{I}}}$ and ${{\mathbb{M}}_{\rm{II}}}$ have dimension $3n-1$ each.

Condition I and flips. If Condition I is satisfied, then the Delaunay mosaic is simplicial. We get a violation when the state trajectory intersects ${{\mathbb{M}}_{\rm{I}}}$. We limit ourselves to discussing what we call the typical case, in which there is only one violation of general position. Indeed, multiple violations correspond to intersections of the $(3n-1)$-dimensional submanifolds. In the typical case, the intersection is an isolated point where the trajectory passes locally from one side of ${{\mathbb{M}}_{\rm{I}}}$ to the other. The corresponding change in the Delaunay mosaic is a flip, of which there are four types. For integers $1\leq b\leq 4$ and $a=5-b$, the $b$-to-$a$ flip replaces $b$ tetrahedra by $a$ tetrahedra; see [7] for details on this operation. The Delaunay mosaic is a simplicial complex geometrically realized in ${{\mathbb{R}}}^{3}$, both before and after the replacement. This implies that the union of the $b$ tetrahedra before the flip be the same as the union of the $a$ tetrahedra after the flip.

To be more formal, let $B$ be the complex consisting of the $b$ tetrahedra and their faces, and let $A$ be the complex consisting of the $a$ tetrahedra and their faces. Then $C=B\cap A$ is the common boundary of both complexes, and the flip substitutes $A\setminus C$ for $B\setminus C$ in the Delaunay mosaic. By construction, either all simplices of $B\setminus C$ and of $A\setminus C$ belong to the alpha complex before and after the flip, or none does. In the latter case, none of the formulas in Theorem 4.1 are affected by the flip. In the former case, some terms in the formulas get replaced by other terms. At the moment of the flip, the replaced and the replacing terms evaluate to the same numerical value, which implies that the weighted mean curvature and its derivative are unchanged and therefore continuous. We can see this by noting that the flip does not alter the combinatorial structure of the space-filling diagram boundary and therefore affects neither the weighted mean curvature nor its derivative.

Condition II and critical simplices. Condition II addresses situations in which the alpha complex changes while the Delaunay mosaic stays the same: we either add one or more simplices in the mosaic to the complex or we remove them from the complex. Here we consider the case in which only one simplex is added or removed. With reference to the discrete Morse theory of Delaunay mosaics outlined in [3], we call this a critical simplex. We distinguish three cases.

Case** ${{\bf{C}}_{\bf{1}}}$:**

the added or removed simplex is an edge of the Delaunay mosaic. At the moment of the change, the spheres whose centers are the endpoints of the edge touch in a single point. Letting $i$ and $j$ be the indices of the spheres, ${\varepsilon}=r_{i}+r_{j}-{\|{x_{i}}-{x_{j}}\|}$ is negative when the spheres are disjoint and positive when they intersect in a circle. For ${\varepsilon}<0$, the edge does not belong to the alpha complex and therefore contributes neither to the weighted mean curvature nor to its gradient. For ${\varepsilon}>0$, the change to the weighted mean curvature caused by the edge depends on the change in area of the space-filling diagram, the change in length of its arcs, and the dihedral angles at the arcs. Focussing on the rough order of the changes, it is not difficult to see that $\Delta{\it area}={\varepsilon}$, $\Delta{\it length}=\sqrt{{\varepsilon}}$, ${\it Angle}=1$, and therefore $\Delta{\it mean}=\sqrt{{\varepsilon}}$ and ${\|{\nabla{\it mean}}\|}=\infty$ at ${\varepsilon}=0$.

The other two cases are similar so we will be brief. In each case, we use ${\varepsilon}$ to parameterize the local motion, with ${\varepsilon}=0$ at the moment of change. We feel free to make further non-essential simplifying assumptions, such as the radii being distinct.

Case** ${{\bf{C}}_{\bf{2}}}$:**

the added or removed simplex is a triangle of the Delaunay mosaic. At the moment of change, the spheres whose centers are the vertices of the triangle meet in a single point. By easy analysis, we get $\Delta{\it area}={\varepsilon}\sqrt{{\varepsilon}}$, $\Delta{\it length}=\sqrt{{\varepsilon}}$, ${\it Angle}=1$, and therefore $\Delta{\it mean}=\sqrt{{\varepsilon}}$ and ${\|{\nabla{\it mean}}\|}=\infty$ at ${\varepsilon}=0$.

Case** ${{\bf{C}}_{\bf{3}}}$:**

the added or removed simplex is a tetrahedron of the Delaunay mosaic. At the moment of change, the spheres whose centers are the vertices of the tetrahedron meet in a single point. By easy analysis, we get $\Delta{\it area}={\varepsilon}^{2}$, $\Delta{\it length}={\varepsilon}$, ${\it Angle}=1$, and therefore $\Delta{\it mean}={\varepsilon}$ and $\Delta{\|{\nabla{\it mean}}\|}=1$ at ${\varepsilon}=0$.

In conclusion, the weighted mean curvature is continuous but its gradient occasionally changes discontinuously; see Table 1 for a convenient summary of the three cases above as well as the six cases to be discussed shortly.

Condition II and non-singular intervals. Here we consider the case in which at least two simplices are added or removed at the same moment. Generically, these simplices form an interval of size $2$, $4$, or $8$; see [3] for details. We encounter six cases.

Case** ${{\bf{N}}_{\bf{01}}}$:**

A vertex together with an incident edge are added or removed. At the moment of change, the sphere centered at the vertex breaks through the surface of the sphere whose center is the other endpoint of the edge. By easy analysis, we get $\Delta{\it area}={\varepsilon}$, $\Delta{\it length}=\sqrt{{\varepsilon}}$, and because $\phi_{ij}$ is proportional to $r_{ij}$ also ${\it Angle}=\sqrt{{\varepsilon}}$. Hence $\Delta{\it mean}={\varepsilon}$ and $\Delta{\|{\nabla{\it mean}}\|}=1$ at ${\varepsilon}=0$.

Case** ${{\bf{N}}_{\bf{02}}}$:**

A vertex together with an incident triangle and its two edges that share the vertex are added or removed. At the moment of change, the sphere centered at the vertex breaks through the circle at which the two spheres centered at the other vertices of the triangle intersect. By easy analysis, we get $\Delta{\it area}={\varepsilon}\sqrt{{\varepsilon}}$, $\Delta{\it length}=\sqrt{{\varepsilon}}$, ${\it Angle}=1$, and therefore $\Delta{\it mean}=\sqrt{{\varepsilon}}$ and ${\|{\nabla{\it mean}}\|}=\infty$ at ${\varepsilon}=0$.

Case** ${{\bf{N}}_{\bf{03}}}$:**

A vertex together with an incident tetrahedron and its three edges and three triangles that share the vertex are added or removed. At the moment of change, the sphere centered at the vertex breaks through the pair of points in which the three spheres centered at the other vertices of the tetrahedron intersect. We get $\Delta{\it area}={\varepsilon}^{2}$, $\Delta{\it length}={\varepsilon}$, ${\it Angle}=1$, and therefore $\Delta{\it mean}={\varepsilon}$ and $\Delta{\|{\nabla{\it mean}}\|}=1$ at ${\varepsilon}=0$.

Case** ${{\bf{N}}_{\bf{12}}}$:**

An edge together with an incident triangle are added or removed. At the moment of change, the circle in which the two spheres centered at the endpoints of the edge intersect breaks through the surface of the sphere centered at the third vertex of the triangle. By easy analysis, we get $\Delta{\it area}={\varepsilon}\sqrt{{\varepsilon}}$, $\Delta{\it length}=\sqrt{{\varepsilon}}$, ${\it Angle}=1$, and therefore $\Delta{\it mean}=\sqrt{{\varepsilon}}$ and ${\|{\nabla{\it mean}}\|}=\infty$ at ${\varepsilon}=0$.

Case** ${{\bf{N}}_{\bf{13}}}$:**

An edge together with an incident tetrahedron and its two triangle that share the edge are added or removed. At the moment of change, the circle in which the two spheres centered at the endpoints of the edge intersect breaks through the circle in which the two spheres centered at the other vertices of the tetrahedron intersect. We get $\Delta{\it area}={\varepsilon}^{2}$, $\Delta{\it length}={\varepsilon}$, ${\it Angle}=1$, and therefore $\Delta{\it mean}={\varepsilon}$ and $\Delta{\|{\nabla{\it mean}}\|}=1$ at ${\varepsilon}=0$.

Case** ${{\bf{N}}_{\bf{23}}}$:**

A triangle together with an incident tetrahedron are added or removed. At the moment of change, the pair of points in which the three spheres centered at the vertices of the triangle intersect break through the surface of the sphere centered at the fourth vertex of the tetrahedron. We get $\Delta{\it area}={\varepsilon}^{2}$, $\Delta{\it length}={\varepsilon}$, ${\it Angle}=1$, and therefore $\Delta{\it mean}={\varepsilon}$ and $\Delta{\|{\nabla{\it mean}}\|}=1$ at ${\varepsilon}=0$.

See again Table 1 for a convenient summary of the results in all cases. We note that the situation in the unweighted case has better continuity properties along ${{\mathbb{M}}_{\rm{II}}}\subseteq{{\mathbb{R}}}^{3n}$ in some cases, but the analysis is more involved. Subtracting ${{\mathbb{M}}_{\rm{I}}}\cup{{\mathbb{M}}_{\rm{II}}}$ from ${{\mathbb{R}}}^{3n}$, we are left with finitely many open cells such that the formula of the weighted mean curvature and of its derivative are both invariant over each cell. All terms in these formulas are continuous over the cell, which implies that both ${\it mean}\colon{{\mathbb{R}}}^{3n}\to{{\mathbb{R}}}$ and $\nabla{\it mean}\colon{{\mathbb{R}}}^{3n}\to{{\mathbb{R}}}^{3n}$ are continuous over the open cell. As argued above, ${\it mean}$ and $\nabla{\it mean}$ are also continuous at states ${\bf x}\in{{\mathbb{M}}_{\rm{I}}}\setminus{{\mathbb{M}}_{\rm{II}}}$, hence they are continuous at all ${\bf x}\in{{\mathbb{R}}}^{3n}\setminus{{\mathbb{M}}_{\rm{II}}}$. We summarize the findings of this section.

Theorem 5.1** (Continuity of Gradient).**

The gradient of the weighted mean curvature of a space-filling diagram of $n$ closed balls in ${{\mathbb{R}}}^{3}$ is continuous provided the state ${\bf x}\in{{\mathbb{R}}}^{3n}$ of the diagram does not belong to ${{\mathbb{M}}_{\rm{II}}}$, which is a $(3n-1)$-dimensional subset of ${{\mathbb{R}}}^{3n}$.

6 Discussion

The main contribution of this paper is the analysis of the derivative of the weighted mean curvature of the space-filling diagram of a set of moving balls. Specifically, we give an explicit description of the gradient of the weighted mean curvature function, which for $n$ spheres is a map from ${{\mathbb{R}}}^{3n}$ to ${{\mathbb{R}}}$. In addition, we characterize the subset of ${{\mathbb{R}}}^{3n}$ at which the derivative violates continuity. In total, this is sufficient information for an efficient and robust implementation of the weighted mean curvature derivative, one that can be added to the inner loop of a molecular dynamics simulation of a physical system. There are several questions about this application that go beyond the scope of this paper:

• •

Is there a connection between the coefficients with which the morphological approach combines the intrinsic volumes [13, 18] and the states at which their gradients are not continuous?

• •

Splitting the mean curvature concentrated along an arc in equal parts is suggested by the Apollonius diagram of the spheres, but splitting it according to the Voronoi diagram is also feasible. Are there physical reasons to prefer one split over the other?

We remark that the formulas in this paper can be easily adapted to splitting the mean curvature according to the Voronoi tessellation of the spheres.

Acknowledgment

The authors of this paper thank Roland Roth for suggesting the analysis of the weighted curvature derivatives for the purpose of improving molecular dynamics simulations and for his continued encouragement. They also thank Patrice Koehl for the implementation of the formulas and for his encouragement and advise along the road. Finally, they thank two anonymous reviewers for their constructive criticism.

Appendix A Notation

Appendix B Definitions and Results

• •

Section 1: Introduction.

• •

Section 2: Background.

• –

Formula 2 (Archimedes’ Theorem).

• –

Formula 2 (Heron’s Theorem).

• –

Proposition 2.1 (Intrinsic Volumes).

• –

Proposition 2.2 (Gradient of Weighted Volume).

• –

Proposition 2.3 (Gradient of Weighted Area).

• •

Section 3: Derivatives.

• –

Lemma 3.1 (Derivative of $\sigma_{i}$).

• –

Lemma 3.2 (Derivatives of $r_{ij}$ and of $\phi_{ij}$).

• –

Lemma 3.3 (Change of Motion).

• –

Lemma 3.4 (Derivative of $\sigma_{ij}$).

• •

Section 4: Gradients.

• –

Theorem 4.1 (Gradient of Weighted Mean Curvature).

• •

Section 5: Violations of Continuity.

• –

Theorem 5.1 (Continuity of Gradient).

• •

Section 6: Discussion.