Manifold topological deep learning for biomedical data

Domain: neuroscience
Published: 2026-04-01
Editorial outcome: major_revision

Abstract
Recently, topological deep learning (TDL), which integrates algebraic topology with deep neural networks, has achieved significant success in processing point-cloud data and has emerged as a promising paradigm in data science. However, TDL has not been extended to differentiable-manifold data, including images, due to the challenges introduced by differential topology. We address this challenge by introducing a manifold topological deep learning (MTDL) framework. To highlight the contribution of Hodge theory, we integrate it into a streamlined convolutional neural network within the MTDL framework. In this framework, original images are represented as smooth manifolds with vector fields that are decomposed into three orthogonal components based on Hodge theory. These components are then concatenated to form an input image for the convolutional neural network architecture. The performance of MTDL is evaluated using the MedMNIST v2 benchmark database, which comprises 717,287 biomedical images from eleven 2D and six 3D datasets. MTDL significantly outperforms other competing methods, extending TDL to a wide range of data on smooth manifolds.

Introduction
Topological deep learning (TDL) is an emerging field that integrates topological methods with deep learning techniques to perform learning tasks such as regression, classification, and representation learning [1]. Unlike traditional black-box deep learning models, TDL models offer greater interpretability by leveraging topological features and representations to explicitly capture the underlying geometric and structural properties of data. Since its introduction in 2017 [2], TDL has rapidly evolved, leading to a diverse set of methods and models. Besides its methodological advances, TDL has been successfully applied in various domains, including biology, chemistry, materials science, neuroscience, and social networks [3,4]. Current TDL models mainly focus on combinatorial data structures, such as point clouds and graphs. Compared to combinatorial data, differentiable manifold data contains richer geometric information and is more suitable for analysis using methods from differential topology, such as differential forms and differential operators. These methods enable the study of continuous, smooth phenomena that cannot be adequately captured through purely combinatorial approaches. Despite this advantage, there are currently no TDL models designed for differentiable manifold data.

Differentiable manifolds, such as curves and surfaces, are ubiquitous in real-world data. For example, object surfaces, and images are all natural examples of differentiable manifolds. Thus, extending TDL to differentiable manifold data is both meaningful and necessary. However, two primary challenges hinder this extension. Firstly, although the images are inherently manifold data, it is nontrivial to rigorously model them as differentiable manifolds while preserving essential differentiable and topological properties. Secondly, designing an efficient model that combines the mathematical methods from differential topology with deep learning poses theoretical and computational challenges.

With the advancements in Topological Data Analysis [5, 6], particularly the remarkable successes achieved by persistent homology [7], several studies have employed topological methods for manifold data analysis. For instance, using simplicial complexes or cubical complexes to model images [8] and using curves or knots to represent amino acid chains [9], then computing persistent homology for data analysis. While these methods are powerful in capturing the topological structures of manifold data across various scales, they exhibit limitations in capturing the smooth differentiable information within the manifold data, which can be addressed by incorporating methods from differential topology, such as vector fields, differential forms, and differential operators. Moreover, although differentiable manifolds have been utilized in manifold topological learning [10], such as in modeling protein-ligand complexes [11], these methods have not yet been extended to deep learning architectures.

Recently, a discrete topology-preserving Hodge theory for differentiable manifolds embedded in Cartesian grids has been introduced [12] and successfully applied to single-cell RNA velocity analysis [13]. This theory provides an efficient way for modeling images as differentiable manifolds since images are naturally embedded in Cartesian grids. On the other hand, the MedMNIST v2 dataset offers a standard and reliable benchmark for evaluating model performance in medical image classification. The MedMNIST v2 dataset contains twelve 2D datasets and six 3D datasets, covering major medical data modalities, the data scale ranges from 100 to 100000, and the task type includes binary, multi-class classification, ordinal regression, and multi-label classification, making it highly suitable for assessing the efficiency, robustness, and generalizability of models [14].

Here, we introduce a Manifold Topological Deep Learning (MTDL) Model using the de Rham-Hodge theory, a landmark of the 20th Century’s mathematics. MTDL integrates the discrete Hodge theory of differential topology and a CNN architecture, providing a novel framework for extending TDL to differentiable manifold data. In the MTDL model, the input image is represented as a discrete differentiable manifold and a vector field defined on this manifold. The Hodge Laplacian theory is then employed to decompose the vector field into three orthogonal components: curl-free, divergence-free, and harmonic parts. These components are concatenated to form a new image representation, which is passed to the CNN architecture for the prediction task. We evaluate MTDL on the MedMNIST v2 dataset, including 717,287 images from eleven 2D datasets and six 3D datasets. MTDL significantly outperforms other models, establishing MTDL as an efficient framework for TDL on differentiable manifold data.

Results
2.1 Overview of MTDL

The discrete Hodge theory on Cartesian grids provides an approach to decomposing an image into three distinct components, each capturing different geometric and topological features. To leverage this theory in TDL, we propose a MTDL model that integrates discrete Hodge theory and a CNN architecture for image classification. The architecture of MTDL is illustrated in Fig. 1. As shown in the figure, the original image is first represented as a discrete manifold on Cartesian grids under normal or tangential boundary conditions (Fig. 1a). This manifold representation establishes a mathematical formalization of the images, serving as the basis for further analysis using Hodge theory, such as using the harmonic spectra of Hodge Laplacian to detect the loop structures of the manifold (Fig. 1c). Subsequently, a vector field that encodes the image information is constructed on the discrete manifolds (Fig. 1b). There are several methods for constructing vector fields from images (section 2 of the Supplementary Information), each method providing a specific perspective on the image’s content and structure. The generated vector field is then decomposed into three orthogonal components by the Hodge decomposition. These components, including curl-free, divergence-free, and harmonic parts, are concatenated to form a multi-channel representation of the decomposed images (Fig. 1d). Finally, the resulting representation is fed into the CNN for image classification (Fig. 1e). This CNN is based on the Transformer encoder architecture by adding a maxpooling operation and replacing the multihead attention and feedforward layers with convolution operations.

2.2 Datasets and Evaluation Protocols

The MedMNIST v2 dataset [14] is an updated version of the original MedMNIST dataset [15]. It is an MNIST-like collection of standardized biomedical images comprising twelve 2D datasets and six 3D datasets that cover primary medical imaging modalities, such as X-ray, Optical Coherence Tomography (OCT), Ultrasound, Computed Tomography (CT), Electron Microscope, and Magnetic Resonance Angiography (MRA). These datasets support a wider range of classification tasks, including binary classification, multi-class classification, ordinal regression, and multi-label classification. Data sizes range from 100 to 100,000 samples. In total, MedMNIST v2 includes 708,069 2D images and 9,998 3D images, with standard train-validation-test splits provided for all datasets.

Among these datasets, BreastMNIST2D is derived from a dataset of 780 breast ultrasound images [16]. The original dataset has been reported to contain certain inconsistencies that could significantly impact model performance [17]. To ensure the validity and reliability of our evaluation, we exclude this dataset and utilize the remaining eleven 2D datasets along with all six 3D datasets for assessing our model’s performance. The image resolutions we used are 224×224 for 2D images and 64×64×64 for 3D images. Further details about the datasets can be found in section 3 of the Supplementary Information.

We use the MedMNIST v2 split training and validation sets to train and select hyperparameters and report the results of the test set. Accuracy (ACC) and Area Under the ROC Curve (AUC) are used as evaluation metrics to ensure a fair comparison with benchmark methods reported in the literature [14, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32]. All benchmark results are taken from their corresponding publications. To enhance the reliability of the results, we repeated the process three times with different random seeds and used the average value as the final performance of our model.

2.3 Overall Performance

Performance comparison of the proposed MTDL model with other state-of-the-art methods on the MedMNIST v2 dataset, in terms of AUC and ACC, is presented in Fig. 2 (detailed values refer to section 9 of the Supplementary Information). Two radar charts are used to show the performance comparison among different models across all 17 datasets for ACC and AUC respectively. As shown in the figure, the polygon corresponding to MTDL model covers the largest area and is situated at the outermost edge of the region occupied by the polygons of all the models, demonstrating its superior overall performance for medical image classification (Fig. 2a). Notably, MTDL demonstrates significant improvements over the second-best models in specific datasets:

• DermaMNIST: AUC improves from 0.937 to 0.962, and ACC improves from 0.780 to 0.836.
• RetinaMNIST: AUC improves from 0.773 to 0.874, and ACC improves from 0.568 to 0.655.
• OrganMNIST3D: AUC improves from 0.995 to 0.999, and ACC improves from 0.912 to 0.952.
• SynapseMNIST3D: AUC improves from 0.866 to 0.951, and ACC improves from 0.820 to 0.931.

Additionally, we compute the average AUC and ACC separately for 2D and 3D datasets, MTDL consistently outperforms all other models for 2D and 3D tasks (Fig. 2b). Specifically, for 2D datasets, MTDL achieves an average AUC of 0.956 and an ACC of 0.868, outperforming the second-best model, which achieves an average AUC of 0.943 and an ACC of 0.846. For 3D datasets, MTDL achieves an average AUC of 0.910 and an ACC of 0.855, compared to the second-best model’s average AUC of 0.901 and ACC of 0.832. A statistical significance analysis of the improvement is provided in Section 6.1 of the Supplementary Information.

Furthermore, we count the frequency of top performance for all models. As shown in the figure, MTDL can outperform all other models in AUC and ACC for 2D and 3D tasks (Fig. 2c). Specifically, for 2D tasks, MTDL achieves the highest AUC and ACC on six tasks, including RetinaMNIST (1,600 samples), DermaMNIST (10,015 samples), BloodMNIST (17,092 samples), OrganCMNIST (23,660 samples), OrganAMNIST (58,850 samples), and OCTMNIST (109,309 samples). This demonstrates its ability to perform effectively on prediction tasks of varying data scales. When considering the top-2 models, MTDL ranks within the top 2 in a frequency of 10/11 for AUC, 8/11 for ACC, and 8/11 for both AUC and ACC, which is significantly better than the second-best model, which ranks within top 2 in 5/11 for AUC, 4/11 for ACC, and 3/11 for both AUC and ACC. For 3D tasks, MTDL ranks best in both AUC and ACC for 2 tasks while no other model achieves the top rank for both metrics on any dataset. Moreover, MTDL ranks in the top 2 with a frequency of 6/6 for AUC, 6/6 for ACC, and 6/6 for both ACU and ACC, compared with the second-best model’s performance of 4/6 for AUC, 4/6 for ACC, and 2/6 for both AUC and ACC.

These results highlight the overall superiority of MTDL in comparison to other state-of-the-art models, demonstrating its effectiveness in handling both 2D and 3D medical image classification tasks.

2.4 Robustness Analysis Across Data Modality, Scale, and Task Type

To assess the robustness and generalizability of model performance, we divide the 17 datasets into groups based on data modality, data scale, and task type (refer to section 4 of the Supplementary Information), and then compare the average performance of all models within each group. To ensure a fair comparison, for each group, MTDL is evaluated only against models that reported results for all datasets in the respective group.

For data modality, we divide the datasets into four groups: Radiology (X-ray, CT, MRA), Microscopy (Pathology, Electron Microscope), Ophthalmology, and Dermatology. The performance comparison between MTDL and other models is shown in Fig. 3a. It can be seen that MTDL consistently outperforms all other models in both AUC and ACC across all groups. Specifically, MTDL achieves an average AUC (ACC) of 0.926 (0.875) for Radiology, compared to the second-best model’s performance of 0.903 (0.836). For Microscopy, MTDL obtains an average AUC (ACC) of 0.973 (0.890), while the second-best model achieves 0.942 (0.829). For Ophthalmology, MTDL gets an average AUC (ACC) of 0.932 (0.772), significantly surpassing the second-best model’s score of 0.869 (0.710). For Dermatology, MTDL obtains an average AUC (ACC) of 0.962 (0.836), compared to the second-best model’s performance of 0.937 (0.780). Notably, MTDL maintains an AUC above 0.930 for all four groups and an ACC above 0.835 for groups except Ophthalmology, the ACC for Ophthalmology is slightly smaller than other groups. We attribute this to the RetinaMNIST dataset within the Ophthalmology group since this dataset contains only 1600 samples. Despite this, MTDL still significantly outperforms other models for this group.

