Modality-agnostic, patient-specific digital twins modeling temporally varying digestive motion

TL;DR

This paper introduces a pipeline to create digital twins of the digestive system to evaluate the accuracy of image registration methods in dynamic, complex regions.

Contribution

A novel pipeline for generating patient-specific digital twins of temporally varying digestive motion for DIR validation.

Findings

Digital twins achieved motion amplitudes and log Jacobian determinants comparable to real patient data.

The pipeline supports detailed quantitative evaluation of DIR performance with spatial and dosimetric metrics.

Dose warping errors were assessed in both low- and high-dose regions for patient-specific error estimation.

Abstract

Objective. Clinical implementation of deformable image registration (DIR) requires voxel-based spatial accuracy metrics such as manually identified landmarks, which are challenging to implement for highly mobile gastrointestinal (GI) organs. To address this, patient-specific digital twins (DTs) modeling temporally varying motion were created to assess the accuracy of DIR methods. Approach. A total of 21 motion phases simulating digestive GI motion as 4D image sequences were generated from static 3D patient scans using published analytical GI motion models through a multi-step semi-automated pipeline. Eleven datasets, including six T2-weighted FSE MRI (T2w MRI), two T1-weighted 4D golden-angle stack-of-stars, and three contrast-enhanced computed tomography scans were analyzed. The motion amplitudes of the DTs were assessed against real patient stomach motion amplitudes extracted from independent 4D MRI datasets using hierarchical motion reconstruction. The patient-specific DTs were then used to assess six different DIR methods using target registration error, Dice similarity coefficient (DSC), and the 95th percentile Hausdorff distance using summary metrics and voxel-level granular visualizations. Finally, for a subset of T2w MRI scans collected from patients treated with magnetic resonance-guided radiation therapy, dose distributions were warped and accumulated to assess dose warping errors (DWEs), including evaluations of DIR performance in both low- and high-dose regions for patient-specific error estimation. Main results. Our proposed pipeline synthesized patient-specific DTs modeling realistic GI motion, achieving mean and maximum motion amplitudes and a mean log Jacobian determinant within 0.8 mm and 0.01, respectively, similar to published real-patient gastric motion data. It also enables the extraction of detailed quantitative DIR performance metrics and supports rigorous validation of dose mapping accuracy prior to clinical implementation. Significance. The developed pipeline enables rigorously testing DIR tools for dynamic, anatomically complex regions facilitating granular spatial and dosimetric accuracies.

Figures

Overview of the pipeline used to generate a patient-specific DT from a 3D abdominal scan. Step 1: Extract the medial axis by skeletonizing and pruning the AI-generated segmentation masks. Step 2: Generate a NURBS surface based on the medial axis. Step 3: Apply peristaltic motion to the target organs (stomach and large bowel), resulting in 21 phases representing different contraction states. Step 4: Compute the DVFs using the original and deformed NURBS surfaces.

Pseudocode illustrating the implementation of the motion model applied to the NURBS surfaces of the stomach and large bowel, as detailed in the methods section part (e).

Top plots (a),(b),(c): comparison of mean and maximum displacement magnitudes (±standard deviation) between ground truth and DT-synthetic stomach deformations across 21 phases in three distinct 4D-MRI datasets. Bottom plots (d),(e),(f): corresponding comparison of the mean log-Jacobian determinant (±standard deviation), reflecting local volumetric changes during deformation. Together, these plots illustrate both the extent (top) and anatomical plausibility (bottom) of the predicted motion fields relative to ground truth.

Example deformed segmentation masks and corresponding DVFs shown for each of the four imaging modalities: CECT, T1w gaSOS, T2w MRI for 3 different time points. Additionally, for the T2w MRI we show a sequence of snapshots of a stomach NURBS surface over four motion phases, illustrating the progression of the applied wave-like motion traversing the organ. The sectional views highlight how the modeled peristaltic deformation propagates along the surface.

For five T2w MRI datasets, we report segmentation performance (DSC and HD95), registration accuracy (TRE), and the DWE across various deep learning and variational DIR methods. Motion was applied separately to the stomach, duodenum, and large bowel in each dataset.

Patient-Specific granular performance of the different DIR methods on patient A where the stomach, duodenum and large bowel have a small overlap with the high radiation zone. (a) Global mean RMSE (mm) binned by motion magnitude (b) stomach, duodenum and large bowel mean RMSE (mm) binned by Gy radiation level. (c) Patient scan visualization with the RMSE overlapped. (d) Patient scan visualization with the dose map overlapped.

Patient-specific granular performance of the different DIR methods on patient B where the stomach, duodenum and large bowel have a big overlap with the high radiation zone. (a) Global mean RMSE (mm) binned by motion magnitude (b) stomach, duodenum and large bowel mean RMSE (mm) binned by Gy radiation level. (c) Patient scan visualization with the RMSE overlapped. (d) Patient scan visualization with the dose map overlapped.

Tables

