Parallel Computation of Alpha Complex for Biomolecules

TL;DR

This paper introduces a GPU-based parallel algorithm for computing the alpha complex of biomolecules, significantly improving performance by up to 50 times over existing methods without needing prior Delaunay triangulation.

Contribution

The novel algorithm computes the alpha complex directly from biomolecular data using a two-stage process, bypassing the need for Delaunay triangulation and leveraging biomolecular spatial properties.

Findings

Up to 50x speedup compared to existing methods

Effective in handling various biomolecular structures

Eliminates the need for prior Delaunay triangulation

Abstract

The alpha complex, a subset of the Delaunay triangulation, has been extensively used as the underlying representation for biomolecular structures. We propose a GPU-based parallel algorithm for the computation of the alpha complex, which exploits the knowledge of typical spatial distribution and sizes of atoms in a biomolecule. Unlike existing methods, this algorithm does not require prior construction of the Delaunay triangulation. The algorithm computes the alpha complex in two stages. The first stage proceeds in a bottom-up fashion and computes a superset of the edges, triangles, and tetrahedra belonging to the alpha complex. The false positives from this estimation stage are removed in a subsequent pruning stage to obtain the correct alpha complex. Computational experiments on several biomolecules demonstrate the superior performance of the algorithm, up to a factor of 50 when compared to existing methods that are optimized for biomolecules.

Tables

Table 1. Table 1 : Runtime comparison of the proposed algorithm with gReg3D on an nVidia GTX 660 Ti graphics card. Timings are reported in milliseconds. % percent \% Simplex refers to the size of the alpha complex as a percentage of the size of the weighted Delaunay triangulation. The last column shows the speedup in runtime of our algorithm over gReg3D . ‘*’ indicates the data was partitioned and processed in chunks. ‘–’ indicates that the code could not execute due to insufficient memory.

$\alpha$	PDB id	$\#$Atoms	$K_{\alpha }$		gReg3D		$\%$Simplex	Speed up
			$\#$Simplices	Time(ms)	$\#$Simplices	Time(ms)
0.0	1GRM	260	932	13	6295	117	14.8	9.0
	1U71	1505	5696	13	40878	115	13.9	11.1
	3N0H	1509	5739	14	41244	137	13.9	10.0
	4HHB	4384	38796	29	150141	193	25.8	6.6
	2J1N	8142	29642	18	227719	229	13.0	12.7
	1K4C	16068	62851	27	446383	347	14.1	12.9
	2OAU	16647	123175	56	466586	344	26.4	6.2
	1AON	58674	262244	65	1650841	879	15.9	13.5
	1X9P*	217920	924086	113	6142811	2555	15.0	22.6
	1IHM*	677040	2713083	277	–	–	–	–
	4CWU*	5905140	23450403	2709	–	–	–	–
	3IYN*	5975700	24188892	2874	–	–	–	–
1.0	1GRM	260	1598	15	6295	117	25.4	7.9
	1U71	1505	10828	17	40878	115	26.5	8.5
	3N0H	1509	10965	30	41244	137	26.6	4.6
	4HHB	4384	65987	86	150141	193	44.0	2.2
	2J1N	8142	58205	30	227719	229	25.6	7.6
	1K4C	16068	118467	52	446383	347	26.5	6.7
	2OAU	16647	199101	159	466586	344	42.7	2.2
	1AON	58674	495683	160	1650841	879	30.0	5.5
	1X9P*	217920	1653778	196	6142811	2555	26.9	13.0
	1IHM*	677040	5058507	605	–	–	–	–
	4CWU*	5905140	44411353	5118	–	–	–	–
	3IYN*	5975700	45790463	5501	–	–	–	–

Table 2. Table 2 : Time spent within different steps of the algorithm. Timings are reported in milliseconds for memory transfer, grid computation, computing potential simplices, and pruning. The last column shows the total time taken for all steps.

$\alpha$	PDB id	$\#$Atoms	$\#$Simplices	Memory	Grid	Potential Simplices			Pruning			Total
						Edges	Tris	Tets	Tets	Tris	Edges	time
0.0	1GRM	260	932	0.8	1.0	2.8	1.1	1.0	1.2	3.1	1.9	13.0
	1U71	1505	5696	0.9	0.8	2.3	1.7	1.0	1.8	2.6	1.9	13.0
	3N0H	1509	5739	0.7	0.9	2.3	1.4	2.5	1.8	2.5	1.5	13.7
	4HHB	4384	38796	1.1	0.8	2.7	2.0	6.0	8.3	4.6	3.7	29.2
	2J1N	8142	29642	1.1	1.2	3.7	1.6	1.3	2.0	3.9	3.2	18.1
	1K4C	16068	62851	1.7	2.0	4.3	1.5	1.6	3.8	7.6	4.3	26.9
	2OAU	16647	123175	2.1	1.3	4.7	4.3	5.8	21.9	9.5	6.0	55.5
	1AON	58674	262244	4.6	2.8	11.1	5.6	4.1	16.7	10.9	9.4	65.2
1.0	1GRM	260	1598	1.2	1.4	2.7	1.2	2.0	1.8	2.7	1.8	14.7
	1U71	1505	10828	0.9	0.8	2.3	1.5	3.4	2.5	3.2	2.5	17.1
	3N0H	1509	10965	0.9	1.4	2.4	1.7	14.6	2.6	3.5	2.9	30.0
	4HHB	4384	65987	1.5	1.9	3.4	5.4	23.7	32.8	12.4	4.8	86.0
	2J1N	8142	58205	1.4	1.1	3.4	2.5	3.6	9.0	4.6	4.5	30.3
	1K4C	16068	118467	2.1	1.8	6.1	3.0	3.4	19.5	8.3	7.9	52.2
	2OAU	16647	199101	3.0	1.7	6.0	10.2	28.3	90.0	12.1	7.5	158.9
	1AON	58674	495683	6.3	2.0	12.4	9.9	12.5	87.9	17.9	11.0	159.9

Table 3. Table 3 : Incorrectly identified simplices of the alpha complex.

$\alpha$	PDB id	$\#$Atoms	$\#$Simplices	$\#$Misclassified Simplices				Error rate
				Edges	Triangles	Tetrahedra	Total
0.0	1GRM	260	932	0	0	0	0	0.0000
	1U71	1505	5696	0	0	0	0	0.0000
	3N0H	1509	5739	0	0	0	0	0.0000
	4HHB	4384	38796	0	0	0	0	0.0000
	2J1N	8142	29642	0	0	0	0	0.0000
	1K4C	16068	62851	15	33	16	64	0.0010
	2OAU	16647	123175	12	21	5	38	0.0003
	1AON	58674	262244	22	39	21	82	0.0003
