Exploratory application of DMD for particle deposition and fluid field in the respiratory tract

TL;DR

This paper explores using DMD to efficiently simulate particle deposition and airflow in the respiratory tract, reducing computational costs while maintaining accuracy.

Contribution

The first demonstration of DMD as a reduced-order model for respiratory tract particle and fluid simulations.

Findings

DMD reconstructed fluid fields with a Mean Relative Error of 12% compared to CFD simulations.

DMD achieved 85% correlation with ground truth in approximating overall particle distribution.

Localized discrepancies in DMD reconstructions reached up to 70% relative error.

Abstract

Simulating particle deposition in the respiratory tract requires high computational effort due to the intricate airway geometry and complex airflow–particle interactions. To address this challenge, this study introduces the first demonstration of Dynamic Mode Decomposition (DMD) as a reduced-order model to infer the trajectories of inhaled particles during a breathing cycle and to evaluate the applicability of DMD as a fluid field interpolator. The periodic nature of respiration and the predominance of sinusoidal boundary conditions make it well-suited for DMD analysis. Three high-fidelity computational fluid dynamics (CFD) simulations were performed under three different inlet volume airflow conditions for the same realistic adult male anthropomorphic phantom respiratory tract model. Reduced-rank DMD reconstructions were compared to the CFD ground truth, yielding a Mean Relative Error (MRE) of 12% in the velocity field. Additionally, a fourth simulation was conducted at an intermediate point to evaluate the interpolation capability of the parametric DMD framework in complex systems. This interpolation resulted in an MRE of 20%, with the reconstructed flow field capturing dominant fluid modes and overall dynamics, though localized discrepancies reached relative errors up to 70%.

While DMD effectively reconstructed fluid fields, preserving mean flow regimes, some deviations were observed in Lagrangian particle tracking, specifically in spatial deposition resolution. However, the method approximated overall particle distribution with an 85% correlation to ground truth and was effective in representing regional deposition patterns across the tracheobronchial tree. These findings support the utility of DMD a computationally efficient approach for fluid field reconstruction and particle transport analysis in respiratory flow simulations.

Figures

Tables

Algorithm 1. Parametric DMD algorithm via reduced (discrete) Koopman operator interpolation (Huhn et al., 2023)

1.	Solve the fluid field for each of the NS breathing flow rates v˙jj=1NS and collect temporal snapshots.
2.	Divide each matrix of snapshots into a shifted matrix Xfluidv˙j+ and a lagged matrix Xfluidv˙j-.
3.	Perform SVD of Xfluidv˙j-:Xfluidv˙j-=UjΣjVjT, for 1≤j≤NS.
4.	Retain r modes and compute the reduced Koopman operator A~ for each of the training breathing flow rates (v˙j), see Eq. (12).
5.	Interpolate the Koopman operator A~θ using the respective operators A~fluidv˙1 and A~fluidv˙2 of the closest breathing flow rates (v˙1,v˙2) to the breathing flow rate of interest (v˙θ:v˙1<v˙θ<v˙2), see Eq. (19).
6.	Perform the eigen-decomposition of A~θ to obtain the reduced eigen-pair (Λθ,Wθ), see Eq. (21).
7.	Interpolate the POD modes Uθ using the respective U1 and U2 of the closest breathing flow rates (v˙1,v˙2) to the breathing flow rate of interest v˙θ:v˙1<v˙θ<v˙2, see Eq. (22).
8.	Construct the DMD modes Φθ=UθWθ, see Eq. (23).
9.	Interpolate the initial amplitudes b0θ using the respective b01 and b02 of the closest breathing flow rates (v˙1,v˙2) to the breathing flow rate of interest (v˙θ:v˙1<v˙θ<v˙2), see Eq.(24).
10.	Reconstruct the full-state solution Xreconstructedθ(t), see Eq. (18).

Equations

Full text

Introduction

Inhalation of aerosolized particles occurs across a range of environments, with elevated exposure risk in occupational settings involving dust, fumes, or hazardous aerosols. When inhaled material contains harmful substances, such as radioactive aerosols or toxic compounds, quantifying deposition patterns in the respiratory tract is necessary for evaluating internal exposure and guiding dosimetric analysis (Bartol et al., 2024). Obtaining high-resolution particle deposition profiles (PDPs) is also relevant in therapeutic applications involving inhaled drug delivery, as these profiles provide spatial information on impacted airway regions and, in the case of radioactive aerosols, inform dose estimation to radiosensitive tissues. Characterizing subject-specific particle distribution within the human airways can improve exposure assessment models, evaluate deposition variability across anatomical geometries, and support individualized simulation frameworks.

