Persistent de Rham-Hodge Laplacians in Eulerian representation for manifold topological learning

1.   Introduction

Recent years have witnessed a fast growth of topological data analysis (TDA) in data science and engineering ^[63]. The growth is driven by the great promise of topological approaches to real-world data that are distinguished from any other statistical, mathematical, physical, and engineering methods ^[11,42]. Typically, TDA offers a multi-scale topological characterization of data, which is the case with persistent homology ^[23,70], a key method employed in TDA. A major feature of persistent homology is its multi-scale analysis, which creates a family of topological spaces from the original data to track the topological persistence, i.e., the lifespan of topological invariants across scales ^[5,27]. The other major feature of persistent homology is its topological description of a space (like connected components, loops, and voids) in terms of topological invariants, such as Betti numbers. As such, persistent homology-based TDA leads to much topological simplification of the geometric information in the data ^[1,20]. Consequently, TDA typically works extremely well for data with intricate complexity ^[59,67]. Unfortunately, for data without geometric complexes, TDA may give rise to an oversimplification of key geometric characteristics, leading to a less competitive approach.

For many years, persistent homology has been used in qualitative analysis, which is somewhat counterintuitive and unproductive for nonexperts. The power of persistent homology was not demonstrated until it was utilized in quantitative and predictive analysis via machine learning algorithms ^[8,39]. Topological deep learning (TDL), coined in 2017 ^[9], was introduced to deal with large and intrinsically complex datasets using both persistent homology and deep neural networks. More recently, simplicial neural networks and other topological neural techniques have been applied in TDL to the design of neural network architecture. TDL has become an emerging paradigm in data science and machine learning ^[49]. However, an increasing concern associated with this rising popularity is whether TDL brings any practical benefit beyond its mathematical elegance. There are many applications where TDL has demonstrated superiority to other competitive methods ^[44]. Perhaps some of the most compelling examples are TDL's dominant wining of D3R Grand Challenges, an annual worldwide competition series in computer-aided drug design ^[45,46], its discovery of the mechanisms of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) evolution ^[16,60], and its successful forecast of emerging dominant SARS-CoV-2 variants BA.2 ^[17] and BA.4/BA.5 about two months in advance ^[15].

It is interesting to understand why TDL (or TDA) was so successful in the aforementioned examples, but was not competitive in many other situations in the literature ^[50]. First, biomolecular data, which is intricately complex in their internal structures^[67], was involved in the above compelling examples. As such, topological simplification was a productive process, whereas TDL leads to the severe loss of crucial geometric information in many other data that is relatively simple in their internal structures. Additionally, it was element-specific persistent homology, rather than the plain persistent homology, that was applied in the above examples. This approach captures physical and biological interactions in the biomolecular data ^[9]. In fact, in the forecast of emerging dominant SARS-CoV-2 variants BA.4/BA.5, persistent Laplacian, rather than persistent homology, was utilized. This happens because persistent homology has many drawbacks or limitations ^[64]. First, the topological invariant extracted from persistent homology is qualitative, rather than quantitative. For example, the barcode from persistent homology does not distinguish a five-number from a six-number ring. Additionally, persistent homology is incapable of dealing with different elements in a point cloud, which is ineffective with the physics and chemistry of (bio)molecular data. Moreover, persistent homology cannot describe non-topological changes, i.e., homotopic shape evolution during the multi-scale (or filtration) analysis. Further, persistent homology is incapable of handling directed networks and digraphs, such as polarization, regulation, and control issues in applications. Finally, persistent homology is unable to characterize structured data, e.g., hypergraphs, directed networks, etc. These challenges call for innovative new topological methods.

To address these challenges, the persistent spectral graph, also known as persistent combinatorial Laplacian or persistent Laplacian (PL), was introduced in 2019 ^[61]. The harmonic spectra of PLs fully recover the topological invariants of persistent homology. However, the nonharmonic spectra of PLs capture the homotopic shape evolution during the multi-scale analysis that cannot be observed with persistent homology. Computational algorithms ^[22,62] and mathematical analysis ^[33,40] of PLs have been reported. In the past few years, much effort has been given to extend persistent Laplacian to further address other limitations of persistent homology ^[28], leading to persistent sheaf Laplacians ^[65], persistent path Laplacians, persistent hypergraph and hyperdigraph Laplacians ^[37], persistent directed flag Laplacians, persistent Mayer Laplacians, and persistent interaction Laplacians ^[64]. PLs have been shown to outperform persistent homology in many applications ^[15,41].

However, defined on point cloud data, neither persistent homology nor PL can directly deal with two other commonly occurring data formats, namely, data on manifolds ^[18], such as electron density ^[68], cryogenic electron microscopy density, and computed tomography images ^[14], and curves embedded in the three-dimensional Euclidean space, such as knots, links, and tangles, and their generalizations ^[30,48]. Multi-scale Gauss link integral ^[54] and evolutionary Khovanov homology have been proposed to deal with embedded curve data ^[55]. Evolutionary Khovanov homology integrates algebraic topology, geometric topology, and metric analysis for the first time. However, effective computational algorithms are needed for this approach to be widely used in practical applications.

To carry out manifold topological analysis of data on manifolds, the evolutionary de Rham-Hodge method was introduced ^[18]. This approach creates a family of multi-scale manifolds with boundaries from a given data and then builds evolutionary Hodge Laplacian operators on the multi-scale manifolds with appropriate boundary conditions. While originated from sharply different topological spaces, evolutionary Hodge Laplacian and PLs share the same algebraic structure and capture topological invariants in their harmonic spectra ^[52]. Case studies have been given to demonstrate evolutionary de Rham-Hodge theory-based manifold topological analysis of data on manifolds ^[18]. However, this approach was based on discrete exterior calculus ^[19,21] or finite element exterior calculus ^[3] in the Lagrangian representation, which is not efficient for multi-scale analysis and machine learning studies. Specifically, the regeneration of the evolving manifolds at different scales with different Lagrangian meshes causes numerical inconsistencies and becomes expensive for practical applications in machine learning studies. This challenge calls for new effective manifold topological analysis approaches for data on manifolds.

The objective of this work is to develop a persistent de Rham-Hodge theory on the Euler representation for manifold topological learning (MTL). To this end, we solve Hodge Laplacians on a pre-designed structure-persevering Cartesian grid for all scales to avoid numerical inconsistency. We construct a natural mapping of differential forms from a manifold with boundary embedded in $\mathbb{R}^3$ to a large manifold, use it to produce persistent cohomology mapping, and construct a persistent Hodge Laplacian with built-in boundary conditions. Our new approach draws on differential geometry, algebraic topology, partial differential equations, metric analysis, and numerical analysis. To give a proof-of-principle demonstration, we pair the proposed persistent de Rham-Hodge Laplacians with an effective machine learning algorithm to predict protein-ligand binding affinities. Based on two benchmark datasets in the Protein Data Bank (PDB), PDBbind v2007 and PDBbind v2016, we show that our MTL model gives rise to cutting-edge performance.

The rest of this paper is organized as follows: Section 2 offers a primer on the de Rham-Hodge theory on manifolds with boundaries; Section 3 presents our discretization for evolutionary de Rham-Hodge theory based on spectrum calculation of Laplacians associated with sublevel sets on Cartesian grids; Section 4 presents our construction for persistent de Rham-Hodge Laplacians both in the continuous setting and for given level set functions on Cartesian grids; Section 5 showcases preliminary studies on the applications of MTL; and Section 6 concludes the paper.

2.   De Rham-Hodge theory

The de Rham-Hodge theory is an advanced mathematical framework that merges ideas from differential geometry, algebraic topology, analysis, and partial differential equations to study the properties of differential forms on smooth manifolds. It plays a crucial role in understanding the topology and geometry of manifolds through differential forms. The de Rham-Hodge theory consists of de Rham cohomology and Hodge theory. The former concerns differential forms, exterior derivative, and cohomology groups, while the latter deals with Riemannian manifolds, Hodge star operator, Hodge Laplacian, and Hodge decomposition.

Let $M$ be an $m$-dimensional smooth, orientable, compact Riemannian manifold with boundary. Denote by $\Omega^k(M)$ the space of all differential $k$-forms on $M$, i.e., the space of all smooth antisymmetric covariant tensor fields on $M$ of degree $k$. The differential $d$, also called exterior derivative, is the unique $\mathbb{R}$-linear mapping from the space of $k$-forms $\Omega^k(M)$ to the space of $(k\!+\!1)$-forms $\Omega^{k+1}(M)$ satisfying the Leibniz rule with respect to the wedge product $\wedge$ and the nilpotent property $dd = 0$. A key property of differential forms is that they can be integrated over any orientable $k$-submanifolds of $M$. For any oriented $(k\!+\!1)$-submanifold $S\subset M$ with boundary $\partial S$, Stokes' theorem, as a generalization of the Newton-Leibniz rule, states that the integral of a differential $k$-form $\omega$ over $\partial S$ is equal to the integral of its differential over $S$, i.e.,