1.0	1GRM	260	1598	0	0	0	0	0.0000
	1U71	1505	10828	0	0	0	0	0.0000
	3N0H	1509	10965	0	0	0	0	0.0000
	4HHB	4384	65987	0	0	0	0	0.0000
	2J1N	8142	58205	0	0	0	0	0.0000
	1K4C	16068	118467	20	34	14	68	0.0006
	2OAU	16647	199101	10	22	10	42	0.0002
	1AON	58674	495683	10	26	21	57	0.0001

Equations

$$\pi (p,b_i)=\Vert p-p_i{\Vert }^2-r_i^2.$$

$$K_{\alpha }=\{{\sigma }^d\in D\text{suchleavevmodenobreak that}\mathsf{Size}({\sigma }^d)\le \alpha \}.$$

Full text

Parallel Computation of Alpha Complexes for Biomolecules

Talha Bin Masood111Corresponding author

  

Tathagata Ray

  

Vijay Natarajan

1 Introduction

The alpha complex of a set of points in $\mathbb{R}^{3}$ is a subset of the Delaunay triangulation. A size parameter $\alpha$ determines the set of simplices (tetrahedra, triangles, edges, and vertices) of the Delaunay triangulation that are included in the alpha complex. It is an elegant representation of the shape of the set of points [16, 18, 14] and has found various applications, particularly in molecular modeling and molecular graphics. The atoms in a biomolecule are represented by weighted points in $\mathbb{R}^{3}$, and the region occupied by the molecule is represented by the union of balls centered at these points. The geometric shape of a biomolecule determines its function, namely how it interacts with other biomolecules. The alpha complex represents the geometric shape of the molecule very efficiently. It has been widely used for computing and studying geometric features such as cavities and channels [25, 26, 10, 30, 33, 28, 23]. Further, an alpha complex based representation is also crucial for accurate computation of geometric properties like volume and surface area [24, 12, 27].

Advances in imaging technology have resulted in a significant increase in the size of molecular structure data. This necessitates the development of efficient methods for storing, processing, and querying these structures. In this paper, we study the problem of efficient construction of the alpha complex with particular focus on point distributions that are typical of biomolecules. In particular, we present a parallel algorithm for computing the alpha complex and an efficient GPU implementation that outperforms existing methods. In contrast to existing algorithms, our algorithm does not require the explicit construction of the Delaunay triangulation.

1.1 Related work

The Delaunay triangulation has been studied within the field of computational geometry for several decades and numerous algorithms have been proposed for its construction [1]. Below, we describe only a few methods that are most relevant to this paper.

A tetrahedron belongs to the Delaunay triangulation of a set of points in $\mathbb{R}^{3}$ if and only if it satisfies the empty circumsphere property, namely no point is contained within the circumsphere of the tetrahedron. The Bowyer-Watson algorithm [4, 34] and the incremental insertion algorithm by Guibas et al. [21] are based on the above characterization of the Delaunay triangulation. In both methods, points are inserted incrementally and the triangulation is locally updated to ensure that the Delaunay property is satisfied. The incremental insertion method followed by bi-stellar flipping works in higher dimensions also [20] and can construct the Delaunay triangulation in $O(n\log n+n^{\left\lceil d/2\right\rceil})$ time in the worst case, where $n$ is the number of input points in $\mathbb{R}^{d}$. A second approach for constructing the Delaunay triangulation is based on its equivalence to the convex hull of the points lifted onto a $(d+1)$-dimensional paraboloid [19].

A third divide-and-conquer approach partitions the inputs points into two or generally multiple subsets, constructs the Delaunay triangulation for each partition, and merges the pieces of the triangulation finally. The merge procedure depends on the ability to order the edges incident on a vertex and hence works only in $\mathbb{R}^{2}$. The extension to $\mathbb{R}^{3}$ requires that the merge procedure be executed first [6]. The divide-and-conquer strategy directly extends to a parallel algorithm [31, 5]. The DeWall algorithm [6] partitions the input point set into two halves and first constructs the triangulation of points lying within the boundary region of the two partitions. The Delaunay triangulation of the two halves is then constructed in parallel. The process is repeated recursively resulting in increased parallelism. Cao et al. [5] have developed a GPU parallel algorithm, gDel3D, that constructs the Delaunay triangulation in two stages. In the first stage, points are inserted in parallel followed by flipping to obtain an approximate Delaunay triangulation. In the second stage, a star splaying procedure works locally to convert non-Delaunay tetrahedra into Delaunay tetrahedra. The algorithm can be extended to construct the weighted Delaunay triangulation for points with weights. Cao et al. report a speed up of up to a factor of $10$ over a sequential implementation for constructing the weighted Delaunay triangulation of 3 million weighted points.

Existing algorithms for constructing the alpha complex [11, 18, 27, 9] often require that the Delaunay triangulation be computed in a first step, with the exception of a recent method that guarantees output sensitive construction under mild assumptions on weights [32] or a possible construction from Čech complexes [2]. Simplices that belong to the alpha complex are identified using a size filtration in a second step. Simplices that belong to the alpha complex are identified using a size filtration in a second step. In the case of biomolecules, only small values of the size parameter are of interest and the number of simplices in the alpha complex is a small fraction of those contained in the Delaunay triangulation. Hence, the Delaunay triangulation construction is often the bottleneck in the alpha complex computation. The key difficulty lies in the absence of a direct characterization of simplices that belong to the alpha complex.

1.2 Summary of results

We propose an algorithm that avoids the expensive Delaunay triangulation computation and instead directly computes the alpha complex for biomolecules. The key contributions of this paper are summarized below:

• •

A new characterization of the alpha complex – a set of conditions necessary and sufficient for a simplex to be a part of the alpha complex.

• •

A new algorithm for computing the alpha complex of a set of weighted points in $\mathbb{R}^{3}$. The algorithm identifies simplices of the alpha complex in decreasing order of dimension without computing the complete weighted Delaunay triangulation.

• •

An efficient CUDA-based parallel implementation of this algorithm for biomolecular data that can compute the alpha complex for a 10 million point dataset in approximately 10 seconds.

• •

A proof of correctness of the algorithm and comprehensive experimental validation to demonstrate that it outperforms existing methods.

While the experimental results presented here focus on biomolecular data, the algorithm is applicable to data from other application domains as well. In particular, the efficient GPU implementation may be used for points that arise in smoothed particle hydrodynamics (SPH) simulations, atomistic simulations in material science, and particle systems that appear in computational fluid dynamics (CFD).

2 Background

In this section, we review the necessary background on Delaunay triangulations required to describe the algorithm and also establish a new characterization of the alpha complex that does not require the Delaunay triangulation. For a detailed description of Delaunay triangulations, alpha complexes, and related structures, we refer the reader to various books on the topic [1, 13, 15].