Table 1.. Details of the image volumes and voxel dimensions used for motion synthesis across 11 patient datasets.

Patient	MR sequence details	Volume size (voxels)	Voxel size (mm)
1	T1w gaSOS	224 × 224 × 96	2.188 × 2.188 × 3.5
2	T1w gaSOS	224 × 224 × 96	2.188 × 2.188 × 3.5
3	CECT	512 × 512 × 201	0.9766 × 0.9766 × 2
4	CECT	512 × 512 × 201	0.9766 × 0.9766 × 2
5	CECT	512 × 512 × 201	0.9766 × 0.9766 × 2
6	T2w MRI	512 × 512 × 50	0.7813 × 0.7813 × 5
7	T2w MRI	448 × 448 × 125	1 × 1 × 2
8	T2w MRI	448 × 448 × 125	1 × 1 × 2
9	T2w MRI	448 × 448 × 125	1 × 1 × 2
10	T2w MRI	448 × 448 × 125	1 × 1 × 2
11	T2w MRI	448 × 448 × 125	1 × 1 × 2

Table 2.. Maximum motion magnitudes for 2 gastric organs (stomach and large bowel) computed from the DVFs generated across the 21 phases of motion for all 11 patients with synthesized motion. Patients 1–2 lacked large bowel data due to a limited field of view, and patients 3–6 were excluded because of inaccurate large bowel segmentations that prevented reliable NURBS surface generation.

		Motion magnitude (mm)
Patient	Full body	Stomach	Large bowel
1	13.2	13.2	—
2	14.2	13.9	—
3	13.2	13.2	—
4	13.5	13.1	—
5	13.5	13.4	—
6	14.3	14.3	—
7	8.7	8.7	8.6
8	8.6	6.6	8.6
9	8.6	8.6	8.6
10	8.7	8.7	7.7
11	8.6	8.6	8.6

Table 3.. Means and standard deviations of motion magnitudes for 2 gastric organs (stomach and large bowel) computed from the DVFs generated across the 21 phases of motion for all 11 patients with synthesized motion. Patients 1–2 lacked large bowel data due to a limited field of view, and patients 3–6 were excluded because of inaccurate large bowel segmentations that prevented reliable NURBS surface generation.

	Motion magnitude (mm)
Patient	Full body	Stomach	Large bowel
1	0.04 ± 0.01	3.05 ± 0.83	—
2	0.03 ± 0.01	2.85 ± 0.9	—
3	0.04 ± 0.02	2.22 ± 0.72	—
4	0.01 ± 0.01	2.32 ± 0.65	—
5	0.06 ± 0.02	3.27 ± 0.96	—
6	0.04 ± 0.02	3.24 ± 0.82	—
7	0.06 ± 0.00	0.92 ± 0.03	0.77 ± 0.03
8	0.08 ± 0.00	2.57 ± 0.01	2.74 ± 0.04
9	0.08 ± 0.00	3.65 ± 0.05	2.96 ± 0.01
10	0.07 ± 0.00	3.57 ± 0.03	2.78 ± 0.05
11	0.09 ± 0.00	3.69 ± 0.01	2.96 ± 0.03

Table 4.. Segmentation accuracy of various variational registration methods applied to T1WI gaSOS, and CECT. Stomach & Duo: duodenum was included in the segmentation mask of the stomach.

		T1w gaSOS		CECT
		Patient 1	Patient 2	Patient 3	Patient 4	Patient 5
Metric	DIR Algorithm	Stomach & Duo	Stomach & Duo	Stomach	Stomach	Stomach
TRE ↓ (mm)	HSOF	2.14 ± 1.31	2.70 ± 1.94	3.06 ± 1.62	3.33 ± 1.93	3.11 ± 2.01
	EVO	3.02 ± 1.70	3.64 ± 2.12	3.88 ± 2.30	3.33 ± 1.90	3.44 ± 2.27
	Elastix	3.40 ± 2.19	2.58 ± 1.45	2.48 ± 1.33	3.48 ± 2.12	3.48 ± 2.41
	Demons	3.26 ± 1.93	3.08 ± 1.94	3.27 ± 1.68	3.60 ± 2.00	3.46 ± 2.05
DSC ↑	HSOF	0.81	0.78	0.99	0.79	0.77
	EVO	0.77	0.81	0.98	0.82	0.81
	Elastix	0.68	0.77	0.98	0.76	0.71
	Demons	0.74	0.8	0.99	0.79	0.77
HD95 ↓ (mm)	HSOF	6.00	7.00	1.00	7.48	8.31
	EVO	8.06	8.77	1.73	6.71	8.06
	Elastix	8.06	7.35	2.24	7.18	8.66
	Demons	8.06	6.16	1.41	7.34	8.83

Full text

Introduction

