Numerical simulations of targeted delivery of magnetic drug aerosols in the human upper and central respiratory system: a validation study

In the present study, we investigate the concept of the targeted delivery of pharmaceutical drug aerosols in an anatomically realistic geometry of the human upper and central respiratory system. The geometry considered extends from the mouth inlet to the eighth generation of the bronchial bifurcations and is identical to the phantom model used in the experimental studies of Banko et al. (2015 Exp. Fluids 56, 1–12 (doi:10.1007/s00348-015-1966-y)). In our computer simulations, we combine the transitional Reynolds-averaged Navier–Stokes (RANS) and the wall-resolved large eddy simulation (LES) methods for the air phase with the Lagrangian approach for the particulate (aerosol) phase. We validated simulations against recently obtained magnetic resonance velocimetry measurements of Banko et al. (2015 Exp. Fluids 56, 1–12. (doi:10.1007/s00348-015-1966-y)) that provide a full three-dimensional mean velocity field for steady inspiratory conditions. Both approaches produced good agreement with experiments, and the transitional RANS approach is selected for the multiphase simulations of aerosols transport, because of significantly lower computational costs. The local and total deposition efficiency are calculated for different classes of pharmaceutical particles (in the 0.1 μm≤dp≤10 μm range) without and with a paramagnetic core (the shell–core particles). For the latter, an external magnetic field is imposed. The source of the imposed magnetic field was placed in the proximity of the first bronchial bifurcation. We demonstrated that both total and local depositions of aerosols at targeted locations can be significantly increased by an applied magnetization force. This finding confirms the possible potential for further advancement of the magnetic drug targeting technique for more efficient treatments for respiratory diseases.


SK, 0000-0002-7568-5513
In the present study, we investigate the concept of the targeted delivery of pharmaceutical drug aerosols in an anatomically realistic geometry of the human upper and central respiratory system. The geometry considered extends from the mouth inlet to the eighth generation of the bronchial bifurcations and is identical to the phantom model used in the experimental studies of Banko et al. (2015 Exp. Fluids 56, 1-12 (doi:10.1007/s00348-015-1966-y)). In our computer simulations, we combine the transitional Reynolds-averaged Navier-Stokes (RANS) and the wall-resolved large eddy simulation (LES) methods for the air phase with the Lagrangian approach for the particulate (aerosol) phase. We validated simulations against recently obtained magnetic resonance velocimetry measurements of Banko et al. (2015 Exp. Fluids 56, 1-12. (doi:10.1007/s00348-015-1966-y)) that provide a full three-dimensional mean velocity field for steady inspiratory conditions. Both approaches produced good agreement with experiments, and the transitional RANS approach is selected for the multiphase simulations of aerosols transport, because of significantly lower computational costs. The local and total deposition efficiency are calculated for different classes of pharmaceutical particles (in the 0.1 µm ≤ d p ≤ 10 µm range) without and with a paramagnetic core (the shell-core particles). For the latter, an external magnetic field is imposed. The source of the imposed magnetic field was placed in the proximity of the first bronchial bifurcation. We demonstrated that both total and local depositions of aerosols at targeted locations can be significantly increased by an applied magnetization

