Development of a realistic human respiratory tract cast representing physiological thermal conditions

Abstract

Toxicological assessment of inhaled aerosols and their effects on respiratory tissue cells requires accurate measurements of particle delivery, properties, evolution, and deposition in the human respiratory tract. These measurements are affected by both anatomical features and physiological factors, such as thermal conditions in the respiratory system. We constructed a segmented cast model, based on a realistic geometry of the human respiratory tract, equipped with features to control the temperature of the air flow. We then evaluated the thermal equilibrium of the air flow through the cast using both experimental measurements and computational fluid dynamics modeling. Uniform temperature distribution in the cast shell was achieved using parallel connection of the cast segments to a thermal bath by flow splitter. Air temperature inside the cast was shown to require <1 min to equilibrate with the temperature of the cast during internal circulation of hot water. Air flow temperature was shown to equilibrate with the cast temperature in the mouth–throat region, and a uniform temperature distribution was achieved in the other segments, resembling the thermal conditions of respiratory flow in the human lung. We showed that the cast successfully represents the physiological thermal conditions of the human respiratory system.

1. Introduction

Toxicological assessment of aerosols is often performed in vivo to assess the biological effects of compounds on living organisms (Höfer et al. 2004). Legal and ethical mandates to reduce animal involvement in toxicological assessments following the 3Rs principle (replacement, reduction, and refinement) (Baumans 2004; Festing and Wilkinson 2007) have promoted the development of in vitro and in silico models to complement in vivo exposure systems (Holmes, Creton, and Chapman 2010). The applicability of such models depends on the extent to which they simulate the living system and represent physiologically relevant scenarios (Davila et al. 1998).

Inhalation exposes lung tissue cells to aerosols consisting of vapor in gas phase and dispersed particles in solid or liquid phases. Particle transport through the respiratory tract and deposition on airway surfaces depends on the particles physical characteristics (Hinds 2012; Srirama et al. 2012). Smaller particles are subject to interactions with the carrier gas molecules through Brownian motion (Gupta and Peters 1985); larger particles are not affected by such interactions, but their inertia and weight can drive them out of the flow stream and trigger deposition on respiratory tract walls (Crane and Evans 1977). The biological effects of particle deposition depend on the delivered dose of the substance, which is typically defined as the received particle mass per unit of time and surface area (Donaldson et al. 2002; Pöschl 2005).

The refinement of in vitro exposure systems offers increased capacity and reliability in studying direct aerosol deposition in cell cultures (Paur et al. 2011; Thorne and Adamson 2013; Majeed et al. 2014). To represent realistic conditions, these systems require adjustment of the aerosol concentration to achieve a delivered aerosol mass that is comparable with in vivo exposure conditions (Mülhopt et al. 2009). The estimated delivered aerosol mass is obtained using analytical and computational models (Comouth et al. 2013; Secondo, Liu, and Lewinski 2017). Recently Lucci et al. (2018) developed a general model, not system specific, which has been thoroughly compared and validated with the in vitro exposure system measurements (Comouth et al. 2013; Fujitani et al. 2015).

Advances in computational fluid dynamics (CFD) have made it feasible to realistically estimate the retention of aerosol in the human respiratory tract under different inhalation regimes (Rostami 2009; Longest and Holbrook 2012; Rygg and Longest 2016). These computational models consider distinct deposition mechanisms such as impaction, diffusion, and gravitational settling, and can be applied for a wide range of particle sizes (Frederix et al. 2017). CFD models have been used with complex morphometric data from the human upper respiratory tract (Rahimi-Gorji et al. 2015; Frederix et al. 2018). In addition, a wide variety of aerosol properties have been modeled for various applications such as cigarette smoke particles (Zhang, Kleinstreuer, and Hyun 2012) and drug-delivery spray particles (Kimbell et al. 2007). Extended CFD models have also been developed to simulate aerosol evolution in the respiratory tract, namely the condensation/evaporation processes of liquid particles (Longest and Hindle 2010; Asgari, Lucci, and Kuczaj 2018). Despite extensive progress in the development of computational models, there are still challenges in model verification and validation under complex experimental conditions (Oldham 2006).

