An adaptable model for growth and/or shrinkage of droplets in the respiratory tract during inhalation of aqueous particles

The site of deposition of pulmonary delivered aerosols is dependent on the aerosol's droplet size distribution, which may change during inhalation. The aim of this study was to develop a freely accessible and adaptable model that describes the growth (due to condensation) and shrinkage (due to evaporation) of inhaled droplets as a function of the distance from the airway wall during various inhalation conditions, for a laminar flow scenario. This was achieved by developing a model with which the evaporation of water from a droplet surface or condensation of water onto the droplet surface can be calculated. This model was then applied to a second model that describes the heat and mass transfer from the airway wall to the inhaled aerosol. The latter was based on the Weibel model. It was found that the growth and shrinkage of inhaled droplets markedly differs, depending on the distance from the airway wall. Droplets near the wall start to grow immediately due to fast water vapor transfer from the wall to the cold inhaled air. This growth continues until the air reaches body temperature and is fully saturated. However, droplets in the center of the airway first evaporate partly, due to a delay in water vapor transfer from the airway wall, before they start to grow. Depending on the conditions during inhalation, the droplet size distribution can widen considerably, which may affect the lung deposition and thereby the efficacy of the inhalation therapy. In conclusion, the model was able to show the effect of the conditions in the respiratory tract on the growth and shrinkage of inhaled droplets during standard inhalation conditions. Future developments can be aimed at expanding the model to include turbulent flow and hygroscopic growth, to improve the accuracy of the model and make it applicable to both droplets of solutions and dry particles. & 2015 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Inhalation of drugs is widely applied in the treatment of asthma and COPD. Furthermore, the pulmonary route is currently targeted for the treatment of diseases that are not directly located or even associated with the respiratory tract. Due to the proximity of the capillaries and blood vessels with the lining in the respiratory tract and the leaky nature of the pulmonary membranes (especially in the alveoli), the absorption of active pharmaceutical ingredients can be extremely rapid (Patton, 1996;Rabinowitz et al., 2006). However, the rate and amount of absorption heavily depend on the location of deposition of the inhaled substance (Labiris & Dolovich, 2003;Usmani, Biddiscombe & Barnes, 2005;Zanen, Go & Lammers, 1994).
During a standard inhalation procedure, the largest particles (either dry powder or droplets, although in this paper only aqueous droplets will be considered) are deposited in the mouth and the back of the throat, meaning their effect is essentially lost. The rest of the particles are deposited throughout the respiratory tract, where the spread of deposition is mainly dependent on the particle size distribution of the liquid aerosol or dry powder. This aerosol's particle size distribution is ideally tailored to obtain deposition at a location where it is desired for a specific application. As a general rule of thumb, particles with an aerodynamic diameter larger than 5 mm are deposited in the throat and smaller than 1 mm are exhaled after inhalation (Labiris & Dolovich, 2003). Anything in between is deposited in the lung, with smaller particles having a higher chance to reach the peripheral parts of the lung and even the alveoli and larger particles being deposited more centrally in the respiratory tract. However, in reality it is more complex.
When particles are inhaled, they are subjected to quickly changing conditions. When the particles exit the inhaler together with the ambient air, conditions surrounding the particles closely match those of the ambient air, which means relatively low relative humidity and room temperature. During the passage through the respiratory tract, the relative humidity quickly increases as well as the temperature. For both dry powder and liquid aerosols, this could affect their aerodynamic diameter, as the increased humidity could cause the dry powder particles to get wet and grow, while droplets in the liquid aerosol could grow due to condensation or shrink due to evaporation. If it would be possible to tell whether, and if so, how much the aerodynamic diameter of inhaled particles is changed while traversing the respiratory tract, it could aid in optimizing formulations for inhalation. By taking into account growth and/or shrinkage of the particles, the deposition location in the respiratory tract could be further tailored for each treatment.
Indeed, a lot of effort has gone into modeling the condensational growth of solid submicron particles during inhalation to allow targeted deposition in the respiratory tract, mainly by Longest et al. (Ferron, 1977;Golshahi et al., 2013;Longest, McLeskey & Hindle, 2010;Longest, Tian, Li, Son & Hindle, 2012;Tian, Longest, Li & Hindle, 2013;Tian, Longest, Su & Hindle, 2011;Worth Longest et al., 2012) but others as well (Broday & Georgopoulos, 2001;Kleinstreuer & Zhang, 2010). In addition, a lot is known about the conditions to which inhaled particles and droplets are subjected when inhaled (Ferron, 1977;Olson, Sudlow, Horsfield & Filley, 1973;Tian, Longest, Su & Hindle, 2011). Most of these studies revolve around the use of computational fluid dynamics (CFD), which is an excellent method to accurately model the complex flow profiles that exist inside the airways. However, such models lack openness and accessibility for researchers who are not directly involved in such fields. The same can be said for more analytical models that use pre-formulated sets of equations that describe, for example, the evaporation rate from a droplet. Furthermore, the complete models are not readily available online. Although simpler models using standard heat and mass balances usually result in less detailed results, their use often allows more openly accessible models to be developed that can be used and adapted by others. In addition, given enough effort, the complexity of the model can be increased to obtain more accurate results. Furthermore, we believe that understanding the basic mechanics underlying such a model can also greatly increase the understanding of the system under consideration.
Therefore, the aim of this study is to investigate the behavior, i.e. growth and shrinkage, of inhaled particles in the respiratory tract with a theoretical model based on standard heat and mass balances. This will be done by primarily focusing on aqueous particles. Although dissolved components such as drugs and excipients will be left out of the model to prevent it from becoming overly complex, the model will be constructed such that it can be expanded by anyone who wishes to do so. Therefore, the model will be developed in an open source platform called GNU Octave and be made available freely. Consequently, a model will be obtained that can form a basis for further expansion and development to create a more complex but also more complete model of dry particle and droplet behavior in the respiratory tract without resolving to CFD. Finally, the clinical relevance of the results of this study towards inhalation practices will also be discussed.