The differential $d$ generalizes and unifies the classical operators in vector calculus, such as gradient $\nabla$, curl $\nabla\times$, and divergence $\nabla\cdot$ in $\mathbb{R}^2$ and $\mathbb{R}^3$. For instance, in $\mathbb{R}^3$, $0$-forms and $3$-forms can be identified with scalar fields, while $1$-forms and $2$-forms can be identified with vector fields. In this case, the differential $d$ corresponds to the gradient operator $\nabla$ when applied to $0$-forms, the curl operator $\nabla\times$ when applied to $1$-forms, or the divergence operator $\nabla\cdot$ when applied $2$-forms. The nilpotent property $dd = 0$ directly leads to the vector field analysis identities $\nabla\times\nabla = 0$ and $\nabla\cdot\nabla\times = 0$.

A differential form $\omega\in\Omega^k(M)$ is called closed if $d\omega = 0$, or exact if there is a $(k\!-\!1)$-form $\zeta\in\Omega^{k-1}(M)$ such that $\omega = d\zeta$. Due to the property $dd = 0$, every exact form is closed. Thus, the differential $d$ links the sequence of the spaces of differential forms on $M$ into a co-chain complex

The $k$-th de Rham cohomology group, denoted by $H^k_{dR}(M)$, is then defined to be the $k$-th homology of this chain complex, i.e., the quotient space of closed $k$-forms modulo the space of exact $k$-forms, i.e.,

The de Rham cohomology, by the de Rham theorem, is naturally isomorphic to the singular cohomology, and thus depends only on the manifold topology.

Let $g$ be a Riemannian metric on $M$ and $\langle\cdot, \cdot\rangle_g$ be the point-wise inner product induced by $g$ on $\Omega^k(M)$. The Hodge star operator $\star$ provides an isomorphism from the space of differential $k$-forms $\Omega^k(M)$ to the space of $(m\!-\!k)$-forms $\Omega^{m-k}(M)$, defined by the following formula:

where $\mu_g$ is the volume form on $M$ induced by $g$. The Hodge $L^2$-inner product on the space of $k$-forms $\Omega^k(M)$ can then be obtained by taking the integral of the formula (2.4)

The codifferential $\delta: \Omega^k(M)\to\Omega^{k-1}(M)$ is defined by

which also has the nilpotent property $\delta\delta = 0$. We call a differential form $\omega\in\Omega^k(M)$ co-closed if $\delta\omega = 0$, or co-exact if there is a $(k+1)$-form $\eta\in\Omega^{k+1}(M)$ such that $\omega = \delta\eta$. The codifferential $\delta$, as the differential $d$, also extends the classical gradient, curl, and divergence in vector calculus. In $\mathbb{R}^3$, it corresponds to $-\nabla\cdot$, $\nabla\times$, and $-\nabla$ when applied to $1$-forms, $2$-forms, and $3$-forms, respectively.

The Hodge Laplacian for differential forms is defined as $\Delta = d\delta + \delta d: \Omega^k(M)\to\Omega^k(M)$. Its kernel, consisting of all differential $k$-forms $\omega$ on $M$ with $\Delta\omega = 0$, is called the space of harmonic $k$-forms. We denote by $\mathcal{H}^k_{\Delta}(M)$ the space of harmonic $k$-forms and by $\mathcal{H}^k(M)$ the space of $k$-forms that are both closed and co-closed, i.e., $\mathcal{H}^k(M) = \ker d\cap \ker\delta$. The latter space $\mathcal{H}^k(M)$, known as the space of harmonic $k$-fields, is in general only a subset of the space of harmonic forms $\mathcal{H}^k(M)\subset\mathcal{H}^k_{\Delta}(M)$, and is infinite-dimensional ^[53]. However, in the case of closed manifolds where $\partial M = \emptyset$, the space of harmonic forms $\mathcal{H}^k_{\Delta}(M)$ reduces to the space $\mathcal{H}^k(M),$ as any harmonic form is both closed and co-closed. The result follows directly from the following formula:

due to the $L^2$-adjointness of the codifferential $\delta$ and the differential $d$ on closed manifolds, i.e., $(d\omega, \eta) = (\omega, \delta\eta)$.

The classical Hodge decomposition theorem for closed manifolds states that the space of differential $k$-forms $\Omega^k(M)$ can be decomposed as

These three subspaces are mutually orthogonal with respect to the inner product (2.5). Moreover, Hodge theorem identifies the harmonic space $\mathcal{H}^k_{\Delta}(M)$ with the $k$-th de Rham cohomology group $H^k_{dR}(M)$, which states that each harmonic form corresponds to exactly one equivalence class in $H^k_{dR}(M)$. Therefore, the harmonic space $\mathcal{H}^k_{\Delta}(M)$ is fully determined by the manifold topology, and is finite-dimensional with its dimension given by the Betti number $\dim\mathcal{H}^k_{\Delta}(M) = \beta_k$.

2.1.   Hodge decomposition for manifolds with boundary

In the presence of a nonempty boundary $\partial M$, the two operators $d$ and $\delta$ are not $L^2$-adjoint, as integration by parts leads to ^[69]

which contains a boundary term that may not vanish, and thus the decomposed subspaces in (2.8) are not orthogonal. However, certain boundary conditions can be enforced, ensuring the adjointness of the differential $d$ and the codifferential $\delta$, thereby inducing an orthogonal decomposition of the space of differential forms.

The most common choices of boundary conditions ensuring the adjointness of $d$ and $\delta$ are the normal (Dirichlet) and tangential (Neumann) boundary conditions. A differential form $\omega\in\Omega^k(M)$ is called normal (Dirichlet) if it gives zero when applied to tangent vectors of the boundary, or tangential (Neumann) if the same holds for its dual $\star\omega$ instead. Denote by $\Omega^k_n(M)$ the set of normal differential $k$-forms and by $\Omega^k_t(M)$ the set of tangential differential forms, i.e.,

Following their definitions, the spaces $\Omega^k_n(M)$ and $\Omega^{m-k}_t(M)$ are isomorphic under the Hodge star operator $\star$, also known as the Hodge duality. Moreover, the differential $d$ preserves the normal boundary conditions, while the codifferential $\delta$ preserves the tangential boundary conditions.

The Hodge-Morrey decomposition ^[43] states that there is a 3-component $L^2$-orthogonal decomposition

The orthogonality of the decomposition directly comes from the adjointness of $\delta$ and $d$ when enforcing the normal or tangential boundary conditions. For $\omega\in\Omega^k(M),$ there is a unique decomposition of $\omega$ given as follows:

where $\alpha_n\in \Omega^{k+1}_n(M)$, $\beta_t\in\Omega_t^{k+1}(M)$, and $\eta\in\mathcal{H}^k(M)$. Note that the potentials $\alpha_n$ and $\beta_t$ are not uniquely determined as all $\alpha_n+d\eta$ and $\beta_t+\delta\gamma$ with any $\eta \in\Omega_n^{k-2}(M)$ and $\gamma\in\Omega_t^{k+2}(M)$ serve as potentials for the same components. However, the issue can be addressed by enforcing gauge conditions, such as

The potentials $\alpha_n$ and $\beta_t$ can then be uniquely determined by the following equations:

by resolving the (finite) rank deficiencies of $\Delta$ under these boundary conditions (Eqs (2.10), (2.11), (2.14), and (2.15)).

Remark 1. In the case that $M$ is a closed manifold, i.e., $\partial M = \emptyset$, both the spaces $\Omega^k_n(M)$ and $\Omega^k_t(M)$ coincide with the space of differential forms $\Omega^k(M)$, and the space of harmonic fields is identical to the space of harmonic forms $\mathcal{H}^k(M) = \mathcal{H}^k_{\Delta}(M)$. The Hodge decomposition (2.12) then reduces to the classical Hodge decomposition (2.8) for closed manifolds.

Remark 2. The Hodge-Morrey decomposition (2.12) in the low-dimensional Euclidean spaces $\mathbb{R}^2$ and $\mathbb{R}^3$, often referred to as the Helmholtz-Hodge decomposition in vector calculus, states that any vector field $\mathbf{v}$ defined on a compact domain can be orthogonality decomposed as

where $f$ is a scalar potential that vanishes on the boundary of the domain, $\mathbf{u}$ is a vector field orthogonal to the boundary, and $\mathbf{h}$ is the harmonic vector field satisfying $\nabla\times \mathbf{h} = 0$ and $\nabla\cdot \mathbf{h} = 0$. The first component $\nabla f$ and the second component $\nabla\times\mathbf{u}$ are often called the curl-free and divergence-free parts of the vector field $\mathbf{v}$ respectively. Note that in the presence of a boundary, the resulting scalar potential $f$ is also called satisfying the normal boundary of $0$-forms, and the vector field $\mathbf{u}$ is called satisfying the tangential boundary condition of $2$-forms, which are direct counterparts of the potentials $\alpha_n$ and $\beta_t$ in (2.12). For a complete correspondence between scalar or vector fields, and differential forms under the normal and tangential boundary conditions, see ^[69].