Introduction
A rapid worldwide urbanization and the current development of mega-cities in Asia is associated with a significant deterioration of the air quality due to a fast increase of industrial and traffic pollution [1][2][3]. Especially alarming are the results of numerous studies directly connecting exposure to the longterm ambient air pollution, with a sudden increase in various respiratory diseases among children [4,5]. Inhalation of medical drug aerosols is applied for treating some of the most common respiratory diseases including asthma and chronic obstructive pulmonary disease (according to a World Health Organization fact sheet published in 2016, more than 235 million people have asthma and more than 64 million have chronic obstructive pulmonary disease). Similarly, due to numerous limitations associated with conventional treatments of lung cancer (which is one of the leading cancer killers in both men and women), alternative treatments based on the pulmonary route of drug administration directly to the lungs have recently been proposed [6][7][8].
Computer simulation-based studies can be a powerful tool in portraying airflow patterns and spatial distributions for different classes of inhaled particles within the human airways. Furthermore, such studies can improve different strategies to deliver inhaled medical drug aerosols to specific predefined locations within the human respiratory system to fight respiratory diseases [9][10][11][12]. Ultimately, such simulations can be applied to patient-specific conditions, providing the best strategies to achieve the most efficient delivery of the medical drug aerosols to diseased sites within the upper and central respiratory system. To achieve this goal, we would like to focus on the subject-and patient-specific human airway geometries.
A realistic extra-thoracic airway model was studied in [13]. The single representative geometry was selected from CT scans of five subjects. The final geometry included mouth cavity, larynx and trachea, without bronchial bifurcations. Different inhalation rates of 15, 30 and 60 l min −1 were considered with particle diameters ranging between 2 and 20 µm. The fluid flow and particle deposition analyses were performed by combining a low-Reynolds number k-ω model of [14] for the airflow and a Lagrangian stochastic trajectory method for the particulate phase. One of the major findings was that the laminar to turbulent transition was sensitive to the geometry of the airway model, especially for lower flow rates, stressing the importance of considering subject-specific geometries. It was shown that the major percentage of total deposition of particles occurred within the mouth cavity. The flow structures inside an idealized human upper airway were analysed in [15]. The lattice Boltzmann method (LBM) was applied to perform a direct numerical simulation study of the airflow patterns, which were compared with the X-hot-wire anemometry experiments of Johnstone et al. [16]. It was demonstrated that the LBM was in good agreement with experimentally observed flow features. The dynamics of airflow in a short inhalation (sniff) for an anatomically realistic geometry (starting from the nose to the second bronchial generation) was computationally studied in [17]. Additionally, the convective transport of a scalar species with relatively low diffusivity (with a Schmidt number of 900) was analysed.
The experimental studies of flow and deposition in models of human airways are necessary for validation of the computer simulations. The local and total deposition measurements in several realistic and throat geometries were performed in [18]. In total, seven representative geometries were studied at flow rates of 30 and 90 l min −1 , where the range of particles was between 3 and 6.5 µm. The geometry of the models was obtained from the magnetic resonance imaging (MRI) scans. The gamma scintigraphy and gravimetry were used to measure particle depositions. The authors proposed a simple equation correlating the total deposition of particles by introducing the equivalent diameter (D mean = 2 √ V/π L, where V is the total cast volume and L is the length of the central sagittal line of the model) and velocity (U mean = QL/V, where Q is the volume flow rate) in the Stokes number, which significantly reduced scattering of the data. A steady inspiratory flow, in an anatomically representative model of the human airways, was experimentally investigated in [19]. Magnetic resonance velocimetry (MRV) was used to measure three components of the mean velocity, where water was used as a working fluid. The subject-specific geometry was manufactured from the CT scan images by three-dimensional (3D) printing (stereolithography). It was concluded that both streamwise (due to streamwise gradients) and lateral dispersion (due to secondary flows) were relevant transport mechanisms. This investigation also stressed the importance of considering the subject-and patient-specific geometries instead of simplified airway configurations. In their follow-up study, time-varying flow conditions, which correspond to breathing regimes during moderate exercises, were also analysed with the same technique [20].
In the present study, in addition to a passive distribution of inhaled aerosols, we would like to analyse the potential of a controlled drug-aerosol delivery in the anatomically realistic human respiratory system, based on the magnetic drug targeting (MDT) concept. Magnetic drug delivery involves a straightforward concept of imposing external magnetic field gradients in the proximity of the disease location to act upon the magnetically responsive carriers. An active pharmaceutical ingredient is attached to these magnetic carriers, making it possible to deliver required therapeutic concentrations at pre-specified locations.
Some pre-clinical and phase I/II clinical in vivo trials of the MDT technique were reported for treatments of advanced solid malignant tumours [21][22][23], liver cancer [24,25], prostate cancer [26] and breast cancer [27,28]. Despite these significant advancements of applying the MDT technique, there are still numerous challenges that need to be solved before extensive use of magnetic drug delivery [29][30][31]. The primary challenge includes a limited exposure of deep organs to sufficiently high magnetic fields to enable an efficient capture of medical drug particles. In our previous work, we demonstrated how computer simulations could contribute to further advancement and optimization of the MDT technique in dealing with localized treatments of cardiovascular diseases [32,33] and brain tumours [34]. In the present contribution, we extend further applications of the MDT technique for potential treatment of respiratory diseases, specifically for localized treatment of a lung cancer.
The motivation of the present research is triggered by the experimental study of the targeted delivery of magnetic aerosol droplets to intact mice lungs, as reported in [35]. It was demonstrated that steering magnetic droplets towards the right lung could be achieved by imposing sufficiently high magnetic field gradients close to the main trachea bifurcation. This magnetic steering enabled an overall fourfold increase of medical drug deposition compared to the neutral situation. In conclusion, Dames et al. [35] suggested further developments and scaling-up of the magnetic steering approach to the human lungs. The primary goal of the present research is to check, by performing a series of numerical simulations, to what extent (if at all) it is possible to achieve a magnetic control of drug aerosol distribution in the human airways.
In this work, we focus on computer simulations of the MDT technique to steer and capture magnetic aerosols in an anatomically realistic geometry of human upper and central airways (that extends from the mouth inlet to the eighth generation of bronchial bifurcations) for which a complete 3D velocity field was measured by the MRV technique [19,20]. First, we performed computer simulations of the airflow inside airways to compare the mean flow features with the MRV measurements. Here, the major challenge lies in the local generation of turbulence caused by the formation of the wall jet in the larynx region. Then we analysed the distribution of various classes of the neutral aerosols inside the different parts of the airway system. Finally, aerosols with magnetic core were introduced and a few scenarios of steering and deposition by imposed external magnetic field gradients were analysed.