Let $B=\{b_{i}\}$ denote a set of balls or weighted points, where $b_{i}=(p_{i},r_{i})$ represents a ball centered at $p_{i}$ with radius $r_{i}$. We limit our discussion to balls in $\mathbb{R}^{3}$, so $p_{i}=(x_{i},y_{i},z_{i})\in\mathbb{R}^{3}$. Further, we assume that the points in $B$ are in general position, i.e., no two points have the same location, no three points are collinear, no four points are coplanar, and no subset of five points are equidistant from a point in $\mathbb{R}^{3}$. Such configurations are called degeneracies. In practice, a degenerate input can be handled via symbolic perturbation [17].

2.1 Simplex and simplicial complex

A $d$-dimensional simplex $\sigma^{d}$ is defined as the convex hull of $d+1$ affinely independent points. Assuming the centres of balls in $B$ are in general position, all $(d+1)$ sized subsets of $B$ form a simplex $\sigma^{d}=(p_{0}^{\sigma},p_{1}^{\sigma},\cdots,p_{d}^{\sigma})$. For simplicity, we sometimes use $b_{i}$ instead of the center $p_{i}$ to refer to points incident on a simplex. For example, we may write $\sigma^{d}=(b_{0}^{\sigma},b_{1}^{\sigma},\cdots,b_{d}^{\sigma})$.

A non-empty strict subset of $\sigma^{d}$ is also a simplex but with dimension smaller than $d$. Such a simplex is called a face of $\sigma^{d}$. Specifically, a $(d-1)$-dimensional face of $\sigma^{d}$ is referred to as a facet of $\sigma^{d}$. A set of simplices $K$ is called a simplicial complex if: 1) a simplex $\sigma\in K$ implies that all faces of $\sigma$ also belong to $K$, and 2) for two simplices $\sigma_{1}$, $\sigma_{2}\in K$, either $\sigma_{1}\cap\sigma_{2}\in K$ or $\sigma_{1}\cap\sigma_{2}=\emptyset$.

2.2 Power distance and weighted Voronoi diagram

The power distance $\pi(p,b_{i})$ between a point $p\in\mathbb{R}^{3}$ and a ball $b_{i}=(p_{i},r_{i})\in B$ is defined as

[TABLE]

The weighted Voronoi diagram is an extension of the Voronoi diagram to weighted points. It is a partition of $\mathbb{R}^{3}$ based on proximity to input balls $b_{i}$ in terms of the power distance. Points $p\in\mathbb{R}^{3}$ that are closer to the ball $b_{i}$ compared to all other balls $b_{j}\in B$ ($j\neq i$) constitute the Voronoi region of $b_{i}$. Points equidistant from two balls $b_{i}$, $b_{j}\in B$ and closer to these two balls compared to other balls constitute a Voronoi face. Similarly, points equidistant from three balls and fours balls constitute Voronoi edges and Voronoi vertices of the weighted Voronoi diagram, respectively. Figure 1(b) shows the weighted Voronoi diagram for a set of 2D weighted points or disks on the plane. Similar to the unweighted case, the Voronoi regions of the weighted Voronoi diagram are convex and linear. However, the weights may lead to a configuration where the Voronoi region of $b_{i}$ is disjoint from $b_{i}$. This occurs when $b_{i}$ is contained within another ball $b_{j}$. Further, the Voronoi region of $b_{i}$ may even be empty.

2.3 Weighted Delaunay triangulation

The weighted Delaunay triangulation is the dual of the weighted Voronoi diagram, see Figure 1(c). It is a simplicial complex consisting of simplices that are dual to the cells of the weighted Voronoi diagram. The following equivalent definition characterizes a simplex $\sigma^{d}$ belonging to a Delaunay triangulation $D$.

Definition 2.1** (Weighted Delaunay Triangulation).**

A simplex $\sigma^{d}$ = $(p_{0}^{\sigma},p_{1}^{\sigma},\cdots,p_{d}^{\sigma}),\,0\leq d\leq 3$, belongs to the weighted Delaunay triangulation $D$ of $B$ if and only if there exists a point $p\in\mathbb{R}^{3}$ such that

DT1:

$\pi(p,b_{0}^{\sigma})=\pi(p,b_{1}^{\sigma})=\cdots=\pi(p,b_{d}^{\sigma})$, and

DT2:

$\pi(p,b_{0}^{\sigma})\leq\pi(p,b_{i})$* for $b_{i}\in B-\sigma^{d}$.*

A point $p$ that satisfies the above two conditions, DT1 and DT2, is called a witness for $\sigma^{d}$. We call a point that minimizes the distance $\pi(p,b_{0}^{\sigma})$ and satisfies both conditions as the closest witness, denoted by $p_{\mathsf{min}}^{\sigma}$. This minimum distance $\pi(p_{\mathsf{min}}^{\sigma},b_{0}^{\sigma})$ is called the $\mathsf{Size}$ of the simplex $\sigma^{d}$. A point that minimizes the distance $\pi(p,b_{0}^{\sigma})$ and satisfies DT1 is called the $ortho\text{-}center$ $p_{\mathsf{ortho}}^{\sigma}$ of simplex $\sigma^{d}$. The distance $\pi(p_{\mathsf{ortho}}^{\sigma},b_{0}^{\sigma})$ is called the $\mathsf{OrthoSize}$ of the simplex $\sigma^{d}$. Clearly, the $\mathsf{Size}$ of a simplex is lower bounded by its $\mathsf{OrthoSize}$. Figure 2 shows the two possible scenarios, namely when $\mathsf{OrthoSize}=\mathsf{Size}$ and $\mathsf{OrthoSize}<\mathsf{Size}$.

2.4 Alpha complex

Given a parameter $\alpha\in\mathbb{R}$, we can construct a subset of the weighted Delaunay triangulation by filtering simplices whose $\mathsf{Size}$ is less than or equal to $\alpha$, see Figures 1(d) and 1(e). The resulting subset, called the alpha complex, is a subcomplex of the Delaunay complex and is denoted $K_{\alpha}$:

[TABLE]

The following equivalent definition characterizes simplices of the alpha complex without explicitly referring to the Delaunay triangulation.

Definition 2.2** (Alpha complex).**

A $d$-dimensional simplex $\sigma^{d}$ = $(p_{0}^{\sigma},p_{1}^{\sigma},\cdots,p_{d}^{\sigma}),0\leq d\leq 3$, belongs to the alpha complex $K_{\alpha}$ of $B$ if and only if there exists a point $p\in\mathbb{R}^{3}$ such that the following three conditions are satisfied:

AC1:

$\pi(p,b_{0}^{\sigma})=\pi(p,b_{1}^{\sigma})=\cdots=\pi(p,b_{d}^{\sigma})$,

AC2:

$\pi(p,b_{0}^{\sigma})\leq\pi(p,b_{i})$* for $b_{i}\in B-\sigma^{d}$, and*

AC3:

$\pi(p,b_{0}^{\sigma})\leq\alpha$* or equivalently, the $\mathsf{Size}$ of $\sigma^{d}$ is at most $\alpha$.*

