Mathematical modeling and numerical simulation of surfactant delivery within a physical model of the neonatal trachea for different aerosol characteristics

ABSTRACT

Surfactant aerosol delivery in conjunction with a noninvasive respiratory support holds the potential to treat neonatal respiratory distress syndrome in a safe manner. The objective of the present study was to gain knowledge in order to optimize the geometry of an intracorporeal inhalation catheter and improve surfactant aerosol delivery effectiveness in neonates. Initially, a mathematical model capable of predicting the aerosol flow generated by this inhalation catheter within a physical model of the neonatal trachea was implemented and validated. Subsequently, a numerical study was performed to analyze the effect of the aerosol liquid droplet size and mass flow rate on surfactant delivery and on the required aerosolization time period. Experimental validation of the mathematical model showed a close prediction of the air axial velocity at the distal end of the physical model, with an absolute error between 0.01 and 0.15 m/s. Furthermore, an admissible absolute error between 0.2 and 2 µm was attained in the prediction of the aerosol mean aerodynamic diameter and mass median aerodynamic diameter in this region. The numerical study highlighted the beneficial effects of generating an intracorporeal aerosol with a mass median aerodynamic diameter higher than 4 µm and a surfactant mass flow rate above 8.93 mg/s in order to obtain effective surfactant delivery in neonates with minimal airway manipulation.

Copyright © 2017 American Association for Aerosol Research

EDITOR:

• Warren Finlay

Introduction

The current therapy for treating neonatal respiratory distress syndrome (RDS) is surfactant replacement therapy, which consists of administering surfactant via intratracheal instillation (Sweet et al. 2013). Nevertheless, this treatment entails several risks due to the necessity of tracheal intubation and mechanical ventilation, leading to the possibility of lung injury and chronic lung disease (Jobe and Bancalari 2001). To date, several alternatives have been tested in order to avoid these undesirable outcomes. For instance, there is a trend toward using noninvasive respiratory support as the primary treatment for RDS instead of giving a surfactant dose (Rojas-Reyes et al. 2012). However, this can be prejudicial if the patient does not respond to the respiratory support and still needs a surfactant dose, since surfactant should be administered as early as possible to maximize its therapeutic efficacy (Carlton et al. 1995, Seidner et al. 1995). Another alternative is to use the INSURE technique, where the patient is intubated during a short period, just to administer the surfactant intratracheally, and then the patient is rapidly extubated to noninvasive respiratory support. In this case, however, it is still necessary to intubate the infant (Dani et al. 2010). Recent advances in surfactant delivery have demonstrated that the application of surfactant via a thin catheter to spontaneously breathing preterm infants receiving continuous positive airway pressure (CPAP) reduces the need for mechanical ventilation (Göpel et al. 2011). This technique requires highly skilled personal and careful manipulation of the airways.

Aerosolized surfactant delivery together with a noninvasive respiratory support is a long pursued alternative to treat neonatal RDS. Several in vivo studies have shown that aerosolized administration of the drug leads to equal or even further improvement of the lung function (Lewis et al. 1991, Schermuly et al. 2000). In addition, some clinical studies have administered surfactant aerosols to neonates (Jorch et al. 1997, Arroe et al. 1998, Berggren et al. 2000, Finer et al. 2010). Unfortunately, only the study conducted by Jorch et al. yielded positive results. Aerosol delivery to neonates is challenging in comparison to adults. The characteristics of the neonatal respiratory system, such as their narrow airways and high respiratory frequency, as well as the high viscosity of the surfactant, and the extracorporeal generation of aerosols with adult devices, account for a low surfactant lung deposition in neonates, less than 1% of the atomized mass (Kölher et al. 2008).