Deformable image registration (DIR) is a critical component of both computed tomography (CT) and magnetic resonance (MR) image-guided adaptive radiotherapy (Rigaud et al 2019, Lowther et al 2022, Meyer et al 2025). It enables automated propagation of organ-at-risk (OAR) segmentations across treatment fractions and facilitates accurate estimation of the radiation dose delivered to both OARs and target structures throughout the treatment course. A key obstacle towards clinical implementation of DIR, particularly for luminal gastrointestinal (GI) organs where motion can range from a few millimeters up to 3 cm (Mostafaei et al 2018, Liu et al 2021, Zhang et al 2021), is the lack of methods to assess DIR accuracy at a patient-specific and granular level.

Target registration error (TRE), which measures the displacement errors of known key points as outlined in TG-132 (Brock et al 2017), provides a granular and unbiased assessment of DIR accuracy. In contrast, segmentation overlap metrics are influenced by the size of evaluated structures, potentially limiting their sensitivity and interpretability. However, TRE requires manually or semi-automatically placed visually discernible keypoints, which are difficult to identify reliably in GI organs due to their large and often unpredictable motion across and within radiation treatment fractions. Keypoints are typically placed in high-contrast regions with large intensity gradients, areas where DIR algorithms tend to perform well, potentially biasing the evaluation. Manually placed keypoints may not reflect registration accuracy in low-contrast soft-tissue regions, leading to an overestimation of DIR performance. Previous studies have highlighted the importance of using digital phantoms to address this limitation and evaluate DIR methods for dose accumulation (Papachristou et al 2024) and even generated ground truths to assess the performance of DIR methods for prostate cancer (Katsoulakis et al 2024). The goal of this study is to develop a patient-specific digital twin (DT) (Bosma et al 2021, 2024) framework to serve as known ground truth for evaluating the accuracy of DIR and accumulated dose estimates in the GI tract.

Population-level digital phantoms simulating respiratory and GI motions have been previously developed (Segars et al 2013, Subashi et al 2023). Prior work by Subashi et al implemented various GI motion patterns including peristalsis, rhythmic segmentation in the stomach, small and large bowel, high amplitude propagating contractions (HAPCs) in the large bowel, and tonic contractions in the GI sphincters using published analytical models within an XCAT computational phantom (2023). Prior studies have also focused on extracting individual motions such as breathing versus stomach contractile motions to extract the amplitude and time scale of such motions from real patients (Johansson et al 2021, Zhang et al 2021). However, while population-based anatomical models are an important tool for advancing research in medical imaging and radiation therapy (Segars et al 1999, Caon 2004, Johnson et al 2009, Cassola et al 2011, Ding et al 2012, Bosca and Jackson 2016), they fall short in capturing many of the complexities of individual patient anatomy.

First, population-level models such as XCAT can generate a wide range of physiologically realistic motions; however, these motions are simulated on generic male and female anatomies. In contrast, patients exhibit nuanced anatomical variations that cannot be fully captured by population-level phantoms. Second, XCAT-based MR phantoms rely on predefined signal intensity values to represent body composition, which cannot be used to evaluate deep learning (DL) DIR methods, as they are sensitive to intensity distribution shifts from training patient scans.

To address the aforementioned limitations of XCAT phantoms, we developed a framework to generate patient-specific DTs simulating temporally varying GI motion. Our DT framework aligns with the definition of DTs set forth by the ecosystem DTs in health (EDITH) Coordination and Support Action funded by the European Commission (EDITH, 2022). EDITH defines DTs for health as computer simulations, incorporating both knowledge-driven and data-driven models to predict clinically relevant quantities otherwise challenging to experimentally measure. In accordance with EDITH, our framework predicts feasible anatomic variations starting from 3D scans for MR guided adaptive treatment planning, where acquiring 4D MRIs are not clinically performed during treatment. In addition, our models used published amplitude and time-scale data derived from real patient imaging (Zhang et al 2021) and apply this motion directly to each patient’s scan.

Our pipeline starts from individual patient imaging data, extracting the various organs using a semi-automated manner, and creating a digitized representation as non-uniform rational B-spline (NURBS) surfaces. A motion pattern is then applied to this model to generate a patient-specific sequence of deformations simulating organ motion. We validate these simulated motion patterns by comparing them to real patient motion data from the literature.

Our key contributions include: (a) development of a patient-specific DT framework that simulates temporally varying GI motion and provides known ground truth for DIR evaluation, (b) demonstration of this framework across multiple imaging modality, including contrast-enhanced CT (CECT), T2-weighted fast spin echo MRI (T2w MRI), and T1-weighted 4D golden-angle stack-of-stars MRI (T1w gaSOS), (c) the first framework, to our knowledge, enabling granular, patient-specific assessment of registration and dosimetric accuracy in the GI tract, and (d) a modality and algorithm independent evaluation framework, demonstrated on both iterative and DL-based DIR methods.

We envision our approach to be useful for performing day-to-day and patient-specific assessments of a preferred DIR method used in a clinic for the MR-Linac (Fallone 2014, Keall et al 2014, Lagendijk et al 2014, Mutic and Dempsey 2014) enabled adaptive radiation treatment or assessing quality of treatment planning simulations.