3 Algorithm

We now describe an algorithm to compute the alpha complex and prove its correctness. The algorithm utilizes the characterizing conditions introduced above. It first identifies the tetrahedra that belong to the alpha complex, followed by the set of triangles, edges and vertices. Figure 3 illustrates the algorithm as applied to disks on the plane.

3.1 Outline

The alpha complex of a point set in $R^{3}$ consists of simplices of dimensions 0–3, $K_{\alpha}=K_{\alpha}^{0}\cup K_{\alpha}^{1}\cup K_{\alpha}^{2}\cup K_{\alpha}^{3}$, where $K_{\alpha}^{d}\subset K_{\alpha}$ is the set of $d$-dimensional simplices in $K_{\alpha}$. We initialize $K_{\alpha}^{d}=\emptyset$ and construct $K_{\alpha}$ in five steps described below:

Step 1:

For $0\leq d\leq 3$, compute the set of all simplices $\sigma^{d}$ such that $\mathsf{Ortho\mathsf{Size}}(\sigma^{d})\leq\alpha$. Let this set be denoted by $\Sigma_{\mathsf{ortho}}=\Sigma_{\mathsf{ortho}}^{0}\cup\Sigma_{\mathsf{ortho}}^{1}\cup\Sigma_{\mathsf{ortho}}^{2}\cup\Sigma_{\mathsf{ortho}}^{3}$.

Step 2:

For all tetrahedra $\sigma^{3}\in\Sigma_{\mathsf{ortho}}^{3}$, check condition AC2 using $p=p_{\mathsf{ortho}}^{\sigma}$. If $\sigma^{3}$ satisfies AC2 then insert it into $K_{\alpha}^{3}$.

Step 3:

Insert all triangles that are incident on tetrahedra in $K_{\alpha}^{3}$ into $K_{\alpha}^{2}$. Let $\Sigma_{\mathsf{free}}^{2}=\Sigma_{\mathsf{ortho}}^{2}-\mathsf{Facets}(K_{\alpha}^{3})$, where $\mathsf{Facets}(K_{\alpha}^{3})$ denotes the set of facets of tetrahedra in $K_{\alpha}^{3}$. For all triangles $\sigma^{2}\in\Sigma_{\mathsf{free}}^{2}$, check condition AC2 using $p=p_{\mathsf{ortho}}^{\sigma}$. If $\sigma^{2}$ satisfies AC2 then insert it into $K_{\alpha}^{2}$.

Step 4:

Insert all edges incident on triangles in $K_{\alpha}^{2}$ into $K_{\alpha}^{1}$. Let $\Sigma_{\mathsf{free}}^{1}=\Sigma_{\mathsf{ortho}}^{1}-\mathsf{Facets}(K_{\alpha}^{2})$, where $\mathsf{Facets}(K_{\alpha}^{2})$ denotes the set of facets of triangles in $K_{\alpha}^{2}$. For all edges $\sigma^{1}\in\Sigma_{\mathsf{free}}^{1}$, check condition AC2 using $p=p_{\mathsf{ortho}}^{\sigma}$. If $\sigma^{1}$ satisfies AC2 then insert it into $K_{\alpha}^{1}$.

Step 5:

Insert all endpoints of edges in $K_{\alpha}^{1}$ into $K_{\alpha}^{0}$. Let $\Sigma_{\mathsf{free}}^{0}=\Sigma_{\mathsf{ortho}}^{0}-\mathsf{Facets}(K_{\alpha}^{1})$, where $\mathsf{Facets}(K_{\alpha}^{1})$ denotes the set of balls incident on edges in $K_{\alpha}^{1}$. For all balls $b_{i}=(p_{i},r_{i})\in\Sigma_{\mathsf{free}}^{0}$, check condition AC2 using $p=p_{i}$. If $p_{i}$ satisfies AC2 then insert it into $K_{\alpha}^{0}$.

Step 1 selects simplices that satisfy AC3. Step 2 recognizes tetrahedra that belong to the alpha complex by checking AC2 using $p=p_{\mathsf{ortho}}^{\sigma}$. Triangle faces of these tetrahedra also belong to $K_{\alpha}$. The other “free” triangles belong to $K_{\alpha}^{2}$ if they satisfy AC2. Step 4 identify edges similarly. First all edge faces of triangles in $K_{\alpha}^{2}$ are inserted followed by those “free” edges that satisfy AC2. Vertices are identified similarly in Step 5.

A notion related to free simplices, called unattached simplices, was introduced by Edelsbrunner [11]. However, the characterization of unattached simplices depends on the fact that they belong to the Delaunay complex.

3.2 Proof of correctness

We now prove that that the algorithm described above correctly computes the alpha complex of the given set of weighted points by proving the following four claims. Each claim states that the set of simplices computed in Steps 2, 3, 4 and 5 are exactly the simplices belonging to the alpha complex. We assume that the input is non-degenerate.

Claim 1**.**

Step 2 computes $K_{\alpha}^{3}$ correctly.

Proof 3.1**.**

For a tetrahedron $\sigma^{3}$, $p_{\mathsf{ortho}}^{\sigma}$ is the only point that satisfies condition AC1. In Step 2 of the proposed algorithm, we check if AC2 holds for $p_{\mathsf{ortho}}^{\sigma}$. If yes, then $p_{\mathsf{ortho}}^{\sigma}$ is a witness for $\sigma^{3}$, i.e., $p_{\mathsf{ortho}}^{\sigma}=p_{\mathsf{min}}^{\sigma}$. Further, since $\mathsf{OrthoSize}(\sigma^{3})\leq\alpha$ and $p_{\mathsf{ortho}}^{\sigma}=p_{\mathsf{min}}^{\sigma}$, we have $\mathsf{Size}(\sigma^{3})\leq\alpha$ thereby satisfying AC3. Therefore, $\sigma^{3}$ belongs to $K_{\alpha}^{3}$ because it satisfies all three conditions.

We now prove that the algorithm correctly identifies the triangles of the alpha complex.

Lemma 3.2**.**

A triangle $\sigma^{2}\in\Sigma_{\mathsf{free}}^{2}$ belongs to $K_{\alpha}^{2}$ if and only if it satisfies AC2 with $p=p_{\mathsf{ortho}}^{\sigma}$.

Proof 3.3**.**