For data scale, we divide the datasets into four groups based on the sample size n of each dataset: G1 (n ⩽10K), G2 (10K < n ⩽50K), G3 (50K < n ⩽100K), and G4 (100K < n). Performance in terms of AUC and ACC for all models is presented in Fig. 3b. MTDL also ranks best in both metrics for all groups. Specifically, MTDL achieves an average AUC (ACC) of 0.930 (0.782) for G1, compared to the second-best model’s performance of 0.884 (0.761). For G2, MTDL scores 0.984 (0.890) compared to the second-best model’s score of 0.965 (0.876). In G3, both MTDL and the second-best model achieve an average AUC of 0.998, but MTDL has a slightly higher ACC (0.956 vs. 0.951). For G4, MTDL obtains 0.931 (0.869) compared to the second-best model’s 0.924 (0.852). MTDL can get an AUC exceeding 0.930 for all four groups and an ACC exceeding 0.860 for groups except G1, which is reasonable because bigger data usually leads to better performance.

For task type, we divide the datasets into four groups based on the number of classes n for each classification task: G1 (n=2), G2 (2< n ⩽5), G3 (5< n ⩽10), and G4 (10< n). The performance for all models is shown in Fig. 3c. MTDL again achieves the best overall performance. Specifically, MTDL achieves an average AUC (ACC) of 0.914 (0.909) for G1 compared to the second-best model’s performance of 0.915 (0.880). For G2, MTDL scores 0.872 (0.709), significantly surpassing the second-best model’s score of 0.823 (0.663). For G3, MTDL obtains 0.975 (0.866) compared to the second-best model’s 0.970 (0.851). In G4, MTDL attains an average AUC (ACC) of 0.993 (0.911) compared to the second-best model’s value of 0.990 (0.882). MTDL achieves strong performance for G1, G3, and G4, with AUC exceeding 0.910 and ACC exceeding 0.880. We think the slightly lower performance on G2 is due to the inclusion of the small-sized RetinaMNIST dataset within this group.

These results demonstrate the superiority, robustness, and generalizability of MTDL across various data scales, data modalities, and task types, indicating its great potential for medical image analysis. It is noteworthy that MTDL has only 0.56M parameters for 2D tasks and 0.75M parameters for 3D tasks, which is significantly smaller than models such as ResNet, GoogleNet, Vision Transformer (ViT), and MedViT. Despite its lightweight architecture, MTDL demonstrates exceptional performance.

2.5 Evaluation on Clinical Data

In MedMNIST v2, most of datasets are derived from clinical sources, i.e., human subjects treated in hospitals and medical centers, such as the German National Center for Tumor Diseases [33], Zhongshan Hospital Affiliated to Fudan University [14], and Guangzhou Women and Children’s Medical Center [34], among others. The majority of source datasets are simply processed through center-cropping and resizing to uniform dimensions for inclusion in MedMNIST v2. Consequently, models based on MedMNIST v2 are generally reliable.

To better understand the clinical applicability of our model for medical image analysis, we need to check the effects of image resizing process on model performance. We utilized the HAM10000 dataset, the original clinical dataset for DermaMNIST, and evaluated the performance of MTDL on it. HAM10000 consists of 10015 dermatoscopic images from different populations, including a representative collection of all important diagnostic categories in the realm of pigmented lesions. Over 50% of lesions in it have been confirmed by pathology, while the remaining cases are validated through either follow-up examinations, expert consensus, or in-vivo confocal microscopy [35].

The images in HAM10000 have the same size of 3×600×450. We center-crop the images to 3×450×450 and then resize into five resolutions: 3×35×35, 3×75×75, 3×150×150, 3×300×300, 3×450×450, with cubic spline interpolation. The performance of MTDL on these groups is shown in Table 1. As seen in the table, MTDL achieves improved performance as the image resolution increases. This indicates that MTDL is capable of extracting more detailed features from higher-resolution inputs, which makes it well-suited for clinical applications where high-resolution images are prevalent. Notably, the best performance is achieved on the largest image resolution (3×450×450), surpassing the results obtained on the DermaMNIST dataset. Specifically, the AUC improves from 0.962 to 0.973 and the ACC improves from 0.836 to 0.863. Even at the lowest resolution (3×35×35), MTDL achieves an AUC (ACC) of 0.943 (0.797), outperforming the best existing models’ performance of 0.937 (0.780). This highlights the robust lower-bound performance of MTDL across varying data resolutions, a critical attribute for addressing real-world clinical challenges.

2.6 Ablation Study

In our proposed MTDL model, the original image is decomposed into distinct orthogonal components, which are then concatenated to form a new composite image, serving as input to the CNN architecture. To evaluate the importance of the Hodge decomposition method, we performed an ablation study by replacing the decomposed images with the original images and denote the resulting model as ImgCNN. We compare the performance of MTDL and ImgCNN on five 2D datasets spanning different data scales, including RetinaMNIST (1,600 samples), PneumoniaMNIST (5856 samples), DermaMNIST (10,015 samples), OrganAMNIST (58,830 samples), and PathMNIST (107,180 samples). Additionally, the comparison extends to two 3D datasets: VesselMNIST3D (binary classification) and FractureMNIST3D (three-class classification). The results are summarized in Table. 2. As shown in the table, MTDL consistently outperforms ImgCNN across 2D and 3D tasks. Notably:

• For RetinaMNIST, AUC improves from 0.838 to 0.874, and ACC improves from 0.608 to 0.655.
• For DermaMNIST, AUC improves from 0.957 to 0.962, and ACC improves from 0.808 to 0.836.
• For VesselMNIST3D, AUC improves from 0.924 to 0.937, and ACC improves from 0.903 to 0.938.

These findings demonstrate the significant potential of the Hodge decomposition approach to enhance medical image representation, enabling improved performance across diverse datasets and classification tasks.

Methods
4.1 Topology-preserving Hodge Decomposition for Images

Hodge decomposition is a fundamental result in differential geometry and algebraic topology, specifically for the analysis of differential forms on Riemannian manifolds. Recently, a discrete topology-preserving Hodge decomposition for manifolds with boundaries on Cartesian grids has been introduced [12]. This method is particularly well-suited for image analysis, as images can be naturally treated as discrete manifolds with boundaries embedded in Cartesian grids.

Let M be an m-dimensional smooth, orientable, compact manifold with boundary ∂M, Ωk(M) represent the space of differential k-forms on M, and d denote the differential (exterior derivative) from k-forms to (k+1)-forms. A differential k-form ω is called closed if dω = 0 and exact if there exists a (k−1)-form ζ such that dζ = ω.

Given a Riemannian metric g on M, let ⋆ be the Hodge star operator that maps k-forms to (m−k)-forms and (·,·) denote the induced Hodge L2 inner product on Ωk(M). The codifferential δ : Ωk(M)→Ωk−1(M) is defined as

δ = (−1)m(k−1)+1⋆d⋆. (1)

A differential k-form ω is called coclosed if δω = 0, and coexact if there exists a (k+1)-form ζ such that δζ = ω. The operators d and δ satisfy the following relationship

(dω,η) = (ω,δη) + ∫∂M ω∧⋆η, (2)

where ω is a (k−1)-form, η is a k-form, and ∧ is the wedge product on differential forms. This implies that d and δ are adjoint if M is a closed manifold, i.e., ∂M = ∅.

The Hodge Laplacian for differential forms is defined as

∆ = dδ+δd. (3)

The Laplacian operator maps k-forms to k-forms. The kernel of ∆ is called the space of harmonic forms. We denote by Hk(M) the space of harmonic k-forms and by Hk(M) the space of k-forms that are both closed and coclosed. We have Hk(M) ⊂ Hk(M).

When M is a closed manifold, i.e., a compact manifold without boundary. The standard Hodge decomposition [37] states that

Ωk(M) = dΩk−1(M) ⊕ δΩk+1(M) ⊕ Hk(M), (4)

where the adjointness of d and δ ensure that these three subspaces are orthogonal with respect to the Hodge L2 inner product.

When M is a manifold with non-empty boundary, the operators d and δ are generally not adjoint, as noted in (2). To ensure their adjointness and consequently achieve an orthogonal decomposition of differential forms, appropriate boundary conditions must be imposed.

Two most commonly used boundary conditions are the normal (Dirichlet) boundary condition and the tangential (Neumann) boundary condition. These conditions define the following subspaces,

Ωk(M) = {ω ∈ Ωk(M) | ω| = 0}, Ωk(M) = {ω ∈ Ωk(M) | ⋆ω| = 0}. (5)
n ∂M t ∂M

The forms in Ωk(M) and Ωk(M) are called normal and tangential respectively.

The Hodge-Morrey decomposition [38] states that

Ωk(M) = dΩk−1(M) ⊕ δΩk+1(M) ⊕ Hk(M). (6)
n t

The exterior derivative d preserves the normal boundary condition and the codifferential δ preserves the tangential boundary condition. As a result, any k-form can be decomposed as the sum of an exact normal form, a coexact tangential form, and a harmonic form that is both closed and coclosed.

ω = dα + δγ + η, (7)
n t

where ω ∈ Ωk(M), α n ∈ Ωk n −1(M), γ t ∈ Ωk t + 1(M), η ∈ Hk(M). When we focus on the compact manifold in Euclidean spaces, the third term Hk(M) in (6) can be further decomposed into three orthogonal components [39], resulting in a five-component decomposition. A more detailed description of the Hodge decomposition can be found in section 1 of the Supplementary Information.

Here we focus on 2D and 3D Cartesian grids, as medical images are typically in these dimensions. The discrete manifold M in Cartesian grids can be given as a sublevel set of a level set function on the grid. We employ the strategy in [40] to determine the boundary of M for two boundary conditions. For normal boundary condition, cells with at least one vertex inside M are included, while for tangential boundary condition, cells with at least one vertex of their dual cells inside M are included. The resulting sets of cells are referred to as the normal support for the normal boundary condition and the tangential support for the tangential boundary condition. These supports can be seen as discrete versions of the manifold with boundary. The discrete differential operators, including exterior derivative, Hodge star, codifferential, and Laplacian can be defined for two boundary conditions accordingly. As in the continuous case, the Kernels of these discrete Laplacians are fully determined by the topology of M. Specifically, the dimension of kerL k,n equals the Betti number β , while the dimension of kerL equals β . The Betti number β quantifies the number of k-dimensional topological features in M: β represents the number of connected components, β the number of loops, and β the number of voids.

Fig. 4 illustrates an example demonstrating the topology-preserving property of the discrete Laplacian. As shown in the figure, a blood cell image is represented as a discrete manifold under boundary conditions (Fig. 4a). This manifold exhibits three distinct loop structures, resulting in a Betti number β of 3. We compute the Laplacian L under the normal boundary condition, and the eigenvectors corresponding to the three zero eigenvalues are displayed. These eigenvectors align precisely with the three loops present in the manifold (Fig. 4b). In Fig. 4c, a vector field (1-form) on the manifold is first derived from the image using the flow-based method described in Supplementary Information. This vector field is subsequently decomposed into three orthogonal components: the curl-free, divergence-free, and harmonic parts. The harmonic component represents the global topological structure of the underlying manifold, whereas the normal and tangential components characterize distinct aspects of the local information. Specifically, the textures of the normal and tangential components exhibit an approximately perpendicular relationship, and the harmonic component appears smoother compared to the other two components.

4.2 CNN Architecture

The CNN we used is based on the Transformer encoder architecture by adding a maxpooling operation and replacing the multihead attention and feedforward layers by convolution operations (Fig. 1e). As illustrated, the decomposed image x is first processed through an initialization block, which consists of a convolutional layer, a batch normalization operation, and a nonlinear ReLU activation function. The initialized image is then passed through a sequence of simple-CNN layers to extract hierarchical features, where each simple-CNN layer is composed of two convolutional blocks followed by a pooling operation. Finally, the extracted features are spatially averaged and fed into a multilayer perceptron (MLP) for classification. More details can be found in section 5.1 of the Supplementary Information.

4.3 Model Implementation Detail

In our implementation, each image is considered as a scalar field on the vertices of a standard Cartesian grid. The discrete manifold is generated through segmentation, where foreground pixels are extracted via thresholding to remove background regions (refer to section 6.3 of the Supplementary Information). We use the grid vertices, edges, faces, and cubes to construct the differential operators and projection operators. We use the gradient-based method (refer to section 2 of the Supplementary Information) to construct the vector field from the images and perform the Hodge decomposition.

The proposed MTDL model is implemented using PyTorch [41] and evaluated on an NVIDIA Tesla V100S GPU. For 2D datasets, the batch size and learning rate are set to 64 and 10−3, respectively, across all tasks. The training process spans 30 epochs for tasks with a sample size smaller than 100,000 and 10 epochs for tasks with a sample size exceeding this threshold. The number of layers in the model is adapted based on the data distribution. For the majority of tasks, a 5-layer structure is employed, with detailed configurations provided in section 5.2 of the Supplementary Information. The hidden channel dimension C and head number h are set to 72 and 4, with the group number g configured as 1 for grayscale images and 3 for colored images.

