The open-source software package of OpenFOAM (www.openfoam.org) is used to numerically solve the particle-laden flow equations inside the nasal airway. The Pressure-Implicit with Splitting of Operators (PISO) method combined with Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) algorithm (pimpleFoam solver) is used for the solution of the Eulerian phase and the Lagrangian Intermediate library is combined with the pimpleFoam solver to enable the solution of the secondary phase in the gas phase (Patankar and Spalding, 1972; Patankar, 1980; Syrakos et al., 2017). The pimpleFoam solver is suitable for LES simulations with the local dynamic k-equation sub-grid scales (SGS) model. The local dynamic k-equation SGS model is used since it covers the laminar-transitional-turbulent region so that it is very suitable for the simulation of the present work which is predominantly laminar with possible transitional and turbulent flow regimes. It is important to note that in the local dynamic k-equation SGS model, C k and C ε is dynamically chosen depending on the time and space.

The diffusion terms are discretized with the least-squares scheme, which is second-order accurate on general unstructured meshes (Martínez et al., 2015). Since upwind-dominated schemes could lead to numerical dissipations in orders of SGS viscosity dissipation and consequently decrease the precision of the numerical simulations (Sweby, 1984), the convection term in the momentum equation requires a proper discretization scheme. In the present work, “Gauss filteredLinear” scheme, which is a low-dissipation second-order central differencing method, is implemented for LES. The convection term of other (scalar) quantities is discretized by the second-order “Gauss limitedLinear” scheme based on Total Variation Diminishing (TVD) (Moukalled et al., 2016). The time derivatives are discretized by the Second-Order Upwind Euler (SOUE) scheme (Wang et al., 2009).

Based on Figure 1, the nasal cavity is divided into three main parts: the wall, inlet, and outlet. The nosepieces are extended to represent the fully developed condition (Shi et al., 2007). The sinusoidal and steady air flow profiles are implemented for the velocity at the inlet by a Dirichlet boundary conditions. The air velocity for both inlet conditions represents an average flow rate of 4 L/min. The velocity profile of the gas phase is described as U = U 0 + U 0 s i n ( 2 π ω t ) , (16) where ω is 45 Hz. OpenFOAM utilizes a phase by face-to-face particle tracking procedure to simulate the motion of the particles (Macpherson et al., 2009). The analytical method is used to integrate particle Lagrangian equations (Minier et al., 2003). Furthermore, the properties of water are assigned to the particles since most of the drugs are aqueous solutions and as a result their densities are close to water density (Edwards et al., 1997). Utilizing the so-called parcel method, 4 million computational particles (parcels) with a total injection mass of 2.5 mg are uniformly and randomly entering through the right nosepiece for a duration of 0.5 s to represent realistic conditions as encountered for nebulizers such as PARI SINUS in contrast to the majority of the research in the literature which only considers a very short injection duration in range of a single time step. Moreover, two-way coupling between particles and gas phase is considered due to the high volume fraction of the particles of 7.5 E-5, which lies above the threshold of 1 E-6 (Ferrante and Elghobashi, 2007; Crowe et al., 2011). Since two-way coupling between air momentum and particles is considered, a particle-related source term is added to the gas phase momentum equation. In order to avoid statistical error, a large number of particles are required to be injected (Tofighian et al., 2019). In the literature typically short injections are studied and one-way coupling is implemented for simulation of drug delivery with different devices (Ahookhosh et al., 2021; Taheri et al., 2021); however, for the simulation of the aerosol delivery with nebulizers which have continuous cloud injection, one-way coupling between both phases does not require a particle-related source term in the momentum equation. The introduction of the source term requires larger statistical power, which requires an injection of a few million particles as compared to the previously used 10,000–50,000 particles (Dastan et al., 2014; Ghahramani et al., 2014; Keeler et al., 2016). Moreover, due to deposition of very few particles in the olfactory region, an injection with large number of particles is required to achieve independency of DE in olfactory region from number of injected particles and consequently to avoid statistical errors in prediction of olfactory epithelium deposition.