Open source software
The model was developed entirely in GNU Octave, which is freely available (Eaton & Community, 2014). Furthermore, the developed model itself is available online as Supplementary data.

Model input parameters
The basis for the input of the models was obtained from the Respimat s soft mist inhaler (Boehringer Ingelheim, Ingelheim, Germany). The inhaled jet that is ejected from the inhaler contains about 15 mg of water and lasts for about 1.2 s. Furthermore, the starting droplet diameter is 4.5 mm (Dalby, Eicher & Zierenberg, 2011). For the simulations, an inhalation flow of 40 L/min was chosen, although for comparison, a fast inhalation with an inhalation flow of 100 L/min will also be considered (Brand, Hederer, Austen, Dewberry & Meyer, 2008). Finally, a room temperature of either 293 K or 303 K, and a relative humidity of either 35, 50 or 70% will be used. Note that all of these parameters can be easily customized in the developed model, if desired.
Using these parameters, the volume of air that is surrounding each droplet can be estimated. To do so, it is both assumed that each droplet will be of equal size and separated by an equal distance from each other and therefore, the volume of air surrounding each droplet will be the same. To simplify the model, the layer of air around each droplet is considered to be a spherical shell instead of cubical. The volume of air around each droplet is then simply calculated by dividing the amount of air inhaled during the ejection of the jet from the inhaler by the amount of droplets generated. Furthermore, the droplet is assumed to move at the same velocity as the surrounding air, which implies that the air surrounding each droplet can be considered as stagnant.

Modeling strategy
The complete model for predicting the behavior of droplets in the respiratory tract consists of two separate models that work together. The first (droplet) model describes the mass transfer of water and energy from and to a water droplet, which will be described in the first section. The second (respiratory tract) model will describe the conditions inside the respiratory tract, and will also be based on mass and energy balances. With the respiratory tract model, the exchange of heat and mass between the airway wall and the inhaled air is first calculated. The changes in relative humidity and temperature are then passed on to the droplet model to calculate the exchange of heat and mass between the droplet and the inhaled air. By keeping the models separated, it is very easy to vary the amount of droplets that are considered during the calculation. The development of the respiratory tract model will be described in further detail in the second section.

