In Silico Optimization of Fiber-Shaped Aerosols in Inhalation Therapy for Augmented Targeting and Deposition across the Respiratory Tract

Motivated by a desire to uncover new opportunities for designing the size and shape of fiber-shaped aerosols towards improved pulmonary drug delivery deposition outcomes, we explore the transport and deposition characteristics of fibers under physiologically inspired inhalation conditions in silico, mimicking a dry powder inhaler (DPI) maneuver in adult lung models. Here, using computational fluid dynamics (CFD) simulations, we resolve the transient translational and rotational motion of inhaled micron-sized ellipsoid particles under the influence of aerodynamic (i.e., drag, lift) and gravitational forces in a respiratory tract model spanning the first seven bifurcating generations (i.e., from the mouth to upper airways), coupled to a more distal airway model representing nine generations of the mid-bronchial tree. Aerosol deposition efficiencies are quantified as a function of the equivalent diameter (dp) and geometrical aspect ratio (AR), and these are compared to outcomes with traditional spherical particles of equivalent mass. Our results help elucidate how deposition patterns are intimately coupled to dp and AR, whereby high AR fibers in the narrow range of dp = 6–7 µm yield the highest deposition efficiency for targeting the upper- and mid-bronchi, whereas fibers in the range of dp= 4–6 µm are anticipated to cross through the conducting regions and reach the deeper lung regions. Our efforts underscore previously uncovered opportunities to design the shape and size of fiber-like aerosols towards targeted pulmonary drug delivery with increased deposition efficiencies, in particular by leveraging their large payloads for deep lung deposition.

where u is the velocity vector,  is the mass density, P is the pressure,  is the dynamic viscosity.

RANS based Realizable k- Turbulence Model
The modeled transport equations in the realizable k- model are defined as follows [34,[35][36]: and 2 1 2 1 3 where 1 max 0.43, , =S , 2 5 In these equations, t  is defined as turbulent dynamic viscosity, k represents turbulent kinetic energy,  is turbulent dissipation, k G represents the generation of turbulence kinetic energy due to the mean velocity gradients. b G is the generation of turbulence kinetic energy due to buoyancy. M Y represents the contribution of the fluctuating dilatation in compressible turbulence to the overall dissipation rate.  Figure S1). Figure S1. Flow profile at the larynx at pick inhalation of different mesh sizes.

Convergence test for the Mesh in the bronchi model
We meshed the bronchial model with 5,6.8,7.5 million elements, and compared between the velocity profile at the centreline and at the first bifurcation (see Figure S2 and Figure S3).

Ellipsoid geometric representation using Euler angles
We begin by reviewing the geometrical definitions of the particle orientation. First, let us build two additional auxiliary coordinate systems. In Figure S4, a schematic diagram of an ellipsoid is presented; we define three coordinate systems, the first (x,y,z), marked green, is the global lab coordinate system. The system (x',y',z'), marked black, is the particle coordinate system, its center coinciding with the particle center of mass, and the z' axis is aligned with the ellipsoid's major axis. Finally, (x'',y'',z''), marked blue, is a coordinate system centered at the particle center of mass, but parallel to the lab coordinate system ("Euler Angles --from Wolfram MathWorld," n.d.). The transformation matrix between the systems (x',y',z') and (x'',y'',z'') is marked as 'A', such that for a vector ' v or a matrix ' M : The 'A' matrix is defined by Euler angles and Euler quaternions (Lin Tian, Ahmadi, Wang, & Hopke, 2012), as cos cos cos sin sin cos sin cos cos sin sin sin sin cos cos sin cos sin sin cos cos cos cos sin sin sin sin cos cos where: we mark the vector N (orange) to be the intersection of the two planes, (x-y) with the plane (x'-y');  is the angle between x and N;  is the angle between z to z';  is the angle between x' and N In Figure S5, we have 3 examples illustrating this idea: Figure S5. Examples of 3 simple cases to illustrate Euler angle.

Euler Quaternions
The use of Euler angles directly in the solution of particle's motion, will lead to singularities (Fan & Ahmadi, 1995) at 0,

 =