We first prove the backward implication, namely if $\sigma^{2}\in\Sigma_{\mathsf{free}}^{2}$ satisfies AC2 with $p=p_{\mathsf{ortho}}^{\sigma}$, then $\sigma^{2}\in K_{\alpha}^{2}$. Note that $p_{\mathsf{ortho}}^{\sigma}$ satisfies AC1 by definition. Further, it satisfies AC2 by assumption and hence $\mathsf{Size}(\sigma^{2})=\mathsf{OrthoSize}(\sigma^{2})$. We also have $\mathsf{OrthoSize}(\sigma^{2})\leq\alpha$ because $\sigma^{1}\in\Sigma_{\mathsf{free}}^{2}\subseteq\Sigma_{\mathsf{ortho}}^{2}$. So, $\mathsf{Size}(\sigma^{2})\leq\alpha$ thereby satisfying AC3. The triangle $\sigma^{2}$ with $p=p_{\mathsf{ortho}}^{\sigma}$ satisfies all three conditions and hence belongs to $K_{\alpha}^{2}$.

We will now prove the forward implication via contradiction. Suppose there exists a triangle $\sigma^{2}\in\Sigma_{\mathsf{free}}^{2}$ that belongs to $K_{\alpha}^{2}$ but does not satisfy AC2 with $p=p_{\mathsf{ortho}}^{\sigma}$. In other words, there exists a ball $b_{i}\in B-\sigma^{2}$ for which $\pi(p_{\mathsf{ortho}}^{\sigma},b_{i})<\pi(p_{\mathsf{ortho}}^{\sigma},b_{0}^{\sigma})$. Let $B_{v}$ denote the set of all such balls $b_{i}$. The set of points that are equidistant from the three balls $(b_{0}^{\sigma},b_{1}^{\sigma},b_{2}^{\sigma})$ corresponding to $\sigma^{2}$ form a line perpendicular to the plane containing $\sigma^{2}$ called the radical axis. Each ball $b_{i}\in B_{v}$ partitions the radical axis into two half-intervals based on whether the point on radical axis is closer to $b_{i}$ or to $b_{0}^{\sigma}$, see Figure 4. Let $I^{+}(b_{i})$ denote the half interval consisting of points that are closer to $b_{0}^{\sigma}$ compared to $b_{i}$. Let $I^{+}(B_{v})$ denote the intersection of all such half intervals $I^{+}(b_{i})$. We have assumed that $\sigma^{2}\in K_{\alpha}^{2}$, so there must exist a closest witness $p_{\mathsf{min}}^{\sigma}$, and it lies within $I^{+}(B_{v})$. Thus, $I^{+}(B_{v})$ is non-empty. In fact, $I^{+}(B_{v})=I^{+}(b_{j})$ for some $b_{j}\in B_{v}$ and $p_{\mathsf{min}}^{\sigma}$ is exactly the end point of $I^{+}(b_{j})$. This implies that $p_{\mathsf{min}}^{\sigma}$ is also a closest witness for the tetrahedron $\sigma^{3}=(b_{0}^{\sigma},b_{1}^{\sigma},b_{2}^{\sigma},b_{j})$. So, $\sigma^{3}$ belongs to $K_{\alpha}^{3}$ and its $\mathsf{Size}$ is equal to $\mathsf{Size}(\sigma^{2})$. However, this means that $\sigma^{2}\notin\Sigma_{\mathsf{free}}^{2}$, a contradiction. So, the forward implication in the lemma is true.

Claim 2**.**

Step 3 computes $K_{\alpha}^{2}$ correctly.

Proof 3.4**.**

If a simplex $\sigma^{3}$ belongs to $K_{\alpha}$ then naturally all of its faces also belong to $K_{\alpha}$. The algorithm includes such triangles into $K_{\alpha}^{2}$ and remove them from $\Sigma_{\mathsf{ortho}}^{2}$ to obtain the set of free triangles $\Sigma_{\mathsf{free}}^{2}$ . It follows directly from Lemma 3.2 that AC2 is a necessary and sufficient condition for a triangle in $\Sigma_{\mathsf{free}}^{2}$ to belong to $K_{\alpha}^{2}$. Hence, Step 3 correctly computes the triangles belonging to $K_{\alpha}^{2}$.

The above arguments need to be extended to prove that the edges of the alpha complex are also correctly identified.

Lemma 3.5**.**

An edge $\sigma^{1}\in\Sigma_{\mathsf{free}}^{1}$ belongs to $K_{\alpha}^{1}$ if and only if it satisfies the condition AC2 with $p=p_{\mathsf{ortho}}^{\sigma}$.

Proof 3.6**.**

First, we assume that $\sigma^{1}\in\Sigma_{\mathsf{free}}^{1}$ satisfies AC2 with $p=p_{\mathsf{ortho}}^{\sigma}$. The point $p_{\mathsf{ortho}}^{\sigma}$ satisfies AC1 by definition. Further, it satisfies AC2 by assumption and hence $\mathsf{Size}(\sigma^{1})=\mathsf{OrthoSize}(\sigma^{1})$. We also have $\mathsf{OrthoSize}(\sigma^{1})\leq\alpha$ because $\sigma^{1}\in\Sigma_{\mathsf{free}}^{1}\subseteq\Sigma_{\mathsf{ortho}}^{1}$. So, $\mathsf{Size}(\sigma^{1})\leq\alpha$ thereby satisfying AC3. The edge $\sigma^{1}$ with $p=p_{\mathsf{ortho}}^{\sigma}$ satisfies all three conditions and hence belongs to $K_{\alpha}^{1}$.

We will prove the forward implication via contradiction. Suppose there exists an edge $\sigma^{1}\in\Sigma_{\mathsf{free}}^{1}$ that belongs to $K_{\alpha}^{1}$ but does not satisfy AC2 with $p=p_{\mathsf{ortho}}^{\sigma}$. In other words, there exists a ball $b_{i}\in B-\sigma^{1}$ such that $\pi(p_{\mathsf{ortho}}^{\sigma},b_{i})<\pi(p_{\mathsf{ortho}}^{\sigma},b_{0}^{\sigma})$. Let $B_{v}$ denote the set of all such balls $b_{i}$. The set of points that are equidistant from the two balls $(b_{0}^{\sigma},b_{1}^{\sigma})$ corresponding to $\sigma^{1}$ form a plane perpendicular to the line containing $\sigma^{1}$ called the radical plane. Each ball $b_{i}\in B_{v}$ partitions the radical plane into two half-planes based on whether the point on the radical plane is closer to $b_{i}$ or to $b_{0}^{\sigma}$, see Figure 5. Let $H^{+}(b_{i})$ denote the half-plane consisting of points that are closer to $b_{0}^{\sigma}$ compared to $b_{i}$. Let $H^{+}(B_{v})$ denote the intersection of all such half-planes $H^{+}(b_{i})$. We have assumed that $\sigma^{1}\in K_{\alpha}^{1}$, so there must exist a closest witness $p_{\mathsf{min}}^{\sigma}$, and it lies within $H^{+}(B_{v})$. Thus, $H^{+}(B_{v})$ is non-empty. In fact, $p_{\mathsf{min}}^{\sigma}$ lies on the envelope of $H^{+}(B_{v})$ because it minimizes the distance to $b_{0}^{\sigma}$. Let $p_{\mathsf{min}}^{\sigma}$ lie on the bounding line corresponding to $H^{+}(b_{j})$ for some $b_{j}\in B_{v}$. This implies that $p_{\mathsf{min}}^{\sigma}$ is also a closest witness for the triangle $\sigma^{2}=(b_{0}^{\sigma},b_{1}^{\sigma},b_{j})$. So, $\sigma^{2}$ belongs to $K_{\alpha}^{2}$, and its size is equal to $\mathsf{Size}(\sigma^{1})$. However, this means that $\sigma^{1}\notin\Sigma_{\mathsf{free}}^{1}$, a contradiction. So, the forward implication in the lemma is also true.

