Homology, Homotopy and Applications, vol.5(2), 2003, pp.19–37

COMPUTING LINKING NUMBERS OF A FILTRATION

HERBERT EDELSBRUNNER and AFRA ZOMORODIAN

(communicated by Gunnar Carlsson)

Abstract We develop fast algorithms for computing the linking number of a simplicial complex within a filtration. We give experimental results in applying our work toward the detection of non-trivial tangling in biomolecules, modeled as alpha complexes.

1. Introduction

In this paper, we develop fast algorithms for computing the linking numbers of simplicial complexes. Our work is within a framework of applying computational topology methods to the fields of biology and chemistry. Our goal is to develop useful tools by researchers in computational structural biology.

1.0.0.1. Motivation and Approach. In the 1980’s, it was shown that the DNA, the molecular structure of the genetic code of all living organisms, can become knotted during replication [1]. This finding initiated interest in knot theory among biologists and chemists for the detection, synthesis, and analysis of knotted molecules [8]. The impetus for this research is that molecules with non-trivial topological attributes often display exotic chemistry. Taylor recently discovered a figure-of-eight knot in the structure of a plant protein by examining 3,440 proteins using a computer program [19]. Moreover, chemical self-assembly units have been used to create catenanes, chains of interlocking molecular rings, and rotaxanes, cyclic molecules threaded by linear molecules. Researchers are building nano-scale chemical switches and logic gates with these structures [2, 3]. Eventually, chemical computer memory systems could be built from these building blocks.

Catenanes and rotaxanes are examples of non-trivial structural tanglings. Our work is on detecting such interlocking structures in molecules through a combinatorial method, based on algebraic topology. We model biomolecules as a sequence of alpha complexes [7]. The basic assumption of this representation is that an alpha-complex sequence captures the topological features of a molecule. This sequence is also a filtration of the Delaunay triangulation, a well-studied combinatorial object, enabling the development of fast algorithms.

The focus of this paper is the linking number. Intuitively, this invariant detects if components of a complex are linked and cannot be separated. We hope to eventually incorporate

Research by both authors is partially supported by ARO under grant DAAG55-98-1-0177. Research by the first author is also partially supported by NSF under grant CCR-97-12088, EIA-99-72879, and CCR-00-86013. Received October 1, 2001, revised November 1, 2002; published on April 22, 2003. 2000 Mathematics Subject Classification: 55N99, 55U10, 55M99. Key words and phrases: Computational geometry and topology, knots, linking number, three-manifolds, filtrations, alpha shapes, algorithms. \begin{matrix}c \\ ◯ \end{matrix} 2003, Herbert Edelsbrunner and Afra Zomorodian. Permission to copy for private use granted.

our algorithm into publicly available software as a tool for detecting existence of interlocked molecular rings.

Given a filtration, the main contributions of this paper are:

(i) the extension of the definition of the linking number to graphs, using a canonical basis,

(ii) an algorithm for enumerating and generating all cycles and their spanning surfaces

within a filtration, (iii) data structures for efficient enumeration of co-existing pairs of cycles in different

components, (iv) an algorithm for computing the linking number of a pair of cycles,

(v) and the implementation of the algorithms and experimentation on real data sets.

Algorithm (iv) is based on spanning surfaces of cycles, giving us an approximation to the linking number in the case of non-orientable or self-intersecting surfaces. Such cases do not arise often in practice, as shown in Section 6. However, we note in Section 2 that the linking number of a pair may be also computed by alternate algorithms. Regardless of the approach taken, pairs of potentially linked cycles must be first detected and enumerated. We provide the algorithms and data structures of such enumeration in (i-iii).

1.0.0.2. Prior work. Important knot problems were shown to be decidable by Haken in his seminal work on normal surfaces [10]. This approach, as reformulated by Jaco and others [13], forms the basis of many current knot detection algorithms. Haas et al. recently showed that these algorithms take exponential time in the number of crossings in a knot diagram [12]. They also placed both the unknotting problem and the splitting problem in NP, the latter being the focus of our paper. Generally, other approaches to knot problems have unknown complexity bounds, and are assumed to take at least exponential time. As such, the state of the art in knot detection only allows for very small data sets. We refer to Adams [1] for background in knot theory.

Three-dimensional alpha shapes and complexes may be found in Edelsbrunner and Mücke [7]. We modify the persistent homology algorithm to compute cycles and surfaces [6]. We refer to Munkres [15] for background in homology theory that is accessible to non-specialists.

1.0.0.3. Outline. The remainder of this paper is organized as follows. We review linking numbers for collections of closed curves, and extend this notion to graphs in ℝ^{3} in Section 2. We describe our model for molecules in Section 3. Extending the persistence algorithm, we design basic algorithms in Section 4 and use them to develop an algorithm for computing linking numbers in Section 5. We show results of some initial experiments in Section 6, concluding the paper in Section 7.