In parallel with CFD modeling, development of in vitro cast models enables researchers to measure the delivered aerosol dose in different regions of the respiratory tract such as trachea (Schlesinger and Lippmann 1976) and under various inhalation flow and aerosol exposure conditions (Chan and Lippmann 1980; Cheng, Zhou, and Chen 1999). The cast models are also used to verify the computational models and assess their accuracy (Longest and Vinchurkar 2007). Some cast models are based on idealized geometries of the tracheobronchial tree (Grgic et al. 2004b), such as the Weibel model (Yeh and Schum 1980). Idealization of the geometry, however, does not address the influence of the anatomical complexity on air flow and aerosol dynamics. Thus, cast models developed for realistic geometries use imaging tools such as computed tomography and magnetic resonance imaging (Grgic et al. 2004b; Lizal et al. 2015). Furthermore, anatomical variability of respiratory tract is taken into account using statistical analysis and geometry optimization (Borojeni et al. 2015). Most investigations to date have concentrated on the anatomical features of the human respiratory tract; less attention has been paid to reproducing physiologically relevant conditions such as temperature and humidity in an in vitro cast.

Studies by Ferron, Haider, and Kreyling (1985), Eisner and Martonen (1989), Daviskas, Gonda, and Anderson (1990) used 1D mathematical models and experiments to evaluate the temperature distribution in the human tracheobronchial system. Ferron, Haider, and Kreyling (1985) approximated temperature and water vapor concentration in the human upper airways solving the 1D transport equations. Daviskas, Gonda, and Anderson (1990) obtained the analytical solution for air flow temperature assuming a cylindrical geometry for different segments of the respiratory tract. Eisner and Martonen (1989) developed a surrogate tracheobronchial system with a multicomponent physical model, in which the average Reynolds number value for the airflow within each respective compartment equals the value within the corresponding bronchial airway generation of Weibel’s model. They used the surrogate model to estimate the temperature in various segments of the surrogate bronchial system. The evaluation of temperature in the respiratory tract is particularly important for estimating liquid particles retention in the respiratory tract, as the particles are subject to influential interactions with the carrier vapor phase and condensation/evaporation processes.

We addressed some of the above-mentioned challenges by developing an in vitro respiratory tract cast model that is based on realistic anatomical features. The cast is designed to control the air flow temperature and maintain it at a constant and realistic level. This thermal-control feature is useful for conducting aerosol evolution and deposition measurements for evolving liquid particles. In this article, we present the cast development procedure, technical features, and validation by a combination of CFD modeling and experimental analysis of internal thermal equilibrium with and without air flow.

2. Materials and methods

The cast development procedure and our measurement methods are outlined in this section comprising the following steps: preparation of respiratory tract geometry, computer-aided design (CAD) of the cast, material selection for 3D printing, rapid prototyping technique, and experimental setup used for the cast assessment.

2.1. Respiratory tract geometry and cast development

Three-dimensional (3D) cast geometry was based on the structure of the human respiratory tract (Zhang, Kleinstreuer, and Hyun 2012). Figure 1a shows the 3D reconstructed model of the human airways up to the ninth generation of the tracheobronchial tree. This geometry was used as a reference model for the CAD of the internal surface. The airway pathway was processed in MeshLab open source 3D processing software (ISTI-CNR, Pisa, Italy) and imported as a stereolithography file into SolidWorks 3D design software (Dassault Systèmes, Vélizy-Villacoublay, France). The external shell surface was constructed by applying an oversizing parameter on the initial geometry, and a Boolean operation yielded the internal shell surface. Channels for water circulation were designed inside the shell with additional Boolean operations. The structure obtained was split into five segments. Water connection interfaces were also designed in SolidWorks and added to each segment of the cast. The modular segmented format was designed to represent the mouth cavity, throat, trachea, first and second branch generations, and third to fifth branch generations of the tracheobronchial tree as shown in Figure 1b. The mouth inlet was designed as a cylindrical extrusion from the mouthpiece segment. Each cast segment has a cylindrical flange at the junction with its neighboring segment for assembly of the cast segments. Placeholders for four pins were designed at each junction to orient the segments correctly (Figure 1b, inset). The seals were designed to guarantee air-tightness at the junctions.