Droplet model
The droplet model was developed as a 1-D model, meaning that the system is completely symmetrical and changes in morphology are not considered.
The primary aim of the droplet model is to calculate the changes in droplet radius as a function of the conditions to which the droplet is exposed inside the respiratory tract. To achieve this, the basic equation shown in Eq. (1) was used to start the calculation.
where r t is the radius of the droplet at time t, m t À 1 is the mass of liquid at the previous time point, Φ m is the mass flow rate from the droplet to the environment (or vice versa when negative), Δt is the time step for the iterative calculation, and ρ is the density of the liquid (1000 kg/m 3 for water, assumed to be constant). Because the mass flow rate (Φ m ) is still unknown, Eq.
(2) (derived from Fick's first law) has to be used to calculate the rate of mass flow from or towards the droplet.
where D is the diffusion coefficient of water molecules in the gas phase, M w is the molecular mass of the liquid, Δp is the pressure difference between the partial water vapor pressure at the surface of the droplet and the partial water vapor pressure in the gas phase, R is the gas constant, T is the temperature at the surface of the droplet (K), r d is the radius of the droplet and r a is the total radius of the droplet and gas phase surrounding the droplet. Here, it is assumed that the partial water vapor pressure at the surface of the droplet is equal to the saturation pressure. Therefore, the pressure gradient, Δp, can be calculated with both the Antoine equation and the ideal gas law (Eqs. (3) and (4), respectively).
where p sat is the saturated vapor pressure (Pa), A, B and C are the Antoine constants (for water these are 10.20, 1730.63, and À39.72, respectively (Yaws & Yang, 1989)), T is the temperature (K), p vap is the partial water vapor pressure (Pa), x is the specific humidity (defined as the mass of water divided by the mass of dry air) and ρ g is the density of the gas phase. Besides the partial water vapor pressure gradient, the above equations can also be used to calculate the relative humidity (RH), which is defined as p vap /p sat Á 100. The temperature in the equations for the mass flow rate (Eq. (2)) and both the vapor pressures (Eqs. (3) and (4)) can be calculated with Eq. (5).
where T t À 1 is the temperature at the previous time point, Q is the heat flow rate, C p is the specific heat capacity of the considered phase and m is the mass. The heat flow rate consists of heat flow through conduction and evaporation, which can be calculated using Eq. (6) (derived from Fourier's law) and Eq. (7), respectively. Note that Eqs.
(2) and (6) are similar due to the fact that the base equations for mass and heat transfer (Fick's law and Fourier's law, respectively) are similar as well.
where λ is the thermal conductivity (W/m/K), ΔT is the temperature difference between the temperature at the surface of the droplet and the temperature at the boundary of the gas layer surrounding the droplet and H evap is the heat of evaporation of the liquid. Note that the heat of evaporation is equal to the negative of the heat of condensation. Whether or not the heat is added to (condensation) or subtracted from (evaporation) the heat of the droplet purely depends on whether mass is flowing towards or from the droplet, respectively. By iterating Eqs.
(1)-(7), the radius of the droplet in time can be calculated as a function of the conditions to which the droplet is subjected.

Respiratory tract model
To simulate the conditions inside the respiratory tract to which the droplet is subjected, a second model was developed. Whereas the droplet model is used to calculate the exchange of heat and water between a droplet and the surrounding gas phase, the respiratory tract model calculates the exchange of heat and water between the airway wall and the inhaled mist. Furthermore, because the dimensions of the airways that a droplet is passing through changes with each generation (sections of the airway between bifurcations, with the trachea being generation zero), these are taken into account as well.
These dimensions are given by the Weibel lung model, and are shown together with the cumulative residence time of the droplets in Table 1 (Weibel, 1963).
In this respiratory tract model, the heat and mass transfer that occurs during passage from the inhaler through the mouth to the beginning of the trachea is not considered. This is done to simplify the model, since the dimensions of the mouth during inhalation are not well defined. The aerosol that is inhaled will be at ambient conditions, meaning room temperature and RH, and will start to heat up the moment it reaches the trachea. During the travel through the airways heat and mass is transferred from the airway wall to the center of the airway, increasing both the temperature and the RH of the air surrounding the droplets. This can then lead to evaporation or condensation from or onto the droplets, until equilibrium is reached or droplets are evaporated. In reality, heat and water will also be transferred to the inhaled air during travel to the trachea, which will affect the generation where the equilibrium is reached Olson, Sudlow, Horsfield & Filley, 1973;Tian, Longest, Su & Hindle, 2011).
For this respiratory tract model, the airway is considered to be a smooth cylinder consisting of concentric cylinders to divide the area inside the airway (see Fig. 1). The content of each of the cylinders are considered to be ideally mixed, meaning that there is no variation in conditions within each cylinder at any time point. This will result in a stepwise Table 1 Respiratory tract dimensions according to the Weibel lung model (Weibel, 1963).

Generation
Diameter ( gradient of temperature and RH from the airway wall to the center of the airway, with the number of steps equal to the number of chosen concentric cylinders. Furthermore, at each bifurcation (when a new generation is reached) the dimensions are changed such that the same number of cylinders will be kept the same, together with the conditions in each ring at that time point. Although in reality the flow might be split through the middle at each bifurcation leaving half the inner cylinder on the outside, in this model the distribution of the concentric cylinders inside the airway does not change when a bifurcation is passed. How this might affect the result will be discussed. In addition, loss of droplets due to impaction is not taken into account. Finally, the mass and heat transfer is modeled for a laminar flow, thereby neglecting the slightly turbulent flow in the first two to four stages and the effect of droplets on the flow (one-way coupling). Although this will impose an error on the final estimated behavior of droplets inside the respiratory tract, the importance and magnitude of this error will also be discussed. The calculations for the respiratory tract model are comparable to the calculations that are done in the droplet model. The mass transfer between the concentric cylinders is calculated using Eq. (8) and the heat conduction is calculated using Eq. (9), derived from Fick's law and Fourier's law, respectively.
where h is the height of the considered layer of mist moving through the respiratory tract (arbitrarily chosen as the diameter of the droplet with its surrounding gas layer), r n and r n þ 1 are the radii of two adjacent concentric cylinders n and n þ1 in the airway, respectively, and Δp and ΔT are the partial water vapor pressure and temperature of cylinder n and n þ1, respectively. Finally, the evaporation of water from the airway wall to the airway is assumed to occur through a small 1 mm boundary layer with 99.5% RH at body temperature (Ferron, 1977), and for the heat transfer to the airway, the thickness of the airway wall decreases from 2.75 mm in the trachea to about 1 mm at generation 11, after which it is assumed to remain constant (Liu, Chen, Tawhai, Hoffman & Sonka, 2009). The thermal conductivity, λ, used for the airway wall is about 0.28 W/m/K (Poppendiek et al., 1964). It is important to realize that the outer cylinder, due to the ideally mixed conditions, will not be at the body temperature from the start.

Combined model
The droplet model and the respiratory tract model are combined to calculate the growth and shrinkage of a droplet in each concentric cylinder in the respiratory tract model. For each droplet, the heat and water transfer from or onto the droplet is calculated. In addition, the heat and water transfer between each concentric cylinder in the airway is also Fig. 1. Schematic overview of the modeled airway with a top view on the left and a side view on the right. Left: Each concentric cylinder transfers mass and heat but there is no variation in conditions inside each cylinder. The arrows depict the transfer of mass and heat from the airway wall through each cylinder to the center of the airway. Right: The airway is considered to be a straight smooth cylinder without bifurcations but with sudden changes in diameter, through which a slice of the inhaled mist is followed. calculated. By adding the heat and water transfer from the airway and the droplet, the total change in temperature and relative humidity in each concentric cylinder and droplet is obtained after each time step. Furthermore, after the simulated time in each generation surpasses the residence time of that specific generation (Table 1), the dimensions of the concentric cylinders are instantly changed to that of the next generation. To keep calculation times limited, while still obtaining a good resolution of the size distribution of droplets inside the airway, the airway was divided into 30 cylinders and thus 30 droplets were modeled between the airway wall and the center of the airway.

Droplet model
A model was developed that describes the heat and mass transfer to and from a liquid droplet. To validate the accuracy of the droplet model alone (without external influences, e.g. as imposed by the respiratory tract model), the output as shown in Table 2 was compared to data found in a psychrometric chart. This chart gives information on the temperature of a droplet dispersed in air with a specific temperature and humidity. When a droplet is exposed to air with a moisture content below the dew point, the droplet will evaporate. Due to the evaporation, the temperature of the droplet will decrease to a specific temperature, the wet bulb temperature, which is dependent on the rate of evaporation. This effect is referred to as evaporative cooling. The rate of evaporation, and consequently the equilibrium temperature of the droplet, is dependent on the relative humidity, the ambient temperature, the type of liquid (in this case only water is considered), and the velocity difference between the ambient air and the droplet. For this case it was assumed that the velocity of the droplet is equal to the velocity of the surrounding air, i.e. there is a stagnant layer of air surrounding the droplet. Furthermore, the droplet diameter at the start of the simulation was 4.5 mm. Finally, the droplet temperature (T w model) was calculated at three different relative humidities, i.e. 0%, 50% and 95%, at an ambient temperature of 293 K and compared to values found in a psychrometric chart at similar conditions (T w chart). It was found that the calculated droplet temperatures were in full agreement with values found in psychrometric charts. Therefore, it can be concluded that the mass and heat transfer in the droplet model is modeled correctly and accurately. Since the same calculations were used for the mass and heat transfer Table 2 Comparison of calculated droplet temperature in stagnant air at 0%, 50% and 95% RH with the wet bulb temperature found in psychrometric charts.

RH (%)
T w model (  between the airway wall and the airway, it can be assumed that the results will also correspond well to the conditions set for this model (i.e. laminar flow).
To further show the application of the droplet model to the growth of droplets due to condensation, a simulation was run for droplets with a diameter of 0.9 mm suspended in air of 37°C with a relative humidity of 99.5, 100, 101 or 104%. The volume of air surrounding the droplets was high enough to negate any effect of evaporation and condensation of the droplets on the relative humidity. These settings were specifically chosen to be able to compare the results to data reported by Longest, McLeskey & Hindle (2010), to determine whether the rate of condensation is correctly modeled (see Fig. 2). Although the values reported by Longest et al. include the effect of hygroscopic growth (starting droplets consisted of a 0.5% w/v albuterol sulfate solution), at higher relative humidity it is expected that growth due to condensation will be high, causing a high dilution of the solution and in turn reducing the effect of hygroscopic growth. At a relative humidity of 104% the growth rate corresponds well with the values reported by Longest et al. for condensational growth of droplets containing albuterol sulfate. When the relative humidity is 101%, the effect of hygroscopic growth is noticeable as the condensational growth rate calculated with our droplet model is slightly lower. However, at a relative humidity of 100% and 99.5% the difference is much more pronounced: droplets without solutes evaporate, whereas values reported by Longest et al. still show growth due to hygroscopic effects.

Droplet growth and shrinkageconditions in the respiratory tract
The airway was divided into 30 concentric cylinders and for each cylinder the behavior of a droplet was modeled. For the default simulation, the ambient air temperature and relative humidity used were 293 K and 35%, respectively. In addition, the starting diameter for each droplet was 4.5 mm. For each of the droplets the surrounding air temperature and relative humidity are shown in time (Fig. 3). Furthermore, the diameter of each of the droplets in time is shown in relation to the generation where the droplets are located during the simulation (Fig. 4). Finally, the droplet size distribution was calculated when the system reached equilibrium and the droplet diameter did not change any further (Fig. 5).
The relative humidity and the temperature close to the wall of the airway increased very rapidly. However, due to the assumption that the flow is laminar, the temperature in the center of the airway does not increase immediately. Only after about 0.075 s (less for locations closer to the airway wall) the temperature in the center of the airway increased to body temperature. Therefore, depending on the location inside the airway, inhaled droplets are subjected to different conditions. A noticeably high relative humidity was found during the initial entry of the colder inhaled air into the trachea of about 240% closest to the wall, implying supersaturation. This quickly decreases again as the temperature increases and the water vapor diffuses to the center of the airway. During the stage where the temperature of the airway has not yet reached body temperature, the simulated relative humidity never decreased below 110%. In fact, a slight increase to about 125% is seen at the same time when the temperature of the air further from the airway wall starts to increase rapidly. This is due to a quick decrease in the diameter of the airway, by which the diffusion rate of water vapor and heat rapidly increases. Another effect that results from a change in the diameter of the airway (and thus the modeled concentric cylinders) is the slight drop in temperature that can be observed at several time points. Here, the heat diffusion towards the center is suddenly increased, decreasing the temperature near the airway wall until a new equilibrium is reached and the temperature increases again. The effect of the supersaturation, i.e. relative humidities higher than 100%, on the behavior of the droplets will be investigated further below in the results section entitled: "Droplet behavior -Maximum relative humidity", as nucleation and formation of separate water droplets which might occur by this supersaturation will affect the results shown in Figs.

3-5.
A clear influence of the location of the droplet in the airway on the droplet behavior was found. Droplets closest to the airway wall grew, while the droplets in the center of the airway first shrunk and after a while started to grow. This is clearly the result of the higher relative humidity and temperature near the airway wall and the lower relative humidity and temperature at the center of the airway. Air near the airway wall is quickly humidified and heated by the airway wall, which causes condensation on the cooler droplets. On the other hand, moisture and heat from the airway wall takes a while to reach the air in the center of the airway (farthest from the airway wall). Therefore, droplets near the center of the airway will first evaporate until the air is saturated, which is when the droplet size reaches a plateau value. Due to evaporative cooling, the temperature of the air will also decrease. When the heat and moisture from the airway wall reaches these air layers, the relative humidity and temperature increases rapidly, resulting in condensation on the droplets. Finally, no large changes in droplet size occur after generation 9-10 because equilibrium was reached, although there is a slight and steady evaporation of the droplets due to the lower relative humidity of the airway wall (99.5%).
When the actual droplet size distribution after reaching equilibrium is considered (Fig. 5), a large widening of the droplet size distribution is found, ending up with a major volume fraction of the droplets larger than the originally inhaled droplets. Under the simulated circumstances (inhaled air with 35% RH, 293 K, 40 L/min) about 63% of the total droplet volume had a diameter above 4.5 mm, which was the original size of the inhaled droplets.

Droplet behaviormaximum relative humidity
In the previous simulation, the maximum relative humidity in the airway was not limited. In reality, however, the relative humidity might be limited by nucleation and formation of separate water droplets (mist). Although in the current model, the formation of these droplets is not considered, the maximum relative humidity can easily be changed. Because the relative humidity has a large influence on the condensation and evaporation rate (and thus droplet growth and shrinking), a simulation was run to compare three relative humidity settings: unlimited (same as in Figs. 3-5), limited to 104% and limited to 100% (Fig. 6). The relative humidity of 104% was added to match that used by Longest et al. Longest, McLeskey & Hindle, 2010), in addition, it was added to serve as a scenario for modest supersaturation, i.e. between no supersaturation (100% RH) and unlimited supersaturation.
When the maximum relative humidity is limited to a maximum of 100% or 104%, a clear overall decrease of the droplet size was found when compared to droplets at unlimited maximum relative humidity. However, it is questionable whether limiting the maximum relative humidity would result in a better representation of conditions during inhalation. It is well known that the degree of supersaturation required for nucleation depends strongly on purity of the liquid and the air. Under ideal circumstances (no impurities), the relative humidity corresponding to the supersaturation required for homogeneous nucleation can be as high as 2200% (Hinds, 1999). The presence of any nuclei, either droplets or other particles, will lower the relative humidity at which condensation occurs. Therefore, in our situation, a lot of the water might actually condensate on the inhaled droplets instead of forming new droplets. In addition, by limiting the relative humidity in the model, a lot of the water that will evaporate from the airway wall to the cooler and drier inhaled air is simply ignored. In other words, all the extra water that evaporates from the airway wall nucleates and does not interact with the droplets in any way. Most likely, in reality some or perhaps most of the moisture may condense onto or agglomerate with the droplets. Finally, the time during which the relative humidity reaches values well above 100% is quite short. The initial burst in relative humidity before being dissipated by transport and condensation is less than 0.02 s (Fig. 3), and the total time before equilibrium is reached with the airway wall is about 0.12 s. Although it is unknown what time scales are involved with the formation of nuclei, it might be possible that the high relative humidities observed in Fig. 3 can be reached for such short duration, even when impurities are present in the inhaled air (besides the droplets themselves). Therefore, it is expected that the results with the unlimited model are more accurate than the model with limitations with respect to RH. However, below both unlimited and limited relative humidity will be considered to better assess the importance of different conditions during inhalation on the behavior of inhaled droplets in the respiratory tract.

Droplet behaviorinlet conditions
Because the conditions at which an inhaler is used can vary significantly, the equilibrium droplet size distribution of droplets inhaled under varying conditions was calculated (Fig. 7). The conditions chosen were an inhaled air temperature of 293 K with a relative humidity of 35%, 50%, and 70%. Furthermore, the influence of inhaled air temperature was determined with an air temperature of 303 K and a relative humidity of 35%. These conditions were simulated with a relative humidity that was either not limited, limited to a maximum of 104%, or limited to a maximum or 100%. Finally, all conditions were simulated during a slow inhalation (40 L/min) and a fast inhalation (100 L/min).  6. Comparison of equilibrium droplet size distribution with relative humidity in the airway being: limited to 100% ( À ), limited to 104% ( À À) and unlimited ( À Á À) relative humidity. The original droplet diameter at the start of the simulation was 4.5 mm, inhaled air with 35% RH, 293 K, inhaled at 40 L/min. Fig. 7. Comparison of equilibrium droplet size distribution after slow inhalation of air (40 L/min, left) or fast inhalation of air (100 L/min, right), with a relative humidity that was not limited (top), limited to a maximum of 104% (middle), or limited to a maximum of 100% (bottom). Simulated inhalation conditions were a relative humidity of 35% ( À ), 50% ( À À) and 70% ( À Á À) at 293 K and 35% at 303 K ( Á Á Á ). The original droplet diameter at the start of the simulation was 4.5 mm.
Although the ambient relative humidity did not have a great effect on the maximum droplet diameter, there was a noticeable effect on the smallest droplet size. When air with a higher relative humidity was chosen, the diameter of the smallest droplets increased. This is caused by a faster saturation of the air closer to the center of the airway when it is not yet heated by the airway wall. Therefore, the droplets evaporate less and the minimum droplet diameter increases.
When the inhaled air temperature was 303 K instead of 293 K, a large overall decrease of the droplet diameter was observed. For droplets closest to the airway wall the higher inhaled air temperature resulted in a lower super saturation of the air with moisture (data not shown). Therefore, the condensation rate onto the droplets was lower than onto droplets in air with a lower temperature. This effect was absent when the relative humidity in the airway was limited to a maximum of 104% and 100%, and the maximum droplet size does not change with a change in inhaled air temperature. Droplets closer to the center of the airway became smaller with increased inhaled air temperature, due to the higher moisture capacity. Air with a fixed relative humidity but a higher temperature will have a higher absolute moisture content before becoming saturated. Therefore, droplets can evaporate more before the air is saturated and equilibrium is reached.
When the inhalation is performed quickly, i.e. at 100 L/min instead of 40 L/min, an overall decrease in droplet size was found. In addition, the effects of changing the inhalation temperature and relative humidity were more pronounced as well. At the same relative humidity and temperature of the inhaled air, the total volume of air surrounding each droplet is increased when the inhalation flow is increased. Therefore, much more water can evaporate before the air is saturated, resulting in smaller droplets near the center of the airway. The effect is smaller near the airway wall, where it takes a bit longer for the air to be saturated and cause slightly slower condensation before equilibrium is reached. If besides the inhalation flow also the temperature or the relative humidity is changed, the effect on droplet size distribution becomes more pronounced as the change in moisture capacity due to both parameters also depends on the volume of air. The same change in temperature from 293 K to 303 K or a change in relative humidity from 35% to 70% will result in a much bigger change in moisture capacity at a high inhalation flow rate than at a lower inhalation flow rate.
When the relative humidity in the respiratory tract is limited to a maximum of 104%, there is an overall decrease in the droplet size distribution compared to the scenario where the relative humidity in the respiratory tract is unlimited. Although droplets in the center of the airway only reduce less than 0.5 mm in size, due to the larger decrease in droplet size near the airway wall, there is a noticeable reduction in the width of the droplet size distribution. Furthermore, there is a decrease of the volume fraction of droplets with a size larger than the originally inhaled droplets of 4.5 mm.
Finally, when the relative humidity was limited to 100%, no growth of droplets is observed. Instead, the droplet diameter does not change when close to the airway wall and decreases when closer to the center of the airway. Compared to the scenario where the relative humidity is limited to 104%, the width of the droplet size distribution remains the same, although the distribution shifts towards smaller droplets with no fraction larger than the original droplet size.

Discussion
It is shown that the developed droplet model can be used to accurately estimate the shrinkage or growth of a droplet by evaporation or condensation, respectively. When this droplet model is coupled to the Weibel model of the lung to obtain the respiratory tract model, estimations can be made regarding the behavior of inhaled droplets in the respiratory tract. Our simulations indicate that whether a droplet grows or shrinks is dependent on where the droplet is situated with regard to the airway wall . When close to the wall, the relatively cold droplet is quickly subjected to an environment of higher temperature and humidity which promotes condensational growth, resulting in larger droplets. However, droplets farther away from the wall will first evaporate before the heat and moisture from the airway wall will reach these droplets. These changes in droplet size distribution happen rather quickly, as it was found that at generation 9-10 the conditions inside the airway will have reached equilibrium and no further change in droplet size occurs. Therefore, there could be a substantial effect of this change in droplet size distribution on the deposition site. Instead of reaching the lower airways, the larger droplets (which are, as shown in the results, already closer to the airway wall) will be deposited higher in the respiratory tract than might be originally intended, thereby reducing the efficacy of the inhalation therapy, while the smaller droplets will deposit deeper in the respiratory tract than originally intended. Knowing the behavior of the droplets inside the respiratory tract can help in countering the effect with either changes in formulation, design of the inhalation device, or inhalation instructions.
Although the idea of droplet growth and shrinkage in the respiratory tract is not new, we have developed a basic model that enables one to estimate the effect of various conditions on inhaled droplets. The model demonstrates the relevance of temperature and relative humidity of the inhaled air and inhalation flow rate to the behavior of inhaled droplets in the respiratory tract. Furthermore, the model is made available freely for anyone to use and adapt. Due to the step by step approach using discrete equations, we hope to enable users to easily understand and expand the model to further improve the accuracy and perhaps complexity. While a more accurate and complete result might be obtained with a full computational fluid dynamics (CFD) model, the simplicity of this model is easier to understand and adapt by researchers that are not involved in the field of CFD.
The conditions at which the aerosol is inhaled determine the magnitude of the effects on changes in size distribution. Both the relative humidity and the temperature of the inhaled air were found to have a noticeable effect on the fraction of droplets smaller than the inhaled droplets, and thus on the deposition in the respiratory tract. With increasing relative humidity, the deposition will occur at the higher regions of the respiratory tract, while with increasing temperature, the deposition will shift toward lower regions of the respiratory tract. Although the conditions during inhalation will change inevitably with location and season, it is something to keep in mind when developing devices for worldwide applications. It could, for example, be worthwhile to look into air conditioning methods that are built into the device to reduce the effects when the deposition site is critical.
Also the rate at which the aerosol is inhaled is of importance. Especially in the case of the Respimat s soft mist inhaler, which ejects the aerosol in a fixed timeframe, the concentration of the droplets in the air mainly influences the evaporation of droplets in the center of the airway. With more diluted droplets (fast inhalation), it takes longer for the air in the center of the airway to saturate, allowing droplets to evaporate more. This in turn results in a smaller droplet size than when the aerosol is slowly inhaled. Therefore, depending on the device, the instructions should clearly specify which type of inhalation is required. Furthermore, the resistance of the inhaler can also play an important role here. A high resistance inhaler will promote slow inhalation, whereas a low resistance inhaler will allow for a fast inhalation, but not necessarily enforce one. Finally, whether or not the inhalation is performed fast or slow will also affect the flow regime inside the respiratory tract. When the air is inhaled slowly, the flow regime will be closer to laminar, also in the upper airways. However, when the air is inhaled faster, the flow will be more turbulent, by which the behavior of the droplets inside the respiratory tract will change. Also worth mentioning is that the same result of a change in inhalation flow rate can be obtained by an equal inverse change of the inhaled amount of droplets of the same size, because the concentration of droplets in the inhaled air will in that case be the same.
Because the concentration of the generated droplets in the inhaled air can affect the behavior of the droplets in the respiratory tract, the aerosol ejection profile from the inhaler and the inhalation flow profile can also play a role. To minimize variations due to growth and shrinkage of droplets during inhalation, a constant inhalation flow and a square ejection profile, where the maximum throughput is reached instantly and also ends instantly, is preferred. However, in reality there will be a gradual increase and decrease in the ejected volume from the device, the droplets at the beginning and at the end of the ejection will be diluted more than droplets ejected in the middle, when the ejected flow is higher. In part, this could be counteracted by the inhalation flow profile, which will also not be square. At the start of the inhalation procedure the inhalation flow will increase until a maximum is reached and then slowly decrease again. This will also affect the behavior of inhaled droplets in the respiratory tract differently, similar to the inhalation flow rate. These notices could also be taken into account when designing a new generation of Adaptive Aerosol Delivery systems (Dhand, 2010), with an even further improved control over lung deposition.
It has to be noted that there are several limitations to the model that should be kept in mind. First of all, the model proposed in this study assumes a fully laminar flow. While this is true for the largest part, the flow in the first few (two to four) generations may be turbulent (Olson, Sudlow, Horsfield & Filley, 1973). This will have an effect on the behavior of the droplets in the respiratory tract. It can be expected that under turbulent conditions, the difference between the center of the airway and the layers close to the airway wall is smaller. Both the temperature and the relative humidity gradients will be smaller and there will be some mixing of the droplets throughout the airway. This should result in a smaller droplet size distribution, as droplets near the airway wall will grow less and droplets near the center of the airway will evaporate less. In fact, in an ideally mixed system, the droplets would all have the same size. This system is the opposite extreme compared to our fully laminar system. However, when the flow becomes laminar after 2-4 generations, the droplet size distribution will most likely become wider again due to the formation of a temperature and relative humidity gradient. It is possible that, due to increased mass and heat transfer under turbulent conditions from the airway wall to the airway, the average droplet size will increase when turbulent conditions are considered in the first few generations. Therefore, it may be worthwhile to expand the model in the future to also include the various flow regimes that are encountered in the subsequent generations of the respiratory tract.
Furthermore, the heat and moisture transferred to the inhaled air during passage from the mouth to the trachea that was omitted could change the location in the respiratory tract where equilibrium is reached. If fully laminar flow were to be assumed for this region as well then the resulting changes in droplet size would most likely change very little, because the conditions would be similar as in the trachea. More specifically, the droplets near the center will shrink to the same size, because the air will be saturated as this also happens in the trachea (Figs. 3 and 4) and only droplets near the wall might have a bit more time to grow. However, equilibrium will be reached sooner, perhaps after only a few generations depending on the inhalation flow. However, it is known that the flow in this region is slightly turbulent and more accurate results would be obtained by taking this into account when expanding the model (Olson et al., 1973).
Another aspect that is not considered in the model, but might be very important under laminar flow conditions, is the influence of the bifurcations. In our model, the airway simply reduces in size (following the Weibel model) but retains the internal distribution of the droplets. However, under fully laminar conditions, the airway is split in half and the droplets at the center of the airway are suddenly close to the airway wall and exposed to conditions that enhance condensational growth. Therefore, the droplets that are originally located in the center will have less time to evaporate before they start to grow after the bifurcation. This could cause an overall increase in the droplet size due to limited evaporation of droplets in the center compared to what is simulated in our model. However, it was also shown that the evaporation of the droplets near the center of the airway reaches equilibrium due to saturation of the air with moisture (Figs. 3 and 4). This equilibrium is under most conditions reached before the end of generation zero, the trachea, thus before any bifurcations are encountered. Therefore, although the droplets may start to grow sooner (upon entering generation one) the effect on the overall droplet size distribution will be small. Furthermore, the hygroscopicity of solutes is not included in the model. A lot is known about the influence of solutes on the growth on both particles and droplets, and it is possible to model these influences (Ferron, 1977;Golshahi et al., 2013;Longest, McLeskey & Hindle, 2010;Tian et al., 2013Tian et al., , 2011Worth Longest et al., 2012). The model without hygroscopicity can be used for conditions where droplet growth is mainly caused by condensation instead of hygroscopicity, i.e. during inhalation of aerosols where the relative humidity tends to reach over 100%. However, for inhalation of solid particles or aerosols with solutes at lower relative humidity there will be an underestimation of droplet growth due to the sizeable contribution of hygroscopicity to particle growth. Because of the stepwise calculation using discrete equations, the addition of such influence should be quite straightforward and is intended to be done during future development. Therefore, although current results will give a slight underestimation of the droplet size due to the absence of solute hygroscopicity (Fig. 2), future expansions of the model can improve the accuracy. In addition to more accurate results for aerosols, this will also make the model more suitable for particle behavior in the respiratory tract. However, current results indicate that the same effects can be expected for dry particles, since both relative humidity and temperature will dictate how much water will be taken up by dry particles. Therefore, it can be assumed that particles near the airway wall will grow more than particles in the center of the airway when the airflow is laminar.
Finally, it is unknown exactly how much the exclusion of nucleation at an RH 4100% will influence the behavior of inhaled droplets inside the respiratory tract. It is clearly shown that the reduction of the maximum relative humidity due to possible nucleation will definitely affect the droplet growth near the airway wall, as can be expected. However, it is unknown how much water will nucleate and form separate droplets that do not coalesce with the inhaled droplets. An experimental setup, as was used by , could help to tune the model to more exact results. Despite the fact that the actual size distribution of the droplets can shift depending on the effect of nucleation at high relative humidity, the practical implications of the results remain unchanged.

Conclusion
The developed model was used to predict the influence of several parameters on the behavior of inhaled aqueous droplets in the respiratory tract. It was found that, under laminar conditions, the location of the droplet with respect to the airway wall determines whether a droplet will grow or shrink. Droplets in the center of the airway will first evaporate until the air is saturated, while droplets near the airway wall will immediately start to grow due to condensation. The extent of the effects is also shown to be largely influenced by the relative humidity and temperature of the inhaled air, and the inhalation flow rate. Understanding how these factors influence the behavior of inhaled droplets can help during formulation, design of the inhalation device, or inhalation instructions. Furthermore, the model can easily be adapted to include phenomena such as turbulent flow, hygroscopic growth, and others that can help to improve the accuracy of the model. In addition, these additions could also make the model applicable to dry powders.