2. Linking Number

In this section, we define links and discuss two equivalent definitions of the linking number. While the first definition provides intuition, the second definition is the basis of our computational approach.

2.0.0.4. Links. A knot is an embedding of a circle in three-dimensional Euclidean space, k : 𝕊^{1} → ℝ^{3}. Two knots are equivalent if there is an ambient isotopy that maps the first to the second. That is, we may deform the first to the second by a continuous motion that does not cause self-intersections. A link l is a collection of knots with disjoint images. A link is separable (splittable) if it can be continuously deformed so that one or more components can be separated from other components by a plane that itself does not intersect any of the components. We often visualize a link l by a link diagram, which is the projection of a link onto a plane such that the over- and under-crossings of knots are presented clearly, We give an example in Figure 1. For a formal definition, see [12].

Figure 1: A link diagram for the Whitehead link. Vertices occur at crossings and are labeled according to the convention in Figure 2.

2.0.0.5. Linking number. A knot (link) invariant is a function that assigns equivalent objects to equivalent knots (links.) Seifert first defined an integer link invariant, the linking number, in 1935 to detect link separability [18]. Given a link diagram for a link l, we choose orientations for each knot in l. We then assign integer labels to each crossing between any pair of knots k, k′, following the convention in Figure 2. Let λ(k, k′) of the pair of knots

Figure 2: The crossing label is +1 if the rotation of the overpass by 90 degrees counterclockwise aligns its direction with the underpass, and −1 otherwise.

be one half the sum of these labels. A standard argument using Reidemeister moves shows that λ is an invariant for equivalent pairs of knots up to sign [1]. The linking number λ(l) of a link l is

λ(l) = ∑_{k≠k′∈l} |λ(k, k′)|.

We note that λ(l) is independent of knot orientations. Also, the linking number does not completely recognize linking. The Whitehead link in Figure 1, for example, has linking number zero, but is not separable.

2.0.0.6. Surfaces. The linking number may be equivalently defined by other methods, including one based on surfaces [17]. A spanning surface for a knot k is an embedded surface with boundary k. An orientable spanning surface is a Seifert surface. Because it is orientable, we may label its two sides as positive and negative. We show examples of such surfaces for the Hopf link in Figure 3.

Figure 3: The Hopf link and Seifert surfaces of its two unknots. Clearly, λ = 1. This link is the 200th complex for data set H in Section 6.

Given a pair of oriented knots k, k′, and a Seifert surface s for k, we label s by using the orientation of k. We then adjust k′ via a homotopy h until it meets s in a finite number of points. Following along k′ according to its orientation, we add +1 whenever k′ passes from the negative to the positive side, and −1 whenever k′ passes from the positive to the negative side. The following lemma asserts that this sum is independent of our the choice of h and s, and it is, in fact, the linking number.

Seifert Surface Lemma.  λ(k, k′)  is  the  sum  of  the

signed intersections between k′ and any Seifert surface for k.

The proof is by a standard Seifert surface construction [17]. If the spanning surface is nonorientable, we can still count how many times we pass through the surface, giving us the following weaker result.

Spanning Surface Lemma. λ(k, k′) (mod 2) is the parity of the number of times k′ passes through any spanning surface for k.

2.0.0.7. Graphs. We need to extend the linking number to graphs, in order to use the above lemma for computing linking numbers for simplicial complexes. Let G = (V, E), E ⊆ \binom{V}{2} be a simple undirected graph in ℝ^{3} with c components G^{i}. Let z_{1}, …, z_{m} be a fixed basis for the cycles in G, where m = |E| − |V | + c. We then define the linking number between two components of G to be λ(G^{i}, G^{j}) = |λ(z_{p}, z_{q})| for all cycles z_{p}, z_{q} in G^{i}, G^{j}, respectively. The linking number of G is then defined by combining the total interaction between pairs of components:

λ(G) = ∑_{i≠j} λ(G^{i}, G^{j}).

The linking number is computed only between pairs of components following Seifert’s original definition. Linked cycles within the same component may be easily unlinked by a homotopy. Figure 4 shows that the linking number for graphs is dependent on the chosen basis. While it may seem that we want λ(G) = 1 in the figure, there is no clear answer Figure 4: We get different λ(G) for a graph G (top) depending on our choice of basis for G^{2}: two small cycles (left) or one large and one small cycle (right.)

in general. We will define a canonical basis in Section 4 using the persistent homology algorithm to compute λ(G) for simplicial complexes.

3. Alpha Complexes

Our approach to analyzing a topological space is to assume a filtration for such a space. A filtration may be viewed as a history of a growing space that is undergoing geometric and topological changes. While filtrations may be obtained by various methods, only meaningful filtrations give meaningful linking numbers. As such, we use alpha complex filtrations to model molecules. The alpha complex captures the connectivity of a molecule that is represented by a union of spheres. This model may be viewed as the dual of the space filling model for molecules [14].