Nevertheless, few works have studied drug delivery devices for neonates. Syedain et al. (2015) evaluated an intrapulmonary aerosol generation device for surfactant delivery in neonates. Holbrook et al. (2015) developed several nebulization devices to administer pharmaceutical aerosols to ventilated infants. However, these devices require the intubation of the infant. For this reason, it is necessary to develop a device intended for neonates in order to deliver drug in a noninvasive and effective way. We have previously investigated the feasibility of using an intracorporeal inhalation catheter (IC) to deliver surfactant and perfluorocarbon aerosol in vitro. Moreover, we also evaluated the feasibility of using this IC to generate intracorporeal aerosols in neonates, positioning the distal tip of the catheter at the beginning of the trachea in a neonatal physical model. We concluded that the design of this device should be optimized to improve the fraction of surfactant deposition in infant lungs. For this purpose, performing a computational fluid dynamics (CFD) study would be a helpful methodology (Longest and Holbrook 2012, Carrigy et al. 2014). For instance, the numerical study of Longest et al. (2013) demonstrated the feasibility of optimizing the design of a nebulizer, by means of CFD simulations, to improve aerosol delivery during nasal high flow therapy in adults. A mathematical model that is able to predict aerosolized drug delivery within the neonatal upper conducting airways would allow us to find adequate aerosol characteristics for improving surfactant delivery in neonates and to modify the design of the device to generate this aerosol.

The aim of the present study is to validate the feasibility of numerically predicting the aerosol flow within a physical model of the neonatal trachea, as well as to perform a numerical study in order to analyze the influence of aerosol characteristics on surfactant delivery within this neonatal physical model and to provide guidelines for optimizing the design of the inhalation catheter.

Mathematical modeling

Two mathematical models were defined by means of CFD techniques to predict the aerosol flow at two different scenarios. On the one hand, the flow beyond the distal end of an intracorporeal IC (AeroProbe, Trudell Medical International, London, ON, Canada) consisting of six outer air lumens and one inner liquid lumen (Murgia et al. 2011) was evaluated. On the other hand, the flow within a physical model of the neonatal trachea (before the branching) was analyzed. It must be emphasized that the distal tip of the IC was placed at the beginning of the trachea of the physical model, representing the position just above the vocal cords, for a noninvasive aerosol delivery condition. Moreover, no external air flow was assumed beyond the flow of the IC, i.e., the respiratory air flow of the neonate was not included. In addition, simplified two-dimensional mathematical models, specifically axisymmetric models, were proposed.

Flow domain and boundary conditions

Two different axisymmetric flow domains were proposed. Domain 1 represents the flow beyond the distal tip of the IC (Figure 1a). It is composed of the internal geometry of the last 2 mm of the IC and the downstream region beyond the IC's distal end. Domain 2 represents the flow within and beyond the trachea (denominated as G0) of the physical model described in the online supplementary information [SI], Figure 1b. Domain 2 is composed of the internal region of the IC, the internal region of the physical model and the downstream region beyond the physical model. IC dimensions were obtained from the experimental measurements provided in the SI, whereas the dimensions of the physical model of the neonatal trachea were obtained from Rozanek and Roubik (2007). A physical model of the neonatal trachea was manufactured in order to carry out this research work because medical experts highlight that once surfactant reaches the trachea, it would be spread towards the distal part of the respiratory system (Enhörning and Robertson 1972).

Figure 1. Computational geometry and boundary conditions for (a) Domain 1 and (b) Domain 2.

An air mass flow rate of 1.13.10⁻⁵ kg/s, corresponding to a driving pressure of 4 bar (Goikoetxea et al. 2014), a stagnation temperature of 293K, a turbulent intensity of 7% and a hydraulic diameter of 0.06 mm were fixed as the inlet boundary condition. The isothermal boundary condition was defined for the walls, with a static temperature of 293K. Besides, trapping condition was considered for the liquid droplets impacting on the wall. Axial-symmetry was applied to the central axis of the IC. A mean null gauge pressure was indicated for downstream boundaries, with an ambient static temperature of 293K, a turbulent viscosity ratio of 2, and a backflow turbulent intensity of 1%.