Materials and methods

DT pipeline for GI motion simulation

2.1.

We developed a semi-automated pipeline to generate patient-specific DTs that simulate temporally varying GI motion. The pipeline produces 4D image sequences by deforming static 3D patient scans using analytical GI motion models (Subashi et al 2023) applied to organ-specific surface representations. Twenty-one motion phases were synthesized to match the temporal resolution of stomach contractions reported by Johanssen et al (2021). The pipeline comprises the following steps (figure 1):

• (a)Organ segmentation: A published AI model called the self-distilled masked image transformer (Jiang et al 2022), was used to segment the liver, stomach, duodenum, spleen, kidneys, small and large bowel from a given patient image set (MR or CT). Segmentations were manually verified by an expert (figure 1. Step 1, row 1).
• (b)Organ-specific skeleton graph extraction: Morphological thinning of binary 3D segmentation masks is performed to obtain a minimally connected sequence of points by making successive passes of the image and removing the identified border pixel, while preventing a break in connectivity of points spanning the organ (Lee et al 1994). This skeleton may contain branches which require further processing (figure 1. Step 1, row 2).
• (c)Organ centerline extraction using the longest medial axis: Starting and ending points are manually selected from each extracted skeleton. From the starting point we build a graph by connecting all directly neighboring points until all points have been explored. Then we apply a breadth-first search algorithm (Lee 1961) to determine the longest connected path graph between the two manually selected points. This path is retained as the organ centerline while all branching structures are discarded (figure 1. Step 1, row 3).
• (d)NURBS surface extraction: For each organ, we sample points along the extracted centerline to create a set of equidistant center points. For each of these points we cast radial vectors perpendicular to the direction of the centerline. These vectors are expanded outward from the centerpoint until they intersect with the boundary of the mask surface. These boundary points on the mask surface served as control points to define the NURBS surface to define the organ’s geometry. Each cross-section along the centerline defined a sectional curve (figure 1. Step 2) For the surface. These NURBS surfaces were then used to synthesize GI motion by applying a motion model at the surface points.
• (e)Multi-phase motion synthesis: Analytical motion models developed by Subashi et al (2023) were applied to the NURBS surfaces to simulate realistic GI motion patterns using traveling sinusoidal waves. These include HAPCs, tonic contractions, and peristalsis (figure 1. Step 3).The overall modeling framework is given by the following equation:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{equation*}P_{_{i,j}}^{\prime} = {P_{i,j}} + F\left( {\overrightarrow {{x_i}} ,t} \right) \cdot {D_s} \cdot {D_t} \cdot \vec d\end{equation*}\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} {P_{i,j}}\end{document} represents the control point \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} j\end{document} at the sectional curve \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} i\end{document} of a given NURBS surface. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} F\left( {\vec i,t} \right)\end{document} is a function modeling the non-dispersive component of the wave, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} {D_s}\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} {D_t}\end{document} represent functions modeling dispersions of the wave in the spatial coordinate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \left( {\overrightarrow {{x_i}} } \right)\end{document} and temporal \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \left( t \right)\end{document} domain respectively and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \vec d\end{document} represents a directional vector to apply a consistent expansion/contraction in the radial direction from the center for each sectional curve.Using this framework, we iterate through each sectional curve of the given NURBS surface and compute the wave and dispersion magnitudes as functions of space and time. In this work, we specifically focused on modeling the peristaltic motion \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ({F_{{\text{PS}}}}\left( {\overrightarrow {{x_i}} ,t} \right)\end{document} ), within the stomach and large bowel. The wave function was modulated via \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} {F_{{\text{PS}}}}\left( {\overrightarrow {{x_i}} ,t} \right)\end{document} and an exponential dispersion function \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} {D_z}\left( u \right)\end{document} . Of note, the DT framework is flexible and can be adapted to incorporate any GI motion patterns and dispersion functions. The peristaltic wave and exponential dispersion function are given by:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{equation*}{F_{{\text{PS}}}}\left( {\overrightarrow {{x_i}} ,t} \right) = { }\frac{1}{{\sqrt 3 }} \cdot A \cdot \sin \left( {2\pi \frac{{{L_i} - c \cdot t}}{\lambda }} \right)\end{equation*}\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{equation*}{D_z}\left( u \right) = {{\text{e}}^{ - \alpha \cdot u}}\end{equation*}\end{document}

where A, c, λ, and α are user-defined parameters representing the amplitude, speed, and wavelength of the peristaltic wave, and the attenuation of amplitude in the spatial or temporal domain, respectively. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} {L_i}\end{document} is the position of the sectional curve along the organ’s length, and u is a generalized variable that corresponds to either \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \overrightarrow {{x_i}} \end{document} or \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} t\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} z\end{document} is an index identifying the spatial \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \left( s \right)\end{document} or temporal \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \left( t \right)\end{document} domain.As illustrated in figure 2, for each motion phase we iterate through each sectional curve (figure 2. Line 3) and each control point in the NURBS (figure 2. Line 8) calculating and applying the deformation to each control point (figure 2. Line 9–15). Each full iteration over the centerline of the NURBS generates one phase. Applying the motion models to the NURBS surfaces yields a sequence of deformed surfaces, which serve as the basis for computing the corresponding deformation vector fields (DVFs).