While computational fluid-particle dynamics (CFPD) has emerged in recent decades as a tool to assess subject-specific PDPs, as demonstrated by Katz et al. (1999), Kleinstreuer and Zhang (2010), Kleinstreuer et al. (2008), CFPD-based approaches remain a computationally intensive solution. Reduced-order models (ROMs) are an alternative tool for accelerating CFPD simulations by identifying and retaining only the dominant coherent structures in the flow. The concept behind ROMs is to reduce the number of degrees of freedom in the governing equations; this enables a rapid reconstruction of the flow field with minimal loss of fidelity. Multiple data-driven techniques operate by decomposing high-dimensional CFD data into a set of spatial modes that evolve with a set of fixed temporal characteristics. In the study by Habibi et al. (2020), a Dynamic Mode Decomposition (DMD)-based ROM was applied to simulations of pulsatile blood flow in a reduced geometry. By retaining 53 DMD modes, the reconstructed solution achieved approximately 99.5% agreement with the corresponding full-order computational fluid dynamics (CFD) results. This high level of agreement illustrates the capacity of DMD to capture the dominant temporal and spatial features in systems characterized by periodic boundary conditions.

In contrast, the study by Holemans et al. (2022) utilized Spectral Proper Orthogonal Decomposition (SPOD) to develop a ROM for large CFD data sets obtained from simulations of an annular swirling jet. By extracting data from two perpendicular 2D planes instead of processing the full 3D data set, the study achieved reduced computational time, achieving a 68-fold acceleration and reduced memory usage. However, this approach incurred a reconstruction error of ~10% for the relatively simple geometry, modeled with high resolution (> 10 million elements) and with a consistent revolving frequency of 22 Hz. Despite the higher error, the method effectively captured the large-scale coherent structures that govern the main flow dynamics. The ability to isolate and preserve these dominant features makes ROMs well-suited to respiratory flow modeling, where computational cost and geometric fidelity are both important for high-throughput or patient-specific simulation efforts.

In a related computational study, Jing et al. (2024) employed Large Eddy Simulation (LES) to analyze inspiratory airflow dynamics in a CT-derived human respiratory tract model under tidal and quasi-steady flow conditions, validated through grid independence tests. Using DMD on LES-derived vorticity fields, Jing et al. extracted coherent flow structures characterized by single frequencies and explicit decay/growth rates, which enabled an accurate identification of transient mechanisms. As a comparison, Proper Orthogonal Decomposition (POD) was applied, revealing its limitations in isolating frequency-specific dynamics due to energy-ranked modes that combine multiple spectral components. The study further demonstrated that DMD outperformed POD as the choice of ROM in respiratory applications, as its ability to resolve temporally evolving modes, such as turbulence in the pharynx-larynx region, was not captured when using POD. This advantage positions DMD as a more effective tool for dynamic analysis, including particle transport studies, where transient flow interactions dictate deposition outcomes.

To distinguish among the three decomposition methods, the work of Schmidt and Colonius (2020) provides a comparative analysis. POD extracts orthogonal spatial modes by performing an eigenvalue decomposition of the entire time-resolved dataset. These modes are ranked by energy and represent the most spatially coherent structures present throughout the simulation. Because POD does not retain temporal evolution, it produces a purely spatial decomposition that captures dominant features without time-dependent behavior.

In contrast, SPOD first applies a Fourier transform to the time-resolved data to isolate individual oscillation frequencies, followed by an eigenvalue decomposition at each frequency to obtain spatial modes. This frequency-domain approach enables SPOD to resolve structures that are coherent in both space and time, preserving the highest-energy modes at each frequency and linking them directly to the system’s temporal dynamics via the Fourier domain.

Finally, DMD seeks a linear operator that advances the system in time, using singular value decomposition (SVD) as part of the procedure. DMD identifies modes that are not constrained to be orthogonal but are optimized to reconstruct the system’s evolution, often capturing growth or decay rates and oscillatory behavior.