As far as the discrete phase is concerned, ten partial injections were simultaneously released in order to represent the droplet size distribution obtained from experimental measurements (Goikoetxea et al. 2014). A uniform aerodynamic diameter (D_a) and the corresponding mass flow rate (m) were indicated for each partial injection (Table 1). The compounds Perfluorodecalin (PFD, C₁₀F₁₈, F2 Chemicals Ltd., Lancashire, UK; density 1.95 g/ml; kinematic viscosity 10.5 cSt) and Poractant-alfa (Curosurf®, Chiesi Farmaceutici S. P. A., Parma, Italy; density 1 g/ml; kinematic viscosity 2.7 cSt) were tested. A total droplet mass flow rate of 8.31 mg/s for PFD and 4.47 mg/s for surfactant (Poractant-alfa) were obtained from experimental measurements made at 4 bar driving pressure (Goikoetxea et al. 2014). The released PFD aerosol had a mass median aerodynamic diameter (MMAD) of 2.29 µm and a mean aerodynamic diameter, D_a, of 1.69 µm, whereas the surfactant aerosol had a MMAD of 3.76 µm and a mean D_a of 1.62 µm. Liquid droplets were released 2 mm beyond the distal tip to reproduce exact droplet size measurement conditions. The droplet release axial velocity (ranging between 139 and 280 m/s) and the static temperature (ranging between 255 and 287K) were determined by taking into account the numerical results of the converged air flow of the present study at the injection point for each turbulence model and flow domain. These high velocity and low temperature values of the air flow at the injection point are due to the fact that the flow reaches nearly sonic conditions at the distal end of the inhalation catheter (Potter and Wiggert 2002).

Table 1. Uniform droplet aerodynamic diameter (D_a) and mass flow rate (m) for each partial injection, for PFD and surfactant aerosols.

	PFD		Surfactant
Injection N°	Da (µm)	m (mg/s)	Da (µm)	m (mg/s)
1	0.524	0.009	0.730	0.004
2	0.785	0.006	1.09	0.009
3	1.17	0.10	1.63	0.05
4	1.76	0.62	2.46	0.24
5	2.64	2.18	3.67	0.71
6	3.94	2.89	5.50	1.42
7	5.91	1.56	8.24	1.01
8	8.84	0.48	12.33	0.56
9	13.24	0.35	18.43	0.40
10	19.82	0.12	27.66	0.06

Data obtained from experimental measurements.

Continuous and discrete phase modeling

The flow of the IC consisted of a primary air flow plus a discrete phase. Due to the compressibility and the turbulence regime of the air flow, the Reynolds and Favre averaged Navier–Stokes equations detailed by Wilcox (1998) were used to solve the continuous flow (Wilcox 1998). Turbulence was modeled by means of four different Eddy Viscosity Models (EVM), i.e., the Standard and Realizable k-ϵ turbulence models as well as the Standard and Shear Stress Transport (SST) k-ω turbulence models (Wilcox 1998). To properly predict flow characteristics close to solid walls where viscous effects are prevalent, a near-wall model of the turbulence is required. In the EVM based on the dissipation rate (ϵ), the two-layer approach was employed, where the low-Reynolds model of Wolfshtein (1969) is combined with the previously mentioned high-Reynolds models through the methodology proposed by Chen and Patel (1988). When the EVM was based on the specific dissipation rate (ω), the transport equation of ω was integrated through the viscous sublayer fixing the value of ω on the wall (Wilcox 1998).

A particle tracking method based on the Euler-Lagrange approach was employed to solve the discrete phase. The drag law for spherical droplets, based on the correlations of Morsi and Alexander (1972), was selected. The dispersion effect of turbulent fluctuations on the discrete phase was predicted by the eddy lifetime model of Gosman and Ioannides (1983).

The prediction of the collisions between parcels was based on the algorithm in O'Rourke (1981). Coalescence between droplet parcels as well as secondary break-up physical models were taken into account to predict the interaction between aerosol droplets. Due to the value of the collisional Weber number (We_c < 100), the Taylor Analogy Breakup (TAB) model was selected to predict the secondary break-up of the discrete phase (Taylor 1963).

Discretization

The spatial discretization of the mathematical model was carried out by generating a mesh of the flow domain using Gambit 2.4.6 software (ANSYS Inc.). After undertaking a mesh refinement study (Goikoetxea 2015), a mesh consisting of 820,140 rectangular elements was selected. Refinements in special regions, such as the symmetry axis, the near wall region and the region next to the distal tip of the IC, were made. Moreover, a temporal discretization was defined by means of a constant time-step of 1·10⁻⁴ s for both phases and a maximum number of 650 iterations per time-step for convergence. Discretization of the differential equations and boundary conditions was done by means of the Finite Volume Method. Finally, the commercial CFD code ANSYS Fluent 14.5 (ANSYS Inc.) was used to run the numerical simulations.

Numerical methods

