- Email: [email protected]
Transport and fate of inhaled particles after deposition onto the airway surface liquid: A 3D numerical study
Computers in Biology and Medicine 117 (2020) 103595
Contents lists available at ScienceDirect
Computers in Biology and Medicine journal homepage: www.elsevier.com/locate/compbiomed
Transport and fate of inhaled particles after deposition onto the airway surface liquid: A 3D numerical study Shayan M. Vanaki a ,∗, David Holmes a ,∗, Kabir Suara a , Pahala Gedara Jayathilake b , Richard Brown a a b
School of Chemistry, Physics and Mechanical Engineering, Queensland University of Technology, Brisbane, Queensland, Australia Department of Oncology, University of Oxford, UK
ARTICLE
INFO
Keywords: Therapeutic drug particles Particle tracking Particle deposition Cilia Periciliary liquid layer Respiratory diseases
ABSTRACT In this paper, a numerical investigation is carried out to provide insights into the fate of inhaled aerosols after their deposition on the lung lining fluid in both healthy and diseased states. Pulmonary drug delivery is a well-known non-invasive route of administration compared to intravenous delivery. Aerosol particles are formulated and used as drug carriers, which are then sent to the airways using aerosol drug delivery devices. This approach is useful for site-specific treatment of lung diseases, treatment of central nervous system (CNS) disorders and a variety of other diseases. Bioavailability of the inhaled therapeutic particles after landing on the airway lining fluid can be significantly altered by the lung muco-ciliary clearance, a process through which hairlike structures known as cilia beat in a harmonised manner and induce the mucus in the proximal direction, leading to an effective clearance of the foreign inhaled particles entrapped by this sticky layer from the airways. Here, we set up a 3D computational model of ciliary arrays interacting with periciliary liquid film (i.e. confined between the epithelium and mucus layer) and a detailed analysis is conducted to better understand the fate of drug nanoparticles that are able to penetrate the mucus. Consistent with clinical findings, we find that the actions of cilia result in a low rate of drug retention and absorption by the pulmonary tissues in healthy lungs. However, under conditions associated with abnormal ciliary beats, the retention time of particles is notably increased at the site of release. Nonetheless, the results associated with some of the ciliary impairments reveal that deposition of drug aerosols on the ciliated cells may still be a significant challenge. These findings have potentially important implications on the modification of therapeutic drug particles to achieve a higher absorption rate.
1. Introduction Computational Fluid Dynamics (CFD) is garnering a great deal of attention in the medical context, as it can be used as a powerful tool to provide valuable insights into physical phenomena occurring inside the human body without the need of expensive clinical experiments and the associated ethical challenges. As such, many researchers have used CFD methods to study the transport of biofluids through various human organs, such as the heart [1], brain [2], and stomach [3]. Lung modelling has also been of interest for many researchers over the years. Modelling inhaled particle deposition in the lung airways has already received a lot of attention from scientists [4]; see e.g. Hofmann [4] and Longest et al. [5] for reviews. However, the clearance process of aerosols from the pulmonary airways has been more challenging to investigate numerically due to the multi-scale, multi-physics, and multi-components involved. In particular, the conductive airways
are protected by a lining fluid called ‘airway surface liquid’ (ASL) which is a multiphase fluid composed of two layers [6]: (i) an aqueous periciliary liquid layer (PCL) that resides on the airway epithelium and contains a carpet of motile microstructures called cilia, and (ii) a nonNewtonian, viscoelastic mucus layer [7] which captures the foreign particles (e.g. dust, PM, bacteria, etc.). The harmonised beating of cilia within the PCL propels the mucus layer continuously towards the trachea and eventually out of the lungs. This fluid–structure interaction (FSI) process is known as ‘muco-ciliary clearance’ which is the primary defence mechanism, with 80%–90% of inhaled particles being discharged from the upper and central lung within 24 h [8,9]. Some works in the literature have simplified such a complex clearance process by using theoretical assumptions, where the transport of particles within the mucus layer is merely estimated at every airway generation, without taking the interaction of cilia with the PCL or mucus layer
∗ Corresponding authors. E-mail addresses: [email protected], [email protected] (S.M. Vanaki), [email protected] (D. Holmes).
https://doi.org/10.1016/j.compbiomed.2019.103595 Received 2 October 2019; Received in revised form 14 December 2019; Accepted 27 December 2019 Available online 7 January 2020 0010-4825/© 2020 Elsevier Ltd. All rights reserved.
Computers in Biology and Medicine 117 (2020) 103595
S.M. Vanaki et al.
Fig. 1. Fate of an aerosol drug after deposition in the lungs. (1) Release of the active nano-sized pharmaceutical ingredient from the micro-sized aerosol particle deposited on the mucus gel. (2) Absorption of the active ingredient across the pulmonary epithelium into the bloodstream. (3) Clearance of the undissolved particle. Source: Figure adapted from [36].
into account; see e.g. [10–15]. Numerical modelling approaches that explicitly account for cilia PCL interaction, including continuum cilia models [16–18] as well as discrete cilia models (comprised of prescribed ciliary beating models [19–28] and FSI models [29–31]) have been developed, potentially allowing for more detailed investigation into the complex physics involved in particle transport. None of the above-mentioned works have investigated the transport of particles across the lung lining layer, which is of importance, particularly in the context of aerosol medicine and pulmonary drug delivery, and the focus has been limited to fluid mixing (e.g. [27,32]) and particle capture by beating cilia (e.g. [33–35]). Dissolution of active pharmaceutical ingredients (or carrier particles) in the ASL and their absorption into the bloodstream would be desirable for pulmonary specialists. However, as shown in Fig. 1, the persistence and subsequent absorption of the inhaled drug carrier in the lung is significantly influenced and reduced by the muco-ciliary clearance process, thus substantially minimising the impact of therapeutic drug particles [36]. Nanoparticles intended for use as drug carriers can permeate through the pores of the mucus network [37] and reach to the underlying PCL, leading to absorption across epithelial and endothelial cells. Very little is known about the fate and transport of these nanodrugs after deposition and penetration inside the mucus layer [36], especially in the presence of an active biological factor such as cilia. The goal of this paper is, therefore, to understand the fate of drug carriers in a flow generated by a periodic array of beating cilia in the lung. Also, investigated is the transport of such particles in diseased lungs, due to abnormal ciliary beatings, compared to that in healthy states. To achieve these objectives, in Section 2, we use a computational model based on the projection method for solving a three-dimensional Navier–Stokes equation and a particle tracking method. In Section 3, we track and quantify the transport of inhaled ultrafine particles in the periciliary liquid film under both healthy and diseased lung conditions. Finally, Section 4 concludes this paper and future numerical challenges and directions are provided.
Fig. 2. Particle tracking verification over a span of 43 s. (a) Present model. (b) Experimental work of Tsorng et al. [38].
the motile cilia and fluid is computed explicitly through the immersed boundary (IB) method, analogous to the one-way FSI model presented in our previous work [39] and hence not repeated here. By feeding the obtained velocity fields from the FSI model into the Lagrangian particle tracking model, a detailed numerical study is then carried out to better understand the fate of drug nanoparticles that are able to penetrate the mucus. 2.1. Particle tracking The motion of finite-inertial spherical particle in fluid can be described by the full four dimensional Maxey–Riley (MR) ordinary differential equation for particle position [40,41]. The inertial particle positions as well as initial particle velocities which are generally unknown are required to integrate the MR equation. However, using the slow manifold approximation, Ref. [42] showed that asymptotic motion of finite-sized particles can be described by a two-dimensional system in the position, that is
2. Numerical methods In the present work, the PCL flow is modelled and solved using a finite-difference projection algorithm, and the interaction between
𝐱̇𝐩 = 𝐮 + 𝜖 2
𝛥𝐮 , 𝛥𝑡
(1)
Computers in Biology and Medicine 117 (2020) 103595
S.M. Vanaki et al.
sizes comparable to cilia movement, the difference between the motion of the fully resolved particles and that of passive tracers is small. To verify the appropriateness of the particle tracking method, we track a particle trajectory in a 3D lid-driven cavity numerically and compare the result with that observed in the experimental work of Tsorng et al. [38]. Since fluid advection is the dominant mechanism for the transport of particles within the cilia-induced PCL flow [27,32], it is still reasonable to have this comparison with an advection dominant cavity problem where Re = 470. The simulated trajectory was obtained by integrating Eq. (3) from the initial location of the observation over the time span of the observation. The comparison is presented in Fig. 2. The simulated and observed trajectories are similar for a major part of the time span with some differences occurring at the end of the tracking. Eq. (3) is solved using the Runge–Kutta (4,5) method, i.e. the MATLAB ode45, to obtain particle trajectories. 2.2. Numerical setting and boundary conditions Fig. 3. Schematic of a row of cilia to illustrate how a ciliated domain can be replicated. For simplicity, it is assumed that the shape of each cilium changes from phase 1 to 4 through one beat cycle. Extra cilia are added to each side of the ciliated domain such that the cilium with the phase number of 4 has the one with phase number of 1 on its right side at the beginning of each beat cycle. This would allow for replicating the ciliated domain and removing the discrete boundary effects. Note that adding extra cilia to any side of the ciliated domain would require extension of the periodic boundary conditions, from 𝐵𝐶𝑖𝑛𝑖𝑡𝑖𝑎𝑙 to 𝐵𝐶𝑓 𝑖𝑛𝑎𝑙 .
For all of the simulation cases presented in this paper, the fluid grid size is ℎ = 0.5 μm, the time step of the computation is 𝛥𝑡 = 1 × 10−6 s, the fluid dynamic viscosity is 𝜇 = 0.01 g∕cm.s, the fluid density is 𝜌 = 1 g∕cm3 , and the number of control points along the cilia is set at 20. As for the boundary conditions, periodic boundaries are applied in the wall-parallel directions. The bottom boundary at the epithelium is treated as a no-slip wall in our simulations consistent with the velocity field at the epithelium boundaries of cultured mouse tracheal cells [45]. Given that the boundary condition at the PCL-mucus interface is influenced by the surface tension force between the mucus layer and PCL as well as the complex fluid rheology of the mucus, it is reasonable to assume the top boundary as a free slip wall in the present mucus-free model. To allow the tracking of the dispersed particles over a longer period of time, we perform the following steps: (i) once the simulations reach the steady state condition, the stable (or the last) cycle is replicated 𝑁𝑐𝑦𝑐𝑙𝑒 times until it reaches at 𝑡𝑒𝑛𝑑 = 𝑁𝑐𝑦𝑐𝑙𝑒 ×𝑇 ; (ii) the physical domain is extended in the 𝑥- and 𝑦- directions through boundary replication. To achieve this and remove any boundary effects, extra rows/columns of cilia with specific phase numbers are added to the each side of the square ciliated domain such that the selected periodic domain for particle tracking does not experience disturbance when it is replicated. Fig. 3, for example, illustrates how a simplified computational domain including a row of cilia can be replicated in the 𝑥- direction, without being influenced by the boundaries. As mentioned in the first section, nanoparticles smaller than 200 nm are able to penetrate into mucus [37] and reach to the PCL. Therefore, in this study, the tracers are tracked after entering the mucus layer. Herein, the tracers are initially (𝑡 = 𝑡0 ) set at an assumed distance of 𝑑 = 160 nm apart from each other within a rectangular patch of 16 × 16 tracers right above the ciliary tips. Note that as soon as the distance between the tracers and the epithelium becomes smaller than 𝑑∕2, the particles are assumed deposited on the bottom wall.
Fig. 4. Schematics of (a) ciliary beat patterns and (b) ciliary packings implemented in the present study.
where 𝐱̇𝐩 denotes the particle velocity, 𝐮 = (𝐱, 𝑡) is the fluid velocity, and 𝛥𝐮 is the total or material acceleration. As in [42], the inertial 𝛥𝑡 parameter 𝜖 can be described as 𝜖=
2 2 𝑎 𝑅𝑒(1 − 𝜌𝑟 ), 9
(2)
where 𝑎 is dimensionless particle diameter normalised by the characteristic length scale (i.e. cilium length 𝐿) of the flow and 𝜌𝑟 is the ratio 𝐿2 of particle density to that of fluid. Given the Reynolds number 𝑅𝑒 = 𝜈𝑇 (where 𝜈 is the fluid kinematic viscosity and 𝑇 is the ciliary beat period) 3 −2 and assuming a particle density of the order of ≈ 10 g∕cm [43] (for nanoparticles used as therapeutic carriers), one would obtain 𝜖 of the order of ≈ 10−8 which is 𝜖 ≪ 1. Hence, Eq. (1) can be reduced to that of a passive fluid particle motion; 𝑑𝐱𝐩 𝑑𝑡
= 𝐮(𝐱(𝑡), 𝑡),
3. Results and discussion After introducing all the simulation cases and their corresponding parameter set, first, the results of particle tracking are presented and the obtained findings are discussed from the clinical point of view. This is then followed by velocity field information related to each case study in order to better understand the physical behaviour of particle transports. We examine the transport of particles released in a cilia-driven periciliary liquid flow for six distinct cases, as listed in Table 1. Case 1 is referred to as ‘standard case’ which is characterised by (i) the 3D beat kinematics reported in [46]; (ii) ciliary length of 7 μm; (iii) beating frequency of 10 Hz [47]; (iv) and triangular ciliary pack [48]; see Fig. 4
(3)
thus, herein, we assume the particles advected are passive, infinitely small, and massless. This can be supported by referring to the work of Buchmann [44] in which it has been observed that for the particle 3
Computers in Biology and Medicine 117 (2020) 103595
S.M. Vanaki et al.
Fig. 5. Particle tracking after 250 beating cycles (∼ 𝑡𝑒𝑛𝑑 = 25 s) for the standard case (i.e. #1 in Table 1). A patch of 16 × 16 tracers (coloured in red and blue in 50:50 ratio) is released at 𝑡 = 0 s located at (0, 0, 8 μm).
Fig. 6. Fate of the released particle tracers after 250 beating cycles for the simulation cases 1–6 of Table 1. A patch similar to the one illustrated in Fig. 5 is released at 𝑡 = 0 s from the origin. Note, where deposition % is not shown it is 0%.
for more details on the beat patterns and the packings between cilia that will be considered in the present work. The model of Sanderson and Sleigh [48] interpreted an area of ciliary activity and reported triangular packing between cilia on the epithelial surface. Dense packing of the cilia, however, makes it difficult to visualise them by means of scanning electron microscopy. Therefore, we suggest another possible packing type; i.e. square pack (case 6) which has commonly been used in muco-ciliary studies (e.g. [26,28,32]). In this study, we refer to cases 1 and 6 as ‘healthy state’. Note that, for these two cases, the ciliary density (i.e. cilium/area taken up for each cilium) is kept constant. On the other hand, cases 2–5 are referred to as ‘diseased state’, which is
associated with acquired and inherited respiratory disorders; see the footnote of Table 1 for more information on common disorders and their associated abnormalities. Taking inspiration from an area of ciliary activity on cultured rabbit tracheal epithelium [48], the initial arrangements of the cilia are chosen such that there is a time lag of 𝑇 ∕12 between successive cilia. The selection of the rabbit lung for the ciliary arrangements is due to its similarities to humans in terms of anatomy, physiology, genetic, and biochemical similarities [58]. Here, a ciliated domain including 96 cilia is assumed. The direction in which the metachronal wave travels is chosen to be diagonal, at an angle of 135◦ to the effective stroke. 4
Computers in Biology and Medicine 117 (2020) 103595
S.M. Vanaki et al.
Fig. 7. Averaged linear distance vs path distance travelled by a patch of particles during 250 cycles. Note that all the values in cases 2–6 have been normalised by the linear distance obtained from the standard case. The error bar is an indication of particle dispersion intensity in the flow field.
For the standard case, initially, a patch of 16×16 tracers (coloured in red and blue in 50:50 ratio) is introduced above the tips of the cilia, and the particles are tracked for 250 beating cycles, as shown in Fig. 5. With a slight deviation towards the 𝑦-direction, the seeded particles have a significant forward transport in the same direction as the effective cilia motion, i.e. the 𝑥-direction. Moreover, the majority of particles have remained above the cilia tip (close to the upper boundary), and only 0.78% of them were deposited after the prescribed period of time. Such particle transport behaviour may be viewed from two different aspects: (i) toxic agents that can pass the first defencive barrier in the healthy lung are more likely to be transported back to the mucus gel and be cleared via the muco-ciliary mechanism. This shows the key role
of ciliary actions in protecting our respiratory system against harmful submicron substances — thus hindering the rate of their deposition onto the epithelium; (ii) the active pharmaceutical ingredients may not easily reach to the epithelium and eventually the blood stream. This finding is of importance given that the pulmonary tissues are a noninvasive absorption route for the air-to-blood drug delivery – not only for local treatment of lung diseases, but also to treat other diseases such as central nervous system (CNS) disorders [59]. The colour-coded particle tracking is also carried out for other simulation cases and compared against the standard case, as shown in Fig. 6. The rate of particle deposition is provided for each case. The lower beating frequency (case 2), which is a common defect 5
Computers in Biology and Medicine 117 (2020) 103595
S.M. Vanaki et al. Table 1 Simulation cases in the present work. Simulation cases Model parameters
#1 (standard) #2 #3 #4 #5 #6
State
Frequency Length Beat pattern
Packing
10 Hz 5 Hza 10 Hz 10 Hz 10 Hz 10 Hz
Triangular ✓ Triangular Triangular Triangular Triangular Square ✓
7 7 4 7 7 7
μm μm 𝛍𝐦b μm μm μm
3D 3D 3D Stiff planarc Wiperd 3D
Healthy Diseased ✓ ✓ ✓ ✓
a
Reduced ciliary beat frequency – e.g. primary ciliary dyskinesia (PCD) [49]; asthma [50]; cigarette smoking [51]; environmental pollutants [52]. b Shorter cilia – e.g. cigarette smoking [53,54]; chronic obstructive pulmonary disease (COPD) [55]. c & d Abnormal ciliary motions – e.g. PCD [56,57].
in PCD [49] and asthma patients [50], cigarette smokers [51], and individuals exposed to severe air pollution [52], may decrease the rate of particle deposition. Additionally, similar to the standard case, most of the particles have a tendency to stay above the cilia tips. Such behaviour is much more significant in case 3, where cilia are short in length. This can be expected because the location where the tracers are released is not directly perturbed by ciliary actions. Therefore, the chance of delivering drugs to the epithelium of COPD patients and smokers, who have shorter cilia (see Table 1), may be notably reduced. Primary ciliary dysfunction (PCD) can cause abnormal ciliary beat patterns in the lungs — such as stiff planar and wiper cilia motions [56,57]. Prescribing these two abnormal cilia beats in the model significantly hindered the transport of particles. For the stiff planar beat (case 4), a significant number of particles can surprisingly move in the negative 𝑥-direction from the release point, with 1.17% of them being deposited. For the stiff planar beat (case 5), the particle tracers do not experience a significant transport in the 𝑥-direction, and move close to one another. From a pharmaceutical point of view, this would prolong the presence of drug particles in the desired site of action. When the cilia were arranged in a square pack in case 6, a deposition rate of 4.69% was observed and the particles experienced a higher mixture compared to the standard case. By looking at healthy state results, it is evident that some measures (such as modification of the particle surface [60]) must be taken to avoid immediate clearance of the particles and to prolong their residence time at the site of absorption. To quantify these statements, we compute both linear and path distances that a patch of particles migrates during 250 cycles, and present results for all the simulation cases in Fig. 7. The linear distance is obtained by the projection of particle location on the 𝑥-axis at every cycle of simulation, as this is the primary direction that particles are cleared from the respiratory tracts or may move to the lower airways. Using particle trajectories, the path distance can be calculated as √ ( ) ( ) ( ) 𝑑𝑥𝑝 2 𝑑𝑦𝑝 2 𝑑𝑧𝑝 2 𝑃 𝑎𝑡ℎ 𝑑𝑖𝑠𝑡𝑎𝑛𝑐𝑒 = + + 𝑑𝑡 (4) ∫ 𝑑𝑡 𝑑𝑡 𝑑𝑡
Fig. 8. Sensitivity tests. The effect of (a) patch size of particles, and (b) different releasing time of particles. 𝑇 denotes the period of beating cycle, and the results have been provided for 100 cycles for the standard case.
velocity profiles are constructed by plotting the time-averaged velocity over the stable cycle at each height 𝑧 of the cube located at the centre of the ciliated domain (see Fig. 9a), which can be written as ⟨𝐮⟩ = (1∕𝑇 𝐴𝐿 )
∫
𝐮𝑑𝐴𝐿 𝑑𝑡,
(5)
where the cube is of cross section 𝐴𝐿 = 8 × 8 μm2 and of height 𝑧 = 8 μm. Notice that because the cilia have abnormal beat patterns in diseased states, variations in the flow field should be analysed threedimensionally. Therefore, three representative planes are chosen to visualise the velocity contours — which are the blue, red, and grey sides of the cube in Fig. 9a. Overall, the 𝑥-component of velocity ⟨𝑢⟩ = (1∕𝑇 𝐴) ∫ 𝑢𝑑𝐴𝑑𝑡 has the highest value in magnitude for all the simulation cases, except for stiff planar and wiper ciliary patterns. For these two, the 𝑦-component of velocity ⟨𝑣⟩ is dominant, with only small changes in ⟨𝑢⟩. This signifies the severity of diseased lungs under such abnormal beating conditions, as the net axial transport of deposited foreign particles could be greatly impaired. However, there would be a high chance for therapeutic particles to be absorbed by the epithelium, due to the high residence time of particles. While the calculated values of ⟨𝑤⟩ (i.e. the 𝑧-component of velocity) were positive for most of the cases, the square pack led to negative values along the 𝑧-axis. This may be surprising, and simply means that the dispersed particles would be transported to the epithelial surface upon deposition in case 6, whereas the particles could experience upwards transport in the other cases. To summarise the findings discussed in this section, we have looked at the variations of the vector fields over one beating cycle for each simulation case, and a schematic diagram of fluid flow has been interpreted and shown on 𝑥-𝑧 plane in Fig. 10. This could qualitatively be used to
It was seen that the linear distance halved with decreasing the ciliary beating frequency from 10 to 5 Hz. Reducing the length of cilia from 7 to 4 μm led to approximately 70% reduction in the linear transport of particles. The results for cases 4 and 5 reveal that the particles have experienced a relatively large displacement, which must be in the 𝑦-direction, as the linear distance has been quite small. Finally, the square pack resulted in a 40% reduction in linear distance as compared with the standard case. Further to these results, two sensitivity tests were performed, and the results showed that the effect of patch size and releasing time of particles on their transport is insignificant, as seen in Fig. 8. The role of abnormal beats in particle transport can be better understood by having the velocity field information. As such, the 𝑥-, 𝑦-, and 𝑧-components of the velocity profile and velocity contour, for each simulation case, are shown in Figs. 9b and 9c, respectively. The 6
Computers in Biology and Medicine 117 (2020) 103595
S.M. Vanaki et al.
Fig. 9. (a) A square column with the cross sectional area of 𝐴𝐿 ∼ 8 × 8 μm2 is shown in the middle of the computational domain through which (b) the time-averaged velocity profiles in the 𝑥-, 𝑦-, and 𝑧-directions are computed. (c) Time-averaged velocity contours of ⟨𝑢⟩, ⟨𝑣⟩, and ⟨𝑤⟩ are shown. All the plots are presented for the simulation cases 1–6 (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.).
7
Computers in Biology and Medicine 117 (2020) 103595
S.M. Vanaki et al.
an advection dominant PCL flow for the healthy case, where the typical ciliary beating frequency is around 10 Hz, the ciliary waves are antiplectic, and the diffusivity of nanoparticles 𝐷 in a watery PCL flow is of order of ≈ 10−12 . Similarly, Ding et al. [32] and Chateau et al. [27] reported that fluid advection is the dominant mechanism for the particle transport inside the PCL. This suggests the necessity of modifying the physico-chemical properties of the nanocarriers to increase their absorption rate. We conclude by noting that the present model is a basic one and effects such as the particle Brownian motion, particle–cilia as well as particle–particle interactions, and particle charges may be implemented in future models. Declaration of competing interest The authors confirm that there are no known conflicts of interest associated with this publication and there has been no significant financial support for this work that could have influenced its outcome. Acknowledgements Shayan M. Vanaki gratefully acknowledges Queensland University of Technology (QUT) for support through the QUTPRA scholarship. The authors appreciatively acknowledge the High Performance Computing Unit of QUT. References [1] S. Khalafvand, E. Ng, L. Zhong, CFD simulation of flow through heart: a perspective review, Comput. Methods Biomech. Biomed. Eng. 14 (01) (2011) 113–132. [2] D.M. Sforza, C.M. Putman, J.R. Cebral, Computational fluid dynamics in brain aneurysms, Int. J. Numer. Methods Biomed. Eng. 28 (6–7) (2012) 801–808. [3] M. Ferrua, R. Singh, Modeling the fluid dynamics in a human stomach to gain insight of food digestion, J. Food Sci. 75 (7) (2010) R151–R162. [4] W. Hofmann, Modelling inhaled particle deposition in the human lung—A review, J. Aerosol Sci. 42 (10) (2011) 693–724. [5] P.W. Longest, K. Bass, R. Dutta, V. Rani, M.L. Thomas, A. El-Achwah, M. Hindle, Use of computational fluid dynamics deposition modeling in respiratory drug delivery, Expert Opinion Drug Deliv. 16 (1) (2019) 7–26. [6] M.R. Knowles, R.C. Boucher, Mucus clearance as a primary innate defense mechanism for mammalian airways, J. Clinic. Investigation 109 (5) (2002) 571–577. [7] D.J. Thornton, J.K. Sheehan, From mucins to mucus: toward a more coherent understanding of this essential barrier, Proc. Amer. Thoracic Soc. 1 (1) (2004) 54–61. [8] M. Geiser, L.M. Cruz-Orive, V.I. Hof, P. Gehr, Assessment of particle retention and clearance in the intrapulmonary conducting airways of hamster lungs with the fractionator, J. Microsc. 160 (1) (1990) 75–88. [9] M. Geiser, W.G. Kreyling, Deposition and biokinetics of inhaled nanoparticles, Particle Fibre Toxicol. 7 (1) (2010) 2. [10] C. Yu, J. Hu, B. Yen, D. Spektor, M. Lippmann, Models for mucociliary particle clearance in lung airways, in: Aerosols: Research, Risk Assessment and Control Strategies, Chelsea, MI, Lewis, 1986, pp. 569–578. [11] P.S. Lee, T.R. Gerrity, F.J. Hass, R.V. Lourenco, A model for tracheobronchial clearance of inhaled particles in man and a comparison with data, IEEE Trans. Biomed. Eng. (11) (1979) 624–630. [12] R. Cuddihy, H. Yeh, Respiratory tract clearance of particles and substances dissociated from particles, in: Inhalat. Toxicol., Springer, 1988, pp. 169–193. [13] W. Hofmann, L. Koblinger, Stochastic model of radon daughter deposition and clearance in human bronchial airways, in: IRPA9: 1996 International Congress on Radiation Protection. Proceedings, vol. 2, 1996. [14] B. Asgharian, W. Hofmann, F. Miller, Mucociliary clearance of insoluble particles from the tracheobronchial airways of the human lung, J. Aerosol Sci. 32 (6) (2001) 817–832. [15] R. Sturm, W. Hofmann, Stochastic modeling predictions for the clearance of insoluble particles from the tracheobronchial tree of the human lung, Bull. Math. Biol. 69 (1) (2007) 395. [16] J. Blake, H. Winet, On the mechanics of muco-ciliary transport, Biorheology 17 (1–2) (1980) 125–134. [17] D. Smith, E. Gaffney, J. Blake, Erratum: A viscoelastic traction layer model of muco-ciliary transport, Bull. Math. Biol. 69 (4) (2007) 1449. [18] P. Kurbatova, N. Bessonov, V. Volpert, H. Tiddens, C. Cornu, P. Nony, D. Caudri, C.W. Group, et al., Model of mucociliary clearance in cystic fibrosis lungs, J. Theoret. Biol. 372 (2015) 81–88.
Fig. 10. Schematic diagram of fluid flow (𝑥-𝑧 plane view) in the simulation cases 1–6. The arrows represent the flow direction, indicating in which directions particles could transport when entering the cilia-induced flow field. The chosen plane is approximately located in the middle of the ciliated domain.
predict the possible directions in which particles might transport when entering the cilia-induced flow fields. 4. Conclusions The major objective of the present work was to better understand the fate of biodegradable nanoparticles, used frequently as drug delivery vehicles, in the presence of the ciliated cells where their transport can be significantly altered by ciliary beats. To achieve this, we set up a computational model based on the projection method combined with the IB and particle track methods to simulate the transport of inhaled particles deposited onto the periciliary liquid layer in healthy and diseased lungs. A detailed study of the particle transport induced by abnormal ciliary beats was performed, and the results showed that impairments such as low ciliary beat frequency and short ciliary length (seen in PCD patients and cigarette smokers) can significantly reduce the rate of transport of therapeutic particles in the proximal direction of the tracheobronchial tree. The configurations corresponding to stiff planar and wiper ciliary motions (seen on PCD patients) were found to greatly reduce the transport of such particles. These findings show that the site-specific delivery of drugs may be achieved in diseased states. However, the beating of cilia creates a barrier for the particles to easily reach to the epithelial cells. Our quantitative results also reveal that the chance of drug retention and absorption by pulmonary tissues may be markedly diminished in healthy lungs if drugs are intended to be delivered into the blood stream by aerosol carriers, specifically in the upper airways. However, such observation might be changed by the diffusive effect of particles in a cilia-induced flow field which has not been implemented in the present model. Nonetheless, based on the flow rate changes during one cycle (shown in [39]) for the present model, it is likely to have 8
Computers in Biology and Medicine 117 (2020) 103595
S.M. Vanaki et al.
[42] G. Haller, T. Sapsis, Where do inertial particles go in fluid flows? Physica D 237 (5) (2008) 573–583. [43] W.S. Cheow, M.L.L. Ng, K. Kho, K. Hadinoto, Spray-freeze-drying production of thermally sensitive polymeric nanoparticle aggregates for inhaled drug delivery: effect of freeze-drying adjuvants, Int. J. Pharm. 404 (1–2) (2011) 289–300. [44] A.L. Buchmann, L.J. Fauci, K. Leiderman, E.M. Strawbridge, L. Zhao, Flow induced by bacterial carpets and transport of microscale loads, in: Applications of Dynamical Systems in Biology and Medicine, Springer, 2015, pp. 35–53. [45] E. Herawati, D. Taniguchi, H. Kanoh, K. Tateishi, S. Ishihara, S. Tsukita, Multiciliated cell basal bodies align in stereotypical patterns coordinated by the apical cytoskeleton, J. Cell Biol. 214 (5) (2016) 571–586. [46] L. Gheber, Z. Priel, Extraction of cilium beat parameters by the combined application of photoelectric measurements and computer simulation, Biophys. J. 72 (1) (1997) 449–462. [47] D.J. Lubkin, E.A. Gaffney, J.R. Blake, A viscoelastic traction layer model of muco-Ciliary transport, Bull. Math. Biol. (ISSN: 1522-9602) 69 (1) (2006) 289. [48] M. Sanderson, M. Sleigh, Ciliary activity of cultured rabbit tracheal epithelium: beat pattern and metachrony, J. Cell Sci. 47 (1) (1981) 331–347. [49] C.M. Rossman, R.M. Lee, J.B. Forrest, M.T. Newhouse, Nasal ciliary ultrastructure and function in patients with primary ciliary dyskinesia compared with that in normal subjects and in subjects with various respiratory diseases, Amer. Rev. Respiratory Disease 129 (1) (1984) 161–167. [50] B. Thomas, A. Rutman, R.A. Hirst, P. Haldar, A.J. Wardlaw, J. Bankart, C.E. Brightling, C. O’callaghan, Ciliary dysfunction and ultrastructural abnormalities are features of severe asthma, J. Allergy Clin. Immunol. 126 (4) (2010) 722–729. [51] A.E. Tilley, M.S. Walters, R. Shaykhiev, R.G. Crystal, Cilia dysfunction in lung disease, Annu. Rev. Physiol. 77 (2015) 379–406. [52] H. Riechelmann, K. Kienast, J. Schellenberg, W. Mann, An in vitro model to study effects of airborne pollutants on human ciliary activity, Rhinology 32 (3) (1994) 105–108. [53] P.L. Leopold, M.J. O’Mahony, X.J. Lian, A.E. Tilley, B.-G. Harvey, R.G. Crystal, Smoking is associated with shortened airway cilia, PLoS One 4 (12) (2009) e8157. [54] J. Hessel, J. Heldrich, J. Fuller, M.R. Staudt, S. Radisch, C. Hollmann, B.-G. Harvey, R.J. Kaner, J. Salit, J. Yee-Levin, et al., Intraflagellar transport gene expression associated with short cilia in smoking and COPD, PLoS One 9 (1) (2014) e85453. [55] H.C. Lam, S.M. Cloonan, A.R. Bhashyam, J.A. Haspel, A. Singh, J.F. Sathirapongsasuti, M. Cervo, H. Yao, A.L. Chung, K. Mizumura, et al., Histone deacetylase 6–mediated selective autophagy regulates COPD-associated cilia dysfunction, J. Clinic. Investigation 123 (12) (2013) 5212–5230. [56] C.M. Rossman, J.B. Forrest, R.M. Lee, M.T. Newhouse, The dyskinetic cilia syndrome: ciliary motility in immotile cilia syndrome, Chest 78 (4) (1980) 580–582. [57] M.A. Chilvers, A. Rutman, C. O’callaghan, Ciliary beat pattern is associated with specific ultrastructural defects in primary ciliary dyskinesia, J. Allergy Clin. Immunol. 112 (3) (2003) 518–524. [58] N.A. Kamaruzaman, E. Kardia, N. Kamaldin, A.Z. Latahir, B.H. Yahaya, The rabbit as a model for studying lung disease and stem cell therapy, BioMed Res. Int. 2013 (2013). [59] P. Noymer, S. Biondi, D. Myers, J. Cassella, Pulmonary delivery of therapeutic compounds for treating CNS disorders, Therapeutic Deliv. 2 (9) (2011) 1125–1140. [60] J.C. Sung, B.L. Pulliam, D.A. Edwards, Nanoparticles for drug delivery to the lungs, Trends Biotechnol. 25 (12) (2007) 563–570.
[19] G.R. Fulford, J.R. Blake, Muco-ciliary transport in the lung, J. Theoret. Biol. 121 (4) (1986) 381–402. [20] W. Lee, P. Jayathilake, Z. Tan, D. Le, H. Lee, B. Khoo, Muco-ciliary transport: effect of mucus viscosity, cilia beat frequency and cilia density, Comput. & Fluids 49 (1) (2011) 214–221. [21] R. Chatelin, P. Poncet, A hybrid grid-particle method for moving bodies in 3D stokes flow with variable viscosity, SIAM J. Sci. Comput. 35 (4) (2013) B925–B949. [22] P. Jayathilake, D. Le, Z. Tan, H.P. Lee, B.C. Khoo, A numerical study of mucociliary transport under the condition of diseased cilia, Comput. Methods Biomech. Biomed. Eng. 18 (9) (2015) 944–951. [23] M. Sedaghat, M. Shahmardan, M. Norouzi, P. Jayathilake, M. Nazari, Numerical simulation of muco-ciliary clearance: immersed boundary-lattice Boltzmann method, Comput. & Fluids 131 (2016) 91–101. [24] R. Chatelin, P. Poncet, A parametric study of mucociliary transport by numerical simulations of 3D non-homogeneous mucus, J. Biomech. 49 (9) (2016) 1772–1780. [25] H. Guo, E. Kanso, A computational study of mucociliary transport in healthy and diseased environments, Eur. J. Comput. Mech. 26 (1–2) (2017) 4–30. [26] R. Chatelin, D. Anne-Archard, M. Murris-Espin, M. Thiriet, P. Poncet, Numerical and experimental investigation of mucociliary clearance breakdown in cystic fibrosis, J. Biomech. 53 (2017) 56–63. [27] S. Chateau, U. D’Ortona, S. Poncet, J. Favier, Transport and mixing induced by beating cilia in human airways, Front. Physiol. 9 (2018) 161. [28] Y.L.R. Quek, K.M. Lim, K.-H. Chiam, Three-dimensional computational model of multiphase flow driven by a bed of active cilia, Comput. & Fluids 170 (2018) 222–235. [29] R.H. Dillon, L.J. Fauci, C. Omoto, X. Yang, Fluid dynamic models of flagellar and ciliary beating, Ann. New York Acad. Sci. 1101 (1) (2007) 494–505. [30] S.M. Mitran, Metachronal wave formation in a model of pulmonary cilia, Comput. Struct. 85 (11–14) (2007) 763–774. [31] S. Lukens, X. Yang, L. Fauci, Using Lagrangian coherent structures to analyze fluid mixing by cilia, Chaos 20 (1) (2010) 017511. [32] Y. Ding, J.C. Nawroth, M.J. McFall-Ngai, E. Kanso, Mixing and transport by ciliary carpets: a numerical study, J. Fluid Mech. 743 (2014) 124–140. [33] R. Ghosh, G.A. Buxton, O.B. Usta, A.C. Balazs, A. Alexeev, Designing oscillating cilia that capture or release microscopic particles, Langmuir 26 (4) (2009) 2963–2968. [34] A. Bhattacharya, G.A. Buxton, O.B. Usta, A.C. Balazs, Propulsion and trapping of microparticles by active cilia arrays, Langmuir 28 (6) (2012) 3217–3226. [35] Y. Ding, E. Kanso, Selective particle capture by asynchronously beating cilia, Phys. Fluids 27 (12) (2015) 121902. [36] C.A. Ruge, J. Kirch, C.-M. Lehr, Pulmonary drug delivery: from generating aerosols to overcoming biological barriers—therapeutic possibilities and technological challenges, Lancet Respiratory Med. 1 (5) (2013) 402–413. [37] X. Murgia, P. Pawelzyk, U.F. Schaefer, C. Wagner, N. Willenbacher, C.-M. Lehr, Size-limited penetration of nanoparticles into porcine respiratory mucus after aerosol deposition, Biomacromolecules 17 (4) (2016) 1536–1542. [38] S. Tsorng, H. Capart, J. Lai, D. Young, Three-dimensional tracking of the long time trajectories of suspended particles in a lid-driven cavity flow, Exp. Fluids 40 (2) (2006) 314–328. [39] S.M. Vanaki, D. Holmes, P.G. Jayathilake, R. Brown, Three-dimensional numerical analysis of periciliary liquid layer: ciliary abnormalities in respiratory diseases, Appl. Sci. 9 (19) (2019) 4033. [40] M. Maxey, The gravitational settling of aerosol particles in homogeneous turbulence and random flow fields, J. Fluid Mech. 174 (1987) 441–465. [41] K.-K. Tio, A. Liñán, J.C. Lasheras, A.M. Gañán-Calvo, On the dynamics of buoyant and heavy particles in a periodic stuart vortex flow, J. Fluid Mech. 254 (1993) 671–699.
9
Contents lists available at ScienceDirect
Computers in Biology and Medicine journal homepage: www.elsevier.com/locate/compbiomed
Transport and fate of inhaled particles after deposition onto the airway surface liquid: A 3D numerical study Shayan M. Vanaki a ,∗, David Holmes a ,∗, Kabir Suara a , Pahala Gedara Jayathilake b , Richard Brown a a b
School of Chemistry, Physics and Mechanical Engineering, Queensland University of Technology, Brisbane, Queensland, Australia Department of Oncology, University of Oxford, UK
ARTICLE
INFO
Keywords: Therapeutic drug particles Particle tracking Particle deposition Cilia Periciliary liquid layer Respiratory diseases
ABSTRACT In this paper, a numerical investigation is carried out to provide insights into the fate of inhaled aerosols after their deposition on the lung lining fluid in both healthy and diseased states. Pulmonary drug delivery is a well-known non-invasive route of administration compared to intravenous delivery. Aerosol particles are formulated and used as drug carriers, which are then sent to the airways using aerosol drug delivery devices. This approach is useful for site-specific treatment of lung diseases, treatment of central nervous system (CNS) disorders and a variety of other diseases. Bioavailability of the inhaled therapeutic particles after landing on the airway lining fluid can be significantly altered by the lung muco-ciliary clearance, a process through which hairlike structures known as cilia beat in a harmonised manner and induce the mucus in the proximal direction, leading to an effective clearance of the foreign inhaled particles entrapped by this sticky layer from the airways. Here, we set up a 3D computational model of ciliary arrays interacting with periciliary liquid film (i.e. confined between the epithelium and mucus layer) and a detailed analysis is conducted to better understand the fate of drug nanoparticles that are able to penetrate the mucus. Consistent with clinical findings, we find that the actions of cilia result in a low rate of drug retention and absorption by the pulmonary tissues in healthy lungs. However, under conditions associated with abnormal ciliary beats, the retention time of particles is notably increased at the site of release. Nonetheless, the results associated with some of the ciliary impairments reveal that deposition of drug aerosols on the ciliated cells may still be a significant challenge. These findings have potentially important implications on the modification of therapeutic drug particles to achieve a higher absorption rate.
1. Introduction Computational Fluid Dynamics (CFD) is garnering a great deal of attention in the medical context, as it can be used as a powerful tool to provide valuable insights into physical phenomena occurring inside the human body without the need of expensive clinical experiments and the associated ethical challenges. As such, many researchers have used CFD methods to study the transport of biofluids through various human organs, such as the heart [1], brain [2], and stomach [3]. Lung modelling has also been of interest for many researchers over the years. Modelling inhaled particle deposition in the lung airways has already received a lot of attention from scientists [4]; see e.g. Hofmann [4] and Longest et al. [5] for reviews. However, the clearance process of aerosols from the pulmonary airways has been more challenging to investigate numerically due to the multi-scale, multi-physics, and multi-components involved. In particular, the conductive airways
are protected by a lining fluid called ‘airway surface liquid’ (ASL) which is a multiphase fluid composed of two layers [6]: (i) an aqueous periciliary liquid layer (PCL) that resides on the airway epithelium and contains a carpet of motile microstructures called cilia, and (ii) a nonNewtonian, viscoelastic mucus layer [7] which captures the foreign particles (e.g. dust, PM, bacteria, etc.). The harmonised beating of cilia within the PCL propels the mucus layer continuously towards the trachea and eventually out of the lungs. This fluid–structure interaction (FSI) process is known as ‘muco-ciliary clearance’ which is the primary defence mechanism, with 80%–90% of inhaled particles being discharged from the upper and central lung within 24 h [8,9]. Some works in the literature have simplified such a complex clearance process by using theoretical assumptions, where the transport of particles within the mucus layer is merely estimated at every airway generation, without taking the interaction of cilia with the PCL or mucus layer
∗ Corresponding authors. E-mail addresses: [email protected], [email protected] (S.M. Vanaki), [email protected] (D. Holmes).
https://doi.org/10.1016/j.compbiomed.2019.103595 Received 2 October 2019; Received in revised form 14 December 2019; Accepted 27 December 2019 Available online 7 January 2020 0010-4825/© 2020 Elsevier Ltd. All rights reserved.
Computers in Biology and Medicine 117 (2020) 103595
S.M. Vanaki et al.
Fig. 1. Fate of an aerosol drug after deposition in the lungs. (1) Release of the active nano-sized pharmaceutical ingredient from the micro-sized aerosol particle deposited on the mucus gel. (2) Absorption of the active ingredient across the pulmonary epithelium into the bloodstream. (3) Clearance of the undissolved particle. Source: Figure adapted from [36].
into account; see e.g. [10–15]. Numerical modelling approaches that explicitly account for cilia PCL interaction, including continuum cilia models [16–18] as well as discrete cilia models (comprised of prescribed ciliary beating models [19–28] and FSI models [29–31]) have been developed, potentially allowing for more detailed investigation into the complex physics involved in particle transport. None of the above-mentioned works have investigated the transport of particles across the lung lining layer, which is of importance, particularly in the context of aerosol medicine and pulmonary drug delivery, and the focus has been limited to fluid mixing (e.g. [27,32]) and particle capture by beating cilia (e.g. [33–35]). Dissolution of active pharmaceutical ingredients (or carrier particles) in the ASL and their absorption into the bloodstream would be desirable for pulmonary specialists. However, as shown in Fig. 1, the persistence and subsequent absorption of the inhaled drug carrier in the lung is significantly influenced and reduced by the muco-ciliary clearance process, thus substantially minimising the impact of therapeutic drug particles [36]. Nanoparticles intended for use as drug carriers can permeate through the pores of the mucus network [37] and reach to the underlying PCL, leading to absorption across epithelial and endothelial cells. Very little is known about the fate and transport of these nanodrugs after deposition and penetration inside the mucus layer [36], especially in the presence of an active biological factor such as cilia. The goal of this paper is, therefore, to understand the fate of drug carriers in a flow generated by a periodic array of beating cilia in the lung. Also, investigated is the transport of such particles in diseased lungs, due to abnormal ciliary beatings, compared to that in healthy states. To achieve these objectives, in Section 2, we use a computational model based on the projection method for solving a three-dimensional Navier–Stokes equation and a particle tracking method. In Section 3, we track and quantify the transport of inhaled ultrafine particles in the periciliary liquid film under both healthy and diseased lung conditions. Finally, Section 4 concludes this paper and future numerical challenges and directions are provided.
Fig. 2. Particle tracking verification over a span of 43 s. (a) Present model. (b) Experimental work of Tsorng et al. [38].
the motile cilia and fluid is computed explicitly through the immersed boundary (IB) method, analogous to the one-way FSI model presented in our previous work [39] and hence not repeated here. By feeding the obtained velocity fields from the FSI model into the Lagrangian particle tracking model, a detailed numerical study is then carried out to better understand the fate of drug nanoparticles that are able to penetrate the mucus. 2.1. Particle tracking The motion of finite-inertial spherical particle in fluid can be described by the full four dimensional Maxey–Riley (MR) ordinary differential equation for particle position [40,41]. The inertial particle positions as well as initial particle velocities which are generally unknown are required to integrate the MR equation. However, using the slow manifold approximation, Ref. [42] showed that asymptotic motion of finite-sized particles can be described by a two-dimensional system in the position, that is
2. Numerical methods In the present work, the PCL flow is modelled and solved using a finite-difference projection algorithm, and the interaction between
𝐱̇𝐩 = 𝐮 + 𝜖 2
𝛥𝐮 , 𝛥𝑡
(1)
Computers in Biology and Medicine 117 (2020) 103595
S.M. Vanaki et al.
sizes comparable to cilia movement, the difference between the motion of the fully resolved particles and that of passive tracers is small. To verify the appropriateness of the particle tracking method, we track a particle trajectory in a 3D lid-driven cavity numerically and compare the result with that observed in the experimental work of Tsorng et al. [38]. Since fluid advection is the dominant mechanism for the transport of particles within the cilia-induced PCL flow [27,32], it is still reasonable to have this comparison with an advection dominant cavity problem where Re = 470. The simulated trajectory was obtained by integrating Eq. (3) from the initial location of the observation over the time span of the observation. The comparison is presented in Fig. 2. The simulated and observed trajectories are similar for a major part of the time span with some differences occurring at the end of the tracking. Eq. (3) is solved using the Runge–Kutta (4,5) method, i.e. the MATLAB ode45, to obtain particle trajectories. 2.2. Numerical setting and boundary conditions Fig. 3. Schematic of a row of cilia to illustrate how a ciliated domain can be replicated. For simplicity, it is assumed that the shape of each cilium changes from phase 1 to 4 through one beat cycle. Extra cilia are added to each side of the ciliated domain such that the cilium with the phase number of 4 has the one with phase number of 1 on its right side at the beginning of each beat cycle. This would allow for replicating the ciliated domain and removing the discrete boundary effects. Note that adding extra cilia to any side of the ciliated domain would require extension of the periodic boundary conditions, from 𝐵𝐶𝑖𝑛𝑖𝑡𝑖𝑎𝑙 to 𝐵𝐶𝑓 𝑖𝑛𝑎𝑙 .
For all of the simulation cases presented in this paper, the fluid grid size is ℎ = 0.5 μm, the time step of the computation is 𝛥𝑡 = 1 × 10−6 s, the fluid dynamic viscosity is 𝜇 = 0.01 g∕cm.s, the fluid density is 𝜌 = 1 g∕cm3 , and the number of control points along the cilia is set at 20. As for the boundary conditions, periodic boundaries are applied in the wall-parallel directions. The bottom boundary at the epithelium is treated as a no-slip wall in our simulations consistent with the velocity field at the epithelium boundaries of cultured mouse tracheal cells [45]. Given that the boundary condition at the PCL-mucus interface is influenced by the surface tension force between the mucus layer and PCL as well as the complex fluid rheology of the mucus, it is reasonable to assume the top boundary as a free slip wall in the present mucus-free model. To allow the tracking of the dispersed particles over a longer period of time, we perform the following steps: (i) once the simulations reach the steady state condition, the stable (or the last) cycle is replicated 𝑁𝑐𝑦𝑐𝑙𝑒 times until it reaches at 𝑡𝑒𝑛𝑑 = 𝑁𝑐𝑦𝑐𝑙𝑒 ×𝑇 ; (ii) the physical domain is extended in the 𝑥- and 𝑦- directions through boundary replication. To achieve this and remove any boundary effects, extra rows/columns of cilia with specific phase numbers are added to the each side of the square ciliated domain such that the selected periodic domain for particle tracking does not experience disturbance when it is replicated. Fig. 3, for example, illustrates how a simplified computational domain including a row of cilia can be replicated in the 𝑥- direction, without being influenced by the boundaries. As mentioned in the first section, nanoparticles smaller than 200 nm are able to penetrate into mucus [37] and reach to the PCL. Therefore, in this study, the tracers are tracked after entering the mucus layer. Herein, the tracers are initially (𝑡 = 𝑡0 ) set at an assumed distance of 𝑑 = 160 nm apart from each other within a rectangular patch of 16 × 16 tracers right above the ciliary tips. Note that as soon as the distance between the tracers and the epithelium becomes smaller than 𝑑∕2, the particles are assumed deposited on the bottom wall.
Fig. 4. Schematics of (a) ciliary beat patterns and (b) ciliary packings implemented in the present study.
where 𝐱̇𝐩 denotes the particle velocity, 𝐮 = (𝐱, 𝑡) is the fluid velocity, and 𝛥𝐮 is the total or material acceleration. As in [42], the inertial 𝛥𝑡 parameter 𝜖 can be described as 𝜖=
2 2 𝑎 𝑅𝑒(1 − 𝜌𝑟 ), 9
(2)
where 𝑎 is dimensionless particle diameter normalised by the characteristic length scale (i.e. cilium length 𝐿) of the flow and 𝜌𝑟 is the ratio 𝐿2 of particle density to that of fluid. Given the Reynolds number 𝑅𝑒 = 𝜈𝑇 (where 𝜈 is the fluid kinematic viscosity and 𝑇 is the ciliary beat period) 3 −2 and assuming a particle density of the order of ≈ 10 g∕cm [43] (for nanoparticles used as therapeutic carriers), one would obtain 𝜖 of the order of ≈ 10−8 which is 𝜖 ≪ 1. Hence, Eq. (1) can be reduced to that of a passive fluid particle motion; 𝑑𝐱𝐩 𝑑𝑡
= 𝐮(𝐱(𝑡), 𝑡),
3. Results and discussion After introducing all the simulation cases and their corresponding parameter set, first, the results of particle tracking are presented and the obtained findings are discussed from the clinical point of view. This is then followed by velocity field information related to each case study in order to better understand the physical behaviour of particle transports. We examine the transport of particles released in a cilia-driven periciliary liquid flow for six distinct cases, as listed in Table 1. Case 1 is referred to as ‘standard case’ which is characterised by (i) the 3D beat kinematics reported in [46]; (ii) ciliary length of 7 μm; (iii) beating frequency of 10 Hz [47]; (iv) and triangular ciliary pack [48]; see Fig. 4
(3)
thus, herein, we assume the particles advected are passive, infinitely small, and massless. This can be supported by referring to the work of Buchmann [44] in which it has been observed that for the particle 3
Computers in Biology and Medicine 117 (2020) 103595
S.M. Vanaki et al.
Fig. 5. Particle tracking after 250 beating cycles (∼ 𝑡𝑒𝑛𝑑 = 25 s) for the standard case (i.e. #1 in Table 1). A patch of 16 × 16 tracers (coloured in red and blue in 50:50 ratio) is released at 𝑡 = 0 s located at (0, 0, 8 μm).
Fig. 6. Fate of the released particle tracers after 250 beating cycles for the simulation cases 1–6 of Table 1. A patch similar to the one illustrated in Fig. 5 is released at 𝑡 = 0 s from the origin. Note, where deposition % is not shown it is 0%.
for more details on the beat patterns and the packings between cilia that will be considered in the present work. The model of Sanderson and Sleigh [48] interpreted an area of ciliary activity and reported triangular packing between cilia on the epithelial surface. Dense packing of the cilia, however, makes it difficult to visualise them by means of scanning electron microscopy. Therefore, we suggest another possible packing type; i.e. square pack (case 6) which has commonly been used in muco-ciliary studies (e.g. [26,28,32]). In this study, we refer to cases 1 and 6 as ‘healthy state’. Note that, for these two cases, the ciliary density (i.e. cilium/area taken up for each cilium) is kept constant. On the other hand, cases 2–5 are referred to as ‘diseased state’, which is
associated with acquired and inherited respiratory disorders; see the footnote of Table 1 for more information on common disorders and their associated abnormalities. Taking inspiration from an area of ciliary activity on cultured rabbit tracheal epithelium [48], the initial arrangements of the cilia are chosen such that there is a time lag of 𝑇 ∕12 between successive cilia. The selection of the rabbit lung for the ciliary arrangements is due to its similarities to humans in terms of anatomy, physiology, genetic, and biochemical similarities [58]. Here, a ciliated domain including 96 cilia is assumed. The direction in which the metachronal wave travels is chosen to be diagonal, at an angle of 135◦ to the effective stroke. 4
Computers in Biology and Medicine 117 (2020) 103595
S.M. Vanaki et al.
Fig. 7. Averaged linear distance vs path distance travelled by a patch of particles during 250 cycles. Note that all the values in cases 2–6 have been normalised by the linear distance obtained from the standard case. The error bar is an indication of particle dispersion intensity in the flow field.
For the standard case, initially, a patch of 16×16 tracers (coloured in red and blue in 50:50 ratio) is introduced above the tips of the cilia, and the particles are tracked for 250 beating cycles, as shown in Fig. 5. With a slight deviation towards the 𝑦-direction, the seeded particles have a significant forward transport in the same direction as the effective cilia motion, i.e. the 𝑥-direction. Moreover, the majority of particles have remained above the cilia tip (close to the upper boundary), and only 0.78% of them were deposited after the prescribed period of time. Such particle transport behaviour may be viewed from two different aspects: (i) toxic agents that can pass the first defencive barrier in the healthy lung are more likely to be transported back to the mucus gel and be cleared via the muco-ciliary mechanism. This shows the key role
of ciliary actions in protecting our respiratory system against harmful submicron substances — thus hindering the rate of their deposition onto the epithelium; (ii) the active pharmaceutical ingredients may not easily reach to the epithelium and eventually the blood stream. This finding is of importance given that the pulmonary tissues are a noninvasive absorption route for the air-to-blood drug delivery – not only for local treatment of lung diseases, but also to treat other diseases such as central nervous system (CNS) disorders [59]. The colour-coded particle tracking is also carried out for other simulation cases and compared against the standard case, as shown in Fig. 6. The rate of particle deposition is provided for each case. The lower beating frequency (case 2), which is a common defect 5
Computers in Biology and Medicine 117 (2020) 103595
S.M. Vanaki et al. Table 1 Simulation cases in the present work. Simulation cases Model parameters
#1 (standard) #2 #3 #4 #5 #6
State
Frequency Length Beat pattern
Packing
10 Hz 5 Hza 10 Hz 10 Hz 10 Hz 10 Hz
Triangular ✓ Triangular Triangular Triangular Triangular Square ✓
7 7 4 7 7 7
μm μm 𝛍𝐦b μm μm μm
3D 3D 3D Stiff planarc Wiperd 3D
Healthy Diseased ✓ ✓ ✓ ✓
a
Reduced ciliary beat frequency – e.g. primary ciliary dyskinesia (PCD) [49]; asthma [50]; cigarette smoking [51]; environmental pollutants [52]. b Shorter cilia – e.g. cigarette smoking [53,54]; chronic obstructive pulmonary disease (COPD) [55]. c & d Abnormal ciliary motions – e.g. PCD [56,57].
in PCD [49] and asthma patients [50], cigarette smokers [51], and individuals exposed to severe air pollution [52], may decrease the rate of particle deposition. Additionally, similar to the standard case, most of the particles have a tendency to stay above the cilia tips. Such behaviour is much more significant in case 3, where cilia are short in length. This can be expected because the location where the tracers are released is not directly perturbed by ciliary actions. Therefore, the chance of delivering drugs to the epithelium of COPD patients and smokers, who have shorter cilia (see Table 1), may be notably reduced. Primary ciliary dysfunction (PCD) can cause abnormal ciliary beat patterns in the lungs — such as stiff planar and wiper cilia motions [56,57]. Prescribing these two abnormal cilia beats in the model significantly hindered the transport of particles. For the stiff planar beat (case 4), a significant number of particles can surprisingly move in the negative 𝑥-direction from the release point, with 1.17% of them being deposited. For the stiff planar beat (case 5), the particle tracers do not experience a significant transport in the 𝑥-direction, and move close to one another. From a pharmaceutical point of view, this would prolong the presence of drug particles in the desired site of action. When the cilia were arranged in a square pack in case 6, a deposition rate of 4.69% was observed and the particles experienced a higher mixture compared to the standard case. By looking at healthy state results, it is evident that some measures (such as modification of the particle surface [60]) must be taken to avoid immediate clearance of the particles and to prolong their residence time at the site of absorption. To quantify these statements, we compute both linear and path distances that a patch of particles migrates during 250 cycles, and present results for all the simulation cases in Fig. 7. The linear distance is obtained by the projection of particle location on the 𝑥-axis at every cycle of simulation, as this is the primary direction that particles are cleared from the respiratory tracts or may move to the lower airways. Using particle trajectories, the path distance can be calculated as √ ( ) ( ) ( ) 𝑑𝑥𝑝 2 𝑑𝑦𝑝 2 𝑑𝑧𝑝 2 𝑃 𝑎𝑡ℎ 𝑑𝑖𝑠𝑡𝑎𝑛𝑐𝑒 = + + 𝑑𝑡 (4) ∫ 𝑑𝑡 𝑑𝑡 𝑑𝑡
Fig. 8. Sensitivity tests. The effect of (a) patch size of particles, and (b) different releasing time of particles. 𝑇 denotes the period of beating cycle, and the results have been provided for 100 cycles for the standard case.
velocity profiles are constructed by plotting the time-averaged velocity over the stable cycle at each height 𝑧 of the cube located at the centre of the ciliated domain (see Fig. 9a), which can be written as ⟨𝐮⟩ = (1∕𝑇 𝐴𝐿 )
∫
𝐮𝑑𝐴𝐿 𝑑𝑡,
(5)
where the cube is of cross section 𝐴𝐿 = 8 × 8 μm2 and of height 𝑧 = 8 μm. Notice that because the cilia have abnormal beat patterns in diseased states, variations in the flow field should be analysed threedimensionally. Therefore, three representative planes are chosen to visualise the velocity contours — which are the blue, red, and grey sides of the cube in Fig. 9a. Overall, the 𝑥-component of velocity ⟨𝑢⟩ = (1∕𝑇 𝐴) ∫ 𝑢𝑑𝐴𝑑𝑡 has the highest value in magnitude for all the simulation cases, except for stiff planar and wiper ciliary patterns. For these two, the 𝑦-component of velocity ⟨𝑣⟩ is dominant, with only small changes in ⟨𝑢⟩. This signifies the severity of diseased lungs under such abnormal beating conditions, as the net axial transport of deposited foreign particles could be greatly impaired. However, there would be a high chance for therapeutic particles to be absorbed by the epithelium, due to the high residence time of particles. While the calculated values of ⟨𝑤⟩ (i.e. the 𝑧-component of velocity) were positive for most of the cases, the square pack led to negative values along the 𝑧-axis. This may be surprising, and simply means that the dispersed particles would be transported to the epithelial surface upon deposition in case 6, whereas the particles could experience upwards transport in the other cases. To summarise the findings discussed in this section, we have looked at the variations of the vector fields over one beating cycle for each simulation case, and a schematic diagram of fluid flow has been interpreted and shown on 𝑥-𝑧 plane in Fig. 10. This could qualitatively be used to
It was seen that the linear distance halved with decreasing the ciliary beating frequency from 10 to 5 Hz. Reducing the length of cilia from 7 to 4 μm led to approximately 70% reduction in the linear transport of particles. The results for cases 4 and 5 reveal that the particles have experienced a relatively large displacement, which must be in the 𝑦-direction, as the linear distance has been quite small. Finally, the square pack resulted in a 40% reduction in linear distance as compared with the standard case. Further to these results, two sensitivity tests were performed, and the results showed that the effect of patch size and releasing time of particles on their transport is insignificant, as seen in Fig. 8. The role of abnormal beats in particle transport can be better understood by having the velocity field information. As such, the 𝑥-, 𝑦-, and 𝑧-components of the velocity profile and velocity contour, for each simulation case, are shown in Figs. 9b and 9c, respectively. The 6
Computers in Biology and Medicine 117 (2020) 103595
S.M. Vanaki et al.
Fig. 9. (a) A square column with the cross sectional area of 𝐴𝐿 ∼ 8 × 8 μm2 is shown in the middle of the computational domain through which (b) the time-averaged velocity profiles in the 𝑥-, 𝑦-, and 𝑧-directions are computed. (c) Time-averaged velocity contours of ⟨𝑢⟩, ⟨𝑣⟩, and ⟨𝑤⟩ are shown. All the plots are presented for the simulation cases 1–6 (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.).
7
Computers in Biology and Medicine 117 (2020) 103595
S.M. Vanaki et al.
an advection dominant PCL flow for the healthy case, where the typical ciliary beating frequency is around 10 Hz, the ciliary waves are antiplectic, and the diffusivity of nanoparticles 𝐷 in a watery PCL flow is of order of ≈ 10−12 . Similarly, Ding et al. [32] and Chateau et al. [27] reported that fluid advection is the dominant mechanism for the particle transport inside the PCL. This suggests the necessity of modifying the physico-chemical properties of the nanocarriers to increase their absorption rate. We conclude by noting that the present model is a basic one and effects such as the particle Brownian motion, particle–cilia as well as particle–particle interactions, and particle charges may be implemented in future models. Declaration of competing interest The authors confirm that there are no known conflicts of interest associated with this publication and there has been no significant financial support for this work that could have influenced its outcome. Acknowledgements Shayan M. Vanaki gratefully acknowledges Queensland University of Technology (QUT) for support through the QUTPRA scholarship. The authors appreciatively acknowledge the High Performance Computing Unit of QUT. References [1] S. Khalafvand, E. Ng, L. Zhong, CFD simulation of flow through heart: a perspective review, Comput. Methods Biomech. Biomed. Eng. 14 (01) (2011) 113–132. [2] D.M. Sforza, C.M. Putman, J.R. Cebral, Computational fluid dynamics in brain aneurysms, Int. J. Numer. Methods Biomed. Eng. 28 (6–7) (2012) 801–808. [3] M. Ferrua, R. Singh, Modeling the fluid dynamics in a human stomach to gain insight of food digestion, J. Food Sci. 75 (7) (2010) R151–R162. [4] W. Hofmann, Modelling inhaled particle deposition in the human lung—A review, J. Aerosol Sci. 42 (10) (2011) 693–724. [5] P.W. Longest, K. Bass, R. Dutta, V. Rani, M.L. Thomas, A. El-Achwah, M. Hindle, Use of computational fluid dynamics deposition modeling in respiratory drug delivery, Expert Opinion Drug Deliv. 16 (1) (2019) 7–26. [6] M.R. Knowles, R.C. Boucher, Mucus clearance as a primary innate defense mechanism for mammalian airways, J. Clinic. Investigation 109 (5) (2002) 571–577. [7] D.J. Thornton, J.K. Sheehan, From mucins to mucus: toward a more coherent understanding of this essential barrier, Proc. Amer. Thoracic Soc. 1 (1) (2004) 54–61. [8] M. Geiser, L.M. Cruz-Orive, V.I. Hof, P. Gehr, Assessment of particle retention and clearance in the intrapulmonary conducting airways of hamster lungs with the fractionator, J. Microsc. 160 (1) (1990) 75–88. [9] M. Geiser, W.G. Kreyling, Deposition and biokinetics of inhaled nanoparticles, Particle Fibre Toxicol. 7 (1) (2010) 2. [10] C. Yu, J. Hu, B. Yen, D. Spektor, M. Lippmann, Models for mucociliary particle clearance in lung airways, in: Aerosols: Research, Risk Assessment and Control Strategies, Chelsea, MI, Lewis, 1986, pp. 569–578. [11] P.S. Lee, T.R. Gerrity, F.J. Hass, R.V. Lourenco, A model for tracheobronchial clearance of inhaled particles in man and a comparison with data, IEEE Trans. Biomed. Eng. (11) (1979) 624–630. [12] R. Cuddihy, H. Yeh, Respiratory tract clearance of particles and substances dissociated from particles, in: Inhalat. Toxicol., Springer, 1988, pp. 169–193. [13] W. Hofmann, L. Koblinger, Stochastic model of radon daughter deposition and clearance in human bronchial airways, in: IRPA9: 1996 International Congress on Radiation Protection. Proceedings, vol. 2, 1996. [14] B. Asgharian, W. Hofmann, F. Miller, Mucociliary clearance of insoluble particles from the tracheobronchial airways of the human lung, J. Aerosol Sci. 32 (6) (2001) 817–832. [15] R. Sturm, W. Hofmann, Stochastic modeling predictions for the clearance of insoluble particles from the tracheobronchial tree of the human lung, Bull. Math. Biol. 69 (1) (2007) 395. [16] J. Blake, H. Winet, On the mechanics of muco-ciliary transport, Biorheology 17 (1–2) (1980) 125–134. [17] D. Smith, E. Gaffney, J. Blake, Erratum: A viscoelastic traction layer model of muco-ciliary transport, Bull. Math. Biol. 69 (4) (2007) 1449. [18] P. Kurbatova, N. Bessonov, V. Volpert, H. Tiddens, C. Cornu, P. Nony, D. Caudri, C.W. Group, et al., Model of mucociliary clearance in cystic fibrosis lungs, J. Theoret. Biol. 372 (2015) 81–88.
Fig. 10. Schematic diagram of fluid flow (𝑥-𝑧 plane view) in the simulation cases 1–6. The arrows represent the flow direction, indicating in which directions particles could transport when entering the cilia-induced flow field. The chosen plane is approximately located in the middle of the ciliated domain.
predict the possible directions in which particles might transport when entering the cilia-induced flow fields. 4. Conclusions The major objective of the present work was to better understand the fate of biodegradable nanoparticles, used frequently as drug delivery vehicles, in the presence of the ciliated cells where their transport can be significantly altered by ciliary beats. To achieve this, we set up a computational model based on the projection method combined with the IB and particle track methods to simulate the transport of inhaled particles deposited onto the periciliary liquid layer in healthy and diseased lungs. A detailed study of the particle transport induced by abnormal ciliary beats was performed, and the results showed that impairments such as low ciliary beat frequency and short ciliary length (seen in PCD patients and cigarette smokers) can significantly reduce the rate of transport of therapeutic particles in the proximal direction of the tracheobronchial tree. The configurations corresponding to stiff planar and wiper ciliary motions (seen on PCD patients) were found to greatly reduce the transport of such particles. These findings show that the site-specific delivery of drugs may be achieved in diseased states. However, the beating of cilia creates a barrier for the particles to easily reach to the epithelial cells. Our quantitative results also reveal that the chance of drug retention and absorption by pulmonary tissues may be markedly diminished in healthy lungs if drugs are intended to be delivered into the blood stream by aerosol carriers, specifically in the upper airways. However, such observation might be changed by the diffusive effect of particles in a cilia-induced flow field which has not been implemented in the present model. Nonetheless, based on the flow rate changes during one cycle (shown in [39]) for the present model, it is likely to have 8
Computers in Biology and Medicine 117 (2020) 103595
S.M. Vanaki et al.
[42] G. Haller, T. Sapsis, Where do inertial particles go in fluid flows? Physica D 237 (5) (2008) 573–583. [43] W.S. Cheow, M.L.L. Ng, K. Kho, K. Hadinoto, Spray-freeze-drying production of thermally sensitive polymeric nanoparticle aggregates for inhaled drug delivery: effect of freeze-drying adjuvants, Int. J. Pharm. 404 (1–2) (2011) 289–300. [44] A.L. Buchmann, L.J. Fauci, K. Leiderman, E.M. Strawbridge, L. Zhao, Flow induced by bacterial carpets and transport of microscale loads, in: Applications of Dynamical Systems in Biology and Medicine, Springer, 2015, pp. 35–53. [45] E. Herawati, D. Taniguchi, H. Kanoh, K. Tateishi, S. Ishihara, S. Tsukita, Multiciliated cell basal bodies align in stereotypical patterns coordinated by the apical cytoskeleton, J. Cell Biol. 214 (5) (2016) 571–586. [46] L. Gheber, Z. Priel, Extraction of cilium beat parameters by the combined application of photoelectric measurements and computer simulation, Biophys. J. 72 (1) (1997) 449–462. [47] D.J. Lubkin, E.A. Gaffney, J.R. Blake, A viscoelastic traction layer model of muco-Ciliary transport, Bull. Math. Biol. (ISSN: 1522-9602) 69 (1) (2006) 289. [48] M. Sanderson, M. Sleigh, Ciliary activity of cultured rabbit tracheal epithelium: beat pattern and metachrony, J. Cell Sci. 47 (1) (1981) 331–347. [49] C.M. Rossman, R.M. Lee, J.B. Forrest, M.T. Newhouse, Nasal ciliary ultrastructure and function in patients with primary ciliary dyskinesia compared with that in normal subjects and in subjects with various respiratory diseases, Amer. Rev. Respiratory Disease 129 (1) (1984) 161–167. [50] B. Thomas, A. Rutman, R.A. Hirst, P. Haldar, A.J. Wardlaw, J. Bankart, C.E. Brightling, C. O’callaghan, Ciliary dysfunction and ultrastructural abnormalities are features of severe asthma, J. Allergy Clin. Immunol. 126 (4) (2010) 722–729. [51] A.E. Tilley, M.S. Walters, R. Shaykhiev, R.G. Crystal, Cilia dysfunction in lung disease, Annu. Rev. Physiol. 77 (2015) 379–406. [52] H. Riechelmann, K. Kienast, J. Schellenberg, W. Mann, An in vitro model to study effects of airborne pollutants on human ciliary activity, Rhinology 32 (3) (1994) 105–108. [53] P.L. Leopold, M.J. O’Mahony, X.J. Lian, A.E. Tilley, B.-G. Harvey, R.G. Crystal, Smoking is associated with shortened airway cilia, PLoS One 4 (12) (2009) e8157. [54] J. Hessel, J. Heldrich, J. Fuller, M.R. Staudt, S. Radisch, C. Hollmann, B.-G. Harvey, R.J. Kaner, J. Salit, J. Yee-Levin, et al., Intraflagellar transport gene expression associated with short cilia in smoking and COPD, PLoS One 9 (1) (2014) e85453. [55] H.C. Lam, S.M. Cloonan, A.R. Bhashyam, J.A. Haspel, A. Singh, J.F. Sathirapongsasuti, M. Cervo, H. Yao, A.L. Chung, K. Mizumura, et al., Histone deacetylase 6–mediated selective autophagy regulates COPD-associated cilia dysfunction, J. Clinic. Investigation 123 (12) (2013) 5212–5230. [56] C.M. Rossman, J.B. Forrest, R.M. Lee, M.T. Newhouse, The dyskinetic cilia syndrome: ciliary motility in immotile cilia syndrome, Chest 78 (4) (1980) 580–582. [57] M.A. Chilvers, A. Rutman, C. O’callaghan, Ciliary beat pattern is associated with specific ultrastructural defects in primary ciliary dyskinesia, J. Allergy Clin. Immunol. 112 (3) (2003) 518–524. [58] N.A. Kamaruzaman, E. Kardia, N. Kamaldin, A.Z. Latahir, B.H. Yahaya, The rabbit as a model for studying lung disease and stem cell therapy, BioMed Res. Int. 2013 (2013). [59] P. Noymer, S. Biondi, D. Myers, J. Cassella, Pulmonary delivery of therapeutic compounds for treating CNS disorders, Therapeutic Deliv. 2 (9) (2011) 1125–1140. [60] J.C. Sung, B.L. Pulliam, D.A. Edwards, Nanoparticles for drug delivery to the lungs, Trends Biotechnol. 25 (12) (2007) 563–570.
[19] G.R. Fulford, J.R. Blake, Muco-ciliary transport in the lung, J. Theoret. Biol. 121 (4) (1986) 381–402. [20] W. Lee, P. Jayathilake, Z. Tan, D. Le, H. Lee, B. Khoo, Muco-ciliary transport: effect of mucus viscosity, cilia beat frequency and cilia density, Comput. & Fluids 49 (1) (2011) 214–221. [21] R. Chatelin, P. Poncet, A hybrid grid-particle method for moving bodies in 3D stokes flow with variable viscosity, SIAM J. Sci. Comput. 35 (4) (2013) B925–B949. [22] P. Jayathilake, D. Le, Z. Tan, H.P. Lee, B.C. Khoo, A numerical study of mucociliary transport under the condition of diseased cilia, Comput. Methods Biomech. Biomed. Eng. 18 (9) (2015) 944–951. [23] M. Sedaghat, M. Shahmardan, M. Norouzi, P. Jayathilake, M. Nazari, Numerical simulation of muco-ciliary clearance: immersed boundary-lattice Boltzmann method, Comput. & Fluids 131 (2016) 91–101. [24] R. Chatelin, P. Poncet, A parametric study of mucociliary transport by numerical simulations of 3D non-homogeneous mucus, J. Biomech. 49 (9) (2016) 1772–1780. [25] H. Guo, E. Kanso, A computational study of mucociliary transport in healthy and diseased environments, Eur. J. Comput. Mech. 26 (1–2) (2017) 4–30. [26] R. Chatelin, D. Anne-Archard, M. Murris-Espin, M. Thiriet, P. Poncet, Numerical and experimental investigation of mucociliary clearance breakdown in cystic fibrosis, J. Biomech. 53 (2017) 56–63. [27] S. Chateau, U. D’Ortona, S. Poncet, J. Favier, Transport and mixing induced by beating cilia in human airways, Front. Physiol. 9 (2018) 161. [28] Y.L.R. Quek, K.M. Lim, K.-H. Chiam, Three-dimensional computational model of multiphase flow driven by a bed of active cilia, Comput. & Fluids 170 (2018) 222–235. [29] R.H. Dillon, L.J. Fauci, C. Omoto, X. Yang, Fluid dynamic models of flagellar and ciliary beating, Ann. New York Acad. Sci. 1101 (1) (2007) 494–505. [30] S.M. Mitran, Metachronal wave formation in a model of pulmonary cilia, Comput. Struct. 85 (11–14) (2007) 763–774. [31] S. Lukens, X. Yang, L. Fauci, Using Lagrangian coherent structures to analyze fluid mixing by cilia, Chaos 20 (1) (2010) 017511. [32] Y. Ding, J.C. Nawroth, M.J. McFall-Ngai, E. Kanso, Mixing and transport by ciliary carpets: a numerical study, J. Fluid Mech. 743 (2014) 124–140. [33] R. Ghosh, G.A. Buxton, O.B. Usta, A.C. Balazs, A. Alexeev, Designing oscillating cilia that capture or release microscopic particles, Langmuir 26 (4) (2009) 2963–2968. [34] A. Bhattacharya, G.A. Buxton, O.B. Usta, A.C. Balazs, Propulsion and trapping of microparticles by active cilia arrays, Langmuir 28 (6) (2012) 3217–3226. [35] Y. Ding, E. Kanso, Selective particle capture by asynchronously beating cilia, Phys. Fluids 27 (12) (2015) 121902. [36] C.A. Ruge, J. Kirch, C.-M. Lehr, Pulmonary drug delivery: from generating aerosols to overcoming biological barriers—therapeutic possibilities and technological challenges, Lancet Respiratory Med. 1 (5) (2013) 402–413. [37] X. Murgia, P. Pawelzyk, U.F. Schaefer, C. Wagner, N. Willenbacher, C.-M. Lehr, Size-limited penetration of nanoparticles into porcine respiratory mucus after aerosol deposition, Biomacromolecules 17 (4) (2016) 1536–1542. [38] S. Tsorng, H. Capart, J. Lai, D. Young, Three-dimensional tracking of the long time trajectories of suspended particles in a lid-driven cavity flow, Exp. Fluids 40 (2) (2006) 314–328. [39] S.M. Vanaki, D. Holmes, P.G. Jayathilake, R. Brown, Three-dimensional numerical analysis of periciliary liquid layer: ciliary abnormalities in respiratory diseases, Appl. Sci. 9 (19) (2019) 4033. [40] M. Maxey, The gravitational settling of aerosol particles in homogeneous turbulence and random flow fields, J. Fluid Mech. 174 (1987) 441–465. [41] K.-K. Tio, A. Liñán, J.C. Lasheras, A.M. Gañán-Calvo, On the dynamics of buoyant and heavy particles in a periodic stuart vortex flow, J. Fluid Mech. 254 (1993) 671–699.
9