While POD and SPOD emphasize spatial compression using orthogonal basis functions, DMD prioritizes reconstruction of time-dependent system behavior (Towne et al., 2018). DMD modes may overlap spatially, but they are selected to best represent the system’s dynamic trajectory.

The present study investigates the reconstruction capacity of the DMD framework and evaluates its use for interpolating flow fields in highly complex and non-linear systems. Such interpolated fields can then be utilized under various particle deposition calculation schemes, such as the direct solution of uncoupled particle transport equations or more sophisticated data-driven methods.

DMD is a powerful data-driven method that extracts coherent structures and dynamics from high-dimensional systems (Kutz & Brunton, 2016). It is based on the Koopman operator theory, which describes the evolution of observables in a dynamical system. The Koopman operator, Eq. (1), $[eqn]$ , acts on a space of observables $[eqn]$ such that:

$$[eqn]$$

where $[eqn]$ represents the state of the system at time $[eqn]$ . DMD approximates the Koopman operator by decomposing the data into modes that evolve exponentially in time, allowing for the reconstruction and prediction of the system’s dynamics. Mathematically, given a sequence of snapshots of the system state $[eqn]$ , DMD, Eq. (2), seeks a linear operator $[eqn]$ that best approximates the dynamics:

$$[eqn]$$

The DMD modes and eigenvalues are then obtained from the eigen-decomposition of $[eqn]$ , providing a low-dimensional representation of the system’s behavior.

This study is structured into three phases. The first phase involved the application of DMD to reconstruct particle positions in the respiratory tract by directly applying DMD to the dynamics of particle positions. The second phase examined the reconstruction power of DMD for the fluid field in the respiratory tract and compared it with high-fidelity simulations. The third phase presents the results of interpolating the Koopman operator to reconstruct a fluid field for which the model was not explicitly trained. The interpolation scheme builds upon the work of Huhn et al. (2023), who demonstrated that interpolating the reduced Koopman operator relaxes smoothness constraints and improves fidelity relative to direct snapshot stacking. This three-phase structure provides a basis for evaluating the performance and limitations of DMD and Koopman-based surrogates relative to high-fidelity models.

Methodology

To reconstruct particle positions and fluid fields, three high-fidelity CFPD simulations were first performed on the same adult male human respiratory geometry, which was reconstructed from CT scans and implemented with two-way coupling between particles and fluid, meaning particles have mass and volume and that particles have feedback on the Navier–Stokes equations. A two-way coupled Lagrangian tracing was performed in StarCCM+ (Software, 2021), using a Lagrangian Multi-phase (LMP) solver.

CFD simulation details

2.1.

The airway model used represents the lower respiratory tract up to the 5th generation of bronchi. This geometry was reconstructed by segmenting the air cavities of a CT scan deck and then arranged sequentially using a pre-trained neural network from Wang et al. (2022) and post-processed for a CFPD simulation using the methodology adopted from Bartol et al. (2024). A schematic of the resulting geometry can be seen in Fig. 1:

The particle distribution used corresponded to a polydisperse log-normal distribution with mean $[eqn]$ and standard deviation $[eqn]$ equal to $[eqn]$ and $[eqn]$ , and a density of $[eqn]$ . This particle distribution relates to a form of aerosolized iodine-131, a radioisotope formed as a byproduct in nuclear fission relevant to nuclear accidents and fallout scenarios, which was previously used by Talaat et al. (2019) for CFPD simulations.

The simulation was performed using the licensed finite volume solver StarCCM+ (Software, 2021). The boundary conditions for the velocity $[eqn]$ were chosen to replicate a representative respiratory cycle of a human subject under heavy-exercise breathing conditions, completing a 2.4-second breathing cycle. The inlet flow was computed as in , where $[eqn]$ is the time of the simulation, and $[eqn]$ is the breathing period, which in this case will be $[eqn]$ . Lastly, $[eqn]$ is the peak air flow in a breathing cycle $[eqn]$ .

The mesh consisted of polyhedral elements, with four boundary layers and a base element size of 1.5 mm. Furthermore, refinement of the cells near the wall region was performed, resulting in an element size of approximately 1.1 mm. This base element size in the bulk flow region was chosen to ensure convergence of the simulations, per previous studies by Bartol et al. (2024). The total number of polyhedral cells resulted in 60,204 and 177,802 mesh vertices. The total volume of the lower respiratory tract model of the subject was 19.82 mL.