Numerical simulations were initially started by solving the continuous phase flow in a steady state. After reaching a converged solution for the continuous phase, discrete phase injections were added and both phases were simulated in a transient state until they converged. A pressure-based solver was selected for the numerical simulation of the mathematical model. The pressure and velocity fields were solved by using a coupled strategy. A courant number between 50 and 200 was used. A least-squares procedure was employed for gradient discretization. Second order upwind schemes were used for pressure interpolation, density, momentum, turbulence and energy. A first order scheme was used for temporal discretization. The under-relaxation factor for the discrete phase source terms was set to 0.1.

A two-way coupling method was employed for the interaction between continuous and discrete phases. The automated combination of an implicit low order tracking scheme and a trapezoidal high order tracking scheme was selected for the discretization of the discrete phase equations. Node-based averaging was employed to smoothly distribute the effect of the discrete phase source terms in the continuous phase. A constant-number parcel release method was used, setting 400 particles per parcel for PFD injections and 100 particles per parcel for surfactant injections to achieve a good balance between the accuracy of the results and the computational cost.

Continuous phase calculations were considered converged when the axial velocity at the exit of the flow domain changed less than 0.1% each 1,000 iterations and the locally scaled residuals were below 1.10⁻⁵. After reaching a converged solution for the air flow and releasing the droplet injections, the locally scaled residual for the continuity equation was limited to below 1.10⁻³. Moreover, the discrete phase mass balance was checked to ensure that the predicted flow time was appropriate in order to start taking droplet samples. Several samples were taken for each set of data to confirm that the obtained results correctly represented the behavior of the converged aerosol flow.

Experimental validation

Both mathematical models were validated by comparing numerical and experimental results. The air axial velocity profile measurements detailed in the SI were taken by means of a hot-wire anemometer with associated software (StreamLine Constant Temperature Anemometer and StreamWare, Dantec Dynamics, Denmark) to validate air flow prediction. Numerical results beyond the distal tip of the IC reached an absolute error of between 0.1 and 2.9 m/s depending on the turbulence model, whereas numerical results beyond the G0 of the physical model reached an absolute error of between 0.01 and 0.15 m/s (Figure 2). The Realizable k-ϵ turbulence model yielded the best prediction of the experimental results. This agreement was correct for both scenarios, i.e., with and without the physical model.

Figure 2. Comparison of experimental and numerical air axial velocity profiles (a) 30 mm beyond the distal tip of the IC and (b) 5 mm beyond the exit of the G0 of the physical model, for a driving pressure of 4 bar.

In addition, droplet size measurements were taken by means of an aerodynamic particle sizer spectrometer with associated software (APS 3321 and Aerosol Instrument Manager, TSI, Brooklyn, NY, USA) to predict droplet size evolution (Goikoetxea et al. 2014). These measurements were taken after determining the optimal distance between the distal tip of the inhalation catheter and the suctioning nozzle of the equipment (35 mm for surfactant and 60 mm for PFD) to ensure that the suctioned number of particles was within the appropriate range for each compound in order to achieve an accurate measurement. Numerical results beyond the distal tip of the IC (Table 2) showed an absolute error of between 0.9 and 1.4 µm for the mean D_a and a value between 0.6 and 1 µm for the MMAD, depending on the aerosolized compound. The most accurate prediction of the aerosol flow beyond the IC was achieved with the Standard k-ϵ turbulence model. Regarding the numerical results beyond the G0 of the physical model, an absolute error of between 1.6 and 2 µm was obtained for the mean D_a and a value between 0.2 and 1.1 µm was obtained for the MMAD. For this scenario, the best flow prediction was attained by means of the Standard k-ω turbulence model.

Table 2. Experimental and numerical results of droplet size characterization 35 mm (surfactant) and 60 mm (PFD) beyond the distal tip of the IC and 10 mm beyond the G0 of the physical model.

	PFD				Surfactant
	Beyond IC		Beyond G0		Beyond IC		Beyond G0
	Exp.	CFD	Exp.	CFD	Exp.	CFD	Exp.	CFD
Mean Da (µm)	1.1	2.5	4.5	2.9	0.9	1.8	0.9	2.9
MMAD (µm)	4.4	3.4	4.9	6.0	6.7	6.1	10.7	10.5

