Aerosol dynamics simulations of the anatomical variability of e-cigarette particle and vapor deposition in a stochastic lung

Electronic cigarette (EC) aerosols are typically composed of a mixture of nicotine, glycerine (VG), propylene glycol (PG), water, acidic stabilizers and a variety of flavors. Inhalation of e-cigarette aerosols is characterized by a continuous modification of particle diameters, concentrations, composition and phase changes, and smoker-specific inhalation conditions, i.e. puffing, mouth-hold and bolus inhalation. The dynamic changes of inhaled e-cigarette droplets in the lungs due to coagulation, conductive heat and diffusive heat/convective vapor transport and particle phase chemistry are described by the Aerosol Dynamics in Containment (ADiC) model. For the simulation of the variability of inhaled particle and vapor deposition, the ADiC model is coupled with the IDEAL Monte Carlo code, which is based on a stochastic, asymmetric airway model of the human lung. We refer to the coupled model as “ IDEAL/ADIC_v1.0 ” . In this study, two different e-cigarettes were compared, one without any acid ( “ no acid ” ) and the other one with an acidic regulator (benzoic acid) to establish an initial pH level of about 7 ( “ lower pH ” ). Corresponding deposition patterns among human airways comprise total and compound-specific number and mass deposition fractions, distinguishing between inhalation and exhalation phases and condensed and vapor phases. Note that the inhaled EC aerosol is significantly modified in the oral cavity prior to inhalation into the lungs. Computed deposition fractions demonstrate that total particle mass is preferentially deposited in the alveolar region of the lung during inhalation. While nicotine deposits prevalently in the condensed phase for the “ lower pH ” case, vapor phase deposition is dominating the “ no acid ” case. The significant statistical fluctuations of the particle and vapor deposition patterns illustrate the inherent anatomical variability of the human lung structure.


Introduction
Due to the chronic exposure to tobacco smoke toxicants, inhalation of mainstream tobacco smoke has been commonly recognized as an important health risk and identified as the main contributor to the occurrence of bronchial carcinomas. Thus the dosimetry of inhaled smoke aerosols in the human respiratory tract is an important component of any toxicological assessment (Asgharian, Price, Yurteri, Dickens, & McAughey, 2014;Pichelstorfer, Hofmann, Winkler-Heil, Yurteri, & McAughey, 2016;Pichelstorfer, Winkler-Heil, & Hofmann, 2013;Robinson & Yu, 2001). In recent years, however, electronic cigarettes (EC) have become increasingly popular because they are considered as a less toxic alternative to conventional smoking due to the absence of combustion aerosols.
Electronic cigarette aerosols are produced by heating of a liquid mixture typically composed of nicotine, vegetable glycerine (VG), propylene glycol (PG), water and a variety of flavors, which are distributed between the condensed and the gas phases. Thus, compared to combustible cigarette aerosols, which contain a multitude of chemical substances (Roemer et al., 2012;St. Charles, McAughey, & Shepperd, 2013), EC aerosol mixtures have significantly fewer constituents (Cheng, 2014), which allows a more realistic description of dynamic changes and deposition of the inhaled aerosol in the human respiratory tract.
Deposition of inhaled EC aerosols throughout the human respiratory tract depends on the physical aerosol characteristics. While vapor is transported very quickly to the system boundaries due to its high diffusivity, particles in the range from about 400 nm to 600 nm can penetrate deeply into the lung and may be deposited there by diffusion, sedimentation or impaction (Hofmann, 2012). Nicotine, a key component associated with addiction, is initially present mainly (>98%; e.g. Lewis, Colbeck, & Mariner, 1995) in the condensed phase. However, simulations have indicated that nicotine is primarily deposited in the lungs in vapor form (Pichelstorfer et al., 2016). Thus, the condensed phase of the aerosol serves as a vector for nicotine during inhalation into the lung. Consequently, a description of nicotine evaporation from the particle phase is required to determine the amount of nicotine deposited in condensed or vapor form. Indeed, if aerosol dynamics are not considered in deposition calculations for e-cigarette droplets, then nicotine mass deposition fractions are significantly underestimated.
A first approach of modeling EC aerosol deposition in the human respiratory tract with the aerosol dynamics model ADiC (Aerosol Dynamics in Containments) was presented in Pichelstorfer et al. (2016). In that paper, the ADIC model was combined with a single-path version of the stochastic lung model and did not yet consider chemical reactions of the particle phase. Model simulations for specific aerosol properties revealed that the contribution of the vapor phase to nicotine mass deposition is significantly higher than that of the condensed phase. Recently, Pithawalla, 2018a, 2018b) reported the application of an aerosol dynamics model to the deposition of e-cigarette constituents in the oral cavity and the human respiratory tract. Model predictions indicated that over 90% of nicotine and propylene glycol were delivered to lung tissue by both aerosol deposition and vapor uptake, while glycerine remained mostly in the condensed phase and hence delivered by particle deposition alone.
In the present study, dynamic changes of inhaled evolving EC aerosol and resulting condensed and gas phase species deposition patterns in the human respiratory tract are simulated by the combination of two models: The aerosol dynamics model ADiC for the simulation of aerosol dynamics mechanisms, and the Monte Carlo deposition model IDEAL (Inhalation, Deposition and Exhalation of Aerosols in the Lungs) for the calculation of particle and vapor transport and deposition in the lungs. The aerosol dynamics model ADiC has been developed in a step-wise fashion from the consideration of coagulation and its effect of inhaled particle deposition in the human respiratory tract (Pichelstorfer et al., 2016;Pichelstorfer et al., 2013) to the inclusion of heat/vapor transport, phase transition (condensation/evaporation) and particle phase chemistry for the analysis of cigarette smoke denuder tube studies (Pichelstorfer & Hofmann, 2015, 2017, and, finally, to the application of the model to particle and vapor phase deposition of nicotine of combustible and electronic cigarette aerosols (Pichelstorfer et al., 2016).
The stochastic lung deposition model IDEAL has been developed for the simulation of particle deposition patterns in the human respiratory tract (Hofmann, 2012;. While the standard version of the IDEAL code refers to normal breathing conditions with continuous inspiration and expiration phases, inhalation of the cigarette aerosols is characterized by three stages: puffing, mouth-hold, and bolus inhalation/exhalation (Pichelstorfer & Hofmann, 2015;Pichelstorfer et al., 2013).
In the case of inhalation of inert particles commonly used in human inhalation studies, inhaled particle sizes remain constant during a complete breathing cycle, and deposition is the only mechanism to reduce the initial particle mass and number concentrations. For comparison, inhaled EC droplets experience -besides deposition -coagulation, condensation of water vapor, evaporation of volatile compounds, and chemical reactions which change particle diameters, number concentrations and particle properties (Pichelstorfer et al., 2016). Moreover, in addition to particle transport in the lungs, it is necessary to simulate the transport and deposition of the vapor as nicotine is deposited in the lungs in both particulate and vapor form (Pichelstorfer et al., 2016).
Thus the objectives of the present study are (i) to combine the ADiC model with the IDEAL code to fully explore the effects of the stochastic airway structure on the dynamic changes of inhaled EC droplets, (ii) to compute resulting number and mass deposition patterns of particle and vapor phase constituents, particularly for nicotine, and (iii) to assess the statistical fluctuations of particle and vapor deposition patterns as a result of the inherent anatomical variability of the human lung structure.

Materials and methods
For the simulation of EC aerosol deposition in the human respiratory tract, three phases have to be distinguished, partly operating in different anatomical regions: 1) puffing in the mouth, 2) mouth-hold in the oral cavity, and 3) inhalation/exhalation in the throat and downstream airways of the lungs. During puffing, the mouth cavity expands and (e-) cigarette aerosol is continuously filling the mouth, where they experience coagulation, phase transition, heat and vapor transport, dilution and mixing, chemical reactions and turbulent particle and vapor deposition. During the subsequent mouth-hold phase, operating mechanisms in the oral cavity are coagulation, phase transition, heat and vapor transport, chemical reactions and diffusional and gravitational deposition. During inhalation and exhalation in the throat and the lungs, coagulation, phase transition, heat and vapor transport, chemical reactions and diffusional, inertial and gravitational deposition are considered. Inhalation of the EC droplets into the lungs is characterized by an aerosol bolus, representing the volume of the oral cavity, followed by the inhalation of aerosol free air during the remaining time of the inhalation phase. These aerosol dynamics mechanisms will be explained in more detail in the framework of the ADiC aerosol dynamics model and the IDEAL deposition model as described in subsequent text.

