Abstract
For many of the one billion sufferers of respiratory diseases worldwide, managing their disease with inhalers improves their ability to breathe. Poor disease management and rising pollution can trigger exacerbations which require urgent relief. Higher drug deposition in the throat instead of the lungs limit the impact on patient’s symptoms. To optimise delivery to the lung, patient-specific computational studies of aerosol inhalation can be used. However in many studies, inhalation modelling does not represent an exacerbation, where the patient’s breath is much faster and shorter. Here we compare differences in deposition of inhaler particles in the airways of a healthy male, female lung cancer and child cystic fibrosis patient. We aimed to evaluate deposition differences during an exacerbation with image-based healthy and diseased patient models. We found that during an exacerbation, particles progressing to the lower airways were distributed similarly to those inhaled during healthy breathing, but fewer in quantity. Throat deposits were halved in the healthy patient compared to the diseased patients under extreme inhalation, due to changes in the detailed shape of the throat. Our results identify that the modelled upper airway must be patient-specific, and an exacerbating profile tested for optimal measurement of reliever inhaler deposition.
1. Introduction
More than one billion people worldwide suffer from asthma, cystic fibrosis, and other chronic respiratory diseases (The Global Asthma Network, 2018), many experiencing distress and anxiety due to restrictions to activities and limited productivity (Dockrell et al., 2007). These limitations are most prominent among the young and elderly (The Global Asthma Network, 2018). One of the largest contributors to the diseased population is asthma, which incurs an annual cost per patient of €1,700 and $3,100 in Europe and the USA, respectively (Nunes et al., 2017) from direct cost of treatment and indirect costs such as work absence or decreased productivity (Katsaounou et al., 2018; Gruffydd-Jones et al., 2019). Similar impacts are induced from cystic fibrosis, a highly common hereditary disease in the UK, USA and Australia (Elborn, 2016). Consistent treatment in alignment with disease management plans are recommended to minimise symptoms (Ring et al., 2015), but adherence is an issue in young patients (McQuaid et al., 2003) and many adults are purposely inconsistent to limit exposure to side-effects such as osteoporosis and cataracts (Dockrell et al., 2007). Even in adherent patients, efficiency of the metered-dose inhaler varies greatly across patients (Clark, 1995) as many (particularly children) experience difficulties in device technique (Usmani, 2019) due to the rapid spray of the drug. The issue of technique (patient breathing pattern and coordination with device actuation) and differences in lung structure are the main influences in drug delivery (Darquenne et al., 2016).
Optimisation of the medication deposition could be achieved through in silico analyses, by providing the clinician information on the local deposition and therapeutic outcome. One available deposition tool is the Multiple-Path Particle Dosimetry (MPPD) model (Anjilvel and Asgharian, 1995; Asgharian et al., 2001). However this calculates deposition of ambient particles (Borghardt et al., 2015), which does not mirror the physics of spray aerosol inhalation (Longest et al., 2008). This issue has been recognised and a commercial counterpart to predict deposition of pharmaceutical aerosols has been developed (Olsson and Bäckman, 2018). However, as the equations used are based on probabilistic 1D equations, similar to that of MPPD, complex fluid phenomena and particle interactions cannot be included. Comparisons between 1D models and computational particle-fluid dynamics (CPFD) deposition led Zhang et al. (2009) to observe significant differences in local deposition, attributed tolocal flow. Flow and particle phenomena can be readily accounted for in CPFD by solving equations governing the transport of air and particles (e.g. see (Sundaresan et al., 2018)). Flow can be solved in airways extracted from medical images and particle properties can represent inhaler particles to produce patient-specific deposition analyses. But with the complexity of the flow regimes in the system creating demanding simulations and the massive number of respiratory patients, the benefit to cost of image-based modelling is an issue. Work exists interpreting pharmaceutical deposition differences in adult patients (van Holsbeke et al., 2018; Poorbahrami and Oakes, 2019; Poorbahrami et al., 2019), but only sparse research models deposition in children. Existing children-focused image-based deposition studies analyse the nasal cavity (Xi et al., 2011, 2012) or central airways (Das et al., 2018; Oakes et al., 2018). To move towards enhanced treatments for all ages, a study considering the impact of the upper airways is needed. We begin this by comparison of deposition in the airways (from mouth to central airways) in three diverse patients ranging from 11 - 49 years old. This can allow simulation of flow created in the upper airways and provide understanding of how this affects deposition throughout the airways, across patients.
Results from these simulations are dependent on the inflow conditions, which here is based on the duration and strength of the patient’s inhalation. Colasanti et al. (2004) published breathing profiles of patients who had recently had an exacerbation (around one hour before measurements were taken) show vastly different shape and magnitude from the inhalation waveform which is typically used in respiratory CPFD studies (Inthavong et al., 2010). A sinusoidal inhalation waveform is used by most existing studies (Oakes et al., 2018; Inthavong et al., 2010; Naseri et al., 2017), which mimics a healthy patient’s tidal breathing. The variation in flow patterns produced by these differences would therefore give a deposition analysis which may not represent the reliever inhaler’s true operating conditions. Different breathing profiles have been utilised (Longest et al., 2012; Khajeh-Hosseini-Dalasm and Longest, 2015), but these were used to represent techniques for different devices (dry-powder compared to metered-dose inhalers), not exacerbating and healthy breathing in a metered-dose inhaler. By applying different inflow conditions, one can understand changes in flow structure during an exacerbation or another desired breathing situation. This knowledge would allow manufacturers to tailor inhaled therapeutics to ensure optimal dosage reaches the desired site in the airways, and to provide clinicians and patients the tools with which to maximise inhalation technique during exacerbations.
In addition, existing MDI deposition research simplifies particle interactions, which includes treating particle-wall collision as complete sticking. This is a coarse approximation of dissipative lubrication forces between the particle and wall (Legendre et al., 2005; Holbrook and Longest, 2013). To treat as fully sticking therefore neglects rebound and could overestimate deposition of high inertia particles. Furthermore, MDI studies often exclude van der Waals forces which can cause small particles agglomerate (Hamaker, 1937), and become more inertial. If the increase in inertia is large enough, this can change particle trajectory or chance of rebound, which, in turn would cause early deposition. How these forces alter deposition and compete against each other should be tested. Applying this to simulations made specific to a patient’s airway structure under extreme and optimal breathing will further understanding of drug transport across different health states.
Therefore, we aimed to (i) evaluate the need for patient-specific domains in future simulations through medical image-based modelling of three diverse patients. To satisfy this we also aimed to (ii) identify the necessary physical effects to produce accurate models of the system (including the mass of drug simulated, van der Waals and particle-wall lubrication forces), and (iii) evaluate variation in deposition during exacerbating breathing.
2. Methods
To answer the research aims above, we evaluated changes in deposition produced by the following variations. Parameters varied included the mass of drug simulated, particle cohesiveness, lubrication forces in particle wall-collisions, and deletion or saving of deposited particles from the system to mimic absorbing particles with mucus layer. The use of a healthy and exacerbation breathing profile (Colasanti et al., 2004) allowed for analysis of the reliever inhaler during an exacerbation. The effects of patient variation were then analysed using this optimised modelling setup, through comparison under realistic inflow conditions in three patients.
2.1. Medical image processing
Three patients were studied retrospectively using computed tomography (CT) (detailed in Table 1). Use of these retrospective images was approved by Heriot-Watt University (ID: 2020-0500-1452). The surface file of the healthy patient’s segmented airways was provided by Dr Filippo Coletti (University of Minnesota). The patients had sufficient variation in age, gender and health status to gather an indication of the benefit of patient-specificity.
Images were processed using 3D Slicer 4.10 (Fedorov et al., 2012). Images were preprocessed by applying isotropic spacing to account for the anisotropic resolution of CT scans (Table 1). For the cancer patient, voxel size was then reduced using a linear interpolation, creating a voxel size of ~ 0.4 mm to allow extraction of smaller airways. This was also performed on the cystic fibrosis patient. To preserve edges of the airway while blurring lung tissue, an anisotropic diffusion filter was then applied (Duan et al., 2019). Parameters used were conductance 3, 5 iterations and step size of 0.06025. These were based on visual comparison of our segmentations produced from different values of conductance used to good success in medical imaging (Behnaz et al., 2010; Sen et al., 2011).
Images were then segmented using a threshold-based region-growing approach to segment areas classified as air without leakage to the background air (Nardelli et al., 2015; Mayer et al., 2004; De Nunzio et al., 2011; Aykac et al., 2003). This semi-automatic process grows the segment from ‘seeds’ which are user declared points within each region, set separately for the airway and surrounding lung tissue. Once each region was grown, labels denoting leaked regions were found by overlaying the segmentation on the scan. Leaked regions were then labelled as airway and the region was grown again until the scan was observed to be acceptable quality. This process was applied to the external airways, right and left lungs separately due to variations in the airway image density in each region (Nardelli et al., 2015). This allowed segmentation to a depth ranging from the fourth to sixth bifurcation level (G4 – G6, Figure 1).
Some of the upper airways were not included in the CT data, we have found this to be common in most clinical CT scans. To account for this we merged the oral cavity of the healthy patient to the throat of the lung cancer patient. This region was then extracted, scaled and joined to the trachea of the cystic fibrosis patient to complete the missing regions. Scaling was performed such that the intersection of the new and existing regions matched in diameter (resulting in a scaling of 0.8 for the added region). This was later repeated using the throat of healthy patient, after being advised such a narrow, obstructed throat is not characteristic of cystic fibrosis, and likely unique to the cancer patient. The regions which were artificially added are shown graphically in Figure 1 within the dashed box. These excluded regions are of large importance in inhaler simulations as a large portion of the dosage is lost within this part of the airway and the turbulence created here is cascaded through the trachea and main bronchus (Banko et al., 2015).
2.2. Mathematical modelling
Here we present the mathematical relations used to determine the physical factors included in our model of the system. We provide the equations governing the fluid and particle solvers in Appendix A. Briefly, particle transport was solved by the discrete element method (DEM) using the particle simulator LIGGGHTS (Kloss and Goniva, 2011). This tracks individual particles’ trajectories by integrating Newton’s equations of motion in time for each particle (Verlet, 1967). Particle collisions were modelled as a linear spring-dashpot system (Cundall and Strack, 1979). Fluid transport through the airways was solved by the volume-filtered mass and momentum continuity equations (Anderson and Jackson, 1967; Capecelatro and Desjardins, 2013) implemented in OpenFOAM v2.2 (Weller et al., 1998). Particle and fluid phases were coupled through a version of the CFDEMcoupling platform (Kloss et al., 2012) modified to benefit from faster two-way coupling by Ozel et al. (2016).
Filtering the fluid transport equations creates unresolved stresses. In our simulations, we consider stress contributions from the gas pressure gradient () and residual stresses from volume filtering of fluid velocity fluctuations (Ru) below the cell size, Δ. Ru is dependent upon the eddy viscosity (μt), a term representing turbulence dissipation into the smaller, unresolved scales. Eddy viscosity is modelled using a dynamic Smagorinsky model (Germano et al., 1991; Lilly, 1992). This models the energy transfer to flow structures smaller than the cell size without a fixed model constant, as this is computed dynamically.
We solve particle drag from the relation derived by Beetstra et al. (2007) for monodisperse particles. This drag law is based on the particle’s Reynolds number, Rep, given as
Where Vr,i is the relative particle velocity, dp is the particle diameter, μf is the viscosity. As well as determining particle drag, Rep also characterises particle inertia through the particle’s Stokes number,
Where D is the airway diameter. Due to their small size, metered-dose inhaler particles have a Stokes number below one, meaning they show high sensitivity to carrier fluid fluctuations induced by the complex structure of the upper airways (Kleinstreuer and Zhang, 2010). This makes the Stokes number the primary parameter to be adjusted when optimising deposition (Kleinstreuer and Zhang, 2010).
Again, due to the particle’s small size, particle-particle cohesion from van der Waals forces may influence deposition. Particle attractive energy due to van der Waals force is determined by the material’s Hamaker constant, A (Hamaker, 1937). To permit a larger timestep, particles’ elastic properties are softened. Therefore the real stiffness (kR) is reduced to a soft stiffness (kS), and the Hamaker constant is amended by the relationship AS = AR (kS/kR)1/2 using the model of Gu et al. (2016a). This reduction of real particle stiffness (kR) to a softer stiffness (kS) has negligible effect on fluid hydrodynamics and particle cohesion (Gu et al., 2016a; Ozel et al., 2017). To determine A, which is not given in literature for the metered-dose inhaler propellant HFA-134A, we evaluate deposition at three Bond numbers, Bo. Bo is provided by Ozel et al. (2017) as
Where is the maximum van der Waals force magnitude occuring when (the minimum separation distance for a particle at its real stiffness). This was varied three orders of magnitude, Bo = 10,100,1000 which was sufficient variation to interpret cohesive differences. This magnitude of variation was chosen as it showed changes in deposition without running a large amount of simulations at finer Bo intervals. This also gave A at Bo = 1000 of the same order of magnitude (A = 10−20 J) to drug particles in the inhaler propellant HFA-227 (Engstrom et al., 2009).
Further complexities arise when considering particle-wall interactions. This is widely treated as a fully plastic collision due to the presence of a respiratory mucus layer (Miyawaki et al., 2012; Zhang et al., 2018; Chen et al., 2012; Naseri et al., 2017). This may occur due to lubrication interactions, which have been experimentally shown to damp collision forces (Legendre et al., 2005) due to the formation of a thin interfacial film during contact. This relation has been shown to follow the expression for solid particles (Legendre et al., 2005). Where subscript R and C are the rebound and precollision velocities, respectively, νT is the terminal velocity of the particle, before it slows due to interaction with the wall and Stcoll = (ρp + ρf)dpνT/9μf (Legendre et al., 2006). e follows a sigmoid trend when plotted against Stcoll, with Stcoll < 10 creating a plastic collision similar to that approximated in modelling of particle-mucus layer interactions. We implemented this relationship to model particle-wall contact force, then when later deleting or freezing deposited particles, use it to determine the cutoff. We use this to model particle-mucus layer collision instead of the typical ‘sticking’ condition (Miyawaki et al., 2012; Zhang et al., 2018; Chen et al., 2012; Naseri et al., 2017). This allows potential for particle rebound after impacting the wall. Therefore, as the particle-wall interaction is better represented, deposition is not over-estimated.
2.3. Simulation configuration
2.3.1. Single-phase validation
We first verified the fluid phase by qualitative comparison to a published experimental study (Banko et al., 2015). The study observed water flow in a 3D printed hollow cast of an adult male patient’s airways (Figure 1a) at a Reynolds number representative of heavy breathing ( corresponding to Q = 1L/s). Water flow was simulated to first replicate the experiment before advancing to a simulation of airflow to show independence of the fluid simulated when velocity is scaled by its Re. Water was simulated using OpenFOAM solver pimpleFoam, then when simulating air we used the CFDEM solver with no particles. We gauged mesh sensitivity by simulating two uniform, hexahedral meshes (5 × 105 and 106 cells). As interaction between the gas and particle phases (namely, drag) is dependent on the size of Δ relative to dp, we made our fine mesh cell size Δ = 50dp to minimise drag overestimation. This is of the same order of magnitude as coarse-fluid grid simulations performed by Radl and Sundaresan (2014). We kept grid size consistent with our first particle size tested (at dp = 10 μm, Δ = 500 μm) throughout the study to minimise excessive computation time when later reducing particle size.
As upstream fluctuations in the experimental apparatus necessitate a turbulent inlet (Sagaut, 2006), we implemented a boundary condition which adds random components to the flow. This was used with an inflow velocity of uwater = 0.167 m/s and uair = 2.677 m/s, with a no-slip condition at the wall. Outlets had a uniform, fixed pressure applied to each bronchi. Although not truly representative of outflow conditions within the lung, a uniform outlet condition has been shown to be acceptable when the bronchi are extended a few diameters in the axial direction to eliminate secondary flow effects (Zhang et al., 2012). This has been shown to resemble more advanced outlet conditions which consider compliance of the lung (Ma and Lutchen, 2006).
2.3.2. Aerosol transport modeling
In multiphase simulations we coupled the single-phase approach above to the DEM solver to track monodisperse particles of dp = 10 and 4 μm. We began with 10 μm, which are on the upper end of those used in metered-dose inhaler research (Kleinstreuer and Zhang, 2010; Watanabe and Watanabe, 2019), to allow faster simulations at this stage to observe variation between patients. This is as timestep size and number of particles (for a desired mass of drug) are dependent upon particle diameter. Particles are classed as deposited and deleted from the system when impacting the wall with low inertia (Stcoll < 20 Legendre et al. (2005)). As CT scan resolution only permitted segmentation to approximately the sixth bifurcation level (Figure 1), particles reaching the end of the bronchial path were classified as reaching the distal airways and deleted. As we this only model a limited amount of bifurcation levels (between four to six), this is a coarse approximation of downstream behaviour and dosage reaching the small airways.
The dosage was released over a period of t = 0.1 s (Ju et al., 2010). The particle motion was discretised in time based on calculation of the collision time of the particles, (Gu et al., 2016a). The van der Waals stiffness scaling (Gu et al., 2016a) was set to AS = AR (kS/kR)1/2 = AR/31.6. Particle lubrication force was modelled using mucus layer viscosity ranging from μmucus = 0.026 − 0.05 Pa · s dependent on disease (Rubin, 2007) and a density of ρmucus = 1000 kg/m3 as it is largely made up of water in the upper airways (Olsson et al., 2011). Parameters used are summarised in Table 2. Simulations were performed on high-performance computers ROCKS of Heriot-Watt and EDDIE of the University of Edinburgh high-performance computers using 28 − 64 CPU cores.
Once the optimal parameters were determined, a time-varying inlet condition was implemented to represent a real breathing cycle (Figure 2). For this we used data provided by Colasanti et al. (2004) whom analysed the breathing profile of patients suffering from COPD and cystic fibrosis shortly after an exacerbation. This was also compared to a healthy inhalation profile. To compare the effect of inhalation waveform, the healthy patient was simulated using both the healthy and exacerbation profile presented by Colasanti et al. (2004). The two diseased patients were then simulated under the exacerbation inhalation profile. To evaluate the effect of the missing throat and mouth of the cystic fibrosis patient, we model this patient using the mouth and throat of the healthy and cancer patient. The magnitude of the child’s inhalation velocity was lowered by half, based on inhalation velocities used in a similar study (Longest et al., 2006). Following this we focus on the two patients of main interest, the child cystic fibrosis patient and the healthy adult. We compare the effect of reducing particle size from 10 μm to 4μm, and alternating the healthy and exacerbation breathing profile.
2.4. Deposition assessment
In all of these cases particle deposition was analysed by regional groupings of the mouth, throat, trachea, main bronchus and bronchi within each lobe of the lung (as used by Asgharian et al. (2001) and van Holsbeke et al. (2018)). We also provided results grouped into external airways (mouth to end of first branch), and internal airways (airways within the lungs). We use this to evaluate preferential deposition within airways’ different regions.
In contrast to the fully plastic particle-wall condition used in literature, our particles may rebound as we resolve collisions over multiple timesteps. We classed particles that had low inertia prior to impacting the wall as deposited. This inertia was based on Equation (4). We compared results where these particles were kept active in the domain to interact with floating particles, or where we simply deleted them. The difference in these was found to be minimal (under 5% in all regions). Therefore we opted to delete them due to superior computational efficiency (reached 0.1 s physical time in 30% faster clock time).
Due to the particle’s ability to slide along the wall in our parameter study simulations, when stuck it would not be completely stationary. To extract particles on the wall for comparison to our sticking condition, deposition was defined when the particle velocity was sufficiently below that of the free-stream (ν = 0.01 m/s, 900 times less than gas velocity in the throat). This velocity cutoff was found through observed comparison of velocities of slowly floating and deposited sliding particles. Particles that were below this threshold were classed as deposited deleted during post-processing.
Due to the locally-acting nature of metered-dose inhalers (Lu et al., 2015), it is important to understand therapeutic distribution and dosage experienced by the patient based on deposition concentration (Solomon et al., 2012). We interpreted this through the dosimetry measure of deposition enhancement factor (DEF) (Balashazy et al., 1999; Longest et al., 2006). We calculate this using the number of particles deposited within a fixed distance (1 mm, area Aconc = π (1 mm)2) of the central point of each wall face. This distance was based on that used by Dong et al. (2019). Other studies have used much narrower radii (Longest et al., 2006; Xi et al., 2012), but using this wider radius can account for particle translocation during the time between deposition and absorption. This is made relative to the global deposition by
We use this instead of the sum of the deposition efficiencies as the defined areas may overlap, and this therefore prevents particles being counted multiple times towards the global average.
3. Results
3.1. Single-phase validation
The sensitivity of aerosol particles to changing flow structures (as St < 1) makes singlephase flowfield validation crucial before considering particles. Turbulence induced in the upper airways is the foundation of the flow structure, creating secondary flows which are responsible for deposition in the trachea (Jin et al., 2007; Kleinstreuer and Zhang, 2010). Therefore, the upper airways are a suitable region for the validation (Figure 3). All results have been normalised by the mean velocity in the trachea (VT,water = 0.22 m/s and VT,air = 3.51 m/s). The strong jet of flow formed at the throat matches the experimental study well in magnitude and structure, capturing the recirculation zones as the throat expands well. Banko et al. (2015) gave the bulk (area-averaged) velocity relative to VT as 1.73 in the glottis, compared to our value of 1.77, producing a 2.3% relative error.
The resemblance in flow structure continues in the lower airways, from both coronal and axial views as given in Figure 4 (see Banko et al. (2015) for further comparison). The slight separation seen at the left bronchus (right side of (a)) agrees well with the experimental data. Contours in the axial region (shown in (b)) agree well in shape, accurately capturing recirculation and asymmetrical flow. However, velocity vectors shown in C-C’ and D-D’ differ from the published results. This likely stems from differences in time-averaging of such a sensitive parameter in highly unsteady flow. Vectors in B-B’ and E-E’ match the experimental data well. Qualitatively comparing numerical and experimental velocity fields (Banko et al., 2015) verifies that the single-phase flow configuration is suitable to capture the flow patterns present. To quantitatively validate the solver, we compare the relative bulk velocity (/VT). Banko et al. (2015) gave this as 1.03 in the trachea, matching our solver exactly. In the left main bronchus, we had a value of 0.67, differing from the experimental results of 0.70 by 4.3%. In the right bronchus the numerical and experimental and 0.87, respectively (relative error 1.2%). These relative bulk velocity comparisons at key cross sections of the airway validate that our solver can reproduce the in vitro respiratory velocimetry measurements of Banko et al. (2015) to an error below 5%.
3.2. Effect of parameter variations
To model inhaler inhalation the simulation’s sensitivity to the dosage simulated, van der Waals forces and particle-wall lubrication forces were evaluated. The influence of each of these were examined by their influence on dosage deposition.
As physical parameters were varied, one can observe that deposition distinctions are mostly minor, as deposition changes across each lobe were all under 1.5% of the total dosage (Figure 5). When increasing the number of particles the only change came at Np = 1,000,000 (Figure 5c), as the throat deposition rose by 2%, in all other regions the difference was less than 1%. For this reason 200,000 particles was chosen for the remaining simulations to reduce run times. Although we have purposely underestimated Np, the number of particles should not be chosen arbitrarily. Instead the sensitivity of the model to this parameter should be included as an early simulation in all studies where the true payload cannot be simulated. When including particle cohesion and particle-wall lubrication (Figure 5h,i), deposition in the throat due to particle-wall lubrication forces rose by 13% of the total dosage. This effect is reduced by 4% when modelling van der Waals forces, likely explained by the agglomeration of particles providing additional inertia (increasing Stcoll), thus reducing the energy lost due to lubrication. This also explains why minimal changes are observed in Figure 5(d – f).
The parameters taken forward to evaluate deposition variance in the diseased patients were: a dosage of 0.126 μg (Np = 200,000 at dp = 10 μm), a Bond number of 1000, and particle-wall lubrication forces. Time-varying breathing profiles were also applied (Figure 2). Although particle count is underestimated, the error is within 5%, for a computation time reduction of 62% (Np = 200,000 took three weeks, whereas Np = 1,000,000 took eight). These parameters were considered to model the particle behaviour accurately, with a small sacrifice made to reduce computational cost.
3.3. Inter-patient variation
Due to the dominance of upper airway deposition, it is difficult to visualise differences in central airway deposition. When viewing deposition as a logarithm the behaviour in the lobes can be analysed with greater ease (Figure 6). Across the five simulations the dosage was mainly deposited in the external airways (Figure 6), deposition here was always greater than half of the payload. Given the size of the particles compared to the distribution used in inhalers, a large deposit here was expected. For the most constricted external airways (Figure 6c,d) the deposition fraction was 0.95 and 0.85, respectively. The lower throat deposit in the child, despite using a scaled replica of the cancer patient’s airways above the trachea, indicating the weaker inhalation is beneficial here. The cystic fibrosis patients’ deposition dropped by 13.4% in the trachea (relative to total dose) when using the dilated airway, however deposition rises elsewhere meant the total drop was 12.5% (Figure 6d, e). There is anexternal airway deposition rise of 1.9,1.7 and 1.5 times when comparing the healthy patient with exacerbation breathing against the cancer and cystic fibrosis patients (constricted and dilated), respectively. Considering only ‘stuck’ particles (non ‘exiting’ particles), the ratio of external to internal airway deposition ranged from 3 to 61.
Both patients with the constricted throat (Figure 1ii, iii) received below 5% of the dosage to the deep airways (Figure 7). When replacing the constricted throat with the dilated throat in the cystic fibrosis patient, dosage in this region increased to 14% of the total payload. In the two healthy patient simulations this number increased to 31 and 19% for healthy and the exacerbation inhalation. Lobar distribution of these particles for healthy and the exacerbation inhalation was mainly directed to the lower lobes (67 and 56% of total exited) and the left:right lung ratio of dosage reaching the distal airways was 1.1 and 1.6. Left:right deep lung ratio was 0.2 for the lung cancer patient, 0.3 for cystic fibrosis patient (i) and 0.91 for cystic fibrosis patient (ii). This shows a noticeable difference between healthy and diseased airways.
3.4. Particle size reduction
As improvements to device efficiency are generally realised through the patient’s breathing and the particle sizes, we applied these to the two main cases of interest (the healthy adult and cystic fibrosis child patients). Particle size was reduced to 4 μm using the exacerbation profile. We applied the healthy breathing profile to both patients with 10 μm particles. Deposition in the external airways gives no clear trend of improvement (Figure 8). We see the change in breathing yields a deposition reduction of 43% for the cystic fibrosis patient, but a rise of 4% in the healthy patient. Similarly reducing the particle size yields a deposition increase of 3% for the cystic fibrosis patient, but a deposition reduction of 56%. This reduction does not appear to impact the dosage reaching the deep lung, as during exacerbating breathing this value does not differ by more than 0.5% for 10 and 4 μm particles.
3.5. Dosimetry assessment
In all simulations, hotspots were observed at the back of the mouth, throat, upper trachea, first and second bifurcations (Figure 9). The healthy patient’s drug deposits were limited to this region under normal breathing (Figure 9a,e), but during an exacerbation the particles form a uniform coating across the airways with a DEF value of 1–5 in the bronchi, 20 at the first bifurcation, and 20–250 in regions of the throat. This is much larger than under healthy breathing, where the in the throat the range is 10–100, and 10 at the first bifurcation. We see less intense hotspots in bronchi of the two diseased patients. At the first bifurcation, DEF = 0.2 and 5 for the cancer and cystic fibrosis patients (i), respectively. The dose is concentrated in the throat at the narrow throat, where in both patients the values range from 20 to 120. Trachea deposition is favoured to the front as the mean DEF here is 1, and 1.6 for the cancer and cystic fibrosis patient, respectively (Figure 9c,d). This is much smaller (the order of 10-3) at the rear of the trachea in these patients (Figure 9g,h). The cystic fibrosis patient with the constricted and dilated throat shows similar local behaviour below the carina, as there are sparse patches of drug deposits. These are narrow and more intense in the cystic fibrosis patient with the dilated throat (DEF = 1 to 13 in the first few bifurcations, decreasing to 0.1 to 1.5 in the distal bronchi), but appear the slightly larger in area in the cystic fibrosis patient with the constricted throat. Mean DEF at the first bifurcation is 5.5 (Figure 9e), only 0.5 higher than the same patient with the constricted upper trachea (Figure 9d).
4. Discussion
In this study we aimed to evaluate deposition in patient-specific models during exacerbating breathing, compared to during healthy breathing. This allowed the investigation to consider short-acting bronchodilators (‘reliever’ inhalers). By investigating this in a small, diverse patient group we found patient-specific domains to be a necessity in future studies. Results showed less drug reaching the deep lung during an exacerbation, but the main differences came in the patients with distinct changes in upper airways (shape and size). The increased upper airway deposition is attributed to the greater constriction and complexity here, particularly in the cancer patient (Figure 1b). From the brief parameter study, this appeared to overrule the effect of particle-particle interactions such as collisions. This shows that patient-specificity in the upper airways is crucial for accurate deposition prediction.
We have seen unexpected reactions from some patients to reducing particle size and applying the healthy breathing profile (Figure 8). Typically smaller particles and steadier breathing reduce mouth and throat deposition as shown in vivo by Usmani et al. (2005). This further emphasises the need for the use of patient-specific domains, as we have seen a small rise in external airway deposition when reducing particle size from 10 to 4μm for the cystic fibrosis patient. Conversely, for the healthy patient we see a external airway deposition rise of over double with this change in particle sizes. This agrees with reduced deposition in central airways for smaller particles (Usmani et al., 2005). Similar trends were seen when comparing breathing profiles in the healthy patient, as there was a small rise in external airway deposition with the healthier breathing profile (but also a rise in deep lung dosage). A potential cause for these abnormal observations could be the unpredictability of the mesoscale interactions such as particle collisions, van der Waals forces and particle-wall lubrication, and their change with the changing flow structures in different patient airways (Figure 5). In contrast to other patient-specific deposition studies that did not consider particle-particle interactions (Choi et al., 2018; Koullapis et al., 2018; Poorbahrami and Oakes, 2019), we chose to use CFD-DEM with particle-particle interactions. This was as cohesive van der Waals forces have been known to influence the transport of micro-scale particles (Gu et al., 2016a). Soft-sphere (treating particles as deformable) deposition studies have been implemented before but only to observe dynamics in simplified airways (Chen et al., 2012; Wang et al., 2017). This allows understanding of particle dynamics in a fixed system. However, our results in Figure 6 have shown, there are changes in flow dynamics across different patients which may alter the influence of particle interactions. For example, in narrow airways the particle concentration will be denser than in a patient with wider airways and cause more cohesion. This is a speculative explanation as our results in Figure 5 show small variation with cohesive forces. The minor influence of cohesion in our simulations is aligned with the small localised changes shown by Wang et al. (2017) when modelling cohesion. As suggested by Islam et al. (2020) the effect of particle-particle interactions should be investigated across many cases to understand their influence on respiratory drug delivery, but our study has provided a first step towards this. To use modelling to inform patient treatments an understanding of different particle forces and their change with drug properties is important.
Dosimetry analysis across the patients show noticeable differences in deposition in the trachea (Figure 9), attributable to transitioning from the fast-flowing, constricted throat to the expanding and curved trachea, of the cystic fibrosis and cancer patients. This agrees with single-phase simulations of Wei et al. (2017) that found flow to vary with orientation, as well as Bates et al. (2016) who showed the effect of pathological trachea curvature on gas pressure and energy loss. Also the change in relative orientation of the trachea to gravity present an additional factor to deposition changes. We also see a more uniform coating in the central airways of the healthy patient during exacerbating breathing, compared to normal breathing. As drug concentration across the airway surface is important in its dissolution into the tissue (Solomon et al., 2012), this provides a deeper understanding of the drug’s absorption than a simple deposition analysis. These results could be input to a model of the particle’s interaction with the tissue at smaller scales, as has been done by Olsson and Bäckman (2018) using 1D deposition data. Using CPFD for this could further inform clinicians in comparing potential treatment and delivery techniques for a patient.
When observing intra-patient drug delivery, the results presented show an uneven distribution across the lobes. There is a clear favour in transport to the lower lobes of the lung (Figure 7). This is similar to observations made by Lambert et al. (2011). Deposition is dominant in the right lung in the diseased patients, agreeing with experimental data (Usmani et al., 2005; De Backer et al., 2008). Patient-specific models could be used to attain a balanced dosage distribution, or target a specific region. This could allow for improved symptom relief through inhaler design informed by knowledge of regions sensitive to local inflammation (Barbu et al., 2011), allowing for more efficient devices. Additionally, cancerous regions of the lung could be targetted for chemotherapy using nebulisers (Tatsumura et al., 1993; Kleinstreuer and Zhang, 2003; Kleinstreuer et al., 2007). Clinicians could therefore predict and tailor the chemotheraputic agent delivery through patient-specific CPFD models to minimise radiation reaching undesired areas of the lung. As course of chemotherapeutic treatment is given over a large period of time (given in three week intervals by Wittgen et al. (2007)), this fits well with a high accuracy, time-consuming modelling method such as CFD-DEM.
Interestingly, the only lobar distribution changes in deep lung delivery occurred in the healthy patient as the left upper lobe received less medication with healthier breathing (Figure 7). Apart from this, the quantity of the dosage reaching distal lung regions increased with the longer, slower inhalation, but the distribution among lobes remained similar. This was also true for the cystic fibrosis patient (Figure 8b). This finding is beneficial as it suggests the drug’s behaviour does not differ significantly in extreme circumstances. However this is limited by a lack of knowledge of ventilation changes during exacerbations. Such information can be obtained by imaging combined with respiratory tests during healthy conditions (De Backer et al., 2010, 2014). To model deposition during an exacerbation, parameters usually obtained from imaging may need to be approximated in future studies. However, CPFD can be used to model exacerbations and avoid any dangerous experiments.
Limitations of this study included the absence of imaging data available for the throat and mouth of the child, meaning this region was taken from the other patients and scaled. This of course limits the specificity of the main deposition site. However by simulating with the two available throat geometries merged to the child’s trachea we can see that the deposition in this region itself is unaffected (only changing 1%). The only differences are observed downstream in trachea deposition (Figure 6) and deep lung dosage (Figure 7). Therefore, to neglect this region completely would harm the accuracy of the results downstream as flow generated in the upper airways heavily impacts downstream behaviour of both the gas and the particles (Figure 4). However, adding this region from adult patient to a child does take into account the maturational effects of adolescence on the airway, which may impact deposition here. In future studies, images containing the mouth and throat are needed to ensure this important section is patient-specific.
An additional limitation is that patients lying down in scans may have slightly different airway shape and orientation when standing or sitting upright (Jan et al., 1994), as is a typical posture when taking inhalers. If using image-based CPFD to recommend treatments, the patient images should be consistent with their typical posture when using their treatment where possible.
Only modelling the small number of airway bifurcation levels visible in the CT scan limits the model to only providing deposition information about the upper and central airways. To understand drug delivery in the smaller more distal airways the particles ‘exiting’ from our model’s outlets could be coupled to analytical 1D models (Koullapis et al., 2019). These models predict deposition based on particle size, estimated airway length and diameter, and flow rate. This would allow for a coarse prediction of drug delivery and efficacy in the targeted small airways.
Assumptions included the simulation of a uniform particle distribution. A uniform distribution is commonly used to research aerosol physics experimentally (Usmani et al., 2005, 2003), although does not directly represent the varying size distribution of inhaler particles (Dolovich, 1991). This is due to the small range of particle sizes (geometric standard deviation below two (Mitchell et al., 2003)) in a real device, meaning St < 1 in all cases, therefore the particles will exhibit high sensitivity to flow changes in both polydisperse and monodisperse flows. Deposition results should not be appreciably different as no transport characteristics are changed by this approximation. This allowed patient morphology differences to be evaluated in a simpler manner, without sacrificing accuracy.
A final modelling limitation is the softening of the particles. Particle softening is standard in DEM simulations to permit a larger DEM timestep as the particle’s collision happens over a longer period due to a lower stiffness. We mitigate an impact of increased van der Waals forces by including the cohesion model of Gu et al. (2016a). This model makes cohesive particle flows independent of their stiffness, and therefore increases our simulation timestep from tDEM = 0.8 ns to 25 ns. This reduction of real particle stiffness (kR) to a softer stiffness (kS) has been shown to have negligible effect on fluid hydrodynamics and particle cohesion (Gu et al., 2016a; Ozel et al., 2017). The use of this model therefore makes the effect of softening particles negligible. This allows us to use DEM to model inhaler particle-particle and particle-wall interactions with feasible computation times.
5. Conclusion
We have performed patient-specific inhaler deposition simulations across three diverse patients during exacerbating inhalation conditions. This has been compared to a healthy control patient and using a healthy and exacerbating breathing profile, showing that during an exacerbation less of the drug reaches the deep lung.
We have applied common means of improving drug delivery such as reduced particle size and a slower, steadier inhalation. Neither of these showed a common reduction in external airway deposition. This also shows the shortcomings in treating diverse populations with generic treatments. This demonstrates the need for personalised airways in respiratory CPFD studies, including the mouth and throat.
In the healthy patient the distribution of particles behaves as expected. This was fairly balanced across each lung, and primarily in the lower lobes of the lung (due to particle inertia). However we found particle distribution to be far less balanced in the diseased patients, as left to right lung deposition ratio was as low as 0.2 in their worst case. We predict this heterogeneity may be furthered upon inclusion of ventilation differences in the study to follow—particularly if patients studied suffer from issues such as mucus plugging or collapsed lung regions.
Data Availability
Simulation data available upon request.
Acknowledgements
For their early involvement and insight into the issue faced by patients, the authors would like to thank Elisabeth Ehrlich and Olivia Fulton. Also Carol Porteous (Patient Public Involvement Advisor) who arranged our contact.
The authors thank Prof. Vicki Stone from Heriot-Watt University for fruitful discussions. The authors also thank Dr Filippo Coletti from University of Minnesota for providing the surface file of the healthy patient’s segmented airways that enabled our model validation and was used to compare deposition across the patients. AO thanks Prof. Sankaran Sundaresan from Princeton University for his support and fruitful discussions. JW thanks Dr Rudolf Hellmuth from Vascular Flow Technologies for fruitful discussions.
JW was funded by an Institution of Mechanical Engineers Postgraduate Masters Scholarship 2018, a Scottish Funding Council Masters fee scholarship 2018, and a Carnegie-Trust for the Universities of Scotland PhD scholarship 2019.
Appendix A. Gas and solid-phase modelling
Simulations in this study solve fluid transport through the respiratory system on a Eulerian grid with Lagrangian particle tracking using the discrete element method (DEM). We solve the volume filtered mass (A.1) and momentum (A.2) continuity equations at each finite-volume cell. Here they are presented in terms of the volume filtered variables, accounting for particle interactions as derived by Anderson and Jackson (1967) and Capecelatro and Desjardins (2013), with overbars denoting filtered terms and bold, lower- and uppercase characters describing vectors and second order tensors, respectively,
Where subscript t is time, ū is the filtered local velocity vector, is the filtered local fluid-pressure, ρf is the fluid density, ϕ is the particle volume fraction, g is gravity (assumed constantly acting downwards at 9.81 N/kg), the force caused by interaction with the discrete phase is − Φd, Ru is the sub-grid stress from filtering, modelled using a dynamic Smagorinsky model (described in Section 2.2). is the filtered fluid stress tensor, composed of the fluid pressure gradient (), the deviatoric viscous stress tensor (labelled below), and an additional term arising from filtering of sub-grid velocity fluctuations (Rμ) (Capecelatro and Desjardins, 2013)
Where I is an identity tensor and Rμ is the term arising from filtering velocity gradients. In this study we dismissed Rμ to be included in a later study, and the deviatoric part of the stress tensor due to its minor influence on gas-solid flows in comparison to Ru (Agrawal et al., 2001). Φd is dependent on the interaction force between particles and fluid (ff→p,i) of all particles within a cell volume () by, (Anderson and Jackson, 1967; Capecelatro and Desjardins, 2013; Ozel et al., 2017). For a particle i, this is related to the filtered fluid stresses and the particle’s drag by , where fd,i is the drag force, taken from Beetstra et al. (2007).
Newton’s equations of motion are used to track particle linear motion (A.4) and angular motion (A.5) (Cundall and Strack, 1979; Capecelatro and Desjardins, 2013), by
Where mi is the mass of a particle, i, ν and ω are the particle’s translational and angular velocity, respectively, fc is the contact force from a particle-particle collision (subscript ij), and particle-wall collision (subscript iw), in the normal and tangential directions shown by sub or superscript n and t, respectively. van der Waals forces are shown with ν. Angular momentum (A.5) from inter-particle collisions depends on outward unit normal vector from particle centre to the point of collision, n, and the tangential contact force . Particle contact forces are solved here using a linear spring-dashpot model (Cundall and Strack, 1979; Capecelatro and Desjardins, 2013),
Where k, is the particle’s spring constant with kt = 2kn/7 (Matuttis et al., 2000), δ is the particle overlap distance, γ is the viscous damping coefficient. γn is calculated from the coefficient of restitution (Gu et al., 2016b). γt = 2γn/7. Effective particle mass is m* = mimj/ (mi + mj), in particle-wall collisions m* = mi, as one radius is assumed infinite (Gu et al., 2016a); μs is the sliding coefficient; tij represents the tangential displacement due to a collision, found from the integral of its velocity component.