The aerosols are generated using a collision-type atomizer (Model 3076, TSI) with filtered air as the gas source, maintained at P = 2 bars. To minimize water vapor content while preventing particle loss, the aerosol stream is directed through two consecutive diffusion dryers (Model 3062, TSI). The resulting airflow rate is Q A G = 1.2 L/min (Figure 1a), measured using a flow meter (Model 4100, TSI). Such airflow rate is about three orders of magnitude higher than the desired bronchiolar flow rate (Q ₀) under quiet breathing conditions (Pedley, 1977) to feed the first generation of the four microchannel models positioned on a single microscope slide (Figures 1a,b). To achieve the required flow conditions, a stepwise air reduction and splitting mechanism is implemented, as detailed in the table of Figure 1c. The first stage of flow reduction is achieved using an assembled array of tubes (all apparatus dimensions are provided in the Supplementary Material (SM) Supplementary Figure S1a and S1b), which gradually reduces the airflow in a manner analogous to airway branching in the lungs, from the bronchi down to the bronchioles. This system, referred to as the “air reducer” (Figure 1a), expels approximately ∼80% of the excess airflow (Q _bleed) through a 90° outlet into a chemical hood. The remaining airflow is symmetrically split to feed two parallel lung-on-chip models (Q _chip1 and Q _chip2). Each of these flows supplies four microchannels (Q _micro) through their respective inlet ports (Figure 1b).

Before entering the first microchannel generation (i.e., Gen. 0, see Figure 1d and top view in Figure 2c) the bulk airflow passes through a small vertical pipe for additional reduction Q _pipe. The desired flow rate into Gen 0, .designated as Q ₀, is then drawn under a constant steady flow into the microchannels via a syringe pump (PHD Ultra, Harvard Apparatus) connected to the suction outlet of each chip (Figure 1b; Supplementary Figure S1a). The constant flow approximation in the model is sufficient to accurately represent physiological breathing, as further justified in the device design section. The final reduction step ensures that Q ₀ is only 6% of Q _micro and is directed to the model inlet to reach the desired quite breathing rate, while the remaining 94% is expelled. The syringe pump simultaneously withdraws a total airflow of Q _S = 4Q ₀. Within each microchannel, the airflow is further subdivided in a manner mimicking bronchiolar bifurcation. The airway-on-chip design features individual airway trees with three bifurcating generations, includes both normal and constricted segments.

Fluorescent polystyrene latex (PSL) particles are used to investigate size-dependent aerosol deposition (ρ _p = 1,050 kg/m³, Fluoromax red and green, fluorescent microspheres, 1% solid, Thermo Scientific). The use of PSL particles is well established in aerosol exposure studies involving airway-on-chip models (Bessler et al., 2023; Elias-Kirma et al., 2020; Fishler et al., 2015) due to their well-characterized physiochemical properties and charge acquisition during aerosolization. To examine the effects of particle size and electrostatic charge on deposition, distinct monodisperse PSL sizes are used: d p = 0.2, 0.5 and 2 μm. While the smallest particles (i.e. 0.2 μm) are known to be significantly affected by Brownian diffusion, the largest ones (i.e. 2 μm) are subject to considerable sedimentation in small airways. In contrast, 0.5 μm particles sit in a size range where airborne transport is primarily driven by convection such that aerosols are exhaled under normal breathing and deposit in minimal amounts due to weak diffusive or gravitational forces (Hofem et al., 2015; Sznitman, 2013). Furthermore, and most relevant to the present study, PSL particles are selected due to their natural tendency to acquire electrostatic charge during the aerosolization process. The “Charged Group” consists of particles exiting the aerosol generator after diffusion drying, with an estimated mean charge of >100 e based on previous experimental data (Forsyth et al., 1998; Whitby and Liu, 1968). In contrast, the “Neutralized Control Group” is obtained by employing an electrical ionizer (Model 1090 MSP, TSI) and ensuring charge neutralization following a Boltzmann charge equilibrium. Under these conditions, the mean charge per particle is expected to be q ∼ 1 e for 0.2 μm and q ∼ 3 e for 2 μm particles (Hinds and Zhu, 2022). Note that these levels are insufficient to induce significant electrostatic deposition, allowing for direct comparison between charged and neutralized aerosols (Osman et al., 2015).