Here, 1 µm particles are considered with an average flow rate of 4 L/min. Moreover, pulsatile aerosol delivery is compared with non-pulsatile (steady) inlet airflow carrying particles with a diameter of 1 µm. As mentioned above, this particle size is used due to its relatively low relaxation time and high diffusivity which means that particles follow the streamlines (low impaction) and reach the olfactory region where they can deposit due to diffusion more than larger particles. This concept also rationalizes the choice of a relatively low flow rate of 4 L/min, since this reduces impaction in bends like vestibule and nasal valve and enhances the residence time and thus diffusive deposition of the particles in the olfactory region. The particles were distributed randomly and uniformly across the inlet plane and the velocities of the particles at the inlet are equal to the instantaneous velocity of the gas phase at the specific time and location of each particle. It is assumed that the particles stick to the wall once they touch it. The particles are followed until they hit the wall or reach the outlet and escape the computational domain. Although aerosol injection ceases after 0.5 s, the computational simulation is continued for 1.98 s to ensure that all particles have either deposited or exited the domain. The simulations are performed using 256 CPUs and the time steps were variably chosen based on the Courant number criteria to be lower than 1 which led to time steps in a range of 10⁻⁵ s to 10⁻⁶ s.

The validity of the present numerical solver was assessed in two steps. In our previous publications (Farnoud et al., 2017; Farnoud et al., 2020a; Farnoud et al., 2020b), the DE in 90° bends was compared to the experimental data of Pui et al. (1987) and other research in the literature (Breuer et al., 2006; Nicolaou and Zaki, 2016; Inthavong, 2019) to ensure that the solver is able to capture the DE of particles correctly. This benchmark is a popular approach to validate the capability of a numerical solver to predict particle deposition in curvature geometries.

In the next step, our model was compared to experimental and numerical data obtained for other nasal passages. Since the majority of the experimental and numerical studies on nasal aerosol delivery have focused on bilateral aerosol delivery, we validated our computational model with experimental DE data from a different nasal passage with bilateral drug delivery condition. Thus, for the current nasal passage geometry with bilateral air flow (drug delivery), open soft palate and delivery through both nostrils, DE is calculated as a function of the characteristic impaction parameter ( IP = d 2 Q ) and compared with previous nasal drug delivery studies in the literature (Cheng et al., 2001; Kelly et al., 2004; Shi et al., 2007; Hsu and Chuang, 2012). Moreover, in our previous study (Farnoud et al., 2020a) the pressure drop in nasal passage was compared to pressure drop in other patient-specific nasal airway models which showed a good agreement between the studies (Kelly et al., 2004; Schroeter et al., 2011). For mono-disperse particles with diameters of 1–30 µm (100,000 particles launched for each case) carried by three airflow rates of 4, 15, 30 L/min entering through both nostrils (regular inhalation), the calculated DEs agree well with experimental data and other computational studies in the literature that used different nasal airways and similar boundary and initial conditions (see Figure 2) (Cheng et al., 2001; Kelly et al., 2004; Shi et al., 2007; Hsu and Chuang, 2012). It is also important to note that for 1 µm at 4 L/min the DE is with 1.25% rather low indicating relatively low particle deposition due to impaction and diffusion. This is consistent with the well-known nasal incapability limit of ca. 10 µm, i.e., for normal breathing only particles smaller than ca. 10 µm can reach the lung (Hinds, 1999).

The result of large eddy simulation (LES) is highly dependent on the quality and resolution of the computational grid. To determine the resolution and quality of a LES simulation various criteria, typically referred to as “LES indices of Resolution Quality” (LES_IQk), have been introduced. One LES_IQk is defined as the ratio of resolved turbulent kinetic energy (TKE) to the total TKE, and can be formulated as: M = k res / ( k res + k sgs ) (17) where k res s is the resolved turbulent kinetic energy and k sgs is the sub-grid scale part of the turbulent kinetic energy. Celik et al. (2005) recommended a range of 0.77–0.85 for M to ensure that the grid resolution is adequate. Pope (2004) also considered grid resolution as sufficient if M is larger than 0.8. It is evident from Figure 3 that M is higher than 0.99 in the entire computational domain, indicating that the grid is finer than typically required. This hyperfine grid resolution is chosen since the fine-grained and complex geometry of the nasal cavity and relatively high flow rates require such a fine mesh.

In the present study, the effect of bi-directional pulsatile inlet airflow on the deposition pattern of the particles in the nasal airway, and specifically the olfactory epithelium region, is investigated. Since the inlet airflow profile directly affects dispersion and deposition of the particles, pulsating and non-pulsating airflow conditions are implemented at the inlet to assess the effect of the pulsation on the airflow pattern inside the nasal airway and specially the olfactory region.