• (f)(f) Extraction of DVFs: From the sequence of deformed NURBS surfaces, we can derive a 3D DVF that captures the deformation for each NURBS by:

• Computing inner-shell vectors using radial interpolation. To capture internal surface deformation within each GI organ, intermediate ‘inner shell’ surfaces were generated between the NURBS outer surface and the organ’s centerline. These shells were parameterized by a radial factor \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} p \in \left[ {0,1} \right]\end{document} , with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} p = 0{\text{ }}\end{document} at the outer surface and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} p = 1\end{document} at the centerline. For each shell defined by p_k_, a new NURBS surface was computed by interpolating between the outer and centerline surfaces. Sampling each shell over a uniform \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \left( {u,{\text{ }}v} \right)\end{document} grid yielded corresponding 3D points in both the original and deformed configurations. These parameters are standard in NURBS surface definitions, where u and v define positions within the 2D parameter space of the surface.

• Deformation vectors: Deformation vectors were then calculated pointwise as the difference between the deformed and original positions: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} {V_{{\text{surf}}}}\left( {{u_i}{ }{v_j},{p_k}} \right) = X{^{^{\prime}}}\left( {{u_i}{ }{v_j},{p_k}} \right) - X\left( {{u_i}{ }{v_j},{p_k}} \right)\end{document}

• Voxelization and smoothing of the DVF: The computed deformation vectors across all \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \left( {u,{\text{ }}v,{\text{ }}p} \right)\end{document} samples were mapped into a fixed 3D voxel grid by assigning each vector to a voxel based on its original undeformed position. When multiple vectors mapped into the same voxel, they were averaged. To improve continuity, a spatial smoothing step was applied wherein zero-motion voxels were iteratively updated by averaging the vectors of neighboring non-zero voxels until the average change in each voxel was less than 0.001 mm. This will generate a DVF that can be applied to the patient’s static scan to deform it in a way that reflects both surface and internal anatomical displacements (figure 1, Step 4).

• (g)Creation of 4D MRI/ 4D CT sequences: By repeating the previous step across the full sequence of deformed NURBS surfaces, we synthesize patient-specific motion that is agnostic to any imaging modality. Applying the corresponding DVFs to the original patient scans produces a 4D sequence of GI motion while preserving the internal contents of the luminal organs (e.g. air, fluid). Sample videos of synthetized motion sequences are included in supplemental material for representative cases.

Overview of the pipeline used to generate a patient-specific DT from a 3D abdominal scan. Step 1: Extract the medial axis by skeletonizing and pruning the AI-generated segmentation masks. Step 2: Generate a NURBS surface based on the medial axis. Step 3: Apply peristaltic motion to the target organs (stomach and large bowel), resulting in 21 phases representing different contraction states. Step 4: Compute the DVFs using the original and deformed NURBS surfaces.

Pseudocode illustrating the implementation of the motion model applied to the NURBS surfaces of the stomach and large bowel, as detailed in the methods section part (e).

For additional details on steps (a) through (d), please refer to the pseudocode provided in the appendix document. Additionally, the code is available upon request to the authors.

Datasets and motion parameters

2.2.

Under institutional review board–approved protocols (IRB 21-129 and HUM00260112), all retrospective patient data were analyzed in accordance with ethical standards. In total, 11 datasets were analyzed, comprising six T2w MRI, two T1w gaSOS, and three CECT scans. Additionally, three independent T1w gaSOS datasets were used solely for quality assurance, with motion models tailored to the motion patterns observed in those specific cases. To the stomach, we applied a traveling wave equation with parameters A = 16 mm, λ = 55 mm, and c = 5 mm s^−1^. In addition, for the five datasets that included the large bowel, a second traveling wave equation A = 16 mm, λ = 40 mm, c = 8 mm s^−1^ was applied to simulate large bowel motion. Large bowel motion was not analyzed for a subset of patients either because of a smaller field of view, as in the cases of patients with liver cancer (patients 1 and 2), or due to insufficient scan lengths that did not capture the full extent of the large bowel (patients 3–6), which prevented the reliable generation of corresponding NURBS surfaces. For all synthesized motions used in the DIR evaluation, α = 0 was chosen to emulate high-magnitude, non-dispersive GI motion. The amplitude and wavelength parameters were chosen to align with the peristaltic motion measured from 10 patients treated on a 1.5 T MR-linac and reported by Subashi et al (2023).

Patient-specific motion quality assurance

2.3.