Mathematical model
In this section, we provide a comprehensive set of equations used to describe the motion of air and the dynamics of the particulate phase. The airflow dynamics equations are based on the Eulerian approach, whereas the dynamics of the particulate phase are represented in the Lagrangian frame of reference.

Airflow
For the airflow, conservation of mass and momentum in a Cartesian coordinate frame of reference were used: and where ρ is the density of fluid and ν its kinematic viscosity. As we are dealing with a flow where a local generation of turbulence can occur, we have to account for the turbulence contributions through the turbulent viscosity (ν t ). In the present study, we apply two different approaches to model turbulence.  [36] adopted here, the following set of four transport equations is used, i.e. the intermittency (γ ), the transition momentum thickness Reynolds number (Re θt ), the TKE (k) and, finally, the specific dissipation rate (ω). The second simulation approach is the large eddy simulation (LES), which requires that the simulation is run in a time-dependent manner, despite the steady inhalation conditions. Here, the numerical resolving of the flow dominates over its modelled part, which is contained only to flow and turbulence structures smaller than the defined filter size (a size of control volume). Here, we adopted the wall-adopting local eddy-viscosity (WALE) subgrid model of [37] to calculate the subgrid contributions.
2.1.1. Reynolds-averaged Navier-Stokes: SST-transition turbulence model The transport equation for the intermittency is formulated as follows: where the main transition source (P γ ) and dissipation (E γ ) terms are defined as where S is the strain rate magnitude, Ω is the vorticity magnitude, while F length , F onset and F turb represent intermittency functions, which are designed to properly capture transition from laminar to turbulent flow regime (for detailed derivations of these functions and their finite forms, see [36]). Finally, the remaining model constants are given as σ γ = 1, c a1 = 2, c a2 = 0.06, c e2 = 50. The transport equation for transition momentum thickness Reynolds number (Re θt ) can be written as The main source term (P θt ) is defined as where the Re * θt is the local value calculated from an empirical correlation, t = 500ν/U 2 is the characteristic local time scale, F θt is the blending function (for the full expression see [36]) and, finally, c θt = 0.03 and σ θt = 2.0 are model coefficients. The transport equation of the TKE (k) is given as Here, P k and Y k are the production and dissipation terms of the turbulence kinetic energy and γ eff is the effective intermittency calculated from: where Re v = y 2 S/ν, R T = k/ω and Re θc are the local strain-rate Reynolds number, the turbulent Reynolds number and the critical Reynolds number where the intermittency first starts to increase in the boundary layer, respectively. Finally, the transport equation for the specific dissipation rate (ω) is given as with P ω and D ω as the production terms and Y ω as the dissipation term. The system of equations (2.1)-(2.9) is finally closed with the expression for the turbulent viscosity, which is written as where a 1 = 0.31 is the model constant and F 2 is the blending function (for definition see [14]).