3.0.0.8. Dual complex. A spherical ball û = (u, U^{2}) ∈ ℝ^{3} × ℝ is defined by its center u and square radius U^{2}. If U^{2} < 0, the radius is imaginary and so is the ball. The weighted distance of a point x from a ball û is πû(x) = ∥x − u∥^{2} − U^{2}. Note that a point x ∈ ℝ^{3} belongs to the ball iff πû(x) ⩽ 0, and it belongs to the bounding sphere iff πû(x) = 0. Let S be a finite set of balls. The Voronoi region of û ∈ S is the set of points for which û minimizes the weighted distance,

Vû = {x ∈ ℝ^{3} | πû(x) ⩽ πv̂(x), ∀v̂ ∈ S}.

The Voronoi regions decompose the union of balls into convex cells of the form û ∩ Vû, as illustrated in Figure 5. Any two regions are either disjoint or they overlap along a shared portion of their boundary. We assume general position: at most four Voronoi regions can have a non-empty common intersection. Let T ⊆ S have the property that its Voronoi regions have a non-empty common intersection, and consider the convex hull of the corresponding centers, σ_{T} = conv {u | û ∈ T }. General position implies that σ_{T} is a d-dimensional simplex, where d = card T − 1. The dual complex of S is the collection of simplices constructed in this manner,

K = {σ_{T} | T ⊆ S, ⋂_{û∈T}(û ∩ Vû) ≠ ∅}.

Figure 5: Union of nine disks, convex decomposition using Voronoi regions, and dual complex.

Any two simplices in K are either disjoint or they intersect in a common face which is a simplex of smaller dimension. Furthermore, if σ ∈ K, then all faces of σ are simplices in K. A set of simplices with these two properties is a simplicial complex [15]. A subcomplex is a subset L ⊆ K that is itself a simplicial complex.

3.0.0.9. Alpha complex. A filtration ordering is an ordering of a set of simplices such that each prefix of the ordering is a subcomplex. The sequence of subcomplexes defined by taking successively larger prefixes is the corresponding filtration. For dual complexes of a collection of balls, we generate an ordering and a filtration by literally growing the balls. For every real number α^{2} ∈ ℝ, we increase the square radius of a ball û by α^{2}, giving us û(α) = (u, U^{2} + α^{2}). We denote the collection of expanded balls û(α) as S(α). If U^{2} = 0, then α is the radius of û(α). If α^{2} < 0, then α is imaginary, and so is the ball û(α). The α-complex K(α) of S is the dual complex of S(α) [7]. For example, K(−∞) = ∅, K(0) = K, and K(∞) = D is the dual of the Voronoi diagram, also known as the Delaunay triangulation of S. For each simplex σ ∈ D, there is a unique birth time α^{2}(σ) defined such that σ ∈ K(α) iff α^{2} ⩾ α^{2}(σ). We order the simplices such that α^{2}(σ) < α^{2}(τ) implies σ precedes τ in the ordering. More than one simplex may be born at a time and such cases may arise even if S is in general position. In the case of a tie, it is convenient to order lower-dimensional simplices before higher-dimensional ones, breaking remaining ties arbitrarily. We call the resulting sequence the age ordering of the Delaunay triangulation.

3.0.0.10. Modeling molecules. To model molecules by alpha complexes, we use representations of molecules as unions of balls. Each ball is an atom, as defined by its position in space and its van der Waals radius. These atoms become the spherical balls we need to define our complexes. Our representation gives us a filtration of alpha complexes for each molecule, as shown in Figure 6. We compute a linking number for each complex in a filtration of m complexes. Let [m] denote the set {1, 2, …, m}. Then, the linking number may be viewed as a signature function λ : [m] → ℤ that maps each index i ∈ [m] to an integer λ(i) ∈ ℤ. For other signature functions for filtrations of alpha complexes, see [5, 7].

Figure 6: Six complexes in the filtration of 42,787 complexes for data set Z in Section 6.

4. Basis and Surfaces

To compute the linking numbers for an alpha complex, we need to recognize cycles, establish a basis for the set of cycles, and find spanning surfaces for the basis cycles. We do so by extending an algorithm we developed for computing persistent homology [6]. We dispense with defining persistence and concentrate on the algorithm and its extension.

4.0.0.11. Homology. We use homology to define cycles in a complex. Homology partitions cycles into equivalence classes using the boundary class of bounding cycles as the null element of a quotient group in each dimension. We use ℤ_{2} homology, so the group operation, which we call addition, is symmetric difference. Addition allows us to combine sets of simplices in a way that eliminates shared boundaries, as shown in Figure 7. Intuitively, non-bounding 1-cycles correspond to the graph notion of a cycle. We need to