To evaluate motion realism, we compared synthetic stomach motion from the DT framework to motion derived from three T1w gaSOS MRI datasets with known stomach contractions. Specifically, we used the previously published known stomach deformation motions, which includes detailed recordings of stomach contractions over 21 distinct phases. Previous studies have used hierarchical MRI reconstructions to model stomach contractions and slow drifts with and without respiratory motion (Zhang et al 2021). Employing the ground truth DVFs derived from three scans with hierarchical motion models as baseline references, we fit traveling sinusoidal waves tailored specifically to each contraction phase. This allowed the DT pipeline to generate synthetic stomach motions closely aligned with the experimentally measured motions. To validate the synthesized motion, we computed the mean and maximum displacements of the ground-truth DVFs and the log of the Jacobian which captures the degree of local non-rigid, non-affine deformations.

Motion analysis

2.4.

For all 21 synthetically generated motion phases across the 11 patient scans, the maximum and mean displacement magnitudes of the DVFs generated by the pipeline were computed for the full body as well as two gastric organs: the stomach and large bowel.

Evaluation of DIR methods using DTs

2.5.

The goal of the experiment was to assess the feasibility of using the DTs to evaluate a range of DIR methods using the synthesized GI motions produced with various imaging modalities for multiple patients. Of note, the goal was not to compare individual methods with respect to one another; all methods were used as is without performing any hyperparameter or other optimization for the evaluated DT datasets.

Six different DIR algorithms were assessed including 4 variational and 2 DL methods. The variational methods included Horn–Schunck optical flow (HSOF), EVolution multimodality (EVO), Elastix with mutual information, and the Iterative Demons algorithm, all of which used iterative optimization using image intensities. Key differences between the methods include the use of brightness consistency and smoothness regularization for motion fields used in HSOF (Horn and Schunck 1981), normalized intensity gradients with smoothness regularization used in EVO (Denis de Sennevill et al 2016), B-spline parameterization of the transformation used in Elastix (Klein et al 2010, Leibfarth et al 2013) and optical flow with Gaussian smoothing in diffeomorphic Demons (Thirion 1998, Vercauteren et al 2009). The two DL DIR methods were VoxelMorph (Balakrishnan et al 2019) and an enhancement of VoxelMorph called progressively refined registration-segmentation (ProRSeg) that uses convolutional long short-term memory networks in the encoder (Jiang et al 2023). Testing datasets used for evaluation were never used for training VoxelMorph and ProRSeg to avoid data leaks. Both VoxelMorph and ProRSeg were trained on different sets of real patient T2w datasets used in Jiang et al (2023) and were only applied to DTs corresponding to the same MR scanning sequence.

Evaluation metrics

2.6.

Registration was performed between the original scan and the phase exhibiting the maximum deformation. Geometric accuracy was computed by measuring organ segmentation accuracy using the DSC and Hausdorff distance at 95th percentile as well as TRE using the surface and inner shell points of a given NURBS surface. Differences in dose accumulation across methods were assessed by simulating dose accumulation using both the ground truth DVF and the DVFs produced by each registration method. Dose deformation was performed using direct dose mapping. DWE comparing the accumulated dose distributions was calculated as:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{equation*}{\text{DWE}} = { }\mathop \sum \limits_{i = 1}^N \frac{{{\text{Accum }}{D_{i,{\text{DIR}}}} - {\text{Accum }}{D_{i,{\text{GT}}}}}}{{{\text{Accum }}{D_{i,{\text{GT}}}}}}.\end{equation*}\end{document}

Granular assessment of dosimetric and displacement errors were computed at the anatomic voxel level and quantified using root mean squared error (RMSE), which then were also visualized as heatmaps for a more detailed evaluation of model accuracy in regions subject to large deformations or high-radiation dose exposure.

Results

Our semi-automated pipeline was applied to generate motion for 11 different patients from various anatomic imaging modalities. The framework takes 5–15 min to synthesize peristaltic motion from a 3D input image, given fixed motion parameters, NURBS surfaces, with the synthesis time varying due to the number of set motion phases. The framework allows to interactively change the motion parameter configurations as well as NURBS surfaces to also visually match known peristaltic patterns. The framework was capable of generating motion for image volumes and voxels with varying sizes, as shown in table 1. The method was applicable to axial and coronal image reconstructions.

Assessing patient-specific motion quality assurance

3.1.

The top row of figure 3 shows a comparison between DT-generated and ‘ground truth’ stomach motions derived from Zhang et al using hierarchical motion modeling (2021), showing both mean and maximum displacement magnitudes. Across all motion phases, the DT-generated motion exhibited mean and maximum displacement differences within 0.8 mm of ground truth. Furthermore, as shown in the bottom row of figure 3, the DT-simulated deformations yielded consistent log-Jacobian means across all motion phases, remaining within 0.01 of the ground truth mean over all 21 phases. This indicates strong agreement in the overall deformation characteristics.