Within both the trachea and the primary bronchi, large eddies and flow separation arise from the abrupt change in direction at the carina ridge. These large-scale turbulent structures are captured by the $[eqn]$ Shear-Stress Transport (SST) Langtry–Menter (LM) RANS turbulence model. This four-equation model resolves transitional and fully developed turbulence and captures relaminarization in distal bronchi. Although detailed information on boundary conditions, choice of turbulence model, and fluid-particle coupling can be found in previous work by Bartol et al. (2024), the following description provides a summary of the boundary conditions used for the model.

At the outlets, the pressure was set to a fixed value, while zero-gradient conditions were applied at the inlets and walls. Turbulence inlet specification was computed with the usual isotropic estimate for the kinetic energy, $[eqn]$ . Because at heavy-exercise conditions, a medium level of turbulence is expected, the intensity $[eqn]$ was set to 4%, and the reference speed $[eqn]$ was the peak inhalation velocity at the inlet for the particular airflow used. The specific dissipation rate $[eqn]$ is then set with $[eqn]$ , using the SST constant $[eqn]$ and a length scale $[eqn]$ equal to the hydraulic diameter of the inlet. At the outlets, both $[eqn]$ and $[eqn]$ were set to the condition of zero pressure gradients.

The All $[eqn]$ wall treatment wall function was used in STAR-CCM+, which blends wall functions and near-wall resolution models depending on the local $[eqn]$ value. In our simulation, the $[eqn]$ values ranged from 0 to 5, falling within the viscous sublayer.

The intermittency $[eqn]$ is set to 1.0 at the inlet, assuming fully turbulent conditions, and uses zero gradients elsewhere in the regions. Lastly, the transition momentum thickness Reynolds number $[eqn]$ is computed from the local turbulence intensity Tu with a more complex piece-wise correlation. Walls and outlets again were set with zero gradients for both $[eqn]$ and $[eqn]$ .

A total of 105,492 particles were injected into the system using 20 randomized injector points located on the inlet surface, corresponding to the trachea inlet. Of the total original particles injected, only 17,839 deposited to the walls.

Three boundary conditions were implemented for the volume flow inlet: a baseline flow rate of $[eqn]$ , followed by configurations with 20% higher $[eqn]$ and 20% lower $[eqn]$ flow rates. An additional simulation with 10% increased flow $[eqn]$ was performed to validate the interpolation capabilities of the Koopman operator.

Data arrangement

2.2.

For particle tracking, snapshots of particle positions for the baseline simulation were stacked vertically to form a third-ordered tensor $[eqn]$ , where:

• 3 is the number of spatial coordinates (x,y,z);
• $[eqn]$ is the number of particles for the DMD study $[eqn]$ ; and
• $[eqn]$ is the number of uniform timesteps (m = 144).

Only the trajectories of particles that deposited to the walls at their last position were used; particles that were inhaled and then exhaled were omitted from the tensor definition.

Similarly, the fluid velocity field was organized into a tensor $[eqn]$ , where:

3 is the number of velocity components (u,v,w); $[eqn]$ is the number of mesh nodes $[eqn]$ ; and $[eqn]$ is the number of uniform timesteps (m = 144).

Exact dynamic mode decomposition

2.3.

The exact DMD algorithm from Kutz and Brunton (2016) was applied to both datasets, $[eqn]$ and $[eqn]$ .

For the decomposition of the particles, the lagged $[eqn]$ and shifted $[eqn]$ tensors were constructed by arranging each state vector in chronological order. Each vector contains the three spatial coordinates of each particle at time $[eqn]$ , represented here as $[eqn]$ , with each column corresponding to a different particle:

$$[eqn]$$

Here $[eqn]$ is the number of snapshots taken for the simulation, which for all purposes of this study was $[eqn]$ .

For the $[eqn]$ tensors corresponding to each of the $[eqn]$ -th breathing flow rates $[eqn]$ , where the index $[eqn]$ refers to the breathing flow rate used for the corresponding simulation. The data is arranged as $[eqn]$ , from that structure, the following tensors were then constructed as:

$$[eqn]$$