Figure 7: Symmetric difference in dimensions one and two. We add two 1-cycles to get a new 1-cycle. We add the surfaces the cycles bound to get a spanning surface for the new 1-cycle.

define a basis for the first homology group of the complex which contains all 1-cycles, and choose representatives for each homology class. We use these representatives to compute linking numbers for the complex.

A simplex of dimension d in a filtration either creates a d-cycle or destroys a (d − 1)- cycle by turning it into a boundary. We mark simplices as positive or negative, according to this action [5]. In particular, edges in a filtration which connect components are marked as negative. The set of all negative edges gives us a spanning tree of the complex, as shown in Figure 8. We use this spanning tree to define our canonical basis. Every time a positive edge Figure 8: Solid negative edges combine to form a spanning tree. The dashed positive edge σ_{i} creates a canonical cycle.

σ_{i} is added to the complex, it creates a new cycle. We choose the unique cycle that contains σ_{i} and no other positive edge as a new basis cycle. We call this cycle a canonical cycle, and the collection of canonical cycles, the canonical basis. We use this basis for computation.

4.0.0.12. Persistence. The  persistence  algorithm  matches  positive

and negative simplices to find life-times of homological cycles in a filtration. The algorithm does so by following a representative cycle z for each class. Initially, z is the boundary of a negative simplex σ_{j}, as z must lie in the homology class σ_{j} destroys. The algorithm then successively adds class-preserving boundary cycles to z until it finds the matching positive simplex σ_{i}, as shown in Figure 9. We call the half-open interval [i, j) the persistence in-

Figure 9: Starting from the boundary of the negative triangle σ_{j}, the persistence algorithm finds a matching positive edge σ_{i} by finding the dashed 1-cycle. We modify this 1-cycle further to find the solid canonical 1-cycle and a spanning surface.

terval of both the homology class and its canonical representative. During this interval, the homology class exists as a class of homologous non-boundings cycles in the filtration. As such, the class may only affect the linking numbers of complexes K_{i}, …, K_{j−1} in the filtration. We use this insight in the next section to design an algorithm for computing linking numbers.

4.0.0.13. Computing canonical cycles. The persistence algorithm halts when it finds the matching positive simplex σ_{i} for a negative simplex σ_{j}, often generating a cycle z with multiple positive edges and multiple components. We need to convert z into a canonical cycle by eliminating all positive edges in z except for σ_{i}. We call this process canonization. To canonize a cycle, we add cycles associated with unnecessary positive edges to z successively, until z is composed of σ_{i} and negative edges, as shown in Figure 9. Canonization amounts to replacing one homology basis element with a linear combination of other elements in order to reach the unique canonical basis we defined earlier. A cycle undergoing canonization changes homology classes, but the rank of the basis never changes.

4.0.0.14. Computing spanning surfaces. For each canonical cycle, we need a spanning surface in order to compute linking numbers. We may compute these by maintaining surfaces while computing the cycles. Recall that initially, a cycle representative is the boundary of a negative simplex σ_{j}. We use σ_{j} as the initial spanning surface for z. Every time we add a cycle y to z in the persistence algorithm, we also add the surface bounded by y to the z’s surface. We continue this process through canonization to produce both canonical cycles and their spanning surfaces. Here, we are using a crucial property of our filtrations: the final complex is always the Delaunay complex of the set of weighted points and does not contain any 1-cycles. Therefore, all 1-cycles are eventually turned to boundaries and have spanning surfaces.

If the generated spanning surface is Seifert, we may apply the Seifert Surface Lemma to compute the linking numbers. In some cases, however, the spanning surface is not Seifert, as in Figure 10. In these cases, we may either compute the linking number modulus 2 by applying the Spanning Surface Lemma, or compute the linking number by alternative methods.

Figure 10: The spanning surface produced for the cycle which is the boundary of a Möbius strip is non-orientable.

5. Algorithm

In this section, we use the basis and spanning surfaces computed for 1-cycles to find linking numbers for all complexes in a filtration. Since we focus on 1-cycles only, we will refer to them simply as cycles.

5.0.0.15. Overview. We assume a filtration K_{1}, K_{2}, …, K_{m} as input, which we alternately view as a single complex undergoing growth. As simplices are added, the complex undergoes topological changes which affect the linking number: new components are created and merged together, and new non-bounding cycles are created and eventually destroyed. We use a basic insight from the last section: a basis cycle z with persistence interval [i, j) may only affect the linking numbers of complexes K_{i}, K_{i+1}, …, K_{j−1} in the filtration, Consequently, we only need to consider basis cycles z′ that exist during some subinterval [u, v) ⊆ [i, j) in a different component than z’s. We call the pair z, z′ a potentiallylinked (p-linked) pair of basis cycles, and the interval [u, v) the p-linking interval.