Large eddy simulation: Wall-adapting local eddy-viscosity subgrid model
Within the LES framework we applied the WALE subgrid model in [37]. This model is selected because of the fact that, in contrast with the standard Smagorinsky model, it does not require any additional empirical wall-damping functions and provides a proper near-wall behaviour of the subgrid turbulent viscosity in the proximity of the wall.
whereŪ i is the resolved velocity, κ is the von Kármán constant, d is the distance to the nearest wall, = V 1/3 is the characteristic local grid scale (grid filter) and, finally, C w = 0.325 is the model constant.

Aerosol/particulate phase
The Lagrangian approach is applied for the particle tracking and for deposition of particles. This approach is based on a simple balance of all forces (F) acting on the particle [38]: where m p is the particle mass and t is the characteristic time. After integrating the particle velocity with respect to time, the particle position is calculated as We next consider the drag, virtual mass, gravity and magnetization forces. The drag force is expressed as where u is the velocity of the fluid phase and f is the drag factor, which depends on the particle Reynolds number, defined as and τ p is the particle relaxation time, defined as In the present study, we assume that the slip correction factor C c = 1. We define the Stokes number as the ratio between the characteristic particle relaxation time to the time of the flow: whereŪ is the characteristic velocity and D is a typical length scale of the fluid. The drag coefficient is specified as a function of the Re p after [38]: The virtual mass (added mass) and gravity forces are evaluated as Finally, the magnetization force can be written as  Figure 1. The front view of the geometry of the patient-specific human airway system that extends from the mouth inlet to the eighth generations of the bronchial bifurcations [11,19,20] (a); details of the numerical mesh used (zoom-ins)(a); a side view of geometry with imposed boundary conditions: a single inlet and multiple outlets (in total 72), where a zero gauge pressure condition is imposed (b); details of the numerical mesh (zoom-ins) at 1-2-3 locations (c), (d) and (e), respectively, as indicated in figure 1a.
where V mp is the volume of the magnetic core of the particle, μ 0 is the magnetic constant, M is the magnetization and H is the auxiliary magnetic field. In the present study, we are dealing with fully magnetically saturated particles for which magnetization reduces to M = M sat H/|H|, where M sat is the saturation magnetization (a material property). We also define characteristic non-dimensional parameters to quantify the particle deposition efficiency. The local deposition efficiency, defined as the dimensionless surface concentration of the particles, is calculated as where is the typical search radius from the centre of the wall control volume (e.g. ξ 1 mm defines that the search radius is = 1 mm), N in is the number of particles injected at the inlet and N wall is the number of the particles captured within the pre-specified search radius ( ) [34].

Numerical method
For the gaseous phase (air), the system of equations (2.1)-(2.12) is numerically solved by using the finite-volume approach for arbitrary complex geometries, Ansys/Fluent Inc. 14  in transport equations are discretized by the second-order central differencing scheme (CDS). For the RANS-SST model, the second-order upwinding differencing scheme is used for convective terms in all transport equations. By contrast, the second-order CDS is used in the LES approach for convective terms. Besides, the fully implicit second-order time integration scheme based on the three consecutive time steps is applied in LES. The SIMPLE algorithm is used for predictor/corrector coupling between the velocity and pressure fields. For the particulate phase, equations (2.13) and (2.14) are integrated by the fourth order Runge-Kutta time integration method. The separate user-defined functions are developed which include calculations of the spatial distribution of the auxiliary magnetic field and resulting magnetization force. The former is based on the integration of Biot-Savart's Law for a current-carrying wire (e.g. [39,40]).