Beyond IC: Domain 1. Beyond G0: Domain 2.

Numerical and experimental results of the cumulative mass distributions in function of the droplet size for both scenarios are illustrated in Figure 3. On the one hand, droplet size curves at the optimal distance from the distal tip of the IC (Domain 1) for each compound are shown. On the other hand, droplet size curves beyond the G0 (Domain 2) are presented. It must be pointed out that the “Droplet release” data correspond to the injected droplet size, which is equal for numerical and experimental data. From this figure, it can be observed that the complex physics involving the aerosol flow was predicted reasonably well regardless of the aerosolized compound and the scenario studied.

Figure 3. Comparison of experimental and numerical cumulative mass based droplet size distribution curves (a) 35 mm (surfactant) and 60 mm (PFD) beyond the distal tip of the IC, and (b) 10 mm beyond the exit of the G0 of the physical model. Beyond IC: Domain 1. Beyond G0: Domain 2.

Furthermore, information about the aerosolized compound delivery was obtained (Table 3). The numerical results of the aerosolized mass delivery were analyzed by studying three different quantities: (1) The mass that was deposited at the G0 walls of the physical model, which is referred to as mass deposited at the physical model, (2) the mass that escaped from the outlet boundaries of the flow domain, referred to as mass directed beyond the G0, and (3) the mass that had still not been deposited or had not yet escaped when the results were analyzed, referred to as mass in suspension. The highest mass fraction of the aerosolized compound (between 39.1 and 63.4%) was deposited within the G0 of the physical model. In addition, a mass fraction of between 10.4 and 23.5%, depending on the aerosolized compound, went beyond the G0 of the physical model. The rest of the aerosolized mass fraction (over 26%) stayed in suspension.

Table 3. Numerical results of the percentage distribution of the atomized mass.

Mass (%)	PFD	Surfactant
Deposited at the physical model	63.4	39.1
Directed beyond the G0	10.4	23.5
In suspension	26.2	37.4

The experimental validation confirmed that the flow predictions obtained by means of the implemented mathematical models closely match the experimental results of the aerosol flow beyond the IC and within the physical model of the neonatal trachea. It is known that the EVM, due to its nature, cannot reproduce the anisotropy of Reynolds stresses, such as that existing near the solid walls. When the turbulent dispersion in the movement of the discrete phase is considered within the context of EVM, it is necessary to correct the calculation of the continuous phase's turbulent fluctuation velocities to take that anisotropy into account. In duct flows, in which the solid walls play an important role and they are the main source of anisotropy, it has been shown that the isotropic hypothesis leads to an overprediction of the deposition velocities of small particles (Wang and James 1999, Matida et al. 2000). Several corrections have been made to consider the effect of the solid walls on the Reynolds stresses. Most of these corrections are functions of the dimensionless wall distance (y⁺) and have been calibrated in duct flows (Matida et al. 2004, Zhang et al. 2005). In this study, this correction has not been introduced in the mathematical model. Because of this, a bigger discrepancy would be expected between the experimental and simulated results of the particle size distribution in the case that considers the G0 than in the case that does not consider G0. However, noticeable differences in the experimental validation have not been found between cases, and therefore it seems that not taking this correction into account in the model does not play a crucial role in the prediction of the analyzed aerosol flow.

Numerical study

A numerical study implemented one of the validated mathematical models in order to test different surfactant aerosols and to obtain practical information about the influence of the aerosol characteristics on surfactant delivery, i.e. on droplet size change and surfactant deposition.

First, the influence of the aerosolized surfactant total mass flow rate was analyzed. Two alternatives to the original surfactant mass flow rate administered by the IC were numerically evaluated: at half and twice the original value, namely 2.23 and 8.93 mg/s, respectively. The results showed that the highest surfactant mass fraction obtained at the distal end of the G0 of the physical model corresponded to droplets with an aerodynamic diameter (D_a) of between 10 and 15 µm for the three cases (Figure 4a). Figure 4 shows mass based droplet size distribution at the end of the distal end of the G0 for all the hypothetical aerosol injections that were simulated. Figure 4a indicates that there was no clear relation between the mass flow rate variation and its influence on droplet size. Although the mean D_a of the aerosol at the distal end of the G0 decreased when a lower mass flow rate was administered, this did not happen in the case of the MMAD (Table 4). Nevertheless, in comparing the results corresponding to the two extreme cases (2.23 mg/s and 8.93 mg/s), it was clearly observed that a higher coalescence between droplets occurred when a higher mass flow rate was released.