Claim 3**.**

Step 4 computes $K_{\alpha}^{1}$ correctly.

Proof 3.7**.**

All edge faces of triangles in $K_{\alpha}^{2}$ naturally belong to $K_{\alpha}^{1}$. Step 4 inserts all edges incident on triangles in $K_{\alpha}^{2}$ into $K_{\alpha}^{1}$ as valid edges and removes them from $\Sigma_{\mathsf{ortho}}^{1}$ to obtain the set of free edges $\Sigma_{\mathsf{free}}^{1}$. It follows directly from Lemma 3.5 that AC2 is a necessary and sufficient condition for an edge $\sigma^{1}\in\Sigma_{\mathsf{free}}^{1}$ to belong to $K_{\alpha}^{1}$. Therefore, Step 4 correctly computes the edges belonging to $K_{\alpha}^{1}$.

Claim 4**.**

Step 5 computes $K_{\alpha}^{0}$ correctly.

Proof 3.8**.**

All vertices incident on $K_{\alpha}^{1}$ naturally belong to $K_{\alpha}^{0}$. Step 5 inserts all such vertices in $K_{\alpha}^{0}$ as valid vertices and removes them from $\Sigma_{\mathsf{ortho}}^{0}$ to obtain the set of free vertices $\Sigma_{\mathsf{free}}^{0}$. Next, the vertices in $\Sigma_{\mathsf{free}}^{0}$ for which the center of the ball $b_{i}=(p_{i},r_{i})$ satisfies AC2 are also inserted into $K_{\alpha}^{0}$. Clearly, these vertices also satisfy AC3 because they belong to $\Sigma_{\mathsf{ortho}}^{0}$. The condition AC1 is not relevant for [math]-dimensional simplices. Therefore, these vertices clearly belong to the alpha complex. Similar to Lemmas 3.2 and 3.5, it is easy to prove that checking for AC2 for $p=p_{i}$ is necessary and sufficient condition to decide whether a vertex in $\Sigma_{\mathsf{free}}^{0}$ belongs to the alpha complex. That is, it is possible to show that vertices in alpha complex that have non-empty Voronoi regions but do not satisfy AC2 for $p=p_{i}$ would be incident on some edge in $K_{\alpha}^{1}$, and therefore must have been already detected by Step 4 and hence can not belong to $\Sigma_{\mathsf{free}}^{0}$. Therefore, Step 5 correctly computes the vertices belonging to $K_{\alpha}^{0}$.

The arguments in the proofs of Lemma 3.2 and Lemma 3.5 are similar. A general result is likely true for $d$-dimensional simplices in $\Sigma_{\mathsf{free}}^{d}$. However, given the focus on alpha complexes in $\mathbb{R}^{3}$, we prefer to state and prove these results specific to lower dimensions. We also prefer to provide individual proofs for edges and triangles because it simplifies the exposition and could also potentially help in the design of improved data structures to accelerate computation of different steps of the algorithm.

4 Parallel algorithm for biomolecules

Although the algorithm as described above is provably correct, a straightforward implementation will be extremely inefficient with a worst-case running time of $O(n^{5})$, where $n$ is the number of weighted points in $B$. This is because Step 1 requires $O(n^{4})$ time to generate all possible tetrahedra. In later steps, we need $O(n)$ effort per simplex to check AC2. However, the input corresponds to atoms in a biomolecule. We show how certain properties of biomolecules can be leveraged to develop a fast parallel implementation.

4.1 Biomolecular data characteristics

Atoms in a biomolecule are well distributed. The following three properties of biomolecules are most relevant:

• •

The radius of an atom is bounded. The typical radius of an atom in a protein molecule ranges between 1Å to 2Å [3]. Further, a protein molecule contains upwards of thousand atoms. So, the radius is small compared to the total size of the molecule.

• •

There is a lower bound on the distance between the centres of two atoms. This is called the van der Waals contact distance, beyond which the two atoms start repelling each other. In the case of atoms in protein molecules, this distance is at least 1Å. This property together with the upper bound on atomic radii ensures that no atom is completely contained inside another. This means that the weighted Voronoi regions corresponding to the atoms in a biomolecule can be always be assumed to be non-empty.

• •

Structural biologists are interested in small values of $\alpha$. The two crucial values are 0Å and 1.4Å. The former corresponds to using van der Waals radius and the latter corresponds to the radius of water molecule, which acts as the solvent.

In the light of the above three properties, we can say that the number of simplices of the alpha complex that are incident on a weighted point (atom) is independent of the total number of input atoms and hence bounded by a constant [22].

4.2 Acceleration data structure

The algorithm will benefit from an efficient method for accessing points of $B$ that belong to a local neighborhood of a given weighted point. We store the weighted points in a grid-based data structure. Let $r_{max}$ denote the radius of the largest atom and assume that the value of the parameter $\alpha$ is available as input. First, we construct a grid with cells of side length $\sqrt{r_{max}^{2}+\alpha}$ and then bin the input atoms into the grid cells. In our implementation, we do not store the grid explicitly because it may contain several empty cells. Instead, we compute the cell index for each input atom and sort the list of atoms by cell index to ensure that atoms that belong to a particular cell are stored at consecutive locations. The cell index is determined based on a row-major or column-major order. Alternatively, a space-filling curves like the Hilbert curve could also be used to order the cells.

After the atoms are stored in grid cells, the alpha complex is computed in two stages. In the first stage, we employ a bottom-up approach to obtain a conservative estimate of the edges, triangles, and tetrahedra belonging to the alpha complex. The false positives from the first stage are removed in a subsequent pruning stage resulting in the correct alpha complex. We describe these two stages in the following subsections.

4.3 Potential simplices