before applying DMD, each snapshot $[eqn]$ was reshaped into a column vector $[eqn]$ , forming the snapshot matrix

$$[eqn]$$

(and similarly for $[eqn]$ ). The lagged tensor is subsequently decomposed using SVD.

$$[eqn]$$ $$[eqn]$$

The matrix $[eqn]$ (and its transpose $[eqn]$ ) contains the spatial orthonormal POD modes of the simulation. The diagonal matrix $[eqn]$ holds the singular values, which represent the energy captured by each corresponding mode. The matrix $[eqn]$ (and its transpose $[eqn]$ ) contains the temporal coefficients that describe how the spatial modes in $[eqn]$ are weighted and combined over time.

The DMD operator $[eqn]$ with $[eqn]$ is approximated as $[eqn]$ by projecting it onto the subspace spanned by the truncated POD modes $[eqn]$ , retaining only a subset of $[eqn]$ dominant modes. The derivation proceeds as follows:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

Where $[eqn]$ , and $[eqn]$ . Combining these components yields:

$$[eqn]$$

The DMD modes $[eqn]$ and eigenpairs $[eqn]$ (eigenvalues) and $[eqn]$ (eigenvectors) are obtained by performing eigen-decomposition on each of the $[eqn]$ :

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

An approximation of the original data, $[eqn]$ , can be reconstructed for a dimensionless time $[eqn]$ , where $[eqn]$ represents the discrete timestep used in both the simulation and the DMD formulation, the last 16 frames of the 2.4s breathing cycle were discarded as the fluid field was nearly zero throughout the geometry, causing numerical instability in the DMD. The reconstruction is obtained by computing the initial amplitudes $[eqn]$ using the pseudoinverse of the previously computed DMD modes $[eqn]$ .

$$[eqn]$$ $$[eqn]$$

This reconstruction is made for each of the simulations by choosing the appropriate $[eqn]$ and $[eqn]$ .

Koopman operator interpolation

2.4.

A linear combination of Koopman operators (Huhn et al., 2023) was performed for flow rate interpolation between boundary conditions. Given two trained operators $[eqn]$ and $[eqn]$ with corresponding breathing flow rates $[eqn]$ and $[eqn]$ , the interpolated operator $[eqn]$ for an intermediate flow rate $[eqn]$ is:

$$[eqn]$$

where $[eqn]$ is the interpolation weight for an intermediate $[eqn]$ determined by:

$$[eqn]$$

Then, performing the eigen-decomposition of the interpolated Koopman operator to obtain the eigenpair.

$$[eqn]$$

The POD modes can be interpolated $[eqn]$ to compute the interpolated DMD modes $[eqn]$ by projecting onto the POD modes since the $[eqn]$ is not available to compute the modes via the exact way as in Eq. (14), and thus computed as follows:

$$[eqn]$$ $$[eqn]$$

Finally, the amplitudes are interpolated as: $[eqn]$ as:

$$[eqn]$$

The reconstructed field $[eqn]$ is then computed using Equation (18) with $[eqn]$ and $[eqn]$ .

All the computations described above follow an order that can be summarized in the following algorithm:

Algorithm 1: Parametric DMD algorithm via reduced (discrete) Koopman operator interpolation (Huhn et al., 2023)

Performance evaluation

2.5.

To evaluate the performance of the DMD reconstruction, three quantitative metrics were used: (1) Mean Relative Error (MRE), to assess the general precision of fluid field reconstructions; (2) the Bhattacharyya correlation coefficient, which measures the similarity between two probability distributions by quantifying their overlap (where a value of 1 indicates a perfect match and 0 indicates completely dissimilar distributions). This metric was used to compare the reconstructed and ground-truth particle deposition distributions; and (3) the maximum local relative error, to identify the worst-case deviations in spatial precision.

These metrics were applied in each phase of the study to quantify the fidelity of the reconstruction. The Bhattacharyya correlation was computed as follows:

$$[eqn]$$

Where $[eqn]$ and $[eqn]$ are the particle marginal deposition distributions at position r, obtained by CFPD simulations and DMD reconstruction respectively.

The CFD and DMD data matrices $[eqn]$ contain the velocity components (u,v,w) across all $[eqn]$ mesh nodes and $[eqn]$ time steps, where each column corresponds to one temporal snapshot. That is,