Figure 4 shows the streamlines of the air with non-pulsatile (constant) inlet flow rate of 4 L/min which are colored by normal velocity (normalized to average air speed at the inlet, i.e., right nostril). The simulations are performed until t = 1.98 s and Figure 4 represents the results at this time. Different regions of the nasal anatomy are shown in the bottom panel of Figure 4 to enable the further discussion based on the regions in the nasal cavity. At the right vestibule, a recirculation zone is observed and swirling flows occur as the air passes through the nasal valve. At the upper parts of the nasal valve (Figure 4. Top left panel) a swirling flow is observed which moves downwards to the inferior right meatus. The air flow from the nasal valve is divided into three meatuses and some streamlines move towards the top of the nasal cavity and pass through the olfactory region with relatively high velocity. At the nasopharynx the airflows from the three meatuses mix and make a 180° turn while entering the left cavity. At the left cavity, the flow moves from the nasopharynx towards the left nostril and uniform streamlines are observed in the entire main nasal airway except for the right vestibule where a recirculation zone is observed. In both right (entering) and left (exiting) vestibules recirculation zones are observed. The air velocity at the left vestibule increases since the upstream airflow encounters a sudden contraction in the nasal valve region.

Similarly, Figure 5 shows the streamlines of the normal velocity in the nasal cavity and olfactory region at T/2 of the last air flow pulse for the pulsatile inlet airflow with average flow rate of 4 L/min and frequency of 45 Hz (oscillation period: T = 22.2 ms). A complex airflow pattern including recirculation zones and swirling flows are observed in the anterior parts of both right and left main nasal airway due to periodical pulses in the airflow. For both inlet airflow conditions, the maximum velocity is observed in the anterior region and is close to 1.4 times the average inlet velocity. In contrast, streamlines with higher velocities pass through the right olfactory region when the non-pulsatile airflow is implemented at the inlet. The pulsatile airflow at the inlet avoids the impingement of the airflow to the top of the cavity after the nasal valve, which results in streamlines with much lower velocity passing the olfactory region.

Since the streamlines only show the flow pattern and not the amount of the air reaching to each region, the airflow fraction in the right and left olfactory regions are calculated. For this purpose, a new variable, the concentration of a passive scalar tracer (C), is introduced in the OpenFOAM code. C is initially 0 in the nasal cavity but changes to 1 when the drug-laden air flow reaches the specific computational cell. Analogous to a gaseous tracer, once all of the initially present air in the nasal cavity is replaced by the inhaled air, C = 1 throughout the nasal cavity. Hence, C reflects the fractional replacement of the initially present air with the inhaled air (drug-laden air fraction). The variable C is calculated in the whole domain including non-convectively ventilated paranasal sinuses and the olfactory region.

Figure 6 shows the contour plots of the variable C and the streamlines of the air flow (at the end of inhalation period t = 1.98 s) colored by variable C in different regions of the nasal cavity specifically in the olfactory region for non-pulsatile inlet airflow condition. It is noteworthy that C is larger than zero even in some regions of the non-ventilated paranasal sinuses (Figure 6, top), which is due to diffusional (not convective) transport of inhaled air into these regions. Moreover, Figure 6 (bottom) depicts the iso-volumes (regions that represent a specific value of a parameter) of variable C for values more than 0.99 which means 99% of these regions are filled with inhaled air indicating that these regions can be readily reached by the inhaled air either via convection or via diffusion. As even the lowest iso-volumes of C = 97.67% in the left nasal valve and vestibule are close to 100%, gas exchange in the entire nasal cavity (including olfactory epithelium region) is highly effective for a 1.98 s inhalation (bottom panel and pie chart in Figure 6). This implies that also olfactorially relevant gas molecules carried by the inhaled air will be able to efficiently reach the olfactory epithelium, which is a prerequisite for the sense of smell to work properly. For fine particles with low relaxation time and high diffusivity, the same statement is valid.