The first stage essentially corresponds to Step 1 of the algorithm described in the previous section. We compute the set $\Sigma_{\mathsf{ortho}}$ of potential simplices for which $\mathsf{OrthoSize}(\sigma^{d})\leq\alpha$. However, for efficiency reasons we process the simplices in the order of increasing dimension. First, we identify edges that satisfy the AC3 condition as described below. Given the size of the grid cell, endpoints of edges that satisfy the condition either lie within the same grid cell or in adjacent cells. So, the grid data structure substantially reduces the time required to compute the list of potential edges $\Sigma_{\mathsf{ortho}}^{1}$. Beginning from this set of edges, we construct the set of all possible triangles and retain the triangles whose $\mathsf{OrthoSize}$ is no greater than $\alpha$, resulting in the set $\Sigma_{\mathsf{ortho}}^{2}$. Finally, we use the triangles in $\Sigma_{\mathsf{ortho}}^{2}$ to construct the list of tetrahedra that satisfy the $\mathsf{OrthoSize}\leq\alpha$ condition. The above procedure works because the $\mathsf{OrthoSize}$ of a simplex is always greater than or equal to the $\mathsf{OrthoSize}$ of its faces. The set of simplices identified in this stage contains all simplices of the alpha complex. False positives are pruned in the second stage described below.

4.4 Pruning

The second stage corresponds to Steps 2-5 of the algorithm and processes the potential simplices in the decreasing order of dimension. This stage checks the characterizing condition AC2 to prune $\Sigma_{\mathsf{ortho}}$ into $K_{\alpha}$. The tetrahedra are processed by checking if any of the input balls are closer to the $ortho\text{-}center$ than the balls incident on the tetrahedron. If yes, the tetrahedron is pruned away. Else, the tetrahedron is recognized as belonging to the alpha complex and inserted into $K_{\alpha}^{3}$. Triangles incident on these tetrahedra also belong to the alpha complex and are inserted into $K_{\alpha}^{2}$ after they are removed from the list of potential triangles $\Sigma_{\mathsf{ortho}}^{2}$. Next, the triangles in $\Sigma_{\mathsf{ortho}}^{2}$ are processed by checking if they satisfy AC2. If yes, they are inserted into $K_{\alpha}^{2}$. Otherwise, they are pruned away. All edges incident on triangles belonging to $K_{\alpha}^{2}$ are inserted into $K_{\alpha}^{1}$ and removed from the set $\Sigma_{\mathsf{ortho}}^{1}$. Next, the edges in $\Sigma_{\mathsf{ortho}}^{1}$ are processed by checking if they satisfy AC2. Edges that satisfy AC2 are inserted into $K_{\alpha}^{1}$ and the others are pruned away. All the vertices in $\Sigma_{\mathsf{ortho}}^{0}$ are directly inserted into $K_{\alpha}^{0}$ without the AC2 check because for biomolecular data we assume that Voronoi regions of all the atoms are non-empty. The check for condition AC2 for each simplex is again made efficient by the use of the grid data structure. Atoms that may violate AC2 lie within the same cell as that containing the $ortho\text{-}center$ or within the adjacent cells. Atoms that lie within other cells may be safely ignored.

4.5 CUDA implementation

We use the CUDA framework [7] and the thrust library [8] within CUDA to develop a parallel implementation of the algorithm that executes on the many cores of the GPU. The grid computation is implemented as a CUDA kernel where all atoms are processed in parallel. The computation of potential simplices and pruning stages are broken down into multiple CUDA kernels and parallelized differently in order to increase efficiency. We now describe the parallelization strategy in brief.

For computing the set of potential edges, the initial enumeration of possible edges incident on an atom is done using the atoms in the corresponding grid cell and its neighbouring cells. This is done per atom in parallel, the thread corresponding to the atom $i$ being responsible for generating the edges $ij,j>i$ to ensure no duplicate edges are generated. Subsequently, the AC3 condition is checked for the edges in parallel to finally generate the list of potential edges $\Sigma_{\mathsf{ortho}}^{1}$. For computing potential triangles $\Sigma_{\mathsf{ortho}}^{2}$, the potential edge list is used as a starting point for the initial enumeration of all possible triangles. This step is also parallelized per atom, the thread $i$ being responsible for generation of triangles of the type $ijk;j,k>i$ if all three edges $ij$, $ik$ and $jk$ are potential edges. The AC3 condition for the triangle is checked next within a separate kernel and parallelized per triangle to generate the potential triangles $\Sigma_{\mathsf{ortho}}^{2}$. A similar strategy is used for computing the set of potential tetrahedra $\Sigma_{\mathsf{ortho}}^{3}$.

The pruning stage is parallelized per tetrahedron, triangle, and edge as required. So, computation of $p_{\mathsf{ortho}}^{\sigma}$ is done in parallel. The grid data structure is again useful in checking for potential violators of the AC2 condition. Only the balls belonging to the grid cell corresponding to $p_{\mathsf{ortho}}^{\sigma}$ or the those in neighbouring grid cells can violate the AC2 condition.

4.6 Handling large data sizes

Typical protein structures consist of up to 100,000 atoms. Our implementation can handle datasets of this size easily for reasonable values of $\alpha$. However, the size of datasets is ever increasing. Protein complexes that are available nowadays may consist of millions of atoms, necessitating smart management of GPU memory while handling such data sets.

We propose two strategies and implement one of them. The first strategy is to partition the grid by constructing an octree data structure and choosing an appropriate level in the octree to create partitions. Each partition together with its border cells can be processed independently of other partitions. So, we can copy one partition and its border to the GPU memory, compute its alpha complex, and copy the results back from GPU to CPU memory. After all the partitions are processed, the list of simplices can be concatenated followed by duplicate removal to generate the final alpha complex.

The second strategy is to partition the sorted list of atoms into chunks of equal sizes and to process each chunk independently. Here, we assume that the complete list of atoms together with the grid data structure fits in the GPU memory. This is a reasonable assumption considering that datasets containing several million atoms can easily fit on modern GPUs, which typically have at least 2GB video memory. Also, the main difficulty in handling large protein structures is managing the large lists of simplices generated within the intermediate steps of the algorithm, when compared to handling the input list of atoms or the output list of simplices. We compute the alpha complex by executing the algorithm in multiple passes. Each pass computes the alpha complex for a single chunk and copies it back to the CPU memory. We have implemented this second strategy and can handle data sizes of up to 16 million atoms on a GPU with 2GB of memory. Results are reported in the next section.

5 Experimental results

We now present results of computational experiments, which demonstrate that the parallel algorithm is fast in practice and significantly better than the state-of-the-art. We also performed runtime profiling to better understand the bottlenecks and effect of the parameter $\alpha$ on the runtime. We present results for $\alpha$ in the range $0.0$ to $2.0$. This range is important for structural analysis of biomolecules as it corresponds to solvent accessible surface of the biomolecule for typical solvent molecules like water (van der Waals radius = 1.4Å). The value $\alpha=0$ corresponds to the van der Waals surface of the biomolecule. All experiments, unless stated otherwise, were performed on a Linux system with an nVidia GTX 660 Ti graphics card running CUDA 8.0 and a 2.0GHz Intel Xeon octa core processor with 16 GB of main memory. The default number of threads per block was set at $512$ for all the CUDA kernels.