Figure 1. (a) The geometry of the human respiratory tract, obtained by computed tomography (Zhang, Kleinstreuer, and Hyun 2012). The red line indicates the branches selected for the in vitro cast up to the fifth generation. (b) CAD model of the in vitro cast. Tract segments are numbered in black; cast outlet connections are numbered in red. The figure inset magnifies a flange between two cast segments with pin placeholders to ensure correct orientation at assembly.

2.2. Water capillaries for temperature control

Water capillaries were embedded in the cast shell for thermal control. These capillaries were designed to circulate water from a thermal water bath at a specified temperature. The water capillaries were built into all cast segments covering a large surface area to facilitate heat exchange with the air flow inside the cast. Water capillaries were connected with inlet/outlet connections at the junctions between neighboring cast segments (see Figure 2). Water circulation between cast segments can be serial, with pipes connecting the inlets/outlets of neighboring segments, or a parallel connection of all segments directly to the thermal water bath using a flow splitter. Figure 2 shows the water capillaries and their connections in a cross-section of the geometry.

Figure 2. CAD model of the cast showing the inlet/outlet connections for water circulation and water capillaries inside the cast shell.

2.3. Material selection

The material selected for 3D prototyping should satisfy several property requirements, including high thermal conductivity and compatibility with various chemical compounds. The thermal conductivity of the material is vital to controlling the temperature of the air flow inside the cast using water circulation. High thermal conductivity of the resin material reduces the time required to reach thermal equilibrium between the circulating water in the cast shell and air flow through the cast. We therefore selected high-temperature resin from Formlabs (Somerville, MA, USA), as it satisfies the required high thermal conductivity and is compatible with 3D printing (Formlabs 2018). The mechanical tensile strength of the resin, 51.1 MPa, makes the printed cast robust to withstand the forces implemented during experiments. The thermal capacity of the material was assessed using differential scanning calorimetry. Thermal assessment was performed for both raw (uncured) and ultraviolet light (UV)-cured samples (Figure 3). The raw sample was polymerized by sharp changes in heat flow during the first measurement cycle (Figure 3a). The UV light-cured sample did not polymerize with increased temperature indicating that it was fully polymerized in UV-curing process. To obtain the thermal capacity values, the measurements had to be calibrated using a well-documented reference material. We used sapphire as the reference. Table 1 lists the calculated values. These values are in the range of heat capacity values for polymers, which is 1.0–2.5 J/g/°C. Finally, the solvent compatibility of the resin with different chemical compounds was tested by extended exposure of the resin to the target compounds reported in Table 2. The resin did not appear to degrade following contact with the chemical compounds, and its weight gain over 24 h of exposure was below 1%.

Figure 3. Differential scanning calorimetry measurements of heat flux against temperature for (a) raw material and (b) ultraviolet light-cured material.

Table 1. Differential scanning calorimetry measurements of the resin heat capacity compared with the reference material (Sapphire) heat capacity.

Temperature (°C)	Heat capacity (J/g/°C)
	Sapphire (reference)	High-temperature resin
37.5	0.80	1.54
47	0.82	1.55
57	0.84	1.56
67	0.85	1.57

Table 2. The solvent compatibility test of the resin exposed to the target chemical compounds.

Solvent	Percentage of mass gain over 24 h of exposure
Acetic acid	<1
Glycerol	<1
Propylene glycol	<1
Acetone	<1

The data in this table reports the mass absorption measurement over 24 h of exposure to each chemical compound.

2.4. 3D printing

3D printing of the cast was conducted using the stereolithography apparatus (SLA) technique with a Miicraft 100X 3D printer (Jena, Germany) for light-sensitive photopolymer resins. The technique is used to produce the cast in a layer-by-layer approach with photopolymerization, in which a 2D light pattern is used to form the molecules chain into polymers, building the 3D object (Hull 1986). The source of emitted light in this technique is a digital light projector. SLA printing provides a minimum print layer thickness of 10 µm indicating the smallest features of the geometry, which can be resolved in 3D printed object. This thickness is finer than the one provided by other technologies, such as fused filament fabrication. The high-temperature resin we selected is compatible with the SLA printing technique, and was used to print the cast (Figure 4).