The space of harmonic fields $\mathcal{H}^k$, in general, is infinite-dimensional, and thus has no direct correspondence with the cohomology of the manifold. However, as noted early ^[69], one can restrict to the space of normal harmonic fields, namely, $\mathcal{H}^k_n(M) = \mathcal{H}^k(M)\cap\Omega^k_n(M)$, and the space of tangential harmonic fields, $\mathcal{H}^k_t(M) = \mathcal{H}^k(M)\cap\Omega^k_t(M)$. As a consequence of the de Rham map, these two subspaces $\mathcal{H}^k_n(M)$ and $\mathcal{H}^k_t(M)$ are fully determined by the topology of $M$: the space of normal harmonic fields $\mathcal{H}^k_n(M)$ is isomorphic to the relative de Rham cohomology $H^k_{dR}(M, \partial M)$, while the space of tangential harmonic fields $\mathcal{H}^k_t(M)$ is isomorphic to the absolute de Rham cohomology $H^k_{dR}(M)$ ^[25]. The two subspaces $\mathcal{H}^k_n(M)$ and $\mathcal{H}^k_t(M)$ are thus finite-dimensional, with dimensions given by the Betti numbers: $\dim\mathcal{H}^{k}_n(M) = \beta_{m-k}$ and $\dim\mathcal{H}^{k}_t(M) = \beta_k$. Furthermore, the kernels of the Hodge Laplacian $\Delta$, when restricted to the space of normal forms $\Omega^k_n(M)$ and the space of tangential forms $\Omega^k_t(M)$ with gauge conditions on the boundary, can be identified to the space of normal harmonic fields and the space of tangential harmonic fields, respectively. Denote by $\Delta_{n}$ and $\Delta_{t}$ the restrictions of the Hodge Laplacian $\Delta$ on the space of normal fields $\Omega^k_n(M)$ satisfying Eq (2.14) and the space of tangential fields $\Omega^k_t(M)$ satisfying Eq (2.15), i.e., $\Delta_n:\Omega^k_n(M)\to\Omega^k(M)$ and $\Delta_t:\Omega^k_t(M)\to\Omega^k(M)$. Then, immediately we have $\ker\Delta_n = \mathcal{H}^k(M)\cap\Omega^k_n(M) = \mathcal{H}^k_n(M)$ and $\ker\Delta_t = \mathcal{H}^k(M)\cap\Omega^k_t(M) = \mathcal{H}^k_t(M)$. The result follows directly from Eq (2.7). These identifications finally enable us to study the topology of the underlying manifold $M$ through the Hodge Laplacians on normal and tangential forms.

Remark 3. In fact, let $\mathcal{H}^k_{co} = \mathcal{H}^k(M)\cap \delta\Omega^{k+1}(M)$ and $\mathcal{H}^k_{ex} = \mathcal{H}^k(M)\cap d\Omega^{k-1}(M)$. The space of harmonic fields $\mathcal{H}^k(M)$ can be further orthogonally decomposed for smooth manifolds

which results in the Hodge-Morrey-Friedrichs decomposition given as follows:

In particular, if $M$ is a compact domain in Euclidean spaces, then there is a unique orthogonal $5$-component decomposition

as the spaces $\mathcal{H}^k_n(M)$ and $\mathcal{H}^k_t(M)$ are $L^2$-orthogonal, instead of just being transversal for compact manifolds in general ^[56]. Due to the correspondence between differential forms and vector fields in the low-dimensional Euclidean spaces, the implementation of this $5$-component Hodge decomposition has been applied and implemented to the study of vector fields for surface triangle meshes, for tetrahedral meshes ^[69] and for regular Cartesian grids ^[58].

As we mainly focus on applications of compact domains in $\mathbb{R}^3$, to study the geometric and topological information of the underlying manifolds, there are eight Laplacians to be considered, which are defined on the spaces of differential $k$-forms with $k = 0, 1, 2, 3$ satisfying either the normal or the tangential boundary conditions. However, thanks to the duality between the space of normal fields and tangential fields, the study of the spectra of these eight Laplacians reduces to that of four Laplacians on one of the two types of boundary conditions, and finally to the singular spectra of three differential operators, applied to differential forms of degree $k = 0, 1, 2, 3$ ^[18]. Further details will be discussed in the next section for the discretization of Laplacians.

3.   Discretization and construction of Laplacians

In this section, we elaborate on the discretization of the Hodge Laplacian and introduce the boundary-induced graph (BIG) Laplacian for compact domains in low-dimensional Euclidean spaces ^[52]. Although the theory works for 2D compact domains, for the remainder of the paper we focus only on compact domains in $\mathbb{R}^3$, as we target mainly 3D applications. We use discrete exterior calculus (DEC) to discretize all differential operators and differential forms on regular Cartesian grids, as it allows for efficient and accurate numerical algorithms relying on just matrix algebra, while keeping the $L^2$ orthogonality between different components in Hodge decomposition. In addition, the constructed discrete differential operators and differential forms in DEC approximate their smooth analogs. For the characterization of the underlying manifold, we choose the Eulerian formulation, where the manifold is given as a sublevel set of a level set function defined on a regular Cartesian grid. Another common way, called the Lagrangian formulation, discretizes the manifold as simplicial meshes, i.e., triangular or tetrahedral meshes in 2D or 3D. The spectrum analysis of the Hodge Laplacians has been discussed in ^[69] for the Lagrangian formulation and in ^[58] for the Eulerian formulation. Compared to the Lagrangian case, the Eulerian representation uses vertices, edges, faces, and cells all fixed in a Cartesian grid, which significantly simplifies the data structures and algorithms. The Hodge stars, in the latter case, are close to rescaled identity matrices. This fact simplifies the study of Hodge Laplacians to that of BIG Laplacians with no Hodge stars involved, and thus leads to algorithms with efficient computations.

3.1.   Discretization on entire grid

Denote by $I_m$ a rectangular $m$-dimensional regular Cartesian grid with $k$-cells oriented according to their alignments with the coordinate axes. The entire grid $I_m$ can be treated as a cell complex tessellating a rectangular domain in $\mathbb{R}^m$, where each $k$-cell is a $k$-dimensional hypercube with edge length $\ell$. A continuous differential $k$-form $\omega$ on $I_m$, following the de Rham map, can be discretized by its integral value over each oriented $k$-cell $\sigma_i$, given as $W^i = \int_{\sigma_i} \omega$ ^[19]. The discrete differential on discrete $k$-forms of the grid $I_m$ is then encoded by a sparse matrix $D^I_k$, which stores the signed incidence between $(k\!+\!1)$-cells and $k$-cells and is given as the transpose of the cell boundary operator $\partial_{k+1}^T$ on $(k\!+\!1)$-cells following from Stokes' theorem $\int_{\sigma} d\omega = \int_{\partial \sigma} \omega$. An illustration of the chain complex formed by boundary operator $\partial$ for a simple grid complex with a single 2D cell can be seen in Figure 1, which is a straightforward generalization of the chain complex on simplicial complexes. Note that the boundary of the boundary of a cell always results in a 0 chain, i.e., $\partial\partial = 0$, whose transpose immediately produces $D^I_{k+1}D^I_k = 0$, thus preserving the nilpotent property in the continuous setting.

The discrete Hodge star establishes a one-to-one correspondence between discrete $k$-forms on the primal grid $I_m$ and discrete $(m\!-\!k)$-forms on its dual grid, given as the translated grid with grid points located at the $m$-cell centers of $I_m$, based on the following formula:

where $\star\sigma_k$ is the dual $(m\!-\!k)$-cell formed by the dual grid points located at the centers of the primal $m$-cells incident to $\sigma_k$. See Figure 2 for an illustration of the correspondences between the primal and dual cells in the Cartesian grid case. Following from the discretization of differential forms, this correspondence leads to a diagonal matrix $S^I_k$ with diagonal entries given by the ratio between the volumes of the dual $(m\!-\!k)$-cells and the primal $k$-cells, $\ell^{m-k}/\ell^k = \ell^{m-2k}$. The associated discrete Hodge $L^2$-inner product (2.5) of two discrete $k$-forms $V_k$ and $W_k$ on grid $I_m$ is then given by

The discrete codifferential, by definition of its smooth counterpart (2.6), can be assembled from the discrete differential and Hodge star operators as $\delta^I_k = (S^I_{k-1})^{-1}D^I_{k-1}S^I_k$. Note that the discrete counterpart of the Hodge Laplacian $\Delta = d\delta + \delta d$ by replacing the differential and codifferential operators results in a nonsymmetric matrix. Instead, we consider the counterpart of $\star\Delta$ as the discrete Hodge Laplacian given by

where the operators are considered to be null for $k < 0$ or $k > m$.

3.2.   Discrete differential forms and operators on $M$