Mach and Koehl describe two techniques for computing alpha complex of biomolecules called AlphaVol and UnionBall in their paper [27]. Both approaches construct the weighted Delaunay triangulation of input atoms first followed by a filtering step to obtain the alpha complex.UnionBall is the state-of-the-art technique for alpha complex computation for biomolecules on multi-core CPU. It uses heuristics and optimizations specific to biomolecular data to improve upon AlphaVol. For biomolecules containing 5 million atoms, AlphaVol takes approximately 8600 seconds for computing the alpha complex, while UnionBall takes approximately 150 seconds. Our method computes the alpha complex in less than 3 seconds for similar sized data, see Table 1.

5.1 Comparison with gReg3D

We are not aware of any available software that can compute the alpha complex directly without first constructing the complete Delaunay triangulation. In order to compare the performance, we chose the state-of-the-art parallel algorithm for computing the weighted Delaunay triangulation in 3D, gReg3D [5]. The CUDA implementation of gReg3D is available in the public domain. Table 1 compares the running times of our proposed algorithm with that of gReg3D for twelve different biomolecules at $\alpha=0$ and $\alpha=1$. As evident from the table, we consistently observe significant speedup over gReg3D. The observed speedup is as high as a factor of 22 for the biomolecule 1X9P at $\alpha=0$, one of the largest molecules in our dataset. Clearly, the speedup goes down for $\alpha=1$ when compared to $\alpha=0$ because of the increased number of simplices in the output alpha complex. We also report the number of simplices in the alpha complex compared to the total number of simplices in the Delaunay triangulation under the column ‘$\%$Simplex’. This makes it clear why the speedup decreases as $\alpha$ is increased from 0 to 1. For example, for the protein 1AON, the fraction of alpha complex simplices increases from $15.9\%$ to $30\%$ as $\alpha$ is increased from 0 to 1. Correspondingly, the speedup decreases from a factor of 13.5 to 5.5. We repeated the experiment on a MS Windows system with an nVidia GTX 980 Ti card running CUDA 8.0 and observed similar speedups. However, the individual runtimes both for our algorithm and for gReg3D were higher on the GTX 980 Ti.

The starred entries in Table 1 are results for execution using the data partitioning approach. This is necessitated because these four large molecules generate large intermediate simplex lists that can not fit into the GPU memory if all the atoms in the molecule are processed at once. We observe that gReg3D is able to successfully compute the Delaunay complex for only one out of these four large molecules and runs out of GPU memory for the remaining three molecules.

5.2 Runtime profiling

The two stages of our parallel algorithm (potential simplices and pruning) are further divided into three steps each, corresponding to the computation of edges, triangles, and tetrahedra respectively. A grid computation step precedes the two stages. We study the computation effort for each of these seven steps of the algorithm. We also report the time spent in memory transfers from CPU to GPU and vice-versa. Thus, we report the split up of the total runtime into eight categories namely, ‘Memory transfer’, ‘Grid computation’, ‘Potential edges’, ‘Potential triangles’, ‘Potential tetrahedra’, ‘AC2 tetrahedra’, ‘AC2 triangles’ and ‘AC2 edges’.

Table 2 summarizes the observed split up of runtime for eight different biomolecules. Figure 6 shows the actual time spent for different steps and Figure 7 shows relative time spent for each step. From these figures, it is clear that the pruning stage consumes the maximum amount of time. The pruning stage involves checking the neighboring balls for violations of the AC2 condition for each simplex. Specially, the tetrahedra-pruning step (red) takes approximately $25\%$ of the total time required for alpha complex computation.

We performed additional experiments to determine the average split up over multiple runs. We computed the relative time spent for each step for different values of $\alpha$ between 0.0 and 2.0. These observations are reported in Figure 8. It is clear that the memory transfers and grid computation combined do not take more than $10\%$ of the total time. The potential-simplices-estimation stage consumes $30\%$ of the time. However, the pruning stage is most expensive, taking up roughly $60\%$ of the computation time. The pruning-tetrahedra step takes up $35\%$ of the time on average. This suggests that this step should be the focus of the optimization efforts in future. It should be noted that proportion of time spent for each step depends on the distribution of atoms in the biomolecule as well as the value of $\alpha$. This explains the significant deviation from the averages as shown by the error bars.

5.3 Effect of the value of alpha

We also performed experiments to observe the runtime performance as the value of $\alpha$ is varied between 0.0 and 2.0. Figures 9 and 10 show the results for the proteins 1K4C and 1AON, respectively. We also show how the number of simplices in the computed alpha complex increases as the value of $\alpha$ is increased. The runtime and total number of simplices follow a near-linear trend. However, increase in time required for pruning, especially for pruning the tetrahedra, is greater than time required for other steps of the algorithm. Note that although both graphs appear linear, this is not guaranteed behavior for other input. The scaling behavior depends on the distribution of the atoms in the molecule and on the range of $\alpha$ values for which the experiment is conducted.

5.4 Numerical issues

The proposed algorithm requires computation of $\mathsf{OrthoSize}$ for each simplex, which in turn requires solving systems of linear equations. These computations require higher precision than is available on the GPU. So, the results may contain numerical errors. These numerical errors ultimately manifest as misclassification of a simplex as belonging to $K_{\alpha}$ or not. We performed extensive experimentation and observed that the alpha complex was correctly computed in several cases. In cases where the results were not correct, the number of false positives and negatives (extra or missing simplices) is extremely small as compared to the total number of simplices in the alpha complex. We observed a worst case error rate of $0.001$ in our experiments, see Table 3. This error rate is tolerable for several applications. If exact computation is required, we could use a tolerance threshold to tag some simplices as requiring further checks, which can in turn be performed on the CPU or GPU using a multi-precision library [joldes2016campary]. The use of multi-precision libraries will adversely affect the computation time. Our implementation can be extended to use such a hybrid strategy. However, it will require additional experimentation to identify appropriate tolerance thresholds and performance optimization. We plan to implement this in future.

6 Conclusions

We proposed a novel parallel algorithm to compute the alpha complex for biomolecular data that does not require prior computation of the complete Delaunay triangulation. The useful characterization of simplices that belong to the alpha complex may be of independent interest. The algorithm was implemented using CUDA, which exploits the characteristics of the atom distribution in biomolecules to achieve speedups of up to a factor of 22 compared to the state-of-the-art parallel algorithm for computing the weighted Delaunay triangulation, and up to a factor of 50 speedup over the state-of-the-art implementation that is optimized for biomolecules. In future work, we plan to further improve the runtime efficiency of the parallel implementation and to resolve the numerical issues using real arithmetic.

Applications of alpha complex outside the domain of biomolecular analysis often require the complete filtration of Delaunay complex. The algorithm as presented here is not best suited for such cases. However, the algorithm may be modified to utilize a previously computed alpha complex to efficiently compute the alpha complex for higher values of $\alpha$. We plan to investigate this extension in future work.