The geometry of the upper and central airways and the numerical mesh used
The geometry used is, apart from the inlet region, identical to that of [11] and which was also experimentally studied in [19,20], figure 1a,b. The CAD geometry was reconstructed from the CT scans with a resolution of 1.25 mm, which was received from [19]. The geometry starts with a part of the mouth and ends with the eighth bronchial generations in the lungs. An Octree-type mesh was selected. The surface mesh was generated first. Then, this smoothed surface mesh was extruded inwards by placing between five (coarse mesh) and 10 (fine mesh) segments to create the prismatic boundary layer. This local mesh refinement is used to more accurately handle possible flow separation and deposition of particles in the near-wall regions. The remaining part of the geometry was then filled using the Delaunay meshing method. Some details of the finest mesh containing approximately 11 × 10 6 control volumes and 10 prismatic segments within the boundary layers are shown in figure 1c-e. All results shown are for the finer mesh.

Imposed boundary conditions
The imposed boundary conditions are also given in figure 1b. The walls are treated as the no-slip type of boundary condition. The inlet velocity profiles were assumed to be uniform with the velocity magnitude obtained from the pre-defined value of the Reynolds number, i.e.  [19]). For the RANS-SST model, the inlet turbulence intensity is assumed to be 5% and that the ratio between the turbulent and molecular viscosity is ν t /ν = 10. The zero gauge pressure boundary condition was applied to all outlets (here 72 outlets were imposed in total).