Prior to exposure (see previous section), polystyrene latex (PSL) particles are suspended in deionized (DI) water at a high concentration (∼10⁷–10⁸ particles/mL) to ensure sufficient deposition data for ensemble statistics within reasonable exposure times (∼30–90 min, depending on particle size and charge conditions) (Fishler et al., 2015). The particle suspension is prepared at a weight concentration of 5%, carefully selected to maintain an aerosol with <1% aggregate formation. This assumption follows from the expected evaporation of water droplets exiting the atomizer, characterized by a count mean diameter of 0.35 μm and a geometric standard deviation of 2 μm, as specified by the manufacturer. The size distribution of the generated aerosol is assumed to follow a log-normal distribution with monodisperse spherical particles, consistent with previous studies (Fishler et al., 2015; Raabe, 1968). To minimize aerosol aggregation due to the high particle concentration, each suspension undergoes 30 min of sonication in an ultrasonic water bath (Elma Elmasonic S10) immediately before experimentation. To prevent particle sedimentation during experiments, a magnetic stirrer continuously operates inside the suspension container.

Following each aerosol exposure experiment, the bronchiole microchannels were examined using an inverted fluorescent microscope (Nikon Eclipse Ti) at varying magnifications (×2, ×10 and ×20) depending on particle size (i.e., 0.2 μm, 0.5 μm and 2 μm respectively), see Figures 1e,f. Using ANDOR Zyla sCMOS camera, we captured the entire model by stitching them together with zero overlap, using large free-shape images function. To reduce the image memory space, we used 4 × 4 binning. The quantification of deposition involved identifying the precise planar locations (x, z) of individual deposited fluorescent particles through digital processing of local intensity maxima using ImageJ software. Subsequently, the particles were counted with correspondence to their two-dimensional (2D) location and the sum of deposited particles were obtained for each airway generation. The deposition fraction was then calculated as the percentage of deposited particles in a specific airway generation out of the total number of deposited particles across the entire model. This approach enables us to determine the spatial distribution of particle deposition within the in vitro model.

The in vitro airway-on-chip platform consists of four planar, symmetric airway path trees spanning four generations 0–3, as schematically shown in Figures 1, 2c. The tree design broadly mimics the deep bronchial airway branches (i.e., bronchioles) in the distal region of the conducting zone, with hydraulic diameters of <1 mm (Figures 1a,b). Each airway segment consists of square cross-sections, with dimensions summarized in Figure 2b. The chosen generation dimensions were based on morphometric measurements representative of typical length and diameter of an adult human lung, corresponding to generations 12–15, following the seminal works of Weibel (1963), Horsfield et al. (1971). These airway dimensions mimic the deep bronchial regions, referring to airways with a diameter smaller than 2 mm, which are associated with constricted diseases (McNulty et al., 2014). While the Weibel A model assumes idealized cylindrical airways, actual distal bronchiolar cross-sections in the human lung are often irregular and deviate from perfect circles (Hogg, 2004; Wang et al., 2020). In our airway-on-chip design, we employed square cross-sections (Figures 2c,d) as a deliberate engineering choice to enable precise microfabrication (Figure 1b), optical access for high-resolution imaging (Figures 1e,f), and consistent channel alignment across multiple generations. Importantly, flow behavior and particle transport in rectangular microchannels remain comparable to those in circular channels when the hydraulic diameter (defined as D _ℎ = 4A/P, where A is the cross-sectional area and P is the wetted perimeter) and flow conditions (Re) are matched. To ensure physiologically relevant flow, we designed each generation as a square channel to have a hydraulic diameter D _ℎ equal to that of the corresponding cylindrical segment in the Weibel model (see Figure 2b). While minor deviations in flow may occur near corners due to reduced velocities, these are localized and do not significantly affect the overall transport and deposition behavior of aerosols. Therefore, the underlying physics remain valid despite the geometric simplification. The first airway (Gen. 0, see Figure 1d and see top view Figure 2c) serves as the model’s inlet, ensuring fully-developed flow conditions over a distance exceeding the anticipated entrance length L _e by more than two orders of magnitude at such low Reynolds numbers ( L e ∼ 0.06 Re D 0 ). These conditions result in fully-developed laminar Poiseuille flow profiles within the model’s bronchioles (see Figures 2c,d).