Figure 4. Three-dimensional cast printed using the SLA technique.

The cast is designed to be used for air flow measurements and aerosol exposure experiments, so the internal surface roughness of the printed cast is an important factor. We measured surface roughness with a Bruker optical profilometer (Billerica, MA, USA); Figure 5 shows the acquired data. The layer-by-layer fabrication of the cast is visible in the profilometry figure. The average surface roughness measured was approximately 0.55 µm and the maximum measured surface roughness was 3.54 µm indicating the final print resolution. The maximum roughness occurs in Y direction, where the profile covers the height of the deposited resin layer. The data indicated that the SLA printing method yielded a smoother cast surface than the fused filament fabrication technique, which offers an approximate average surface roughness of2 µm.

Figure 5. Surface roughness image (left) and measurements by optical profilometer (right). The surface pattern indicates the layer-by-layer fabrication by SLA.

2.5. Temperature and flow measurement setup

Figure 6 shows the experimental setup used to control air flow through the cast and to measure the surface and internal temperatures. Air flow rates at the cast outlets were fixed using red-y gas flow controllers (Vgtlin Instruments GmbH, Aesch, Switzerland). The flow controllers were connected to a KNF N 920 G vacuum pump (KNF Neuberger GmbH, Freiburg, Germany). The air flow rate at the inlet connection was measured using a TSI 4100 flowmeter (TSI Inc., Shoreview, MN, USA), to ensure that the inlet flow equaled the total flow pulled from the outlets; the equivalence of the inlet and outlet flow rates indicates the air-tightness of the cast shell. The water capillaries in the cast shell were connected to an Isotemp water bath (Thermo Fisher Scientific, Waltham, MA, USA) controlling the water temperature. Water capillaries in the cast segments were connected in parallel using a flow splitter. Temperature was measured by two techniques. A Fluke Ti300 infrared thermal camera (Fluke Corporation, Everett, WA, USA) was used to measure the external surface temperature of the cast with and without internal air flow. To measure air flow temperature, we inserted a thermocouple probe through the inlet connection inside the mouth cavity and along the central axis. Temperature was measured inside the cast at different locations to map the temperature gradient along the path of the air flow.

Figure 6. Flow and temperature measurement setup.

2.6. Computational model

We used a CFD model to simulate the air flow distribution in the cast and obtain the air flow rates through the outlet branches. The calculated flow rates from simulations are used in the experimental setup to adjust the flow controllers. The CFD model is then used to simulate the air flow heat transfer in the cast. After validation by experimental temperature measurements at various points inside the cast, the model can give a full picture of thermal equilibrium in the cast. The model includes the continuity equation for mass, Navier–Stokes equations for momentum, and energy equation for temperature. These equations are: (1) $${\partial }_t\rho +\nabla ·\left(\rho \mathbf{u}\right)=0,$$(1) (2) $${\partial }_t\left(\rho \mathbf{u}\right)+\nabla ·\left(\rho \mathbf{uu}\right)=-\nabla p-\left(\nabla ·\mathit{\tau }\right),$$(2) (3) $$c_p\left[{\partial }_t\left(\rho T\right)+\nabla ·\left(\rho \mathbf{u}T\right)\right]=\nabla ·\left(\kappa \nabla T\right)-\left(\mathit{\tau }:\nabla \mathbf{u}\right)+{\mathrm{D}}_{t}p.$$(3)

In this set of equations, ρ is the gas density, $\mathbf{u}$ represents the gas velocity, p is the pressure, t is the time, and T represents the temperature. The material properties of the air including μ as the dynamic viscosity, c_p as the heat capacity at constant pressure, and κ as the heat conductivity are temperature dependent. $\mathit{\tau }$ denotes the rate of the strain tensor defined as $\mathit{\tau }=-\mu \left[\nabla \mathbf{u}+{\left(\nabla \mathbf{u}\right)}^{\mathrm{T}}\right]+\frac{2}{3}\mu \left(\nabla ·\mathbf{u}\right)\mathbf{I}$ with $\mathbf{I}$ as the identity tensor.