On the adequacy of the numerical mesh used
Because of the expected turbulent character of the flow within the geometry and specified inspiration rate considered, it is important to check to what extent the generated numerical mesh is adequate for the low-Reynolds turbulence closure (within the RANS framework) or for the even more demanding (regarding the spatial mesh resolution) LES approach. The first criterion is evaluated from the low-Re RANS-SST results by plotting the non-dimensional wall distance calculated as y + = u τ · y n /ν, with the friction velocity u τ = √ τ wall /ρ and y n as the wall-normal coordinate. The wall-shear stress (WSS) is calculated as τ wall = |μ(∂u/∂y) wall |. The contours of the non-dimensional wall distance (y + n ) are shown in figure 2a-c. It can be seen that the maximum value of y + n ≈ 2 can be observed at particular locations, whereas the criterion that y + n ≤ 1 holds for a large part of the simulated domain. This confirms that the numerical mesh is sufficiently refined in the near-wall region. Similarly, in our preliminary study [33], we also estimated the ratio between mesh-based and Kolmogorov length scales (l ratio = V 1/3 mesh /η K ). Here, the V mesh is the volume of the characteristic mesh element, whereas the Kolmogorov length scale is evaluated as η K = (ν 3 /ε) 1/4 (ε from the low-Re RANS-SST is used). For a well-resolved LES, it is to be specified that this ratio is between 5 and 10. We found that apart from a relatively narrow range of 0.13 ≤ z ≤ 0.15 m, this condition was also satisfied. It can be concluded that the carefully designed locally refined numerical mesh with prismatic boundary layers is appropriate for the well-resolved LES.  as well as its overall shape. Compared to experiment, the RANS-SST simulation slightly underpredicts the strength of the jet in 0.135 ≤ z ≤ 0.145 m region. Consequently, the region downstream from the laryngeal jet also shows some differences. It can be concluded that good agreement between simulations and experiment is obtained concerning predictions of the laryngeal jet behaviour. Two selected isosurfaces of the turbulent kinetic energy (TKE) are shown in figure 6. Note that experimental data do not provide this quantity, so only simulations results are compared. The distribution of the TKE nicely portrays the local generation of turbulence within the laryngeal jet and its relatively rapid decay in the downstream region. This distribution corresponds to regions where the most dense population of vortical structures is observed, as shown in figure 4. It can be seen that both simulations capture similar regions. The LES results (figure 6a) indicate stronger contributions of the velocity fluctuations compared to the RANS (figure 6b) results, but the overall agreement is good.
Next, we move closer to the regions just after the first bifurcation; figure 7. Here, we compare measured and simulated (LES) non-dimensional z-velocity components. The velocity contours show good agreement for the incoming jet in the 0.25 ≤ z ≤ 0.3 m region. Then, flow splits asymmetrically with a distinctly higher intensity of velocity in the left branch. In the left branch, the LES shows slightly higher velocity compared to measurements in the z > 0.34 m region. In the right branch after bifurcation, the flow behaviour is well captured. This includes the recirculation zone in the 0.32 ≤ z ≤ 0.33 m, as well as the jet deflection towards the left/lower wall within the right branch.
We next analyse distributions of the velocity and vorticity in selected horizontal cross sections (i.e. in the x-y planes perpendicular to the main flow direction), figures 8 and 9. Here, we focus mainly on the pharynx and post-laryngeal regions ( figure 1a,b), where the flow exhibits the most complex behaviour, which includes helical motion (swirling), separation, strong acceleration and recirculation. The four pre-selected planes (plane (1), plane (3), plane (4) and plane (6)) are defined with z = 0.1, 0.135, 0.146 and 0.161 m (figure 5), respectively. It can be seen that a good agreement between the experiment and simulations (LES) is obtained at practically all locations. All experimentally captured (figure 8a,c,e,g) salient flow features are also visible in simulations ( figure 8b,d,f,h). There is just a slight disagreement in plane (6), where the simulation shows a somewhat more extended left-side-oriented jet-like structure compared to the experiment. The contours of the non-dimensional vertical vorticity (ω * z = ω z · D/U ref , ω z = −(∂u x /∂y − ∂u y /∂x)) in identical planes (plane (1)-plane (6)) are shown in figure 9. Note that the vorticity is a more sensitive parameter relative to the velocity magnitude because it contains derivatives of both velocity components in the particular plane. It can be seen that overall behaviour is relatively well captured in simulations, but there are some minor differences at specific locations. Differences in the proximity of walls can be explained by the lower spatial resolution of the measurements, which is particularly important in the near-wall regions.  (1) plane (3) plane (4) plane (     Breuer et al. [42], R 0 = 5.6 Cohen Stuart et al. [40], R 0 = 5.6 present-random inlet, R 0 = 5.6 present-uniform inlet, R 0 = 5. 6 Pui et al. [41], steel, R 0 = 5. 6 Pui et al. [41], steel, R 0 = 5. 7 Pui et al. [41], glass, R 0 = 5.7 (b) (a) Figure 14. Benchmark solutions for the 90 • bend geometry validated against numerical simulations of [40,42] and experiments of [41].   x (m) Finally, based on all the presented comparisons between the water model experiment and simulations (figures 5-13), it is concluded that generally a good agreement is obtained at all locations considered. The LES results proved to be in a slightly better agreement with experiment than the RANS-SST. Despite this, considering the total computation costs between the LES and RANS-SST used here (the LES are at least O(10 2 ) more computationally expensive due to a necessity to perform the time-dependent simulations even for the steady inspiration conditions) and that the RANS-SST model properly captured practically all of the most salient flow features in the upper human airways, we decided to use the RANS-SST model in the context of the multiphase simulations that will include Lagrangian tracking of injected aerosols.

Multiphase simulations of aerosol transport and deposition
Prior to full multiphase simulations of aerosol distribution in the human upper and central respiratory system, we performed a series of simulations in simplified geometries for which experimental and numerical results of other authors are available. The first configuration is a 90 • bent tube with characteristic curvature ratio (R 0 = R bent /R tube = 5.6 giving a Dean number of De = 423) and with a fully developed laminar inflow (Re D = 1000) with non-magnetic particles, previously experimentally studied in [41] and numerically in [40,42]. Here, we study effects of the generated secondary motion on the particle distributions. The sketch of the geometry under study as well as specification of different classes of particles (with a particle diameter range 20 ≤ d p ≤ 85 µm) are shown in figure 14. It is shown that the total deposition efficiency agrees very well with experiments and with other numerical studies for all simulated classes of particles, confirming the accuracy of the present approach. In the second simplified geometry case, we address flow and particle deposition in a double bifurcation geometry, identical to [43]. Despite simplifications, this geometry contains many flow features occurring in real human airways. Here, the air is a working fluid and particles have a diameter of d p = 10 µm and a density of ρ p = 1060 kg m −3 . It can be seen that both velocity distributions and cumulative deposition efficiency agree very well with results of [43] ( figure 15). We conclude that the present results are in good agreement with the two above-mentioned benchmark studies in predicting both integral and cumulative deposition efficiencies, and as a next step, we simulate the geometry of the human respiratory system  (as shown in figure 1). It is noted that here we combine the RANS-SST model for the airflow and the Lagrangian tracking for the particulate phase. The particles are uniformly released at the inlet plane. As a very diluted concentration of the aerosols is inhaled, we assume one-way coupling between particles dynamics and airflow. Furthermore, the particle-particle interactions are neglected. This makes it possible to simultaneously simulate various classes of particles in a single computational run. The particle-wall interactions are assumed to be fully inelastic. This is due to the presence of the adhesive mucous layer surrounding the airway walls. To  figure 16. These plots are generated by recording locations where the aerosol particles are crossing the specified horizontal planes. It can be seen that for the smallest particles (d p = 0.1 µm) strongest scattering is observed, whereas the larger particles (d p = 3 and 5 µm) are confined within significantly smaller regions. For all classes of particles, a strong asymmetry in distribution is caused by the presence of substantial helical airflow pattern within the trachea. The contours of the local deposition efficiency (ξ ) with a search radius of = 1 mm, for different classes of the neutral (nonmagnetic) particles (d p = 0.1, 1 and 5 µm), are shown in figure 17. It can be seen that the particles deposit primarily at the back of the pharynx (for all d p ) and at bifurcations (for particles with smaller d p = 0.1 and 1 µm). The local concentration of deposited particles at the back of the pharynx increases with increasing a particle diameter; figure 17d-f.
In the next step, we apply an external magnetic field and investigate its impact on the behaviour of magnetic particles. As a proof of concept study, we consider two situations, which are characterized with different distances from the point-like magnetic source as illustrated in figure 18. The magnetic source is placed at the right side from trachea just above the first bifurcation; figure 18a. The characteristic distances between the magnetic source and the wall of the trachea are 1 cm and 10 cm, respectively.  Figure 19. The contours of the local deposition efficiency (ξ 1 mm ) for two classes of the magnetic core-shell particles with d p = 1 µm ((a)-(c)) and 5 µm ((d)-(f ))-both with d m /d p = 0.84. No magnetic field ((a) and (d)), the magnetic field in case with r mag = 1 cm ((b) and (e)) and r mag = 10 cm ((c) and (f )) distance from the magnetic source, as shown in figure 18.
field in the present generation of MRI scanners in clinical applications. A colour map of the spatial distribution of the imposed magnetic field for these two distances are plotted in figure 18b,c, respectively. For the two situations, the strength of the electric current through the point-like (wire) source is the same, resulting in a magnetic field strength along the trachea of 2 T and 0.2 T, for distances of r mag = 1 cm and r mag = 10 cm, respectively. In the present investigation, we define two classes of magnetic particles. We introduce a characteristic ratio of the magnetic-core (d m ) and total (d p ) diameters, defined as d * mp = d m /d p . When d * mp = 1, we have fully magnetic particles, whereas, for d * mp = 0, we have magnetically neutral particles. For the magnetic drug delivery, we introduce particles containing the magnetic core. The shell is a mixture of a medical drug and a carrier. For the carrier, we adopt the biodegradable poly(lactide-co-glycolide), commonly  referred to as the PLGA. The resulting density of such a shell particle can be calculated as As the shell contains both magnetic carrier and medical drug, its density depends on the medical drug loading. The latter is expressed as a mass percentage of the medical drug in the shell layer. Finally, the shell density is calculated as where f l is the fractional loading of the medical drug within the shell layer (based on the mass). In the present study, the magnetic core is made from the iron oxide-maghemite (γ Fe 2 O 3 ), with a density of ρ = 4860 kg m −3 [44], magnetic susceptibility χ m = 3 [45] and saturation magnetization M sat = 3.9 × 10 5 A m −1 [46]. The shell carrier is PLGA with a characteristic density of 1300 kg m −3 . The fractional loading of the medical drug is taken to be f l = 30%, which corresponds to an anti-tubercular drug, with a typical density of 1610 kg m −3 [47]. The resulting local deposition efficiency of the shell-core magnetic particles is shown in figure 19. For illustration, we selected the core-shell particles with d p = 1 µm and 5 µm, and with d * mp = d m /d p = 0.84. Compared to the neutral situation shown in figure 19a,d, contours indicate a significant enhancement of the particle deposition efficiency in the proximity of the magnetic sources, as shown in figure 19b,c,e,f. With increasing distance between the magnetic source and the targeted regions, from r mag = 1 cm to a more realistic distance of r mag = 10 cm for the entirely non-invasive treatment (figure 19c,f ), the local deposition efficiency is still increased in comparison to the neutral case.
The total deposition efficiency (which is defined as a ratio of the number of the particles deposited along the walls and the number of the injected particles) for the neutral and for two classes of the magnetic particles (d * mp = 1 and 0.84, respectively) for the previously analysed two different locations of the magnetic sources, is plotted in figure 20. The experimental results of [9] are also added for a qualitative comparison. No significant effect of the magnetization can be observed in the range St ≤ 5 × 10 −4 . For the St > 0.1 range, practically all particles found are deposited along the walls. The dependence observed for the non-magnetic particles shows a qualitatively similar behaviour to the results of [9]. The magnetic modulation is the most effective in the 5 × 10 −4 ≤ St ≤ 10 −1 range. This range corresponds to particles with a diameter approximately between 0.3 and 5 µm. For example, at St = 0.05, the total deposition efficiency has increased from 23% for the neutral particles to 46% and 85% for the magnetic particles for the magnetic source distances of r = 10 cm and 1 cm, respectively. Similar analysis for St = 0.005 shows that the efficiency increased from a rather low 2% for neutral particles to 5% and 15% for magnetic sources 2 and 1, respectively. It can be noted that relatively small differences are obtained between the fully magnetic (solid lines) and core-shell particles (dashed lines) confirming potentials of the MDT technique to deliver medical drugs coated around a magnetic core. This kind of information, together with maps of the local concentrations of the deposited particles, can be used for improvements of the existing generation of magnetically neutral medical aerosols, as well as for the design of a new generation of medical drugs with the magnetic core.

Conclusion
We applied numerical simulations to describe the mean airflow patterns and distributions of aerosols in an anatomically realistic geometry of the human airways under the steady inspiration condition with a flow rate of 60 l min −1 . The geometry considered extended from the mouth inlet to the eighth bronchial generations. The results obtained with RANS-SST and LES-WALE were compared against MRV experiments under identical conditions. A good agreement between simulations and measurements was obtained at all locations considered. Despite the fact that a slightly better agreement with experiments was obtained with the LES compared to the RANS, the Lagrangian tracking of the particulate phase was performed in conjunction with the RANS-SST model due to significantly lower computational costs. We covered an extensive range of particle sizes, ranging from 0.1 to 10 µm, without and with the magnetic core. For the latter, to mimic as realistically as possible the structure of the pharmaceutical medical drugs, we have selected some of the recently proposed aerosol carriers based on the shell-core concept. We found that the local and total deposition efficiencies can be significantly enhanced by activation of the magnetization force. The most effective enhancement was observed in the 5 × 10 −4 ≤ St ≤ 10 −1 range, which corresponds to particles with a diameter approximately between 0.3 and 5 µm.
Data accessibility. The data supporting the findings of this study are available at the Dryad Digital Repository (http://dx.doi.org/10.5061/dryad.0jt43) [48].