Focusing on p-linked pairs, we get an algorithm with three phases. In the first phase, we compute all p-linked pairs of cycles. In the second phase, as shown in Figure 11, we compute the linking numbers of such pairs. In the third and final phase, we aggregate these contributions to find the linking number signature for the filtration.

for each p-linked pair z_{p}, z_{q} with interval [u, v) do Compute λ = |λ(z_{p}, z_{q})| ; Output (λ, [u, v))

               endfor.

Figure 11: Linking number algorithm.

Two  cycles  z_{p}, z_{q}  with  persistence  intervals  [i_{p}, j_{p}),

[i_{q}, j_{q}) co-exist during [r, s) = [i_{p}, j_{p}) ∩ [i_{q}, j_{q}). We need to know if these cycles also belong to different components during some sub-interval [u, v) ⊆ [r, s). Let t_{p,q} be the minimum index in the filtration when z_{p} and z_{q} are in the same component. Then, [u, v) = [r, s) ∩ [0, t_{p,q}). If [u, v) ≠ ∅, z_{p}, z_{q} are p-linked during that interval. In the remainder of this section, we will first develop a data structure for computing t_{p,q} for any pair of cycles z_{p}, z_{q}. Then, we use this data structure to efficiently enumerate all pairs of p-linked cycles. Finally, we give an algorithm for computing λ(z_{p}, z_{q}) for a p-linked pair of cycles z_{p}, z_{q}.

5.0.0.16. Component history. To compute t_{p,q}, we need to have a history of the changes to the set of components in a filtration. There are two types of simplices that can change this set. Vertices create components and are therefore all positive. Negative edges connect components. We construct a binary tree called component tree recording these changes using a union-find data structure [4]. The leaves of the component tree are the vertices of the filtration. When a negative edge connects two components, we create an internal node and connect it to the nodes representing these components, as shown in Figure 12. The component tree has size O(n) for n vertices, and we construct it in time

Figure 12: The union-find data structure (left) has vertices as nodes and negative edges as edges. The component tree (right) has vertices as leaves and negative edges as internal nodes.

O(nA^{−1}(n)), where A^{−1}(n) is the inverse of the Ackermann’s function which exhibits insanely slow growth. Having constructed the component tree, we find the time two vertices w, x are in the same component by finding their lowest common ancestor (lca) in this tree. We utilize Harel and Tarjan’s optimal method to find lca’s with O(n) preprocessing time and O(1) query time [11]. Their method uses bit operations. If such operations are not allowed, we may use van Leeuwen’s method with the same preprocessing time and O(log log n) query time [20].

5.0.0.17. Enumeration. Having constructed the component tree, we use a modified union-find data structure to enumerate all pairs of p-linked cycles. We augment the data structure to allow for quick listing of all existing canonical cycles in each component in K_{i}. Our augmentation takes two forms: we put the roots of the disjoint trees, representing components, into a circular doubly-linked list. We also store all existing cycles in each component in a doubly-linked list at the root node of the component, as shown in Figure 13. When components merge, the root x_{1} of one component becomes the parent of the root x_{2}

Figure 13: The augmented union-find data structure places root nodes in the shaded circular doubly-linked list. Each root node stores all active canonical cycles in that component in a doubly-linked list, as shown for the darker component.

of the other component. We concatenate the lists stored at the x_{1}, x_{2}, store the resulting list at x_{1}, and eliminate x_{2} from the circular list in O(1) time. When cycle z_{p} is created at time i, we first find z_{p}’s component in time O(A^{−1}(n)). Then, we store z_{p} at the root of the component and keep a pointer to z_{p} with simplex σ_{j}, which destroys z_{p}. This implies that we may delete z_{p} from the data structure at time j with constant cost.

Our algorithm to enumerate p-linked cycles is incremental. We add and delete cycles using the above operations from the union-find forest, as the cycles are created and deleted in the filtration. When a cycle z_{p} is created at time i, we output all p-linked pairs in which z_{p} participates. We start at the root which now stores z_{p} and walk around the circular list of roots. At each root x, we query the component tree we constructed in the last subsection to find the time t when the component of x merges with that of z_{p}. Note that t = t_{p,q} for all cycles z_{q} stored at x. Consequently, we can compute the p-linking interval for each pair z_{p}, z_{q} to determine if the pair is p-linked. If the filtration contains P p-linked pairs, our algorithm takes time O(mA^{−1}(n) + P), as there are at most m cycles in the filtration.

5.0.0.18. Orientation. In the previous section, we showed how one may compute spanning surfaces s_{p}, s_{q} for cycles z_{p}, z_{q}, respectively. To compute the linking number using our lemma, we need to orient either the pair s_{p}, z_{q} or z_{p}, s_{q}. Orienting a cycle is trivial: we orient one edge and walk around to orient the cycle. If either surface has no selfintersections, we may easily attempt to orient it by choosing an orientation for an arbitrary triangle on the surface, and spreading that orientation throughout. The procedure either orients the surface or classifies it as non-orientable. We currently do not have an algorithm for orienting surfaces with self-intersections. The main difficulty is distinguishing between two cases for a self-intersection: a surface touching itself and passing through itself, as shown in Figure 14.