Compared to the case of simplicial or polygonal meshes, where the projection matrices to the interior can be straightforward to implement with the boundary elements explicitly labeled, modeling the manifold $M$ as the volume bounded by a level set surface leads to delicate computation of the projection matrices. Note that the boundary of $M$ using grid representation typically intersects with boundary $k$-cells instead of being its supersets. We restrict the computation to relevant cells by implementing the two types of boundary conditions through the inclusion or exclusion of the entire $k$-cells. We use the strategy as in ^[58] for the computation of projection matrices for each type of boundary condition: for the normal boundary condition, we include all cells if at least one of its vertices is inside or on the boundary of $M$, while for the tangential boundary condition, we include all cells with at least one of the vertices of the corresponding dual cells is inside or on the boundary. We refer to the former set of cells as the normal support and the latter as the tangential support. In contrast to the mesh case, it is important to note that neither the normal nor the tangential support is necessarily a superset of the other. See Figure 3 for one example showing the distinction of these two supports for $1$-forms.

In the computation of the discrete Hodge star operators, it is essential to consider and incorporate the boundary conditions. Following the procedure in ^[58], we keep the dual cell volumes and adjust the primal cell volumes for normal boundary conditions, and do conversely for tangential boundary conditions with the primal cell volumes kept and the dual cell volumes changed. To be specific, when dealing with normal (resp., tangential) boundary conditions, we only compute the volume of the region of the primal (resp., dual) $k$-cells within the boundary $\partial M$ for the denominator (resp., numerator) of the ratio in the discrete Hodge star matrix, and leave the dual (resp., primal) cell volumes in the numerator (resp., denominator) unchanged. Each unaltered $k$-cell has a $k$-volume of $\ell^k$. In addition, for numerical stability, we do not alter the volume of outside primal $k$-cells, and perturb the level set function evaluated at primal/dual grid points to have an absolute value above $\epsilon = 10^{-5}\ell$, which ensures well-behaved fractional $k$-volumes. We denote by $S^I_{k, n}$ and $S^I_{k, t}$ the diagonal Hodge star matrices defined on the entire grid $I^m$ corresponding to the normal and tangential boundary conditions, respectively.

The projection matrix to the corresponding support, for each type of boundary condition, can be constructed from the identity matrices by eliminating the rows corresponding to $k$-cells outside the support. Denote by $P_{k, n}$ the projection matrix for $k$-cells onto the normal support and by $P_{k, t}$ the one onto the tangential support. We then obtain a new set of differential and Hodge star operators for $M$:

The nilpotent property $D_{k+1, n}D_{k, n} = 0$ and $D_{k+1, t}D_{k, t} = 0$ still holds for both boundary conditions due to $D^I_{k+1}D^I_k = 0$ and the following observations:

The discrete Hodge $L^2$-inner products of the two types of discrete $k$-forms on the manifold $M$ for these two boundary conditions are then given by

whose domains are the discrete $\Omega^{k}_n(M)$ and the discrete $\Omega^{k}_t(M)$, respectively. Finally, we assemble the two types of discrete Hodge Laplacians as in the mesh case:

The null spaces of these discrete Hodge Laplacians, as in the continuous case, are fully determined by the topology of the underlying manifold $M$, since they only depend on the differential and projection matrices. The dimension of the kernel of $L_{k, n}$ is given by the Betti number $\beta_{m-k}$, while the dimension of the kernel of $L_{k, t}$ is given by $\beta_{k}$. Here, the Betti number $\beta_k$ presents directly the number of $k$-dimensional holes on the manifold $M$. For instance, $\beta_0$ gives the number of connected components, $\beta_1$ gives the number of tunnels, and $\beta_2$ provides the number of closed cavities, respectively. The spectra of these Laplacians, in addition, could be used to study the geometric information of the manifold. It is known that the nonzero eigenvalues of the Laplacians provide rich insights into the shape of a manifold. For instance, the Fiedler value, defined as the smallest nonzero eigenvalue of a graph Laplacian, describes connectivity. As another example, the multiplicity of eigenvalues can reveal certain symmetries of the shape.

Remark 4. The two types of discrete Hodge Laplacians (3.9) not only provide rich geometrical and topological information of the underlying manifold, but also play a central role in the computation of the discrete Hodge decomposition (2.22) of differential forms for compact domains in 2D and 3D Euclidean spaces. In particular, they can be utilized, by resolving the rank deficiencies, to compute the potentials of the decomposed components in Hodge decomposition on normal or tangential support satisfying the corresponding boundary conditions. In addition, as the kernel sizes of Laplacians are finite, their eigenvectors corresponding to $0$ eigenvalues, for each $k$, form a basis for the space of normal or tangential harmonic fields.

Note that the discrete Hodge stars in the Eulerian setting are almost identical to rescaled identity matrices. Therefore, the computations of the Hodge Laplacian can be further simplified by replacing the Hodge stars with identity matrices, leading to the definition of the BIG Laplacians as follows:

The BIG Laplacians were introduced in ^[52] for bounded domains to facilitate the comparison and contrast of the Hodge Laplacians and the combinatorial Laplacians. They preserve the Hodge Laplacian's capability to perform differential calculus but also retain the discrete nature of combinatorial Laplacians. The convergence of the spectra of the BIG Laplacians to Hodge Laplacians has been discussed in ^[52], showing that the spectra of (3.11) converge to those of Hodge Laplacians up to a scaling value $\ell^{-2}$ when enforcing the boundary conditions. This scaling value $\ell^{-2}$ is exactly the ratio between the missing scaling factor $\ell^{m-2(k+1)}$ in $L_k$ and the missing factor $\ell^{m-2k}$ of $S_k$. As the BIG Laplacians produce results similar to those obtained from the discrete Hodge Laplacians with less computation, they can also be used to study the geometric and topological information of the underlying manifolds.

Note that the dual grid is also a Cartesian grid staggered with the primal grid by a replacement of $\ell/2$ in all three axial directions of the Cartesian coordinates. For the study of the spectra of these Laplacians, one only needs to implement one type of boundary condition, for instance, the normal boundary condition, as $L_{k, n}$ defined on the primal grid with normal boundary conditions is equivalent to $L_{m-k, t}$ defined on its dual grid with tangential boundary conditions.

3.3.   Topology-preserving construction of Laplacians

Preserving the topological structure is a major characteristic of the present work. However, preserving topological structure in various Laplacian operators is a nontrivial job. In this section, we present the detailed construction of topology-preserving Laplacians.

Note that, on the grid, the Hodge Laplacians and the BIG Laplacians are of the same sparsity patterns. For simplicity in exposition when discussing the spectrum analysis of the Laplacians, we let $L_k$ be a generic Laplacian matrix of the form

Here, the Laplacian $L_k$ can be interpreted, under choices of boundary conditions and Hodge star accuracy, as either a Hodge Laplacian or BIG Laplacian (with $S_k$ set to identity) under tangential or normal boundary condition. The eigenvalues and eigenvectors of $L_k$ can be solved by considering the generalized eigenvalue problem

where $\lambda$ is an eigenvalue and $W$ is the associated eigenvector. To analyze the results, we perform the following transformation in the space of discrete forms: $\bar D_k = S_{k+1}^{1/2}D_kS_k^{-1/2}$, $\bar L_k = S_k^{-1/2}L_k S_k^{-1/2}$, and $\bar W = S_k^{1/2}W$. Rewriting the formulas above yields a simplified form of the Laplacian

and a regular eigenvalue problem:

Note that the property $\bar{D}_{k}\bar{D}_{k-1} = 0$ is preserved. As the nonzero eigenvalues of $\bar{D}_k^T\bar{D}_k$ and $\bar{D}_{k}\bar{D}_{k}^T$ for each $k$ are the same, given by the squared nonzero singular values of the discrete differential $\bar{D}_k$, and each Laplacian $\bar L_k$ is just the combination of $\bar{D}_k^T\bar{D}_k$ and $\bar{D}_{k-1}D_{k-1}^T$, the entire spectrum of the Laplacians can thus be studied through the singular values of discrete differentials. Let

be the singular value decomposition of $\bar D_k$, where $U_{k+1}$ and $V_k$ are orthogonal matrices and $\Sigma_k$ is a rectangular diagonal matrix with diagonal values given by the singular values of $\bar D_k$. It follows immediately from $\bar{D}_k\bar{D}_{k-1} = 0$ that

Therefore, the columns of $V_k$ corresponding to nonzero singular values of $\bar{D}_k$ are orthogonal to columns of $U_k$ associated with nonzero singular values of $\bar{D}_{k-1}$. In addition, it follows from

that the spectrum of $\bar{L}_k$ is given by the union of squared nonzero singular values of $\bar{D}_k$, $\bar{D}_{k-1}$, and $0$, with the multiplicity of $0$ given by the $k$-th Betti numbers. The columns of $U_k$ and $V_k$ corresponding to nonzero singular values, together with the set of harmonic forms, span the entire space of differential $k$-forms.

