- Email: [email protected]
Topological analysis of particle transport in lung airways: Predicting particle source and destination
Journal Pre-proof Topological analysis of particle transport in lung airways: Predicting particle source and destination Ali Farghadan, Filippo Coletti, Amirhossein Arzani
PII: DOI: Reference:
S0010-4825(19)30362-2 https://doi.org/10.1016/j.compbiomed.2019.103497 CBM 103497
To appear in:
Computers in Biology and Medicine
Received date : 6 September 2019 Revised date : 9 October 2019 Accepted date : 10 October 2019 Please cite this article as: A. Farghadan, F. Coletti and A. Arzani, Topological analysis of particle transport in lung airways: Predicting particle source and destination, Computers in Biology and Medicine (2019), doi: https://doi.org/10.1016/j.compbiomed.2019.103497. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Ltd.
Journal Pre-proof
The author hereby declares no conflicts of interest.
Jo
urn al P
re-
pro
Amirhossein Arzani Assistant Professor Cardiovascular Biomechanics Lab https://www2.nau.edu/cardio-lab/ Mechanical Engineering Department Northern Arizona University
of
Conflict of Interest Statement
0
Jo
FTLE
6
33
map
Destination
Graphical Abstract (for review)
urn al P
re-
pro
of
Journal Pre-proof
1
tL = 0.0s
tL = 0.2s
tL = 0.4s
Color-coded outlets
Journal Pre-proof *Revised manuscript (clean) Click here to download Revised manuscript (clean): Farghadanetal19_revise_notred.pdf
Filippo Coletti2
Amirhossein Arzani1
pro
Ali Farghadan1
of
Topological analysis of particle transport in lung airways: predicting particle source and destination
1
Jo
urn al P
Correspondence: Amirhossein Arzani, Northern Arizona University, Flagstaff, AZ, 86011 Phone: +1-9285230541 Email: [email protected]
re-
Department of Mechanical Engineering, Northern Arizona University, Flagstaff, AZ, United States 2 Department of Aerospace Engineering and Mechanics, University of Minnesota, Minneapolis, MN, United States
1
Journal Pre-proof
Abstract
Jo
urn al P
re-
pro
of
Particle transport in lung airways can induce respiratory disease and play a vital role in aerosol drug delivery. Herein, we present dynamical systems features that influence airflow and particle transport in the tracheobronchial trees. Computational fluid dynamics (CFD) was used to solve for unsteady airflow in a patient-specific model. Particle tracking simulations were performed for micron-size particles. The destination map that connects the particle final location to the initial location and injection time was constructed. Finite-time Lyapunov exponent (FTLE) fields were calculated to identify inertial Lagrangian coherent structures (ILCS), topological features that act as separatrices. Our results demonstrated that these topological features control the destination map at the trachea. The temporal evolution of ILCS influenced the sensitivity of particle transport fate to injection time, whereas the emergence of new ILCS with an increased integration time controlled transport to different generations of airways. Additionally, particles starting at the ILCS were shown to mostly deposit at the airway walls. Finally, an innovative source inversion strategy was introduced to integrate the MaxeyRiley equation backward in time and identify the origin of dispersed particles. Our study explores novel dynamical systems tools that improve our understanding of particle transport and deposition in the airways and could be used to guide future targeted drug delivery studies. Keywords: airflow, computational fluid dynamics, FTLE, Lagrangian coherent structures, source inversion, particle tracking, aerosol.
Journal Pre-proof
1
Introduction
Jo
urn al P
re-
pro
of
The essence of particle transport in the respiratory system is well known. Harmful substances such as dust, fumes, and gaseous pollutants are vastly suspended in the environment with a variety of shapes and sizes ranging from ultrafine (smaller than 1 µm) to large-size matters (larger than 100 µm) [5]. Detrimental substances have either immediate effect and cause minor disorders such as dizziness or can cause serious long-term damage to the lungs including asthma or cancer [19, 48]. In the case of many of the developed chronic lung diseases, prescribing an effective treatment is necessary. Overall, for both medical and pathological purposes, studying particle transport provides diagnostic and prophylactic measures for adverse respiratory health conditions. To study particle transport, computational fluid dynamics (CFD) coupled with particle dynamics simulations are widely employed [23, 35, 25, 56, 38]. The dichotomous network of airways represents one of the complicated dynamical systems in the human body where chaotic flow and advection occurs [49]. Several characteristics contribute to complex airflow patterns such as lung geometry with numerous planar and non-planar conducting tubes, unsteady nature of the breathing cycle, asymmetric tapering down into distal branches, and other physiological factors [30, 6, 51]. All of these features are known to influence particle transport and deposition patterns [7, 12, 31, 20]. In practice, Lagrangian approach and in particular solving the Maxey-Riley (MR) equation [33], which governs the dynamics of small spherical inertial particles, has been shown to be an effective method to simulate particle transport in lung airways [32, 36, 27]. This equation is extensively studied in a broad range of applications for all types of particles such as aerosols, neutrally buoyant particles, and bubbles [21, 46]. In the context of aerosol transport in lung airways, this equation is used to answer two important questions: First, what is the fate of each particle entering the upper airways? Second, what spatial spot and temporal time-point should particles be released to reach a specific desired location? The first question is of interest regarding harmful inhaled particles, while the latter question is one of the major problems involved in targeted drug delivery. To answer the aforementioned questions, particle destination maps have been proposed [24, 54, 26]. Namely, particles are densely seeded at the inlet of an airway model and integrated based on the MR equation. Subsequently, the destination of each particle (typically based on the outlets of the model) is identified and mapped back to the initial location and time-point of release. Such destination maps are used to study the effect of spatial location [59] and intra-cycle time-point of release [13] on particle transport and fate. Additionally, the effect of particle diameter (or more generally the Stokes number) on destination maps has been studied [58]. While destination maps are useful in identifying the fate of injected (or incoming) particles, they do not provide a mechanistic explanation of particle transport results and are often difficult to interpret. Moreover, if one is interested in targeted drug delivery to a specific location of the lung airways, a large number of particles are often required to be released and tracked such that the localized region of interest is sufficiently sampled. In this study, we propose novel dynamical systems tools to overcome these issues and improve our understanding of particle transport in the conducting airways. First, we perform a topological analysis of fluid flow and particle transport using finite-time Lyapunov exponent (FTLE) fields and Lagrangian coherent structures (LCS) [42]. LCS mark boundaries of flow with distinct dynamical behavior and could be considered as invisible surfaces embedded in fluid flow that act as transport barriers. Therefore, they provide a template to study transient and chaotic flows. Herein, we compute the extension of this concept to inertial particles (inertial LCS) [40, 46]. Our group and others have used FTLE and LCS to characterize transport in cardiovascular flow applications [43, 4, p. 3
Journal Pre-proof
2.1
Methods Airflow simulations
urn al P
2
re-
pro
of
3, 41, 14]. In respiratory flows, only one group has studied the significance of FTLE in the airways and illustrated the correlation between FTLE and destination maps [45, 44]. However, the effect of temporal FTLE variations and integration time on particle distribution as well as the effect of LCS on particle deposition was not revealed. In this study, we provide further insight into how FTLE maps could fully characterize particle transport in the tracheobronchial tree. Finally, we introduce a novel source-inversion technique to identify the initial position of dispersed particles with very little computational effort. Namely, we use the recent slow manifold technique developed by Haller and Sapsis [17] to integrate the MR equation backward in time and identify the initial location of release for any desired set of particles. The manuscript is organized as follows. First, we describe our CFD and particle dynamics frameworks for a patient-specific model under unsteady flow conditions. Next, we explain FTLE and inertial LCS calculation as well as source inversion using the slow manifold method. We discuss how these techniques could be used to study the spatiotemporal distribution of particle destination maps. We demonstrate how the integration time in FTLE calculation and the resulting LCS correlate with particle distribution to different generations of bifurcating airways. Finally, we demonstrate how backward integration of the MR equation could be used to identify the particle injection location, which is of high interest in targeted drug delivery.
The inhalation and exhalation simulation is performed on the same model employed in an invitro study [22]. The inlet of the model is extended at the trachea. A parabolic, Newtonian g ), incompressible (ρf = 1.2 × 10−3 cmg 3 ), unsteady airflow is imposed at (µf = 1.81 × 10−4 cm.s the inlet. The airflow is assumed to be sinusoidal (Q = Qmax Sin(2πf qt)) with peak Reynolds ρf V¯max din ) of 1500 and Womersley number (W o = d2in 2πf ) of 2.6, where νf is number (Remax = µf νf
Jo
air kinematic viscosity, din = 1.6 cm is the inlet diameter and f = 15 respiration is the respiratory min frequency where each breathing cycle takes 4 seconds. The corresponding cross-section averaged L peak velocity is V¯max = 141.4 cm and peak flow rate is Qmax = 17.1 min . The nondimensional s Remax and Wo numbers belong to a normal breathing rate in an adult. The flow is assumed laminar throughout the tracheobronchial tree. No-slip boundary condition is prescribed at the wall. For the inhalation, resistance boundary conditions are imposed at the outlets and are tuned iteratively such that the flow rate distributions are 14, 5, 32, 17, 32 percent in right upper lobe (RUL), right middle lobe (RML), right lower lobe (RLL), left upper lobe (LUL), and left lower lobe (LLL), respectively. These values follow the lobar flow rate distributions obtained by multidetector computed tomography in a cohort of healthy non-smokers [52]. For the exhalation, a plug velocity profile is assigned to each terminal branch matching the unsteady exhalation flow rate. The simulation time step for the whole reparation cycle is 2 × 10−4 s and the results are found to be time-step independent with this selection of time-step. An overview of the 3D model and boundary conditions are shown in Fig. 1. The mesh consisted of 5.3M quadratic tetrahedral elements with five boundary layers created with SimVascular [50]. The boundary layer meshing is used to capture the near-wall flow phenomena. This second-order mesh is equivalent to a first-order mesh with ∼42M linear tetrahedral elements, resulting in an extremely fine spatial resolution compared to other highresolution studies in the literature (mostly less than 15M elements [39, 38, 55]). Unsteady velocity data are obtained using Oasis [34], an open-source Navier-Stokes solver (based on the finite element p. 4
Journal Pre-proof
method) with minimum numerical dissipation. The second respiratory cycle (t = 4 to 8s) is saved for post-processing. The focus of this study is on the inhalation phase and the exhalation phase simulation is performed to provide more realistic flow conditions at the beginning of inhalation.
Exhalation Phase
Vinh(t)
PTrachea
Trachea seeded surface
PRUL RRUL
RLUL
17% Vexh(t) RRML
PLLL
14% Vexh(t)
pro
PLUL
32% Vexh(t)
re-
PRLL
32% Vexh(t)
(b)
urn al P
(a)
RML seeded surface
5% Vexh(t)
PRML
RRLL
RLLL
Seeded surfaces
of
Inhalation Phase
(c)
Figure 1: The patient-specific lung airway geometry used in this study is displayed. The CFD boundary conditions for (a) inhalation and (b) exhalation phases of the breathing cycle are shown. (c) The location of a thin box at the trachea where particles are launched and FTLE fields are calculated is shown. Also, the location of seeded particles at the entrance of right middle lobe (RML) where source inversion study is performed is displayed.
2.2
Particle tracking
Particles are tracked within the domain of interest throughout the inhalation phase. The governing equation of spherical particle-laden flows (the MR equation) is integrated in a Lagrangian approach
Jo
dv(x(t)) 3 Du(x, t)) 18µf 1 = (ρp − ρf )g + ρf − 2 (v(x(t)) − u(x, t))) , (ρp + ρf ) 2 dt 2 Dt dp
(1)
where g is the gravity, facing downward resembling an upright standing, v and u are particle and air velocities, respectively. µf is fluid viscosity, and ρp and ρf are particle and air densities, respectively. Stokesian drag is assumed, which is justified by the very small Reynolds number based on the particle diameter and particle velocity relative to the surrounding fluid (Rep ∼ O(10−3 )). Particle-particle interactions and back-reaction on the carrier fluid are ignored. Faxen correction and Basset history terms are also neglected. We also ensured that Brownian motion was ignorable (mainly due to the domination of advection over diffusion) for microscale particle transport in the tracheobronchial tree. Particles are seeded at a trachea cross-section right before the model extension with a velocity equal to the fluid velocity at releasing location as shown in Fig. 1c. Particles are launched at different time-points and uniformly from all spots (∼770,000 location points) on the trachea surface p. 5
Journal Pre-proof
Destination maps
pro
2.2.1
of
(the upper green slice in Fig. 1c). Two extreme limits of microparticles typically entering the tracheobronchial tree are simulated, assuming particle diameter of dp = 1 and 10 µm. The bulk portion of large particles (over 10 µm) cannot reach the trachea [11], and nanoscale particles are ρ d2 w not considered in this study. The Stokes numbers (St = 18µp fpdin , w is the mean inhalation velocity at the trachea) at the inlet corresponding to these diameters are St = 1.8 × 10−4 and 1.8 × 10−2 , respectively. Particles are integrated until they either leave the computational domain, deposit on the bronchial wall (the distance between particle center and wall is less than particle radius), or the inhalation cycle is completed.
2.3
re-
Destination maps are a useful tool that give a better understanding of particle fate. To obtain the destination map, first, particle tracking is performed as explained above. Each outlet is tagged with an index (n = 1, 2, ...). For each particle tracked, the outlet that it exits the domain is determined. Subsequently, the outlet index is assigned to the initial release location of the particle at the trachea. The resulting map of indices at the trachea is called the destination map. All deposited particles are assigned a zero index (n = 0).
Finite-time Lyapunov exponent (FTLE) and inertial Lagrangian coherent structures (ILCS)
urn al P
Lagrangian coherent structures (LCS) are invisible material surfaces embedded in 3D flows (lines in 2D flows) that have distinct attracting or repelling properties [16]. The attracting LCS attract nearby particles in their basin of attraction and could be sometimes visualized in the lab using experimental visualization techniques. On the other hand, the repelling LCS repel nearby trajectories and act as transport barriers. In this study, we only compute the repelling LCS and demonstrate their influence on particle distribution. While mathematically more accurate techniques have been proposed for identification of LCS [15], the finite-time Lyapunov exponent (FTLE) field [42] remains a computationally robust approach in identifying LCS. FTLE is utilized to characterize coherent structures in a wide range of studies [8, 9, 1, 41]. FTLE measures the exponential rate of stretching of initially closely seeded particles. The ridges of the FTLE field (loosely speaking, local maxima) identify the repelling or attracting LCS under forward and backward time integration, respectively. To calculate the FTLE field, a densely distributed set of particles is seeded in a region of interest and advected via integration of the MR equation. The FTLE field is subsequently computed as [42]
Jo
σ(x0 ; t0 , t) =
1 p ln λmax C(x0 ; t0 , t) , T
(2) |
where λmax C(x0 ; t0 , t) is the maximum eigenvalue of C(x0 ; t0 , t), C(x0 ; t0 , t) = ∇Ftt0 · ∇Ftt0 is the right Cauchy-Green strain tensor, Ftt0 : x(t0 ) 7→ x(t) is the flow map (maps the initial particle location to the final location), T = t − t0 is the integration time, and the FTLE field is calculated as a scale of the largest eigenvalue of the Cauchy-Green strain tensor. The FTLE field (σ(x0 ; t0 , t)) is a scalar field and in unsteady flows it can vary depending on the initial time-point of particle release (t0 ). The division over integration time (T ) in Eq. 2 could be dropped to avoid sensitivity of FTLE values to particles leaving the computational domain early [4]. The FTLE field is computed for a thin box in the same location where particles are initially seeded. The calculation is only confined to a thin box at the inlet because we are interested in the trajectory of particles released at this region. Over 1.5M particles with diameters of 1 and 10 µm p. 6
Journal Pre-proof
2.4
pro
of
are launched in six time-points of the inhalation phase. Forward time integration is performed, and therefore the ridges of the FTLE field identify the repelling LCS. Additionally, since we are tracking inertial particles (instead of ideal tracers) the inertial Lagrangian coherent structures (ILCS) [40, 46] are identified. The main difference between ILCS and LCS is in the flow maps. Constructing flow maps by solving the MR equation identifies the ILCS, whereas the LCS are identified by integrating ideal tracers along velocity pathlines. To gain a better understanding of coherent structures inside the lung, the FTLE field is calculated for a 3D box spanning the trachea to the carina region. Calculating FTLE for a large 3D box is computationally costly since millions of particles (in this example, 27M particles densely seeded in the 3D domain) need to be integrated by solving the MR equation.
Source inversion
re-
The Maxey-Riley (MR) equation, mathematically, could be integrated forward and backward in time. However, it has been proven that backward integration of the MR equation for locating the origin of dispersed particles (source inversion) is quite challenging. In fact, for large-scale timedependent flows, even a high-precision integration scheme with very fine time-step blows up when the MR equation is integrated backward in time. To overcome this instability, Haller and Sapsis [17] proposed the idea of slow manifolds under unsteady flow conditions, which is briefly explained here. We can rewrite the MR equation (Eq. 1) by dividing both sides by the coefficient of particle acceleration (ρp + 12 ρf ) as
18µ
urn al P
Du(x, t)) dv(x(t)) = (1 − 1.5R)g + 1.5R − KR(v(x(t)) − u(x, t))) , dt Dt
(3)
2ρ
Jo
f where K = d2 f and R = ρf +2ρ . This equation is singular at the limit of infinitesimally small p p particle diameters. Hence, integrating the MR equation backward in time is numerically impossible when KR >> 1. The term −KR(v(x(t)) in equation 3 causes a quick and inevitable blowup of the numerical solution to this equation backward in time. Using high-order numerical integration schemes or smaller time step would only delay the blowup process, making backward MR integration impractical for large-scale time-dependent problems. Time-dependent invariant slow manifolds (M ), which are three-dimensional global attractors, are proven to exist in the six-dimensional phase space of the MR equation [17]. M can be written as the following form: Du(x, t)) 2 − g + O( ) , (4) M = (x, v) : v = u + (1.5R − 1) Dt
ρ d
where = ρfp win St quantifies the significance of inertia. Note that in the original paper [17], = St R R was written in non-dimensional form and a different definition of St number was used. In order to recover the starting position of an inertial particle without struggling with the instability of solving the full MR equation backward in time, one must first project the final location of a dispersed particle onto the slow manifold along its velocity direction. Then the governing particle motion equation on slow manifolds should be integrated from the final time to the initial time [17]. A first-order approximation to the particle motion on slow manifolds is given by Du(x, t)) x˙ = u(x, t) + (1.5R − 1) − g + O(2 ) , (5) Dt p. 7
Journal Pre-proof
3
pro
of
which is the lowest-order truncation without numerical divergence. The zero limit of makes Equation 5 equivalent to the governing equation of ideal tracers. Sufficiently large values of , however, might lead to inaccuracy. For microscale particle tracking in lungs << 1, and Equation 5 is easily integrable backward in time with good accuracy. Using this method, the initial positions of dispersed particles can be recovered with an error of O(2 ). In this study, as an example of the utility of this method in lung airways, the source of particles reaching the right-middle lobe (RML) have been investigated. The same approach can be utilized for any other lobes or branches. Figure 1c shows the surface of interest at the RML entrance. Particles were released uniformly from over 62,500 spots seeded on this surface. They were injected 380 times during the inhalation cycle (uniformly from 2.0 to 0.1s), thus overall, ∼23.75M particles were tracked and integrated backward with a time step of 5 µs.
Results
3.1
FTLE maps
re-
The corresponding FTLE, destination, and source inversion maps are presented in the following sections.
3.1.1
urn al P
In this section, variations in the FTLE maps at the trachea are found for various cases such as different integration times, several particle injection time-points, and two particle diameters. The results are compared to the destination maps. Effect of temporal integration and particle diameter
3.1.2
Jo
In Figure 2, FTLE plots associated with 1 and 10 µm diameter particles are displayed. Particles are released at the start of inhalation and FTLE fields are calculated for several integration times (T). As the integration time increases, more ILCS emerge corresponding to more distal branches. That is, additional ILCS are created with an increase in integration time, since particles have more time to pass through new generations of bifurcations, and therefore experience separation. The longest integrated FTLE fields (rightmost column) remains qualitative invariant by extending integration time for both diameters. The near-wall region indicates high FTLE values. The reason behind these high values is the high shear rates in these regions that make the initially close particles far apart from each other in time. Moreover, some portion of near-wall particles deposits on the immediate curved wall below the trachea which magnifies these artificially high FTLE values. Finally, the figure shows that FTLE ridges (the ILCS) are robust to variations in particle diameters between 1 and 10 µm. Effect of injection time
FTLE maps are calculated for 1 µm particle diameter launched in different time-points during inhalation. As seen in Figure 3, with an increase in integration time (T), ILCS appear sooner (nearby particles separate sooner) when particles are released later in the breathing cycle. For a fixed value of the integration time, variations in the FTLE field for different launch/injection times (tL ) show the temporal evolution of the FTLE and ILCS, which demonstrate the sensitivity of particle fate to time-point of injection. The FTLE field varies with increasing the injection time up to tL = 0.4s, and is relatively insensitive to further increase in tL . p. 8
pro
of
Journal Pre-proof
3.1.3
Destination maps
re-
Figure 2: The effect of integration time (T) and particle diameter on the FTLE fields is shown. The first and second row correspond to 1 and 10 µm particle diameters, respectively. The same FTLE color bar range is used for both rows. The emergence of new FTLE ridges (ILCS) is shown with an increase in integration time.
3.1.4
urn al P
The particle destination maps are compared with ILCS for 1 µm particle size as shown in Figure 4. The deposited particles are excluded from the destination maps. A very strong correlation between the ILCS and FTLE field and the corresponding destination map is seen. Each region bound inside ILCS lines (with lower FTLE values) corresponds to a different tag in the destination map. A very close match exists between ILCS lines and boundaries lying on the destination maps for the three injection time-points shown (which are representative of the considered values of tL within the inhalation phase). Deposition maps
ILCS surfaces computed from forward FTLE fields are repellers and act as transport barriers. Therefore, virtually no particle is expected to pass from one side of ILCS to the other. The ILCS partition transport to different branches and this leads to a very interesting scenario in which particles right on the ILCS deposit on the wall. As seen in Figure 5, the original location of ultimately deposited particles is right on top of the ILCS. The same scenario occurs for different injection times as shown in the figure. 3D FTLE visualization
Jo
3.1.5
As particles are injected at the trachea, it suffices to calculate FTLE fields for a thin slice at the same location where the corresponding ILCS surfaces are visualized as lines (shown in the above figures). However, FTLE calculation throughout an extended region of the airways could be performed to visualize the 3D ILCS surfaces. This is done in Figure 6 where the 3D FTLE and ILCS patterns are shown. The results are shown for 1 µm particle diameter and 0.24s integration time when the first ILCS emerge at the trachea. One can see how the transport barrier represented by the ridge of high FTLE extends down the trachea, and into the main carina.
p. 9
T = 0.14s
T = 0.24s
T = 0.28s
T = 0.12s
T = 0.14s
T = 0.02s
T = 0.06s
T = 0.08s
0
T = 0.04s
T = 0.06s
T = 010s
T = 0.08s
urn al P
T = 0.02s
T = 0.16s
T = 0.02s
T = 0.02s
T = 0.04s
T = 0.04s
T = 0.20s
re-
FTLE
6 T = 0.02s
T = 0.32s
T = 0.06s
T = 0.06s
T = 0.08s
T = 0.08s
T = 0.12s
T = 0.10s
tL = 0.0s
T = 0.40s
pro
T = 0.02s
of
Journal Pre-proof
tL = 0.2s
T = 0.40s
tL = 0.4s
T = 0.40s
tL = 0.6s T = 0.40s
tL = 1.0s T = 0.10s
T = 0.40s
tL = 1.2s T = 0.10s
T = 0.40s
3.2
Jo
Figure 3: The temporal evolution of FTLE fields for several integration times is shown. For each integration time (T), FTLE fields are calculated for different injection times (tL ) to demonstrate the temporal evolution of FTLE. All calculations are performed for 1 µm particle diameter. The same FTLE color bar range is used for all results.
Source inversion
By launching particles seeded at RML entrance and tracking them backward in time, particles reach the trachea in a range of times. Some particles never find their way back to the trachea as they are seeded in a region that is not sampled by the forward MR equation. That is, no particle injected at the trachea reaches these locations, indicating that such locations are virtually impossible to reach under the present flow conditions and range of injection time-points. In order to create a source map of the target RML at the trachea in a specific time (tref ), only particles that arrive at the trachea within tref ± 5ms range are kept. That is, particles at the RML are integrated backward in time and the union of particles reaching the trachea at the p. 10
Journal Pre-proof
FTLE
6
of
33
pro
map
Destination
0
1
tL = 0.0s
tL = 0.2s
tL = 0.4s
Color-coded outlets
urn al P
re-
Figure 4: The comparison between FTLE and destination maps are shown for dp = 1 µm. The top row shows the FTLE field with T = 0.40s integration time and the bottom row is the corresponding destination maps. The three panels correspond to 0.0, 0.2, and 0.4s injection times from left to right, respectively. The destination map color bar is enumerated from 1 to 33. 1-12 are left lower lobe (LLL), 13-18 are left upper lobe (LUL), 19-26 are right lower lobe (RLL), 27-29 are right middle lobe (RML), and 30-33 are right upper lobe (RUL). The blue color family (from blue to green) corresponds to the left side and the red color family (from yellow to red) corresponds to the right side of the lung. desired time range is identified as the source map (defined at the trachea). This time-thresholding process is carried out for all 380 series of backward injections. Figure 7 exhibits the source maps versus RML destination maps for tref = 0.2, 0.4 and 1.0 second time-points of the respiration cycle. The similarity between the corresponding map pairs confirms the qualitative accuracy of the source inversion method. To further quantify the accuracy of the method, the full MR equation is integrated forward in time for particles released from the source maps of various tref values. Table 1 shows that between 90% to 99% of particles reach the target RML.
Discussion
Jo
4
In this study, a dynamical systems approach was used to study particle transport in lung airways. The FTLE field was calculated and a close connection between the identified ILCS and destination map was demonstrated. Destination map could be created for any cross-section at the trachea to visualize the fate of seeded particles. The ridges of the FTLE field (the ILCS) mark the boundaries of different destinations and provide a mechanistic explanation for the destination map patterns. We also, for the first time, demonstrated the utility of integrating the MR equation backward in time using slow manifolds to create source maps. The source map demonstrates the spatiotemporal origin of dispersed particles. A very good match between the source and destination maps was observed. The connection between FTLE and destination maps in lung airways was previously demonp. 11
re-
pro
of
Journal Pre-proof
urn al P
Figure 5: The FTLE map (left panel) and the initial location of deposited particles (middle panel) are shown for dp = 1 µm. The two fields are superimposed on top of each other and shown in the right panel. Three different injection times are shown in this figure. T = 0.4s is selected when the FTLE patterns become independent of T.
Jo
strated [45, 44]. Our study provides several additional insights. First, we have demonstrated a relationship between the integration time used in calculating FTLE, the resulting ILCS, and particle transport to different generations of bifurcations. Namely, as we increase the integration time, the first emerging ILCS at the trachea marks the boundary between particles separating at the first bifurcation. As the integration time is further increased, new ILCS emerge that correspond to the next generation of airways and bifurcations. We have also demonstrated the temporal evolution of the FTLE and ILCS. That is, by injecting new sets of particles at new time-steps the temporal evolution of these fields could be studied. The temporal evolution of ILCS was shown to control the changes in the destination maps with respect to the time-point of particle release (Fig. 4). This is valuable information since it allows one to track the evolution of ILCS during inspiration to study how sensitive the destination map is to particle release time-point. Interestingly, at the beginning of inspiration when the flow rate is lower, the temporal evolution of ILCS is more pronounced, whereas with an increase in flow rate the ILCS minimally changes in time. These results imply that releasing particles at the beginning of inspiration is more sensitive to the choice of particle release time compared to later stages of inhalation. Finally, we demonstrated minor changes in ILCS results when the particle diameter changed from 1 to 10 µm. Since the Stokes numbers are within a small range (O(10−4 ) to O(10−2 )), the bulk of particles follow the fluid flow, and therefore the FTLE patterns look very similar. We verified that qualitative LCS patterns were identical between the microparticles investigated in this study and tracers (results not shown). This means that the destination maps closely follow the streamline distribution from the trachea to the outlets. Quantitative differences, however, were observed between the FTLE values along the ILCS (inertial p. 12
re-
pro
of
Journal Pre-proof
urn al P
Figure 6: (a) A 3D volume-rendered view of the FTLE field is shown. The red surfaces correspond to the ILCS for dp = 1 µm. (b) The region of interest where 3D FTLE calculation is performed is illustrated. (c) 2D FTLE slices at the trachea (top) and right below the carina (bottom) are shown. (d) Iso volumes of FTLE field with a value of 3.5 are displayed at the same locations. An integration time of T = 0.24s is used in this figure.
Destination map
Source map
tref = 0.2s
tref = 0.4s
tref = 0.6s
tref = 1.0s
Jo
Figure 7: The first row is RML destination map obtained from post-processing of traditional aerosol bolus tracking released at tref = 0.2, 0.4, and 1.0s time-points of inhalation. The second row displays the source map with the corresponding injection times by computing backward MR integral (source inversion). particles) and LCS (ideal tracers). The repelling ILCS mark the boundaries of the flow with high stretching. Particles originating at the different sides of these structures, by definition, separate exponentially fast. In the upper airways, this corresponds to bifurcating vessels that create such separation. The Lagrangian coherent structures are invariant manifolds where particles at the LCS surfaces remain on these structures
p. 13
Journal Pre-proof
of
Table 1: The number of particles at the trachea obtained from backward integration of particles seeded at RML (second column) and the number of particles starting at the source map that reach the RML from solving the full MR equation forward in time (third column) are listed. The success rate is the ratio of the third column to the second column and quantifies the accuracy of the source inversion process. The rows are associated with various injection times. tref = 0s refers to the inhalation starting point.
0.0 0.2 0.4 0.6 1.0 1.2
No. source particles at the trachea from RML (backward integration) 2904 64500 94316 104949 120884 134790
No. particles reaching RML from source map (forward integration) 2659 59595 89701 102351 119114 133408
re-
tref (s)
pro
Source map verification
Success rate percentage 91.6% 93.1% 95.1% 97.5% 98.5% 99.0%
Jo
urn al P
while being advected with the flow [42]. These observations, collectively, imply that particles originating from the ILCS should deposit at the bifurcation regions. Indeed, our results confirmed this conjecture where particles originating at the ILCS ended up depositing at the airway wall (Fig. 5). Therefore, the FTLE and ILCS not only determine the final particle distribution but also accurately identify the origin sites of particles that end up depositing on the bronchial walls. Interestingly, prior studies have shown that release positions of particles at the trachea depositing at the airway wall resemble 2D lines (with structures very similar to our ILCS patterns) [29]. Our numerical results are also consistent with the recent experiments of Amili et al. [2]. These authors showed how, when inertial particles are injected in a successively bifurcating vessel, their fate is strongly dependent on the release location and weakly affected by the flow rate through the distal branches. We also introduced a novel method to create source map with backward in time integration of the MR equation. We presented an example of how the MR equation could be integrated backward in time using the slow manifold technique to overcome the numerically ill-conditioned backward MR equation. The results demonstrate that the source map fully covers the destination map for the target lobe for different reference injection times. In addition to the excellent qualitative agreement, more than 90% of the particles that were integrated backward in time reached the expected location of release. The slow manifold technique is an approximation to the backward MR equation, and therefore perfect accuracy is not expected [47]. Nevertheless, this approach presents an excellent method for identification of the source of any desired particle. For example, instead of launching millions of particles at the trachea to identify particles passing through a specific airway vessel, particles could be launched only at the desired location and integrated backward in time. This could improve the spatial resolution of the destination map by significantly reducing the computational cost. The methods presented in this study have important physiological implications. Inhaled aerosols are commonly used in respiratory disease drug delivery [37]. It is imperative to increase the dosage of the drug reaching the diseased site to increase drug delivery efficacy and reduce harmful side p. 14
Journal Pre-proof
urn al P
re-
pro
of
effects [18, 53]. Construction of destination maps at the trachea to guide spatial and temporal drug injection strategy is the first step to improve efficacy [57, 49]. In this study, we demonstrated the important role that FTLE and ILCS play in these maps. Namely, our novel topological analysis could be used to distinguish between the boundaries of the destination map corresponding to different generations of airway bifurcations. Additionally, the degree of the temporal evolution of ILCS highlights the sensitivity of drug transport results with respect to injection time. The ILCS surfaces can potentially appear, disappear, or move with the flow during inhalation, which highlight the importance of an appropriate temporal release strategy. Finally, source inversion using the slow manifold method could be used to identify drug release location with minimal computational effort. The present study has a few limitations. Particle-particle interactions were ignored in our study. The inlet velocity was assumed to be parabolic although the laryngeal jet can create asymmetric flow profiles. Information from the upstream flow could be considered in future studies. The lower generations of tracheabronchial trees move relatively little, which was not modeled. We have not modeled potential cycle-to-cycle variations in respiratory flow patterns. This work is limited to spherical particles, whereas other types of particles could exhibit different trajectory and results. Our study is limited in its parameter space (just one Re and one Wo number). This limitation is modulated by the fact that this study is aimed to demonstrate the value of the tools, rather than the effect of varying the flow regimes. The results are sensitive to the lung geometry and boundary conditions, thus we expect our results to be different from patient to patient, however, the correlation between FTLE and destination maps is expected to preserve. Finally, the pharmaceutic particles are normally administered by the oral route in clinical practice. Therefore, it would be more practical if destination maps were found at the oronasal cavity. In conclusion, we have presented novel topological tools from dynamical systems theory, which can provide valuable insight into inertial particle transport patterns in lung airways. Our proposed tools assist in coming up with a more accurate spatial and temporal drug release strategy. Future work should extend our framework to the oronasal cavity where more complex and turbulent flow is present. Inertial particle transport in complex, chaotic flows can lead to rich and unique physics [10] (e.g. transient chaos [28]) that has received very little attention in respiratory flows [49].
Conflict of interest None.
Acknowledgement
Jo
The authors would like to thank Kamran Poorbahrami and Jessica Oakes for fruitful discussions related to particle dynamics in lung airways and Milad Habibi for discussions related to slow manifolds.
References
[1] M. R. Allshouse and T. Peacock. Lagrangian based methods for coherent structure detection. Chaos: An Interdisciplinary Journal of Nonlinear Science, 25(9):097617, 2015. [2] O. Amili, J. Golzarian, and F. Coletti. In vitro study of particle transport in successively bifurcating vessels. Annals of Biomedical Engineering, pages 1–13, 2019. p. 15
Journal Pre-proof
[3] A. Arzani, A. S. Les, R. L. Dalman, and S. C. Shadden. Effect of exercise on patient specific abdominal aortic aneurysm flow topology and mixing. International Journal for Numerical Methods in Biomedical Engineering, 30(2):280–295, 2014. [4] A. Arzani and S. C. Shadden. Characterization of the transport topology in patient-specific abdominal aortic aneurysm models. Physics of Fluids, 24(8):081901, 2012.
pro
of
[5] B. Asgharian, O. T. Price, M. Oldham, L. Chen, E. L. Saunders, T. Gordon, V. B. Mikheev, K. R. Minard, and J. G. Teeguarden. Computational modeling of nanoscale and microscale particle deposition, retention and dosimetry in the mouse respiratory tract. Inhalation Toxicology, 26(14):829–842, 2014. [6] A. J. Banko, F. Coletti, C. J. Elkins, and J. K. Eaton. Oscillatory flow in the human airways from the mouth through several bronchial generations. International Journal of Heat and Fluid Flow, 61:45–57, 2016.
re-
[7] W. D Bennett, M. S. Messina, and G. C. Smaldone. Effect of exercise on deposition and subsequent retention of inhaled particles. Journal of Applied Physiology, 59(4):1046–1054, 1985. [8] G. Boffetta, G. Lacorata, G. Redaelli, and A. Vulpiani. Detecting barriers to transport: a review of different techniques. Physica D: Nonlinear Phenomena, 159(1-2):58–70, 2001.
urn al P
[9] A. E BozorgMagham, S. D. Ross, and D. G. Schmale III. Real-time prediction of atmospheric Lagrangian coherent structures based on forecast data: An application and error analysis. Physica D: Nonlinear Phenomena, 258:47–60, 2013. [10] J. H. E. Cartwright, U. Feudel, G. K´arolyi, A. De Moura, O. Piro, and T. T´el. Dynamics of finite-size particles in chaotic fluid flows. In Nonlinear Dynamics and Chaos: Advances and Perspectives, pages 51–87. Springer-Verlag Berlin Heidelberg, 2010. [11] Y. Cheng, Y. Zhou, and B. T. Chen. Particle deposition in a cast of human oral airways. Aerosol Science & Technology, 31(4):286–300, 1999. [12] C. C. Daigle, D. C. Chalupa, F. R. Gibb, P. E. Morrow, G. Oberd¨orster, M. J. Utell, and M. W. Frampton. Ultrafine particle deposition in humans during rest and exercise. Inhalation Toxicology, 15(6):539–552, 2003.
Jo
[13] P. Das, E. Nof, I. Amirav, S. C. Kassinos, and J. Sznitman. Targeting inhaled aerosol delivery to upper airways in children: Insight from computational fluid dynamics (CFD). PloS One, 13(11):e0207711, 2018. [14] A. Farghadan and A. Arzani. The combined effect of wall shear stress topology and magnitude on cardiovascular mass transport. International Journal of Heat and Mass Transfer, 131:252– 260, 2019. [15] G. Haller. A variational theory of hyperbolic Lagrangian coherent structures. Physica D: Nonlinear Phenomena, 240(7):574–598, 2011. [16] G. Haller. Lagrangian coherent structures. Annual Review of Fluid Mechanics, 47:137–162, 2015. p. 16
Journal Pre-proof
[17] G. Haller and T. Sapsis. Where do inertial particles go in fluid flows? Physica D: Nonlinear Phenomena, 237(5):573–583, 2008. [18] D. R. Hess, P. Anderson, R. Dhand, J. L. Rau, G. C. Smaldone, and G. Guyatt. Device selection and outcomes of aerosol therapy: Evidence-based guidelines. Chest, 127(335):335–371, 2005.
of
[19] T. Higenbottam, T. Siddons, and E. Demoncheaux. The direct and indirect action of inhaled agents on the lung and its circulation: lessons for clinical science. Environmental Health Perspectives, 109(suppl 4):559–562, 2001.
pro
[20] W. Hofmann. Modelling inhaled particle deposition in the human lung–a review. Journal of Aerosol Science, 42(10):693–724, 2011. [21] H. Homann and J. Bec. Finite-size effects in the dynamics of neutrally buoyant particles in turbulent flow. Journal of Fluid Mechanics, 651:81–91, 2010.
re-
[22] S. Jalal. Velocity & Vorticity Transport In 3D-Printed Idealized & Realistic Human Airways Using Magnetic Resonance Velocimetry (MRV) & Particle Image Velocimetry (PIV). PhD thesis, University of Minnesota, 2019. [23] C. Kleinstreuer and Z. Zhang. Airflow and particle transport in the human respiratory system. Annual Review of Fluid Mechanics, 42:301–334, 2010.
urn al P
[24] C. Kleinstreuer, Z. Zhang, and J. F. Donohue. Targeted drug-aerosol delivery in the human respiratory system. Annu. Rev. Biomed. Eng., 10:195–220, 2008. [25] A. V. Kolanjiyil and C. Kleinstreuer. Computationally efficient analysis of particle transport and deposition in a human whole-lung-airway model. part I: Theory and model validation. Computers in Biology and Medicine, 79:193–204, 2016. [26] A. V. Kolanjiyil, C. Kleinstreuer, and R. T. Sadikot. Computationally efficient analysis of particle transport and deposition in a human whole-lung-airway model. part II: Dry powder inhaler application. Computers in Biology and Medicine, 84:247–253, 2017. [27] P. G. Koullapis, L. Nicolaou, and S. C. Kassinos. In silico assessment of mouth-throat effects on regional deposition in the upper tracheobronchial airways. Journal of Aerosol Science, 117:164–188, 2018.
Jo
[28] T. Lai, Y. C.and T´el. Transient chaos: complex dynamics on finite time scales, volume 173. Springer New York, 2011. [29] Z. Li, C. Kleinstreuer, and Z. Zhang. Particle deposition in the human tracheobronchial airways due to transient inspiratory flow patterns. Journal of Aerosol Science, 38(6):625–644, 2007. [30] Z. Li, C. Kleinstreuer, and Z. Zhang. Simulation of airflow fields and microparticle deposition in realistic human lung airway models. part I: Airflow patterns. European Journal of MechanicsB/Fluids, 26(5):632–649, 2007. [31] Z. Li, C. Kleinstreuer, and Z. Zhang. Simulation of airflow fields and microparticle deposition in realistic human lung airway models. part II: Particle transport and deposition. European Journal of Mechanics-B/Fluids, 26(5):650–668, 2007.
p. 17
Journal Pre-proof
[32] P. W. Longest, G. Tian, R. L. Walenga, and M. Hindle. Comparing MDI and DPI aerosol deposition using in vitro experiments and a new stochastic individual path (SIP) model of the conducting airways. Pharmaceutical Research, 29(6):1670–1688, 2012. [33] M. R. Maxey and J. J. Riley. Equation of motion for a small rigid sphere in a nonuniform flow. The Physics of Fluids, 26(4):883–889, 1983.
of
[34] M. Mortensen and K. Valen-Sendstad. Oasis: A high-level/high-performance open source navier–stokes solver. Computer Physics Communications, 188:177–188, 2015.
pro
[35] A. Naseri, O. Abouali, P. Ghalati Farhadi, and G. Ahmadi. Numerical investigation of regional particle deposition in the upper airway of a standing male mannequin in calm air surroundings. Computers in Biology and Medicine, 52:73–81, 2014. [36] J. M. Oakes, A. L. Marsden, C. Grandmont, S. C. Shadden, C. Darquenne, and I. E. VignonClementel. Airflow and particle deposition simulations in health and emphysema: from in vivo to in silico animal experiments. Annals of Biomedical Engineering, 42(4):899–914, 2014.
re-
[37] U. Pison, T. Welte, M. Giersig, and D. A. Groneberg. Nanomedicine for respiratory diseases. European Journal of Pharmacology, 533(1-3):341–350, 2006. [38] K. Poorbahrami and J. M. Oakes. Regional flow and deposition variability in adult female lungs: A numerical simulation pilot study. Clinical Biomechanics, 66:40–49, 2019.
urn al P
[39] S. Qi, B. Zhang, Y. Teng, J. Li, Y. Yue, Y. Kang, and W. Qian. Transient dynamics simulation of airflow in a ct-scanned human airway tree: more or fewer terminal bronchi? Computational and Mathematical Methods in Medicine, 2017:1–14, 2017. [40] T. Sapsis and G. Haller. Inertial particle dynamics in a hurricane. Journal of the Atmospheric Sciences, 66(8):2481–2492, 2009. [41] S. C. Shadden and A. Arzani. Lagrangian postprocessing of computational hemodynamics. Annals of Biomedical Engineering, 43(1):41–58, 2015. [42] S. C. Shadden, F. Lekien, and J. E. Marsden. Definition and properties of Lagrangian coherent structures from finite-time Lyapunov exponents in two-dimensional aperiodic flows. Physica D: Nonlinear Phenomena, 212(3-4):271–304, 2005.
Jo
[43] S. C. Shadden and C. A. Taylor. Characterization of coherent structures in the cardiovascular system. Annals of Biomedical Engineering, 36(7):1152–1162, 2008. [44] B. Soni and S. Aliabadi. Large-scale CFD simulations of airflow and particle deposition in lung airway. Computers & Fluids, 88:804–812, 2013. [45] B. Soni, D. Thompson, and R. Machiraju. Visualizing particle/flow structure interactions in the small bronchial tubes. IEEE Transactions on Visualization and Computer Graphics, 14(6):1412–1427, 2008. [46] M. Sudharsan, S. L. Brunton, and J. J. Riley. Lagrangian coherent structures and inertial particle dynamics. Physical Review E, 93(3):033108, 2016. [47] W. Tang, G. Haller, J. Baik, and Y. Ryu. Locating an atmospheric contamination source using slow manifolds. Physics of Fluids, 21(4):043302, 2009. p. 18
Journal Pre-proof
[48] K. Tor´en, I. A. Bergdahl, T. Nilsson, and B. J¨arvholm. Occupational exposure to particulate air pollution and mortality due to ischaemic heart disease and cerebrovascular disease. Occupational and Environmental Medicine, 64(8):515–519, 2007. [49] A. Tsuda, F. S. Henry, and J. P. Butler. Particle transport and deposition: basic physics of particle kinetics. Comprehensive Physiology, 3(4):1437–1471, 2011.
of
[50] A. Updegrove, N. M. Wilson, J. Merkow, H. Lan, A. L. Marsden, and S. C. Shadden. Simvascular: An open source pipeline for cardiovascular simulation. Annals of Biomedical Engineering, 45(3):525–541, 2017.
pro
[51] T. Van de Moortele, U. Goerke, Chris H. Wendt, and F. Coletti. Airway morphology and inspiratory flow features in the early stages of chronic obstructive pulmonary disease. Clinical Biomechanics, 66:60–65, 2017.
re-
[52] T. Van de Moortele, C. H. Wendt, and F. Coletti. Morphological and functional properties of the conducting human airways investigated by in vivo computed tomography and in vitro mri. Journal of Applied Physiology, 124(2):400–413, 2017. [53] J. Vestbo, S. S. Hurd, A. G. Agust´ı, P. W. Jones, et al. Global strategy for the diagnosis, management, and prevention of chronic obstructive pulmonary disease: GOLD executive summary. American Journal of Respiratory and Critical Care Medicine, 187(4):347–365, 2013.
urn al P
[54] J. Xi, X. A. Si, J. W. Kim, E. Mckee, and E. B. Lin. Exhaled aerosol pattern discloses lung structural abnormality: a sensitivity study using computational modeling and fractal analysis. PloS One, 9(8):e104682, 2014. [55] X. Xu, Y. Shang, L. Tian, W. Weng, and J. Tu. Inhalation health risk assessment for the human tracheobronchial tree under pm exposure in a bus stop scene. Aerosol and Air Quality Research, 19(6):1365–1376, 2019. [56] B. Zhang, S. Qi, Y. Yue, J. Shen, C. Li, W. Qian, and J. Wu. Particle disposition in the realistic airway tree models of subjects with tracheal bronchus and COPD. BioMed Research International, 2018:1–15, 2018. [57] Z. Zhang and C. Kleinstreuer. Effect of particle inlet distributions on deposition in a triple bifurcation lung airway model. Journal of Aerosol Medicine, 14(1):13–29, 2001.
Jo
[58] Z. Zhang, C. Kleinstreuer, J. F. Donohue, and C. S. Kim. Comparison of micro-and nano-size particle depositions in a human upper airway model. Journal of Aerosol Science, 36(2):211–233, 2005. [59] Z. Zhang, C. Kleinstreuer, and C. S. Kim. Aerosol deposition efficiencies and upstream release positions for different inhalation modes in an upper bronchial airway model. Aerosol Science & Technology, 36(7):828–844, 2002.
p. 19
PII: DOI: Reference:
S0010-4825(19)30362-2 https://doi.org/10.1016/j.compbiomed.2019.103497 CBM 103497
To appear in:
Computers in Biology and Medicine
Received date : 6 September 2019 Revised date : 9 October 2019 Accepted date : 10 October 2019 Please cite this article as: A. Farghadan, F. Coletti and A. Arzani, Topological analysis of particle transport in lung airways: Predicting particle source and destination, Computers in Biology and Medicine (2019), doi: https://doi.org/10.1016/j.compbiomed.2019.103497. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Ltd.
Journal Pre-proof
The author hereby declares no conflicts of interest.
Jo
urn al P
re-
pro
Amirhossein Arzani Assistant Professor Cardiovascular Biomechanics Lab https://www2.nau.edu/cardio-lab/ Mechanical Engineering Department Northern Arizona University
of
Conflict of Interest Statement
0
Jo
FTLE
6
33
map
Destination
Graphical Abstract (for review)
urn al P
re-
pro
of
Journal Pre-proof
1
tL = 0.0s
tL = 0.2s
tL = 0.4s
Color-coded outlets
Journal Pre-proof *Revised manuscript (clean) Click here to download Revised manuscript (clean): Farghadanetal19_revise_notred.pdf
Filippo Coletti2
Amirhossein Arzani1
pro
Ali Farghadan1
of
Topological analysis of particle transport in lung airways: predicting particle source and destination
1
Jo
urn al P
Correspondence: Amirhossein Arzani, Northern Arizona University, Flagstaff, AZ, 86011 Phone: +1-9285230541 Email: [email protected]
re-
Department of Mechanical Engineering, Northern Arizona University, Flagstaff, AZ, United States 2 Department of Aerospace Engineering and Mechanics, University of Minnesota, Minneapolis, MN, United States
1
Journal Pre-proof
Abstract
Jo
urn al P
re-
pro
of
Particle transport in lung airways can induce respiratory disease and play a vital role in aerosol drug delivery. Herein, we present dynamical systems features that influence airflow and particle transport in the tracheobronchial trees. Computational fluid dynamics (CFD) was used to solve for unsteady airflow in a patient-specific model. Particle tracking simulations were performed for micron-size particles. The destination map that connects the particle final location to the initial location and injection time was constructed. Finite-time Lyapunov exponent (FTLE) fields were calculated to identify inertial Lagrangian coherent structures (ILCS), topological features that act as separatrices. Our results demonstrated that these topological features control the destination map at the trachea. The temporal evolution of ILCS influenced the sensitivity of particle transport fate to injection time, whereas the emergence of new ILCS with an increased integration time controlled transport to different generations of airways. Additionally, particles starting at the ILCS were shown to mostly deposit at the airway walls. Finally, an innovative source inversion strategy was introduced to integrate the MaxeyRiley equation backward in time and identify the origin of dispersed particles. Our study explores novel dynamical systems tools that improve our understanding of particle transport and deposition in the airways and could be used to guide future targeted drug delivery studies. Keywords: airflow, computational fluid dynamics, FTLE, Lagrangian coherent structures, source inversion, particle tracking, aerosol.
Journal Pre-proof
1
Introduction
Jo
urn al P
re-
pro
of
The essence of particle transport in the respiratory system is well known. Harmful substances such as dust, fumes, and gaseous pollutants are vastly suspended in the environment with a variety of shapes and sizes ranging from ultrafine (smaller than 1 µm) to large-size matters (larger than 100 µm) [5]. Detrimental substances have either immediate effect and cause minor disorders such as dizziness or can cause serious long-term damage to the lungs including asthma or cancer [19, 48]. In the case of many of the developed chronic lung diseases, prescribing an effective treatment is necessary. Overall, for both medical and pathological purposes, studying particle transport provides diagnostic and prophylactic measures for adverse respiratory health conditions. To study particle transport, computational fluid dynamics (CFD) coupled with particle dynamics simulations are widely employed [23, 35, 25, 56, 38]. The dichotomous network of airways represents one of the complicated dynamical systems in the human body where chaotic flow and advection occurs [49]. Several characteristics contribute to complex airflow patterns such as lung geometry with numerous planar and non-planar conducting tubes, unsteady nature of the breathing cycle, asymmetric tapering down into distal branches, and other physiological factors [30, 6, 51]. All of these features are known to influence particle transport and deposition patterns [7, 12, 31, 20]. In practice, Lagrangian approach and in particular solving the Maxey-Riley (MR) equation [33], which governs the dynamics of small spherical inertial particles, has been shown to be an effective method to simulate particle transport in lung airways [32, 36, 27]. This equation is extensively studied in a broad range of applications for all types of particles such as aerosols, neutrally buoyant particles, and bubbles [21, 46]. In the context of aerosol transport in lung airways, this equation is used to answer two important questions: First, what is the fate of each particle entering the upper airways? Second, what spatial spot and temporal time-point should particles be released to reach a specific desired location? The first question is of interest regarding harmful inhaled particles, while the latter question is one of the major problems involved in targeted drug delivery. To answer the aforementioned questions, particle destination maps have been proposed [24, 54, 26]. Namely, particles are densely seeded at the inlet of an airway model and integrated based on the MR equation. Subsequently, the destination of each particle (typically based on the outlets of the model) is identified and mapped back to the initial location and time-point of release. Such destination maps are used to study the effect of spatial location [59] and intra-cycle time-point of release [13] on particle transport and fate. Additionally, the effect of particle diameter (or more generally the Stokes number) on destination maps has been studied [58]. While destination maps are useful in identifying the fate of injected (or incoming) particles, they do not provide a mechanistic explanation of particle transport results and are often difficult to interpret. Moreover, if one is interested in targeted drug delivery to a specific location of the lung airways, a large number of particles are often required to be released and tracked such that the localized region of interest is sufficiently sampled. In this study, we propose novel dynamical systems tools to overcome these issues and improve our understanding of particle transport in the conducting airways. First, we perform a topological analysis of fluid flow and particle transport using finite-time Lyapunov exponent (FTLE) fields and Lagrangian coherent structures (LCS) [42]. LCS mark boundaries of flow with distinct dynamical behavior and could be considered as invisible surfaces embedded in fluid flow that act as transport barriers. Therefore, they provide a template to study transient and chaotic flows. Herein, we compute the extension of this concept to inertial particles (inertial LCS) [40, 46]. Our group and others have used FTLE and LCS to characterize transport in cardiovascular flow applications [43, 4, p. 3
Journal Pre-proof
2.1
Methods Airflow simulations
urn al P
2
re-
pro
of
3, 41, 14]. In respiratory flows, only one group has studied the significance of FTLE in the airways and illustrated the correlation between FTLE and destination maps [45, 44]. However, the effect of temporal FTLE variations and integration time on particle distribution as well as the effect of LCS on particle deposition was not revealed. In this study, we provide further insight into how FTLE maps could fully characterize particle transport in the tracheobronchial tree. Finally, we introduce a novel source-inversion technique to identify the initial position of dispersed particles with very little computational effort. Namely, we use the recent slow manifold technique developed by Haller and Sapsis [17] to integrate the MR equation backward in time and identify the initial location of release for any desired set of particles. The manuscript is organized as follows. First, we describe our CFD and particle dynamics frameworks for a patient-specific model under unsteady flow conditions. Next, we explain FTLE and inertial LCS calculation as well as source inversion using the slow manifold method. We discuss how these techniques could be used to study the spatiotemporal distribution of particle destination maps. We demonstrate how the integration time in FTLE calculation and the resulting LCS correlate with particle distribution to different generations of bifurcating airways. Finally, we demonstrate how backward integration of the MR equation could be used to identify the particle injection location, which is of high interest in targeted drug delivery.
The inhalation and exhalation simulation is performed on the same model employed in an invitro study [22]. The inlet of the model is extended at the trachea. A parabolic, Newtonian g ), incompressible (ρf = 1.2 × 10−3 cmg 3 ), unsteady airflow is imposed at (µf = 1.81 × 10−4 cm.s the inlet. The airflow is assumed to be sinusoidal (Q = Qmax Sin(2πf qt)) with peak Reynolds ρf V¯max din ) of 1500 and Womersley number (W o = d2in 2πf ) of 2.6, where νf is number (Remax = µf νf
Jo
air kinematic viscosity, din = 1.6 cm is the inlet diameter and f = 15 respiration is the respiratory min frequency where each breathing cycle takes 4 seconds. The corresponding cross-section averaged L peak velocity is V¯max = 141.4 cm and peak flow rate is Qmax = 17.1 min . The nondimensional s Remax and Wo numbers belong to a normal breathing rate in an adult. The flow is assumed laminar throughout the tracheobronchial tree. No-slip boundary condition is prescribed at the wall. For the inhalation, resistance boundary conditions are imposed at the outlets and are tuned iteratively such that the flow rate distributions are 14, 5, 32, 17, 32 percent in right upper lobe (RUL), right middle lobe (RML), right lower lobe (RLL), left upper lobe (LUL), and left lower lobe (LLL), respectively. These values follow the lobar flow rate distributions obtained by multidetector computed tomography in a cohort of healthy non-smokers [52]. For the exhalation, a plug velocity profile is assigned to each terminal branch matching the unsteady exhalation flow rate. The simulation time step for the whole reparation cycle is 2 × 10−4 s and the results are found to be time-step independent with this selection of time-step. An overview of the 3D model and boundary conditions are shown in Fig. 1. The mesh consisted of 5.3M quadratic tetrahedral elements with five boundary layers created with SimVascular [50]. The boundary layer meshing is used to capture the near-wall flow phenomena. This second-order mesh is equivalent to a first-order mesh with ∼42M linear tetrahedral elements, resulting in an extremely fine spatial resolution compared to other highresolution studies in the literature (mostly less than 15M elements [39, 38, 55]). Unsteady velocity data are obtained using Oasis [34], an open-source Navier-Stokes solver (based on the finite element p. 4
Journal Pre-proof
method) with minimum numerical dissipation. The second respiratory cycle (t = 4 to 8s) is saved for post-processing. The focus of this study is on the inhalation phase and the exhalation phase simulation is performed to provide more realistic flow conditions at the beginning of inhalation.
Exhalation Phase
Vinh(t)
PTrachea
Trachea seeded surface
PRUL RRUL
RLUL
17% Vexh(t) RRML
PLLL
14% Vexh(t)
pro
PLUL
32% Vexh(t)
re-
PRLL
32% Vexh(t)
(b)
urn al P
(a)
RML seeded surface
5% Vexh(t)
PRML
RRLL
RLLL
Seeded surfaces
of
Inhalation Phase
(c)
Figure 1: The patient-specific lung airway geometry used in this study is displayed. The CFD boundary conditions for (a) inhalation and (b) exhalation phases of the breathing cycle are shown. (c) The location of a thin box at the trachea where particles are launched and FTLE fields are calculated is shown. Also, the location of seeded particles at the entrance of right middle lobe (RML) where source inversion study is performed is displayed.
2.2
Particle tracking
Particles are tracked within the domain of interest throughout the inhalation phase. The governing equation of spherical particle-laden flows (the MR equation) is integrated in a Lagrangian approach
Jo
dv(x(t)) 3 Du(x, t)) 18µf 1 = (ρp − ρf )g + ρf − 2 (v(x(t)) − u(x, t))) , (ρp + ρf ) 2 dt 2 Dt dp
(1)
where g is the gravity, facing downward resembling an upright standing, v and u are particle and air velocities, respectively. µf is fluid viscosity, and ρp and ρf are particle and air densities, respectively. Stokesian drag is assumed, which is justified by the very small Reynolds number based on the particle diameter and particle velocity relative to the surrounding fluid (Rep ∼ O(10−3 )). Particle-particle interactions and back-reaction on the carrier fluid are ignored. Faxen correction and Basset history terms are also neglected. We also ensured that Brownian motion was ignorable (mainly due to the domination of advection over diffusion) for microscale particle transport in the tracheobronchial tree. Particles are seeded at a trachea cross-section right before the model extension with a velocity equal to the fluid velocity at releasing location as shown in Fig. 1c. Particles are launched at different time-points and uniformly from all spots (∼770,000 location points) on the trachea surface p. 5
Journal Pre-proof
Destination maps
pro
2.2.1
of
(the upper green slice in Fig. 1c). Two extreme limits of microparticles typically entering the tracheobronchial tree are simulated, assuming particle diameter of dp = 1 and 10 µm. The bulk portion of large particles (over 10 µm) cannot reach the trachea [11], and nanoscale particles are ρ d2 w not considered in this study. The Stokes numbers (St = 18µp fpdin , w is the mean inhalation velocity at the trachea) at the inlet corresponding to these diameters are St = 1.8 × 10−4 and 1.8 × 10−2 , respectively. Particles are integrated until they either leave the computational domain, deposit on the bronchial wall (the distance between particle center and wall is less than particle radius), or the inhalation cycle is completed.
2.3
re-
Destination maps are a useful tool that give a better understanding of particle fate. To obtain the destination map, first, particle tracking is performed as explained above. Each outlet is tagged with an index (n = 1, 2, ...). For each particle tracked, the outlet that it exits the domain is determined. Subsequently, the outlet index is assigned to the initial release location of the particle at the trachea. The resulting map of indices at the trachea is called the destination map. All deposited particles are assigned a zero index (n = 0).
Finite-time Lyapunov exponent (FTLE) and inertial Lagrangian coherent structures (ILCS)
urn al P
Lagrangian coherent structures (LCS) are invisible material surfaces embedded in 3D flows (lines in 2D flows) that have distinct attracting or repelling properties [16]. The attracting LCS attract nearby particles in their basin of attraction and could be sometimes visualized in the lab using experimental visualization techniques. On the other hand, the repelling LCS repel nearby trajectories and act as transport barriers. In this study, we only compute the repelling LCS and demonstrate their influence on particle distribution. While mathematically more accurate techniques have been proposed for identification of LCS [15], the finite-time Lyapunov exponent (FTLE) field [42] remains a computationally robust approach in identifying LCS. FTLE is utilized to characterize coherent structures in a wide range of studies [8, 9, 1, 41]. FTLE measures the exponential rate of stretching of initially closely seeded particles. The ridges of the FTLE field (loosely speaking, local maxima) identify the repelling or attracting LCS under forward and backward time integration, respectively. To calculate the FTLE field, a densely distributed set of particles is seeded in a region of interest and advected via integration of the MR equation. The FTLE field is subsequently computed as [42]
Jo
σ(x0 ; t0 , t) =
1 p ln λmax C(x0 ; t0 , t) , T
(2) |
where λmax C(x0 ; t0 , t) is the maximum eigenvalue of C(x0 ; t0 , t), C(x0 ; t0 , t) = ∇Ftt0 · ∇Ftt0 is the right Cauchy-Green strain tensor, Ftt0 : x(t0 ) 7→ x(t) is the flow map (maps the initial particle location to the final location), T = t − t0 is the integration time, and the FTLE field is calculated as a scale of the largest eigenvalue of the Cauchy-Green strain tensor. The FTLE field (σ(x0 ; t0 , t)) is a scalar field and in unsteady flows it can vary depending on the initial time-point of particle release (t0 ). The division over integration time (T ) in Eq. 2 could be dropped to avoid sensitivity of FTLE values to particles leaving the computational domain early [4]. The FTLE field is computed for a thin box in the same location where particles are initially seeded. The calculation is only confined to a thin box at the inlet because we are interested in the trajectory of particles released at this region. Over 1.5M particles with diameters of 1 and 10 µm p. 6
Journal Pre-proof
2.4
pro
of
are launched in six time-points of the inhalation phase. Forward time integration is performed, and therefore the ridges of the FTLE field identify the repelling LCS. Additionally, since we are tracking inertial particles (instead of ideal tracers) the inertial Lagrangian coherent structures (ILCS) [40, 46] are identified. The main difference between ILCS and LCS is in the flow maps. Constructing flow maps by solving the MR equation identifies the ILCS, whereas the LCS are identified by integrating ideal tracers along velocity pathlines. To gain a better understanding of coherent structures inside the lung, the FTLE field is calculated for a 3D box spanning the trachea to the carina region. Calculating FTLE for a large 3D box is computationally costly since millions of particles (in this example, 27M particles densely seeded in the 3D domain) need to be integrated by solving the MR equation.
Source inversion
re-
The Maxey-Riley (MR) equation, mathematically, could be integrated forward and backward in time. However, it has been proven that backward integration of the MR equation for locating the origin of dispersed particles (source inversion) is quite challenging. In fact, for large-scale timedependent flows, even a high-precision integration scheme with very fine time-step blows up when the MR equation is integrated backward in time. To overcome this instability, Haller and Sapsis [17] proposed the idea of slow manifolds under unsteady flow conditions, which is briefly explained here. We can rewrite the MR equation (Eq. 1) by dividing both sides by the coefficient of particle acceleration (ρp + 12 ρf ) as
18µ
urn al P
Du(x, t)) dv(x(t)) = (1 − 1.5R)g + 1.5R − KR(v(x(t)) − u(x, t))) , dt Dt
(3)
2ρ
Jo
f where K = d2 f and R = ρf +2ρ . This equation is singular at the limit of infinitesimally small p p particle diameters. Hence, integrating the MR equation backward in time is numerically impossible when KR >> 1. The term −KR(v(x(t)) in equation 3 causes a quick and inevitable blowup of the numerical solution to this equation backward in time. Using high-order numerical integration schemes or smaller time step would only delay the blowup process, making backward MR integration impractical for large-scale time-dependent problems. Time-dependent invariant slow manifolds (M ), which are three-dimensional global attractors, are proven to exist in the six-dimensional phase space of the MR equation [17]. M can be written as the following form: Du(x, t)) 2 − g + O( ) , (4) M = (x, v) : v = u + (1.5R − 1) Dt
ρ d
where = ρfp win St quantifies the significance of inertia. Note that in the original paper [17], = St R R was written in non-dimensional form and a different definition of St number was used. In order to recover the starting position of an inertial particle without struggling with the instability of solving the full MR equation backward in time, one must first project the final location of a dispersed particle onto the slow manifold along its velocity direction. Then the governing particle motion equation on slow manifolds should be integrated from the final time to the initial time [17]. A first-order approximation to the particle motion on slow manifolds is given by Du(x, t)) x˙ = u(x, t) + (1.5R − 1) − g + O(2 ) , (5) Dt p. 7
Journal Pre-proof
3
pro
of
which is the lowest-order truncation without numerical divergence. The zero limit of makes Equation 5 equivalent to the governing equation of ideal tracers. Sufficiently large values of , however, might lead to inaccuracy. For microscale particle tracking in lungs << 1, and Equation 5 is easily integrable backward in time with good accuracy. Using this method, the initial positions of dispersed particles can be recovered with an error of O(2 ). In this study, as an example of the utility of this method in lung airways, the source of particles reaching the right-middle lobe (RML) have been investigated. The same approach can be utilized for any other lobes or branches. Figure 1c shows the surface of interest at the RML entrance. Particles were released uniformly from over 62,500 spots seeded on this surface. They were injected 380 times during the inhalation cycle (uniformly from 2.0 to 0.1s), thus overall, ∼23.75M particles were tracked and integrated backward with a time step of 5 µs.
Results
3.1
FTLE maps
re-
The corresponding FTLE, destination, and source inversion maps are presented in the following sections.
3.1.1
urn al P
In this section, variations in the FTLE maps at the trachea are found for various cases such as different integration times, several particle injection time-points, and two particle diameters. The results are compared to the destination maps. Effect of temporal integration and particle diameter
3.1.2
Jo
In Figure 2, FTLE plots associated with 1 and 10 µm diameter particles are displayed. Particles are released at the start of inhalation and FTLE fields are calculated for several integration times (T). As the integration time increases, more ILCS emerge corresponding to more distal branches. That is, additional ILCS are created with an increase in integration time, since particles have more time to pass through new generations of bifurcations, and therefore experience separation. The longest integrated FTLE fields (rightmost column) remains qualitative invariant by extending integration time for both diameters. The near-wall region indicates high FTLE values. The reason behind these high values is the high shear rates in these regions that make the initially close particles far apart from each other in time. Moreover, some portion of near-wall particles deposits on the immediate curved wall below the trachea which magnifies these artificially high FTLE values. Finally, the figure shows that FTLE ridges (the ILCS) are robust to variations in particle diameters between 1 and 10 µm. Effect of injection time
FTLE maps are calculated for 1 µm particle diameter launched in different time-points during inhalation. As seen in Figure 3, with an increase in integration time (T), ILCS appear sooner (nearby particles separate sooner) when particles are released later in the breathing cycle. For a fixed value of the integration time, variations in the FTLE field for different launch/injection times (tL ) show the temporal evolution of the FTLE and ILCS, which demonstrate the sensitivity of particle fate to time-point of injection. The FTLE field varies with increasing the injection time up to tL = 0.4s, and is relatively insensitive to further increase in tL . p. 8
pro
of
Journal Pre-proof
3.1.3
Destination maps
re-
Figure 2: The effect of integration time (T) and particle diameter on the FTLE fields is shown. The first and second row correspond to 1 and 10 µm particle diameters, respectively. The same FTLE color bar range is used for both rows. The emergence of new FTLE ridges (ILCS) is shown with an increase in integration time.
3.1.4
urn al P
The particle destination maps are compared with ILCS for 1 µm particle size as shown in Figure 4. The deposited particles are excluded from the destination maps. A very strong correlation between the ILCS and FTLE field and the corresponding destination map is seen. Each region bound inside ILCS lines (with lower FTLE values) corresponds to a different tag in the destination map. A very close match exists between ILCS lines and boundaries lying on the destination maps for the three injection time-points shown (which are representative of the considered values of tL within the inhalation phase). Deposition maps
ILCS surfaces computed from forward FTLE fields are repellers and act as transport barriers. Therefore, virtually no particle is expected to pass from one side of ILCS to the other. The ILCS partition transport to different branches and this leads to a very interesting scenario in which particles right on the ILCS deposit on the wall. As seen in Figure 5, the original location of ultimately deposited particles is right on top of the ILCS. The same scenario occurs for different injection times as shown in the figure. 3D FTLE visualization
Jo
3.1.5
As particles are injected at the trachea, it suffices to calculate FTLE fields for a thin slice at the same location where the corresponding ILCS surfaces are visualized as lines (shown in the above figures). However, FTLE calculation throughout an extended region of the airways could be performed to visualize the 3D ILCS surfaces. This is done in Figure 6 where the 3D FTLE and ILCS patterns are shown. The results are shown for 1 µm particle diameter and 0.24s integration time when the first ILCS emerge at the trachea. One can see how the transport barrier represented by the ridge of high FTLE extends down the trachea, and into the main carina.
p. 9
T = 0.14s
T = 0.24s
T = 0.28s
T = 0.12s
T = 0.14s
T = 0.02s
T = 0.06s
T = 0.08s
0
T = 0.04s
T = 0.06s
T = 010s
T = 0.08s
urn al P
T = 0.02s
T = 0.16s
T = 0.02s
T = 0.02s
T = 0.04s
T = 0.04s
T = 0.20s
re-
FTLE
6 T = 0.02s
T = 0.32s
T = 0.06s
T = 0.06s
T = 0.08s
T = 0.08s
T = 0.12s
T = 0.10s
tL = 0.0s
T = 0.40s
pro
T = 0.02s
of
Journal Pre-proof
tL = 0.2s
T = 0.40s
tL = 0.4s
T = 0.40s
tL = 0.6s T = 0.40s
tL = 1.0s T = 0.10s
T = 0.40s
tL = 1.2s T = 0.10s
T = 0.40s
3.2
Jo
Figure 3: The temporal evolution of FTLE fields for several integration times is shown. For each integration time (T), FTLE fields are calculated for different injection times (tL ) to demonstrate the temporal evolution of FTLE. All calculations are performed for 1 µm particle diameter. The same FTLE color bar range is used for all results.
Source inversion
By launching particles seeded at RML entrance and tracking them backward in time, particles reach the trachea in a range of times. Some particles never find their way back to the trachea as they are seeded in a region that is not sampled by the forward MR equation. That is, no particle injected at the trachea reaches these locations, indicating that such locations are virtually impossible to reach under the present flow conditions and range of injection time-points. In order to create a source map of the target RML at the trachea in a specific time (tref ), only particles that arrive at the trachea within tref ± 5ms range are kept. That is, particles at the RML are integrated backward in time and the union of particles reaching the trachea at the p. 10
Journal Pre-proof
FTLE
6
of
33
pro
map
Destination
0
1
tL = 0.0s
tL = 0.2s
tL = 0.4s
Color-coded outlets
urn al P
re-
Figure 4: The comparison between FTLE and destination maps are shown for dp = 1 µm. The top row shows the FTLE field with T = 0.40s integration time and the bottom row is the corresponding destination maps. The three panels correspond to 0.0, 0.2, and 0.4s injection times from left to right, respectively. The destination map color bar is enumerated from 1 to 33. 1-12 are left lower lobe (LLL), 13-18 are left upper lobe (LUL), 19-26 are right lower lobe (RLL), 27-29 are right middle lobe (RML), and 30-33 are right upper lobe (RUL). The blue color family (from blue to green) corresponds to the left side and the red color family (from yellow to red) corresponds to the right side of the lung. desired time range is identified as the source map (defined at the trachea). This time-thresholding process is carried out for all 380 series of backward injections. Figure 7 exhibits the source maps versus RML destination maps for tref = 0.2, 0.4 and 1.0 second time-points of the respiration cycle. The similarity between the corresponding map pairs confirms the qualitative accuracy of the source inversion method. To further quantify the accuracy of the method, the full MR equation is integrated forward in time for particles released from the source maps of various tref values. Table 1 shows that between 90% to 99% of particles reach the target RML.
Discussion
Jo
4
In this study, a dynamical systems approach was used to study particle transport in lung airways. The FTLE field was calculated and a close connection between the identified ILCS and destination map was demonstrated. Destination map could be created for any cross-section at the trachea to visualize the fate of seeded particles. The ridges of the FTLE field (the ILCS) mark the boundaries of different destinations and provide a mechanistic explanation for the destination map patterns. We also, for the first time, demonstrated the utility of integrating the MR equation backward in time using slow manifolds to create source maps. The source map demonstrates the spatiotemporal origin of dispersed particles. A very good match between the source and destination maps was observed. The connection between FTLE and destination maps in lung airways was previously demonp. 11
re-
pro
of
Journal Pre-proof
urn al P
Figure 5: The FTLE map (left panel) and the initial location of deposited particles (middle panel) are shown for dp = 1 µm. The two fields are superimposed on top of each other and shown in the right panel. Three different injection times are shown in this figure. T = 0.4s is selected when the FTLE patterns become independent of T.
Jo
strated [45, 44]. Our study provides several additional insights. First, we have demonstrated a relationship between the integration time used in calculating FTLE, the resulting ILCS, and particle transport to different generations of bifurcations. Namely, as we increase the integration time, the first emerging ILCS at the trachea marks the boundary between particles separating at the first bifurcation. As the integration time is further increased, new ILCS emerge that correspond to the next generation of airways and bifurcations. We have also demonstrated the temporal evolution of the FTLE and ILCS. That is, by injecting new sets of particles at new time-steps the temporal evolution of these fields could be studied. The temporal evolution of ILCS was shown to control the changes in the destination maps with respect to the time-point of particle release (Fig. 4). This is valuable information since it allows one to track the evolution of ILCS during inspiration to study how sensitive the destination map is to particle release time-point. Interestingly, at the beginning of inspiration when the flow rate is lower, the temporal evolution of ILCS is more pronounced, whereas with an increase in flow rate the ILCS minimally changes in time. These results imply that releasing particles at the beginning of inspiration is more sensitive to the choice of particle release time compared to later stages of inhalation. Finally, we demonstrated minor changes in ILCS results when the particle diameter changed from 1 to 10 µm. Since the Stokes numbers are within a small range (O(10−4 ) to O(10−2 )), the bulk of particles follow the fluid flow, and therefore the FTLE patterns look very similar. We verified that qualitative LCS patterns were identical between the microparticles investigated in this study and tracers (results not shown). This means that the destination maps closely follow the streamline distribution from the trachea to the outlets. Quantitative differences, however, were observed between the FTLE values along the ILCS (inertial p. 12
re-
pro
of
Journal Pre-proof
urn al P
Figure 6: (a) A 3D volume-rendered view of the FTLE field is shown. The red surfaces correspond to the ILCS for dp = 1 µm. (b) The region of interest where 3D FTLE calculation is performed is illustrated. (c) 2D FTLE slices at the trachea (top) and right below the carina (bottom) are shown. (d) Iso volumes of FTLE field with a value of 3.5 are displayed at the same locations. An integration time of T = 0.24s is used in this figure.
Destination map
Source map
tref = 0.2s
tref = 0.4s
tref = 0.6s
tref = 1.0s
Jo
Figure 7: The first row is RML destination map obtained from post-processing of traditional aerosol bolus tracking released at tref = 0.2, 0.4, and 1.0s time-points of inhalation. The second row displays the source map with the corresponding injection times by computing backward MR integral (source inversion). particles) and LCS (ideal tracers). The repelling ILCS mark the boundaries of the flow with high stretching. Particles originating at the different sides of these structures, by definition, separate exponentially fast. In the upper airways, this corresponds to bifurcating vessels that create such separation. The Lagrangian coherent structures are invariant manifolds where particles at the LCS surfaces remain on these structures
p. 13
Journal Pre-proof
of
Table 1: The number of particles at the trachea obtained from backward integration of particles seeded at RML (second column) and the number of particles starting at the source map that reach the RML from solving the full MR equation forward in time (third column) are listed. The success rate is the ratio of the third column to the second column and quantifies the accuracy of the source inversion process. The rows are associated with various injection times. tref = 0s refers to the inhalation starting point.
0.0 0.2 0.4 0.6 1.0 1.2
No. source particles at the trachea from RML (backward integration) 2904 64500 94316 104949 120884 134790
No. particles reaching RML from source map (forward integration) 2659 59595 89701 102351 119114 133408
re-
tref (s)
pro
Source map verification
Success rate percentage 91.6% 93.1% 95.1% 97.5% 98.5% 99.0%
Jo
urn al P
while being advected with the flow [42]. These observations, collectively, imply that particles originating from the ILCS should deposit at the bifurcation regions. Indeed, our results confirmed this conjecture where particles originating at the ILCS ended up depositing at the airway wall (Fig. 5). Therefore, the FTLE and ILCS not only determine the final particle distribution but also accurately identify the origin sites of particles that end up depositing on the bronchial walls. Interestingly, prior studies have shown that release positions of particles at the trachea depositing at the airway wall resemble 2D lines (with structures very similar to our ILCS patterns) [29]. Our numerical results are also consistent with the recent experiments of Amili et al. [2]. These authors showed how, when inertial particles are injected in a successively bifurcating vessel, their fate is strongly dependent on the release location and weakly affected by the flow rate through the distal branches. We also introduced a novel method to create source map with backward in time integration of the MR equation. We presented an example of how the MR equation could be integrated backward in time using the slow manifold technique to overcome the numerically ill-conditioned backward MR equation. The results demonstrate that the source map fully covers the destination map for the target lobe for different reference injection times. In addition to the excellent qualitative agreement, more than 90% of the particles that were integrated backward in time reached the expected location of release. The slow manifold technique is an approximation to the backward MR equation, and therefore perfect accuracy is not expected [47]. Nevertheless, this approach presents an excellent method for identification of the source of any desired particle. For example, instead of launching millions of particles at the trachea to identify particles passing through a specific airway vessel, particles could be launched only at the desired location and integrated backward in time. This could improve the spatial resolution of the destination map by significantly reducing the computational cost. The methods presented in this study have important physiological implications. Inhaled aerosols are commonly used in respiratory disease drug delivery [37]. It is imperative to increase the dosage of the drug reaching the diseased site to increase drug delivery efficacy and reduce harmful side p. 14
Journal Pre-proof
urn al P
re-
pro
of
effects [18, 53]. Construction of destination maps at the trachea to guide spatial and temporal drug injection strategy is the first step to improve efficacy [57, 49]. In this study, we demonstrated the important role that FTLE and ILCS play in these maps. Namely, our novel topological analysis could be used to distinguish between the boundaries of the destination map corresponding to different generations of airway bifurcations. Additionally, the degree of the temporal evolution of ILCS highlights the sensitivity of drug transport results with respect to injection time. The ILCS surfaces can potentially appear, disappear, or move with the flow during inhalation, which highlight the importance of an appropriate temporal release strategy. Finally, source inversion using the slow manifold method could be used to identify drug release location with minimal computational effort. The present study has a few limitations. Particle-particle interactions were ignored in our study. The inlet velocity was assumed to be parabolic although the laryngeal jet can create asymmetric flow profiles. Information from the upstream flow could be considered in future studies. The lower generations of tracheabronchial trees move relatively little, which was not modeled. We have not modeled potential cycle-to-cycle variations in respiratory flow patterns. This work is limited to spherical particles, whereas other types of particles could exhibit different trajectory and results. Our study is limited in its parameter space (just one Re and one Wo number). This limitation is modulated by the fact that this study is aimed to demonstrate the value of the tools, rather than the effect of varying the flow regimes. The results are sensitive to the lung geometry and boundary conditions, thus we expect our results to be different from patient to patient, however, the correlation between FTLE and destination maps is expected to preserve. Finally, the pharmaceutic particles are normally administered by the oral route in clinical practice. Therefore, it would be more practical if destination maps were found at the oronasal cavity. In conclusion, we have presented novel topological tools from dynamical systems theory, which can provide valuable insight into inertial particle transport patterns in lung airways. Our proposed tools assist in coming up with a more accurate spatial and temporal drug release strategy. Future work should extend our framework to the oronasal cavity where more complex and turbulent flow is present. Inertial particle transport in complex, chaotic flows can lead to rich and unique physics [10] (e.g. transient chaos [28]) that has received very little attention in respiratory flows [49].
Conflict of interest None.
Acknowledgement
Jo
The authors would like to thank Kamran Poorbahrami and Jessica Oakes for fruitful discussions related to particle dynamics in lung airways and Milad Habibi for discussions related to slow manifolds.
References
[1] M. R. Allshouse and T. Peacock. Lagrangian based methods for coherent structure detection. Chaos: An Interdisciplinary Journal of Nonlinear Science, 25(9):097617, 2015. [2] O. Amili, J. Golzarian, and F. Coletti. In vitro study of particle transport in successively bifurcating vessels. Annals of Biomedical Engineering, pages 1–13, 2019. p. 15
Journal Pre-proof
[3] A. Arzani, A. S. Les, R. L. Dalman, and S. C. Shadden. Effect of exercise on patient specific abdominal aortic aneurysm flow topology and mixing. International Journal for Numerical Methods in Biomedical Engineering, 30(2):280–295, 2014. [4] A. Arzani and S. C. Shadden. Characterization of the transport topology in patient-specific abdominal aortic aneurysm models. Physics of Fluids, 24(8):081901, 2012.
pro
of
[5] B. Asgharian, O. T. Price, M. Oldham, L. Chen, E. L. Saunders, T. Gordon, V. B. Mikheev, K. R. Minard, and J. G. Teeguarden. Computational modeling of nanoscale and microscale particle deposition, retention and dosimetry in the mouse respiratory tract. Inhalation Toxicology, 26(14):829–842, 2014. [6] A. J. Banko, F. Coletti, C. J. Elkins, and J. K. Eaton. Oscillatory flow in the human airways from the mouth through several bronchial generations. International Journal of Heat and Fluid Flow, 61:45–57, 2016.
re-
[7] W. D Bennett, M. S. Messina, and G. C. Smaldone. Effect of exercise on deposition and subsequent retention of inhaled particles. Journal of Applied Physiology, 59(4):1046–1054, 1985. [8] G. Boffetta, G. Lacorata, G. Redaelli, and A. Vulpiani. Detecting barriers to transport: a review of different techniques. Physica D: Nonlinear Phenomena, 159(1-2):58–70, 2001.
urn al P
[9] A. E BozorgMagham, S. D. Ross, and D. G. Schmale III. Real-time prediction of atmospheric Lagrangian coherent structures based on forecast data: An application and error analysis. Physica D: Nonlinear Phenomena, 258:47–60, 2013. [10] J. H. E. Cartwright, U. Feudel, G. K´arolyi, A. De Moura, O. Piro, and T. T´el. Dynamics of finite-size particles in chaotic fluid flows. In Nonlinear Dynamics and Chaos: Advances and Perspectives, pages 51–87. Springer-Verlag Berlin Heidelberg, 2010. [11] Y. Cheng, Y. Zhou, and B. T. Chen. Particle deposition in a cast of human oral airways. Aerosol Science & Technology, 31(4):286–300, 1999. [12] C. C. Daigle, D. C. Chalupa, F. R. Gibb, P. E. Morrow, G. Oberd¨orster, M. J. Utell, and M. W. Frampton. Ultrafine particle deposition in humans during rest and exercise. Inhalation Toxicology, 15(6):539–552, 2003.
Jo
[13] P. Das, E. Nof, I. Amirav, S. C. Kassinos, and J. Sznitman. Targeting inhaled aerosol delivery to upper airways in children: Insight from computational fluid dynamics (CFD). PloS One, 13(11):e0207711, 2018. [14] A. Farghadan and A. Arzani. The combined effect of wall shear stress topology and magnitude on cardiovascular mass transport. International Journal of Heat and Mass Transfer, 131:252– 260, 2019. [15] G. Haller. A variational theory of hyperbolic Lagrangian coherent structures. Physica D: Nonlinear Phenomena, 240(7):574–598, 2011. [16] G. Haller. Lagrangian coherent structures. Annual Review of Fluid Mechanics, 47:137–162, 2015. p. 16
Journal Pre-proof
[17] G. Haller and T. Sapsis. Where do inertial particles go in fluid flows? Physica D: Nonlinear Phenomena, 237(5):573–583, 2008. [18] D. R. Hess, P. Anderson, R. Dhand, J. L. Rau, G. C. Smaldone, and G. Guyatt. Device selection and outcomes of aerosol therapy: Evidence-based guidelines. Chest, 127(335):335–371, 2005.
of
[19] T. Higenbottam, T. Siddons, and E. Demoncheaux. The direct and indirect action of inhaled agents on the lung and its circulation: lessons for clinical science. Environmental Health Perspectives, 109(suppl 4):559–562, 2001.
pro
[20] W. Hofmann. Modelling inhaled particle deposition in the human lung–a review. Journal of Aerosol Science, 42(10):693–724, 2011. [21] H. Homann and J. Bec. Finite-size effects in the dynamics of neutrally buoyant particles in turbulent flow. Journal of Fluid Mechanics, 651:81–91, 2010.
re-
[22] S. Jalal. Velocity & Vorticity Transport In 3D-Printed Idealized & Realistic Human Airways Using Magnetic Resonance Velocimetry (MRV) & Particle Image Velocimetry (PIV). PhD thesis, University of Minnesota, 2019. [23] C. Kleinstreuer and Z. Zhang. Airflow and particle transport in the human respiratory system. Annual Review of Fluid Mechanics, 42:301–334, 2010.
urn al P
[24] C. Kleinstreuer, Z. Zhang, and J. F. Donohue. Targeted drug-aerosol delivery in the human respiratory system. Annu. Rev. Biomed. Eng., 10:195–220, 2008. [25] A. V. Kolanjiyil and C. Kleinstreuer. Computationally efficient analysis of particle transport and deposition in a human whole-lung-airway model. part I: Theory and model validation. Computers in Biology and Medicine, 79:193–204, 2016. [26] A. V. Kolanjiyil, C. Kleinstreuer, and R. T. Sadikot. Computationally efficient analysis of particle transport and deposition in a human whole-lung-airway model. part II: Dry powder inhaler application. Computers in Biology and Medicine, 84:247–253, 2017. [27] P. G. Koullapis, L. Nicolaou, and S. C. Kassinos. In silico assessment of mouth-throat effects on regional deposition in the upper tracheobronchial airways. Journal of Aerosol Science, 117:164–188, 2018.
Jo
[28] T. Lai, Y. C.and T´el. Transient chaos: complex dynamics on finite time scales, volume 173. Springer New York, 2011. [29] Z. Li, C. Kleinstreuer, and Z. Zhang. Particle deposition in the human tracheobronchial airways due to transient inspiratory flow patterns. Journal of Aerosol Science, 38(6):625–644, 2007. [30] Z. Li, C. Kleinstreuer, and Z. Zhang. Simulation of airflow fields and microparticle deposition in realistic human lung airway models. part I: Airflow patterns. European Journal of MechanicsB/Fluids, 26(5):632–649, 2007. [31] Z. Li, C. Kleinstreuer, and Z. Zhang. Simulation of airflow fields and microparticle deposition in realistic human lung airway models. part II: Particle transport and deposition. European Journal of Mechanics-B/Fluids, 26(5):650–668, 2007.
p. 17
Journal Pre-proof
[32] P. W. Longest, G. Tian, R. L. Walenga, and M. Hindle. Comparing MDI and DPI aerosol deposition using in vitro experiments and a new stochastic individual path (SIP) model of the conducting airways. Pharmaceutical Research, 29(6):1670–1688, 2012. [33] M. R. Maxey and J. J. Riley. Equation of motion for a small rigid sphere in a nonuniform flow. The Physics of Fluids, 26(4):883–889, 1983.
of
[34] M. Mortensen and K. Valen-Sendstad. Oasis: A high-level/high-performance open source navier–stokes solver. Computer Physics Communications, 188:177–188, 2015.
pro
[35] A. Naseri, O. Abouali, P. Ghalati Farhadi, and G. Ahmadi. Numerical investigation of regional particle deposition in the upper airway of a standing male mannequin in calm air surroundings. Computers in Biology and Medicine, 52:73–81, 2014. [36] J. M. Oakes, A. L. Marsden, C. Grandmont, S. C. Shadden, C. Darquenne, and I. E. VignonClementel. Airflow and particle deposition simulations in health and emphysema: from in vivo to in silico animal experiments. Annals of Biomedical Engineering, 42(4):899–914, 2014.
re-
[37] U. Pison, T. Welte, M. Giersig, and D. A. Groneberg. Nanomedicine for respiratory diseases. European Journal of Pharmacology, 533(1-3):341–350, 2006. [38] K. Poorbahrami and J. M. Oakes. Regional flow and deposition variability in adult female lungs: A numerical simulation pilot study. Clinical Biomechanics, 66:40–49, 2019.
urn al P
[39] S. Qi, B. Zhang, Y. Teng, J. Li, Y. Yue, Y. Kang, and W. Qian. Transient dynamics simulation of airflow in a ct-scanned human airway tree: more or fewer terminal bronchi? Computational and Mathematical Methods in Medicine, 2017:1–14, 2017. [40] T. Sapsis and G. Haller. Inertial particle dynamics in a hurricane. Journal of the Atmospheric Sciences, 66(8):2481–2492, 2009. [41] S. C. Shadden and A. Arzani. Lagrangian postprocessing of computational hemodynamics. Annals of Biomedical Engineering, 43(1):41–58, 2015. [42] S. C. Shadden, F. Lekien, and J. E. Marsden. Definition and properties of Lagrangian coherent structures from finite-time Lyapunov exponents in two-dimensional aperiodic flows. Physica D: Nonlinear Phenomena, 212(3-4):271–304, 2005.
Jo
[43] S. C. Shadden and C. A. Taylor. Characterization of coherent structures in the cardiovascular system. Annals of Biomedical Engineering, 36(7):1152–1162, 2008. [44] B. Soni and S. Aliabadi. Large-scale CFD simulations of airflow and particle deposition in lung airway. Computers & Fluids, 88:804–812, 2013. [45] B. Soni, D. Thompson, and R. Machiraju. Visualizing particle/flow structure interactions in the small bronchial tubes. IEEE Transactions on Visualization and Computer Graphics, 14(6):1412–1427, 2008. [46] M. Sudharsan, S. L. Brunton, and J. J. Riley. Lagrangian coherent structures and inertial particle dynamics. Physical Review E, 93(3):033108, 2016. [47] W. Tang, G. Haller, J. Baik, and Y. Ryu. Locating an atmospheric contamination source using slow manifolds. Physics of Fluids, 21(4):043302, 2009. p. 18
Journal Pre-proof
[48] K. Tor´en, I. A. Bergdahl, T. Nilsson, and B. J¨arvholm. Occupational exposure to particulate air pollution and mortality due to ischaemic heart disease and cerebrovascular disease. Occupational and Environmental Medicine, 64(8):515–519, 2007. [49] A. Tsuda, F. S. Henry, and J. P. Butler. Particle transport and deposition: basic physics of particle kinetics. Comprehensive Physiology, 3(4):1437–1471, 2011.
of
[50] A. Updegrove, N. M. Wilson, J. Merkow, H. Lan, A. L. Marsden, and S. C. Shadden. Simvascular: An open source pipeline for cardiovascular simulation. Annals of Biomedical Engineering, 45(3):525–541, 2017.
pro
[51] T. Van de Moortele, U. Goerke, Chris H. Wendt, and F. Coletti. Airway morphology and inspiratory flow features in the early stages of chronic obstructive pulmonary disease. Clinical Biomechanics, 66:60–65, 2017.
re-
[52] T. Van de Moortele, C. H. Wendt, and F. Coletti. Morphological and functional properties of the conducting human airways investigated by in vivo computed tomography and in vitro mri. Journal of Applied Physiology, 124(2):400–413, 2017. [53] J. Vestbo, S. S. Hurd, A. G. Agust´ı, P. W. Jones, et al. Global strategy for the diagnosis, management, and prevention of chronic obstructive pulmonary disease: GOLD executive summary. American Journal of Respiratory and Critical Care Medicine, 187(4):347–365, 2013.
urn al P
[54] J. Xi, X. A. Si, J. W. Kim, E. Mckee, and E. B. Lin. Exhaled aerosol pattern discloses lung structural abnormality: a sensitivity study using computational modeling and fractal analysis. PloS One, 9(8):e104682, 2014. [55] X. Xu, Y. Shang, L. Tian, W. Weng, and J. Tu. Inhalation health risk assessment for the human tracheobronchial tree under pm exposure in a bus stop scene. Aerosol and Air Quality Research, 19(6):1365–1376, 2019. [56] B. Zhang, S. Qi, Y. Yue, J. Shen, C. Li, W. Qian, and J. Wu. Particle disposition in the realistic airway tree models of subjects with tracheal bronchus and COPD. BioMed Research International, 2018:1–15, 2018. [57] Z. Zhang and C. Kleinstreuer. Effect of particle inlet distributions on deposition in a triple bifurcation lung airway model. Journal of Aerosol Medicine, 14(1):13–29, 2001.
Jo
[58] Z. Zhang, C. Kleinstreuer, J. F. Donohue, and C. S. Kim. Comparison of micro-and nano-size particle depositions in a human upper airway model. Journal of Aerosol Science, 36(2):211–233, 2005. [59] Z. Zhang, C. Kleinstreuer, and C. S. Kim. Aerosol deposition efficiencies and upstream release positions for different inhalation modes in an upper bronchial airway model. Aerosol Science & Technology, 36(7):828–844, 2002.
p. 19