There are considerable similarities and some differences between pulsatile flow (Figure 7) and non-pulsatile flow profile (Figure 6) when comparing the passive scaler tracer contours. Air exchange in the olfactory region and the entire nasal cavity is virtually identical and almost complete (>97%) for both cases (bottom panel and pie chart in Figures 6, 7). One notable difference is the lower ventilation (lower C values) in the left cavity specifically at the nasal valve and vestibule region (Figure 7, bottom right) when pulsatile airflow is used. Moreover, the pattern of the streamlines after the nasal valve is more straight for the non-pulsatile condition. However, in the pulsatile airflow simulations, the streamlines are oriented into all directions due to the disturbance caused by pulsatile airflow pattern.

A more detailed view on the olfactory region is presented in Figure 8. The fractional air flow in/through the left and right olfactory region at the end of the modeling period (t = 1.98 s) is 0.703 and 0.653% (of inhaled particle dose) independent of pulsation state, respectively, which is due to the slight differences in geometry and volume (left: 337 mm³; right: 331 mm³). Hence, from this perspective the conditions for relatively uniform left/right aerosol deposition in the olfactory epithelium region are favorable. This will be discussed in more detail in the next section.

Four million computational particles (parcels) with diameters of 1 µm are introduced for 0.5 s from the right nostril (at the velocity of the air at the nostril) and carried by non-pulsatile and pulsatile airflows (4 L/min) to the left nostril, which led to deposition of the particles in different regions of the cavity or to their exit through the left nostril. These bi-directional flow simulations are performed for a duration of 1.98 s to ensure that all particles have deposited or exited the nasal domain. Particle deposition patterns and DEs in the right and left nasal passages and olfactory regions are visualized and quantitatively evaluated at 1.98 s as presented below.

Figure 9 illustrates the spatial deposition pattern of 1 µm particles for either non-pulsatile (top) or pulsatile (bottom) inlet airflow in the nasal airway and olfactory region (highlighted). The main deposition occurs in the right nasal cavity and specifically in the regions with uneven surfaces which cause changes in the direction of the flow and consequently deposition of particles due to inertial impaction.

Quantitative analysis of our model reveals that total (1 µm) particle deposition in the nose is 50.9 and 34.4% for aerosol delivery with non-pulsatile and pulsatile bi-directional 4 L/min air flow, respectively (Figure 10). This is significantly larger than the corresponding value for bilateral inhalation (1.25%) indicating that during normal inhalation there is almost no impaction in the nasal cavity, while for bi-directional flow the sudden change in flow direction in the nasopharynx region induces significant impaction of 1 µm particles. In bilateral aerosol delivery (injection from both nostrils) particles mainly escape the nasal passage from the nasopharynx to the pharynx. However, during bi-directional delivery, particles enter from one nostril and exit from the other one which leads to substantial enhancement of DE in the nasal cavity.

For further validation of particle deposition modeling beyond what has already been presented above, it is instructive to compare our CFD-derived total (1 µm) particle deposition in the nose (Figure 10) with empirical in vitro and other modeling data on bi-directional nasal aerosol deposition. Kleven et al. (2005) experimentally and numerically found total nasal DEs of 63.1 ± 18.4 and 45.4% for 3.5 µm particles bi-directionally inhaled with a constant flow rate of 6 L/min. Moreover, Xi et al. (2017) reported DE of 24 ± 6% for bi-directional nasal drug delivery with a constant inlet airflow rate of 6 L/min, particle size of 3.2 µm. In another in vitro bilateral aerosol delivery study with a PARI SINUS nebulizer (pulsatile flow, frequency 44.5 Hz), particle MMD of 3.79 ± 0.03 µm, and airflow rate of 6 L/min, DE for the adult nasal model was reported as 80.74 ± 10.17% (Hosseini and Golshahi, 2019). The large range of DEs in the literature are presumably an indication of the strong dependence of nasal aerosol deposition on nasal geometry and experimental conditions. The DE values of this study are in the range of the DEs reported in previous experimental and numerical studies on pulsatile and bi-directional flow. Flow pulsation not only reduces total aerosol deposition by a factor of 1.48 (from 50.9 to 34.4%), it also mitigates the non-uniformity of right-to-left cavity deposition from 1.78-fold to 1.33. Thus, pulsation reduces impaction-induced deposition in both cavities, but more pronounced in the right cavity. Similar but even more pronounced trends were observed for the olfactory region where pulsation reduced olfactory particle deposition from 0.48 to 0.12% of the inhaled aerosol (3.9-fold reduction). Even more striking, the right-to-left olfactory DE ratio was reduced from 624 for non-pulsatile to 53.2 and pulsatile airflow, respectively, i.e., pulsation reduces the right-left non-uniformity by more than one order of magnitude, but it is still severe. It is noteworthy that pulsation increased particles deposited in the left olfactory region by a factor of 4.4, albeit still on a very low level (<0.01%). The main reason for this extremely large difference between left and right olfactory deposition can be seen in Figures 4, 5. For the non-pulsatile flow, streamlines impinge on the top part of the cavity after the nasal valve, and the air passes with relatively high-speed through the right olfactory region.