Aerosol dynamics mechanisms
The dynamic changes of inhaled cigarette smoke particles during puffing, mouth-hold and within the lungs and during inspiration and expiration are described by the aerosol dynamics model ADiC, which considers coagulation, conductive heat transport, diffusive vapor transport, phase transition (i.e. condensation and evaporation of volatiles), dilution and mixing, and chemical reactions.
The dynamic changes of the number concentration of particles with size i, n i , as a function of time t, ∂n i /∂t, in a given morphometric respiratory element, such as mouth, throat, airways and alveoli of the human lung can be described by the discrete general dynamics equation (GDE) (Hidy, 1984): where q is the convection related velocity, D p is the binary diffusion coefficient of the particle in the carrier gas, and q F is the external force related velocity. The first term on the right side of the equation represents convection related changes in number concentration, the second term diffusion, the third, fourth and fifth terms phase transition, coagulation and nucleation, respectively, and the last term the influence of external forces, such as physical deposition mechanisms. In the present case of high particle concentrations, the nucleation term can safely be neglected as particles, potentially formed by nucleation, quickly diffuse onto the surface of larger particles. The different effects of the dynamic changes caused by coagulation, phase transition, heat/vapor transport, dilution/mixing and chemical reactions have already been described in more detail in previous publications on combustion and EC aerosol deposition in the human respiratory tract (Pichelstorfer et al., 2016;Pichelstorfer et al., 2013) and in denuder tubes (Pichelstorfer & Hofmann, 2015, 2017, as well as for the hygroscopic growth of NaCl particles (Winkler-Heil, . However, for the sake of the reader, the salient features of these dynamic processes and their computational realizations are briefly described below.

Coagulation
For sufficiently high particle number concentrations, as in the case of inhaled EC droplets (of the order of 10 9 particles per cm 3 ), coagulation quickly reduces number concentration, increases particle diameters and changes the composition of the particles (Pichelstorfer et al., 2013). The dynamic changes of the number concentration of particles with size i, n k , as a function of time t, dn i /dt, is described by an equation proposed by Smoluchowski (1917): where β ik is the coagulation coefficient of colliding particles i and k, which depends on the coagulation mechanism considered, and s is the number of different particle sizes. In the case of spherical liquid particles, two droplets may collide due to their relative motion, forming a new spherical particle with a larger diameter. Relative motion between two particles can occur by thermal motion, gravitational settling, laminar shear, turbulent flow, electrical charge and inertial effects at airway bifurcations (Pichelstorfer et al., 2013). For the simulation of EC droplets, only coagulation driven by thermal motion has been considered, as this mechanism is dominating coagulation for the given situation. Since the size distribution of e-cigarette particles falls into the transition regime between the continuum (small Knudsen numbers) and the kinetic regime (high Knudsen numbers), the approximation derived by Fuchs (1964) for the thermal coagulation coefficient was applied. For a detailed description of the different mechanisms causing coagulation and of the mathematical simulation of the dynamic changes of the initial particle size distribution (PSD), the reader is referred to the paper of Pichelstorfer et al. (2013).

Heat and vapor transport
Heat and vapor transport affect the temperature, the vapor concentration and, indirectly, the properties of the condensed phase in an airway volume. Due to differences in temperature, heat can be transported from the airway walls to the air volume if the temperature of the inhaled air is lower than the body temperature, such as inhalation of normal outside air. However, heat can also be transferred from the air volume to the walls if the temperature of the inhaled air is higher than the body temperature, as in the case of inhalation of heated EC droplets. Likewise, airway walls can also be a sink or source of vapor concentration. Due to differences in relative humidity between the inhaled air and the very humid environment of the respiratory tract, water vapor is transported from the airway walls to the inhaled air volume. However, in the case of volatile species produced by evaporation from the EC droplets, vapor transport from the air volume to the walls depends on the concentration gradient between airborne and surface vapor concentration.
The physical mechanisms determining heat and vapor transport are, besides convective transport, conductive heat transfer and diffusive vapor transport. The transport of heat and vapor across system boundaries and heat and vapor exchange in the case of laminar flow through a circular tube can be described by the equation of Jacob (1950) for heat conduction as well as vapor transport by diffusion: where u either stands for the temperature or the vapor concentration, α = κ/(c p ρ) (where ρ is the gas density, κ is the heat conductivity and c p is the specific heat capacity of the gas) in case of heat transport and α = D (where D is the diffusion coefficient of the vapor species in the carrier gas) in case of diffusional vapor transport, respectively. In addition to the diffusive vapor transport described by equ.
(3), vapor transport in cylindrical airways is further determined by convective transport and axial diffusion (Scherer, Shendalman, Greene, & Bouhuys, 1975), both depending on the flow velocity in a given airway. A detailed description of the mathematical solution of equ.
(3) for heat and vapor transport in cylindrical airways can be found in the paper of Pichelstorfer and Hofmann (2015).

Phase transition
Condensation of water molecules on inhaled particles in the humid atmosphere of the oral cavity and subsequent lung airways as well as condensation and evaporation of other volatile compounds on or from inhaled particles are caused by differences in partial vapor pressures. C Barrett and F Clement (1988) give the steady-state mass flux for a spherical geometry as: where d p is the particle diameter, p is the total pressure, M v is the molar mass, D is the binary diffusion coefficient of the condensing species in the carrier gas, R is the general gas constant, T ∞ , p v,∞ and p v,dp are the ambient and the vapor pressure at the particle surface, respectively. The particle-vapor phase transition affects the diameter of the particle, either increasing diameters by condensation or decreasing diameters by evaporation of volatile compounds, as well as the composition of the particles. In the present model, the interaction between particle and vapor phase is based on the formulations proposed by Vesala, Kulmala, Rudolf, Vrtala, and Wagner (1997) and Mattila, Kulmala, and Vesala (1996), i.e. assuming stationary vapor and temperature profiles and a quasi-stationary growth and shrinkage processes.
The phase/transition model includes the effect of mixing of volatiles and solubles within the particle on the vapor pressure (Raoult's law) and of the curvature of the particle surface on the vapor pressure (Kelvin correction), along with the effect of nonideality of nicotine-water mixture (including protonation of nicotine in the condensed phase). Corrections for the transition regime are taken into account. Note that mass transfer by phase transition is always accompanied by a flux of latent heat either heating (condensation) or cooling (evaporation) the condensed phase which hampers the process.
For a detailed description of the mathematical formulations describing phase transition, the reader is referred to the paper of Pichelstorfer, Winkler-Heil, & Hofmann (2015).

Equilibrium chemistry
The particle phase chemistry sub-model, which is strongly linked to the phase transition model by affecting the properties of the volatile species, considers acid-base reactions. Note that we assume the ions to be non-volatile (Pichelstorfer & Hofmann, 2015). Reaction rates are neglected since we assume equilibrium chemistry, i.e. forward and reverse reactions of reactants and products are equal. The acid donates a proton and is thereby converted to its conjugate base and a base is converted to its conjugate acid by accepting a proton: where [HA] is the molar concentration of the acid HA and H + , A -, B, BH + are hydronium ion, conjugate form of acid, base and conjugate form of base, respectively. These assumptions allow the construction of a polynomial of H + which can be solved numerically. For a detailed description of the mathematical formulations describing chemical reactions and how they are solved, the reader is referred to the paper of Pichelstorfer and Hofmann (2017).

Dilution, mixing
During puffing, a constant and continuous (time resolution higher than 10 -4 s) flow of inhaled air feeds fresh EC particles to the continuously aging aerosol in the oral cavity. Thus mixing between the residual air in the mouth and the inhaled air volume during puffing and subsequent mouth-hold leads to a mixing of the previously inhaled particle concentration and the incoming fresh aerosol. In addition to a change of the temperature in the oral cavity, mixing affects both, the condensed as well as the gas phase.
Note that instantaneous mixing is assumed, i.e. the incoming aerosol is immediately mixed with the aerosol already present in the mouth. During the mixing process, the air volume in the mouth continuously increases with time as the oral cavity gradually expands.

Particle and vapor deposition
Deposition of inhaled particles in the human respiratory tract reduces not only their number concentration but also affects their size-distribution because of the size-sensitivity of operating physical deposition mechanisms. Likewise, exchange of vapors between the gas phase and the airway walls alters the concentration of the vapors and may thus trigger phase transition. Compared to normal inhalation with constant inhalation and exhalation phases, inhalation of cigarette and EC aerosols is characterized by three distinct inhalation phases, puffing, mouth-hold and bolus inhalation, followed by a constant exhalation phase. As a result of these phases, particle deposition models distinguish between deposition in the mouth during puffing and mouth-hold, and throat and lung deposition of an aerosol bolus. For comparison, deposition of the exhaled bolus is based on a constant flow rate in all regions of the respiratory tract.

Mouth deposition model
Deposition of inhaled particles in the mouth consists of two phases. The first phase is the puffing phase, where fresh aerosol is continuously added to the gradually aging aerosol in the mouth and mixing with the aerosol already there. In the second phase, the aerosol ages in the oral cavity and deposition is caused by diffusion and sedimentation in stagnant air ).

Puffing phase -turbulent mixing deposition.
During the puffing phase, the aerosol is drawn from the e-cigarette and enters the gradually expanding mouth cavity. As a result the aerosol is continuously mixed with fresh EC-aerosol. The optimal solution for the description of air flow in the mouth during puffing would require the application of computational fluid dynamics (CFD) methods. However, there are three main difficulties, which makes CFD simulations difficult to apply, namely (i) the expanding mouth cavity, and (ii) the fact that there is only an inflow but no outflow. Further, the authors are not aware of any available CFD package representing all relevant aerosol dynamics processes, which strongly impact on the aerosol at this stage.
Thus, and in order to provide a computational less extensive and generally applicable approximation, a model introduced by (Lai and Nazaroff, 2000) describing the deposition in a turbulently mixed chamber is used. The model is based on the assumption of a completely mixed volume separated from a surface by a boundary layer. It considers near surface turbulences to describe the effects of gravitational settling, Brownian and turbulent diffusion on particle deposition. Note that the expressions for diffusional mass transport may also be applied to describe conductive heat and vapor transport. Lai and Nazaroff (2000) provide solutions for containments (i) composed of horizontal and vertical walls, and (ii) spherical containments. In the present work we applied solution (i). See Appendix A for details.
Recent experimental (Montigaud et al., 2021;Zhou et al., 2020) and computational (Ashgarian et al., 2018a;Haghnegahdar, Feng, Chen, & Lin, 2018;Feng, Kleinstreuer, Castro, & Rostami, 2016) works described the regional deposition to be in the range from 2% to 10% in the oral cavity. These values are dependent on the conditions (e.g.: puffing volume, puffing time, composition of the e-liquid). In the present work, we assume the regional mass deposition in the extra-thoracic region to be roughly 10%. The changes to the PSD by deposition in the oral cavity are calculated applying the chamber deposition formulation described above scaled to the total loss of roughly 10%.
Note that vapor deposition utilizes the same formalism but without parameterization.

Mouth-hold phasestagnant air deposition.
Particle deposition in the oral cavity during mouth-hold was computed for diffusion and gravitational settling, with diffusion being the dominant deposition mode (Pichelstorfer et al., 2013). Radial diffusive transport was simulated in a volume-equivalent sphere based on the mean square displacement of the particles. Deposition by sedimentation was calculated for a volume-equivalent cuboid, applying the particle settling velocity given by Hinds (1995).

Throat deposition model
Since practically all deposition events in the oropharyngeal region due to inertial impaction are caused by the 90 o bend leading to the pharyngeal and laryngeal regions, but not in the oral cavity, the semi-empirical deposition equation for the whole oropharyngeal region derived by Cheng (2003) was applied for deposition in the throat. Since diffusional deposition in the mouth is already considered by the mouth deposition model described above, deposition of submicron EC droplets and related vapor species in pharynx and larynx, which are approximated by an 8 cm long and 2 cm diameter cylindrical channel (Cheng, Zhou, & Chen, 1999;ICRP, 1994;Stapelton, Guentsch, Hoskinson, & Finlay, 2000), is simulated by the equations proposed by Ingham (1984) for laminar developing flow.
Upon exhalation, the semi-empirical equations for deposition by diffusion and impaction in the oropharyngeal region, i.e. throat plus mouth, are applied (Cheng, 2003).

Lung deposition model
Particle deposition in the airways of the human lung is computed with the IDEAL Monte Carlo code (Hofmann, 2011;Hofmann, Winkler-Heil, & Balásházy, 2006;Winkler-Heil, Ferron, & Hofmann, 2014), which simulates the random walk of individual inhaled and exhaled particles through an asymmetric, stochastic airway model of the human lung . As is current practice, the location of individual airways within the lung relative to the trachea is characterized by an assigned airway generation number, typically starting with generation 0 (ICRP, 1994) for the trachea (1994) and increasing with penetration into the lung. Consequently, the branching network of the whole human lung is represented by a serial succession of cylindrically-shaped airway generations (ICRP, 1994) or by a sequence of Y-shaped airway bifurcations (Koblinger & Hofmann, 1985).
The stochastic airway geometry of the lung is based on statistical analyses of the morphometric data of bronchial (Raabe, Yeh, Schum, & Phalen, 1976) and acinar (Haefeli-Bleuer & Weibel, 1988) airways (Koblinger & Hofmann, 1985). For the deposition calculations, the measured bronchial airway dimensions (Raabe et al., 1976) obtained at total lung capacity (TLC) were scaled down to a functional residual capacity (FRC) of 3300 cm 3 (ICRP, 1994) by assuming a constant scaling factor. However, no scaling procedure was applied to the acinar data of Haefeli-Bleuer and Weibel (1988) as they already represent FRC dimensions. With the exception of trachea, main and lobar bronchi, which are assumed to be rigid structures, an additional scaling procedure was applied to diameters and lengths to account for the linearly increasing and decreasing lung volume during inspiration and expiration of tidal volume V T , approximated by an average constant lung volume of FRC + V T /2.
In the standard version of the IDEAL code for continuous breathing, particles are inhaled at random times during the inspiration phase, which leads to a distribution of path lengths within the lungs, ranging from the most distal airways, for particles inhaled at the very beginning of the inhalation phase, to only the trachea, for particles inhaled towards the end of the inspiration phase. In the present study, however, the inhalation phase is restricted to the length of the aerosol bolus given by the volume of the oral cavity and the flow velocity. Furthermore, asymmetric flow splitting at asymmetrically branching airway bifurcations, based on distal lung volumes, leads to different paths of inspired particles. Both random inhalation times and random paths of inhaled particles lead to highly variable deposition fractions in individual airways. In the present version of the IDEAL code, particle paths are terminated by a pre-specified number of simulations. The equations determining deposition efficiency applied in the present study, are listed in Appendix B. This stochastic particle deposition model has been validated by comparison with the experimental data of Heyder, Gebhart, Rudolf, Schiller, and Stahlhofen (1986) for small and larger particles and of Schiller, Gebhart, Heyder, Rudolf, and Stahlhofen (1988) for ultrafine particles, and with the fitted deposition data provided by the ICRP (1994) report.
In contrast to particle deposition, which is determined by the IDEAL code, vapor deposition is calculated by the ADIC code. Vapor deposition in individual bronchial airways or airway segments was considered by calculating the convective and diffusive transport, along and perpendicular to the main axis, of vapor molecules in a laminar flow field in cylindrical tubes. A detailed description of the equations applied can be found in Pichelstorfer and Hofmann (2017).

Implementation of the ADiC model into the IDEAL model
The initial version of the IDEAL deposition model , addresses the fate of single, i.e. non-interacting, particles, randomly selected from the inhaled particle size distribution (PSD), on their way through a full breathing cycle until they are either deposited or exhaled (Lagrangian formulation). Thus, except for the loss of individual particles by deposition, the inhaled the inhaled PSD does not change its shape. However, in the case of inhaled EC droplets, condensation or evaporation can modify the diameters of the inhaled EC droplets, depending on the surrounding vapor concentration. Furthermore, due to high particle concentrations, the diameter of a particle can also increase by coagulation, which is as a function of the aerosol concentration. As a result of these dynamic processes, the initial size distribution prior to the onset of inhalation is no longer constant and the calculation of the correct diameter in a given site requires information on particle and vapor concentrations at that site. Thus instead of tracking individual particles, the random path of an elemental air volume (bolus) containing information on particle size distribution and concentration as well as on vapor concentration is simulated, thereby adding an Eulerian element to the Lagrangian random path model (Pichelstorfer et al., 2013).
In this combined Lagrangian/Eulerian deposition model, the initial size distribution, vapor composition and size-dependent particle properties are continuously modified by coagulation, condensation, evaporation, chemical changes and deposition. Because of the Lagrangian nature of the deposition model, these simultaneously occurring processes are treated independently in a consecutive fashion. To minimize the effect of considering these mechanisms consecutively rather than simultaneously, airways are divided into a number of tubular segments or time intervals. At present, the length of a given individual time interval is determined by a maximum Fig. 1. Schematic description of the task-distribution in IDEAL/ADiC_v1.0 for simulating lung deposition. On the left hand side, a random (indicated by probability p(x)) lung structure is generated by IDEAL based on morphometric data and its variation. The aerosol, traveling along a random path (orange) through this structure, undergoes dynamic changes, described by ADiC, as depicted by the flow chart on the right hand side. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.) change of the characteristic parameters of a process (i.e. particle number concentration or vapor phase concentration of a condensable substance) by 0.1% in all particle diameter intervals. More details on these characteristic or control parameters can be found in Pichelstorfer and Hofmann (2015). Particle deposition probabilities are then calculated for a sequence of airway segments until the particle size distribution has passed the whole airway bifurcation.
While the initial version of the IDEAL deposition code has been used for the simulation of the transport and deposition of inert particles for a variety of applications, the current version of the IDEAL deposition model for the inhalation of dynamic aerosols, extended by the implementation of the ADiC model, also considers the transport and deposition of different vapor species.
A graphical description of the functionality of the IDEAL/ADiC_v1.0 and which tasks are associated with former IDEAL and ADiC is provided by Fig. 1. A stochastically generated lung structure is provided by IDEAL based on morphometric data and its distribution. The random trajectory along which the aerosol travels through this structure is also provided by IDEAL. This trajectory is divided into intervals which are characterised by mean residence time (t res ) and a description of the boundary conditions (B). Changes to the aerosol during t res are then calculated by ADiC considering dynamic processes as described in section 2.1. The processes are calculated consecutively. However, a control parameter "C" takes care that changes by each process do not exceed a set limit by adjusting the integration time step for aerosol dynamics (Δt) and thus allowing for quasi-simultaneous solution of all processes. More information on this method can be found in Pichelstorfer and Hofmann (2015).
The whole process of generating a stochastic lung structure, determining a random path and calculating the dynamic changes along the path is repeated N times (N is typically between 10 3 -10 4 ). Results are averaged to bridge from considering a single path to considering the whole lung. Note that the deposition fraction of particles of any size are determined by IDEAL for the respective time interval prior to calculating the dynamic changes. Further, the description of dynamics and deposition is fully described by ADiC for the ET region (see section 2.2).
The combined ADiC/IDEAL has been partly verified for a specific combination of several mechanisms by comparison with experimental data on total deposition of cigarette smoke particles, i.e. considering the effect of heat and vapor transport, phase transition and coagulation (Pichelstorfer et al., 2013) and on hygroscopic growth and deposition of NaCl particles, i.e. simulating heat and vapor transport and phase transition (Winkler-Heil et al., 2017). Thus, the excellent agreement with these experimental data represents a reliable baseline for the simulation of the deposition of inhaled EC droplets. However, note that these earlier combinations of ADiC/IDEAL were designed to simulate very specific cases only: in the work by Pichelstorfer et al. (2013), only the coagulation sub-model was coupled with the IDEAL code to consider the effect of coagulation on regional aerosol deposition in the human lung. In Pichelstorfer et al. (2016), the IDEAL code was used solely to generate a single path representation of the human lung, including geometry, deposition probability and mean residence time of the aerosol in the different regions of the lung. In the work Winkler-Heil et al. (2017), the IDEAL code was coupled with a specific (for water condensation on NaCl particles) version of the ADiC phase transition sub-model. Thus, we consider the present work the first with fully coupled IDEAL and ADiC models; and refer to it as "IDEAL/ADiC_v1.0". The latter allows to exploit full functionality of both individual models.

Computational methods
In the present study, the ADiC model was implemented into the IDEAL code. The IDEAL code serves to generate the lung geometry, i.e. the boundary conditions of the system, and determine the deposition of particulate matter. The ADiC model is then used to calculate aerosol dynamics, namely coagulation, phase transition, particle phase chemistry, conductive/convective heat transport, diffusive/convective vapor transport, the interaction of the vapor phase with the boundaries of the system and the deposition of the particle phase in the extrathoracic (ET) region. This section outlines the procedures how the models are merged to solve the equations for transport, aerosol dynamics and deposition of particles and vapors for the ET region and the stochastically generated lung geometry.

The ET region
In the ET region, the ADiC model is exclusively applied to describe aerosol dynamics as well as interaction of the aerosol with the system boundary during three stages, i.e. puffing, mouth-hold and inhalation of the aerosol bolus through pharynx and larynx. Phase transition, particle phase chemistry and coagulation are calculated in a similar manner for all three stages as they can be described independently from the boundaries of the system. Particle deposition, heat and vapor transport are described by solving different sets of equations for each stage in order to account for different boundary conditions (see section on particle and vapor deposition).
During the puffing phase, the mouth expands and the aerosol is drawn into the mouth cavity at a certain rate which may be time dependent. Here we apply a model describing deposition of particles, transport of vapor molecules and heat transport for a turbulently mixed chamber. Note that a scaling parameter is applied to fit the calculated particle mass deposition to experimental data. Heat transport and diffusional particle and vapor transport along a gradient between the center of the mouth towards the walls is applied during the mouth hold period. During inhalation, deposition of particles is calculated by applying empirical deposition formulas. Heat and vapor exchange is calculated by assuming a simplified cylindrical geometry in the inhalation phase. In a previous work it could be shown that no significant error in temperature and water humidity arises by this assumption (Winkler-Heil et al., 2017).

The lung
The simulation of the lung starts with the stochastic generation of a path through the lung by the IDEAL code. This process is described in detail in section 2.2.3. The inhaled aerosol travels along that path during inhalation and exhalation. Each lung generation is divided into several, individual sections, typically consisting of one half of the airway lengths, or even less (Winkler-Heil et al., 2017). First, the IDEAL code calculates the deposition probabilities for particles of all sizes within a section. Then, the particle concentration is reduced accordingly before the ADiC model is applied to determine the change in vapor phase composition (by phase transition and vapor transport from/towards the walls), gas temperature (latent heat from phase transition and heat transport form/towards the wall) and particle size distribution (phase transition, particle phase chemistry, coagulation). Since the effects of changes in aerosol composition, calculated by ADiC, on particle deposition, calculated by IDEAL, are quite small within a given section, these processes are solved quasi-simultaneously.

Solving the equations
The equations are solved in a box (note, the present aerosol dynamics model is a zero dimensional or box model), traveling along a trajectory generated by the IDEAL code, which is also used to calculate particle deposition. It further provides boundary conditions of the respective section where the aerosol is presently residing (wall temperature, wall saturation ratio for vapors, length and diameter section which is assumed to be of cylindrical shape) and the residence time within the section. The aerosol is assumed to be uniformly distributed across the geometry cross section which leaves time as the only unknown. The time changes of vapor and particle phase caused by aerosol dynamics are determined by solving the respective differential equations by means of Eulerian forward integration. The processes are calculated in a consecutive fashion. However, as the integration time steps are automatically adopted to assure the desired accuracy (typically changes to a quantity are limited to 0.1% during an integration time step), the processes are calculated quasi-simultaneously. A detailed description of how the equations are solved can be found elsewhere (Pichelstorfer & Hofmann, 2015).

Description of the particle size distribution
As the aerosols considered in this work experience considerable changes due to coagulation as well as condensation, the particle size distribution can neither be represented by a sectional model nor by a moving average diameter model. The first model is appropriate to describe processes such as coagulation, nucleation or deposition (Gelbard & Seinfeld, 1980), while showing considerable errors, namely numerical diffusion, when describing condensation (Gelbard, 1990). Moving average models are well suited to describe condensational changes, as they allow for continuous particle size changes, but produce a huge number of new particle sizes during coagulation events. Thus we apply a moving distribution sectional structure which features sections containing normalized particle size distributions (Pichelstorfer & Hofmann, 2015). It allows for the simultaneous calculation of coagulation and growth, while keeping the number of sections constant and avoiding numerical diffusion.

Model inputs
The determination of inputs for computer simulations that allow for a generally representative assessment of local effects is not trivial for mainly three reasons: Firstly, there exist plenty of e-liquid brands and formulations although they are mostly composed of the same substances, which are glycerol, propylene glycol, water, nicotine and some flavor additives (Geiss, Bianchi, Barahona, & Barrero-Moreno, 2015;Oldham, Zhang, Rusyniak, Kane, & Gardner, 2018;Sundahl, Berg, & Svensson, 2017). The composition of the aerosol directly affects its potential to grow by condensation or shrink by evaporation. As a consequence, condensed and vapor phase deposition are indirectly affected by the composition of the e-liquid.
Secondly, a recent study (Robinson, Hensel, Morabito, & Roundtree, 2015) conducted under natural conditions reports a significant inter-subject variability regarding the puffing parameters, i.e. puffing topography, volume and time, respectively. Mean puffing ranged from 0.9 to 6.9 s with relatively low standard deviation for each subject. Puff volumes in the range from 29 to 388 mL were reported, indicating that "taking a puff" in some cases means "taking a breath" as the puffed volume often exceeded the volume of the oral cavity which is about 71.2 ± 15 mL for men and 55.4 ± 13.4 mL for women (Nascimento, Cassiani, & Dantas, 2012). Thus, the typical smoking scheme (puffingmouth holdinhalationbreath holdexhalation) as reported in the literature (Guillerm & Radziszewski, 1978) is not applicable to all consumers of e-cigarettes. Further, volume flow rates during puffing between 24 and 102 mL s -1 were observed. The fact that characteristic parameters of the aerosol (particle number concentration and median diameter, respectively) depend on the puffing flow rate and are different for each device further complicates the characterization of e-cigarette aerosol (Ingebrethsen, Cole, & Alderman, 2012).
Finally, the determination of the PSD generated by electronic devices is extremely challenging as the individual aerosol droplets are on the one hand fully volatile (i.e. dilution would lead to evaporation) and on the other hand appear in high concentrations (i.e. no dilution results in considerable coagulation) when they enter the oral cavity (Ingebrethsen et al., 2012;Oldham et al., 2018;Pankow, 2017).
There are many studies focusing on the determination of the PSD generated from electronic cigarettes (Belka, Lizal, Jedelsky, Jicha, & Pospisil, 2017;Fuoco, Buonanno, Stabile, & Vigo, 2014;Ingebrethsen et al., 2012;Pratte, Cosandey, & Goujon-Ginglinger, 2016). As input for this modelling study we focused on the ones applying impactors (Oldham et al., 2018;Steven, Chen, Serban, & Stephen, 2015;Sundahl et al., 2017) as this method does not require dilution or any other manipulation of the aerosol. There, reported mass median diameters range from about 0.5 μm to 1.1 μm and the geometric standard deviations are in the range from 1.5 to 2.1. Oldham et al. (2018) considered the difference between mass loss in the cartridge and the mass gathered on the impactor stages and the filter pad. They found 74%-91% mass being detected. This quantification is of big importance regarding the fact that the aerosol is fully volatile. This is why we decided to use the aerosol size parameters they determined (geometric standard deviation σ g = 1.95 and mass median aerodynamic diameter MMAD = 1.02 μm, which corresponds to a median count number diameter of about 0.27 μm, assuming a log-normal distribution of the aerosol particle diameter).
Further, we chose an average composition of the e-liquids investigated by Oldham et al. (2018), i.e. 50.6% VG; 35.5% PG; 11.0% water and 2.9% nicotine. Tayyarah and Long (2014) found that the aerosol composition is similar to the bulk liquid composition. Thus we assume that the aerosol has the same composition as the bulk liquid. The initial composition of condensables in the gas phase is set to be in thermodynamic equilibrium with the particles of mass median diameter of the aerosol. Unfortunately, none of the above-mentioned studies applying impactors provide information on particle number concentrations. Thus, in order to choose an adequate aerosol mass loading, we considered the study by Gillman, Kistler, Stewart, and Paolantonio (2016) on the effect of power level on the mass yield of e-cigarette devices. They evaluated five different devices with power levels from 5 Watt to 25 Watt, which produced mass loadings from 1.5 to 27 mg per puff. Most loadings are found in the range of 2 mg-10 mg per puff. Based on this data we chose a moderate mass loading of 5.5 mg per puff as simulation input. The mass loading combined with the shape of the particle size distribution allows for determination of the particle concentration: 1.6 x 10 9 cm -3 . The puff volume was set to 55 mL per puff, as this value was used by Oldham et al. (2018) when they determined the particle size distribution.
Related breathing conditions were: 3.5 s puff flow, i.e. continuous inhalation of 55 cm 3 into the oral cavity (Robinson et al., 2015), 1 s mouth-hold, and 2.5 s inhalation into the lungs with a tidal volume of 750 cm 3 , consisting of the inhalation of the 55 cm 3 bolus from the mouth during the first 0.167 s, followed by aerosol-free air in the remaining 2.366 s, and 2.5 s exhalation.

Limitations
The primary objective of the present calculations is to investigate the effect of the anatomical variability of the human lung on the resulting deposition patterns of e-cigarette droplets. The stochastic whole-lung model IDEAL is based on the morphometric measurements of the bronchial tree of two adult lungs (Raabe et al., 1976) and corresponding measurements of the alveolar structure of several selected acini (Haefeli-Bleuer & Weibel, 1988). This anatomical variability described by the stochastic lung model comprises lung volume (FRC) as well as airway dimensions, branching and gravity angles. Although other published anatomical whole lung models refer to different lung structures (e.g. Weibel, 1963;Horsfield & Cumming, 1968;ICRP, 1994), they provide only average values for the corresponding airway dimensions. Thus the stochastic lung model is currently the only model which provides information on the stochasticity and asymmetry of the lung structure, which is the very focus of this study. Although the anatomical variability of the extrathoracic region, i.e. size and shape of the oral cavity and pharynx-larynx dimensions differ among several subjects (Cheng, 2003;Feng et al., 2018), typical values from the published literature were adopted for the present simulations, thereby underestimating the effect of extrathoracic anatomical variability on deposition.
The physiological variability of puffing volume, puffing time, mouth-hold time, tidal volume and breathing cycle times, i.e. inhalation, breath-hold and exhalation times, and constant or cyclic flow rate further increases the variability of the predicted particle and vapor deposition patterns. For example, significant inter-subject variations of puffing parameters and smoking schemes (Robinson et al., 2015), as well as individual variations of tidal volume and breathing frequency (ICRP, 1994) have been reported in the literature. Instead, typical values from the published literature have been selected for the calculations. Thus the variability of the breathing parameters is not captured in the present study, again underestimating the variability of the predicted deposition patterns. Thus, in summary, the limitations of the model in terms of anatomical and physiological variabilities may lead to a higher variability of the calculated particle and vapor deposition patterns among bronchial and acinar airway generations than predicted in this study.
Particle deposition in the IDEAL deposition model is determined by analytical deposition equations for specified idealized flow profiles in straight cylindrical tubes and spherical alveoli. While the transport of individual particles through the airway system is characterized by a one-dimensional path, related analytical deposition equations describe the loss of particles by diffusion, impaction and sedimentation in all three dimensions. For comparison, Computational Fluid and Particle Dynamics (CFPD) models allow the simulation of random paths and deposition of individual particles in three dimensions under realistic flow conditions. For example, CFPD simulations have been published on flow and particle deposition in simple mouth-throat models (Chen, Feng, Zhong, & Kleinstreuer, 2017;Grgic, Finlay, Burnell, & Heenan, 2004) or in a realistic cast of the human upper respiratory tract (Nordlund et al., 2017). For deposition in the bronchial tree, Zhang, Kleinstreuer, & Kim (2009) presented a complete bronchial airway model consisting of a sequence of parallel adjustable triple bifurcation units based on Weibel's (1963) lung model. Computed generational deposition fractions agreed very well with predictions obtained with analytical deposition equations, although considerable differences were observed for local deposition patterns. A more realistic bronchial airway model was developed by Tian et al. (2011) andLongest, Tian, Delvadia, andHindle (2012) based on the Yeh and Schum (1980) model, considering also branching asymmetry and gravity angles in different lobes, coupled to a mouth-throat model. Finally, Kolanjiyil and Kleinstreuer (2017) developed a whole-lung airway model, consisting of a subject-specific upper airways from nose/mouth to about airway generation 3, which are connected to adjustable symmetric triple bifurcation units. These bifurcation units are in series and parallel in one plane throughout the whole bronchial tree and the pulmonary region, with expanding and contracting alveoli attached to the pulmonary airway generations. Predicted results of total and regional deposition fractions for 3 μm particles reveal agreement within ±10% with the experimental data of Heyder et al. (1986). It may be noted here, that simulations with the IDEAL code revealed also excellent agreement with these experimental data (Hofmann, 2011).
For defined inhalation conditions, deposition in the human lung depends on the geometry of the lung and the physics of flow and particle deposition. While CFPD models allow realistic flow and particle transport conditions, they are based on a simplified description of the lung structure, e.g. all airways lie in a plane. In contrast, one-dimensional whole lung models, such as the IDEAL model, rely on simplified descriptions of airflow and particle deposition, which is compensated, however, by a realistic airway geometry. Although results of both the IDEAL model (Hofmann, 2011) and the CFPD whole lung model (Kolanjiyil & Kleinstreuer, 2017) agree with the experimental results of Heyder et al. (1986), it must be noted that only the stochastic lung model allows to study the effect of anatomic airway variability on deposition, which is the very objective of the present study.
Due to mixing in the alveoli, differences between inhaled and exhaled flow profiles and axial diffusion, a small fraction of inhaled droplets will still remain airborne in the lung at the end of the exhalation phase. Thus in the case of multiple breaths, these particles will be drawn deeper into the lung during the next inhalation phase until they are either deposited or exhaled. The present model IDEAL/ADiC_v1.0 is capable of describing the effect of multiple breaths. However, in the present work we considered a single breathing cycle and thus neglected the effect of multiple breaths.

Results
In this study, two different e-cigarettes were compared: (1) an e-cigarette formulation without any acid added ("no acid"). In this case the pH is dominated by nicotine resulting in a basic solution and a very low degree of nicotine protonation, and (2) a formulation which adds an acidic regulator (about 1% benzoic acid of e-liquid mass) to establish a pH level of roughly 7 for the bulk substance. This case is referred to as "lower pH".
Subsequent results illustrate the evolution of the aerosol during a complete breathing cycle and the deposition of particulate and vapor phase substances in the human respiratory system, with special focus on their variability caused by the stochastic generation of the lung paths. Throughout this work, variability is expressed by reporting the mean value ± 1 standard deviation (SD). As the mean values -1 SD typically results in negative values, we chose to use "0" as the lower boundary in such cases. Reported generational mean values should be interpreted with care: Mean values in the upper lung generations are computed from a large number (>1000) of individual simulations. However, with increasing generation number, the number of simulations in these generation decreases, as not all initially generated aerosol trajectories penetrate deeper into the lung. Thus, the deposition fractions in the more peripheral airway generations are more strongly influenced by the stochastics of the simulation process. Note that all relative changes of the quantities given below are related to their initial values when entering the lung.
Changes in particle number and diameter as well as deposition fractions for each constituent in the mouth listed in Table 1 are determined primarily by coagulation, with notable contributions by deposition and phase transition. Indeed, about 60% of the particle number is lost due to diffusion driven coagulation in the ET region, while roughly 11% of the condensed mass is deposited (see VG mass, which can be considered non-volatile during inhalation and exhalation). Water starts to condense on the droplets as the relative humidity increases to about 87%, before the aerosol enters the lung. At the same time, some nicotine evaporates in the "no acid" case (about 4%) while it remains in the droplet for the "lower pH" case, where most (roughly 90%) of the nicotine is protonated at a pH level of 7 and thus literally non-volatile.
As a result of mainly coagulation and condensation of water, the median diameter of the aerosol increases by 230% while remaining in the ET region during puffing and mouth-hold. No statistical variation of the results is reported for the ET region as the stochastics of path selection is only considered for the lung.
The evolution of particle diameter and the distribution of number and mass deposition fraction among the airways of the human lung and their variability during inhalation and exhalation are illustrated in the panels of Fig. 2.
Particle number deposition during a complete breathing cycle is plotted in panel A, exhibiting a distinct maximum in lung generations 15 to 25. Note that particle number loss by coagulation is implicitly considered in the deposition calculations. During inhalation, number deposition peaks around generation 20 to 23, while this peak is shifted to slightly lower generation numbers during exhalation. These deposition calculations further illustrate the considerable variability of the deposition fractions among bronchial and alveolar airways. For example, standard deviations at maximum deposition alone are about twice as high as the corresponding mean values. This trend is even stronger for lower lung generations, where the mean values + SD are up to five times higher than the mean values. This indicates that some paths lead to higher deposition fractions from lung generation 5 onward.
Panel B shows the VG relative mass deposition in liquid form as a function of lung generation. As VG features a very low vapor pressure under the given conditions (see Table 2, which shows that there is negligible vapor phase deposition of VG), we chose it as an indicator for condensed mass (analogous to tar for combustion-type cigarettes). While mean deposition patterns very much resemble particle number deposition patterns, their variations are different. At low lung generation numbers, during inhalation, the number deposition shows higher variation than particle mass. This can be explained by considering the effect of coagulation on the particle Table 1 Summary of the changes to the inhaled aerosol in the extrathoracic region: Deposition fractions (df), particle number concentration (N), and diameter (dp). number available. The evolution of the median particle diameter in the airways of the human lung during inhalation and exhalation is illustrated in panel C. The initial strong increase of the median diameter in the oral cavity (not included in the plot; see Table 1) is mainly caused by coagulation. It removes the smallest particles (roughly 70% of all particles are lost due to coagulation in the oral cavity) and increases the size of larger ones. In addition, hygroscopic growth also contributes to the size increase in the ET region. During inhalation, water condenses while the other substances evaporate from the droplets and either deposit at the lung surface or condense on other particles. The growth of larger droplets at the expense of smaller ones (i.e. based on the Kelvin effect) can be studied with the present model Table 2 Generation-integrated deposition fractions for particle number N and all components of the e-liquid in condensed (liq) and vapor (vap) form, respectively, during inhalation (inh) and exhalation (exh). Abbreviations: nicotine (Nic), vegetable glycerol (VG), propylene glycol (PG), benzoic acid (BA).  Fig. 3. Evolution of the relative nicotine mass, i.e. relative to the initial aerosol puff (panels A and B) and of particle phase pH (panels C and D) during a complete breathing cycle, distinguishing between the two cases "lower pH" (panels A and C) and "no acid" (panels B and D).
formulation. However, its contribution is small compared to the addition of water. Similar is the effect of coagulation on growth, as most of the coagulation takes place during puffing and mouth-hold. Note, however, that this should not be generalized, as in the absence of a mouth-hold phase, the significant contribution by coagulation may reshape the aerosol size distribution during inhalation. The maximum median diameter, which is found around lung generation 27, is the result of several processes. The main contribution is the condensation of water vapor and simultaneous evaporation of propylene glycol. Further, size selective deposition and coagulation also play a role. During exhalation, the median droplet size slowly decreases, mainly due to evaporation of propylene glycol and an accompanied loss of water vapor. The loss of nicotine plays a minor role due to its very limited amount (initially 2.9%). Further, deposition of large particles in the alveolar region reduces the median diameter. Fig. 3 shows the evolution of nicotine mass relative to the initial aerosol puff and of particle phase pH during a complete breathing cycle, respectively, which illustrate significant differences between the two cases "lower pH" and "no acid". In the "lower pH" case, nicotine features a much lower effective vapor pressure in the condensed phase. The more nicotine evaporates, the lower is the pH (panel C) and, at the same time, the nicotine volatility. Thus, nicotine is deposited mainly in its condensed form. During inhalation, on average not more than 9% vaporous nicotine deposition was calculated and less than 1% during exhalation. This is contrasted by 75% vaporous nicotine deposition when simulating inhalation of "no acid" aerosol. Further, the pH drop, although much more expressed, has less consequences as nicotine is almost unprotonated (>90%) at a pH of 9.2. To summarize, nicotine deposits prevalently in the condensed phase for the "lower pH" case, while vapor-phase deposition is dominating the "no acid" case. However, note that this is characteristic for the present setup and may not be generalized without carrying out further calculations.
Distributions of the nicotine deposition fractions, relative to the initial nicotine content, of the condensed (liquid) and vapor phases among the airways of the human lung during inhalation and exhalation are plotted in Fig. 4 for the two cases "lower pH" and "no acid".
In the upper part of the lung, the condensed phase nicotine deposition first evolves in a similar fashion for both scenarios (see panels A and B). However, the vapor phase nicotine deposition is dominating the "no acid" scenario (see Table 2). Thus, although the nonvolatile condensed phase deposition is very similar for both scenarios (see VG in Table 2), the liquid phase nicotine deposition decreases strongly with increasing lung generation number as the nicotine is depleted in the condensed phase due to evaporation.
Distributions of the propylene glycol (PG) deposition fractions, relative to the initial PG content are depicted by Fig. 5. It considers condensed (liquid) and vapor (gas) phase deposition among the airways of the human lung during inhalation and exhalation, Fig. 4. Distribution of the nicotine deposition fractions, relative to the initial nicotine content, of the condensed (liquid) (panels A and B) and vapor phases (panels C and D) among the airways of the human lung during inhalation and exhalation, distinguishing between the two cases "lower pH" (panels A and C) and "no acid" (panels B and D). distinguishing between the two cases "lower pH" and "no acid". These calculations demonstrate that PG deposition rates in condensed and vaporous form are almost identical. Although minor variations are caused by the stochastics of lung path creation, they cancel out when looking at generation-integrated values (see Table 2).
The results of the deposition calculations for particle number and all constituents of the EC aerosol, i.e. nicotine, glycerine, propylene glycol and benzoic acid, in both condensed (liquid) and vapor form, are summarized in Table 2 for both the "lower pH" and the "no acid" case and both inhalation and exhalation. Total deposition fractions were obtained by integration of the computed deposition patterns over all airway generations. In summary, several differences as well as similarities can be observed between the "lower pH" and the "no acid" cases: (i) Number deposition fractions for inhalation and exhalation are hardly affected by the choice of the pH value.
(ii) Significant differences can be observed for nicotine liquid inhalation and exhalation, and nicotine vapor inhalation, while only minor differences exist for nicotine vapor exhalation. (iii) For comparison, only minor differences exist between the liquid and gas phases of glycerine and propylene glycol for both inhalation and exhalation phases. (iv) Benzoic acid liquid and vapor deposition fractions follow the behavior of nicotine liquid and vapor for both inhalation and exhalation.

Conclusions and discussion
In summary, two different e-cigarette formulations, "lower pH" and "no acid", were applied for the simulation with the coupled ADiC/IDEAL model system to calculate the evolution of the EC aerosol during puffing, inhalation and exhalation as well as its interaction with the lung. Deposition of the volatile components of the e-liquid were determined for condensed as well as vaporous form. Due to the nature of the IDEAL code, taking into account several thousands of stochastically generated individual lung paths, average values as well as corresponding statistical fluctuations were simulated.
Evolution of both aerosols is very similar with regard to particle diameter, total mass deposited or evolution of particle number. In contrast, the deposition of nicotine is very different for the two cases and seems to be strongly governed by the evolution of particle phase pH. In the "lower pH" case, starting at an arbitrary pH level of about 7, nicotine deposition is mainly found to be caused by nicotine-containing particles depositing on the lung surface. However, the formulation that does not contain any acid contains nicotine that is much more volatile. As a result, nicotine deposits predominantly in vaporous form in this case.
In the case of tobacco burning cigarettes and especially for electronic cigarettes, nicotine is believed to be the most addictive substance (WHO, 2015). According to several studies, e.g. Lauterbach, Bao, Joza, and Rickert (2010), that deal with the amount of so-called free-base nicotine (non-protonated nicotine in the particle phase) in mainstream cigarette smoke, only a small fraction of nicotine is uncharged. As a result the nicotine activity coefficient is considerably smaller than 1. This suggests that only a small fraction of nicotine is available for evaporation from the particle phase. However, Armitage, Dixon, Frost, Mariner, and Sinclair (2004) found that a considerable fraction of nicotine deposits in the human lung in vaporous form. Thus, transformation of protonated to non-protonated nicotine is a key process for understanding nicotine delivery to the human lung and requires the description of particle phase chemistry. Although the concept of particle phase pH has been discussed controversially for cigarette and e-cigarette aerosols, it has been shown that the application of this concept may close the gap between simulation results and observations in denuder tube experiments (Pichelstorfer & Hofmann, 2015).
The physical mechanisms determining the fate of inhaled EC aerosols in the human respiratory tract adopted in the present study are practically the same as those assumed in Asgharian et al. (2018a,b). However, despite considering more or less the same physical processes describing deposition and the same morphometric lung data-set, specific differences of the present model are: (i) the use of the full size distribution entering the mouth upon puffing, compared to the assumption of an average diameter; (ii) the application of a stochastic, asymmetric multiple path lung model, compared to a deterministic single path model; (iii) the solution of some of the model equations by numerical methods, e.g. Monte Carlo techniques, compared to the solution of differential equations.
Aerosol dynamics simulations reveal that the particle diameter does not remain constant, especially during puffing and inhalation and, to some extent also during exhalation. Growth caused by condensation of water vapor is not taking place instantaneously (Winkler-Heil et al., 2017). The present simulations suggest that condensational growth significantly changes the particle diameter (and thus also composition) during the puffing, mouth-hold and inhalation phases. The reason is the large potential of the considered e-liquids to grow under humid conditions. Thus dynamic changes of the particle diameters primarily during the inhalation phase must be considered in calculations of cigarette aerosol deposition in general, and specifically in e-cigarettes (Pichelstorfer et al., 2016).
For most of the quantities shown in the figures, the variation is about 2-3 times of their average value. For example, the maximum of the average particle number deposition per generation is about 2%, while the average + 1 SD is about 6%. This also applies, though slightly less, for vapor phase deposition. This finding is of importance when considering toxicological effects. Surface deposition levels in certain regions of the lung may therefore be considerably higher (or lower) than suggested by their average values. These significant statistical fluctuations of the deposition patterns illustrate the inherent anatomical variability of the human airway structure within a given lung (intrasubject variability) and between different lungs (intersubject variability).
The present work applies a 1-dimensional solution for solving the particle transport. This approach requires the assumption of an idealized flow field in the lung. On the other hand, it allows for a very detailed description of the lung itself and relatively simple implementation of the description of aerosol dynamics. Limitations resulting from this approach are discussed in section 3 "Limitations".
Finally, we want to note that the present work does not aim to provide a look-up table for deposition fractions or dosimetry of ecigarette aerosol in general. As discussed in sections "2.3.5 Model inputs" and "3. Limitations", pronounced variability is observed, however, not accounted for in the present work, for many input parameters other than lung geometry. Thus, the present work focuses on the variability caused by anatomical variation and represents a case study with regard to other potentially varying parameters. Further studies are needed to take into account variations of other parameters like puffing duration or e-cigarette liquid formulation.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments
This work was supported by the Austrian Science Funds (FWF) under grant J-4241 Schrödinger Programm.

Deposition in turbulently mixed chamber
First order loss coefficients are calculated using the following expression (Lai & Nazaroff, 2000): where v dv , v du , v dd are the deposition velocities on a vertical surface, an upward horizontal surface and a downward horizontal surface, Fig. 5. Distribution of the propylene glycol (PG) deposition fractions, relative to the initial PG content, of the condensed (liquid) (panels A and B) and vapor phases (panels C and D) among the airways of the human lung during inhalation and exhalation, distinguishing between the two cases "lower pH" (panels A and C) and "no acid" (panels B and D).
respectively, and A v , A u and A d are the respective surface areas. Note that deposition velocities are dependent on the friction velocity (u*) which is the only parameter applied in the model. It characterizes the nature and the intensity of near surface turbulent flow: where ν is the kinematic viscosity of air, U is the speed parallel to the surface and dU dy ⃒ ⃒ ⃒ ⃒ y=0 is the average velocity gradient which can be derived for a given freestream velocity U ∞ and a plate length L (see, e.g. Schichtling, 1979): where ρ a is the density of the gas.
In the present work, we used the length of the mouth cavity to describe the plate length and the mean freestream velocity was calculated by: where A cig is the e-cigarette cross section and V mouth is the mouth volume. For an explicit description of how to derive the deposition velocities see (Lai and Nazaroff, 2000).

Equations determining deposition efficiency
Deposition of particles by Brownian motion in upper bronchial airways was determined by the empirical equation proposed by Cohen and Asgharian (1990) for laminar developing flow. In more peripheral bronchiolar and acinar airways, Ingham's (1975) equation for diffusion under parabolic flow conditions was applied, while for expiratory flow a uniform flow profile due to secondary flows was assumed for all airways (Ingham, 1975). Deposition by inertial impaction in bronchial airways was calculated according to Cai and Yu (1988), assuming a uniform profile in upper bronchial airways, a parabolic profile in distal bronchiolar and acinar airways, and again a uniform profile upon expiration. In the trachea, the effect of the laryngeal jet was considered by the empirical equation proposed by Chan, Schreck, and Lippmann (1980). Finally, deposition by sedimentation was computed by the equations proposed by Finlay (2001) for well-mixed plug flow in upper bronchial airways, Poiseuille flow in bronchiolar airways, laminar plug flow in alveolated airways, and well-mixed plug flow in all conducting airways during exhalation. Werner Hofmann graduated from the University of Vienna in 1973 with a PhD in Nuclear Physics. He then joined the Institute of Physics at the University of Salzburg, where he is still actively pursuing his research activities after his retirement as a professor in 2010. During his academic career he spent two sabbaticals at the University of Nebraska and Duke University in North Carolina. He still holds adjunct professor positions at the University of Nova Gorica, Slovenia and the Queensland University of Technology, Brisbane. So far he has authored about 280 publications in peer revised journals and numerous contributions to proceedings of international conferences.