Figure 14: A surface self-intersection viewed from its side. We cannot resolve it as the surface touching or passing through itself.

5.0.0.19. Computing λ. We now show how to compute λ(z_{p}, z_{q}) for a pair of p-linked cycles z_{p}, z_{q}, completing the description of our algorithm in Figure 11. We assume that we have oriented s_{p}, z_{q} for the remainder of this subsection.

Let the star of a vertex u St u be the set of simplices containing u as a vertex. We subdivide the complex via a barycentric subdivision by connecting the centroid of each triangle to its vertices and midpoints of its edges, subdividing the simplices accordingly. This subdivision guarantees that no edge uv will have both ends on a Seifert surface unless it is entirely contained in that surface. We note that this approach mimics the construction of regular neighborhoods for complexes [9].

For a vertex u ∈ s_{p}, the edge property guaranteed by subdivision enables us to mark each edge uv ∈ St u, v ∉ s_{p} as positive or negative, depending on the location of v with respect to s_{p}. We show an example of this marking in Figure 15. After marking edges, we

Figure 15: Edges uv ∈ St u, u ∈ s_{p}, v ∉ s_{p} are marked + or − depending on where they end relative to the oriented Seifert surface s_{p}.

walk once around z_{q}, starting at a vertex not on s_{p}. If such a vertex does not exist, then λ(z_{p}, z_{q}) = 0. Otherwise, we create a string S_{p,q} of + and − characters by noting the marking of edges during our walk. S_{p,q} has even length as we start and end our walk on a vertex not on s_{p}, and each intersection of z_{q} with s_{p} produces a pair of characters, as Figure 16: Starting at v, we walk on z_{q} according to its orientation. Segments of z_{q} that intersect s_{p} are shown, along with their contribution to S_{p,q} = “ + + + + + − − −”. We get λ(z_{p}, z_{q}) = −1.

shown in Figure 16. If S_{p,q} is the empty string, z_{q} never intersects s_{p} and λ(z_{p}, z_{q}) = 0. Otherwise, z_{q} passes through s_{p} for pairs +− and −+, corresponding to z_{q} piercing the positive or negative side of s_{p}, respectively. Scanning S_{p,q} from left to right in pairs, we add +1 for each occurrence of −+, −1 for each +−, and 0, for each ++ or −−. Applying the Seifert Surface Lemma in Section 2, we see that this sum is λ(z_{p}, z_{q}).

5.0.0.20. Computing λ mod 2. If neither of the spanning surfaces s_{p}, s_{q} of the two cycles z_{1}, z_{2} is Seifert, we may still compute λ(z_{1}, z_{2}) mod 2 by a modified algorithm, provided one surface, say s_{p}, has no self-intersections. We choose an orientation on s_{p} locally, and extend it until all the stars of the original vertices are oriented. are oriented. This orientation will not be consistent globally, resulting in pair of adjacent vertices in s_{p} with opposite orientations. We call the implicit boundary between vertices with opposite orientations a flip curve, as shown in bold in Figure 17. When a cycle segment crosses the

Figure 17: The bold flip curve is the border of s_{p}^{+} and s_{p}^{−}, the portions of s_{p} that are oriented differently. S_{p,q} = “ + + − − − + − − − ”, so counting all +’s, we get λ(z_{p}, z_{q}) mod 2 = 3 mod 2 = 1. flip curve, orientation changes. Therefore, in addition to noting marked edges, we add a + to the string S_{p,q} every time we cross a flip line. To compute λ(z_{p}, z_{q}) mod 2, we only count +’s in S_{p,q} and take the parity as our answer.

If s_{p} is orientable, there are no flip curves on it. The contribution of cycle segments to the string is the same as before: +− or −+ for segments that pass through s_{p}, and ++ and −− for segments that do not. By counting +’s, only segments that pass through s_{p} change the parity of the sum for λ. Therefore, the algorithm computes λ mod 2 correctly for orientable surfaces. For the orientable surface on the right in Figure 16, for instance, we get λ(z_{p}, z_{q}) mod 2 = 5 mod 2 = 1, which is equivalent to the parity of the answer computed by the previous algorithm.

5.0.0.21. Remark. We are currently examining the question of orienting surfaces with self-intersections. Using our current methods, we may obtain a lower bound signature for λ by computing a mixed sum: we compute λ and λ mod 2 whenever we can to obtain the approximation. We may also develop other methods, including those based on the projection definition of the linking number in Section 2.

6. Experiments