Approximating a constant inhalation condition may be assumed for modeling airflow in the deep lung regions during quiet breathing, as the Womersley number (Wo = fL ² ρ _f /μ _f ), which is a dimensionless expression of the oscillatory airflow frequency (f ∼ 0.25 Hz for quite breathing) in relation to viscous effects, remains well below unity in our bronchial region of interest (Wo ≪ 1), indicating quasi-steady flow conditions that can be approximated as steady flow conditions (Ménache et al., 2008). However, it is important to note that while the airflow can be considered steady, the absence of exhalation in the model may lead to overestimating particle retention, as particles that would otherwise exit the lungs during exhalation are not accounted for. The Gen. 0 segment terminates in a bifurcation with a 25° split angle. Beyond this point, the airway tree bifurcates symmetrically across all generations, maintaining a consistent 37° branching angle along the symmetry line (Ménache et al., 2008). The model branches into two distinct pathways: one representing normal anatomical bronchioles and the other featuring a 30% constriction, designed to mimic smaller airways or regions affected by obstructive pulmonary diseases (Figure 2). Airflow regulation is achieved using a syringe pump connected to the model outlet port (Figure 1b). The present design enables the investigation of electrostatic deposition by incorporating a conductive coating, while electrical grounding prevents charge buildup (Bessler et al., 2023).

Models were fabricated using 3D-printed molds printed using a PRUSA SL1 3D printer and 3DM-ABS Orange Tough Resin), see Supplementary Material Supplementary Figures S1c and S1d. The models were filled with polydimethylsiloxane (PDMS), inspired by recent microfluidic lung airway models (Bessler et al., 2023; Elias-Kirma et al., 2020). PDMS, a well accepted material for microfluidic devices, was selected for its optical transparency and high flexibility, allowing for easy molding into complex shapes and structures. The printing settings included a layer height of 25 μm and an exposure time of 10 s. Post printing, the molds were cleaned with isopropyl alcohol, and to solidify the resin the molds were cured at 60 °C for 1 h under ultraviolet light (Formlabs Form Cure FH-CU-01 curer). After pouring the PDMS, the PDMS-filled molds were placed in desiccators for 1 hour to release trapped air bubbles. The PDMS models were peeled after 24 h of resting within the mold under ambient conditions. The outlet port for the syringe insertion (Figure 1b) was created using a 2 mm biopsy punch (Miltex, 3331).

The electrical conductivity of lung parenchymal tissue closely approximates that of saline water σ H 2 O ∼ 1 S/m (Osman et al., 2015). This high conductivity allows electrical charges within lung tissue to redistribute rapidly, facilitating the manifestation of image charge phenomena (Melandri et al., 1983). In contrast, glass has a significantly lower conductivity ( σ glass ∼ 10⁻¹¹ S/m), which necessitates an alternative approach to achieve a physiologically relevant conductive environment.

To address this limitation, we applied a conductive indium tin oxide (ITO) coating to our microscope slide, following previous methodologies (Bessler et al., 2023). The ITO pattern was selectively deposited only in regions corresponding to the airway microchannels (see Supplementary Figure S2) to avoid interfering with the adhesion of the PDMS to the glass. The bonding process was performed via a 1 min plasma treatment (ETP, INC. Model BD-20). The ITO coating, provided by Huizhou Konshen Glass Co., Ltd., had a sheet resistivity of approximately R s ∼ 10   Ω / m 2 , significantly enhancing the surface conductivity to σ ITO ∼10² S/m (see Supplementary Material Supplementary Table S1 for details on sheet resistivity). This elevated conductivity enables the manifestation of image charge effects, more faithfully replicating the electrostatic interactions expected in lung tissue. Notably, the electrostatic relaxation time scale ( τ e ∼ 1 pico s) in real lung tissue is orders of magnitude shorter than the characteristic timescale of a particle ( τ p ∼   1 μs). Thus, achieving a sufficiently conductive environment in our model mimics more faithfully electrostatic behavior. Additional details on the relevant timescales are provided in the Supplementary Material (see Supplementary Material Conductive Layering of the Microchannel). To maintain a consistent electrical setup across experiments and prevent excess charge buildup during exposure assays, the coated conductive glass was electrically grounded using a custom external connector wired to the inner surface of the microfluidic channel (See Supplementary Figure S1a).