In the case that $\dim(M) = 3$, for each type of boundary condition, we have four Laplacians of different degrees in total $k = 0, 1, 2, 3$:

Due to the aforementioned discussion on the spectrum of Laplacians and the duality of the normal and tangential boundary conditions, the spectral analysis of all Laplacians can be reduced to the singular spectra analysis of the three discrete differentials $\bar D_0$, $\bar D_1$, and $\bar D_2$ with one type of boundary condition. Note that the numerical evaluation of the singular values of these differentials, in the simplicial mesh case, may differ for the two types of boundary conditions, as the degrees of freedom (DoF) for normal $k$-forms and tangent $m-k$ forms are different. However, in the Cartesian representation, they are strictly equivalent to each other by shifting the grid in all directions of the axis by $\ell/2$, so long as $M$ is at least one grid spacing away from the boundary of the grid.

For the computation of the spectra of the Laplacians, we choose the normal boundary condition. The spectra of all Laplacians $\bar{L}_{k, n}$ for compact domains in $\mathbb{R}^3$ can be finally decomposed into three distinct parts: the squared singular values of the gradient of tangential scalar fields, denoted by $T$, the squared singular values of the gradient of normal scalar fields, denoted by $N$, and the squared singular values of the curl of tangential curl fields, denoted by $C$.

4.   Persistent de Rham-Hodge Laplacians

In this section, we present the construction of the persistent de Rham-Hodge Laplacian on differentiable manifolds, which is based on the filtration of manifolds induced by varying a single parameter (the filtration parameter). The spectra of Laplacians carry rich topological and geometric information of a manifold. Essentially, a single manifold does not provide enough information in practical applications like feature extraction for machine learning analysis. As such, instead of studying just a single manifold, one could examine the spectra of a family of manifolds by adjusting the filtration parameter. The spectra of the Laplacians from this family of manifolds could provide much more information than by considering just one, as the topology and geometry could change for different parameters. This single-parameter family of manifolds, called the evolution of manifolds, was first introduced in ^[18] based on tetrahedral meshes. We briefly recap the background.

The formal definition of the evolving manifold is given by a one-parameter family of immersions $F_c = F(\cdot, c)$ with $F: B\times [a, b]\to N$ being a smooth map, where $B$ is called the base manifold, $N$ is the ambient manifold, and $c\in[a, b]$ is a real parameter within the interval. In practice, the most common way to define the evolution of manifolds without specifying $B$ is through a level set function by adjusting the isovalues. Given a function $f: N\to [a, b]$, then in our case, we consider the sublevel sets $M = \{x\in N\, |\, f(x)\leq c\}$ with the boundary given as $\partial M = \{x\in N\, |\, f(x) = c\}$ for $c\in [a, b]$. A sequence of manifolds can then be obtained by considering evenly distributed isovalues of the function $f$ with the inclusion map

where each $M_l$ is given as the sublevel set corresponding to $c_l$ with $a\leq c_0 < c_1 < \cdots < c_s \leq b$. To ensure that $M_l$ is a manifold, we assume that the function $f$ is a Morse function on $N$, and none of the $c_l$'s corresponds to a critical value of the function $f$, i.e., $f^{-1}(c_l)$ does not contain any critical points. This is always possible as the set of Morse functions on a compact manifold is dense in the space of smooth functions, and their critical points are isolated, nondegenerate, and finite for compact manifolds. Thus we can always perturb any input function slightly to avoid critical isovalues in $\{c_l, l = 0, 1, \cdots s\}$. In addition, we assume that for each $l$, $M_{l, l+1} = \overline{M_{l+1}\backslash M_{l}} = \{x\in N|\, f(x)\in[c_l, c_{l+1}]\}$ contains at most one critical point, which can be realized by refining the parameter sequence. Note that both $M_l$ and $M_{l, l+1}$ are compact. By Morse theory, if $M_{l, l+1}$ contains no critical points, $M_l$ is diffeomorphic to $M_{l+1}$. The retraction from $M_{l+1}$ to $M_l$ can be easily constructed by considering a flow along the gradient of the function. As $M_{l+1}$ is homotopic to $M_l$ in this case, there is no topological change happening between $(c_l, c_{l+1})$. For the other case when there is exactly one critical point in $M_{l, l+1}$, the manifold $M_{l+1}$ is homotopic to $M_l$ with a $k$-cell attached, where $k$ is the index of the critical point, defined to be the dimension of the largest subspace on which the Hessian $\operatorname{Hess}(f)(x)$ is negative definite. The topological change of the sublevel sets occurs precisely at the critical values of the level set function. Depending on the type of the critical points, i.e., local minimum, saddle points, and local maximum, the topology changes in different ways. In general, a local maximum has the full index $m$, a local minimum has index $0$, while saddle points have indices strictly between 0 and $m$. In the case of $\mathbb{R}^3$, the occurrences of minima and maxima correspond to the birth of the $0$-th generators and the death of the $2$nd homology generators, respectively, while the occurrences of $1$-saddle points correspond to the birth of $1$st homology generators or the death of the $0$-th homology generators, and those of $2$-saddle points correspond to the birth of $2$nd homology generators or the death of $1$st homology generators.

4.1.   Persistent harmonic forms

As the de Rham complex depends on the topology, it can also be extended to the filtration of manifolds. Due to the duality of the normal and tangential boundary conditions, without loss of generality, one may focus on the space of normal differential forms. Given $M_l\hookrightarrow M_{l+1}$, we then need to construct a map from the space of normal $k$-forms $\Omega^k_n(M_l)$ to the space of normal $k$-forms $\Omega^k_n(M_{l+1})$, which extends each normal $k$-form on $M_l$ to a normal $k$-form on $M_{l+1}$. Let $\omega\in\Omega^k_n(M_l)$. The idea is to utilize the boundary condition of $\omega$ on $M_l$ and extend the forms $\omega|_{\partial M_l}$ to exact normal forms on the domain $M_{l, l+1}$ with certain boundary conditions on $\partial M_{l, l+1} = \partial M_l\cup \partial M_{l+1}$. Then, the combination $\overline{\omega}$ defines a normal $k$-form on the manifold $M_{l+1}$. Note however that $\delta \overline{\omega}$ is only $0$ in $M_l,$ so the extension of $\omega \in \ker \delta$ may no longer be in $\ker \delta$ on $M_{l+1}.$

To be specific, we consider the biharmonic equation $\Delta^2\zeta = \Delta(\Delta\zeta) = 0$ on $M_{l, l+1}$ with both Dirichlet and Neumann boundary conditions to ensure the smoothness of $d \zeta$ with $\omega$ through $\partial M_{l}$. Note that $d \zeta$ satisfies the normal boundary condition on $M_{l, l+1}$. Let $\overline{\omega}$ be the extension of $\omega$ on $M_{l+1}$ with $\overline{\omega} = \omega$ on $M_{l}$ and $\overline{\omega} = d \zeta$ on $M_{l, l+1}$. It follows that $\overline{\omega}\in\Omega^k_n(M_{l+1})$ as it satisfies the normal boundary condition $\overline{\omega}|_{\partial M_{l+1}} = \zeta|_{\partial M_{l+1}} = 0$.

While the biharmonic equation produces a smooth extension, in practice, it is more efficient to consider the harmonic extension with the boundary condition $\Delta\zeta = 0$ with the boundary condition of $\star d\zeta = \star \omega$ on $\partial M_l$ and the typical normal form boundary condition on $\partial M_{l+1}$. The solution, by ^[53], Theorem 3.4.10], is unique. The resulting $\bar\omega$ is continuous but non-smooth as $\delta \bar\omega$ may lead to a Dirac distribution on $\partial M_l$ when $M_{l, l+1}$ induces a topological change. For instance, for a harmonic normal 1-form $\omega$ on a spherical shell $M_l$ with $M_{l+1}$ turning into a solid ball, the biharmonic extension would create a uniform divergence $\delta\bar\omega$ in $M_{l, l+1},$ whereas the harmonic extension creates a thin layer of nonzero divergence $\delta\bar\omega$ near the part of $\partial M_l$ around the cavity in the middle. Thus, the harmonic extension serves the same purpose in reducing the kernel of $\delta$.

Denote by $\mathcal{I}_{l, 1}$ the map from $\Omega^k_n(M_l)$ to $\Omega^k_n(M_{l+1})$ sending $\omega$ to $\overline{\omega}$. Note that $(d\circ\mathcal{I}_{l, 1})(\omega)$ is $0$ on $M_{l, l+1}$ and thus the same as the extension of the differential of a normal form $d\omega$ on $M_{l, l+1}$, i.e., $d\circ\mathcal{I}_{l, 1} = \mathcal{I}_{l, 1}\circ d$. It follows that there is a commutative diagram

where the horizontal direction gives the de Rham complex and the vertical direction shows the filtration-induced extensions.

Next, we introduce the $p$-persistent Hodge Laplacian. Let $\mathcal{I}_{l, p} = \mathcal{I}_{l+p-1, 1}\circ...\circ \mathcal{I}_{l, 1}$, which then gives an extension map from the space of normal forms on $M_l$ to the space of normal forms on $M_{l+p}$. We have the following commutative diagram:

Here $d_l, \delta_l$ denotes the differential and codifferential on $\Omega^k(M_l)$, $d_{l+p}, \delta_{l+p}$ denotes the differential and codifferential on $\Omega^k(M_{l+p})$, respectively, and $\mathcal{R}_{l, p}$ is the projection of differential forms in $\Omega^k_n(M_{l+p})$ to the space spanned by the harmonic extensions followed by the restriction to $M_{l}$. Let $\tilde{\delta}_{l, p} = \delta_{l+p}\circ\mathcal{I}_{l, p}$ and $\tilde{d}_{l, p} = \mathcal{R}_{l, p}\circ d_{l+p}$. By the construction of the extension, we have $(\tilde{\delta}_{l, p}\omega, \eta) = (\omega, \tilde{d}_{l, p}\eta)$, i.e., $\tilde{\delta}_{l, p}$ are adjoint to $\tilde{d}_{l, p}$. We then define the $p$-persistent Hodge Laplacian operator $\Delta^p_{n, l}: \Omega^k_n(M_l)\to\Omega^k_n(M_{l})$ as follows:

It is easy to see that when $p = 0$, the $p$-persistent Hodge Laplacian gives exactly the usual Hodge Laplacian $\Delta_{n, l}: \Omega^k_n(M_l)\to\Omega^k_n(M_{l})$ restricted to the space of normal forms. We then define the $p$-persistent normal harmonic fields as the kernel of the $p$-persistent Hodge Laplacian $\mathcal{H}_n^{k, p} = \ker\Delta^p_{n, l}$, which can be identified with the space $\ker\tilde{\delta}_{l, p}\cap\ker d_{l}$. Note that by the extension construction and $\mathcal{R}_{l, p}\circ\mathcal{I}_{l, p} = \operatorname{Id}$, one can see that $\ker\tilde{\delta}_{l, p}\subset\ker\delta$ gets smaller as $p$ increases, which confirms that fewer cohomology generators persist longer.

4.2.   Discretization of $p$-persistent de Rham cohomology

The regular Cartesian grid allows one to define persistent graph Laplacian on manifolds in the same way as persistent graph Laplacian ^[61]. It also allows defining persistent Hodge Laplacian in a consistent way, with the inclusion of nontrivial Hodge stars.

Recall that the discrete differential $k$-forms can be seen as a $k$-co-chain, i.e., a linear mapping from the chain space $\mathcal{C}_k$ to $\mathbb{R}$ that sends a $k$-chain $c_k = \sum_i a_i\sigma_i$ to $\int_{c_k}\omega = \sum_i a_iW_i$, where $W_i = \int_{\sigma_i}\omega$ is the integral of a smooth $k$-form $\omega$ over the $k$-cell $\sigma_i$.

By varying the isovalue of the level set function $f$, we can get a sequence of cell complexes given as nested sequences of sub-cell complexes of $K$ satisfying the normal boundary conditions.

See Figure 4 for an example of such a nested sequence of sub-cell complexes in a 2D Cartesian grid. Denote by $\mathcal{C}^k(K_l)$ the space of discrete $k$-forms on sub-complex $K_l$ with $0\leq l\leq s$. Note that $K_l\subset K_{l+1}$. A discrete $k$-form on $K_l$ can be easily extended to $K_{l+1}$ by solving the discrete Laplace equation with the above boundary conditions for values on every $k$-cells in $K_{l, l+1} = \mathrm{Cl} (K_{l+1}\backslash K_l)$, the closure of the difference complex. We denote this extension map as $I_{l, 1}: \mathcal{C}^k(K_l)\to \mathcal{C}^k(K_{l+1})$ and by $I_{l, p} = I_{l+p-1, 1}\circ I_{l+p-2, 1}\circ \cdots \circ I_{l, 1}: \mathcal{C}^k(K_{l})\to\mathcal{C}^k(K_{l+p})$ the extension mapping from the space of discrete $k$-forms on $K_l$ to the space of discrete $k$-forms on $K_{l+p},$ which may also be constructed directly by solving the Laplace equation on $K_{l, l+p} = \mathrm{Cl}(K_{l+p}\backslash K_l)$. With this extension mapping, the space of discrete $k$-forms on $K_{l}$ can be seen as a subspace of discrete $k$-forms on $K_{l+p}$.

A sequence of the discrete de Rham co-chain complexes can be defined as follows:

where $D_l^k: \mathcal{C}^{k+1}(K_l)\to \mathcal{C}^{k}(K_l)$ denotes the discrete differential operator, and $\delta_l^k: \mathcal{C}^k(K_l)\to \mathcal{C}^{k-1}(K_l)$ denotes the discrete codifferential operator on $K_l$.

To define the persistent discrete Hodge Laplacian, we construct the discrete counterparts of $\tilde{d}_{l, p}$ and $\tilde{\delta}_{l, p}$ in the previous section.

Denote by $\delta_{l, p}^{k+1, n}:\mathcal{C}^{k+1}(K_l)\to\mathcal{C}_{l, p}^{k}$ the operator given as $\delta_{l, p}^{k, n} = \delta_{l+p}^k I_{l, p}^{k, n}$, where $\delta_{l+p}^{k, n}$ is the previously defined discrete operator for $K_{l+p}$ and $I_{l, p}^{k, n}$ is the discrete harmonic extension operator defined next. Assuming $K_{l, l+p}$ contains few $k$-cells, the harmonic extension is then constructed by the linear system $L^{k-1, n}_{K_{l, l+p}} \zeta = 0,$ and shifting all $\star d\zeta$ values in the overlap of supports of $K_l$ and $K_{l, l+p}$ to the righthand side and replacing them with a rescaling of $\star \omega$ based on the $k$-volume within each support. More specifically, the resulting system is $\tilde{L}^{k-1, n}_{K_{l, l+p}}\tilde \zeta = -S^{k-1, n}\delta^k_{\partial K_l} \omega,$ where $\tilde{L}^{k-1, n}_{K_{l, l+p}}$ is the Laplace operator applied to a form $\tilde\zeta$ defined on $K_{l, l+p}\backslash \partial K_l,$ and $\delta^k_{\partial K_l}$ is the boundary codifferential operator that uses the values of $\omega$ on $\partial K_l$ to evaluate the neighboring $(k\!-\!1)$-cells in $K_{l, l+p}\backslash \partial K_l.$

The resulting extension operator

where $\operatorname{Id}_{K_l}$ is the identity matrix in $K_l$ up to a rescaling in the boundary, provides the combination of $\omega$ in $K_l$ and $d\tilde\zeta$ in $K_{l, l+p}\backslash \partial K_l,$ when applied to $\omega.$ The matrix corresponding to $I_{l, p}$ is dense for rows corresponding to cells in $K_{l, l+p}$ but diagonal for rows corresponding to cells in $K_l.$ Note that $\delta_{\partial K_l}$ is not necessarily $0$ for co-closed $\omega$, but is $0$ for co-exact $\omega.$

The adjoint operator of $\delta_{l, p}^{k+1, n}$ defines $D^k_{l, p}$. In the following, we drop most of the subscripts for clarity. Recall that $(\omega, \tilde{d} \eta) = (\tilde{\delta} \omega, \eta)$ can be discretized as

with $W$ and $E$ as discrete versions of $\omega$ and $\eta$. Thus, $\tilde{D} = S^{-1} I_{l, p}^T S D,$ from which we may recognize the restriction operator as $R = S^{-1}I_{l, p}^T S$. This restriction operator can be seen as the $L_2$-projection onto the space formed by all harmonic extensions from $\Omega_n^K(M_l).$

Note that in this case, immediately $\delta_{l, p}^k\delta_l^{k+1} = 0$ since the extension operator will generate $\tilde\zeta = 0$ for any co-exact form $\omega = \delta\beta$ on $K_l$ as the righthand side of the associated linear system essentially corresponds to $\delta\delta\beta = 0$. From the adjoint version, we have $D_l^{k}D_{l, p}^{k-1} = 0,$ and thus the following commutative diagram

The discrete $p$-persistent Hodge Laplacian is then given as follows:

and the discrete $p$-persistent BIG Laplacian is

We now present some examples of evolving manifolds and show results for the spectral analysis of their persistent Laplacians. In particular, we focus on the changes of Betti numbers $\beta_0$, $\beta_1$, and $\beta_2$ and the first nonzero eigenvalues $\lambda^T_1, \lambda^C_1$ and $\lambda^N_1$ of the 0-persistent BIG Laplacians in the set T, C, and N, respectively, as introduced in Section 3.3. Four models are considered, including the Bimba model, the kitten model, a genus-3 model, and a four-ball model. For each model, we show on the top row snapshots of evolving manifolds at five evenly spaced isovalues in a chosen interval, and on the bottom row the changes in Betti numbers and the first nonzero eigenvalues $\lambda^T_1, \lambda^C_1$, and $\lambda^N_1$. All the evolving manifolds are generated using isovalues of the signed distance function (SDF) from the original surface model, given as the 0-isosurface of the SDF. As we show below, these values from the evolution of manifolds provide rich information rather than considering just a single manifold. The discontinuity of these variables indicates the topological changes occurring during the evolution process, and the monotonicity of these nonzero eigenvalues reveals the geometric changes.