The main reason for this extremely large difference between left and right olfactory deposition can be seen in Figures 4, 5. For the non-pulsatile flow, streamlines impinge on the top part of the cavity after the nasal valve, and the air passes with relatively high-speed through the right olfactory region. In contrast, due to the oscillating nature of the pulsatile flow, streamlines after the nasal valve are moderated and mixed resulting in lower air speed than that of the non-pulsating flow when passing through the right olfactory region. Only for the right nasal cavity, one can observe some streamlines going into the olfactory region, but none are going into the left olfactory region. Since reduced air flow towards the epithelium results in reduced impaction of particles, pulsation reduces particle deposition for all regions where particle deposition is mostly affected by convective transport (impaction) of particles.

For drug delivery purposes, one is typically not only interested in DE but also in the aerosol (drug) dose deposited onto the epithelium, since the local drug dose per area determines the therapeutic efficacy of a drug. Typically, only a fraction of the inhaled aerosol is due to drug, but for simplicity we will assume a 100% drug loading for the following discussion. Moreover, we assumed that an aerosol (drug) concentration of 75 g/m³ is inhaled at a flow rate of 4 L/min for 1.98 s yielding a total inhaled drug dose of 2.5 mg.

The contour plots of Figure 11 illustrate the deposited drug dose represented as mass per area (epithelium) throughout the entire nasal cavity and specifically in the right and left olfactory regions when non-pulsating inlet airflow is used. The results show a patchy deposition pattern with the highest deposited dose regions in the anterior parts of the right nasal cavity specifically in the nasal valve and the upper parts after the valve. Less dose is observed in the nasal valve and vestibule of the left cavity compared to the right side. To gain a better understanding on the regional distribution of the deposited dose in the right and left olfactory regions, a close-up view of contour plots of dose is presented which shows a significantly higher dose in the right olfactory compared to the left olfactory epithelium region. The posterior nasal regions receive a much lower dose than the anterior region. Therefore, in separate figures these regions with corresponding dose values are shown. In the posterior region most of the deposition is observed where the flow direction changes or vortexes are formed, which is consistent with enhanced inertial impaction of the particles in regions where the air flow changes direction.

Figure 12 is analogous to Figure 11, except that pulsatile aerosol delivery (4 L/min, 45 Hz) is presented. Pulsation of the airflow has almost completely avoided the deposition of particles along the direction of the nosepiece. Therefore, the dose in the regions that reside along the direction of the nosepiece is much lower compared to the condition that non-pulsatile (steady) aerosol delivery is used. In contrast to non-pulsatile aerosol delivery, there are no high dose regions in the upper parts after the nasal valve. However, similar to the non-pulsating aerosol delivery, the regions with highest doses are observed at the vestibule and nasal valve. Analogous to Figure 11 the dose in regions after the nasal valve are separately shown with their corresponding values. Pulsation of the airflow reduces the extreme accumulation of particles in the regions that stand as obstacle perpendicular to the flow direction.

Analogous to Figure 10, average dose per area of epithelium in the nasal cavity is reduced due to pulsation. Figure 13 shows that there is almost no dose (<0.2 ng/mm²) deposited in the left olfactory region irrespective of pulsation conditions, while ca. 54–602-fold higher doses per area (9—36 ng/mm²) are delivered to the right olfactory region with pulsatile and constant flow conditions respectively. However, pulsatile airflow reduces the dose by a factor of ∼3.95 in the right olfactory region.