For 3D datasets, the batch size and learning rate are set to 16 and 10−3, respectively, for all tasks. The training process involves 10 epochs for the FractureMNIST dataset and 20 epochs for the remaining datasets. Similar to the 2D case, the number of layers is determined by the data distribution, with further details available in section 5.2 of the Supplementary Information. The hidden channel dimension C and head number h are set to 64 and 4, with a group number g of 1 since all the 3D images are grayscale.

The model is optimized using the AdamW optimizer [42] with a weight decay of 10−5, and a one-cycle learning rate scheduler employed [43]. For the ChestMNIST2D task, a multi-label classification problem, the Binary Cross-Entropy with Logits is used as the loss function, while Cross-Entropy Loss is applied for all other tasks.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The manuscript introduces Manifold Topological Deep Learning (MTDL) and demonstrates competitive performance on MedMNIST v2. The idea of extending topological deep learning from combinatorial data to differentiable manifolds via Hodge decomposition is interesting. However, the paper suffers from issues in methodological clarity, fairness of comparisons, rigor of evaluation, and substantiation of interpretability claims. The paper requires major revisions before it can be considered further.\nStrengths:\n1. The manuscript successfully proposes a novel framework, MTDL, extending TDL from combinatorial data (like graphs and point clouds) to differentiable manifold data, such as images.\n2. The use of Hodge decomposition to factor an image's vector field into three orthogonal components (curl-free, divergence-free, and harmonic parts) is mathematically elegant. This method captures distinct geometric and topological features, offering the potential for enhanced interpretability.\n3. MTDL achieves competitive or superior results across a diverse range of tasks categorized by dimensionality (2D vs. 3D), data modality, data scale, and task type. Moreover, the model is highly efficient, with a lightweight architecture containing only 0.56M and 0.75M parameters for 2D and 3D tasks, respectively.\nMajor Concerns:\n1. Figure 1a contradicts the definitions of boundary conditions. According to the definition, the normal boundary condition (NBC) should generate a larger support set than the tangential boundary condition (TBC), but Figure 1a shows the opposite. While supplementary Figure 2 is consistent with the definition.\n2. The authors state that topological deep learning models “offer greater interpretability.” However, the paper does not provide any analysis or demonstration of interpretability. The new representation of images as vector fields is introduced, but its meaning and diagnostic value are not justified.\n3. In the ablation study, the authors replace the decomposed images (curl-free, divergence-free, and harmonic parts) with the original images. It remains unclear which of the three modalities contributes most to performance. The study does not clarify whether each component alone is weaker than the original image but the combination yields a gain.\n4. According to Figure 4c, curl-free and divergence-free components appear very similar to each other and to the raw vector field. For images, especially medical images where most pixels are far from the boundary, NBC and TBC differ only at edge regions. Consequently, curl-free and divergence-free maps are nearly indistinguishable. This raises doubts about the necessity of including all three channels. Maybe simply augmenting the original image with the harmonic part would suffice? An ablation explicitly comparing the role of these three channels is needed.\n5. MTDL concatenates three components, thereby increasing the number of input channels and the effective model capacity. The ImgCNN ablation lacks sufficient description and parameter counts are not reported. It is essential to ensure the total number of parameters is equal. Otherwise, the performance gain cannot be attributed specifically to Hodge decomposition.\n6. The BIG Laplacian is introduced without its full name or a clear derivation. The rationale for simplifying the Laplacian into this form, and its numerical and geometric consequences, is insufficiently discussed. It remains unclear whether the BIG Laplacian preserves orthogonality, which is a central property of Hodge decomposition. The conditions under which it is equivalent or a good approximation to the full Laplacian should be explicitly analyzed.\n7. The paper emphasizes that MTDL is a “lightweight” model based on small parameter counts. However, the Hodge decomposition requires additional computation. The extra computational cost is not reported.\n8. MTDL shows competitive results but is not consistently superior. The statistical significance of improvements is not demonstrated. In addition, the comparison set relies on some outdated baselines. More recent state-of-the-art methods from the past two years should be included for the comparison to have empirical relevance.\n9. The so-called “Transformer-encoder-induced convolution layers” replace multi-head attention and the feed-forward network with group convolution and 1×1 convolution. This is essentially a CNN residual block, not a Transformer. The terminology is misleading. Either rename the block to reflect its actual nature, or implement true attention for comparison.\n10. Operator construction lacks transparency. Equations (8)–(11) define several operators, but their explicit matrix construction is not shown. Providing a concrete example using a small toy data would improve reproducibility and make the work accessible to readers unfamiliar with TDL.", "reviewer_2": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_3": "Co-Review of NCOMMS-25-20703-T; Nature Communications: by Zakeri A & Frangi A F\nTitle: Manifold Topological Deep Learning for Biomedical Data\nSummary:\nThe manuscript introduces Manifold Topological Deep Learning (MTDL), which represents images as discrete manifolds, constructs a vector field on that manifold, applies a discrete Hodge decomposition into curl-free / divergence-free / harmonic components, concatenates those channels, and feeds the result into a CNN. The model is evaluated on MedMNIST v2 (11 2D + 6 3D datasets) and reported to outperform prior methods.\nThe idea of converting images to vector fields and apply a discrete Hodge decomposition before feeding the components to a CNN model for classification is promising. However, the paper needs further clarification on methodological details, additional ablation studies, and rigorous evaluation. Detailed comments below outline specific improvements needed before publication.\nOur view is that the manuscript needs a Major Revision:\nMethodology:\n1. Explain how the harmonic component h is computed.\n2. Provide Pseudo-code for the proposed methodological pipeline.\n3. Clarify how the manifold M and its boundary ∂M are defined for each dataset. In practical settings, computing normal and tangential boundary conditions requires a segmentation mask or at least a binary indicator of the region of interest. In Section 4.3.1, the segmentation step is described only as thresholding but the choice of threshold (fixed, adaptive, or dataset-specific), preprocessing (e.g., normalization or filtering), and handling of ambiguous or noisy background regions can significantly influence the resulting manifold boundary ∂M, which in turn affects the computed boundary conditions and the subsequent Hodge decomposition. I recommend the authors (1) explicitly describe the segmentation procedure per dataset, (2) discuss its influence on the topological representation, and (3) provide an ablation or sensitivity analysis to assess robustness to segmentation variations in the results.\n4. The paper would strongly benefit from an explicit pseudocode summarizing the full implementation workflow, including how boundary conditions, discrete operators, and vector-field decompositions are computed. Although the Supplementary Material gives substantial mathematical detail, a procedural description (with solver specifications and preprocessing steps) is essential for reproducibility.\nResults:\n5. Vector-field construction is under-specified in the main text. The paper mentions several methods and refers to the Supplementary Information, but the main text does not make clear which construction is used for each experiment, what hyperparameters were used, and how sensitive results are to that choice. Since the vector field is the core representation, an ablation is necessary: different vector-field constructions (flow-based, gradient-based, others) and their effect on final performance.\n6. Discretization errors and orthogonality checks:\nProvide results showing the discrete decomposition is numerically orthogonal.\nReport: [see formula in attachment]\nProvide these for representative images or in a table for a dataset with mean ± std.\n7. The paper concatenates components as extra channels. Provide an ablation showing which components contribute most to performance (curl-free / divergence-free / harmonic).\n8. Comparisons / fairness of baselines:\n9. The paper reports substantial improvements across multiple MedMNIST tasks; however, it is unclear whether the baseline models were trained and tuned under identical experimental conditions (e.g., data augmentation, optimizer settings, and hyperparameter search strategies). Please clarify whether the baseline results were all adopted from the original publications or reproduced by the authors. If any were re-implemented, details on the hyperparameter parity and training cost should be provided. Without this, it's hard to know whether gains derive from topology or from training differences.\n10. Statistical analysis: While the reported results on MedMNIST v2 are promising, the paper lacks statistical analysis. Since many datasets are small and stochastic training introduces non-negligible variance, it is important to report mean ± standard deviation across multiple runs and perform statistical significance testing. A paired t-test or Wilcoxon test within each dataset, and a Friedman test with Nemenyi post-hoc across datasets, would clarify whether the observed improvements are statistically significant rather than due to random variation. Reporting confidence intervals and effect sizes (e.g., Cohen’s d or Cliff’s δ) would further strengthen the empirical evidence.\n11. In Figure 3, the range of the AUC or ACC values on the x-axis should be consistent across all plots to facilitate clearer and quicker visual comparison between results.\n12. Provide preprocessing time per sample (vector field + decomposition), GPU memory/CPU usage, and training time per epoch for representative datasets (2D and 3D).\nDiscussions:\n13. Invariance/equivariance: discuss which geometric transforms of the input (translation, rotation, scaling) preserve each decomposed channel and whether the CNN architecture exploits or breaks these invariances.\n14. Physical/semantic interpretation: Can we link the decomposed components to expected image structure per modality (e.g., curl-free ~ edge/gradient information; divergence-free ~ texture swirl)? so readers can interpret classifier reliance on components.\n15. How the method might generalize to real-world biomedical images, e.g., high-resolution CT/MRI, anisotropic voxels, or dynamic sequences? MedMNIST is small and standardized; clinical data might pose challenges.\n16. A discussion (or ideally an experiment) comparing MTDL to surface-based or mesh-based features could be helpful. This might help readers understand when grid-based manifold representation is preferred vs direct mesh representation or potential ways these approaches could complement each other in biomedical imaging applications.\n17. The paper could benefit from explicit discussion of scenarios where MTDL might fail."}, "author_response": "Response to Reviewers\nAnswers to Reviewer 1’ Comments\nThe manuscript introduces Manifold Topological Deep Learning (MTDL) and demon- strates competitive performance on MedMNIST v2. The idea of extending topological deep learning from combinatorial data to differentiable manifolds via Hodge decomposition is interesting. However, the paper suffers from issues in methodological clarity, fairness of comparisons, rigor of evaluation, and substantiation of interpretability claims. The paper requires major revisions before it can be considered further.\nStrengths:\n1. The manuscript successfully proposes a novel framework, MTDL, extending TDL from combinatorial data (like graphs and point clouds) to differentiable manifold data, such as images.\n2. The use of Hodge decomposition to factor an image’s vector field into three orthog- onal components (curl-free, divergence-free, and harmonic parts) is mathematically elegant. This method captures distinct geometric and topological features, offering the potential for enhanced interpretability.\n3. MTDL achieves competitive or superior results across a diverse range of tasks catego- rized by dimensionality (2D vs. 3D), data modality, data scale, and task type. Moreover, the model is highly efficient, with a lightweight architecture containing only 0.56M and 0.75M parameters for 2D and 3D tasks, respectively.\nAnswer: We thank the positive comments and insights of the reviewer.\nMajor Concerns:\nComment : 1) Figure 1a contradicts the definitions of boundary conditions. According to the definition, the normal boundary condition (NBC) should generate a larger support set than the tangential boundary condition (TBC), but Figure 1a shows the opposite. While supplementary Figure 2 is consistent with the definition.\nAnswer: We thank the reviewer for the careful review. You are correct that the left part should correspond to the tangential boundary condition, while the right part should correspond to the normal boundary condition. We have revised the figure accordingly.\nComment : 2) The authors state that topological deep learning models “offer greater interpretability.” However, the paper does not provide any analysis or demonstration of inter- pretability. The new representation of images as vector fields is introduced, but its meaning and diagnostic value are not justified.\nAnswer: We thank the reviewer for the insightful comments. The mechanistic in- terpretability of our model arises from the explicit decomposition of the image information into three mathematically meaningful components, rather than relying on the original im- ages. In our model, each image is represented as a vector field, which is then decomposed into a curl-free component, divergence-free component, and a harmonic component. These three components have clear geometric and physical meanings. Specifically, a vector field is curl-free means it does not have the local rotational behavior. A divergence-free vector field has neither sources nor sinks. Source regions correspond to positive divergence, where local intensity or structural patterns spread outward, while sink regions correspond to nega- tive divergence, where patterns converge inward. The harmonic component is both curl-free and divergence-free, representing the smoothest global information associated with topological traits.\nFor the diagnostic value in medical image analysis, the curl-free component can be used to capture the lesion boundaries and margins, as the directions of the strongest intensity transition usually correspond to tissue boundaries. For example, in mammography, malignant lesions often show sharper and more irregular boundaries compared with the benign lesions. Therefore, the curl-free component can be used to distinguish the malignant from benign cases. The divergence-free component captures the vorticity and local rotating behavior in the vector field. Consequently, it can be used to capture the irregular textures, such as the swirling patterns and the heterogeneous internal textures due to necrosis and fibrosis. Especially, this will be useful for malignant tumors since the swirling textures are common in malignant tumors such as breast cancer, brain tumor, and hepatocellular carcinoma. The harmonic component reflects the global topological properties (cycles) of the image and provides robust descriptors insensitive to noises and small perturbations. For a more comprehensive model interpretability about how the decomposed components relate to visual patterns such as edges, textures, and global structures, additional interpretability analysis would be required. We regard this as an important direction for future work.\nWe have added these discussions to section 6.4 of the revised supplementary information.\nComment : 3) In the ablation study, the authors replace the decomposed images (curl- free, divergence-free, and harmonic parts) with the original images. It remains unclear which of the three modalities contributes most to performance. The study does not clarify whether each component alone is weaker than the original image but the combination yields a gain.\nAnswer: We thank the reviewer for the constructive comments. We have added an ablation study in which we evaluate the model performance by using different combinations of the components of the decomposition, the results have been included in section 6.2 of the revised supplementary information.\nFor the reviewer’s convenience, we also summarize these results in Table S1. In this analysis, we consider five 2D datasets of varying scales to provide a comprehensive evaluation:\nRetinaMNIST (1,600 samples), PneumoniaMNIST (5856 samples), DermaMNIST (10,015 samples), OrganAMNIST (58,830 samples), and PathMNIST (107,180 samples). We have\nTable S1: Performance (ACC & AUC) Comparison between the original images and decomposition images.\nCurl: the curl-free component, Div: the divergence-free component, Har: the harmonic component. Img: the original image.\nRetina Pneumonia Derma OrganA Path\nCurl 0.628 0.910 0.820 0.953 0.914\n0.867 0.986 0.957 0.998 0.995\nDiv 0.610 0.885 0.725 0.887 0.700\n0.854 0.967 0.890 0.990 0.938\nHar 0.485 0.793 0.633 0.887 0.847\n0.685 0.921 0.788 0.989 0.976\nCur+Div 0.635 0.910 0.827 0.956 0.914\n0.868 0.982 0.963 0.998 0.995\nCur+Har 0.635 0.911 0.826 0.955 0.918\n0.867 0.984 0.962 0.998 0.995\nDiv+Har 0.615 0.829 0.731 0.906 0.784\n0.843 0.960 0.894 0.993 0.963\nCur+Div+Har 0.655 0.910 0.836 0.956 0.920\n0.874 0.986 0.962 0.998 0.996\nImg 0.608 0.885 0.808 0.955 0.902\n0.838 0.978 0.957 0.998 0.987\nImg+Curl 0.625 0.886 0.821 0.955 0.905\n0.848 0.977 0.958 0.998 0.993\nImg+Div 0.630 0.872 0.810 0.956 0.902\n0.834 0.972 0.959 0.998 0.993\nImg+Har 0.630 0.864 0.810 0.959 0.914\n0.835 0.970 0.958 0.998 0.995\nthe following observations:\n1. The curl-free component contributes most to the model performance compared with the divergence-free and harmonic parts.\n2. Using only the curl-free component can achieve performance comparable to or better than using the original images.\n3. The combination of all three components consistently gives the best overall performance for all tested datasets, which is reasonable since each component captures distinct and complementary information, and their combination provides a more complete represen- tation of the image features.\n4. Augmenting the original images with any single component of the decomposition (curl- free, divergence-free, or harmonic) improves performance compared with using only the original images on several datasets, such as RetinaMNIST, DermaMNIST, and PathMNIST. These findings demonstrate that each component captures specific features that cannot be directly extracted from the original images.\nComment : 4) According to Figure 4c, curl-free and divergence-free components appear very similar to each other and to the raw vector field. For images, especially medical images where most pixels are far from the boundary, NBC and TBC differ only at edge regions. Consequently, curl-free and divergence-free maps are nearly indistinguishable. This raises doubts about the necessity of including all three channels. Maybe simply augmenting the original image with the harmonic part would suffice? An ablation explicitly comparing the role of these three channels is needed.\nAnswer: We thank the reviewer for the constructive comments. In response to Com- ment 3), we have included an ablation study that evaluates the model’s performance using different combinations of the components in the decomposition. As shown in Table S1, the curl-free component contributes the most to the model performance, and the combination of all three components yields the best overall performance across all evaluated datasets.\nComment : 5) MTDL concatenates three components, thereby increasing the number of input channels and the effective model capacity. The ImgCNN ablation lacks sufficient description and parameter counts are not reported. It is essential to ensure the total number of parameters is equal. Otherwise, the performance gain cannot be attributed specifically to Hodge decomposition.\nAnswer: We thank the reviewer for pointing out this. It is true that the MTDL model receives three input channels for each original channel, as the Hodge decomposition automat- ically divides an image into three components. Consequently, for an image with n channels, MTDL uses 3n channels while ImgCNN uses n channels. However, we want to emphasize that this difference affects only the initial convolution layer. As shown in Figure 1e of the main text, both MTDL and ImgCNN employ an initialization convolution layer with 72 out- put channels regardless of the number of input channels. Consequently, after this first layer, the rest of the network architecture is identical for both models, and all hyperparameters and training settings are exactly the same. Therefore, the downstream model capacity remains matched.\nMoreover, Table S1 shows that using only the curl-free component already yields con- sistently better performance than using the original images on RetinaMNIST, Pneumoni- aMNIST, DermaMNIST, and PathMNIST datasets. Note that the curl-free component has the same channel number with the original image. This observation further demonstrates that the performance improvement arises from the informative features introduced by the Hodge decomposition rather than from increased input channel.\nTo have a direct comparison, we conducted an additional ablation study in which the original images were augmented by either zero-padding or by adding blurred low- and high- frequency components such that the total number of input channels matched that of the MTDL representation. Specifically, the low-frequency and high-frequency components were obtained using Gaussian blurring and subsequently combined with the original images. These aug- mented images were then fed into the same CNN architecture used for MTDL. The results are shown in Table S2. It can be seen that the original images actually achieve compara- ble, and in most cases better performance than the zero-padded augmented images. Although augmenting the inputs with low- and high- frequency components yields better performance than simply adding zero-valued channels, the improvement remains substantially below that achieved by the full MTDL representation. This confirms that simply increasing the number of input channels does not improve performance and that the gains achieved by MTDL is from the meaningful decomposition provided by the Hodge decomposition rather than from the increased channel count.\nTable S2: Performance (ACC & AUC) Comparison between the original images and the augmented im- ages. Img: the original image. Img+zero: augmenting the original images with zero-padded channels.\nImg+low+high: augmenting the original images with both low- and high-frequency components.\nRetina Pneumonia Derma OrganA Path\nImg 0.608 0.885 0.808 0.955 0.902\n0.838 0.978 0.957 0.998 0.987\nImg+zero 0.603 0.880 0.801 0.958 0.905\n0.823 0.979 0.955 0.998 0.978\nImg+low+high 0.612 0.888 0.817 0.958 0.912\n0.837 0.977 0.957 0.998 0.986\nMTDL 0.655 0.910 0.836 0.956 0.920\n0.874 0.986 0.962 0.998 0.996\nComment : 6) The BIG Laplacian is introduced without its full name or a clear derivation. The rationale for simplifying the Laplacian into this form, and its numerical and geometric consequences, is insufficiently discussed. It remains unclear whether the BIG Laplacian preserves orthogonality, which is a central property of Hodge decomposition. The conditions under which it is equivalent or a good approximation to the full Laplacian should be explicitly analyzed.\nAnswer: We thank the reviewer for the insightful comments. We have added the full name of the Boundary-Induced Graph (BIG) Laplacian and cited the corresponding reference accordingly. The BIG Laplacian can be derived by replacing the discrete Hodge star with the identity matrix. In the discrete setting, the Hodge star is a diagonal matrix whose entries correspond to the ratio between the volumes of dual cells and primal cells. As a result, only boundary cells yield non-unit diagonal entries, whereas interior cells have value 1. As the grid resolution increases, the number of interior cells becomes dramatically larger than that of boundary cells. Therefore, the identity matrix is an efficient approximation to the Hodge star, and thus the BIG Laplacian serves as an efficient approximation to the Hodge Laplacian.\nIn our implementation, we compute the Hodge decomposition using a direct approach.\nSpecifically, for the decomposition\nω = dα+δβ +h\nwe firstly compute the potentials α and β, and then obtain the harmonic component by\nh = ω −dα−δβ\nActually, this direct approach does not guarantee strict orthogonality of the three components at the discrete level, and orthogonality is recovered only in the continuous limit.\nDespite the lack of exact mathematical orthogonality in the discrete case, the three com- ponents exhibit strong numerical orthogonality in practice. To quantify this, we computed the normalized pairwise inner products between the curl-free, divergence-free, and harmonic components across five representative datasets: RetinaMNIST (1,600 samples), Pneumo- niaMNIST (5,856 samples), DermaMNIST (10,015 samples), OrganAMNIST (58,830 sam- ples), and PathMNIST (107,180 samples). The normalized inner product for two components α,β is defined as\n| < α,β > |\n(α,β) = .\n||α|| ||β||\n2 2\nWe report the mean and standard deviation of these values across all samples for each dataset,\nwhich is shown in Table S3. The results indicate that the three components exhibit strong numerical orthogonality across all evaluated datasets.\nExact orthogonality is not easily achievable in discrete Hodge decompositions because the discrete differential operators generally do not satisfy all properties of the continuous theory.\nHowever, practical applications usually focus on the numerical approximation of orthogonal- ity. Since our numerical errors of orthogonality are near zero, the three components already provide independent and complementary features. For machine learning applications, this level of numerical orthogonality is entirely sufficient, as it ensures that the extracted features are non-redundant and do not interfere with each other.\nWe have added these discussions to section 1.3 of the revised supplementary information.\nTable S3: Numerical analysis of the orthogonality among the three components in the Hodge decomposition\nacross the RetinaMNIST, PneumoniaMNIST, DermaMNIST, OrganAMNIST, and PathMNIST datasets.\nThe mean and standard deviation values of the normalized inner product of all samples within each dataset\nare reported.\nMethod Retina Pneumonia Derma OrganA Path\nmean std mean std mean std mean std mean std\n(dα,δβ) 7e-20 4e-17 2e-19 2e-17 4e-20 3e-17 8e-20 3e-17 7e-20 2e-17\n(dα,h) 2e-13 9e-14 3e-14 2e-14 6e-14 1e-14 3e-14 3e-14 2e-14 1e-14\n(δβ,h) 2e-17 5e-17 2e-18 8e-17 7e-18 2e-16 5e-18 8e-17 2e-17 6e-17\nComment : 7) The paper emphasizes that MTDL is a “lightweight” model based on small parameter counts. However, the Hodge decomposition requires additional computation. The extra computational cost is not reported.\nAnswer: We thank the reviewer for pointing out this. In our experiments, all 2d\nimages are n×224×224, with n = 3 for color images and n = 1 for gray images. All 3d\nimages are grayscale with a size 64×64×64. The Hodge decomposition is computationally\nefficient: 50 grayscale 2d images can be completed within 2 minutes, 50 color 2d images\ncan be completed within 4 minutes, and 50 3d images can be completed within 30 minutes.\nWe have added the running time for the Hodge decomposition in section 7 of the revised\nsupplementary information.\nComment : 8) MTDL shows competitive results but is not consistently superior. The\nstatistical significance of improvements is not demonstrated. In addition, the comparison set\nrelies on some outdated baselines. More recent state-of-the-art methods from the past two\nyears should be included for the comparison to have empirical relevance.\nAnswer: We thank the reviewer for the constructive comments.\nFor statistical significance analysis of improvement, the baseline methods report only\nsingle-run results. Consequently, we conducted the Friedman test followed by the Nemenyi\npost-hoc test. This procedure is standard in multi-dataset benchmark evaluations, particu-\nlarly in settings where only single-run baseline results are available. Specifically, we used the\nscipy.stats.friedmanchisquare package to compute the p-values on 11 MedMNIST2D datasets\nand 6 MedMNIST3D datasets, respectively. For the 2D evaluation, we included only those\nbaseline models that reported results on all 11 MedMNIST2D datasets. Similarly, for the 3D\nevaluation, we restricted the comparison to baselines that provided results on all 6 MedM-\nNIST3D datasets. The resulting p-values are 1.38e-14 (2D) and 5.32e-10 (3D), indicating\nsignificant performance differences among the compared models. Therefore, we provide Crit-\nical Difference (CD) diagrams for both the 2D and 3D cases. The CD value is the minimum\nrank difference required for two models to be considered significantly different under the Ne-\nmenyi test. As shown in Figure S1, our MTDL model achieves the best overall performance\nwith the lowest average rank among all compared methods. Under the computed CD value,\nMTDL is significantly better than most compared methods for both 2D and 3D cases. This\nanalysis has been added to section 6.1 of the revised supplementary information.\nFor comparison with recent state-of-the-art methods, our manuscript was originally sub-\nmitted early this year, and at that time we included all publicly available models evaluated\non this database. To provide an updated and comprehensive comparison, we conducted a\nthorough literature search for works published in 2025 that also utilized this database. We\nsummarize these newly identified results in Table S4 to present an up-to-date benchmark\ncomparison. It can be seen that our model achieves the best overall performance among all\nthese methods. Notably, many of the recently published models are evaluated on only a subset\nrather than the full set of datasets while our model is evaluated on all the datasets. These\nnew results have been added into section 9 of the revised supplementary information.\n[1] Halder, A., Dalal, A., Gharami, S., Wozniak, M., Ijaz, M. F., & Singh, P. K.\nFigure S1: The Critical Difference (CD) diagrams illustrating the performance of our model across 11\nMedMNIST2D datasets and 6 MedMNIST3D datasets.\nTable S4: Performance comparison between our model and recent models on the MedMNIST2D dataset.\nThe benchmark results are taken from the original papers. The best results are in bold. “-” indicates that\nthe corresponding metric was not reported.\nMethods Path Chest Derma OCT\nAUC ACC AUC ACC AUC ACC AUC ACC\nFRDEM [1] - - - - - 0.783 - -\nAIGC [2] - - - 0.83 - -\nMedKAFormer [3] 0.985 0.882 - - 0.919 0.748 0.974 0.810\nLPT [4] - 0.941 - - - 0.787 - 0.832\nSGDGCN [5] 0.987 0.907 0.928 0.854 0.919 0.771 0.933 0.783\nPSNAS-Net-ave [6] - - - - - 0.772 - -\nSRE-CONV [7] - 0.822 - 0.948 - 0.757 - 0.768\nIMET [8] - - - - - - 0.960 0.803\nNQNN [9] - 0.752 - - - 0.827 - 0.591\nMTDL 0.996 0.920 0.793 0.948 0.962 0.836 0.989 0.888\nMethods Pneumonia Retina Blood Tissue\nAUC ACC AUC ACC AUC ACC AUC ACC\nMedKAFormer [3] 0.961 0.904 0.729 0.535 0.998 0.965 - -\nLPT [4] - - - 0.575 - 0.967 - -\nSGDGCN [5] 0.963 0.918 0.782 0.550 0.994 0.940 0.801 0.613\nPSNAS-Net-ave [6] - 0.899 - - - 0.956 - -\nSRE-CONV [7] - 0.906 - 0.523 - 0.960 - 0.712\nIMET [8] 0.895 0.902 - - - - - -\nNQNN [9] - - - - - 0.737 - -\nMTDL 0.986 0.910 0.874 0.655 0.999 0.988 0.945 0.721\nMethods OrganA OrganC OrganS\nAUC ACC AUC ACC AUC ACC\nAIGC [2] - 0.890 - - - -\nMedKAFormer [3] 0.995 0.922 0.991 0.899 0.986 0.788\nLPT [4] - - - 0.906 - 0.796\nSGDGCN [5] 0.977 0.850 0.989 0.880 0.969 0.784\nPSNAS-Net-ave [6] - - - 0.913 - -\nSRE-CONV [7] 0.770 0.769 - - - -\nNQNN [9] - 0.814 - 0.799 - -\nMTDL 0.998 0.956 0.996 0.928 0.978 0.809\n(2025). A fuzzy rank-based deep ensemble methodology for multi-class skin cancer classifica- tion. Scientific Reports, 15(1), 6268.\n[2] Wu, B., Huang, J., & Duan, Q. (2025). Real-time intelligent healthcare enabled by federated digital twins with aoi optimization. IEEE Network.\n[3] Wang, G., Zhu, Q., Song, C., Wei, B., & Li, S. (2025). MedKAFormer: When Kolmogorov-Arnold Theorem Meets Vision Transformer for Medical Image Representation. IEEE Journal of Biomedical and Health Informatics.\n[4] Bai, Y., Bai, L., Yang, X., & Liang, J. (2025). Label-Semantic-Based Prompt Tuning for Vision Transformer Adaptation in Medical Image Analysis. IEEE Transactions on Circuits and Systems for Video Technology.\n[5] Singh, A., Eising, C., Denny, P., & van de Ven, P. (2025). Improving GNNs for Image Classification: Addressing Homophily Challenges. IEEE Open Journal of the Computer Society.\n[6] Zechen, Z., Xuelei, H., Fengjun, Z., & Xiaowei, H. (2025). PSNAS-Net: Hybrid gradient-physical optimizationfor efficient neural architecture search in customized medical imaging analysis. Expert Systems with Applications, 288, 128155.\n[7] Du, Y., Zhang, J., Zeevi, T., Dvornek, N. C., & Onofrey, J. A. (2025, April). Sre-conv: Symmetric rotation equivariant convolution for biomedical image classification. In 2025 IEEE 22nd International Symposium on Biomedical Imaging (ISBI) (pp. 1-5). IEEE.\n[8]Singh, R., &Guggilam, S.(2025). IterativeMisclassificationErrorTraining(IMET): An Optimized Neural Network Training Technique for Image Classification. IEEE Access.\n[9] Rahman, M., & Zhuang, J. (2025, September). NQNN: Noise-Aware Quantum Neural Networks for Medical Image Classification. In International Conference on Medical Image Computing and Computer-Assisted Intervention (pp. 433-442). Cham: Springer Nature Switzerland.\nComment : 9) The so-called “Transformer-encoder-induced convolution layers” re- place multi-head attention and the feed-forward network with group convolution and 1×1 convolution. This is essentially a CNN residual block, not a Transformer. The terminol- ogy is misleading. Either rename the block to reflect its actual nature, or implement true attention for comparison.\nAnswer: We thank the reviewer for pointing this out. We have renamed the block as a simple-CNN layer in the revised manuscript.\nComment : 10) Operator construction lacks transparency. Equations (8)–(11) define several operators, but their explicit matrix construction is not shown. Providing a concrete example using a small toy data would improve reproducibility and make the work accessible to readers unfamiliar with TDL.\nAnswer: We thank the reviewer for pointing out this. In section 1.4 of the revised supplementary information, we have added pseudocode for computing the discrete differential operators and performing the Hodge decomposition. We have also included a simple example in section 1.5 as follows:\nConsider the manifold in Figure S2, it is a 2D disk in the Cartesian grid. The resulting\nFigure S2: A 2d disk and its discrete manifold representation under normal and tangential boundary condi- tions.\ndiscrete manifolds under the normal and tangential boundary conditions are displayed on the left and right sides, respectively. Specifically, the discrete manifold under the normal boundary condition consists of 3 vertices, 8 edges, and 4 faces, whereas the discrete manifold under the tangential boundary condition consists of 4 vertices, 4 edges, and 1 face.\nWe assign an ordering to the primal cells as follows:\n• vertices: 1, 2, 3, 4, 5, 6, 7, 8, 9\n• edges: 12, 23, 45, 56, 78, 89, 14, 47, 25, 58, 36, 69\n• faces: 1245, 2356, 4578, 5689\nWe use the construction of the 0-th Laplacian matrix under the tangential boundary condition L as an example to illustrate the computation. For this computation, we need\n0,t\nto first compute the discrete differential D on the entire grid. Then compute the projection\n0\nmatrix P and P for the tangential boundary condition. Using these projection matrices,\n1,t 0,t\nthe discrete differential under the tangential boundary condition D can be derived by D =\n0,t 0,t\nP D PT . After this, we need to compute the discrete Hodge star S′ on the entire grid.\n1,t 0 0,t 1,t\nThen the discrete Hodge star on the manifold can be derived by S = P S′ PT . Having\n1,t 1,t 1,t 1,t\nthese, the Laplacian can be computed by L = DT S D . We give the details as follows:\n0,t 0,t 1,t 0,t\n1. Let ∂ be the boundary matrix from edges to vertices on the entire grid, then D = ∂T\n1 0 1\n \n−1 0 0 0 0 0 −1 0 0 0 0 0\n 1 −1 0 0 0 0 0 0 −1 0 0 0 \n \n 0 1 0 0 0 0 0 0 0 0 −1 0 \n \n \n0 0 −1 0 0 0 1 −1 0 0 0 0\n \n \n∂ =  0 0 1 −1 0 0 0 0 1 −1 0 0 , D = ∂T\n1   0 1\n 0 0 0 1 0 0 0 0 0 0 1 −1\n \n 0 0 0 0 −1 0 0 1 0 0 0 0 \n \n 0 0 0 0 1 −1 0 0 0 1 0 0 \n \n0 0 0 0 0 1 0 0 0 0 0 1\n2. The projection matrices are obtained by selecting the rows of the identity matrix cor- responding to the cells included in the discrete manifold. Thus, we have\n   \n0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0\n0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0\nP =  , P =  \n0,t 0 0 0 0 1 0 0 0 0 1,t 0 0 0 0 0 0 0 0 1 0 0 0\n   \n0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0\n3. The discrete differential D can be computed as follows\n0,t\n \n−1 1 0 0\n 0 0 −1 0\nD = P D PT =  \n0,t 1,t 0 0,t  0 0 0 0\n \n0 −1 0 0\n4. The discrete Hodge star under tangential boundary condition on the entire grid S′ is\n1,t\na diagonal matrix in which only the entries associated with edges 23, 56, 25, and 36\nare nonzero and all remaining diagonal entries are zero. Each nonzero diagonal entry\nrepresents the ratio between the volume of the dual cell and that of the corresponding\nprimal cell. Under the tangential boundary condition, we follow the strategy described\nin paper “Hodge Decomposition of Vector Fields in Cartesian Grids”, where the primal\ncell volumes are kept unchanged while the dual cell volumes are adjusted. That is, these\nnonzero entries are given by the volume of the intersection between each dual cell and\nthe manifold. Assuming each grid edge has unit length, the intersection of the disk and\nthe dual edge of 23 is a line segment of length 0.24, so is 36. The intersection of the\ndisk and the dual edge of 56 is a line segment of length 0.64, so is 25. Hence, we have\n \n0 0 0 0 0 0 0 0 0 0 0 0\n0 0.24 0 0 0 0 0 0 0 0 0 0\n \n0 0 0 0 0 0 0 0 0 0 0 0\n \n0 0 0 0.64 0 0 0 0 0 0 0 0\n \n \n0 0 0 0 0 0 0 0 0 0 0 0\n \n \n0 0 0 0 0 0 0 0 0 0 0 0\nS′ =   (1)\n1,t 0 0 0 0 0 0 0 0 0 0 0 0\n \n0 0 0 0 0 0 0 0 0 0 0 0\n \n0 0 0 0 0 0 0 0 0 0.64 0 0 0\n \n \n0 0 0 0 0 0 0 0 0 0 0 0\n \n \n0 0 0 0 0 0 0 0 0 0 0.24 0\n0 0 0 0 0 0 0 0 0 0 0 0\n5. The discrete Hodge star under the tangential boundary condition on the manifold is as\nfollows:\n \n0.24 0 0 0\n 0 0.64 0 0 \nS = P S′ PT =   (2)\n1,t 1,t 1,t 1,t  0 0 0.64 0 \n \n0 0 0 0.24\n6. The Hodge Laplacian can be computed as\n \n0.24 −0.24 0 0\n−0.24 0.48 0 0\nL = DT S D =   (3)\n0,t 0,t 1,t 0,t  0 0 0.64 0\n \n0 0 0 0\nThe eigenvalues of L are 0, 0.09, 0.63, and 0.64. The single zero eigenvalue actually\n0,t\ncorresponds to the Betti number β = 0 of the manifold. If we replace the Hodge star\n0\nmatrix with the identity matrix, we get the BIG Laplacian LB . The BIG Laplaican\n0,t\nLB also has exactly one zero eigenvalue since it preserves the topology of the manifold.\n0,t\nAnswers to Reviewer 2 and Reviewer 3’ Comments\nComment : 1)Summary: ThemanuscriptintroducesManifoldTopologicalDeepLearn- ing (MTDL), which represents images as discrete manifolds, constructs a vector field on that manifold, applies a discrete Hodge decomposition into curl-free / divergence-free / harmonic components, concatenates those channels, and feeds the result into a CNN. The model is evaluated on MedMNIST v2 (11 2D + 6 3D datasets) and reported to outperform prior methods.\nThe idea of converting images to vector fields and apply a discrete Hodge decomposition before feeding the components to a CNN model for classification is promising. However, the paper needs further clarification on methodological details, additional ablation studies, and rigorous evaluation. Detailed comments below outline specific improvements needed before publication.\nOur view is that the manuscript needs a Major Revision:\nAnswer: We appreciate the summary from the reviewers of our work and thank the\nreviewers for the positive comments.\nMethodology:\nComment : 2) Explain how the harmonic component h is computed.\nAnswer: We thank the reviewers for pointing out this. For the decomposition\nω = dα +δβ +h,\nn t\nwhere ω ∈ Ω1,α ∈ Ω0,β ∈ Ω2,h ∈ kerd∩kerδ, we first determine the potentials α and\nn n t t n\nβ , and then compute h as\nt\nh = ω −dα −δβ .\nn t\nNote that α and β are not unique, as α +dη and β +δγ with any η ∈ Ω0 and γ ∈ Ω2 can\nn t n t n t\nserve as valid potentials for the first two terms. We solve this issue by imposing the gauge\nconditions. Specifically, we restrict\nα ∈ kerδ ∩Ω0, β ∈ kerd∩Ω2\nn n t t\nUnder these conditions, we have\nδω = δdα = (δd+dδ)α = ∆α, dω = dδβ = (dδ +δd)β = ∆β (4)\nt t t\nThus, α and β can be uniquely determined from the linear system (4) by resolving the rank\nn t\ndeficiency of ∆.\nComment : 3) Provide Pseudo-code for the proposed methodological pipeline.\nAnswer: We thank the reviewers for the constructive comments. In the discretization,\nthree key things need to be considered: the discrete manifold, the discrete differential, and\nthe discrete Hodge star. Once these three components are properly defined, all remaining\ndiscrete operators, such as the Laplacian and the Hodge decomposition, can be assembled or\nimplemented following the same definitions of their continuous counterparts.\nWe provide Algorithm 1 for constructing the discrete manifold under normal and tan-\ngential boundary conditions, Algorithm 2 for computing the projection matrices associated\nwith each boundary condition, Algorithm 3 for building the k-th boundary matrix ∂ on the\nk\nentire Cartesian grid, and Algorithm 4 for computing the discrete Hodge star matrix S′\nk,n\nand S′ . Using these four algorithms, we obtain the k-th boundary matrix ∂ , the discrete\nk,t k\nHodge star S′ ,S′ , and the projection matrices P and P corresponding to the normal\nk,n k,t k,n k,t\nand tangential boundary conditions, respectively.\nBased on these matrices, we can derive the k-th discrete differential operators D and\nk,n\nD , the discrete codifferentials δ and δ , and the discrete Hodge Laplacians L and\nk,t k,n k,t k,n\nL for the two types of boundary conditions. Finally, using these discrete operators, we can\nk,t\ncompute the Hodge decomposition, which is presented in Algorithm 5.\nWe have added these algorithms to section 1.4 of the revised supplementary information.\nComment : 4) Clarify how the manifold M and its boundary ∂M are defined for each\ndataset. In practical settings, computing normal and tangential boundary conditions requires\na segmentation mask or at least a binary indicator of the region of interest. In Section 4.3.1, the segmentation step is described only as thresholding but the choice of threshold (fixed, adaptive, or dataset-specific), preprocessing (e.g., normalization or filtering), and handling of ambiguous or noisy background regions can significantly influence the resulting manifold boundary ∂M, which in turn affects the computed boundary conditions and the subsequent Hodge decomposition. I recommend the authors (1) explicitly describe the segmentation pro- cedure per dataset, (2) discuss its influence on the topological representation, and (3) provide\nAlgorithm 1: Computing Discrete Manifold\nInput: A manifold M on the Cartesian grid I.\nOutput: M : Discrete manifold under the normal boundary condition; M : discrete\nn t\nmanifold under the tangential boundary condition.\nSet M = ∅, M = ∅;\nn t\nfor cell σ in the grid I do\nif one vertex of σ is in M then\nadd σ into M\nn\nLet ⋆σ be the dual cell of σ;\nif one vertex of ⋆σ is in M then\nadd σ into M\nt\nreturn (M ,M )\nn t\nAlgorithm 2: Computing Projection Matrix\nInput: A discrete manifold M with normal or tangential boundary condition on\nb\nthe Cartesian grid I.\nOutput: P : the k-th projection matrix\nk,b\nLet n be the number of k-cells of I;\nk\nSet P be the n×n identity matrix;\nk,b\nfor k-cell σ in the grid I do\nk\nif σ is not in M then\nk b\nRemove the row corresponding to σ in P\nk k,b\nreturn P\nk,b\nan ablation or sensitivity analysis to assess robustness to segmentation variations in the re- sults.\nAnswer: We thank the reviewers for the insightful suggestions.\nFor the segmentation strategy to determine the manifold region from medical images,\nwe carefully examined all 17 datasets. As shown in Figure S3, the majority of the datasets,\nwith the exceptions of PathMNIST, DermaMNIST, and BloodMNIST, contain images with a\nblack background. Note that all datasets share a consistent pixel range from 0 to 255, where\nthe black background corresponds to a pixel value of 0. Based on this observation, we adopt\na unified and straightforward segmentation rule: pixels with intensity strictly greater than\ncutoff of 0 are treated as belonging to the manifold. This ensures consistent preprocessing\nacross all datasets.\nFor the impact of segmentation parameter to the manifold representation, a very small\ncutoff can inadvertently include noise, while a large cutoff may discard meaningful informa-\ntion. Using a cutoff of 0 can guarantee that no potentially informative signal is removed.\nIn addition, our vector-field representation is obtained using the gradient operator. There-\nfore, regions containing low-amplitude noise or approximately constant nonzero background\nAlgorithm 3: Computing Boundary Matrix\nInput: A Cartesian grid I.\nOutput: ∂ : the k-th boundary matrix from k-cells to (k −1)-cells\nk\nLet n be the number of k-cells of I;\nk\nSet ∂ be the n ×n zero matrix;\nk k−1 k\nfor i = 1 to n do\nk\nLet σ be the i-th k-cell of I;\nfor j = 1 to n do\nk−1\nLet τ be the j-th (k −1)-cell of I;\nif τ is a face of σ then\nCompute the incidence sign s(τ,σ) ∈ {−1,+1} according to the\norientation;\nSet ∂ (j,i) = s(τ,σ);\nk\nreturn ∂\nk\nvalues produce near-zero gradients, contributing negligible vector-field information. As a re- sult, selecting 0 as the segmentation threshold is both empirically effective and theoretically reasonable.\nFor the impact of segmentation parameter to model performance, we add an ablation\nstudy by varying the cutoff values across 0, 5, 10, 15, 20, 25, 30. Experiments were per-\nformed on four representative datasets: RetinaMNIST, PneumoniaMNIST, DermaMNIST,\nand OrganSMNIST. The results are summarized in Table S5. It can be seen that model\nperformance remains stable across this range of thresholds. This indicates that the proposed\nmethod is robust to the choice of segmentation cutoff, and that the default setting of 0 is\nreliable.\nWe have added these discussions to section 6.3 of the revised supplementary information.\nComment : 5) The paper would strongly benefit from an explicit pseudocode summariz- ing the full implementation workflow, including how boundary conditions, discrete operators, and vector-field decompositions are computed. Although the Supplementary Material gives substantial mathematical detail, a procedural description (with solver specifications and preprocessing steps) is essential for reproducibility.\nAnswer: We thank the reviewers for the constructive suggestions. We have added detailed\npseudocode to section 1.4 of the revised supplementary information for computing the\ndiscrete manifold with normal and tangential boundary conditions, constructing the discrete\ndifferential operators, and performing the Hodge decomposition. For convenience, we also\nprovide the pseudocode here as Algorithm 1, Algorithm 2, Algorithm 3, Algorithm 4, and\nAlgorithm 5.\nResults:\nComment : 6) Vector-field construction is under-specified in the main text. The pa- per mentions several methods and refers to the Supplementary Information, but the main\nAlgorithm 4: Computing Discrete Hodge Star\nInput: A manifold M on a Cartesian grid I.\nOutput: S′ : the k-th discrete Hodge star matrix under the normal boundary\nk,n\ncondition on the grid I; S′ : the k-th discrete Hodge star matrix under\nk,t\nthe tangential boundary condition on the grid I.\nLet n be the number of k-cells of I;\nk\nSet S′ be the n ×n identity matrix;\nk,n k k\nSet S′ be the n ×n identity matrix;\nk,t k k\nfor i = 1 to n do\nk\nLet σ be the i-th k-cell of I;\nLet ⋆σ be the dual cell of σ;\nSet S′ (i,i) = volume(⋆σ) ;\nk,n volume(σ∩M)\nSet S′ (i,i) = volume(⋆σ)∩M;\nk,t volume(σ)\nreturn (S′ ,S′ )\nk,n k,t\ntext does not make clear which construction is used for each experiment, what hyperparam- eters were used, and how sensitive results are to that choice. Since the vector field is the core representation, an ablation is necessary: different vector-field constructions (flow-based, gradient-based, others) and their effect on final performance.\nAnswer: We thank the reviewers for the insightful comments. The equation (15) in the\n“Decomposed Image Generation” section of the main text is the formula we used to construct\nthe vector field in our MTDL model for all 17 datasets. This formulation corresponds to a\nspecial case of the gradient-based method presented in the supplementary information with\nparameters s = t = 1.\nTo further assess the impact of different vector field construction strategies, we have\nadded an ablation study comparing three distinct methods: the gradient-based method, the\nflow-based method, and the RGB-based method, evaluated on the RetinaMNIST, Pneumo-\nniaMNIST, DermaMNIST, and OrganSMNIST datasets. For the gradient-based method,\nwe examine three parameter settings: (s,t) = (1,1), (1,2), and (2,1). For the RGB-based\nmethod, we construct vector fields using channel pairs (R,G), (R,B), and (G,B). The results\nof all methods are shown in Table S6. As shown in the table, the gradient-based method with\n(s,t) = (1,1) achieves the best performance on the RetinaMNIST, PneumoniaMNIST, and\nDermaMNIST datasets, whereas the RGB-based method yields the best performance on the\nOrganSMNIST dataset. These results show that the optimal vector field construction strategy\ncan vary with dataset characteristics. This could enhance the performance of our method and\nbe further explored in the future study.\nWe have added these to section 2.4 of the revised supplementary information.\nComment : 7) Discretization errors and orthogonality checks: Provide results showing\nthe discrete decomposition is numerically orthogonal. Report:\nAlgorithm 5: Computing Hodge Decomposition\nInput: Vector field ω\nOutput: The decomposition (dα ,δβ ,h) where ω = dα +δβ +h\nn t n t\nUse Algorithm 1 to get the discrete manifold (M ,M );\nn t\nUse Algorithm 2 and (M ,M ) to compute the projection matrix P ,P\nn t k,n k,t\nUse Algorithm 3 to compute the boundary matrix ∂ ;\nk\nUse Algorithm 4 to compute the Hodge star S′ ,S′ ;\nk,n k,t\nSet D = ∂T ;\nk k+1\nSet the discrete differential D = P D PT , D = P D PT ;\n0,n 1,n 0 0,n 1,t 2,t 1 1,t\nSet the discrete Hodge star S = P S′ PT ,S = P S′ PT ;\nk,n k,n k,n k,n k,t k,t k,t k,t\nSet the discrete Laplacian L = DT S D , L = S D S−1DT S ;\n0,n 0,n 1,n 0,n 2,t 2,t 1,t 1,t 1,t 2,t\nSet δ = S−1DT S ;\n1,n 0,n 0,n 1,n\nSolve L α = δ ω to get α ;\n0,n n 1,n n\nSolve L β = D ω to get β ;\n2,t t 1,t t\nSet δ = S−1DT S ;\n2,t 1,t 1,t 2,t\nSet h = ω −D α −δ β ;\n0,n n 2,t t\nreturn (D α ,δ β ,h)\n0,n n 2,t t\n||v −(d ϕ+δψ +h)|| | < d ϕ,δψ > |\no 2 0\nerr = ,err =\nrecon ortho\n||v|| ||d ϕ|| ||δψ||\n2 0 2 2\nProvide these for representative images or in a table for a dataset with mean ± std.\nAnswer: We thank the reviewers for the constructive suggestions.\nFor reconstruction error, in our computation, we first compute the curl-free component\ndϕ and the divergence-free component δψ, and then the harmonic component is obtained by\nh = v −dϕ−δψ\nAs a result, the reconstruction error is always 0.\nFor orthogonality error, we computed the normalized pairwise inner product between the\ncurl-free, divergence-free, and harmonic parts on five representative datasets: RetinaMNIST\n(1,600 samples), PneumoniaMNIST (5856 samples), DermaMNIST (10,015 samples), Or-\nganAMNIST (58,830 samples), and PathMNIST (107,180 samples). The normalized inner\nproduct for two components α,β is defined as\n| < α,β > |\n(α,β) = .\n||α|| ||β||\n2 2\nWe report the mean and standard deviation of these values across all samples for each dataset,\nwhich is shown in Table S7. The results indicate that the three components exhibit strong\nnumerical orthogonality across all evaluated datasets.\nWe have added these discussions to section 1.3 of the revised supplementary information.\nComment : 8) Comparisons / fairness of baselines:\nAnswer: We thank the reviewers for the comments. Our manuscript was originally\nsubmitted earlier this year, so the recently published related work has not been included.\nTo provide an updated and comprehensive comparison, we conducted a thorough literature\nsearch for works published in 2025 that also utilized this database. We summarize these newly\nidentified results in Table S8 to present an up-to-date benchmark comparison. It can be seen\nthat our model achieves the best overall performance among all these methods. Notably,\nmany of the recently published models are evaluated on only a subset rather than the full set\nof datasets while our model is evaluated on all the datasets. These new results have been\nadded into section 9 of the revised supplementary information.\n[1] Halder, A., Dalal, A., Gharami, S., Wozniak, M., Ijaz, M. F., & Singh, P. K.\nFigure S3: Illustration of the 11 2D datasets and 6 3D datasets in the MedMNIST database.\n(2025). A fuzzy rank-based deep ensemble methodology for multi-class skin cancer classifica- tion. Scientific Reports, 15(1), 6268.\n[2] Wu, B., Huang, J., & Duan, Q. (2025). Real-time intelligent healthcare enabled by\nfederated digital twins with aoi optimization. IEEE Network.\n[3] Wang, G., Zhu, Q., Song, C., Wei, B., & Li, S. (2025). MedKAFormer: When\nKolmogorov-Arnold Theorem Meets Vision Transformer for Medical Image Representation.\nIEEE Journal of Biomedical and Health Informatics.\n[4] Bai, Y., Bai, L., Yang, X., & Liang, J. (2025). Label-Semantic-Based Prompt\nTuning for Vision Transformer Adaptation in Medical Image Analysis. IEEE Transactions\non Circuits and Systems for Video Technology.\n[5] Singh, A., Eising, C., Denny, P., & van de Ven, P. (2025). Improving GNNs\nfor Image Classification: Addressing Homophily Challenges. IEEE Open Journal of the\nComputer Society.\n[6] Zechen, Z., Xuelei, H., Fengjun, Z., & Xiaowei, H. (2025). PSNAS-Net: Hybrid\ngradient-physical optimizationfor efficient neural architecture search in customized medical\nimaging analysis. Expert Systems with Applications, 288, 128155.\n[7] Du, Y., Zhang, J., Zeevi, T., Dvornek, N. C., & Onofrey, J. A. (2025, April).\nSre-conv: Symmetric rotation equivariant convolution for biomedical image classification. In\n[8]Singh, R., &Guggilam, S.(2025). IterativeMisclassificationErrorTraining(IMET):\nAn Optimized Neural Network Training Technique for Image Classification. IEEE Access.\n[9] Rahman, M., & Zhuang, J. (2025, September). NQNN: Noise-Aware Quantum\nNeural Networks for Medical Image Classification. In International Conference on Medical\nImage Computing and Computer-Assisted Intervention (pp. 433-442). Cham: Springer\nNature Switzerland.\nComment : 9) Provide Pseudo-code for the proposed methodological pipeline.\nAnswer: We thank the reviewers for the constructive suggestions. We have added the\npseudocode in the revised supplementary information.\nComment : 10) The paper reports substantial improvements across multiple MedM-\nNIST tasks; however, it is unclear whether the baseline models were trained and tuned under\nidentical experimental conditions (e.g., data augmentation, optimizer settings, and hyperpa-\nrameter search strategies). Please clarify whether the baseline results were all adopted from\nthe original publications or reproduced by the authors. If any were re-implemented, details\non the hyperparameter parity and training cost should be provided. Without this, it’s hard to\nknow whether gains derive from topology or from training differences.\nAnswer: We thank the reviewers for pointing this out. All baseline results reported\nin our manuscript are directly taken from their original publications. The MedMNIST v2\ndatabase provides clearly defined training, validation, and testing splits, and we have strictly\nfollowed these official splits to ensure a fair comparison with existing baseline models. We\nhave emphasized this point in section 2.2.2 of the revised manuscript.\nComment : 11) Statistical analysis: While the reported results on MedMNIST v2 are\npromising, the paper lacks statistical analysis. Since many datasets are small and stochas-\ntic training introduces non-negligible variance, it is important to report mean +− standard\ndeviation across multiple runs and perform statistical significance testing. A paired t-test\nor Wilcoxon test within each dataset, and a Friedman test with Nemenyi post-hoc across\ndatasets, would clarify whether the observed improvements are statistically significant rather\nthan due to random variation. Reporting confidence intervals and effect sizes (e.g., Cohen’s\nd or Cliff’s δ) would further strengthen the empirical evidence.\nAnswer: We thank the reviewers for the valuable suggestions.\nFor mean and standard deviation, all reported results of our MTDL model are the aver-\naged values across three independent runs with different random seeds. We have added the\ncorresponding std value in section 9 of the revised supplementary information.\nFor the paired t-test and Cohen’s d analysis, the existing baseline models only provide\nsingle performance values without reporting variance across multiple runs. Therefore, we\ncannot directly perform paired t-tests or compute Cohen’s d between our model and these\nbaselines, because both analyses require multiple runs for each method to estimate variance.\nFor the Friedman test, we use the scipy.stats.friedmanchisquare package to compute the\np-values on 11 2D datasets and 6 3D datasets, respectively. For the 2D analysis, we only\nTable S5: Performance comparison (AUC and ACC) across different segmentation cutoff parameters on\nrepresentative datasets.\nCutoff Retina Pneumonia Derma OrganS\nAUC ACC AUC ACC AUC ACC AUC ACC\n0 0.874 0.655 0.986 0.910 0.962 0.836 0.978 0.809\n5 0.872 0.667 0.986 0.922 0.962 0.843 0.978 0.809\n10 0.870 0.655 0.985 0.915 0.964 0.832 0.978 0.802\n15 0.878 0.643 0.983 0.914 0.962 0.842 0.977 0.806\n20 0.874 0.673 0.985 0.917 0.963 0.838 0.978 0.808\n25 0.875 0.658 0.986 0.914 0.961 0.831 0.977 0.802\n30 0.870 0.643 0.985 0.910 0.965 0.839 0.978 0.803\nTable S6: Performance (ACC & AUC) Comparison of different vector field generation methods, including\nthe gradient-based method with parameter (s,t), the flow-based method, and RGB-based method.\nMethod Retina Pneumonia Derma OrganS\nAUC ACC AUC ACC AUC ACC AUC ACC\nGradient(1,1) 0.874 0.655 0.986 0.910 0.962 0.836 0.978 0.809\nGradient(1,2) 0.864 0.643 0.982 0.904 0.958 0.815 0.979 0.815\nGradient(2,1) 0.854 0.613 0.983 0.899 0.953 0.821 0.979 0.818\nFlow 0.815 0.563 0.973 0.901 0.931 0.785 0.976 0.792\nRGB 0.803 0.585 0.976 0.867 0.952 0.816 0.981 0.818\ninclude baseline models that report results on all 11 MedMNIST2D datasets. Similarly, for\nthe 3D analysis, we only include baselines that report results on all 6 MedMNIST3D datasets.\nThe resulting p-values are 1.38e-14 (2D) and 5.32e-10 (3D), indicating significant perfor-\nmance differences among the compared models. Therefore, we provide Critical Difference\n(CD) diagrams for both the 2D and 3D cases. Figure S4 shows the CD diagram compar-\ning 12 models across the 11 MedMNIST2D datasets using the Nemenyi post-hoc test at the\nsignificance level α = 0.05. Our MTDL model achieves the best overall performance with\nthe lowest average rank among all compared methods. The CD = 3.55 shown at the bottom\nindicates the minimum rank difference required for two models to be considered significantly\ndifferent under the Nemenyi test. Under this threshold, MTDL is significantly better than\nmost compared methods. Similarly, Figure S5 shows the results for the 6 3D datasets, where\nour MTDL model also achieves the best overall performance with the lowest average rank\nand is significantly better than most other methods.\nWe have added these analysis into section 6.1 of the revised supplementary information.\nComment : 12) In Figure 3, the range of the AUC or ACC values on the x-axis\nshould be consistent across all plots to facilitate clearer and quicker visual comparison between\nresults.\nAnswer: We thank the reviewers for pointing this out. We have revised the figure\naccordingly.\nComment : 13) Provide preprocessing time per sample (vector field + decomposition),\nGPU memory/CPU usage, and training time per epoch for representative datasets (2D and\n3D).\nTable S7: Numerical analysis of the orthogonality of the three components in the Hodge decomposition\non the RetinaMNIST, PneumoniaMNIST, DermaMNIST, OrganAMNIST, and PathMNIST datasets. The\naverage and standard deviation values of the normalized inner product of all samples within each dataset\nare reported.\nMethod Retina Pneumonia Derma OrganA Path\nmean std mean std mean std mean std mean std\n(dϕ,δψ) 7e-20 4e-17 2e-19 2e-17 4e-20 3e-17 8e-20 3e-17 7e-20 2e-17\n(dϕ,h) 2e-13 9e-14 3e-14 2e-14 6e-14 1e-14 3e-14 3e-14 2e-14 1e-14\n(δψ,h) 2e-17 5e-17 2e-18 8e-17 7e-18 2e-16 5e-18 8e-17 2e-17 6e-17\nAnswer: We thank the reviewers for the suggestions. In our experiments, all 2D\nimages are of size n × 224 × 224, where n = 3 for RGB images and n = 1 for grayscale\nimages. All 3D images are grayscale volumes with a resolution of 64×64×64. The vector\nfield generation and Hodge decomposition procedures are computationally efficient. When\nimplemented on a single CPU core, processing 50 grayscale 2D images requires approximately\n2 minutes, while processing 50 RGB 2D images requires around 4 minutes. For 3D images,\nthe computation for 50 samples can be completed within 30 minutes. For model training,\nour MTDL framework is implemented in PyTorch and executed on an NVIDIA Tesla V100S\nGPU with 32 GB of memory. Representative training times for selected 2D and 3D datasets\nare reported in Table S9. As shown, the overall training process is also computationally\nefficient.\nThese runtime analysis have been added into section 7 of the revised supplementary\ninformation.\nDiscussions:\nComment : 14) Invariance/equivariance: discuss which geometric transforms of the\ninput (translation, rotation, scaling) preserve each decomposed channel and whether the CNN\narchitecture exploits or breaks these invariances.\nAnswer: We thank the reviewers for the constructive suggestions. The vector field\ngeneration method used in MTDL is translation-equivariant. However, it is not rotationally\nequivariant because the underlying Cartesian grid is not rotationally symmetric, nor is it\nscale-equivariant, since the grid spacing does not vary under geometric scaling. The Hodge\ndecomposition exhibits the same property: it is translation-equivariant, but not rotationally\nTable S8: Performance comparison between our model and recent models on the MedMNIST2D dataset.\nThe benchmark results are taken from the original papers. The best results are in bold. “-” indicates that\nthe corresponding metric was not reported.\nMethods Path Chest Derma OCT\nAUC ACC AUC ACC AUC ACC AUC ACC\nFRDEM [1] - - - - - 0.783 - -\nAIGC [2] - - - 0.83 - -\nMedKAFormer [3] 0.985 0.882 - - 0.919 0.748 0.974 0.810\nLPT [4] - 0.941 - - - 0.787 - 0.832\nSGDGCN [5] 0.987 0.907 0.928 0.854 0.919 0.771 0.933 0.783\nPSNAS-Net-ave [6] - - - - - 0.772 - -\nSRE-CONV [7] - 0.822 - 0.948 - 0.757 - 0.768\nIMET [8] - - - - - - 0.960 0.803\nNQNN [9] - 0.752 - - - 0.827 - 0.591\nMTDL 0.996 0.920 0.793 0.948 0.962 0.836 0.989 0.888\nMethods Pneumonia Retina Blood Tissue\nAUC ACC AUC ACC AUC ACC AUC ACC\nMedKAFormer [3] 0.961 0.904 0.729 0.535 0.998 0.965 - -\nLPT [4] - - - 0.575 - 0.967 - -\nSGDGCN [5] 0.963 0.918 0.782 0.550 0.994 0.940 0.801 0.613\nPSNAS-Net-ave [6] - 0.899 - - - 0.956 - -\nSRE-CONV [7] - 0.906 - 0.523 - 0.960 - 0.712\nIMET [8] 0.895 0.902 - - - - - -\nNQNN [9] - - - - - 0.737 - -\nMTDL 0.986 0.910 0.874 0.655 0.999 0.988 0.945 0.721\nMethods OrganA OrganC OrganS\nAUC ACC AUC ACC AUC ACC\nAIGC [2] - 0.890 - - - -\nMedKAFormer [3] 0.995 0.922 0.991 0.899 0.986 0.788\nLPT [4] - - - 0.906 - 0.796\nSGDGCN [5] 0.977 0.850 0.989 0.880 0.969 0.784\nPSNAS-Net-ave [6] - - - 0.913 - -\nSRE-CONV [7] 0.770 0.769 - - - -\nNQNN [9] - 0.814 - 0.799 - -\nMTDL 0.998 0.956 0.996 0.928 0.978 0.809\nFigure S4: The Critical Difference (CD) diagram over 11 MedMNIST2D datasets.\nor scale-equivariant.\nFortheCNNarchitectureemployedinourmodel, theconvolutionallayersaretranslation-\nequivariant, and the use of global average pooling at the final stage guarantees translation-\ninvariance at the classifier level. Therefore, the overall MTDL model achieves translation-\ninvariance. This analysis suggests that our framework could be further extended to incorpo-\nrate additional forms of equivariance or invariance, such as rotational or scale equivariance.\nComment : 15) Physical/semantic interpretation: Can we link the decomposed com-\nponents to expected image structure per modality (e.g., curl-free edge/gradient information;\ndivergence-free texture swirl)? so readers can interpret classifier reliance on components.\nAnswer: We thank the reviewers for the valuable comments. In our model, each image\nis first represented as a vector field, which is subsequently decomposed into a curl-free com-\nponent, a divergence-free component, and a harmonic component. These components have\nwell-established geometric and physical interpretations. Specifically, a curl-free vector field\nexhibits no local rotational behavior. A divergence-free vector field has neither sources nor\nsinks. Source regions correspond to positive divergence, where local intensity or structural\npatterns spread outward, while sink regions correspond to negative divergence, where patterns\nconverge inward. The harmonic component is both curl-free and divergence-free, capturing\nthe smoothest global properties present in the vector field.\nIn medical image analysis, the curl-free component can be used to capture the lesion\nboundaries and margins, as the directions of maximal intensity change typically align with\nanatomical edges. For example, in mammography, malignant lesions often show sharper and\nmore irregular boundaries compared with the benign lesions. Therefore, the curl-free compo-\nnent can be used to distinguish the malignant from benign cases. The divergence-free compo-\nnent captures the vorticity and local rotating behavior in the vector field. Consequently, it can\nbe used to capture the irregular textures, such as the swirling patterns and the heterogeneous\ninternal textures due to necrosis and fibrosis. Especially, this will be useful for malignant tu-\nmors since the swirling textures are common in malignant tumors such as breast cancer, brain\nFigure S5: The Critical Difference (CD) diagram over 6 MedMNIST3D datasets.\ntumor, and hepatocellular carcinoma. The harmonic component reflects the global topological\nproperties of the image and provides robust descriptors insensitive to noises and small pertur-\nbations. For a more comprehensive understanding of how the decomposed components relate\nto visual patterns such as edges, textures, and global structures, additional interpretability\nanalysis would be required. We regard this as an important direction for future work.\nComment : 16) How the method might generalize to real-world biomedical images,\ne.g., high-resolution CT/MRI, anisotropic voxels, or dynamic sequences? MedMNIST is\nsmall and standardized; clinical data might pose challenges.\nAnswer: We thank the reviewers for pointing out this. For high-resolution real-world\nbiomedical images, we have an ablation study in Section 2.2.5 of the main text. Specifically,\nwe evaluate our model on the clinical HAM10000 dataset, which contains 10,015 dermato-\nscopic images of uniform size 3×600×450. We center-crop each image to 3×450×450\nand then resize it to five different resolutions: 3 × 35 × 35, 3 × 75 × 75, 3 × 150 × 150,\n3×300×300, and 3×450×450. Our MTDL model is then evaluated on each resolution.\nHere we list the results in Table S10. It can be seen that the model performance improves\nconsistently as the input resolution increases, indicating that MTDL is capable of extracting\nricher and more detailed features from higher-resolution inputs. Notably, even at the lowest\nresolution 3×35×35, MTDL achieves an AUC (ACC) of 0.943 (0.797), demonstrating that\nthe model maintains strong and stable performance across a wide range of image resolutions.\nTable S10: Performance of MTDL over different image reslutions on the clinical dataset HAM10000.\nResolution 3×35×35 3×75×75 3×150×150 3×300×300 3×450×450\nAUC 0.943 0.947 0.958 0.970 0.973\nACC 0.797 0.803 0.822 0.856 0.863\nFor dynamic sequences, Hodge decomposition can be directly applied to analyze the tem-\nporal evolution of the three components and thereby characterize dynamic patterns. Alterna-\ntively, representative slices from the sequence can be selected, and Hodge decomposition can\nbe performed independently on each slice. The resulting components can then be concate-\nnated to form a multi-time, multi-channel image representation, which can be combined with\na CNN to perform prediction tasks.\nComment : 17) A discussion (or ideally an experiment) comparing MTDL to surface-based or mesh-based features could be helpful. This might help readers understand when grid-based manifold representation is preferred vs direct mesh representation or potential ways these approaches could complement each other in biomedical imaging applications.\nAnswer: We thank the reviewers for the thoughtful comments. The Hodge theory em-\nployed in this work is formulated on a Cartesian grid. This corresponds to the Eulerian\ndiscretization of the continuous Hodge theory. An alternative is the Lagrangian discretiza-\ntion, in which the Hodge operators are defined on triangular or tetrahedral meshes. The\nchoice between Eulerian and Lagrangian approaches is usually determined by the nature of\nthe underlying data. For example, if a well-defined surface segmentation is available, the\nsurface can be naturally represented by a triangular mesh, making the Lagrangian method\ndirectly applicable. In contrast, when the data are given on a regular grid, such as 2D or 3D\nbiomedical images, the Eulerian method can be applied without additional processing. It is\nimportant to note that Lagrangian methods typically require accurate surface extraction and\nmesh construction, which is not feasible when a clean and anatomically meaningful surface\nis not present or cannot be robustly segmented.\nFor biomedical image analysis, the Eulerian formulation is more suitable compared to the\nLagrangian formulation for two reasons: (1) The primary data modality in clinical practice\nconsists of images that include background regions, noise, and heterogeneous tissue bound-\naries. Applying a Lagrangian Hodge framework would require an additional segmentation\nstep to extract surface meshes, which introduces methodological complexity and potential\nsegmentation-induced variability. (2) Biomedical images are inherently defined on a Carte-\nsian grid, providing a natural domain for Eulerian discrete differential operators. This nat-\nural grid structure makes Eulerian Hodge decomposition both computationally efficient and\nwell-aligned with the data representation in medical image analysis.\nComment : 18) The paper could benefit from explicit discussion of scenarios where\nMTDL might fail.\nAnswer: We thank the reviewers for pointing out this. (1) MTDL relies on the con-\nstruction of image-derived vector fields and their Hodge decomposition. Therefore, when the\nimage lacks meaningful intensity gradients or is dominated by noise, the curl-free, divergence-\nfree, and harmonic components may not produce informative representations. For example,\nimages with extremely low contrast or strong speckle noise may lead to unstable decompo-\nsitions. (2) MTDL is based on Cartesian grids. Consequently, for tasks where the studied\nobject is inherently a high-resolution surface manifold, such as the brain surface, mesh-based\nor surface-based Hodge decomposition may be more suitable. In such cases, applying MTDL\ndirectly to surface data may fail to capture important intrinsic geometric features."}, {"round": 2, "reviewer_comments": {"reviewer_1": "I thank the authors for the detailed response and the extensive additional experiments. Overall, I’m satisfied with the response. I only have a few remaining concerns.\nRegarding Comment 2 on interpretability. I noticed that another reviewer raised a related concern in Comment 15. While I accept the abstract explanations, I was hoping to see some concrete evidence in the specific application. For example, the authors mention that “the directions of maximal intensity change typically align with anatomical edges”, “in mammography, malignant lesions often show sharper and more irregular boundaries compared with benign lesions”, and “swirling patterns and heterogeneous internal textures arise due to necrosis and fibrosis”. I do not expect all of these claims to be demonstrated, but even one or two illustrative examples would significantly improve the clarity and persuasiveness of the interpretability claims.\nFinally, after reading the rebuttal to another reviewer, I agree that the lack of rotation and scale equivariance appears to be an inherent limitation of the proposed method. It would be helpful if this limitation could be properly addressed, or at least stated more explicitly in the discussion.\nMoreover, a few other things to certainly correct for:\n- Introduction exemplifies that the 'DNA chains' as differentiable manifolds. I think this is incorrect and must be toned down - differentiable manifolds can only approximate these at best.\n- Is the original image fully recoverable by the decomposition, e.g. by gradient reconstruction etc.? Which information do we lose particularly? For which data is this not acceptable? Given that only some medical data is used, such a discussion would be helpful.", "reviewer_2": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_3": "Thank you to the authors for the thoughtful and thorough review. I am pleased that all of my previous concerns have been fully addressed.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "Response to Reviewers\nAnswers to Reviewer 1’ Comments\nI thank the authors for the detailed response and the extensive additional experiments.\nOverall, I’m satisfied with the response. I only have a few remaining concerns.\nAnswer: We thank the reviewer for the positive comments.\nComment : 1) Regarding Comment 2 on interpretability. I noticed that another re-\nviewer raised a related concern in Comment 15. While I accept the abstract explanations, I\nwas hoping to see some concrete evidence in the specific application. For example, the au-\nthors mention that “the directions of maximal intensity change typically align with anatom-\nical edges”, “in mammography, malignant lesions often show sharper and more irregular\nboundaries compared with the benign lesions”, and “swirling patterns and heterogeneous internal\ntextures arise due to necrosis and fibrosis”. I do not expect all of these claims to be demon-\nstrated, but even one or two illustrative examples would significantly improve the clarity and\npersuasiveness of the interpretability claims.\nAnswer: We thank the reviewer for the comments. Mammography refers to X-ray\nimaging of the breast, whereas ultrasound is a different imaging modality. Here we use the\nultrasound image to illustrate the general observation that malignant breast lesions often\nexhibit more irregular boundaries compared with the benign cases. As shown in Figure S1, the\nground-truth annotation clearly delineates the lesion boundary, where the malignant case\npresents a more irregular contour.\nComment : 2) Finally, after reading the rebuttal to another reviewer, I agree that the\nlack of rotation and scale equivariance appears to be an inherent limitation of the proposed\nmethod. It would be helpful if this limitation could be properly addressed, or at least stated\nmore explicitly in the discussion.\nAnswer: We thank the reviewer for the comments. Due to the page limit, a discussion\nof this issue has been given to Section 1.6 of the revised Supplementary Information.\nComment : 3) Moreover, a few other things to certainly correct for: Introduction\nexemplifies that the ’DNA chains’ as differentiable manifolds. I think this is incorrect and\nmust be toned down - differentiable manifolds can only approximate these at best.\nAnswer: We thank the reviewer for the comments. The term “DNA chains” has been\nremoved in the revised manuscript.\nobtained from Al-Dhabyani et al., Data in Brief (2020), 28:104863.\nComment : 4) Is the original image fully recoverable by the decomposition, e.g. by\ngradient reconstruction etc.? Which information do we lose particularly? For which data is\nthis not acceptable? Given that only some medical data is used, such a discussion would be\nhelpful.\nAnswer: We thank the reviewer for pointing out this. In our model, the original image\nI is first transformed into a discrete vector field V(i,j) = (x ,y ) via\ni,j i,j\nI(i+1,j)−I(i−1,j)\nx = ,\ni,j\n2\n(1)\nI(i,j +1)−I(i,j −1)\ny = .\ni,j\nsition is algebraically exact for the discrete vector field, meaning that V can be reconstructed\nfrom its components within numerical precision.\nGiven a prescribed value at one reference point I(i ,j ), the remaining values I(i,j)\n0 0\ncan be computed from the vector field V with formula (1). However, since the absolute value\n2\nat the reference point is unknown, the reconstructed image may differ from the original by\nan additive constant. Therefore, the original image can be recovered only up to a constant\noffset.\nIt is true that there are situations in which the computation of Hodge decomposition\nmay fail. In particular, when the imaged object contains extremely thin structures, accurate\nresolution of these features requires a substantially refined discretization grid. Although, in\nprinciple, a sufficiently fine grid allows for meaningful decomposition of a fixed-size object,\nin practice this significantly increases computational cost. Consequently, Hodge decomposition\nmay become computationally impractical for data dominated by very thin structures.\nAnswers to Reviewer 3’ Comments\nComment : Thank you to the authors for the thoughtful and thorough review. I am\npleased that all of my previous concerns have been fully addressed.\nAnswer: We are happy that the reviewer is satisfied with our revision."}]

Decision letter