The results for the Bimba model are presented in Figure 5 with an isovalue interval $[0, 0.2]$. As there is no topological change happening in the evolution process, all Betti numbers $\beta_0$, $\beta_1$, and $\beta_2$ remain constant, and $\lambda^T_1, \lambda^C_1$, and $\lambda^N_1$ are continuous throughout the whole process. Both $\lambda^C_1$ and $\lambda^N_1$ decrease as the isovalue increases.

Figure 6 illustrates the results for the kitten mode with one tunnel formed by its tail. The isovalue interval $[0, 8]$ is considered. One can see all variables are continuous during the evolution process except that $\beta_1$ and $\lambda^C_1$ both drop at the same isovalue, where $\beta_1$ changes from $1$ to $0$. This happens due to the disappearance of the tunnel. In addition, $\lambda^T_1$ increases at the beginning, and then slows down its rate of increase at the isovalue after the tunnel disappears, and $\lambda^C_1$ and $\lambda^N_1$ decrease during the evolution process.

Note that there are also tunnels in the evolving manifolds for the genes-3 model, as we expected, and a similar phenomenon can also be observed in Figure 7 for the change of the Betti numbers and the first nonzero eigenvalues. The isovalue interval $[0.1, 4]$ is considered for this model. The disappearance of the three tunnels leads to a drop of $\beta_1$ from $3$ to $0$ and also a drop of $\lambda^C_1$. $\lambda^T_1$ initially increases, and then changes its behavior to decrease after the tunnels vanish. The evolution process results in a decrease in $\lambda^C_1$ and $\lambda^N_1$, just as in the previous two models.

The evolving process of the four-ball model with isovalue interval $[2, 3.84]$ (see Figure 8) leads to discontinuities in all Betti numbers and the first nonzero eigenvalues. As the four separate components merge in the evolution, $\beta_0$ changes from $4$ to $1$, along with a drop in $\lambda^T_1$ at the same isovalue. In addition, $\beta_1$ increases from $0$ to $3$ due to the appearance of three tunnels when the merge happens and then decreases to $0$ after the disappearance of all tunnels. The nonzero eigenvalue $\lambda^C_1$ has a drop that occurs when the tunnel vanishes, however, it is continuous when the tunnels are formed. This suggests that the continuity of $\lambda^C_1$ is only related to the death but not the birth of tunnels. One can also observe a slowdown in the rate of change of $\lambda^T_1$ following the disappearance of all tunnels. As the isolate increases further, a cavity occurs in the manifold, resulting in an increase of $\beta_2$ from $0$ to $1$ and finally a decrease from $1$ to $0$ after the cavity disappears. This topological change can also be observed in $\lambda^N_1$, where $\lambda^N_1$ becomes non-differentiable.

As illustrated by these models, changes in Betti numbers $\beta_0$, $\beta_1$, and $\beta_2$ and the first nonzero eigenvalues $\lambda^T_1, \lambda^C_1$, and $\lambda^N_1$ not only reflect the changes in topology, but also characterizes the changes in geometry for the evolution of manifolds. The rich information revealed by these variables leads to potential applications in various topological data analysis tasks.

5.   Proof-of-principle experimentation

In this section, we carry out a proof-of-principle experimental demonstration of the proposed persistent de Rham-Hodge theory-based MTL. In this approach, the problem is defined on manifolds with boundaries. Appropriate boundary conditions are implemented to match actual topological dimensions. The resulting persistent Hodge Laplacians are solved to deliver the corresponding series of eigenvectors and eigenvalues at various scales. In this approach, we use these eigenvalues for machine learning predictions of protein-ligand binding affinity. The binding affinity describes the strength of protein-ligand interactions for each protein-ligand complex and has been a popular subject in machine learning studies ^[6,66]. The feature extraction algorithm, i.e., the algorithm for computing eigenvalues of Laplacians, is implemented in MATLAB, while the machine learning model is employed using a gradient boosting regressor (GBR) module from Scikit-learn.

We consider two benchmark datasets, PDBbind-v2007 and PDBbind-v2016 ^[38,51], to demonstrate the effectiveness of our framework in capturing the topological features of protein-ligand complexes. The datasets can be downloaded from http://pdbbind.org.cn/. These two PDBbind datasets provide collections of biomolecular complexes in PDB with experimentally a measured binding affinity for each protein-ligand complex, and are commonly used in various studies such as drug-discovery or molecular recognition, etc ^[10,38]. We aim to build a machine learning model, by utilizing the topological and geometric features of the protein-ligand complexes generated using our persistent Hodge Laplacian (PHL) framework as inputs, for predicting the protein-ligand binding affinities.

The biomolecular complexes in each PDBbind dataset are organized into three sets, including a general set, a refined set and a core set, with each set being a superset of the next ^{[7,24,36,41,57]}. In our experiments, for each dataset, we use the refined set, excluding the core set, to train the predictive model for the binding affinities of the protein-ligand complexes in the core set. The PDBbind-v2007 dataset contains a total of 1,300 complexes with 1,090 in the refined set and 195 in the core set, while the PDBbind-v2016 dataset has a total of 4,057 complexes with 3,767 in the refined set and 290 in the PDBbind core set.

5.1.   Element specific discrete to continuum mapping

The original datasets contain atomic names and coordinates, which are the so-called point cloud data. To generate manifold representations, we carry out the discrete to continuum mapping using the flexibility and rigidity index ^[47]. To compute the topological feature of each protein-ligand complex for the machine learning model, we use the element-specific approach ^[9]. Specifically, we consider the pairwise interactions between element types that are commonly found in proteins and ligands, including Hydrogen (H), Carbon (C), Nitrogen (N), Oxygen (O), and Sulfur (S) in proteins, and Hydrogen (H), Carbon (C), Nitrogen (N), Oxygen (O), Sulfur (S), Phosphorus (P), Fluorine (F), Chlorine (Cl), Bromine (Br), and Iodine (I) in ligands. These interactions result in a total of $50$ pairs of atom types for each protein-ligand complex ^[9]. However, due to the absence of H in most proteins, we reduce the number of atom pairs to $40$ in practice, ignoring the element H in all proteins. These 40 atom pairs, formed by atom types {C, N, O, S} in proteins, and atom types {H, C, N, O, S, P, F, Cl, Br, I} in ligands, along with their $xyz$ coordinates, are used to generate the topological features for each protein-ligand complex. In this paper, all atom-pair complexes are determined by a cutoff distance $12\mathring{\text{A}}$ from the ligand.

Let $\{ \mathbf{x}^{\alpha}_i, i = 1, \cdots, s\}$ be the location coordinates of all $s$ atoms in an atom pair, where $\alpha$ denotes the atom type of the atom either in the protein or in the ligand. For this atom pair, a level set function can then be obtained by considering the negative sum of Gaussian density functions defined at the $xyz$ coordinates of all atoms, given as

where $||\mathbf{x} - \mathbf{x^{\alpha}_i}||$ is the Euclidean distance from position $\mathbf{x}$ to the location $\mathbf{x^{\alpha}_i}$ of the $i$-th atom, $\tau$ is a scalar value, and $r^{\alpha}_i$ is the van der Waals radius of the $i$-th atom, determined by the atom type $\alpha$. Given an isovalue $c$, the sublevel set

defines a compact manifold in $\mathbb{R}^3$ with its boundary given by the isosurface $\partial M = \{\mathbf{x} \, | \, \rho(\mathbf{x}, \tau) = c \}$. A filtration of a manifold for the atom pair can then be obtained by choosing a list of evenly spaced isovalues of this level set function (5.1). Let $c_1 < c_2 < \cdots < c_s$ be such isovalues. We have their corresponding sublevel sets given as follows:

where $M_i$ is the compact manifold associated to isovalue $c_i$. In Figure 9, we present one example of the resulting filtration of manifolds at $3$ different isovalues for atom pair OH in protein-ligand complex 4tmn. Note that the function (5.1) is a special case of the flexibility rigidity index (FRI) density function ^[47], which has been shown computationally stable in converting discrete point cloud representations to continuous embeddings, and been used for generating protein boundary surfaces ^[18] and interactive manifolds ^[47]. Therefore, one can also make other reasonable choices of FRI density functions to generate the filtration of manifolds.

5.2.   Machine learning feature extraction