$$[eqn]$$

where $[eqn]$ are the CFD and DMD state vectors at the $[eqn]$ th time step.

The relative reconstruction error for each snapshot was computed as

$$[eqn]$$

where each $[eqn]$ quantifies the relative $[eqn]$ -error between the DMD reconstruction and the ground-truth CFD field at time step $[eqn]$ .

The mean relative error reported in Table 1 corresponds to the arithmetic average of the individual snapshot errors:

$$[eqn]$$

Finally, the maximum relative error in all of time-steps, was calculated as:

$$[eqn]$$

Results

The results are organized according to the three phases of the study: Lagrangian particle reconstruction, fluid field reconstruction, and interpolation using the Koopman operator. For each phase, performance was assessed using the defined quantitative metrics, MRE, Bhattacharyya correlation, and local relative error, alongside qualitative visualization of DMD-reconstructed fields. Additionally, computational efficiency was compared between DMD-based and full CFD simulations.

Phase I results: Lagrangian reconstruction

3.1.

For each dataset, $[eqn]$ particle trajectories were tracked until they impinged on the wall, $[eqn]$ particles. The radial position $[eqn]$ of every deposited particle $[eqn]$ was recorded. In Fig. 2, a visual comparison between the CFPD (blue) and DMD (orange) predicted trajectories was presented. The Lagrangian reconstruction analysis revealed two key aspects of DMD performance. Spatial tracking, as in Fig. 2, demonstrated that DMD successfully captured the overall trajectory patterns but accumulated positional errors in regions of complex flow, evidenced by the diverging paths in the smaller bronchi branches. The directional consistency suggested that phase relationships were preserved despite local deviations. It is important to note that the DMD model does not explicitly evaluate whether a particle remains within or exits the geometric domain, but instead reconstructs trajectories based on the spatial and temporal information captured in the DMD modes. The geometry is implicitly reflected in the particle paths, rather than being directly encoded or constrained. As a result, some deviations from the original trajectories. To observe this phenomenon, three individual particle trajectories sampled from the original tensor, were reconstructed using DMD, as shown in Fig. 3. The results show that the reconstructed paths may extend beyond the original geometry, since they do not follow the exact particle trajectory at each time step. In Fig. 4, a histogram with uniform bin width $[eqn]$ was constructed and normalized:

$$[eqn]$$

where $[eqn]$ is the count in bin $[eqn]$ .

The same binning and normalization were applied to (a) the Euler–Lagrange CFD data and (b) the DMD-predicted particle field, which was sampled with identical binning resolution.

The Position $[eqn]$ axis in Fig. 4, where $[eqn]$ denotes the radial distance from the carina–centroid axis. The marginal probability-density function (PDF) $[eqn]$ is defined such that,

$$[eqn]$$

i.e. $[eqn]$ gives the probability that a particle deposits between $[eqn]$ and $[eqn]$ . The probability density $[eqn]$ can take values ≫ 1, meaning high particle concentration in narrow bins or ≪ 1, meaning low particle concentration in those bins.

Since only particles that ultimately adhered to the airway walls were considered in the analysis, all reconstructed trajectories correspond to deposited particles, with their final positions indicating their deposition sites. In the study, no particles that had escaped through the outlets or inlet were considered. Another analysis focused on how well the DMD model captured the particle distribution based on a Euclidean distance metric, Eq. (31), with respect to a fixed point ( $[eqn]$ located in the Carina) used for both the ground truth particle simulation and the DMD reconstruction. The distribution analysis, shown in Fig. 4, demonstrated a strong agreement with the ground truth, with DMD reproducing the overall deposition density profile in all radial positions with a Bhattacharyya correlation coefficient of 85%. The presented method captured peak deposition in the carinal ridge, while showing slight underprediction in the smaller bronchi or more distant regions, likely due to the fact that particles traveling farther down the airway follow more convoluted trajectories. In conclusion, the presented DMD approach prove more effective for predicting statistical deposition patterns than for tracking individual particle trajectories, completing the reconstruction in under 7 min using full-ranked Koopman operators $[eqn]$ , whereas the ground truth simulation required approximately 60 h to model a complete 2.4-second breathing cycle.

$$[eqn]$$

Phase II results: Fluid field reconstruction

3.2.