Top plots (a),(b),(c): comparison of mean and maximum displacement magnitudes (±standard deviation) between ground truth and DT-synthetic stomach deformations across 21 phases in three distinct 4D-MRI datasets. Bottom plots (d),(e),(f): corresponding comparison of the mean log-Jacobian determinant (±standard deviation), reflecting local volumetric changes during deformation. Together, these plots illustrate both the extent (top) and anatomical plausibility (bottom) of the predicted motion fields relative to ground truth.

Modality agnostic temporally varying motion synthesis

3.2.

Figure 4 depicts the synthesized motions with the corresponding DVFs for all three analyzed imaging modalities, CECT, T1w gaSOS, and T2w MRI. The deformed masks and 2D image slices taken at different time points along the traveling wave sequence propagating through the stomach are also shown. The maximum stomach motion shown in table 2 ranged from 8.56 mm to 14.34 mm, while the large bowel exhibited motion magnitudes between 7.69 mm and 8.64 mm, indicating the ability to simulate motion at varying amplitudes. The mean and standard deviation in the motion magnitudes for the gastric organs and the whole body are summarized in table 3, which showed stomach mean displacements ranging from 0.92 mm to 3.69 mm, and the large bowel from 0.77 mm to 2.96 mm.

Example deformed segmentation masks and corresponding DVFs shown for each of the four imaging modalities: CECT, T1w gaSOS, T2w MRI for 3 different time points. Additionally, for the T2w MRI we show a sequence of snapshots of a stomach NURBS surface over four motion phases, illustrating the progression of the applied wave-like motion traversing the organ. The sectional views highlight how the modeled peristaltic deformation propagates along the surface.

Evaluating DIR performance using the DT framework

3.3.

DIR was computed using various methods between the phase 0 (static input patient 3D scan) and the phase with the largest amplitude displacement with respect to phase 0. Variational DIR methods were applied to the CECT, T1w gaSOS, and T2w MRI, whereas the DL-DIRs trained with T2w MRI were only applied to T2w MRI datasets.

Table 4 shows the accuracy metrics computed for the variational methods applied to the CECT and T1w gaSOS scans. The DIR methods were similarly accurate for T1w gaSOS compared to CECT using TRE (2.14 ± 1.31 mm to 3.64 ± 2.12 mm versus 2.48 ± 1.33 mm to 3.60 ± 2.00 mm), HD95 (6.00 mm to 8.77 mm versus 1 mm to 8.83 mm), and DSC (0.68–0.81 versus 0.71–0.99). HD95 and DSC were computed for two CECT patients and three T1w gaSOS cases; however, standard deviation values are not reported due to the limited sample size.

In the case of T2w MRI from 5 patients with pancreatic cancer, the DIR methods showed slightly higher TRE 4.14 ± 2.23 mm to 5.14 ± 2.51 mm for the stomach and duodenum compared to T1WI gaSOS and CECT images, a higher DSC from 0.84 ± 0.02**–**0.92 ± 0.02, and lower HD95 ranging from 2.71 ± 0.51 mm to 4.4 ± 0.18 mm. The registration accuracies for the large bowel were slightly lower than the stomach and duodenum with TRE ranging from 3.06 ± 0.48 mm to 4.35 ± 0.29 mm, DSC between 0.89 ± 0.02–0.95 ± 0.01, and HD95 between 1.94 ± 0.65 mm and 3.82 ± 0.44 mm, respectively.

In addition, the same 5 patients also had radiation treatment dose maps, which were used to calculate the DWE for the same three organs. As shown in figure 5, the mean DWE ranged from 6.68 ± 1.45% to 9.80 ± 1.83% for stomach and duodenum and a smaller error of 3.88 ± 1.00% to 6.72 ± 1.81% for the large bowel. The results also show variation in the performance of the different methods across the different datasets and modalities, thus providing an approach to assess the relative merits of the various methods for segmentation, registration, and dose warping on an individual patient level.

For five T2w MRI datasets, we report segmentation performance (DSC and HD95), registration accuracy (TRE), and the DWE across various deep learning and variational DIR methods. Motion was applied separately to the stomach, duodenum, and large bowel in each dataset.

Patient-specific granular errors

3.4.

Capability to assess granular (or voxel-level) error visualization on a patient level is demonstrated for two representative patients, patient A (figure 6) where organs undergoing motion occurred in the low-radiation dose regions (10–30 Gy) (see figure 6(d)) and patient B (figure 7) with the same organs located in the high radiation dose regions (exceeding 30 Gy) (see figure 7(d)). Global RMSE is visualized with respect to increasing motion magnitudes to assess errors as a function of motion (figures 6(a) and 7(a)). Motion magnitude was binned at fixed intervals ranging from a minimum of 0 mm to a maximum of 8.65 mm. Patient specific analysis showed that the mean RMSE ranged from 0.7 mm in low-motion regions (0–1 mm) to 5 mm in high-motion regions (>8 mm) for the two patients. RMSE was also computed within the GI organs undergoing motion to assess impact of motion on the accuracy with respect to the radiation dose delivered to the organ, which showed a range from 0.10 mm to 2.00 mm, with the highest errors occurring mostly in low-dose regions for Patient A, but a higher error of 1.00 mm to 2.5 mm in the high-dose region for Patient B (figures 6(b) and 7(b)).