As seen from Figures 11, 12 and analogous to DE (Figure 10), pulsating flow reduces the right-left imbalance of deposited dose (per area) for both the entire cavity and the olfactory epithelium region. This is due to both reduced and enhanced aerosol deposition in the right and left half of the nose (both cavity and olfactory region), respectively. Albeit the DE in the right cavity (50.9 and 34.4%) is 68-fold and 162-fold higher than that in the right olfactory region for non-pulsating and pulsating flow, respectively (Figure 10), the average local dose (mass per area) in the right cavity is 1.05-fold and 2.5-fold higher than that in the right olfactory region for non-pulsating and pulsating conditions, respectively. This is mainly due to the ca. 64-fold lower area of the olfactory region as compared to nasal cavity. On the other hand, the local dose in the olfactory region is ca. 100-fold lower than that in the left cavity independent of pulsation conditions.

In summary, bidirectional inhalation of 1 µm aerosol provides relatively high nasal DE (>34.4%) and about the same average local (surface-specific) drug dose in the olfactory epithelium (ca. 36 ng/mm²) as observed in the nasal cavity for the right half of the nose (inhalation side). This is not the case for the left half of the nose, where the local aerosol dose in the olfactory regions is ca. 100-fold lower than in the nasal cavity. Hence, pulsatile flow did reduce the right-left imbalance in terms of DE, but not to levels where similar right-left local doses could be observed. Owing to the intricate nasal geometry the spatial aerosol deposition was highly non-uniform with very pronounced dose hotspots of up to 1,000 ng/mm² even in the olfactory epithelium region (Figures 11, 12). Bidirectional aerosol delivery substantially enhanced the total nasal deposition efficiency for 1 µm particles carried by steady inhalation rate of 4 L/min compared to the conventional bilateral aerosol delivery method [50.9 versus ca. 1% (Shi et al., 2007; Shang et al., 2015)]. Similarly, bidirectional aerosol delivery enhanced the olfactory deposition efficiency compared to bilateral aerosol delivery with 1 µm particles [0.48 versus 0.1% (Si et al., 2013)].

The main reason for this extremely large difference between left and right olfactory deposition can be seen in Figures 4, 5. For the non-pulsating flow, streamlines impinge on the top part of the cavity after the nasal valve, and the air passes with relatively high-speed through the right olfactory region. In contrast, due to the oscillating nature of the pulsatile flow, streamlines after the nasal valve are moderated and mixed resulting in lower air speed than that of the non-pulsatile (constant) flow when passing through the right olfactory region. Only for the right nasal cavity, one can observe some streamlines going into the olfactory region, but none are going into the left olfactory region. Since reduced air flow towards the epithelium results in reduced impaction of particles, pulsation reduces particle deposition for all regions where particle deposition is mostly affected by convective transport (impaction) of particles. This is the case for all regions, except the left olfactory epithelium region where convection is virtually absent. Thus, we assume that in this region the oscillatory-flow-induced particle deposition, which has also been reported for non-ventilated paranasal sinuses (Möller et al., 2014), accounts for this finding. This almost complete lack of convective transport into the olfactory region (at least for the left side) can be interpreted as a “protective” mechanism of the olfactory nerve. This nerve is supposed to be sensitive to adsorbed smell-related molecules. Since these molecules reach the olfactory region quite effectively via diffusion (as indicated by large C values in this region, which is indicative of the diffusion of air molecules), the lack of convective transport protects the olfactory region from dust particles, which might deposit onto the olfactory tissue and impair the sensitivity of the olfactory nerve. This hypothesis would have to be substantiated by investigating nasal cavities from other individuals.

Hence, targeting of 1 µm particles to the olfactory tissue with bi-directional flow appears to be difficult. While pulsatile aerosol flow leads to a 3-fold enhanced deposition (DE) in the left olfactory epithelium region, there is a 3.9-fold DE reduction in the right one. The net effect is a more uniform aerosol deposition in the left and right olfactory region, but there is still a ca. 50-fold higher DE in the right olfactory region (inlet side). This suggests that the olfactory region is either not sufficiently conducive to acoustic resonance (which is possibly, since it is completely open on one side of the cavity, rather than fully enclosed with a small opening as required for a Helmholtz resonator) or its resonance frequency is too far away from 45 Hz. We consider the former more likely than the latter.

The effects of different parameters such as morphological differences in various individuals, different inlet flow rates, and particle sizes could play significant role in optimization of particle deposition in olfactory epithelium which should be optimized in future studies.