1 Introduction

Automated medical image analysis with fast, accurate, and insightful results holds substantial significance for clinicians during practice. These analysis results have a wide range of clinical applications, possible usage including disease recognition (Hu et al. 2018), quantitative analysis (Hu et al. 2019), treatment planning (Tian et al. 2022), and prognostic prediction (Zhou et al. 2022), etc. Numerous computer-aided systems (CAD) have been proposed to assist clinicians in accomplishing medical image analysis tasks. Conventional methods used in these analysis tasks are mainly based on handcrafted low-level (SIFT or HOG) that have limited representation capability, resulting in unsatisfactory performance.

During the last decade, the thriving of deep neural networks (DNNs) has changed the paradigm of image analysis from low-level hand-crafted to high-level learning-based ones. The DNNs have achieved enormous success in plenty of applications (LeCun et al. 2015), including object recognition (Krizhevsky et al. 2012), natural language processing (NLP) (Bahdanau et al. 2014), speech recognition (Graves et al. 2013), etc. The essential of DNNs is the implementation of a highly non-linear function F through layer-wise transformations. The F accomplish the mapping \(x \rightarrow y\) with learnable parameter \(W^{l}\), i.e., \(y=F(x;W^l)\), where \(l \in \{1,2,...,L\}\). One of the most simple connectionisms is the fully connected (FC) layer, where each neuron between two layers is connected through connection weight. The computational principle of the FC layer can be formulated as \(a^{l+1}=f( W^{l}a^{l})\), with \(W^{l} \in {\mathbb {R}}^{n_{l+1} \times n_{l}}\). The symbol \(a^l\) denotes the feature representation of layer l with respect to the input, and f can be composed of the non-linear activation function and normalization operator (Ioffe and Szegedy 2015).

The FC layer establishes fundamental interconnections between two layers, inspiring diverse novel connectionism approaches, including the convolutional (LeCun et al. 1998), recurrent (Williams and Zipser 1989), and self-attention (Vaswani et al. 2017) approaches. The characteristic of shared local connection inherent in the convolutional layer can extract local features efficiently, making it an ideal choice for image-related tasks (Krizhevsky et al. 2012). Besides the convolutional layer, the deliberately designed recurrent connection (Hochreiter and Schmidhuber 1997) can describe the long-term relationship among the time steps, which is appropriate for sequential-related tasks. Moreover, the self-attention embedded in the transformer architecture (Vaswani et al. 2017) shows superiority in modeling the contextual information, breaking state-of-the-art records in various tasks. First applied in the NLP-related tasks, later leveraged in the image-related analysis with the proposition of vision transformation (ViT) (Dosovitskiy et al. 2020), the transformer has unified the analysis tasks of different modals.

While various novel DNNs have been proposed to address the analysis tasks, the DNNs are criticized for their weak safety, poor generalization, and lack of explainability. Regarding safety, studies have reported that DNNs are vulnerable to adversarial examples. In other words, a small perturbation to the input can lead to false output with high probability (Goodfellow et al. 2014). Formally, for a given sample x and its corresponding prediction y generated by DNNs. The model is considered safe if, for every \(\epsilon > 0\), there exists \(\delta > 0\) such that when \(\left\| x - x^{\prime } \right\| < \delta\), it follows that \(\left\| y - y^{\prime }\right\| < \epsilon\). Nevertheless, literature has observed that current layer-wise DNNs are unsafe (Zheng et al. 2016). This lack of safety significantly increases the deployment risk of DNNs, despite the promising performance reported in the in-house test dataset. Regarding generalization, DNNs often fail on out-of-distribution (OOD) data, attributed to the independent and identically distributed (IID) assumption of the datasets used in the training and test phases. The limited generalization ability of DNNs may restrict their practical application in medical image analysis tasks, as input data may originate from different centers, vendors, and diseases (Campello et al. 2021). For the explainability, the DNNs have been condemned for a long time due to their inherent black-box computational principle. It is difficult to explain the knowledge contained in the DNNs during their internal layer-wise forward computation. Limitations in the aforementioned three aspects impede the performance of DNNs in real-world scenarios.

Fig. 1
figure 1

Illustration of the trajectory of perturbed input \(x+\delta\) in the NODE

Full size image

Fortunately, the neural ordinary differential equations (NODE) (Chen et al. 2018), a specialized model that describes the DNNs from the viewpoint of the dynamical system, provides a potential solution for the above limitations. The NODE bridges the gap between the DNNs and ODE, changing the paradigm of DNNs from a discrete limited number of layers to continuous unlimited ones. Unlike well-known DNNs with explicit architecture, the NODE implicitly maps the input to the output. The NODE has been observed with higher non-linearity, clearer dynamics behavior, and stronger fitting capacity when compared with the discrete model. For the safety, the NODE evolves the features along with the trajectory embedded in the ODE, alleviating the impact of the noise or attacks contained in the input. Figure 1 provides an example of the trajectory of a disrupted sample in the NODE. For the disrupted sample \(x + \delta\), the vanilla DNNs may output the false prediction. However, the NODE can evolve the input’s representation along the trajectory defined in the dynamical system, resulting in a correct representation. Regarding the generalization, the advantage of NODE in safety also contributes to improving the generalization capability of DNNs. It has been reported that the DNNs embedded with NODE have better robustness compared with the conventional discrete DNNs (Anumasa and Srijith 2021). For the explainability, a variant of NODE termed neural memory ordinary differential equation (nmODE) shows that the ODE possesses strong explainability of the internal representation. By visualizing the internal representation of the ODE, Zhang (Yi 2023) shows that the samples belonging to the same class display high similarity, indicating that the status of the sample in the ODE is explainable.

Recognizing the strength of NODE mentioned above, literature has leveraged NODE and its variants to tackle the medical image analysis tasks, which are characterized by high-noise, volumetric, and multi-modality properties. Numerous encouraging results have been reported in the literature. However, there is a lack of a comprehensive review of these applications, impeding an in-depth understanding of NODE’s working principles and potential applications. To address this limitation, we systematically collected and summarized NODE-related papers in the field of medical image analysis. To the best of our knowledge, this is the first review of NODE in medical image analysis. To cover the field as broadly as possible, the applications are divided into five categories: segmentation, reconstruction, registration, disease prediction, and data generation. For each application, we review the primary principles of the work and summarize the improvements of NODE compared to control methods. Moreover, we also highlight the potential limitations of NODE and give possible research directions in the discussion section.