. Therefore, to solve this problem we use Euler quaternions to visually describe the particle position and orientation, as defined in the following formula (Feng & Kleinstreuer, 2013): We note here that quaternions must always satisfy the condition Finally, in terms of the Euler quaternions, the rotation matrix 'A' now reads: For the completeness of this section, we develop a simple expression to convert the orientation of the particle represented in Euler quaternions, to a simple direction vector representation. This conversion is particularly useful when visualizing calculated behavior of such particles. We imagine the particle to be a unit vector in coordinate system (x',y',z'): We remind that the z' direction is the direction of the ellipsoid's major axis. Using the transformation matrix, the rotated particle would then be Force Balance An ellipsoid particle, moving in a shear flow at a low particle Reynolds number is subjected to the gravitational force ( ), the hydrodynamic drag force ( ) and the lift force ( Subscripts p and f represent particle and fluid respectively. The mass an ellipsoid particle used in equation (14) is given by: Where p a is the particle's semi-minor axis (the semi-major axis will be marked as

Hydrodynamic Drag Force
The hydrodynamic drag force is defined as where f  is the viscosity of the fluid, and '' K is the resistance tensor accounting for a particle's shape. The resistance tensor for an ellipsoid in the particle coordinate system, ' K , is given by: Since the particle coordinate system is also the particle principal coordinate system, this matrix is diagonal. The transformation of this matrix to the lab global system (in which the force balance is solved) is In fact, the resistance tensor in the particle coordinate system ' K is constant, and the drag force is actually changing due to rotation of the particle, through the rotation matrix A. Therefore, the drag force takes the final form

Particle Lift Force
We assume that the particle is small enough that a linear shear flow can be locally approximated around the particle. The lift force is then given by where B is the transformation matrix of velocity gradients, added here to create a simpler to implement sum representation. This matrix is given by (Yu Feng, 2013)

Torque Balance
Previously, the importance of solving the orientation of the ellipsoid fiber in order to solve the force balance was described. Thus, it is of special importance to solve the torque equation of this particle as well. The Euler rotation equations, i.e. the balance of torques, reads I  I I  T  dt  d  I  I I  T  dt  d  I  I I  T We note that all values here are in the particle coordinate system (x',y',z'). ω is the angular velocity, and T is the torque. The principal mass moment of inertia ( )

Hydrodynamic Torque
Under the assumption that the particle is sufficiently small, and the flow around it can be approximated as a linear shear flow, the hydrodynamic torques read where D is the deformation tensor, and W is the spin tensor In Eq. (29), the fluid velocity gradient matrix is in the reference frame of the particle, i.e. it needs to be rotated from the lab frame. This can be done using Eq. (7): Once the torques are solved (Eq. (26)) we are able to integrate the angular velocity and calculate the change in the orientation and find the new quaternions:

Particles Characteristics
Table S1 contain a detailed description of the particles geometry and flow characteristics, including d -is the diameter of a sphere with the equivalent volume, a and b are the minor and major diameters of the fiber, is the fiber equivalent volume diameter given by Shapiro and Goldenberg (1993): = 2 √ ln( +√ 2 −1) , 0 is the relaxation time calculated as 0 = 2 18 , Where particle density = 1000 3 ⁄ and air dynamic viscosity = 1.26 − 6 ⁄ , and is the stokes number of fibers, calculated as = 0 0 , where 0 is the maximal velocity during peak inhalation and is the average alveoli diameter. As the AR increase decrease, as well as the relaxation time and the Stk number of the fibers, meaning they are more convected, reach the flow velocity and follow stream lines.

Converging the number of particles simulated in the upper airways model
We analysed the DE results of different number of randomized particles. According to the results presented in Figure S6, negligible difference was found between simulating 1,500 particles of the same group and simulating 3,000 particles of the same group. Figure S6. DE results of different number of particles, shown in each plot. Negligible difference was found between 1,500 and 3,000 particles results.

Feng & Kleinstreuer(2013)
We compared our upper airways deposition efficiencies to the work of Feng & Kleinstreuer [18] and found a match between the current study results of fibers of AR=3 and AR=10 with = 1.83 . Our results did not match the results of AR=30 (see Figure S8).

Validation of the bronchi model
We compared our bronchial airways deposition efficiencies to the work of Koullapis et al. [23] and found a good match between the current study results of spheres to the range of particle sizes simulated (see Figure S9).