We propose a canonical model featuring rectangular channels connected in series according to Weibel dimensions mimicking our in vitro model (Supplementary Figure S5). In each channel, the streamwise flow (u _f ) is approximated using the laminar Hagen–Poiseuille velocity profile (see Supplementary Figures S6 and S7, and derivations for flow in a square cross-section chanel) at low Reynolds number (Re ∼ 3), corresponding approximately to quiet breathing of a healthy person at these deep lung generations. The numerical simulations are initialized with an aerosol starting at a vertical distance of height y t = 0 = 0.025 D 0 , chosen to align with the image charge method ( y ≪ D 0 . Additionally, this distance from the tissue corresponds to approximately 10% of the annular volume in an ideal cylindrical airway (see Supplementary Figure S8 and derivation in the Supplementary Material), making it a widely accepted threshold for assessing whether electrostatic deposition is significant, as supported by previous studies (Finlay, 2021). This value will later be chosen as the characteristic length (L _c = y ₀ = 0.025D). The initial velocity is set based on the corresponding position in the Hagen–Poiseuille flow profile. Following each iteration, based on the time interval Δt, a 1D random Brownian step ( ±   2 D d i f f Δ t was added to each calculated location, for each axis separately. Here, D _diff represents the Stokes–Einstein diffusion coefficient. To ensure convergence, the time step was set to 0.1 ms (Δt = 0.1 ms). To validate this choice, we performed simulations with progressively smaller time steps, ranging from 10 ms down to 1 μs, confirming that the solution converges (see Supplementary Figure S8a). Additionally, to account for the inherent error in the Explicit Euler method, where terms of O (Δt ²) and higher are neglected, we tracked the local truncation error (LTE) at each iteration. The resulting global truncation errors (GTE) remained around single microns and are summarized in the Supplementary Material (Supplementary Figure S8b; Supplementary Table S3). Additionally, to validate our numerical approach, we compared the results with an analytical solution derived from scaling laws. As the time step (Δt) decreases, both solutions converge well (see Supplementary Material Supplementary Figure S8 for comparison). Since our time intervals are much larger compared to the examined particle relaxation time (τ _p ≪Δt), the particle’s velocity instantaneously follows the surrounding flow velocity ( x ˙ = u f ).

Non-dimensional analysis is helpful for understanding aerosol transport and its leading physical mechanisms, in particular when electrostatic forces are involved. This approach expresses the dominant mechanisms in terms of characteristic length (L _c ∼ D), velocity (U _c ∼Q/A) and time (t _c ∼ L _c /U _c ) scales relevant to the bronchioles. We revisit the governing mechanisms of aerosol deposition in our model using dimensionless parameters (Sznitman, 2022; Tsuda et al., 2013). The deposition of small (approximately <0.5 µm) particles in the bronchioles occurs mainly via Brownian motion, characterized by the inverse Péclet number: Pe ⁻¹ = D _diff /U _c L _c . For larger particles (>1 µm), bronchiolar deposition is primarily governed by gravitational sedimentation, quantified by the gravity number H = v _ter /U _c , where v _ter is the settling velocity given by v _ter = τ _p g, and τ _p represents the particle relaxation time. Concurrently, particle inertia is captured by the particle Stokes number: Stk = τ _p U _c /L _c .

We then extend the discussion to when electrostatic charge is present. When Coulomb forces arising from particle–wall electrostatic interactions introduce an additional dimensionless group relevant to represent the Induced Charge versus the airflow convection (Finlay, 2001): Inc = (Bk _e /U _c ) (q/L _c )². By applying dimensional analysis to Newton’s second law, while omitting the noise introduced by Brownian motion (which can be superimposed subsequently), we obtain the non-dimensional governing equations for particle dynamics, Equations 3, 4, presented below: y ^ : H g ′ − Inc − y ˙ ′ = Stk  y ¨ ′ (3) x ^ : − x ˙ r e l ′ = Stk  x ¨ ′ (4) where x ˙ r e l ′ is the particle’s normalized horizontal velocity relative to the flow. Non-dimensional analysis in the deep bronchioles reveals that impaction is negligible (Stk ≪ Inc, H), and primarily occurs in the upper lung regions (Stahlhofen et al., 1989). As a result, the governing Equations 3 and 4, can be simplified by neglecting the inertial term, reducing the overall complexity of the system to Equations 5 and 6. y ^ : m g − k e q 2 4 y 2 − y ˙ B = 0 (5) x ^ : x ˙ = u f y (6)

Whereas the horizontal particle trajectory follows a straightforward path, showcasing how the particle velocity aligns with that of the surrounding flow (Equation 6), the non-linear equation in the vertical direction (Equation 5) remains more complex and requires further analysis.

We define the characteristic length scale y E as the “Equilibrium Distance” from the wall, where the electrostatic force balances gravity, i.e., F e y E ≡ F g (see Figure 3A). It is important to note that this length scale is specific to each particle’s size and charge levels and is given by y E = k e q 2 / 4 m g . To better characterize electrostatic deposition, we define the normalized distance from the tissue as: η = y t / y E where the initial condition is given by η 0 = y 0 / y E . Hence, electrostatic deposition against gravity is expected only within the range 0 ≤ η < 1 . A careful simplification of Equation 6 (detailed in the Supplementary Material, in section “Analytical and Mathematical Derivation”) leads to the following form (Equation 7): y ^ : η = η 0 + v t e r t y E − 1 2 ln η − 1 η + 1 · η 0 + 1 η 0 − 1 (7)

This non-dimensional equation is no longer a differential equation, making it significantly easier to solve and hence analyze particle behavior over time. We can write it dimensionally as: y ^ : y t = y 0 + v t e r t ⏞ free fall − y E 2 ln y t − y E y t + y E · y 0 + y E y 0 − y E   ⏞ Electrostatic  (8)

This Equation 8 describes displacement through two distinct terms. The first term on the left represents the classic free-fall sedimentation over time, while the nonlinear logarithmic term captures the opposing electrostatic attraction counteracting gravitational pull. Despite this clear separation, the Equation 8 is implicit in y t meaning an explicit expression for the particle’s position as a function of time is not straightforward. Yet, for very small particles, highly charged ones, or in a near proxinity to the tissue, gravitational settling is negligible (H≪Inc), the equilibrium distance y E extends toward infinity in comperison to y t . Under these conditions η , η 0 ≪ 1 , we can expand the logarithmic term in a Taylor series around zero up to the fourth order ( O y 5 ), leading to a simplified explicit results presented in Equation 9: y t = y 0 3 − 3 4 k q 2 B t 1 3 (9)

This result offers a simplified yet direct analytical expression for the vertical displacement of the particle trajectory, providing a fundamental framework to understand how electrostatic charge influences deposition dynamics under Stokes drag.

An alternative approach to Equation (8) can be formulated using a non-dimensional framework. By considering the deposition condition y t f = 0 , the equation reduces to involve only two non-dimensional parameters: Π 1 = v t e r t f / y E , representing the free-fall distance normalized by the equilibrium distance, and Π 2 = y 0 / y E (analogous to η 0 ). This reformulation leads to Equation 10 Π 1 = 1 2 ln Π 2 + 1 Π 2 − 1 − Π 2 (10)

Equation 10 provides a generalized expression for the final deposition time t f , accommodating any combination of charge, particle size and shape across all normalized initial conditions (Figure 3B). Upon reorganizing, the final deposition time is explicitly given in Equation 11: t f = y E v t e r 1 2 ln y 0 + y E y 0 − y E − y 0 (11)

We now turn our attention to the largest particles, i.e. 2 μm. Recalling that H·Pe ∼ 235, these particles, when neutralized, are governed almost entirely by gravitational pull. In this regime, the planar orientation of our model becomes particularly relevant, as gravity acts in a fixed direction, and the absence of 3D airway curvature may influence the observed deposition patterns for clinical interpretation of the results. In order to distinguish the electrostatic forces from gravitational effects in our experiments, we conducted experiments with the model turned upside-down such that the ITO-coated surface is at the top, hence electrostatic may be interpreted as a “lifting force” against the downward pull of gravity. As expected, in the neutralized case, no deposition occurred on the top even when we extended the exposure time (see missing bars in Figures 4c,d). To validate this, as discussed in the Supplementary Material, we repeated the experiment with the model flipped to align with gravity and observed deposition biased by gravity direction, both with and without charge yields significant deposition which is hard to distinguish (see Supplementary Material Supplementary Figure S4c and S4d). This confirms the significant role of gravity in the bronchiolar region, as suggested previously (Ma and Darquenne, 2012). As expected, gravitational deposition occurs only in the direction of the gravitational pull. This can be visually shown in our numerical analysis of particle trajectories (see Figure 5d, Inc/H ∼ 0.01) where particles do not reach deposition at the top of the model when electrostatic forces and diffusion are absent. The particle follows its free-fall trajectory (see Figure 5d, dashed blue line) away from the tissue, with only slight deviations caused by weak diffusion. Its airborne time is approximately half a second (0.48 s), which is shorter compared to the 0.2 µm diffusive particle (1.200 s). This difference arises because the free-falling particle enters closer to the channel centerline, leading it to gain horizontal velocity as it moves into the high-speed region of the parabolic Poiseuille flow profile, resulting in a faster outward convection.

When these same 2 μm particles are charged, they have a chance to deposit. A distinct electrostatic-driven screening effect emerges leading to deposition within the bronchial model (Figure 4b). This suggests that electrostatic forces are enough to counteract gravity, encouraging greater bronchiolar deposition, even in regions where gravitational settling is unlikely to occur. We recall that the characteristic length scale y E represents the “Equilibrium Distance” from the wall at which the electrostatic force balances gravity. This length scale is crucial for understanding particle dynamics under the influence of both electrostatic forces and gravitational effects. Mathematically, it serves as a so-called separatrix, distinguishing between two distinct scenarios: (i) deposition due to electrostatic attraction versus (ii) convective transport through the channel to the next-generation. Most likely, in our case, this results in y E   ∼ 25   μ m (Forsyth et al., 1998; Whitby and Liu, 1968). Based on our in-vitro deposition findings (Figures 4c,d), we can infer that regions where y t < y E are populated with aerosols that are most likely deposited due to electrostatic contributions. In comparison to the screening observed for the 0.2 µm particles, the gradient of the screening is less intense here. This can be attributed to the fact that while Inc⋅Pe is not dependent on particle size, the ratio Inc/H ∝ d _p ⁻³ is highly size-dependent. For the larger 2 µm particles, this strong dependence reduces the effectiveness of electrostatic deposition, as their increased dimensions hinders the influence of electrostatic forces. This is clearly shown and aligns with an in vitro study (Bessler et al., 2023) and an additional upper airway CFD study, where Inc/Stk ∝ d _p ⁻³ follows the same dependency (Koullapis et al., 2016). As a result, more substantial deposition in proximal airways occurs for smaller particles, reinforcing the weight of electrostatic forces for smaller aerosols (Koullapis et al., 2016). Accordingly, a highly significant interaction between particle size and charge (p < 0.0001, for 2 μ m β∼21% for Gen 1, see Supplementary Material Supplementary Table S6 for more), indicating that the effect of particle size on deposition depends on charge, and vice versa. Hence, charge does not influence all particle sizes uniformly, certain sizes are more affected by electrostatic forces than others.

Although electrostatics generally increased overall targeted bronchiolar deposition, it is important to note that across all tested combinations of charge states (neutralized and charged) and particle sizes (0.2, 0.5, and 2 µm), no significant increase in deposition was observed within the constricted regions compared to the normal airways. This finding indicates that, within the limited range of particle sizes and electrostatic charge levels explored in this study, we did not identify conditions that preferentially enhance deposition in constricted regions. Specifically, for 2 µm particles, deposition in the constricted airway was consistently about half of that observed in the normal airway, regardless of charge. These results underscore the ongoing challenge of achieving efficient, targeted aerosol delivery to the bronchioles, especially in distinguishing deposition behavior between constricted and non-constricted geometries. Future studies should explore a broader range of particle sizes, electrostatic charge levels, and airway constriction severities, as this work scope focused mainly on the effect of electrostatics in the bronchiole region and consider only one case of a mild constriction. Given that electrostatic forces become increasingly dominant at smaller length scales, we hypothesize that a promising window of opportunity remains to be uncovered—one that could enable selective targeting of diseased airways through optimized charge and size combinations.

In Figure 5b, we examine the particle trajectory under conditions where electrostatic and gravitational forces are of comparable magnitude (Inc ∼ H), as suggested by Equation (8). The simulation reveals that charged particles can ascend against gravity and deposit on airway surfaces that would otherwise remain unaffected. This behavior aligns with our in vitro findings (Figure 4a), where charged particles exhibited enhanced deposition compared to their neutralized counterparts. Notably, since the release height y 0 is close to the equilibrium distance ( y E   ∼ 25 %   o f   D 0 ), the particle spends most of its airborne time hovering near this equilibrium position. However, beyond a critical threshold, the electrostatic force becomes strong enough to overcome gravitational pull entirely, and the particle trajectory converges with the analytical prediction of Equation (9) (green dashed line). As a result, the final velocity of the particle moving upward against gravity surpasses its free-fall terminal velocity by three orders of magnitude ( v t e r ≪ v f ). In summary, with clinical relevance, electrostatic forces expand the effective bronchial tissue area available for aerosol delivery by enabling particles to reach the top surfaces throughout the airway cross-section, and not only those who favor gravity direction. Hence, electrostatics demonstrate substantial potential for enhancing drug delivery efficiency, particularly in regions where conventional mechanisms fall short.

As a final step, we demonstrated the non-dimensional dependency (Figure 2b) that applies to all deposition scenarios involving electrostatic forces versus gravity. The result of our theoretical non-dimensional plot perfectly aligns with the deterministic trajectory and the one predicted by the numerical simulation, yielding a deposition time of t f = 0.239 s in both approaches. This framework serves as a valuable and comprehensive summary, capturing the interplay between particle shape, size, charge, and initial conditions through two dimensionless parameters: y ₀/y _E and v _ter t _f/y _E . This relationship provides an effective tool for estimating the deposition time t _f , which can be further related to the convective timescale. This representation also accurately predicts deposition times for our specific cases, such as 0.2 and 2 µm particles, showing strong agreement with our numerical solutions. Notably, deposition time approaches infinity when particles are uncharged or if they are initially positioned near the equilibrium stagnation point y E . The initial height y 0 plays a crucial role in determining deposition time. In constricted airways, particles are more likely to start from a lower initial height, smaller distance leading to shorter deposition periods. Additionally, as observed, the particle’s terminal velocity v t e r , influenced by its shape, remains a key factor in deposition dynamics. While our simulations assume spherical particles, it is important to note that alternative particle shapes could exhibit significantly different behaviors (Kleinstreuer et al., 2008; Kleinstreuer and Zhang, 2010). For example, long, straight fibers interact differently with the drug medium (Asgharian and Yu, 1988; Harris and Timbrell, 1975; Shachar-Berman et al., 2018). Asbestos fibers, for instance, with lengths of 50–200 μm, are known to deposit extensively in the lungs and cause severe health conditions (Timbrell, 1965; Donaldson et al., 1989). In the context of electrostatic, in vivo studies in rodents have found much higher acinar deposition in charged asbestos fibers compared to neutralized ones (Davis et al., 1988). This suggests that combining non-spherical shaped particles with electrostatic calls for future investigation (Bessler and Sznitman, 2024). While most accepted approaches neglect the role of electrostatic effects, a growing body of evidence suggests that they play a significant role in altering deposition (Kwok and Chan, 2009).