The mathematical equations were solved with the AeroSolved (AeroSolved 2018) open source software developed in the OpenFOAM framework. Time integration was performed using a first-order Euler implicit scheme. The discretization of flow velocity and temperature was performed with a second-order upwind scheme. All gradient and Laplacian terms were discretized with a second-order accuracy central scheme.

The computational domain for solving the equations was the internal space of the cast geometry (Figure 1b). The computational mesh generated for the domain consisted of 1,264,893 cells. We tested the independence of the numerical results against discretization by refining the mesh to 2,376,540 cells. The flow rates and temperature values obtained using the refined mesh were altered from the coarse mesh calculations by <10%. We used a constant flow boundary condition at the inlet, 1.5 L/min, with fixed room temperature value, 22 °C. Considering a fully developed flow at the cast outlets in steady-state condition, a fixed pressure boundary was applied. A no-slip condition for the flow velocity was applied at the walls. The temperature of the wall was kept constant at 37 °C.

3. Results and discussion

3.1. CFD simulations and flow rate measurements

We performed CFD calculations to measure the air flow velocities and streamlines for a given flow rate at the inlet connection. Figure 7 shows the air flow streamlines color-coded by velocity through the respiratory tract geometry. The inlet air flow rate was 1.5 L/min. Maximum velocity occurred in the jet flow through the inlet connection, decreasing as the air moved through the cast. The density of the streamline in each branch is an indication of the flow rate; Table 3 lists the flow rate in each outlet branch. The outlet flow rate was highest at the two largest outlet connections, numbered 1 and 3. The calculated flow rates were used to adjust the flow controllers in the experimental setup (Figure 6).

Figure 7. Flow streamlines inside the cast color-coded by air flow velocity.

Table 3 Flow rates calculated in each outlet branch of the cast for an inlet flow rate of 1.5 L/min.

Outlet branch number	1	2	3	4	5	6
Flow rate (L/min)	0.77	0.12	0.51	0.06	0.02	0.02

The outlet branches are indicated according to their numbering in Figure 1b.

3.1.1. Temperature measurements

The temperature of the cast was controlled by circulation of hot water in the capillaries inside the cast shell. These water capillaries were connected to a thermal bath, providing stable temperature. Figure 8 shows the surface temperature of the cast, measured using an infrared thermal camera. The temperature of the thermal bath was kept at 40 °C, but the cast remained at an average temperature of 37 °C with ∼1 °C of variation between the cast segments. The temperature change from the thermal bath to the cast (40 to 37 °C) was due to heat loss as the water flowed through the connecting pipes. The temperature variation along the cast was minimized by parallel-connecting the cast segments to the thermal bath using a flow splitter.

Figure 8. External surface temperature of the cast, measured using an infrared thermal camera. Hot water is circulated in the cast shell capillaries using the connected tubes.

The air flow temperature is influenced by the heat exchange with the cast shell. We have measured the air temperature in the cast lumen using a thermocouple inserted through the inlet connection. Figure 9 shows the air temperature measured versus time. The measurement was started without the circulating hot water in the cast and, therefore, the initial temperature was in equilibrium with the ambient conditions. During the measurement, the water circulation was started at t = 120 s and, consequently, air temperature started to increase. The time constant of the heating process can be estimated referring to the lumped heat capacity analysis, in which the convective thermal process is modeled as a first-order system (Lewis, Nithiarasu, and Seetharamu 2004). With this analysis, the time constant of the process equals to the time that the air temperature reaches 63% of the initial temperature difference between air and the cast. This time equals to t = 45 s. Thus, the dynamic response of the air temperature to the changes implemented in the cast temperature can be expected in the order of a minute.

Figure 9. Air temperature measurements inside the cast. Hot water circulation commenced at t = 120 s. The time constant of the heating process, τ, equals the time required to reach 63% of the initial temperature difference between the air and the cast.

The temperature measurement depicted in Figure 9 was performed for stationary conditions without air flow inside the cast, which would change the thermal equilibrium in the cast. Figure 10 shows the temperature distribution inside the cast for an air flow rate of 1.5 L/min, calculated using CFD simulations. The air flow temperature at the inlet was 22 °C (laboratory temperature), and it equilibrated with the cast temperature of 37 °C in the throat region. The thermal equilibration depends on the temperature differences between the air flow and the cast, but also on the flow residence time in the cast, which is indicated by the inlet air flow rate. The lowest temperature was observed for the jet flow entering the cast through the inlet connection. In the mouth cavity, the air velocity was reduced, increasing residence time for conductive heat transfer from the cast to the air flow.

