Salbutamol Transport and Deposition in the Upper and Lower Airway with Different Devices in Cats: A Computational Fluid Dynamics Approach

Simple Summary Administration of inhaled salbutamol via metered-dose inhalers can effectively treat bronchoconstriction. Different devices are used for the delivery of this drug in cats, either in the hospital or at home, for long-term treatment. Effective drug administration may depend on the drug delivery device as well as patient cooperation. By using non-invasive computational fluid dynamics techniques, the impact of these devices on the deposition and transport of salbutamol particles in the cat airways was simulated and assessed. The results confirm a variable drug distribution depending on the device used. The percentage of particles reaching the lung was reduced when using spacers and increased when applied directly into an endotracheal tube. Abstract Pressurized metered-dose inhalers (pMDI) with or without spacers are commonly used for the treatment of feline inflammatory airway disease. During traditional airways treatments, a substantial amount of drugs are wasted upstream of their target. To study the efficiency of commonly used devices in the transport of inhaled salbutamol, different computational models based on two healthy adult client-owned cats were developed. Computed tomographic images from one cat were used to generate a three-dimensional geometry, and two masks (spherical and conical shapes) and two spacers (10 and 20 cm) completed the models. A second cat was used to generate a second model having an endotracheal tube (ETT) with and without the same spacers. Airflow, droplet spray transport, and deposition were simulated and studied using computational fluid dynamics techniques. Four regions were evaluated: device, upper airways, primary bronchi, and downstream lower airways/parenchyma (“lung”). Regardless of the model, most salbutamol is deposited in devices and/or upper airways. In general, particles reaching the lung varied between 5.8 and 25.8%. Compared with the first model, pMDI application through the ETT with or without a spacer had significantly higher percentages of particles reaching the lung (p = 0.006).