The fluid field reconstruction phase demonstrated DMD’s ability to capture the main flow features under varying inlet conditions. While DMD successfully resolved recirculation zones and dominant flow patterns (Fig. 5), maximum errors arose in regions where flow separation occurs because of the abrupt change in the geometry; these are typically locations with steep velocity gradients (e.g., bronchial branching). A worst-case analysis is presented in Fig. 6, revealed localized discrepancies near flow bifurcations (flow-splitting regions), where error magnitudes scaled with flow complexity. It is noteworthy that errors never exceeded 15% at any single-cell level across all inlet flow conditions.

Notably, although DMD preserved main flow features, such as flow separation in high-curvature regions, vortical recirculation, and near-stagnant regions, across all volume flow rates, it also introduced progressive damping of the overall kinetic energy. This smoothing effect, observed consistently across cases, probably originated in the rank reduction step in DMD, where the less dominant modes were truncated to increase computational efficiency, resulting in a total reconstruction time of ~10 min, including file generation.

Phase III results: Koopman interpolation

3.3.

The Koopman operator interpolation highlighted the challenges of projecting nonlinear system dynamics onto a finite-dimensional subspace. By combining modes from 100% and 120% input volume flow rate simulations, the method preserved topological features of the flow field, such as the primary recirculation zones in the modeled parts of the tracheobronchial region (G0 to G5), including the yellow region in Fig. 7, and bulk flow direction. However, the large interpolation interval (20%) struggled to resolve the geometric complexity of the domain, leading to two critical issues: (1) Progressive damping of the mean velocity (loss of kinetic energy fidelity) and (2) inaccurate reconstruction of recirculation regions in bronchial flow-splitting zones (bifurcations). While the mean velocity profile maintained correlation, localized errors concentrated in regions sensitive to volume flow variations, notably bifurcation zones and boundary layer near walls. This suggested that while Koopman interpolation preserved essential nonlinear behavior through its infinite-dimensional projection, practical applications in respiratory flows require training data with smaller parameter intervals to resolve geometrically induced flow instabilities. This necessity reduces the original advantage of the method of reducing computational burden through sparse sampling, as a finer parameter resolution requires more expensive simulations due to a smaller time step.

Table 1 provides a summary of the maximum absolute relative local error and the MRE in all snapshots for each of the conditions studied previously together with the interpolation case and a localized measurement of the error on the walls and outlets of the model. All reconstructions were performed using full-ranked (r = 144) Koopman operators.

The interpolation scheme was also applied to the deposition distributions from simulations at 100% and 120% input volume flow rates. Its performance was comparable to the direct reconstruction results shown in Fig. 4, achieving a Bhattacharyya correlation coefficient of 80%. The corresponding results are presented in Fig. 8.

Discussion

This study investigated the application of DMD in fluid and particle simulations as a reduced-order model, i.e., as a framework in which simulations can be reconstructed at a fraction of the computational cost.

In the first part of this study (Section 3.1), the sensitivity of DMD to trajectory reconstructions was demonstrated. The distribution of particles deposited on the walls at the end of the simulation is strongly correlated with the ground truth, although the exact positions and trajectories are not well reconstructed. This behavior is likely attributable to the lack of coupling among particles, despite their overall path similarity. Each row of the input tensor corresponds to an individual particle and is treated independently, without enforcing inter-particle constraints. As shown in Fig. 3, the reconstructed trajectories diverge from the ground-truth paths, and in Fig. 2, several particles are observed outside the boundaries of the original three-dimensional geometry. In contrast, when simulating fluid fields (Sections 3.2 and 3.3), the momentum equations tightly couple all the nodes (rows) of the input matrix, thereby enabling the DMD framework to better capture the system dynamics and achieve a lower reconstruction error.

Because training excluded non-deposited trajectories, the DMD cannot predict whether a particle exits the domain or remains in suspension; it only localizes deposited particles. Estimation of deposition probability/fraction by size or release location should be provided by a complementary model (e.g., ICRP/MPPD or a dedicated classifier), after which our DMD yields spatial localization for the predicted depositors.

Most data-driven approaches are unbounded by governing equations, making them prone to divergence when the input lies outside the vicinity of the training set. However, the results in Section 3.3 show that DMD produces outcomes consistent with neighboring datasets. A 10% variation in the initial conditions results in an increase in mean reconstruction error by comparable margin, with localized errors increasing by nearly sevenfold.