In this section, we present some experimental timing results and statistics which we used to guide our algorithm development. We also provide visualizations of basis cycles in a filtration. All timings were done on a Micron PC with a 266 MHz Pentium II processor and 128 MB RAM running Solaris 8.

6.0.0.22. Implementation. We have implemented all the algorithms in the paper, except for the algorithm for computing λ mod 2. Our implementation differs from our exposition in three ways. The implemented component tree is a standard union-find data structure with the union by rank heuristic, but no path compression [4]. Although this structure has an O(n log n) construction time and an O(log n) query time, it is simple to implement and extremely fast in practice. We also use a heuristic to reduce the number of p-linked cycles. We store bounding boxes at the roots of the augmented union-find data structure. Before enumerating p-linked cycles, we check to see if the bounding box of the new cycle intersects with that of the stored cycles. If not, the cycles cannot be linked, so we obviate their enumeration. Finally, we only simulate the barycentric subdivision.

6.0.0.23. Data. We have experimented with a variety of data sets and show the results for six representative sets in this section. The first data set contains points regularly sampled along two linked circles. The resulting filtration contains a complex which is a Hopf link, as shown in Figure 3. The other data sets represent molecular structures with weighted points. In each case, we first compute the weighted Delaunay triangulation and the age ordering of that triangulation. The data points become vertices or 0-simplices. Table 1 gives the sizes of the data sets, their Delaunay triangulations, and age orderings. We show renderings of specific complexes in the filtration for data set K in Figure 18.

# simplices of dimension d

0  1  2  3  total
H  100 1,752  3,240 1,587  6,679
G  318 2,322  3,978 1,973  8,591
M 1,001 7,537 13,018 6,481 28,037
Z 1,296 11,401 20,098 9,992 42,787
K 2,370 17,976 31,135 15,528 67,009
D 7,774 60,675 105,710 52,808 226,967

Table 1: H defines a Hopf link. G is Gramicidin A, a small protein. M is a protein monomer. Z is a portion of a periodic zeolite structure. K is a human cyclin-dependent kinase. D is a DNA tile.

Figure 18: Complex K_{8168} of K has two components and seventeen cycles. The spanning surfaces are rendered transparently.

6.0.0.24. Basis. Table 2 summarizes the basis generation process. We distinguish the two steps of our algorithm: initial basis generation and canonization. We give the number of basis cycles for the entire filtration, which is equal to the number of positive edges. We show the effect of canonization on the size of the cycles and their spanning surfaces in Table 3. Note that canonization increases the size of cycles by one or two orders of magnitude. This is partially the reason we try to avoid performing the link detection if possible.

6.0.0.25. Links. In Table 4, we show that our component tree and augmented trees are very fast in practice to generate p-linked pairs. We also show that our bounding box heuristic for reducing the number of p-linked pairs increases the computation time negligibly. The heuristic is quite successful, moreover, in reducing the number of pairs we have to check for linkage It eliminates 99.8% of the candidates for dataset Z, for example, as shown in Table 5. The differences in total time of computation reflect the basic structure of the datasets. Dataset D has a large computation time, for instance, as the average size of the p-linked surfaces is approximately 264.16 triangles, compared to about 1.88 triangles

time in seconds

generate canonize total # cycles
H  0.08  0.04 0.12  1,653
G  0.08  0.03 0.11  2,005
M  0.28  0.20 0.48  6,537
Z  0.46  0.46 0.92  10,106
K  0.72  1.01 1.73  15,607
D  2.63  2.94 5.57  52,902

Table 2: Time to generate and canonize basis cycles, as well as their number.

avg cycle length avg surface size
before after before after
H  3.06 51.03  1.06 63.04
G  3.26 13.02  1.38 52.28
M  3.29 34.18  1.33 71.18
Z  4.71 25.33  3.26 117.81
K  3.48 67.87  1.62 166.70
D  3.46 39.94  1.81 158.99

Table 3: Average number of edges per cycle and number of triangles per spanning surface, before and after canonization.

for dataset K, and about 1.73 triangles for dataset M.

6.0.0.26. Discussion. Our initial experiments demonstrate the feasibility of the algorithms for fast computation of linking. The experiments fail to detect any links in the protein data, however. This is to be expected, as a protein consists of a single component, the primary structure of a protein being a single polypeptide chain of amino acids. Links, on the other hand, exist in different components by definition. We may relax this definition easily, however, to allow for links occuring in the same component. We have implementations of algorithms corresponding to this relaxed definition. Our future plans include looking for links in proteins from the Protein Data Bank [16]. Such links could occur naturally as a

tree alg heur links
H 0.01 0.00 0.00 0.01
G 0.00 0.01 0.02 0.02
M 0.03 0.06 0.06 0.23
Z 0.04 0.07 0.07 0.13
K 0.06 0.13 0.16 0.36
D 0.27 0.56 0.82 8.22