This paper is organized as follows. Section 2 illustrates the working principle of NODE and its variants nmODE, aiming to explain the fundamental mechanism of NODE and clarify the relationship between the NODE and DNNs. Section 3 presents the applications of NODE in medical image analysis tasks, including segmentation, reconstruction, registration, disease prediction, and data generation. Section 4 discusses the strengths and limitations of NODE during the applications. Finally, Sect. 5 summarize the review.

2 Deep neural networks and ordinary differential equation

Despite the fact that the NODE and DNNs belong to distinct continuous and discrete models, the two models are closely related. This section begins by demonstrating the origination of the NODE, followed by the introduction of the recently proposed nmODE, which exhibits superior theoretical and practical properties when compared to the vanilla NODE.

2.1 Neural ordinary differential equation

We begin with the demonstration of the residual module that forms the basis of NODE. The residual module proposed in the ResNet (He et al. 2016) is one of the most prevalent modules in DNNs. It is motivated by the degradation problem during the training of DNNs, i.e., the deeper DNNs have higher training and test errors compared with the shallow ones. The degradation problem limits the depth of DNNs and poses challenges for DNNs with very deep depth, even the depth contributes to the non-linearity. To tackle this limitation, ResNet (He et al. 2016) proposes to incorporate the shortcut connection into the DNNs, resulting in the residual paradigm that has been widely used in building DNNs’ architectures. The residual module can be formulated as follows:

$$\begin{aligned} a^{l+1} = a^{l} + f(a^l;\, W), \end{aligned}$$
(1)

where the \(a^l\) in the right side of Eq. (1) denotes the shortcut connection that facilitates the gradient backpropagation during training. The introduced shortcut connections make it possible to construct DNNs with significant depth, allowing for models with 1000 or more layers.

Considering the superscript layer l in Eq. (1) as the time step t and interpreting the feature \(a^l\) as the representation y(t), we arrive at the following equation:

$$\begin{aligned} y(t+1) - y(t) = f\left( y(t);\, W\right) . \end{aligned}$$
(2)

Taking the limit as t approaches infinity on the left side of Eq. (2), we obtain the NODE formulated as:

$$\begin{aligned} {\dot{y}}(t)=f \left( y(t);\,W \right) . \end{aligned}$$
(3)

So far, we have completed the transformation between the discrete residual module and autonomous NODE. The right-hand side of Eq. (3) is implemented by DNNs, which constitute the name of Neural ODE. The solution to Eq. (3) can be obtained by solving the following integration:

$$\begin{aligned} y(t_1)=y(t_0) + \int _{t_0}^{t_1}f \left( y(t);\,W \right)\, dt. \end{aligned}$$
(4)

Given the initial value \(y(t_0)\), we then can obtain the solution \(y(t_1)\) by using the ODE solver to solve the Eq. (4). Figure 2 demonstrates the working principle of NODE, which can be regarded as a specialized layer inserted into the DNNs. Suppose the representation in layer l of DNNs is denoted as \(a^l\), then the vanilla NODE uses \(a^l\) as the initial point \(y(t_0)\). By feeding the \(y(t_0)\) into the NODE described in Eq. (3) and using ODE solver to accomplish the integration, we can obtain the solution \(y(t_1)\), which can be regarded as the representation in layer \(l+1\).

Fig. 2
figure 2

Illustration of the architecture of NODE

Full size image

Modern DNNs are typically optimized through automatic differentiation. For the NODE, there are two ways to accomplish the optimization, namely the backpropagation or the adjoint (Chen et al. 2018) approaches, which are known as discretize-then-optimize and optimize-then-discretize (Kidger 2021) methods, respectively. Suppose there are T iteration steps in the ODE solver, the backpropagation approach can compute precise gradients and has the advantage of speed, but it consumes \({\mathcal {O}}(T)\) memory since it needs to store all the intermediate variables for backpropagation. Nevertheless, the adjoint approach has a constant memory cost of \({\mathcal {O}}(1)\), but it is slower compared to the backpropagation and introduces numerical discretization errors.

One significant characteristic of ODE is that the integral curves described in Eq. (4) do not intersect. Formally speaking, suppose \(y_1(t)\) and \(y_2(t)\) are two solutions of ODE with different initial values of the identical ODE function, then \(y_1(t) \ne y_2(t)\) for all \(t \in [t_0, t_1]\). The non-intersecting property implies that the integral curves of the starting point and its neighborhood are restricted. If a sample is correctly classified by a NODE, then the prediction of its perturbed version would be bounded. This intrinsic regularization in NODE contributes to increasing the robustness, which is absent from the discrete layer-wise DNNs.

Many variants have been proposed since the prevalence of NODE. ANODE (Dupont et al. 2019) is proposed to augment the dimensions of NODE to learn complex mapping. Gholami et al. (2019) observe that the NODE trained with adjoint method (Chen et al. 2018) may not converge due to inconsistent gradients for small time step sizes. The authors address the underlying problem by incorporating the checkpointing method while keeping the same computational cost as the NODE. Solving the NODE requires evaluating the differential equation multiple times, which is referred to the number of function evaluations (NFE). It is desirable to decrease the NFE without compromising accuracy, which can prompt the computational efficiency of NODE for reality applications. To reduce the NFE, Kelly et al. (2020) simplify the trajectories of NODE by introducing a differentiable regularization term composed of the Kth-order of the internal state with respect to the time. Kidger et al. (2021) use a seminorm to replace the frequently used \(L^2\) norm to reduce the NFE.

2.2 Neural memory ordinary differential equation

Fig. 3
figure 3

Illustration of the nmODE, whose inputs comprise two parts, including the initial internal input \(y(t_0)\) and external input \(\gamma\)

Full size image

One critical aspect of NODE is its input, which is frequently initialized by the output of the previous layer, e.g., the \(y(t_0)\) shown in Fig. 2. Dupont et al. (2019) points out that the NODE cannot represent some specialized function, for example, the function \(g(-1)=1\) and \(g(1)=-1\). The reason lies in the fact that the ODE’s trajectory cannot cross each other. Nevertheless, this restriction can be solved by reorganizing the inputs of NODE, i.e., regarding the feature as external inputs (\(\gamma\)) while using fixed the initial internal input (\(y(t_0)\)) (Yi 2023). For example, the one-dimensional ODE

$$\begin{aligned} {\dot{y}} = -y - x \end{aligned}$$
(5)

can represent the function g. The solution of the above ODE is

$$\begin{aligned} y(t) = \left( y(0) + x\right) e^{-t} - x, \end{aligned}$$
(6)

where we can verify the representation of g by rearranging the input as external input x.

This architecture proposed by Yi (2023) is named neural memory ordinary differential equation (nmODE), which contains two types of neurons, including the learning neurons and memory neurons. The separation of external inputs and internal inputs indicates that learning and memory are divided, where learning only happens in the learning neurons and the memory neurons attempt to capture the feature’s characteristics through ODE. The architecture of nmODE is illustrated in Fig. 3. Formally, the framework of nmODE can be formulated by the following equations:

$$\begin{aligned} {\left\{ \begin{array}{ll} &{} {\dot{y}}(t)=f \left( y(t),\gamma \right) \\ &{} \gamma =g(a^l;\,W). \end{array}\right. } \end{aligned}$$
(7)

Similar to NODE described in Eq. (3), the nmODE outputs the solution \(y(t_1)\) by providing the initial internal input \(y(t_0)\) and external input \(\gamma\). Both the NODE and nmODE can be integrated into DNNs as specialized layers. The most significant difference between NODE and nmODE lies in the learning location. In NODE, learning is integrated into the ODE, where the internal initial value y(0) is assigned as the output from the previous layer \(a^l\). Nevertheless, the nmODE separates the learning from the ODE, specifically the transformation \(g(a^l; W)\) that generates the external input \(\gamma\) for the ODE. The nmODE provides a versatile architecture for implementing implicit non-linear mapping. Yi (2023) introduces a novel implementation of Eq. (7) as follows:

$$\begin{aligned} {\dot{y}}(t) = -y(t) + \sin ^{2} \left( y(t)+\gamma \right) . \end{aligned}$$
(8)

It is widely acknowledged that the attractor in a dynamic system is associated with memory. Theoretically, the nmODE shown in Eq. (8) has been proved that it has one and only one global attractor (Yi 2023), indicating that the model possesses enhanced memory properties. Experiments on natural and medical image recognition tasks show that the nmODE helps to improve the model’s accuracy.

3 Applications in the medical image analysis

3.1 Applications review

Given the strength of NODE and its variants in nonlinear mapping, they have been widely utilized for various analysis tasks, including medical image analysis. To summarize their applications in the medical image analysis tasks, we have thoroughly searched the Web of Science and Google Scholar by using the keywords neural ordinary differential equation and medical image. Based on the search results, we organize the related works into five categories, including medical image segmentation, medical image reconstruction, medical image registration, disease prediction, and medical data generation. Illustrations of related works within each category are presented in the following five subsections.

3.2 Applications in medical image segmentation

Fig. 4
figure 4

Architecture of U-Net with ODE. The orange circle represents the potential applied location of ODE, including the encoder, decoder, and shortcut connections

Full size image

Medical image segmentation plays an important role in clinical practice, contributing to disease diagnosis, treatment planning, and prognostic analysis, etc. Starting from the breakthrough introduced by DNNs (Krizhevsky et al. 2012), especially the deep convolutional neural networks (DCNNs), the learning-based paradigm that extracts high-level abstract features has replaced the conventional hand-crafted low-level ones. The DNNs and their variants have become the de facto choice to accomplish the medical image segmentation tasks. The fully convolutional network (FCN) (Long et al. 2015) is the first attempt to build an end-to-end pixel-level semantic segmentation network. Among the variants of FCN, the U-Net (Ronneberger et al. 2015), a powerful segmentation network with succinct architecture, has been widely utilized to tackle varied medical image segmentation tasks. The composition of U-Net includes three parts, including an encoder, a decoder, and multiple shortcut connections, which is shown in Fig. 4. The encoder is designed to extract abstract features through multiple stages, where each stage consists of a convolutional layer, batch normalization operator, non-linear activation function, and pooling layer. The decoder has a similar architecture to the encoder, except that the pooling layer is replaced by the upsampling or interpolation ones. The deliberately designed shortcut connections contribute to the information flow between the encoders and decoders. Numerous variants of U-Net have been proposed, such as the Attention U-Net (AttU-Net) (Oktay et al. 2018), Recurrent Residual U-Net (R2U-Net) (Alom et al. 2018), and no new U-Net (nnU-Net) (Isensee et al. 2021) etc.

Fig. 5
figure 5

Various modalities used in the medical image segmentation tasks modeled by the NODE. The four images in the first row from left to right are T1, T1 CE, T2, and FLAIR

Full size image

Since the introduction of NODE, plenty of works have incorporated it into the U-Net, aiming to increase the segmentation accuracy. Figure 5 summarizes the various medical image modalities used in segmentation tasks within the literature. The first work using NODE to accomplish the medical image segmentation task can be found in Pinckaers and Litjens (2019), which is termed U-NODE, a combination of U-Net and NODE. The U-NODE applies the NODE to the encoder and decoder parts of U-Net to segment whole-slide histology images of the colon (Sirinukunwattana et al. 2017). The motivation of U-NODE is to leverage the continuous depth property of NODE to dynamically adjust the receptive field of U-Net, for a better accommodation of the target with varied sizes. Experimental results reveal that the U-NODE exhibits superiority in the segmentation noise of large targets compared with the vanilla U-Net. Following the U-NODE, Li et al. (2021) use NODE with U-Net to segment blood cells from microscopy images. Considering the computational cost associated with the NODE, Li et al. tested the applied location of NODE in the U-Net to balance the computational efficacy and segmentation accuracy. A similar work can be found in Sadique et al. (2022), where Sadique et al. employ the NODE in U-Net with context encoding (Rahman et al. 2021) to segment multiple tumor tissues in the BraTS 2021 dataset (Baid et al. 2021), which contains magnetic resonance imaging (MRI) dataset with different parametrics, including the T1, T1 contrast enhancement (T1 CE), T2, and fluid-attenuated inversion recovery (FLAIR). Evaluation results on the out-of-sample MRI dataset differ from the BraTS showing that the NODE helps improve the generalization capability of the segmentation network.

Unlike the Sadique et al. (2022) that stacks the multi-parametric MRI along the channel dimension as input, Yang et al. (2023) utilizes a shared U-Net with NODE to process different modalities independently. The final tumor segmentation result is obtained through a weighted summation of the four branches corresponding to the four MRI modalities. Moreover, the authors deliberately designed a metric termed accumulative contribution curve to inspect the contribution of each modality to the segmentation of different tumor types, including the enhancing tumor (ET), tumor core (TC), and whole tumor (WT). Experimental results reveal that not all modalities have a significant impact on the segmentation performance, indicating a reduced number of MRI modalities for automatic segmentation. This reduction is meaningful for both the clinical and computational modeling, contributing to efficient clinical practice and simplified model architecture. Besides the brain tumor, Ru et al. (2023) tackle the breast tumor segmentation with a proposed model named Att-U-NODE that integrates the attention mechanism (Oktay et al. 2018; Woo et al. 2018) with U-Net applied with NODE. Two datasets are used in the experiments, including the breast ultrasound and MRI images. The attention module is applied to the shortcut connections of the U-NODE, aiming to enhance the representation ability of abstract features.

While the fully-supervised U-Net with NODE delivers promising performance on various segmentation tasks, acquiring large-scale dense pixel-level segmentation datasets is quite expensive, especially in the context of medical imaging. Few-shot segmentation (FSS) (Wang et al. 2019) provides a potential solution that accomplishes segmentation for novel classes in a query image trained on a few annotated support images. However, the robustness of the FSS against the adversarial attack is unclear, which motivates (Pandey et al. 2023) to combine FSS with NODE to build a data-efficient segmentation model that is both accurate and safe. The method proposed by Pandey et al. (2023) is named regularized prototypical neural ordinary differential equation (R-PNODE), designed to leverage the intrinsic properties of NODE for FSS of organs in computed tomography (CT). For the paired query and support images denoted as \((x_q, x_s)\), the authors first obtain their corrupted version by applying Gaussian noise that is formulated as \((x_q^G, x_s^G)\). Three losses are proposed in the R-PNODE, with two designed for the query image and the third one for the support image. For the query image, suppose the predictions of the clean and Gaussian versions from the NODE are denoted as \(y(x_q, t_1)\) and \(y(x_q^G, t_1)\), then R-PNODE minimizes their distance between the groundtruth, which is named as consistency losses. Regarding the support image, R-PNODE maximizes the cosine similarity between the prediction of the clean support image and its related Gaussian one, which is named cluster loss. The cluster loss contributes to a robust representation of the input, aiding the FSS model in defending against adversarial attacks.

Recently, Hu et al. (2023) empirically observe that the conventional DCNNs-based segmentation model is sensitive to the noise in the training dataset, revealing inherent limitations in the robustness of DCNNs. Nevertheless, the U-Net with nmODE (Yi 2023), a novel variant of NODE, helps to alleviate the impact of the noise. The continuous nmODE and discrete DCNNs belong to two different types of DNNs, implying the distinct knowledge contained in the networks. The discrete DCNNs own strength in speed, but are vulnerable to the noise in the training dataset; conversely, the NODE is superior in the aspect of robustness, but with the cost of increased computational complexity. To balance the speed and robustness, Hu et al. (2023) leverage the knowledge distillation (KD) (Hinton et al. 2015) framework to transfer the knowledge from continuous nmODE to discrete DCNNs, simultaneously enhancing the model’s robustness and efficiency. The main idea of the proposed method can be formulated as the following equation:

$$\begin{aligned} {\mathcal {L}}= -(1-\lambda )\log {p_k^s} - \lambda \sum _{c=1}^{C}p_c^{t} \log {p_c^{s}}, \end{aligned}$$
(9)

where the two items in the right-hand side represent the cross entropy and KL divergence balanced by the hyper-parameter \(\lambda\). Superscript s and t denote the student (discrete U-Net) and teacher (continuous nmODE) models, respectively. The \(p_k^s\) represents the prediction from the discrete model regards the target class k. The \(p_c^{t}\) and \(p_c^{s}\) represent the predictions of the nmODE and the U-Net model for the class c, respectively. Experimental results on 18 organs-at-risk (OAR) segmentation tasks demonstrate that nmODE with knowledge distillation exhibits improved robustness compared to ODE-based models while also mitigating the computational efficiency limitation.

The NODE and its variants presented above belong to the first-order method, while Cheng et al. (2023) propose a second-order approach. They theoretically proved and empirically verified that the second-order method has faster convergence, higher robustness, and less sensitivity to noise compared to the first-order ones. The second-order NODE can be formulated as follows:

$$\begin{aligned} {\left\{ \begin{array}{ll} &{} y''(t)=f \left( y(t), y'(t); \,\theta _f \right) \\ &{} y'(t) = g \left( y(t);\, \theta _g \right) , \\ \end{array}\right. } \end{aligned}$$
(10)

where \(\theta _g\) and \(\theta _f\) represent the parameters of the first-order and second-order NODE, respectively. Cheng et al. (2023) compare characteristics of different ODE solvers, and show that the RK4 method is appropriate for the second-order segmentation model. Extensive experiments on six segmentation datasets verify the effectiveness of the second-order NODE model.

Fig. 6
figure 6

Segmentation maps of different time points predicted by the NODE (Lachinov et al. 2023). The upper and lower present the targets related to the Geographic Atrophy and Alzheimer’s disease, respectively. The figures are from Lachinov et al. (2023)

Full size image
Table 1 Comparison of segmentation’s performance between the NODE and control method
Full size table

Despite the static 2D or 3D medical images, the longitudinal imaging that captures multiple scans of a patient along with the time helps the doctor to observe the progress of the disease. The computational method that analyzes the target changes by using longitudinal imaging is named disease progression models (DPMs) (Lachinov et al. 2023). While the recurrent neural networks (RNNs) (Zhang 2013) show strength in modeling the time-related dataset, they are designed for forecasting the fixed intervals only, which fails in analyzing the time-variant dataset. NODE that models the dynamic physical process provides a potential solution for the DPMs. Lachinov et al. (2023) leverage the NODE to forecast the segmentation map at arbitrary time points. The first visit image is processed by a CNN to extract the feature representation \(\xi\), which is used as the initial condition \(y(t_0)\) for the NODE that evolves over time to form the trajectory of the segmentation map. The trained NODE can then be sampled at arbitrary time points to obtain the segmentation map. Moreover, the authors also propose a temporal Dice loss designed for the entire stack of visits. Experimental results on the optical coherence tomography (OCT) of geographic atrophy and MRI of Alzheimer’s disease demonstrate the advantage of NODE in modeling the time-variant dataset. Figure 6 presents the segmentation map of varied time points predicted by the NODE. Table 1 summarizes the comparison between the NODE and control methods in various segmentation tasks.

3.3 Applications in medical image reconstruction

Target reconstruction from medical images is another important task in medical image analysis. For example, reconstructing cortical surfaces from MRI contributes to quantitative analysis of the brain status of different ages. Nevertheless, the reconstructed cortical surface mesh requires to be a closed manifold, which is topologically equivalent to a 2-sphere. This means that the reconstruction results should be without any self-intersection cases, i.e. without any holes or handles. Despite the progress achieved by classical neuroimage analysis software, such as FreeSurfer (Fischl 2012), BrainSuite (Shattuck and Leahy 2002), HCP pipeline (Glasser et al. 2013), two notable limitations restrict their applications (Ma et al. 2022). The first is time consumption. The mentioned conventional methods are time-consuming since several computationally intensive geometric or image processing algorithms are involved. For example, the FreeSurfer requires more than 4 h to accomplish the reconstruction of the cortical surface for a single patient. The second limitation is the generalization capability. The brain images from neonatal and adults significantly varied in intensity values, shapes, and sizes, which pose challenges to the robustness of reconstruction methods.

Fig. 7
figure 7

Comparison of the reconstructed cortical surface from the U-Net, NMF, and CorticalFlow (Lebrat et al. 2021). The figures are from Lebrat et al. (2021)

Full size image

Recently, end-to-end deep learning methods have been widely undertaken to solve the cortical surface reconstruction task. For instance, FastSurfer (Henschel et al. 2020) first leverages the DCNNs to segment the brain tissue, which is later reconstructed to the cortical surface by adopting a fast spherical mapping. DeepCSR (Cruz et al. 2021) uses the occupancy fields and signed distance functions (SDFs) to predict the implicit cortical surface. These implicit methods require topology correction to guarantee the non-intersection property of the reconstructed surface. Besides the implicit methods, Voxel2Mesh (Wickramasinghe et al. 2020) and PialNN (Ma et al. 2021) belong to the explicit approaches that directly reconstruct the cortical surface without post-processing. These explicit approaches show superiority in speed but lack theoretical guarantees to prevent self-intersections of the surface.

Faced with the aforementioned limitations of both classical and DCNNs-based methods, Lebrat et al. (2021) propose the CorticalFlow, a NODE-based approach that reconstructs the cortical surface from input MRI images. Three steps are included in the overall process. The first step involves generating the template mesh from the input MRI volume using the third-party libraries, including the JIGSAW (Engwirda and Ivers 2016) and MeshLab (Cignoni et al. 2008). In the second step, 3D U-Net (Ronneberger et al. 2015) is employed to generate a dense 3D flow field from the MRI volume, with moderate GPU memory consumption. The third step involves feeding both the predicted flow field and template mesh input into an ODE-based model named diffeomorphic mesh deformation (DMD) to compute the diffeomorphic mapping for each vertex. Moreover, a multi-stage framework with an increased number of vertices is proposed, progressively refining the predicted cortical surface from low to high resolution.

The most important aspect of CorticalFlow is the continuous DMD module, which theoretically guarantees that the predicted cortical surface is non-intersecting. The DMD module is solved using the forward Euler method. Figure 7c compares the reconstructed cortical surface from the U-Net (Ronneberger et al. 2015), neural mesh flow (NMF) (Gupta 2020), and the CorficalFlow (Lebrat et al. 2021). It is clear that the U-Net produces irregular mesh, as shown in Fig. 7a, which can be attributed to the absence of geometric regularization of the U-Net. The NMF is another NODE-based method designed for 3D mesh generation (Gupta 2020). However, Lebrat et al. (2021) found that the NMF conditioned on the global feature descriptor of the input image failed to provide the details of the cortical surface, as observed in Fig. 7b. Nevertheless, CorticalFlow uses the DMD module to optimize the local feature descriptor of each vertex and employs three stages to gradually increase the number of processed vertices. Figure 7c shows that the CorticalFlow generates anatomically plausible cortical surfaces with clear boundaries.

Fig. 8
figure 8

The diffeomorphic flow of the CortexODE from the leftmost white matter to the rightmost pial surface (Ma et al. 2022). The figures are from Ma et al. (2022)

Full size image

Motivated by the promising reconstruction results reported by the CortexFlow, Ma et al. (2022) have further proposed a NODE-based cortical surface reconstruction framework named CortexODE. The CortexODE comprises three steps. In contrast to the CortexFlow, which uses the 3D U-Net to predict the dense flow field, CortexODE initially segments the raw white matter from the input MRI volume. In the second step, the raw white matter is transformed into a signed distance function (SDF), subsequently processed by a topology correction algorithm to obtain an initial closed manifold surface. In the final third step, the initial surface is deformed to the white matter and pial surface through two NODEs, respectively. Figure 8 illustrates the diffeomorphic flow modeled by the CortexODE from the leftmost of the white matter surface to the rightmost pial surface.

Besides the Cortical surface reconstruction task modeled by the NODE, reconstructing the fully sampled medical image (e.g., MRI) from the undersampled ones contributes to balancing the efficiency and effectiveness of the imaging process. It is acknowledged that a major challenge of MRI lies in its slow data acquisition process, which is attributed to the hardware restriction. To expedite the process of MRI acquisition, it is common to acquire an undersampled k-space, leading to aliasing artifacts in the resultant image domain. The reconstruction of high-fidelity images from the undersampled k-space data holds great significance in the clinical utilization of MRI.

The MRI acquisition can be formulated by the following equation:

$$\begin{aligned} x = Ey + \epsilon , \end{aligned}$$
(11)

where y denotes the fully sampled MRI image to be reconstructed, and x represents the observed undersampled k-space. \(\epsilon\) is the noise during the imaging process. E represents the Fourier transform and undersampling that convert the fully sampled MRI into the undersampled k-space. Thus, the MRI reconstruction task can be described as the following optimization problem:

$$\begin{aligned} \underset{y}{\textrm{argmin}} \frac{1}{2} \left\| x - Ey \right\| +R(y), \end{aligned}$$
(12)

where the left term measures the distance between the reconstructed and observed MRI image. The right one is the regularization term. The problem shown in Eq. (12) can be solved using the gradient descend algorithm:

$$\begin{aligned} y(n+1) = y(n) - \alpha \left[ E^T \left( Ey(n) - x \right) + \nabla R \left( y(n) \right) \right] , \end{aligned}$$
(13)

where the variable n denotes the iteration step and \(\alpha\) is the learning rate. The Eq. (13) can also be abbreviated as:

$$\begin{aligned} y(n+1) - y(n) = f(y(n), x; \,W). \end{aligned}$$
(14)

W represents the learnable parameters during the optimization. Similar to the equation described in Eq. (2), by considering the step on the left side of Eq. (14) to be infinitely small, we arrive at the NODE in Eq. (3). So far, we can see that the MRI reconstruction task is closely related to the NODE, which can be solved by off-the-shelf solvers.

Among variants of the ODE solver, the Euler solver is the first-order one, which iteratively approximates the integral of ODE through a fixed step size h:

$$\begin{aligned} y(n+1) = y(n) + h * f(y(n), t(n)). \end{aligned}$$
(15)

Advanced solvers, such as the Runge–Kutta (RK) methods with higher orders, can achieve better precision than the Euler approach. The computational principle of the general RK methods with stage s can be formulated by the following equations:

$$\begin{aligned} {\left\{ \begin{array}{ll} &{} y(n+1) = y(n) + \sum \limits_{i=1}^{s}{a_i F_i} \\ &{} F_1 = h * f(y(n), t(n)) \\ &{} F_i = h * f(y(n) + \sum \limits_{j=1}^{i-1}b_{ij}F_j, t(n) + c_i h), i=2,...,s, \end{array}\right. } \end{aligned}$$
(16)

where \(a_i, b_{ij},\) and \(c_i\) are fixed coefficients. Instead of the fixed step size and coefficients, Chen et al. (2020) propose integrating knowledge of the solvers into the network, letting the network learn to adjust the step size and coefficients according to the requirements of tasks. The gradient checkpoint method (Chen et al. 2016) is used to reduce the GPU memory consumption. Experimental results show that the learned ODE solver achieves better reconstruction accuracy compared with the fixed one. Figure 9 compares the reconstruction results from U-Net and ReconODE, where it is clear that the error of ReconODE is smaller than that of U-Net. Similar work for the MRI and cone-beam CT reconstruction can be found in Yazdanpanah et al. (2019); Thies et al. (2022). Table 2 summarizes the comparison between the NODE and control methods in various reconstruction tasks.

Fig. 9
figure 9

Comparison of reconstructed MRI from U-Net (Ronneberger et al. 2015) and ReconODE (Chen et al. 2020). The figures are from Chen et al. (2020)

Full size image
Table 2 Comparison of reconstruction’s performance between the NODE and control method
Full size table

3.4 Applications in medical image registration

Deformable registration computes a dense correspondence between two images, which plays an important role in various medical image analysis tasks, e.g., surgical planning and radiogenomics analysis. The term deformable presents the nonlinear and dense characteristics of the transformation. Let \(y_{\textrm{mov}}\) and \(y_{\textrm{fix}}\) denote the moving and fixed images, where the registration task aims to find the deformation field \(\phi _{W}\) parameterized by W to convert the moving image similar to the fixed one (Dalca et al. 2018). The registration problem can be formulated as:

$$\begin{aligned} {\hat{W}} = \underset{W}{\textrm{argmin}} {\mathcal {L}}(y_{\textrm{mov}} \circ \phi _{W}, y_{\textrm{fix}}) + {\mathcal {R}}(W), \end{aligned}$$
(17)

where the term \(y_{\textrm{mov}} \circ \phi _{W}\) denotes the transformed moving image. Function \({\mathcal {L}}\) measures the similarity between the transformed moving and fixed images. \({\mathcal {R}}(W)\) denotes the regularization term to ensure the smoothness of the transformation. The architecture of the registration method between the moving and fixed images can be found in Fig. 10.

Fig. 10
figure 10

The architecture of learning-based registration of the moving and fixed images

Full size image

Even the conventional registration methods, such as elastic-type models (Bajcsy and Kovačič 1989), B splines (Rueckert et al. 1999), or dense vector fields (Thirion 1998), have achieved impressive registration performance, they are criticized for the slow registration process due to the time-consuming iterative non-linear optimization algorithms. The computational burden hinders their clinical translation. Faced with this challenge, both supervised and unsupervised learning-based methods have been proposed to accomplish the registration task. The learning-based approach requires one forward pass of the model, which is intrinsically faster than the conventional methods. For the supervised approach, the registration model maps the moving and fixed images to an output deformation. The groundtruth deformation field is often obtained through conventional registration tools, which would introduce biases. The unsupervised registration method that does not require the groundtruth deformation field is more practical than the supervised one (Dalca et al. 2018).

Xu et al. (2021) propose an ODE-based registration method named Multi-Scale ODE Network (MS-ODENet), marking the pioneering application of NODE to medical image registration. The authors model the registration optimization problem in Eq. (17) as a continuous dynamics process using NODE. Throughout the experiments, the authors observed that solving NODE for the image registration problem demands a large number of NFE, significantly increasing the training and inference time costs. To address the computational efficiency problem, they propose a multi-scale approach that transforms the input image into L different resolutions. The \((l-1)\)-th version is generated by down-sampling the l-th with a factor of 2 on the x, y, and z axes. The integral interval of NODE, suppose as [0, T], is divided into L segments corresponding to the number of resolutions. In each segment, the NODE is solved at the corresponding resolution, forming the basis of the multi-scale ODE. The MS-ODENet offers two advantages compared to the vanilla one. First is the time cost. The time consumption for the function evaluation of the low-resolution input is much smaller than the high-resolution one. The second advantage is convergence. The search space of the low-resolution input for image registration is greatly reduced, leading to a faster convergence speed.

In addition to the multi-scale approach, the authors propose pretraining the feature extractor. This is motivated by the fact that the anatomical structures of the same patient remain consistent across different contrast images. The authors utilize the image-to-image translation framework (Lee et al. 2020) to separate the content from the style features. Subsequently, the robust pretrained content encoder is transferred to the registration network.

Wu et al. (2022) propose a NODE-based Optimization (NODEO) framework that regards the voxels in the moving image as a high-dimensional dynamical system whose trajectory represents the deformation field. An important contribution of the NODEO lies in the implementation of the regularization term in Eq. (17), which is composed of three parts. The first is the negative Jacobian determinants of the transformed voxels, which aim to reduce the folds in the transformation. The second is the regularization of the magnitude of the velocity field, which is equivalent to penalizing the energy of the transformation. The third is the spatial gradient of the voxels after the deformation, which contributes to the spatial smoothness of the transformed voxels.

Fig. 11
figure 11

Comparison of registration results from VoxelMorph (Balakrishnan et al. 2019) and R2Net (Joshi and Hong 2023). The figures are from Joshi and Hong (2023)

Full size image

Joshi et al. propose a residual registration network named R2Net (Joshi and Hong 2023), which is an unsupervised registration framework that leverages residual connections to generate the desired deformation field. The framework is composed of three components. The first one involves generating the initial velocity field (\(v_0\)) by using the U-Net with the input of the moving and fixed images. In the second step, the generated initial velocity field is integrated into the deformation field (\(\phi\)) with a module termed Lipschitz Continuous ResNet (LC-ResNet), contributing to obtaining the diffeomorphic deformation. During implementation, spectral normalization (Miyato et al. 2018) is applied to each convolutional layer in the residual blocks to guarantee the Lipschitz continuity. Moreover, shared or unshared convolutional weights are considered, resulting in the stationary velocity field or time-varying velocity fields, respectively. In the third step, chunk and downsampling, two frequently used data reduction strategies are adopted, leading to the local and global two branches. For the chunking branch, multiple velocity fields are merged to maintain the size the same as the input. For the downsampling branch, the downsampled velocity field is upsampled back to the original resolution. The two velocity fields are then concatenated to generate the deformation field using the LC-ResNet. Figure 11 compares the registration results from VoxelMorph (Balakrishnan et al. 2019) and R2Net (Joshi and Hong 2023), clearly showing that R2Net achieves a smaller error than VoxelMorph. Table 3 summarizes the comparison between the NODE and control methods in various registration tasks.

Table 3 Comparison of registration’s performance between the NODE and control method
Full size table

3.5 Applications in disease prediction

Fig. 12
figure 12

The architecture of the NODE-based model for the disease prediction tasks

Full size image

The disease prediction task is a crucial aspect of medical image analysis, encompassing two primary domains. The first domain involves identifying the presence of a disease by analyzing medical images, often using sequences or static images in clinical scenarios. The second one focuses on predicting the progression of the disease by utilizing various medical imaging modalities to model morphological or pathobiological changes over time.

Recent advances in DNNs offer substantial promise in addressing these disease prediction tasks. Notably, RNNs have gained popularity for their ability to handle time-dependent patterns in patients’ progression. The recurrent connections allow the network to process input data at each time step, continuously updating its hidden state, which acts as a memory mechanism to model the patterns in sequential data. However, original RNNs face the vanishing gradient problem, where gradients become exceedingly small when modeling long sequences. To overcome this challenge, long short-term memory (LSTM) networks (Memory 2010) and gated recurrent Units (GRUs) (Cho et al. 2014) have emerged as prominent architectures. These architectures incorporate effective gating mechanisms, efficiently mitigating the vanishing gradient problem.

RNNs and their variants are effective in modeling discrete-time dynamical systems. However, clinical data often introduces uncertainty in the time interval between data points, e.g., incomplete data or irregular acquisition times. To address this issue, previous studies either removed several data points to maintain evenly spaced data or employed missing value imputation techniques to generate possible values for intermediate data. However, these methods frequently introduce errors between the generated data and the real values, impacting the prediction results.

In contrast, NODE, the powerful continuous dynamical system, can effectively model data with varied time intervals. Unlike the RNNs, NODE is insensitive to inconsistencies in data points, making it an ideal choice to address the irregular time intervals in clinical data. Figure 12 illustrates the architecture of the NODE-based model for disease prediction tasks. Qiu et al. (2023) propose a gram-based attentive neural ordinary differential equations network (Gram-AODE) method. This approach aims to model the underlying dynamics of eyeballs’ trajectory points and learn the relations of continuous changes. Specifically, the authors first transform the video into a time series of pupil coordinates. The Gram matrix is then employed to convert the time-series data into feature images. After downsampling, the feature images are further augmented by the NODE. The dynamically modeled features by NODE are then processed through a multi-head attention block to output the categories.

The prediction of disease progression is a common requirement in analyzing brain diseases, such as Alzheimer’s disease (AD). Recently, Jeong et al. (2022) have proposed a NODE-based method to tackle challenges in modeling the progression of AD using longitudinal data characterized by irregular and incomplete observations. The workflow involves four steps: transforming input data into the Cholesky space, leveraging recurrent gated recurrent units (RGRUs) for dynamic modeling, employing neural manifold ODEs for continuous modeling, and estimating missing values through a proposed missing value estimation module. The NODE is utilized in their method to model time-series data and estimate the hidden state of missing values.

Hao et al. (2020) aim to utilize NODE to infer the propagation dynamics of amyloid pathology, which is considered the primary pathological event in AD. The method involves constructing a brain network graph where each node represents a cortical or subcortical gray matter region, and edges denote connectivity strength. The authors formulate the NODE to describe the propagation of pathology, incorporating a diffusion constant and a Laplacian-based heat diffusion equation. Additionally, they introduce a function learned through DNNs to modulate the diffusion process, enabling the prediction of pathology burden at multiple time points. The proposed model is evaluated on a neuroimaging dataset, demonstrating promising results.

Cai et al. (2023) propose a NODE-based brain state recognition neural network (OSRNet) to explore NODE’s performance on neuroimaging data. The network comprises three main components: feature Elucidation, a NODE Block, and a fully-connected classifier. Feature Elucidation involves reducing the dimensionality of high-dimensional data by mapping it to a Riemannian manifold. The method employs a generic mapping layer and feature clarification operations to obtain low-dimensional functional connectivity matrix signatures. NODE is used to model brain dynamics. A Log-Euclidean operator is applied to project symmetric positive definite matrices onto a vector space, and an ODE solver is used to simulate continuous-time dynamics. Moreover, a hybrid Euler method is proposed for numerical ODE solving. At each time step, it calculates the system’s state update using the explicit Euler method to capture the slope of the ODE.

In light of the existing unlabeled data situation, Zeghlache et al. (2023) propose an innovative fusion of longitudinal self-supervised learning (LSSL) and NODEs to enhance the modeling of disease progression. The authors introduce a Siamese-like LSSL variant designed to assess the role and significance of the reconstruction term within the LSSL framework. Furthermore, they contextualize the application of NODEs in disease progression modeling, drawing parallels with instances such as predicting Alzheimer’s disease progression and understanding patient dynamics with COVID-19 medication. The introduced LSSL-NODE model extends the application of predicting disease progression. This extension involves formulating an initial value problem through the latent representation of the first image in a consecutive pair, allowing for a more comprehensive exploration of disease dynamics.

Wen (2020) introduces a pioneering approach to modeling the temporal dynamics of resting-state functional magnetic resonance imaging (rsfMRI) data by leveraging NODEs. The key innovation lies in using NODEs as the foundation for video prediction, allowing the compression of spatial-temporal trajectories into latent representations. Specifically, Wen et al. employ NODE to model the spatial-temporal characteristics of rsfMRI data. They fit the observed rsfMRI data into the network to predict future fMRI data based on input trajectories. The authors claim this framework facilitates the prediction of future brain maps and enables data interpolation between two given trajectories. Notably, the utilization of NODE addresses the challenge of degradation in prediction quality over time points, presenting a promising solution to the inherent uncertainty in video prediction problems. Table 4 summarizes the comparison between the NODE and control methods in various disease prediction tasks.

Table 4 Comparison of prediction’s performance between the NODE and control method
Full size table

3.6 Application in medical data generation

Medical data generation refers to the simulation of medical data using synthetic methods to replace or complement real ones. Research in this field involves using generated data to simulate various scenarios in the medical domain, including medical images, clinical data, electronic health records, etc. With the advent of DNNs, learning-based methods have emerged as a prominent tool for data generation.

Medical data involves complex biological and medical dynamic systems, such as disease progression and treatment response. Accurately modeling these systems is challenging, especially when multiple variables and influencing factors are involved. Existing methods attempting to simulate such dynamic systems require a substantial number of learnable parameters and computational resources. To address this issue, NODE has been applied in the field of data generation for the following reasons. First is continuous dynamics modeling. NODE is designed to model continuous dynamics. This makes it well-suited for tasks where the data evolves continuously over time. Second is implicit function representation. NODE implicitly represents the dynamics through trajectory. This implicit representation can capture long-term dependencies and intricate relationships within the data, enabling NODE to generate coherent and realistic sequences. Third is interpolation and extrapolation. Through its continuous representation, NODE can both interpolate between observed data points and extrapolate to generate data at unobserved time points. This is valuable for tasks where generating data at arbitrary time steps is necessary.

Recent endeavors (Hong et al. 2023; Salvador et al. 2023; Wendland et al. 2022) have attempted to incorporate NODE into the medical data generation domain. Hong et al. (2023) present a NODE-based method for forecasting pharmacokinetic parameters in dynamic brain positron emission tomography (PET) imaging. Traditional methods, such as interpolation or extrapolation, require sufficient time-activity curve (TAC) samples throughout the entire acquisition, which is not always practical due to patient motion and noise. In contrast, the proposed NODE method predicts TAC in extended time frames by mimicking the analytical modeling in a data-driven manner. The method predicts complete dynamic images based on limited early-frame images, demonstrating promising results in both simulated and clinical \(^{18}\)F-PI-2620 brain PET data. The potential applications of this method include reducing scan time, minimizing radiation dose, improving imaging quality, and better evaluating drug efficacy.

Wendland et al. (2022) introduced a hybrid, multimodal approach called Multimodal Neural Ordinary Differential Equations (MultiNODEs), which enables the generation of lifelike synthetic patient trajectories over time. This feature facilitates the seamless interpolation and extrapolation of clinical studies. MultiNODEs is proficient in incorporating both static and longitudinal data, implicitly handling missing values. By using MultiNODEs, organizations can generate synthetic patient trajectories that capture the characteristics of real patient data, enabling them to conduct analyses and make decisions based on realistic data without compromising individual privacy or data-sharing restrictions. The authors illustrate the capabilities of MultiNODEs by applying them to actual patient-level data from two independent clinical studies. Simultaneously, they demonstrate its effectiveness with simulated epidemiological data of infectious diseases, emphasizing its potential to address the issue of restricted data collection and sharing among diverse organizations. Table 4 compares the generation performance between the MultiNODEs and control method.

Salvador et al. (2023) have proposed a scientific machine-learning methodology for creating a digital twin model of the heart. The digital twin model is a virtual model based on mathematical models and real-time data. It can simulate and predict the behavior and performance of a physical system. Utilizing latent neural ordinary differential equations (LNODEs), the authors capture the pressure-volume dynamics of a heart failure patient by training a system with 43 model parameters. These model parameters are used to describe cardiac electrophysiology, active and passive mechanics, and cardiovascular fluid dynamics. LNODEs efficiently represent the 3D-0D model in a latent space via a feedforward fully connected network. Remarkably, the model achieves 300x real-time numerical simulations of the cardiac function on a standard laptop processor (Table 5).

Table 5 Comparison of data generation between the NODE and control method
Full size table

4 Discussion

From the above review of NODE in the medical image analysis tasks, it can be found that NODE offers a powerful tool to address the challenges in conventional DNNs. The strengths of NODE can be categorized into practical and theoretical aspects. Practically, the NODE exhibits higher nonlinearity than conventional discrete DNNs, which is a key factor contributing to the successful applications of DNNs. The NODE is closely related to the residual connections in ResNet (He et al. 2016), where the latter is a milestone in the numerous architectures of DNNs that makes it possible to construct the network with very deep depth. Considering the layer in residual connections as the time step and approaching this step to be infinitely small, we arrive at the NODE. Therefore, NODE can be regarded as an advanced ResNet with infinite depth, characterized by higher nonlinearity and stronger representation ability. Theoretically, NODE derives benefits from the theoretical characteristics of ODE, particularly the non-intersecting theorem that ensures the robustness of NODE. Moreover, the theoretical advantage of NODE over the discrete DNNs contributes to the diffeomorphic transformation in the registration and reconstruction of medical images.

4.1 Computational limitation of NODE

While the NODE exhibits notable advantages over conventional DNNs, it is essential to acknowledge its inherent weaknesses. One prominent limitation is its computational cost. It has been observed that the model embedded with NODE is considerably slower than the conventional DNNs during both the training and inference phases (Hu et al. 2023). This slowness is attributed to the iterative computation required for the ODE solver’s integration, referred to as NFE. The introduced computational cost challenges the clinical deployment of NODE-based models, particularly in scenarios that necessitate timely analysis results.

4.2 NODE-based foundation model

In terms of future works, potential research directions include the following aspects. The first one is the NODE-based foundation model. The foundation model (Moor et al. 2023) contributes to the exceptionally rapid development of the generative large language model (LLM). By training the foundation model on a massive and diverse dataset through the self-supervised paradigm, the model exhibits impressive performance on various generative tasks. The NODE can be regarded as a specialized layer with higher nonlinearity and infinite depth. It is worthwhile to endow the foundation model with the NODE, which may further enhance the model’s capacity.

4.3 Pre-training and fine-tuning of NODE

The second one is NODE-based pre-training. The pre-training and fine-tuning paradigm has become the de facto standard for utilizing DNNs, particularly in medical image analysis tasks that often suffer from limited annotated datasets. Due to the superior memory properties of NODEs compared to discrete DNNs, transferring pretrained NODE-based medical image analysis models to downstream tasks may effectively reduce the need for annotated datasets and ease the difficulty of fine-tuning.

4.4 Balance of accuracy and efficiency

The third one is the balance between the computational cost and the accuracy of NODE-based models. Relaxing the error tolerance of NODE during the integration would expedite efficiency but at the cost of decreased accuracy and vice versa. Thus, finding a promising solution to balance efficiency and accuracy would significantly broaden the applications of NODE-based models.

5 Conclusion

This paper comprehensively reviews the applications of NODE-based models in the medical image analysis tasks. Numerous applications are broadly categorized into five types, including segmentation, reconstruction, registration, disease prediction, and data generation. Various novel NODE-based models have been proposed to address the challenges in these analysis tasks, delivering impressive performance compared to conventional discrete DNNs. Additionally, we summarize the strengths and weaknesses inherent in NODE-based models and suggest possible future research directions by integrating the evolution of modern artificial intelligence. We hope this review will provide insights into solutions for the medical image analysis tasks and capture researchers’ attention to the significant potential of NODE.