In the study by Habibi et al. (2020), the authors applied DMD to a pulsatile blood flow system (an oscillating system similar to a breathing cycle). However, they imposed a boundary condition at the pressure inlet consisting of a sum of multiple overlapping complex exponentials (i.e., sines and cosines), making their system similarly intricate. They validated their results against the Womersley solution, which served as the ground truth. Notably, (Habibi et al., 2020) found that reducing the system to only 10 modes resulted in an error of less than 10%, comparable to the error obtained in this study, as shown in Table 1, with the exception of the interpolation case, which involves additional numerical error sources.

In the work by Jing et al. (2024), the authors concluded that DMD is a particularly suitable framework for reconstructing respiratory flows, as it captures single-frequency dynamics. This ability allows DMD to isolate the components comprising the entire simulation and to offer a different perspective on particle transport and deposition. These conclusions are consistent with the results presented in Sections 3.1 and 3.2. Together, these findings reinforce the potential for reduced-order dynamic mode representations to approximate high-fidelity simulations in inhalation dosimetry applications, while also underscoring the importance of spatial and temporal resolution constraints in ensuring model robustness.

Conclusion

This comprehensive analysis of DMD in respiratory flow modeling revealed its dual strengths and limitations across three applications. The method achieved 85% Bhattacharyya correlation with ground-truth particle deposition patterns and showed strong agreement in fluid field reconstruction within trained parameter regimes, significantly reducing computational costs. DMD was particularly effective at extrapolating dominant flow features and preserving large-scale structures such as recirculation zones, and it demonstrated robustness in its reconstruction capabilities. As a comparison, the execution time for a full particle–fluid dynamics model with two-way coupling, utilizing 24 Dual Intel Xeon Gold 6226 CPUs @ 2.7 GHz, DDR4-2933 MHz DRAM, and 32 GB of RAM, required ~60 h to simulate a complete 2.4-second breathing cycle on a mesh with 60,204 cells and 177,802 mesh vertices. Conversely, a DMD-based interpolation of both particle trajectories and the fluid field, performed using 8 CPUs and 32 GB of RAM, achieved the same reconstruction in ~10 min. These findings position DMD as a promising reduced-order modeling framework for efficiently approximating respiratory flow and particle dynamics, supporting future applications where rapid, geometry-specific simulation is needed. However, its ability to interpolate intermediate inlet conditions, such as the untrained +10% inlet flow case, highlighted inherent constraints when applying DMD to complex geometries. The 20% parameter interval between training points led to localized errors exceeding 70% and a mean relative error across all time steps of ≈ 20%, underscoring that while DMD effectively projects nonlinear systems into linear Koopman spaces, geometrically induced flow complexities demand tighter parameter spacing for reliable interpolation. As a result, its applicability as an interpolator in sparsely sampled parametric simulations is limited.

Nonetheless, direct particle interpolation between the high-fidelity simulations at 100% and 120% inlet volume flow rates performed comparably to the direct DMD reconstruction of the full high-fidelity simulation, achieving a Bhattacharyya correlation coefficient of 80%.

While the proposed DMD and Koopman-based methods demonstrate promising reconstruction capabilities, several limitations remain. First, the simulations are limited to a single geometry and breathing condition, constraining generalizability. Second, particle–particle interactions were not modeled beyond two-way coupling with the fluid, potentially affecting accuracy in high-density aerosol environments. Third, the interpolation framework assumes linearity in the reduced-order space, which may not hold for all dynamic regimes. Future work will address these limitations through broader parametric variation, incorporation of subject-specific airway geometries, and development of non-linear reduced-order models and multi-resolution approaches.

These findings suggest that DMD’s current utility lies primarily in regime-specific reconstruction rather than broad parametric prediction in anatomically intricate domains. Future implementations could combine DMD with adaptive parameter sampling or physics-informed constraints to address these gaps, potentially enabling real-time analysis of patient-specific aerosol deposition without sacrificing critical flow dynamics. Incorporating such strategies may extend the applicability of reduced-order modeling frameworks like DMD in clinically relevant respiratory simulations, where computational efficiency and spatial resolution are both required to capture physiologically meaningful flow and deposition patterns.