Table 4: Time in seconds to construct the component tree, and enumerate p-linked pairs (alg), p-linked pairs with intersecting bounding boxes (heur), and links.

alg  heur links
H  1  1  1
G  112  0  0
M 16,503 14,968  0
Z 169,594  308  0
K  12,454 11,365  0
D  98,522 4,448  0

Table 5: Number of p-linked pairs (alg), p-linked pairs with intersecting bounding boxes (heur), and links.

result of disulphide bonds between different residues in a protein.

7. Conclusion

In this paper, we develop algorithms for finding the linking numbers of a filtration. We give algorithms for computing bases of 1-cycles and their spanning surfaces in simplicial complexes, and enumerating co-existing cycles in different components. In addition, we present an algorithm for computing the linking number of a pair of cycles using the surface formulation. Our implementations show that the algorithms are fast and feasible in practice. By modeling molecules as filtrations of alpha complexes, we can detect potential nontrivial tangling within molecules. Our work is within a framework for applying topological methods for understanding molecular structures.

7.0.0.27. Acknowledgments. We thank David Letscher for discussions during the early stages of this work. We also thank Daniel Huson for the zeolite dataset Z, and Thomas La Bean for the DNA tile data set D.

Homology, Homotopy and Applications, vol. 5(2), 2003 36

References

[1] Adams, C. C. The Knot Book: An Elementary Introduction to the Mathematical Theory of Knots. W. H. Freeman and Company, New York, NY, 1994.

[2] Bissell, R. A., Córdova, E., Kaifer, A. E., and Stoddart, J. F. A checmically and electrochemically switchable molecular shuttle. Nature 369 (1994), 133– 137.

[3] Collier, C. P., Wong, E. W., Belohradský, Raymo, F. M., Stoddart, J. F., Kuekes, P. J., Williams, R. S., and Heath, J. R. Electronically configurable moleculear-based logic gates. Science 285 (1999), 391–394.

[4] Cormen, T. H., Leiserson, C. E., and Rivest, R. L. Introduction to Algorithms. The MIT Press, Cambridge, MA, 1994.

[5] Delfinado, C. J. A., and Edelsbrunner, H. An incremental algorithm for Betti numbers of simplicial complexes on the 3-sphere. Comput. Aided Geom. Design 12 (1995), 771–784.

[6] Edelsbrunner, H., Letscher, D., and Zomorodian, A. Topological persistence and simplification. In Proc. 41st Ann. IEEE Sympos. Found. Comput. Sci. (2000), pp. 454–463.

[7] Edelsbrunner, H., and Mücke, E. P. Three-dimensional alpha shapes. ACM Trans. Graphics 13 (1994), 43–72.

[8] Flapan, E. When Topology Meets Chemistry : A Topological Look at Molecular Chirality. Cambridge University Press, New York, NY, 2000.

[9] Giblin, P. J. Graphs, Surfaces, and Homology, second ed. Chapman and Hall, New York, NY, 1981.

[10] Haken, W. Theorie der Normalflächen. Acta Math. 105 (1961), 245–375.

[11] Harel, D., and Tarjan, R. E. Fast algorithms for finding nearest common ancestors. SIAM J. Comput. 13 (1984), 338–355.

[12] Hass, J., Lagarias, J. C., and Pippenger, N. The computational complexity of knot and link problems. J. ACM 46 (1999), 185–211.

[13] Jaco, W., and Tollefson, J. L. Algorithms for the complete decomposition of a closed 3-manifold. Illinois J. Math. 39 (1995), 358–406.

[14] Leach, A. R. Molecular Modeling, Principles and Applications. Pearson Education Limited, Harlow, England, 1996.

[15] Munkres, J. R. Elements of Algebraic Topology. Addison-Wesley, Redwood City, California, 1984.

[16] RCSB. Protein data bank. http://www.rcsb.org/pdb/.

[17] Rolfsen, D. Knots and Links. Publish or Perish, Inc., Houston, Texas, 1990.

[18] Seifert, H. Über das Geschlecht von Knoten. Math. Annalen 110 (1935), 571–592.

[19] Taylor, W. R. A deeply knotted protein structure and how it might fold. Nature 406 (2000), 916–919.

[20] van Leeuwen, J. Finding lowest common ancestors in less than logarithmic time. Unpublished report.

This article may be accessed via WWW at http://www.rmi.acnet.ge/hha/ or by anonymous ftp at

ftp://ftp.rmi.acnet.ge/pub/hha/volumes/2003/n2a2/v5n2a2.(dvi,ps,pdf)

Herbert Edelsbrunner edels@cs.duke.edu Department of Computer Science, Duke University, Durham, and Raindrop Geomagic, Research Triangle Park, NC.

Afra Zomorodian afra@cs.stanford.edu

Department of Computer Science, University of Illinois, Urbana, IL.