Figure 10. Air temperature distribution inside the cast with a flow rate of 1.5 L/min, calculated by CFD simulation. The air flow temperature at the inlet was T = 22 °C, and the cast wall temperature was maintained at T = 37 °C.

We measured the temperature of the air flow inside the cast experimentally using an inserted thermocouple. The temperature was measured inside the mouth cavity by inserting the thermocouple probe through the inlet connection and along the central axis (Figure 11). Ambient air entering the cast through the inlet connection increased in temperature, reaching a temperature just below that of the cast wall (T = 37 °C). This measurement is in good agreement with the CFD prediction (Figure 11), validating the CFD model of thermal equilibrium inside the cast (Figure 10). We compared our results with the analytical solution obtained by Daviskas, Gonda, and Anderson (1990) for the temperature distribution in the tracheobronchial system and we have a variation of <1 °C. The variation in the calculated temperatures can be explained by the fact that Daviskas, Gonda, and Anderson (1990) obtained their empirical coefficients by fitting their equation to experimental data for inspiration flow rates of 15 and 100 L/min. Therefore, using their obtained equation for other flow rates would result in some error.

Figure 11. Air flow temperature measured using a thermocouple inserted through the inlet connection. The air flow rate is 1.5 L/min. The temperature was measured along the line (x) and compared with CFD predictions.

4. Conclusion

We developed a human respiratory tract cast model with realistic geometry that included the first five generations of the tracheobronchial tree. The cast is equipped with water capillaries controlling the temperature of the air flow to represent the physiological conditions in the human lung. The cast design, comprising five separable segments, is suitable for aerosol exposure studies as it facilitates the measurement of particle deposition within the respiratory tract. The resin material selected can withstand high temperatures without polymerizing, and the average surface roughness of the printed cast’s internal surface did not exceed 0.55 µm. Surface characteristics have been shown to be important for the local deposition of particles in airway models (Holbrook and Longest 2013).

The temperature control feature of the cast was tested for uniform temperature distribution in the cast shell and maximum variation of 1 °C was obtained using parallel connection of the cast segments to the thermal bath by flow splitter. We are also able to simulate the physiological temperature gradient with the temperature increasing from throat to deeper segments of the tracheobronchial system (Ferron, Haider, and Kreyling 1985; Eisner and Martonen 1989), since the water capillaries in different segments of the cast can be connected to heat sources with various temperatures.

The air temperature inside the cast was measured without air flow, and required <1 min to equilibrate with the temperature of the cast during circulation of hot water. This relatively short response time justifies the selection of the resin material and will be useful for inhalation studies with different thermal conditions. This time scale is still much longer than one inhalation cycle time (4 s), and hence does not allow the control of the thermal conditions within one breathing cycle. This can be a place for future improvement to reduce the response time to order of a second for more detailed investigations. Air flow rates at the casts outlet branches were calculated by CFD simulations and used to adjust the settings for the experimental measurement. Air temperature inside the cast was also measured with air flow, and was in good agreement with CFD predictions. We then used the validated computational model to assess the thermal equilibrium in the full cast. We showed that the flow temperature equilibrates with the cast temperature in the mouth–throat region, and a uniform temperature distribution is achieved in the other segments. This resembles the thermal conditions of respiratory flow in the human lung.

The cast was developed to obtain an anatomically realistic respiratory tract model that imitates physiological lung thermal conditions and is suitable for aerosol exposure studies. This cast will facilitate evaluation of size evolution and deposition of liquid particles following condensation/evaporation events under varying temperature conditions. Further studies and cast development toward the integration of features controlling air flow humidity to represent physiological conditions in the respiratory tract are warranted to promote in vitro research in inhalation toxicology and aerosol physics.

Disclosure statement

The authors declare no competing financial interests.

Additional information

Funding

The research described here was funded by Philip Morris Products S.A., Switzerland (part of Philip Morris International group of companies).