Figure 4. Numerical results of mass based droplet size distribution curves at the distal end of the G0 of the physical model. (a) Influence of injected mass flow rate. (b) Influence of injected droplet size.

Table 4. Numerical results of droplet size characterization just at the distal end of the G0 of the physical model for different surfactant aerosols.

	Droplet mass flow rate			Droplet size
	4.47 mg/s^*	2.23 mg/s	8.93 mg/s	Original^*	2 μm monodisperse
Mean Da (µm)	2.7	2.5	3.8	2.7	3.0
MMAD (µm)	9.8	8.7	9.4	9.8	7.1

* Original surfactant aerosol generated by the IC.

In terms of surfactant delivery, this increase in coalescence between droplets yielded a lower mass quantity in suspension (between 8.4 and 18.6% lower) (Table 5). As a consequence, a higher mass fraction was either deposited in the wall of the physical model or directed beyond the G0. Regarding the aerosol cone angle, it increased by an average of 14.5% when a higher surfactant mass flow rate was aerosolized (Figure 5).

Table 5. Numerical results of the percentage distribution of the atomized mass for different surfactant aerosols.

	Droplet mass flow rate			Droplet size
Mass (%)	4.47 mg/s^*	2.23 mg/s	8.93 mg/s	Original^*	2 μm monodisperse
Deposited at the physical model	39.1	25.0	46.7	39.1	45.3
Directed beyond the G0	23.5	19.0	24.3	23.5	10.3
In suspension	37.4	56.0	29.0	37.4	44.4

* Original surfactant aerosol generated by the IC.

Figure 5. Numerical results of the behavior of surfactant aerosols within the physical model and the aerosol cone angle (α) beyond the G0 for different mass flow rates and droplet sizes.

Secondly, the influence of surfactant aerosol droplet size was evaluated. Taking into account that the MMAD of the original surfactant aerosol produced by the IC was 3.76 µm and its mean D_a was 1.62 µm, a monodisperse aerosol with a MMAD of 2 µm was numerically studied. Although the mean D_a slightly increased when a monodisperse aerosol of 2 µm was administered, the surfactant mass fraction obtained at the distal end of the G0 of the physical model was transported by smaller droplets, mainly by droplets with a D_a of between 5 and 10 µm (Figure 4b). Consequently, the aerosol MMAD obtained at the distal end of the G0 resulted in a 27.5% decrease. With regard to surfactant delivery, the mass fraction directed beyond the G0 when a monodisperse aerosol of 2 µm was released was reduced by 56.2%. In this case, the value of the aerosol cone angle was decreased by 60.2%.

Discussion

The practical application of the main findings of the present work must be highlighted. As a result of this study, a mathematical model that is able to predict the aerosol flow within a physical model of the neonatal trachea was achieved, and information about the influence of several aerosol characteristics on surfactant delivery was gained. These findings, which will be discussed below, lead to design guidelines that can optimize the geometry of the inhalation catheter and improve drug delivery effectiveness in neonates, which will eventually make it possible to treat neonatal RDS by applying noninvasive respiratory support during intracorporeal yet noninvasive surfactant aerosol delivery.

A primary outcome of the current study, namely that it is feasible for the IC to generate surfactant and PFC aerosols with an adequate size in a context of intracorporeal aerosol generation, is in accordance with previous studies (Murgia et al. 2011; Goikoetxea et al. 2014). Focusing on the surfactant results, an aerosol with a mean D_a of 2.9 µm and a MMAD of 10.5 µm was obtained beyond the G0 of the physical model. In addition, a significant percentage of the aerosolized surfactant (23.5%) went directly beyond the G0 of the physical model. What is more, the aerosolized mass fraction that finally goes towards the neonatal distal airways could reach a higher amount (over 60%) given that medical experts support the idea that if a tracheal deposition of the surfactant is achieved, the surfactant will spread correctly towards the distal airways (Enhörning and Robertson 1972). This amount could be even higher, because the mass that stays in suspension (37.4%) would either be deposited at the G0 wall or directed beyond the G0 during the therapy.