Introduction
Asthma and chronic bronchitis are common inflammatory airway disorders in cats. In people with asthma and chronic obstructive pulmonary disease (COPD; most commonly rat [19], and monkey [20] without clear clinical veterinary applications or just for improving human medicine and health. Each type and formulation of the drug and delivery device must be optimized for the species and respiratory disorder of interest. The deposition and transport of salbutamol varies and depends on particle size, inspiratory flow velocity, and devices design used, among others. CFD technique was used to understand pMDIs performance, reduce drug waste, and improve the delivery of the aerosol particles to the lung considering different inhalers and devices [14]. In addition, the airflow varies throughout the respiratory system. In the upper airways, the flow is turbulent (Reynolds number > 2000), whereas as the airways narrow in the lower zones, the flow becomes laminar. This will also influence the behavior of the particles [7]. The CFD tool allows these flow regimes to be simulated.
In this study using healthy cats, our first aim was to assess the efficiency of pMDIs of salbutamol using different spacers, masks, and an ETT comparing particle deposition and transport in the respiratory tract of cats with numerical models. For reaching this goal, CT-based CFD models were created for determining particle deposition and/or transport in four regions: device, upper airways (nose to the trachea), primary bronchi, and downstream lower airways/parenchyma (herein referred to as "lung").
Our hypothesis was that by using CFD to model pMDI delivery of salbutamol in healthy cats, most particle deposition would take place in the drug delivery devices and the upper airways. We also hypothesized that particle transport would be enhanced in cats with an endotracheal tube compared with a spacer chamber or a pre-oxygenation mask. This study would help to better understand conditions for optimal pMDI drug delivery in healthy cats as a steppingstone to the investigation of cats with inflammatory airway disease such as asthma and chronic bronchitis.

Animals
Two client-owned cats presented to the teaching veterinary hospital (Centre hospitalier universitaire vétérinaire d'Alfort-Chuva) between September 2019 and July 2020 with informed owner consent were enrolled in this study. One cat, without intubation, was under sedation for a head CT for otitis media and a second cat was intubated (ETT 4 ID) for a front limb CT which included the thorax. The weights of the cats enrolled were 4.25 and 4.5 kg, respectively. CT images were selected for geometrical reconstruction and further simulations.

Ethics
All procedures were conducted as part of normal veterinary clinical practice with owner's consent form (Art. R242-48, Ordre National de Vétérinaire) and approved by the Clinical Research Ethical Committee of the École Nationale Vétérinaire d'Alfort (ENVA), France (Number: 2020-05-30).

Pressurised Metered Dose Inhaler (pMDI) and Devices Selection, CT Images, and Design
Nine different possible combinations of devices were geometrically reconstructed based on CT images or manufactured data (see Figure 1).
Using the CT study from the cats, a geometrical reconstruction was performed and then used for further analysis. The selected devices were two types of pre-oxygenation masks with spherical ( Figure 1a) and conical (Figure 1d) shapes. Then, a spacer of 10 cm length and 4 cm diameter and a spacer of 20 cm length and same diameter were connected to the two previous masks (Figure 1b,c,e,f). For these CT images, no cats were required. Finally, a CT image of an intubated cat was selected with the aim of posteriorly creating a numerical model for simulating a direct administration of a pMDI through an endotracheal tube (Figure 1g) or through the two spacers (Figure 1h,i). Pressurized metered-dose inhaler (pMDI) and device models for computational fluid dynamics simulations: (a) pMDI connected to a spherical preoxygenation mask, (b) pMDI connected to a 10 cm spacer and a spherical preoxygenation mask, (c) pMDI connected to a 20 cm spacer and a spherical preoxygenation mask, (d) pMDI connected to a conical preoxygenation mask, (e) pMDI connected to a 10 cm spacer and a conical preoxygenation mask, (f) pMDI connected to a 20 cm spacer and a conical preoxygenation mask, (g) endotracheal tube (ETT) 4 mm ID, trachea, and main bronchus, (h) pMDI connected to a 10 cm spacer and an ETT 4 mm ID, trachea, and main bronchus, (i) pMDI connected to a 20 cm spacer and an ETT 4 mm ID, trachea, and main bronchus. Computational CT-based models with Ansys: a-f, non-intubated cat model; g-i, intubated cat.

Sedation, Multidetector Computed Tomography (MDCT) Protocol and Image Analysis
Cats fasted for 12 h prior to anesthesia but had free access to water until premedication was administered. Sedation or general anesthesia was adapted to each cat. Cats were spontaneously breathing and placed in sternal recumbency with the head elevated and the neck fully extended. Non-contrast-enhanced MDCT examinations were carried out using a 64-detector-row CT system (Brilliance 64; Philips, Amsterdam, The Netherlands) before the programmed CT study. The CT images were obtained using a matrix of 768 × 768, tube voltage of 120 kV, tube current of 196 mA, and a display field of view of 35 cm and a pitch of 0.5. Images for the study were acquired first from the nostrils to the most caudal border of the lungs [21]. One-millimeter-thick images were reconstructed using a high-resolution algorithm. At the conclusion of the CT examination, the cat was supervised until completely recovered.
Images were reviewed using commercial medical imaging software, DICOM (Digital Imaging and Communications in Medicine, DICOM) viewer (Horos v.1.1.7., 64-bit, Horos TM , Brooklyn, NY, USA) using a lung window (window width (WW): 1600; window level (WL): −550). A board-certified radiologist performed the evaluation of the images for major abnormalities involving the airways.

Geometrical Reconstruction and Numerical Discretization
The DICOM files derived from the CT studies were imported into the image-based geometry reconstruction software (MIMICS, Materialise Software, Leuven, Belgium; Figure 2).
Manual reconstruction of the upper airway (defined as the nasal cavity, nasopharynx, oropharynx, larynx, and trachea) and primary bronchi (also called principal or mainstem bronchi) geometry was conducted for the non-intubated cat.
By means of the commercial software (Ansys IcemCFD, v.20, Ansys Inc., Canonsburg, PA, USA), device volumes and airways were filled with tetrahedral elements that composed the computational mesh in which the airflow governing equations, droplet formation, and droplet trajectories were solved. In Figure 3, a grid is shown for the spacer of 10 cm and the spherical mask using sections and corresponding images.
Finally, two different masks of conical and spherical shape and two lengths of sacers as described previously were connected to the cat nose of the reconstructed model. The devices were created and attached to the patient-specific cat model in the commercial software package Rhinoceros (release 5, Robert McNeel and Associates, Seattle, WA, USA). The final model of the non-intubated cat is depicted in Figure 3.  Representation of the numerical discretization in a non-intubated feline CT-based model of pressurized metered-dose inhaler (pMDI) drug delivery. The airways volumes are filled with three-dimensional elements (tetrahedrons). In this example, a pMDI is connected to a 10 cm spacer and a conical preoxygenation mask.
The same procedure was used for the manual reconstruction of the intubated cat. In this case, only the trachea and the first generations of the respiratory tract were segmented, and the tracheal tube and the inhaler were separately generated and attached in Rhinoceros. The models are represented in Figure 1. In this case, the geometry is much less complex with respect to that of the other cat as it starts directly from the trachea bypassing all the nasalto-laryngeal airways tract. Thus, the numerical grid contains far fewer elements resulting in less expensive computational costs. The volume was filled with three-dimensional elements (tetrahedrons; see Table 1), and the total number of elements varied depending on the presence or absence of spacers and the type of spacer. The number of tetrahedrons of the cat geometry within the non-intubated cats was approximately 18 million. The number of elements for the 10 cm and 20 cm length spacers ranged from 3 to 6 million, respectively. The number of elements for the spherical and conical masks ranged from 2.5 to 5 million, respectively. Prior to final computations, a mesh independence study was carried out in order to assess the dependence of the results on the grid size. Details of the study are provided in Appendix A.

CFD Analysis
Once the devices and airways volumes were filled with tetrahedrons, the resultant numerical grids of the cats were imported into a simulation software package (Ansys CFX, v.20, Ansys Inc., Canonsburg, PA, USA). This software solves the Navier-Stokes equations that describe flow motion in different conditions within the geometrical grids using numerical algorithms. In particular, the Ansys CFX software adopts the finite volume method. The exact mathematical formulation and the solving algorithms used by Ansys CFX are provided in the software manual (Ansys, 2020).
The peak inspiratory/expiratory flow of the non-intubated cat was 110 mL/s [22] and of the intubated cat, 30 mL/s. The flow of the intubated cat was obtained by means of a thermal anemometry ("hot wire") flow sensor (Dräger Julian Anaesthesia Machine, Lübeck, Germany). The peak inspiratory flow was imposed at the top of the pMDI (see Figure 4). The flow was considered turbulent, and the k-ω model was used. An initial turbulence intensity value of 5% was adopted.
A respiratory cycle of 3 s (1 s inspiration, 2 s expiration) was selected and 7 respiratory cycles for a total of 21 s using a time step of 0.01 s were computed. Fluid properties, air density of 1.185 kg/m 3 , and viscosity of 1.83 10 -5 Pa·s [21] were used. The pMDI modeling included a dose of 100 µg for puff, an initial droplet diameter of 10 µm, a velocity of 150 m/s, and a spray angle of 20 • [23,24]. The pMDI puff was applied just before the first inhalation waveform starts in order to simulate clinical conditions.
The pMDI used in the model consists of a canister connected to a metering valve capable of producing required dosages of 25 to 100 µg. The user controls the release of the drug through the actuator-nozzle that generates a measured liquid. We used a nozzle orifice of 0.5 mm diameter and located at the center line of the canister in the center of the mask or the spacer (Figure 4). The injected particles were then tracked through the geometric regions until they found one of the three specific conditions: (1) they collided and were trapped on the mask or spacer tube and/or on the airway walls, (2) they escaped from the domain through one of the outlet geometry faces, or (3) they continued in suspension in the flow. Suspension means that particles are still in a hold-up, and they could subsequently be exhaled or impact an airway wall in further breathings. The numerical models allow computing flow velocity and structures inside each cat's airways and devices. The structure of the flow was depicted using 3D streamlines, while the flow intensity was represented using a heatmap.
All computations were performed on an Intel 9 workstation and parallelized on 8 processors. The computational costs of every single simulation are given in Table 1. Further details of the particle modeling and independent study are described in Appendices A and B [25][26][27].

Data Analysis
Airflow streamlines (peak inspiratory flow structures) and droplet deposition were represented using heat maps and described qualitatively for the various conditions in both models. Percentages of particles deposited to the devices, to the upper airways, and transported to the lung as previously defined for all conditions in each model, were reported. Further data analysis used IBM ® SPSS ® Statistics software version 23 (Chicago, IL, USA). A Shapiro-Wilk test was used to test for the normal distribution of percentages of particles deposited to the aforementioned three regions and in suspension for further comparisons. A t-test for independent values was used to compare the means and SD of the delivery percentages of particulates to devices, upper airways, to the lung, and in suspension, between the two models (i.e., the non-intubated cat and intubated cat). Values of p < 0.05 were considered statistically significant with a 95% of confidence interval (CI).

Results
Using CT scans from two cats, nine different combinations of the model devices were simulated. Airflow streamlines, which represent the flow direction depicted with the intensity of the velocity (red = high velocity, dark blue = low velocity) of all simulation models, are represented in Figure 5 at peak flow during inspiration. In the case of the spherical mask, the presence of the spacer seemed to have less influence on the flow. The recirculating airflow structure appeared similar independent of the presence of the spacer (Figure 5a-c). This may have been due to the geometry of the mask. The spherical mask was shorter than the conical mask hence producing a flow recirculation for its proximity to the cat nose in all the cases. In general, the presence of the spacers greatly influenced the flow in the masks, as visible in Figure 6, especially in the case of the conical mask (Figure 5d-f). The flow recirculation appearing in the conical mask in the presence of the spacer was largely reduced when the spacer was not attached. In fact, the air jet generated at the entrance of the conical mask was generated at the entrance of the two spacers when these were connected, and hence the flow velocity reduced, and its structure accordingly modifies. The flow structure inside the spacer appeared similar within the models and independent of its length. After the skewed inlet of the inhaler ring, the airflow velocity in the spacer experimented with a sudden expansion at the spacer entrance, generating a local recirculation near the spacer walls. The airflow velocity increased again at the entrance of the mask and in the oral cavity and further accelerated in the trachea because of the constriction in the larynx. At this stage, the flow was completely turbulent.
The distribution of inhaled airflow was heterogeneous within the cat nasal cavity due to the complex geometry of the airways. The initial high-speed flow in these regions was slowed down in the dorsal meatus ( Figure 5), progressively becoming a low-speed flow in the ventral regions (blue regions of Figure 5). Then, the flow was directed to the nasopharynx and laryngeal regions, where local acceleration occurs due to the physiological constriction represented by the larynx. Finally, the flow moved forward and divided to the lower airways through the bifurcations and generations from the carina.
The airflow in the intubated cat was more homogeneous in comparison to that of the non-intubated cat (Figure 5g-i). As the inhalation proceeded from the endotracheal tube, the airflow entered directly to the trachea and distributed to the airways with a lower velocity as the cat was anesthetized. No influence on the flow patterns was visible when the spacers were added to the tracheal tube. However, important variations were expected regarding the particle inhalation.
The total of particles in percentages (%) reaching the lung after seven respiratory cycles airways or remained in suspension after the simulation are represented in Figure 6 (I:E ratio 1:2, inspiratory time 1 s; see Table 2), attached to the device and/or upper. Particle deposition percentage was computed as the ratio of deposited and injected particles times 100. These percentages are summarized for the inhaler, the spacers, and the masks in Table 3. They represent the amount of drug which is not able to reach the lungs. Those particles deposited on the muzzle or in the upper airways (nasal turbines, oropharynx, larynx, and trachea) are represented in Table 4. In the Supplementary Materials, the unsteady behavior of the droplets injected though the salbutamol pMDI are visualized during the breathing cycles (see Videos S1-S4).
The particle deposition in the devices was 2.6 times higher when a conical mask was used compared with a spherical mask alone. The drug deposition was twice when using a 20 cm spacer with a conical mask and 0.7 times higher when a 10 cm spacer was used compared to when the spherical mask was attached to them. However, the deposition was 1.6 times higher on the muzzle when using the spherical mask with respect to the conical one. Finally, 66.41% of the particles were deposited before reaching the upper airways in the conical mask and 48.96% in the spherical mask when they were used alone. When using these masks with a spacer of 10 cm that, as commented, is specific to cats, the particle deposition in the regions before the upper airways was found to be 22.65% and 26.85% (CM and SM, respectively). When a human pediatric spacer of 20 cm was used, these percentages were 41.35% (CM) and 30.15% (SM).
Mean and standard deviation (±SD) of the percentage of particles deposited into the devices, upper airways, and reaching the lungs (exit the computational model) when using a mask or an ETT with or without a spacer are represented in Table 5. Table 5. Percentage of particle distribution and deposition comparing when using a mask or an ETT with or without a spacer.

Particles
Group Mean and standard deviation (±SD) of the percentage of particles deposited into the devices, upper airways, and reaching the lung when a spacer was added or not are represented in Table 6.
In Figure 7, the salbutamol droplets deposition patterns are represented on the geometrical model surfaces. The pattern provided a qualitative picture of the regions affected by high or low salbutamol concentration. Visible droplets tend to deposit in non-uniform ways. There was a tendency of particles to deposit especially in the upper mask region (see Figure 7a,b). When the spacer was introduced into the model, the deposition on the mask was slightly reduced. A considerable amount of salbutamol was concentrated on the spacer surface (see Figure 7c-f). Salbutamol deposition on the muzzle was enhanced when using the spherical mask. In any case, salbutamol loss due to deposition was not significantly different between particles attached to the devices or upper airways (p= 0.428 and p = 0.089, respectively).

Discussion
The present study aimed to assess the transport, distribution, and deposition of salbutamol 10 µm particles in the upper and lower airways of cats by means of the CFD technique. The pMDI application through the ETT, with or without a spacer, had significantly higher percentages of particles reaching the lung compared with the non-intubated model. In humans, similar studies were used to improve aerosol delivery techniques with similar aims, and the use of the CFD technique is widely known [15,23,28,29]. However, compared with humans, the anatomy of cats' upper airways is more complex and thus challenging for this technique. Studies assessing deposition of aerosol clinical therapies in small animals are limited and based on the use of Scintigraphy in dogs [30]. The use of pertechnetate ( 99m Tc) delivered via a spacer and face mask apparatus in conscious cats has been published [31], but this technique employs the use of radioisotope and gamma-rays, which is more invasive and may be toxic for the patient. Less invasively, Leemans et al., 2009, compared the duration of action and the efficacy of ß2-agonists in cats with induced bronchoconstriction using barometric whole-body plethysmography (BWBP) [32]. The advantage of this numerical tool is that it is a non-invasive technique.
We simulated the most commonly prescribed bronchodilator at home and used in clinics, salbutamol, with a concentration of 100 µg per puff (Ventolin, GlaxoSmithKline plc (GSK), Brentford, UK). It is the same product with exactly the same dosage as the one used in humans (average weight of 70 kg). It is unknown what plasma level is required to achieve effective bronchodilatation in cats, and toxicity trials with aerosolized formulations in cats are missing. Leemans et al., 2009, are the only group that studied the effects and justified a dose for salbutamol pMDI. Their study recommended a single dose of 100 µg with a peak effect of around 15 min post administration and inducing an antispasmodic effect in the airways for up to 4 h. Two-or four-fold increments of drug dosages slightly improve bronchodilatory effects. Similar results were reported in humans [33,34]. Consequently, in this study, a single dose of 100 µg was simulated, with a fixed peak inspiratory flow for the non-intubated and intubated model. Treatment recommendations varied between the time of application and/or the number of breaths after activation of the pMDI [34]. Some companies recommend allowing the cat to breathe through the mask for 10 to 15 s or to take 5 to 10 breaths after activation of the pMDI (Trudell Medical International, London, ON, Canada). In this study, seven respiratory cycles for a total of 21 s (physiological I:E ratio of 1:2) were imposed since we wanted to allow enough time to make sure that all the particles were distributed into the model. Even so, in most cases, there was still a percentage of particles that remained in suspension, except for the direct pMDI endotracheal tube (ETT) administration. In this case, after less than one complete respiratory cycle (after 1.1 s), all particles were already distributed. It is also worth mentioning that for the conical mask, only 2% of the particles remained in suspension. Hence with ETT and conical mask alone, particles were transported faster.
The peak inspiratory/expiratory flow was lower in the simulated intubated cat (30 mL/s) compared with the non-intubated cat (110 mL/s). Significant greater percentages reaching the lung were found in this study when applying pMDI directly in the ETT (19.7-25.8%) compared with the other techniques (5.75-16.2%; p = 0.006), as preconized in our hypothesis. In humans, the serum concentrations of an inhaled drug are greater in intubated patients as there is no oropharyngeal deposition, no enteral absorption, and best lung deposition [35]. One-third of the total aerosol output per pMDI puff is expected in smaller ETT size (4 mm ID as was used in the current study) compared with >5 mm ID ETT [36]. Side effects after ß2-agonist administration directly in the ETT during general anesthesia were described in horses. Specifically, sinus tachycardia, premature ventricular complexes or/and hypotension [37], and sweating [38] were observed after the administration of aerosolized salbutamol. As there are no reports in cats, the exact plasma concentration at which bronchodilation is expected versus at which side effects will appear is unknown, and thus this practice should be employed with caution. Side effects in cats may include tremors, central nervous system (CNS) excitement, vomiting, mydriasis, and/or dizziness [39].
Flow velocities and patterns affect particle behavior and hence their deposition and transport [40]. It was suggested that the use of a spacer might reduce particle velocity from the pMDI and may allow a better distribution [8]. Other studies found that an add-on spacer with a pMDI may result in a considerable reduction in the number of respirable particles available to the patient [41]. The clinical significance of this effect is not well established in humans [7,42]. Although a spacer was specifically designed for cats (AeroKat ® , Trudell Medical International, London, ON, Canada), some veterinarians and/or clients still choose the option of pediatric/baby spacers. These spacers are usually longer and may or may not have a one-way valve (valved holding chambers, VHC). The dose delivered may vary considerably between spacers, and this has to be considered when changing from one spacer to another [43]. For this reason, we have included both types of spacers in the current study. Our results suggest that no benefit is produced by the use of one or the other type of spacer ( Figure 5), independently of the type of mask. The percentage of particles reaching the lung was similar in the presence and in the absence of a spacer (p = 0.757). This may be due to the flow structures inside the masks, which revealed high recirculation. The Global Initiative for Asthma (GINA) recommends the use of pMDI and a dedicated spacer/VHC with facemask for children aged 4 years and younger and a pMDI plus a dedicated spacer/VHC with mouthpiece for children aged between 4 and 6 years. In comparison, for adult humans, no mask is added to the spacer to maximize the spacer's function of reducing particle velocity yet improving particle transport [23]. No guidelines or consensus around this topic exists in veterinary medicine. In this study, more particles were deposited in the masks when the spacer was not added, probably as the spray nozzle was nearer (see Tables 3 and 4). In any case, more particles remained in suspension when using a spacer (p = 0.037). That is, even after 21 s, there were still particles that had not reached any structure or the lung. It is possible that when using a spacer, more time is needed for the salbutamol to fully reach the lungs. However, the total number of particles transported to the lungs was similar between all groups.
Particles reaching the lung using pMDI in humans vary considerably between studies. From 11 to 14%, drug deposition was reported in ambulatory patients using radiolabeled aerosols [34]. Other studies suggest that less than 10% of the inhaled bronchodilator finally reaches the target area when using pMDI alone [44,45], and from 10% to 38% with the use of the pMDI plus a spacer [46]. In a computational study using a salbutamol pMDI with a spacer attached to a human upper airway model and a flow velocity of 30 L/minute, 52.9% of particles traveled to the lung [23]. In the present study, the percentage of particles reaching the lung was higher than 10% in most of the cases and varied between 5.75 and 25.8% depending on the device used. Most of the studies considered particles sizes from 1 to 10 µm [29,40]. In humans, it was reported that the deposition also depends on the particle size [29]. In particular, while nanoparticles deposition tends to increase when the particle size decreases (<0.1 µm), deposition of microparticles tends to increase when the size increases (>1 µm).
Reduced lung deposition percentages will be expected in clinical situations, where cats do not tolerate excessive handling, particularly when respiratory distress is present. This is in comparison to the current study's models in which the devices fit perfectly with the cat's anatomy, and all the particles are utilized. The deposition fractions summarized in Table 3 suggest that spherical masks may enhance particle transport into the lungs compared to conical masks. However, in the presence of the spherical mask, the deposition in the rostral area, on the muzzle, and in the nasal passage is higher than in the presence of the conical mask (see Table 4). Additional studies, including a greater number of cat models, are needed for a better comparison between masks.
Other factors that may modify particle lung deposition when using pMDI aerosol in ventilated humans are the ventilator mode, settings, and circuit [47], and synchronization of pMDI application with breaths, among other factors [48]. In addition, the therapeutic response may be influenced by the patient's airway anatomy and disease severity [49]. Further studies are needed to determine all these factors in cats.
The results that can be extracted from the computational analyses have to be considered carefully. The computational models, in general, are affected by unavailable limitations. The main limitation of this study was that an ideal and perfect fitting mask with or without a spacer was simulated, where 100% of the particles are sent to the cat airways. This may not reflect reality. Another limitation concerned the limited number of animals. Further studies including additional animals are necessary. However, the aim of the study was to compare different devices and not different cat anatomies. For this reason, only two cats were used to develop the models. Finally, only particles of 10 µm size were simulated. Smaller particle diameters should also be investigated as their behavior strongly depends on size, as was demonstrated in humans [29].

Conclusions
Using CT scans to develop feline models and CFD to investigate salbutamol transport and deposition, this study determined that most of the particles deposited in the devices and/or upper airways before reaching the lung. Direct administration of pMDI through the ETT has the largest and fastest lung deposition compared with the rest of the devices. The use of a pediatric spacer versus a specifically designed spacer for cats did not impact droplets transport. Further studies using a larger number of cats, patient-specific airflows, and different particle sizes are warranted to investigate which delivery devices may have better and safer performance.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/ani11082431/s1, Video S1: Unsteady simulation of the particle deposition and transport within the non-intubated cat model with spherical mask. Video S2: Unsteady simulation of the particle deposition and transport within the non-intubated cat model with conical mask. Video S3: Unsteady simulation of the particle deposition and transport within the non-intubated cat model with spherical mask and spacer of 10 cm length. Video S4: Unsteady simulation of the particle deposition and transport within the non-intubated cat model with conical mask and spacer of 10 cm length.

Data Availability Statement:
The data generated during the current study are not publicly available because they are part of a national research project. However, some data could be available from the corresponding author on reasonable request.

Acknowledgments:
The authors acknowledge the statistical support provided by V. Herrería Bustillo. The support of the Institute of Health Carlos III (ISCIII) through the CIBER-BBN initiative is highly appreciated.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A Appendix A.1. Governing Equations
The fluid dynamics model presented in this study is based on the work presented by Kleinstreuer and co-workers [23] and Inthavong and co-workers [24] in humans, and it is well known and widely used for numerical simulations of particle transport and deposition in the human airways. For this reason, here we provide only a summary of the main features of the model, and additionally, validation with the literature data is presented.
The fluid dynamics of the cat breathing subjected to the turbulence and specific boundary conditions were solved simultaneously with the droplet formation and particle trajectories. The latter was solved using Newton's second law governed by drag and gravity as the major forces for dispersed droplets in low-density gas flow [23]. For solving the viscous and incompressible airflow in the cat airways, the well-known Reynolds averaged Navier-Stokes equations were used.
where v is the averaged fluid velocity vector, ρ is the gas density, and p is modified pressure, v eff is the effective kinematic viscosity where v is gas kinematic viscosity, and v t is turbulence kinematic viscosity). Between the kinematic and the turbulent kinematic viscosities exists the relation v t = k/ω. The turbulence was modeled using the Wilcox k-ω model [26]. This model predicts the turbulence using two partial differential equations for the variable k and ω, where k represents the turbulence kinetic energy and ω is the rate of dissipation of the turbulence kinetic energy. Both equations are provided below: where the term P k that represents the turbulent production can be obtained as: P k = υ t ∇v· ∇v + (∇v) T and α = 5/9, β = 0.075, β = 0.09 and σ k = σ ω = 2 are turbulent constants [26].

Appendix A.2. Model of the Salbutamol Spray
The droplet spray was simulated with gas-liquid interactive Enhanced Taylor Analogy Breakup (ETAB) model [23]. With this model, a particle exposed to gas flow is subjected to large drag force by the surrounding gas if a velocity difference between particle and gas exists. The drag force produces droplet formations accordingly to the Weber number (ratio of the inertia force to surface tension): where ρ is the air density, v is the relative velocity between air and particle, r is the particle radius, and σ is the coefficient of liquid surface tension. This non-dimensional ratio connects droplet deformation and the liquid surface tension. Particle trajectories are determined by CFX using Newton's second law in the presence of drag force and gravity. This equation can be expressed as: where v p is the particle velocity, m p is the particle mass, rho is the gas density, g is the gravitational force, c Dp is the drag coefficient, and d p is the particle diameter. The drag coefficient c Dp can be obtained with the relation: where C slip is the Cunningham correction factor and c D can be obtained as follows: where Re p is the particle Reynolds number given by: where v is the kinematic viscosity of the gas. Particles of 10 µm were injected through the nozzle orifice with the aforementioned velocity and the described spray properties and modeling.

Appendix A.3. Mesh and Particle Independence Study
Due to the complex geometry of the cat's upper airways, unstructured tetrahedral computational grids were created using the commercial software Ansys ICEM CFD, Version 20.0 (ANSYS Inc., Canonsburg, PA, USA). Tetrahedral elements filled the 3D airways volumes with a variable number and dimensions of elements. For this reason, to ensure independent results from the selected mesh element size, a mesh independence study was carried out. For this study, we considered the biggest geometry volume that is the cat model with the conical mask and the spacer of 20 cm. The model volume was discretized with an increasing number of elements, and by using each generated mesh, a CFD solution was obtained. We tested several different grids from 10 to 50 million elements. After the analysis of the velocity profiles at different geometrical cross-sections, it was clearly seen that from a grid size of around 30 million elements, the results were grid-independent.
The final model contains about 15 million tetrahedrons in the cat airways, about 6 million tetrahedrons in the spacer, about 5 million tetrahedrons in the mask, and about 4 million for the remaining model objects. As the mesh independence study estimates the necessary regional size, the other geometries used the same discretization technique. Of course, depending on the presence or absence of the spacer, and its length, and on the mask and its type, the global grid size changed. In other words, all geometries were previously subdivided into regions in which the value of the tetrahedron can be specified. For obtaining all the grids in all the analyzed cases, we used the same values obtained in the independent grid for each region of the considered models. As the mesh independence study was carried out on the biggest geometry (the spacer of 20 cm and the conical mask are the biggest devices), the grid sizes of the rest of the models were smaller. In particular, it is worth reporting that the grid size of the spherical mask was about 2.5 million tetrahedrons, and the spacer of 10 cm was composed of about 3 million tetrahedrons. The element sizes per each model that were used for the presented simulations are reported in Table 1.
A particle number-independence study was also conducted in order to verify that 2000 injected particles were statistically sufficient for describing the salbutamol transport and deposition. Since every single particle of 2000 injected is representative of a certain number of particles, it should be proved that the selected number correctly describes the physical phenomenon in the study. The geometry selected with this aim is the same used for the mesh independence study (spacer of 20 cm and conical mask). The CFD simulation was repeated using the same grid size but releasing a different number of particles. Particles exiting the geometries and depositing on the different regions of the airways were systematically collected, increasing the injected number to 200,000. No relevant changes in percentage were observed.

Validation of the Numerical Methodology with Literature Data
To validate the result obtained with the presented simulations, the computational methodology used for creating the in silico models needs to be validated. With this aim, and because no information is available in cats, spherical microparticle deposition was studied in a human upper airways geometry and compared with experimental and computational literature data (see below). For this, we used the typical idealized "path lung model" [26] widely known and used in the human computational aerosol technique. The comparison, at different breathing rates, for different particle diameters, and for the mentioned model are shown in Figure A1. In particular, particle diameters of 1, 5, and 10 µm were injected into the model through the mouth at light, normal, and heavy conditions (i.e., 15, 30, and 60 L/min). Steady inspiratory flows were computed, and particles deposited or exited from the model were collected. The results show that the model set-up is correct as it is capable of reproducing the results presented by Zhang et al. [29]. Additionally, the human upper airways in the presence of the pMDI and of the aerosol spray were simulated. The deposition fractions within the upper airways and the fraction of particles reaching the lung were compared with the results obtained by Yousefi et al. [24]. The latter presents a model of the airflow, droplet spray transport, and deposition in a pMDI attached to a human upper airway model, considering de Ventolin as a propellant. The model described in the present work includes a human oral cavity described in the previous Section B1 and a trachea truncated at its bottom. A pMDI model, as described by Kleinstreuer et al. [23], was attached to the mouth or to the model entrance (see Figure A2).
The following boundary conditions were given to the model: an inspiratory flow rate of Q = 30 L/min and a turbulence intensity value of 5% for the turbulence kinetic energy was applied to the canister ring (see Figure A2); a pressure condition was applied to the outlet (bottom section of the trachea); a perfect absorption to the model walls for the particles. At the inhaler nozzle, a droplet initial diameter of 10 µm was assigned with a spray velocity of 110 m/s and a spray angle of 10 • [24]. Finally, the used Ventolin density was 1230 kg/m 3 with an actuation dose of 100 µg [24]. Again, deposited particles were collected at the oral tract (mouth + soft palate + pharynx + larynx + trachea), and particles traveling to the lung were computed. The comparison with the data presented by Yousefi et al. demonstrates that again the methodology used in the present work is correct, and the model is capable of reproducing the literature results. Figure A2. Computational human model with human upper airway tract and pMDI used for validation (a) and regional deposition fractions (DF; in %) computed and displayed in comparison with the literature data (Yousefi et al., 2017) for inspiratory flows Q of 30 L/min (b). In (c), the velocity distribution inside the oral airways is represented by means of velocity streamlines.