Patient-Specific granular performance of the different DIR methods on patient A where the stomach, duodenum and large bowel have a small overlap with the high radiation zone. (a) Global mean RMSE (mm) binned by motion magnitude (b) stomach, duodenum and large bowel mean RMSE (mm) binned by Gy radiation level. (c) Patient scan visualization with the RMSE overlapped. (d) Patient scan visualization with the dose map overlapped.

Patient-specific granular performance of the different DIR methods on patient B where the stomach, duodenum and large bowel have a big overlap with the high radiation zone. (a) Global mean RMSE (mm) binned by motion magnitude (b) stomach, duodenum and large bowel mean RMSE (mm) binned by Gy radiation level. (c) Patient scan visualization with the RMSE overlapped. (d) Patient scan visualization with the dose map overlapped.

Figures 6(c), 7(c), 6(d) and 7(d) show a visualization of voxel-wise RMSE within the organs undergoing motion and the dose maps, respectively to provide a visual representation of the errors for individual patients.

Discussion

In this work, we extended the concept of population-level digital phantoms to develop patient-specific DTs that model temporally varying GI motion. Our semi-automated pipeline starts from AI automated organ segmentations, which then are used to generate peristaltic motion of varying amplitudes and time scales for luminal organs such as the stomach, duodenum and the large bowel. Our analysis indicates that the synthesized stomach motions can closely match known reference values, both in displacement magnitude and deformation characteristics producing motion magnitudes consistent with those reported by Zhang et al (2021). Thus, the framework showed initial feasibility to generate realistic motion for three different GI organs such as stomach, duodenum, and large bowel, with differing geometries across multiple common radiological imaging modalities,

Our framework is applicable to multiple anatomic imaging modalities and demonstrated feasibility to evaluate multiple variational and two different DL DIR methods. In addition to evaluating registration, our framework can be easily extended to evaluate dose warping accuracy summarizing errors for individual organs as well as on a granular level to assess accuracy variations on voxel-level. As a result, our DT framework enables voxel-wise visualization of registration errors, facilitates analysis of error patterns across motion regimes (e.g. low vs high motion), and supports individualized assessment of how registration inaccuracies can affect radiation dose. To our knowledge, this is the first comprehensive simulation and evaluation of DIR using patient-specific DTs across multiple imaging modalities for GI luminal organs.

Whereas, previous efforts focused on modeling respiratory and cardiac motion, ours focused on modeling the digestive motion (Segars et al 1999, 2013, 2018, 2019). One prior work by Subashi et al demonstrated the ability to generate MR-like digital phantoms incorporating a range of GI motion types including peristalsis, slow and fast gastric contractions, and HAPCs (Subashi et al 2023). However, all aforementioned prior works synthesized digital phantoms modeling population-level anatomy built on the generic adult male and female XCAT models. A limitation with modeling population-level motion is that it does not represent individual patient anatomy variations. Our work, for the first time, addresses the key issue of modeling patient-specific variations by creating patient-specific DTs of gastric motion.

Another limitation of population level modeling using XCAT requires simulation of MR images using fixed signal intensities for each organ, that can create a domain shift for assessing DL registration methods. In our work, synthesis starts from the original MRI, the synthesized motions are also created on MRI, hence allowing to evaluate DL DIR methods.

Finally, GI organs undergo substantial and arbitrary motion that varies from patient to patient despite common motion mitigation strategies such as pneumatic compression belts (Mostafaei et al 2018, Liu et al 2021, Zhang et al 2021, Alam et al 2022). Our approach allows to vary and create a variety of motion amplitudes and rigorously evaluate DIR methods under various GI motion amplitudes.

Limitations of our current framework include the lack of support for modeling respiratory motion as the focus of this work was isolated GI motion. We also excluded small bowel motion simulation because it is generally segmented in the clinic as a ‘bowel-bag’, making the extraction of NURBS to model motion along the tubular region difficult. Nevertheless, our approach can potentially be extended to additional organs and other types of motions involving such organs. Finally, GPU acceleration is not currently implemented, and manual steps such as parameter selection increase computation time. This design prioritizes flexibility, allowing the generation of motions with user-defined amplitudes. While the current pipeline is intended for offline analysis, incorporating GPU acceleration and predefined parameters could reduce processing time and enable near real-time patient-specific DT generation.

Conclusion

We developed a semi-automated DT pipeline to generate realistic GI temporally varying motion in the stomach, duodenum and large bowel from multiple anatomic imaging modalities. Our framework showed capability to generate motions within ranges seen in real patients, indicating feasibility to evaluate multiple DIR methods. Our framework enables evaluating dose warping and registration errors in a granular voxel-wise manner for individualized patient-level analysis, suitable for rigorous analysis required for clinical deployment.