Experimental validation of these results confirmed that a mathematical model able to numerically predict the evolution and behavior of intracorporeal aerosols in neonatal upper conducting airways was obtained. This model permits to numerically predict the turbulent and compressible character of the air flow produced by the IC, obtaining a correct prediction of the axial velocity field and its interaction with the aerosol droplets. It also allows predicting droplet size change within and beyond the physical model, confirming that the results of the complex physics of the flow were modeled correctly enough. Based on these results, the suitability of CFD calculations to advance in this research line was corroborated.

One of the aims of the current work was to run a numerical study that implemented the validated mathematical model in order to analyze the influence of several aerosol characteristics on intracorporeal aerosol delivery and, therefore, to provide design optimization guidelines for the IC. On the one hand, the numerical study indicated that a higher coalescence between droplets was obtained for the highest aerosolized mass flow rate (8.93 mg/s). Droplets tended to be either deposited on the wall or directed beyond the domain more quickly when a higher coalescence between droplets occurred. This could be desirable because when the surfactant mass amount that stays in suspension is low, so is the mass that would be exhaled. In conclusion, in order to improve surfactant delivery in neonates, it would be recommendable to generate a surfactant aerosol with a high droplet mass flow rate; that is, an aerosol that results in a high coalescence between droplets.

As far as droplet size is concerned, it mainly affected the MMAD value at the distal end of the G0, in addition to the surfactant delivery. When the MMAD of the generated aerosol was lowered from 3.76 to 2 µm, the MMAD of the aerosol obtained at the distal end of the G0 had a value that was 27.5% lower. However, the surfactant mass fraction that was directed towards the distal respiratory airways decreased by 11.2%. On the other hand, an important aspect that should also be taken into account is the aerosolization time period required to administer a 200 mg/kg Curosurf dose to treat an infant with RDS (Valls-i-Soler et al. 1997; Speer et al. 2013). For a 1,000 g preterm this would be equivalent to administering a 2.5 ml dose of Curosurf, which implies an aerosolization time period of 5, 10 or 19 min, depending on the aerosolized surfactant mass flow rates that were analyzed in this study, i.e., 8.93, 4.47 and 2.23 mg/s, respectively. Considering that surfactant should be administered to the infant as early as possible to maximize its beneficial effects (Seidner et al. 1995), it would be preferable to generate an aerosol that requires a relatively short aerosolization time period in order to administer the necessary surfactant dose. However, it is necessary to meet the balance between delivering the dose quickly but also gradually, in order to prevent rapid physiological changes in the neonate and avoid airway obstruction. Rey-Santano et al. (2013) concluded that a surfactant aerosol delivery time period of between 10 and 15 min is beneficial for lung response and achieves a better hemodynamic response than the rapid bolus in preterm lambs with RDS. It must be highlighted that, according to our findings, we can modulate the aerosolization parameters in order to achieve a gradual surfactant delivery. In fact, based on the current optimization results, a 5 min aerosolization time period is proposed for the therapy, but it could be modulated to adjust the treatment as needed.

In conclusion, it would be recommendable to generate an aerosol with a MMAD higher than 4 µm and a total mass flow rate higher than 8.93 mg/s in order to increase the surfactant mass fraction that is directed towards the distal respiratory airways and to reduce the time required to administer the dose. However, this study has limitations that must be acknowledged. On the one hand, considering the difficulties and drawbacks of obtaining a real image of the trachea of a neonate, a simplified geometry of the neonatal trachea was implemented. In addition, the respiratory cycle of the infant as well as the humidity and temperature conditions of the respiratory airways were not accounted for in the mathematical model. Further, the possibility of a wall-film formation due to the impact of droplets on the wall of the G0 of the physical model and its consequent sliding to the downstream region were not considered in the numerical model. Bearing these limitations in mind, the knowledge gained from the present study may be beneficial for the optimization of the IC design and, consequently, for obtaining an effective and gradual surfactant delivery method with minimal manipulation of the airway that could be combined with a noninvasive respiratory support, such as nasopharyngeal CPAP.

Acknowledgments

The authors would like to acknowledge the support of the Cátedra Fundación Antonio Aranzábal-Universidad de Navarra.