In the computation of the Laplacians, one can ideally choose a common Cartesian grid such that it contains all manifolds of interest for all protein-ligand complexes, which ensures that all Laplacians are computed consistently, making their spectra comparable for different complexes. However, as atoms are spread out in the space for different atom pairs, we need to use a sufficiently large grid with a fine resolution for accurate computation of Laplacians, which significantly increases the computational load. Instead, we consider, for each type of atom pairs, a fixed Cartesian grid, regardless of the types of protein-ligand complexes. This approach also ensures that the topological features are comparable for different protein-ligand complexes, as all spectra are computed in the same grid for all atom pairs of the same type. For simplicity, we choose a fixed grid spacing for all Cartesian grids across different atom pairs, and in our experiments, the grid spacing is fixed to be $0.549\mathring{\text{A}}$ for both PDBbind datasets.

We consider $9$ evenly spaced isovalues in the interval $[-0.5, -0.001]$ for all level set functions, which provide $9$ compact manifolds for each atom pair. Note that the level set function (5.1) is always less than $0$ and approaches $0$ as the norm of $\mathbf{x}$ increases. This interval is chosen as isovalues greater than $-0.001$ result in no change on the $0$-th Betti number $\beta_0$ of manifolds for most atom pairs, and isovalues smaller than $-0.5$ leads to high computational cost, as finer grids are necessary to resolve those isosurfaces. To ensure that the computation of Laplacians is accurate and no topological information is missing due to numerical errors caused by low resolution, we require that at least $8$ grid cells of the Cartesian grid are contained in each connected component of a manifold. We compute, for each manifold, the BIG Laplacian $L_{3, n}$ under the normal boundary condition, for which the number of its $0$ eigenvalues gives the 0-th Betti number $\beta_0$. We then use the $0$-th Betti number $\beta_0$ and the first $k$ nonzero eigenvalues of $L_{3, n}$, as the topological feature for the manifold. These $k+1$ features for each of the $9$ compact manifolds for each atom pair, amount to $(k+1)\!\times\!9\!\times\!40$ topological features for each protein-ligand complex. While we only used $9$ isovalues within this interval for generating the manifolds in our experiments, more isovalues can be considered, which gives a filtration of more manifolds for each atom pair, and finally leads to more topological features for each protein-ligand complex.

The spectra of the $0$-th Laplacian, which in our case corresponds to $L_{3, n}$ under the normal boundary condition, have proven effective and successful in many machine learning tasks ^{[10,38,41,61]}. While the Laplacians of other orders could also be used for generating more topological features, we utilize, in this preliminary test, only the spectra of $L_{3, n}$ as features for the protein-ligand complexes in the machine learning model due to the computation efficiency. The results, as shown in Section 5.4, indicate that these features are sufficient to validate our framework in the machine learning task for predicting the protein-ligand binding affinities.

5.3.   Machine learning algorithm

As a subset of artificial intelligence (AI), machine learning focuses on developing algorithms and models that allow computers to learn from and make decisions based on data. Instead of being explicitly programmed to perform a task, machine learning models identify patterns in data, enabling them to improve their performance over time. New algorithms have been constantly proposed ^[34,35], including the topology-enabled transformer ^[12]. The machine learning models for predicting protein-ligand binding affinities often fall into two categories depending on the type of input data: complex-based or sequence-based models. The complex-based methods are trained using features obtained from the 3D protein-ligand complexes, while the sequence-based models learn from the one-dimensional protein sequences and the ligand-simplified molecular-input line-entry system (SMILES) strings. In our experiments, besides the topological features from the 3D protein-ligand complexes, we incorporate protein-ligand features obtained from sequence-based models to build consensus models. To be specific, we make use of the recent pretrained transformer protein language model Evolutionary Scale Modeling-2 (ESM-2) ^[32], and the pretrained Transformer-CPZ model ^[13] for generating the protein and ligand features, respectively, and use their concatenation as inputs for the binding affinity prediction. Here CPZ represents the union of three data setscontains ChEMBL ^[26], PubChem ^[31], and ZINC ^[29]. The residue embeddings from the last layer of the pretrained ESM-2 model esm.pretrained.esm2_t33_650M_UR50D are used as the protein features, while the embeddings from the last layer of the pretrained Transformer-CPZ model chembl27_pubchem_zinc_512 are used as the ligand features.

With the topological features and the embedding features obtained from ESM-2 and Transformer-CPZ, we employ the GBR module from Scikit-learn 1.4.2 for predicting the protein-ligand binding affinities. We then use the consensus prediction from these models as the final results. The GBR parameters used in our experiments are: n_estimators = 10,000, max_depth = 5, min_samples_split = 5, learning_rate = 0.005, loss = squared_error, subsample = 0.5, and max_features = sqrt. Changing these parameters does not result in significant differences. To address the randomness from the machine learning algorithm, we repeat each modeling process 20 times with different random seeds, and use the average of predictive results. The Pearson correlation coefficients (PCC) are used as the evaluation metric to assess the performance of our proposed models.

5.4.   Experimental results

The number of topological features for each protein-ligand complex, as in Section 5.2, is given by $(1+k)\!\times\!9\!\times\!40$, where $k$ denotes the number of the first $k$ nonzero eigenvalues of the Laplacians. To find the optimal parameter $k$ leading to the best performance of predictive modules, we carry out the five-fold cross-validation on the training set of each PDBbind dataset with varying values of $k$ based on the average of PCC values. The results indicate that the optimal PCC values for the PDBbind-v2007 and PDBbind-v2016 training sets can be achieved when $k = 5$ and $k = 10$, respectively. For the PDBbind-v2007 training set, the PCC value is $0.709$ and the root mean squared error (RMSE) value is $2.049$, while for the PDBbind-v2016 training set, the PCC value is $0.748$ and the RMSE value of $1.812$. These choices of $k$ result in a total of 2,160 topological features for each protein-ligand complex in the PDBbind-v2007 dataset and 3,960 topological features for each protein-ligand complex in the PDBbind-v2016 dataset. These topological features, along with the concatenated protein-ligand features from ESM-2 and Transformer-CPZ, are then used as inputs of the gradient-boosting regressor for binding affinity prediction.

In Table 1, we report the average PCC values and the average RMSE of our models on the test set for each PDBbind dataset using only the topological features from PHL, the model using only the transformer features (TF), and the consensus module using both types of features. With the incorporation of topological features, one can see a significant improvement in PCC values when using the proposed consensus model for each dataset, compared to the model using only TF features. The best performance is achieved when using the consensus model, yielding a PCC value of $0.826$ with RMSE given as $1.954$ for PDBbind-v2007 and $0.849$ with RMSE $1.728$ for PDBbind-v2016. In addition, we present the Pearson correlation coefficients obtained from our model and those in the previous studies, with results from ^{[7,10,36,41,57]}. As illustrated in Figure 10, our model outperforms all the other models for the two PDBbind datasets. These results demonstrate the utility and effectiveness of our method in capturing the topological features.

6.   Conclusions

Although there has been tremendous success of topological data analysis TDA ^[44,45], particularly, topological deep learning TDL on point cloud data ^[9,49], there are few methods for the topological analysis of data on manifolds or manifold topological analysis ^[19]. To fill this gap, we presented a new method, persistent Hodge Laplacian PHL in the Eulerian representation, for manifold topological learning MTL of real-world data on manifolds. PHL differs from existing state-of-the-art TDA methods on point clouds in the sense that the proposed PHL is defined on manifolds, for which the traditional TDA methods do not work. Additionally, PHL extends our earlier evolutionary de Rham-Hodge theory constructed on the Lagrangian representation ^[18] to the Eulerian representation, which avoids numerical inconsistency over multi-scale manifolds. We offer two discrete Hodge stars that mimic the continuous operator and developed both a continuous theory for mapping of normal forms across manifolds in a filtration to enable persistent cohomology analysis and the associated topology-persevering discrete construction on Cartesian grids. A proof-of-principle test on two benchmark datasets validates our MTL model, highlighting its simplicity and promise for the predictions of data on manifolds.

The popularity of TDA is facilitated by effective software packages, such as JavaPlex ^[2], Perseus ^[42], Ripser ^[4], etc. The further development of efficient PHL software is an important task. The computational efficiency has not been studied in this work. Algorithm acceleration and parallel and graphics processing unit (GPU) architecture are to be explored. Further experimental validations of manifold topological learning are also needed.

Author contributions

Zhe Su: designed the project, wrote the code and draft, analyzed data, and revised the manuscript; Yiying Tong: conceptualized and designed the project, analyzed data, wrote the draft, revised the manuscript, and acquired funding; Guo-Wei Wei: conceptualized and supervised the project, wrote the draft, revised the manuscript, and acquired funding. All authors have read and approved the final version of the manuscript for publication.

Acknowledgments

This work was supported in part by NIH grants R01GM126189, R01AI164266, and R35GM148196, National Science Foundation grants DMS2052983, DMS-1761320, and IIS-1900473, NASA grant 80NSSC21M0023, Michigan State University Research Foundation, and Bristol-Myers Squibb 65109. The authors thank Dr. Hongsong Feng for his help on transformer embeddings.

Conflict of interest

The authors declare no competing interests.