Comparative effects of a candidate modified-risk tobacco product Aerosol and cigarette smoke on human organotypic small airway cultures: a systems toxicology approach

Abstract

Using an in vitro human small airway epithelium model, we assessed the biological impact of an aerosol from a candidate modified-risk tobacco product, the tobacco heating system (THS) 2.2, to investigate the potential reduced risk of THS2.2 aerosol exposure compared with cigarette smoke. Following the recommendations of the Institute of Medicine and the Tobacco Product Assessment Consortium, in which modified-risk tobacco products assessment should be performed in comparison with standard conventional products, the effects of the THS2.2 aerosol exposure on the small airway cultures were compared with those of 3R4F cigarette smoke. We used a systems toxicology approach whereby elucidation of toxic effects is derived not only from functional assay readouts but also from omics technologies. Cytotoxicity, ciliary beating function, secretion of pro-inflammatory mediators and histological assessment represented functional assays. The omics data included transcriptomic and miRNA profiles. Exposure-induced perturbations of causal biological networks were computed from the transcriptomic data. The results showed that THS2.2 aerosol exposure at the tested doses elicited lower cytotoxicity levels and lower changes in the secreted pro-inflammatory mediators than 3R4F smoke. Although THS2.2 exposure elicited alterations in the gene expression, a higher transcriptome-induced biological impact was observed following 3R4F smoke: The effects of THS2.2 aerosol exposure, if observed, were mostly transient and diminished more rapidly after exposure than those of 3R4F smoke. The study demonstrated that the systems toxicology approach can reveal changes at the cellular level that would be otherwise not detected from functional assays, thus increasing the sensitivity to detect potential toxicity of a treatment/exposure.

---

Introduction

Studies of the human small airway region have relied on autopsy specimens or tissue samples removed by surgery.¹ The small airway is a region of the distal lung with bronchioles smaller than 2 mm in diameter, without cartilage.^2,3 Modern technologies, such as high-resolution computed tomography and helium magnetic resonance imaging, are available for studying the distal lung compartment. However, such techniques are limited to assessing the functional and morphological changes, and not suitable for a mechanistic investigation (e.g., cellular and molecular assessment).¹ Mechanistic evaluations are often conducted using animal models to assess the overall effects on the whole organism. Efforts to combine bioinformatics methods to identify the best animal models that represent specific populations of patients with respiratory diseases have been proposed.⁴ However, because of species differences, results obtained from animal studies cannot fully reflect the pathophysiology of human diseases.^5,6

Alternatively, in vitro cellular models can overcome the inter-species differences because such models can be derived from human cells. Despite the availability and convenience of traditional two-dimensional (2D) lung culture systems, they cannot mimic the spatial organization (cell–cell interaction and cell polarity) of the human airway epithelia.⁷ Advances in tissue engineering have enabled the development of more sophisticated cellular systems, e.g., three dimensional (3D) organotypic culture systems and organ-on-a-chip models, that better mimic the organization and functionality of the human tissue counterparts, than conventional 2D cellular systems.^7,8 One of the benefits of the 3D organotypic airway culture systems is that they are cultured at an air–liquid interface; thus, they can be directly exposed to inhalable gases and aerosol (e.g., smoke, aerosols, airborne particles, or nanoparticles). This setup eliminates the need to generate liquid smoke- or aerosol-fractions that are required for an exposure study using 2D submerged airway cells. Usually, a solution containing cigarette smoke (CS) fraction or CS extract, which does not represent whole CS, is prepared and then diluted in the culture media as the exposure. Three dimensional organotypic airway cultures can be exposed to whole CS at the apical side. This setup closely resembles the actual situation in the respiratory epithelium during smoking: CS passes across the smoker's respiratory epithelium where it comes in contact with the apical side of the cultures.

Recently, an air–liquid interface 3D organotypic small airway model became available (SmallAir™ ⁹). The model has a structure of a pseudostratified “tissue” layer composed of basal, ciliated, goblet, and club cells. This structure is similar to that of the human small airway epithelia. The culture is grown on a transparent membrane and is supplied with culture medium underneath.⁹ Using this newly developed organotypic model of the human small airway epithelium, we conducted the present study to compare the biological impact of an aerosol generated from a candidate modified-risk tobacco product (MRTP)¹⁰—the tobacco heating system (THS) 2.2¹¹—with that of a smoke generated from the 3R4F reference cigarette. Readers are referred to a separate publication for more detailed information on the THS2.2 product.¹¹

To assess the biological impact of the exposures, we followed a systems toxicology approach that involves integration of large biological data sets (e.g., omics data) with conventional toxicology readouts.¹² Different from traditional toxicology where visible adverse effects of a compound are often sought (e.g., cellular injury or cell death following high-dose exposure), systems toxicology aims to delineate cellular molecular changes in a system under stress or perturbation. The system's responses following a subtoxic level of exposure that is more relevant to the real-life situation is investigated.¹³ This concept is also highlighted under the toxicity testing in the 21^st Century vision and strategy¹⁴ and offers a “phenotypic anchoring” at the lower doses of toxicants, which “can help to explain a toxicant's mechanism of action […] before histopathological changes were seen […]”.¹⁵ Under such conditions, cellular and molecular responses following exposures will bring a mechanistic understanding of the exposure-induced impacts rather than findings that merely reflect the already-visible adverse effects.^16,17 Accordingly, the molecular mechanistic assessment to deduce the impact of THS2.2 aerosol exposure in the present study was performed in comparison with the benchmark 3R4F smoke-induced impact at a subtoxic concentration (approximately 0.15 mg nicotine per L 3R4F smoke). Cultures were also exposed to a higher dose of 3R4F smoke to demonstrate the progression of toxicity in the human organotypic small airway cultures following smoke exposure. We hypothesized that the cellular and molecular changes ensued following THS2.2 aerosol exposure compared with those following 3R4F smoke at subtoxic doses could provide some indication of adverse outcomes.

Materials and methods

Organotypic human small airway cultures

The organotypic human small airway culture model SmallAir™ was purchased from Epithelix (Plan-Les-Ouates, Geneva, Switzerland). All SmallAir™ cultures used in this study were reconstituted from primary human small airway epithelial cells of bronchiolar origin from the same donor, a healthy non-smoking 55 years-old female. Even though the use of cells from a single donor would only capture a donor-specific response, the decision was made to reduce the impact of donor-to-donor variability. The cells were grown in Transwell® inserts (with a diameter of 6.5 mm) and maintained in 12-well culture plates. The SmallAir™ cultures—fully differentiated upon arrival—were cultured at the air–liquid interface at 37 °C (5% CO₂, 90% humidity) for 12 days in a 12-well culture plate before the experiments. The cultures were maintained in SmallAir™ culture medium (0.7 mL per well) provided by the supplier, with a medium change every 2–3 days. After exposure, the medium was not changed until it was collected for various endpoint measurements (up to 72 h post-exposure).

Reference cigarette smoke and tobacco-heated aerosol

Mainstream CS was generated from 3R4F reference cigarettes, purchased from the University of Kentucky, Lexington, KY, USA.²¹ Mainstream test aerosol was generated from a candidate MRTP, the heat-not-burn-based technology THS2.2 (Philip Morris international R&D, Neuchâtel, Switzerland). The characteristics of THS2.2, including the specification of the product components and how the product operates, were reported in a separate publication.¹¹ The yields of harmful and potentially harmful constituents (HPHCs) in the THS2.2 mainstream aerosol were described in a previous publication;²² the influence of tobacco blends in the concentrations of HPHCs in the THS2.2 mainstream aerosol was reported in a separate publication.²³ In the present study, 3R4F cigarettes and THS2.2 sticks were conditioned according to ISO standard 3402.²⁴

Exposure setup

The smoke and aerosol were generated according to the Health Canada Intense smoking protocol (55 mL puff over 2 s, twice per min²⁵) with an 8 s pump exhaust time. Each 3R4F cigarette was smoked to a standard butt length (approximately 35 mm), and each THS2.2 stick was aerosolized for a total of 12 puffs per stick. Two independent 30-port carousel smoking machines (SM2000; Philip Morris, International), each used for 3R4F or THS2.2, were connected to a dedicated Vitrocell® 24/48 exposure system (Vitrocell Systems GmbH, Waldkirch, Germany) (Fig. 1, upper panel). The Vitrocell® 24/48 exposure system is equipped with a Dilution/Distribution module, in which fresh air can be added into each row to dilute the smoke or aerosol. During an exposure, smoke or aerosol is passed through the Dilution/Distribution module and distributed via trumpets (by negative pressure) into the Cultivation Base Module,²⁶ where biological culture systems were placed and exposed to smoke or aerosol at their apical sides.

Fig. 1 Experimental procedure and biological endpoints tested. Upper panel illustrates the smoking machines and Vitrocell® 24/48 exposure system used for exposing cultures to 3R4F smoke or THS2.2 aerosol. Lower left panel illustrates 3D organotypic cultures and the various biological endpoints assessed in the study. Lower right panel shows the number of experimental repetitions conducted for 3R4F smoke and THS2.2 aerosol exposure.

This study comprised three experimental phases/repetitions (Fig. 1, lower right panel). Each experimental phase was conducted within a week. For each experimental phase, a new batch of the human small airway cultures was obtained from the supplier, and three independent exposure runs were performed (three runs for the 3R4F smoke exposure and three runs for the THS2.2 aerosol exposure, paired with their corresponding air-exposed controls) (Fig. 1, upper panel). The apical sides of the cultures were washed on the Friday preceding to the week of exposure to rinse and standardize the mucus quantity across samples. Because the exposed samples were always paired with the air-exposed controls, a plausible day-to-day variability could be minimized. This design resulted in a total of nine replicate samples for each of the endpoints and exposure conditions. The organotypic small airway cultures were exposed to smoke or aerosol for 28 min. This duration was selected according to a previous finding regarding the sensitivity of organotypic bronchial cultures following CS exposure: a 28 min exposure to 3R4F smoke induced the highest concentration of secreted matrix metalloproteinase (MMP)-1 in bronchial epithelial organotypic cultures.²⁷ Various endpoints were assessed following the exposures (Fig. 1, lower left panel). Endpoints were not measured between the 4 h and 24 h post-exposure time points because sample collections and processing was practically not feasible. For gene expression (mRNA and miRNA profiling), samples were not collected immediately after exposure (0 h post-exposure) because the gene expression alterations were expected to be less robust at this post-exposure time point: a previous study reported that a low enrichment score was found in 3R4F CS-exposed samples 0.5 h post-exposure.²⁷ Furthermore, a previous study showed that the CS-induced nuclear factor erythroid 2-related factor 2 (Nrf2) promoter activation in bronchial epithelial organotypic cultures at 18 h post-exposure and 24 h post-exposure were not markedly different.²⁸ Nrf2 orchestrates cellular defense against cellular-damaging compounds, and is known to be upregulated upon CS exposure.²⁹

A set of small airway cultures was exposed to two dilutions of 3R4F smoke and to air, simultaneously, in one exposure plate (Fig. 1, upper panel):

• A 3R4F smoke dilution (7%) corresponding to a nicotine concentration of around 0.14 mg nicotine per L, termed the “3R4F (0.14)” group,

• A 3R4F smoke dilution (13%) corresponding to a nicotine concentration of around 0.26 mg nicotine per L, termed the “3R4F (0.26)” group,

• A 0% 3R4F smoke exposure, representing the air-exposed controls for the 3R4F exposure, term the “3R4F (Air)” group.

The concentrations of 3R4F smoke were selected based on previous studies using nasal and bronchial organotypic cultures^18–20 and used as the benchmark exposures to investigate the impact of THS2.2 aerosol exposure on the small airway cultures. The 7% 3R4F smoke dilution, a subtoxic level of 3R4F smoke exposure, was shown previously to elicit a no-observed-adverse-effect-level where sufficient cellular and molecular alterations were detected without culture damage.^18–20 The 13% 3R4F smoke dilution, a toxic level of 3R4F smoke exposure, was previously shown to induce overt effects, i.e., culture damage.^18–20 A pilot concentration range study was also performed to confirm that similar findings would be observed in small airway cultures (ESI Fig. 1†).

The nicotine yield from one THS2.2 stick is approximately 30% lower than the yield from one 3R4F cigarette.^11,22 Because THS2.2 is “designed to significantly reduce or eliminate the formation of HPHCs in the inhaled aerosol while preserving as much as possible the taste, sensory experience, nicotine delivery profile and ritual characteristics of cigarette”,¹¹ nicotine was used as the standard compound to compare the biological impact of THS2.2 aerosol and 3R4F smoke in this study. A set of cultures was exposed to THS2.2 aerosol dilutions at nicotine concentrations matched to those of the diluted 3R4F smoke (at least for the two concentrations) and to air, simultaneously, in one exposure plate (Fig. 1, upper panel):

• A THS2.2 aerosol dilution (14%) corresponding to a nicotine concentration of around 0.14 mg nicotine per L, termed the “THS2.2 (0.14)” group,

• A THS2.2 aerosol dilution (24%) corresponding to a nicotine concentration of around 0.30 mg nicotine per L, term the “THS2.2 (0.30)” group,

• A THS2.2 aerosol dilution (31%) corresponding to a nicotine concentration of around 0.45 mg nicotine per L, term the “THS2.2 (0.45)” group,

• A 0% THS2.2 aerosol exposure, representing the air-exposed controls for the THS2.2 exposure, term the “THS2.2 (Air)” group.

Nicotine measurement in the trapped smoke/aerosol

Previous correlation data (not shown) were used to estimate the corresponding nicotine concentrations in a given diluted 3R4F smoke or THS2.2 aerosol; the correlation between nicotine concentrations in the trapped smoke and dilutions of 3R4F smoke had been established before.²⁶ For the present study, the actual concentrations of nicotine in each of the 3R4F smoke or THS2.2 aerosol dilutions were also determined within each week of the experimental phase, to ensure that the target concentrations were appropriately achieved. For this, a trapping experiment was conducted for each of the dilution conditions. Each diluted 3R4F smoke or THS2.2 aerosol was trapped in EXtrelut® 3NT columns (Merck Millipore, Billerica, MA, USA) as described previously.²⁶ The EXtrelut® 3NT column was placed at the end (i.e., the exhaust) of the first row of the Dilution/Distribution module of the Vitrocell® 24/48 exposure system. For this reason, the smoke or aerosol was not distributed to the trumpets—under which culture models would be located—but trapped in the EXtrelut® columns in the first row. Therefore, these trapping experiments are not feasible when the small airway cultures are exposed. Nicotine concentrations were measured from the eluted samples using gas chromatography-flame ionization detection as previously described.²⁶

Histology processing

The histological samples were obtained only from cultures harvested at the 48 h and 72 h post-exposure time points. We hypothesized that morphological alteration would occur at a later time point after exposure and after molecular changes took place, as reported in another study.³⁰ After three rinses with phosphate-buffered saline (PBS), the culture was fixed for 2 h in freshly prepared 4% paraformaldehyde, and then collected from the insert for paraffin embedding using a Leica ASP300S tissue processor (Leica Biosystems Nussloch GmbH, Nussloch, Germany). Sections of 5 μm thickness were obtained using a microtome and mounted on glass slides. The slides were subsequently transferred to an automated slide stainer (Leica ST5020) for staining with hematoxylin (Merck Millipore) and eosin (Sigma-Aldrich, St Louis, MO, USA) (H&E), and alcian blue (Sigma-Aldrich) (AB). The stained slides were then covered with glass coverslips using a Leica CV5030 fully automated coverslipper. Digital microscopic images were generated using a Hamamatsu NanoZoomer 2.0 slide scanner (Hamamatsu Photonics, K.K., Hamamatsu City, Japan).

Adenylate kinase (AK) release assay

AK activity in the basolateral medium of the cultures was measured from different cultures at various time points post-exposure using the ToxiLight™ bioassay kit (Lonza, Basel, Switzerland) according to the manufacturer's instructions. The luminescence signal was measured using a FluoStar Omega reader (BMG Labtech GmbH, Ortenberg, Germany). For each of the three experimental repetitions (i.e., each culture batch), the values of the luminescence signal were normalized to the mean value of the positive and negative control samples as previously described;^19,31 the formula is given in ESI Materials and methods 1.† For positive controls, triplicate samples within a batch were treated with Triton X-100 (at a 1% final concentration) and used to derive the value of 100% cytotoxicity. For negative controls, triplicate untreated samples were used. The mean values (from the total three experimental repetitions) of the normalized relative luminescence units were then reported in the figures as percentages.

Ciliary beating functionality

Video recordings of the cultures were taken before exposure, immediately after exposure (0 h), and 24 h, 48 h, and 72 h post-exposure using a digital high-speed video camera (Sony CXD V60, Sony, Tokyo, Japan). Ciliary beating videos were not generated 4 h post-exposure because first, we hypothesized that the readouts would be similar to that taken 0 h post-exposure (as we observed in another study, ciliary beating and mucocilliary transport at 2, 4, 6, and 24 h post-exposure to a 7 min smoke exposure were largely comparable, data not shown) and second, because measurements at 4 h post-exposure time point would overlap with the sample collections and be logistically challenging. For the video recordings, the camera was connected to an inverted microscope system (Leica DMi8). Images were taken at a rate of 90 frames per second. The ciliary beating functionality was evaluated by four measures: the weighted frequency, the uniformity of the detected frequency, the active area, and the power of the detected signal (fast Fourier transformation [FFT]). Analyses were conducted on a total of 512 video frames recorded from the center of the insert surface. For each pixel, the mean of the 512 frames was subtracted, and an FFT and an approximate Bartlett's Kolmogorov–Smirnov test were performed on the pixel intensity. The weighted frequency was calculated as follows: the mean of the dominant frequency detected was weighted by its FFT power magnitude for each video if the pixel was active (p ≤ 0.001) and its dominant frequency was in the range of 0–20 Hz. The uniformity of the detected frequency was calculated using the mean of the Kolmogorov–Smirnov statistic; i.e., the maximum difference of the normalized cumulative FFT spectrum and the uniform cumulative distribution function. The active area was defined as the proportion of pixels that showed an unadjusted Bartlett's Kolmogorov–Smirnov p-value ≤0.001. The strength of the ciliary beating signal was finally estimated as the sum of the FFT power spectrum in the range of 2.5–20 Hz.

Measurement of secreted pro-inflammatory mediators

For each of the medium samples collected at the given post-exposure time point (24 h, 48 h, and 72 h post-exposure), an independent culture was used (i.e., a cross-section sample, the medium was not changed for the entire duration of the post-exposure period). Based on previous observations (data not shown), the levels of secreted mediators 4 h post-exposure were low (i.e., comparable to the levels of the air control with high variability); therefore, their concentrations were not measured at this time point. Multi-analyte profiling (MAP) of pro-inflammatory mediators secreted was performed using commercially available Milliplex panels (Merck Millipore) with Luminex® xMAP® technology (Luminex, Austin, TX, USA)-based analysis according to the manufacturer's instructions. Samples were analyzed on a FlexMap3D® equipped with xPONENT software v.4.2 (Luminex). The following analytes (mediators) were measured: chemokine (C–C motif) ligand (CCL) 5 and CCL-20; colony-stimulating factor (CSF) 2 and CSF-3; chemokine (C–X–C motif) ligand (CXCL) 1 (also known as GROα) and CXCL-10; epidermal growth factor (EGF); interleukin(IL)-1α, IL-1β, IL-6, and IL-8; vascular endothelial growth factor (VEGF) α; TNFα; soluble intercellular adhesion molecule (sICAM) 1; MMP-1 and MMP-9; tissue inhibitor of metalloproteinase (TIMP) 1, and thymic stromal lymphopoietin (TSLP). As positive controls (to assess the capacity of the cultures to respond to stimuli) from each of the experimental repetitions, three cultures were treated for 24 h with a combination of 10 ng mL⁻¹ TNFα and 10 ng mL⁻¹ IL-1β in PBS added to the basolateral medium; results are presented in ESI Table 1.† As negative controls, from each experimental repetition, three cultures were treated for 24 h with PBS in the basolateral medium.

RNA/microRNA (miRNA) isolation and array analyses

Total RNA, including miRNA, was isolated from the small airway epithelium culture at 4, 24, 48, and 72 h post-exposure, using a previously published method.^19,20,31,32 For the mRNA array, 100 ng of total RNA were reverse-transcribed to cDNA using an Affymetrix® HT 3′-IVT PLUS kit (Affymetrix, Santa Clara, CA, USA). The cDNA was labeled and amplified to complementary RNA (cRNA). The fragmented and labeled cRNA was hybridized to a GeneChip® Human Genome U133 Plus 2.0 Array (Affymetrix) in a GeneChip® Hybridization Oven 645 (Affymetrix) according to the manufacturer's instructions. Arrays were rinsed and stained on a GeneChip® Fluidics Station FS450 DX (Affymetrix) using the Affymetrix® GeneChip® Command Console® Software (AGCC software v-3.2, protocol FS450_0001). The RNA integrity number (RIN) values of the 215 samples (three experimental repetitions) were distributed between 6.8 and 10 (mean: 9.27).

For the miRNA array, a FlashTag™ Biotin HSR kit (Affymetrix) was used to label the miRNA. Two hundred nanograms of total RNA containing low molecular-weight RNA were biotinylated and hybridized to miRNA arrays version 4.0 (Affymetrix) in a GeneChip® Hybridization Oven 645 (Affymetrix) according to the manufacturer's instructions. Arrays were rinsed and stained on a GeneChip® Fluidics Station FS450 DX (Affymetrix) using the Affymetrix® GeneChip® Command Console® Software (AGCC software v-3.2, protocol FS450_0002).

Finally, the arrays were scanned using a GeneChip® Scanner 3000 7G (Affymetrix). Raw images from the scanner were saved as DAT files. The AGCC software automatically gridded the DAT file image and extracted probe cell intensities into a CEL file.

Processing raw CEL files from the mRNA microarray

The raw CEL files were background-corrected, normalized, and summarized using frozen-robust multi-array analysis.³³ Background correction and quantile normalization were used to generate microarray expression values from all arrays passing quality controls (QC), and were performed using the custom CDF environment HGU133Plus2_Hs_ENTREZG v16.0,³⁴ as previously described in greater detail.¹⁹

Analysis of mRNA data

For each experimental factor combination item, concentration, and post-exposure duration, a model to estimate the treatment effect was fitted with limma,³⁵ by including the covariate exposure run as a blocking variable to account for the pairing during an exposure run (exposed vs. air control). The p-values for each computed effect were adjusted across genes using the Benjamini–Hochberg false discovery rate (FDR) method.³⁶ Differentially expressed genes (DEGs) were defined as a set of genes whose FDR was <0.05. The mRNA array dataset is available in the Arrays Express repository (ID: E-MTAB-6098).

Network perturbation amplitude (NPA) analysis of transcriptomic data

The NPA methodology was described in greater detail in a previous publication.³⁷ Briefly, the methodology aims to contextualize transcriptome profiles (treated vs. control, or exposed vs. air control) by combining the alterations in gene expression into differential node values, i.e., one value for each node of a causal network model.³⁸ Relevant network models used for the analysis in this study are listed in ESI Table 2.† The differential node values are determined by fitting procedures inferring the values that best satisfy the directionality of the causal relationships contained in the network model (e.g., positive or negative signs). NPA scores carry a confidence interval accounting for experimental variation, and the associated p-values are computed. In addition, companion statistics, derived to inform the specificity of the NPA score to the biology described in the network models, are reported as *O and K* if their p-values fall below the threshold of significance (0.05). A network is considered significantly affected by exposure if the three values (the p-value for experimental variation, *O, and K*) are below 0.05.³⁷

A systems-wide metric for biological impact, the biological impact factor (BIF),^39,40 summarizes the impacts of the exposure on the cellular system into a single (absolute) number, thus enabling a simple and high-level evaluation of the treatment effects across multiple time points. Calculating the BIF requires the collection of all applicable hierarchically structured network models (ESI Table 2†), and involves aggregating the NPA values of the individual networks.

Processing raw CEL files from the miRNA microarray

The 215 delivered CEL files were read using the oligo package in the Bioconductor suite of microarray analysis tools for the R statistical software environment.^41–43 QC of the miRNA raw data was performed as previously described,¹⁹ using the arrayQualityMetrics package,⁴⁴ and resulted in the exclusion of nine CEL files. Normalized probe-level data were obtained by applying robust multi-array normalization and summarized at the probe-set level with the median polish method.⁴⁵ Using the annotation provided by Affymetrix and the latest miRNA nomenclature according to miRBase v21,⁴⁶ only the probe sets pertaining to human miRNA were kept in the expression matrix. Additionally, the probe sets that were not available on Affymetrix GeneChip® miRNA 3.0 arrays were not considered, to maintain compatibility with other published studies using organotypic cultures.^19,20,31,32 Only the miRNA probe sets with significantly higher intensity values than their matched background probes must be considered as “detected”.⁴⁷ A p-value threshold of 0.01 was selected to determine the detection calls based on Wilcoxon tests. If a miRNA probe set was detected in more than 50% of the samples in at least one sample group, it was kept for further analysis; otherwise, it was discarded. This process leads to a final expression matrix containing 206 columns and 594 rows, corresponding to the accepted samples and filtered miRNA probe sets, respectively. The miRNA array dataset is available in the Arrays Express repository (ID: E-MTAB-6004).

Analysis of miRNA data

For each comparison (exposed vs. air control, at a given stimulus, concentration, and post-exposure time point), a submatrix was extracted from the global expression matrix by keeping only those samples belonging to the corresponding treatment or control groups, as well as the miRNA probe sets that were detected in more than 50% of the samples in at least one of the two sample groups. A linear model for differential expression was applied to the resulting submatrices using the moderated t statistic implemented in the limma package.⁴⁸ The models included an additional variable to take into account the exposure runs. Subsequently, adjusted p-values were obtained following multiple testing corrections with the Benjamini–Hochberg FDR.³⁶ miRNAs below the FDR threshold of 0.05 were considered differentially expressed. The miRNA expression changes are displayed as a heatmap figure, in which they are sorted according to their pattern; the clustering was done based on the “affinity propagation” algorithm implemented in the R “APCluster” package.⁴⁹

Statistical analysis

Statistical analysis was performed with SAS software version 9.2 (SAS Institute, Cary, NC, USA) on the following data: nicotine measurements, luminescence signals of the AK assay, ciliary beating function, and the Luminex-based measurement of secreted mediators. Mean values and standard error of the means are reported (unless otherwise specified). Comparisons of an exposed sample and its air control (i.e., the paired-sample from the same exposure run) were conducted using a paired t-test. When applicable (i.e., for comparing the relative cytotoxicity in the THS2.2-exposed samples with relative toxicity in the 3R4F-exposed samples), the comparison was conducted after subtracting the values of the corresponding air controls (i.e., the paired samples). Next, the comparison was performed using a t-test corrected for non-equal variance (Satterthwaite correction). Numerical values from secreted mediator analysis (Luminex assay) were transformed using the natural log transformation.

Results

Nicotine concentrations in the diluted smoke/aerosol

The dilutions of smoke/aerosol that would match the target nicotine concentrations were estimated from previous experiments as reported before.²⁶ In the present study, to ensure that the target nicotine concentrations were met according to the study design (see the Materials and methods section and Fig. 1), samples of the diluted 3R4F smoke and THS2.2 aerosol were taken throughout the study period, i.e., within the three experimental repetition weeks.

Table 1 shows that the target nicotine concentrations were properly achieved from the specified applied dilutions of 3R4F smoke or THS2.2 aerosol. For example, the actual nicotine concentrations in the 7% dilution of 3R4F smoke were 0.14 mg L⁻¹, 0.16 mg L⁻¹, and 0.16 mg L⁻¹ in the first, second, and third experimental repetitions, respectively. These values were aligned with the target nicotine concentration (i.e., 0.15 mg nicotine per L). Similarly, the actual nicotine concentrations in the 13% dilution of THS2.2 aerosol were 0.11 mg L⁻¹, 0.17 mg L⁻¹, and 0.15 mg L⁻¹ in the first, second, and third experimental repetitions, respectively. These concentrations corresponded to the target nicotine concentration of 0.14 mg L⁻¹. The results confirmed that this approach allowed for a comparison of THS2.2 aerosol and 3R4F smoke at similar nicotine concentrations.

Table 1 Target and actual nicotine concentrations in 3R4F smoke or THS2.2 aerosol dilutions

	Experimental repetitions	Experiment 1	Experiment 2	Experiment 3
3R4F (0.15)	Target nicotine concentration (mg L^−1)	0.15	0.15	0.15
	Applied dilution in Vitrocell® (vol/vol)	7%	7%	7%
	Actual nicotine measured in EXtrelut® (mg L^−1)	0.14 ± 0.00 (N = 3)	0.16 ± 0.01 (N = 3)	0.16 ± 0.01 (N = 3)
3R4F (0.26)	Target nicotine concentration (mg L^−1)	0.26	0.26	0.26
	Applied dilution in Vitrocell® (vol/vol)	13%	13%	13%
	Actual nicotine measured in EXtrelut® (mg L^−1)	0.22 ± 0.01 (N = 3)	0.23 ± 0.01 (N = 3)	0.32 ± 0.02 (N = 3)
THS2.2 (0.14)	Target nicotine concentration (mg L^−1)	0.14	0.14	0.14
	Applied dilution in Vitrocell® (vol/vol)	13%	13%	13%
	Actual nicotine measured in EXtrelut® (mg L^−1)	0.11 ± 0.01 (N = 3)	0.17 ± 0.02 (N = 3)	0.15 ± 0.00 (N = 3)
THS2.2 (0.30)	Target nicotine concentration (mg L^−1)	0.30	0.30	0.30
	Applied dilution in Vitrocell® (vol/vol)	24%	24%	24%
	Actual nicotine measured in EXtrelut® (mg L^−1)	0.29 ± 0.01 (N = 3)	0.34 ± 0.02 (N = 3)	0.29 ± 0.02 (N = 3)
THS2.2 (0.45)	Target nicotine concentration (mg L^−1)	0.45	0.45	0.45
	Applied dilution in Vitrocell® (vol/vol)	31%	31%	31%
	Actual nicotine measured in EXtrelut® (mg L^−1)	0.45 ± 0.01 (N = 3)	0.48 ± 0.01 (N = 3)	0.41 ± 0.02 (N = 3)

---

Effects on cytotoxicity and culture morphology

The exposure effect was first assessed by measuring AK release as a marker of cytotoxicity and then by evaluating the morphology of the exposed small airway tissue cultures. At the earliest time point after exposure (4 h, Fig. 2A), AK release into the culture medium did not differ among the air, smoke, and aerosol exposure groups. The relative cytotoxicity levels were mainly observed in samples exposed to 3R4F (0.26) smoke; cytotoxicity increased with post-exposure duration (up until approximately 30% relative toxicity compared with the 100% cytotoxicity level in the Triton X-treated samples). Following 3R4F (0.15) smoke exposure, a noticeable increase in relative cytotoxicity was observed only 72 h post-exposure, although the difference with levels in the air controls did not reach statistically significance. THS2.2 aerosol exposure at any concentration was not associated with alterations in relative cytotoxicity at all post-exposure time points tested. One statistically significant difference was observed between the samples exposed to THS2.2 (0.30) aerosol and air (at 72 h post-exposure); however, this change was considered not biologically relevant because the relative toxicity level was only 0.34%.

Fig. 2 Cytotoxicity following exposure. Mean cytotoxicity levels evaluated by an adenylate kinase (AK) release assay at various time points post-exposure. AK levels were normalized relative to the positive and negative controls (see Materials and methods section). Nicotine concentrations in 3R4F smoke or THS2.2 aerosol are indicated for each group (mg L⁻¹, x-axis). ↔ indicates p ≤ 0.05 compared with the corresponding air controls. # indicates a p ≤ 0.05 difference from 3R4F smoke exposure at a similar nicotine concentration. In the case of THS2.2 (0.45), 3R4F (0.26) was used as the comparison group.

Histological sections from the exposed cultures were only collected at the 48 h and 72 h post-exposure time points for morphology evaluation after H&E/AB staining (Fig. 3). We assumed that alterations in culture morphology did not occur immediately after exposure; therefore, histological samples were not obtained at the earlier post-exposure time points (4 h and 24 h). The air-exposed and unexposed samples (i.e., incubator controls, data not shown) exhibited similar morphology regarding the thickness of the pseudostratified epithelium, fraction of AB-positive cells, and presence of cilia. Compared with the samples exposed to air, the samples exposed to 3R4F (0.15) smoke had exhibited reduced cilia numbers and increased frequency of empty spaces between cells (see arrowheads in Fig. 3), which could be attributed to a lower cell–cell adherence, and induced detachment at the suprabasal layer. Squamous cells and apoptotic cells were also detected in the 3R4F (0.15) smoke-exposed samples. When the small airway tissue cultures were exposed to the highest concentration of 3R4F smoke (0.26 mg nicotine per L), more pronounced damage was observed 48 h (Fig. 3A) and 72 h (Fig. 3B) after exposure. In contrast, THS2.2 aerosol-exposed samples (at all concentrations tested) did not exhibit any apparent morphological alterations compared with the air-exposed controls at either post-exposure time point (Fig. 3).

Fig. 3 Culture morphology 48 h and 72 h post-exposure. Representative images of hematoxylin and eosin (H&E) and alcian blue (AB)-stained small airway culture sections observed 48 h (A) and 72 h (B) after exposure. Arrowheads indicate empty spaces between cells.

Effects on ciliary beating function

The coordinated movement (beating) of cilia serves to evacuate mucus, which traps toxicants, odorants, and other particulates from the respiratory tract.⁵⁰ In asthmatic patients⁵¹ and in smokers,⁵² this defense mechanism—known as mucociliary clearance—is perturbed. Similar observations have been recorded in in vitro airway models exposed to whole CS or CS condensate.^20,53,54 The small airway tissue cultures have functional cilia at their apical side, as reported by the supplier,⁹ which we also observed (Fig. 3).

Fig. 4 shows that compared with air exposure, the 3R4F smoke exposure, immediately after exposure (0 h), was linked to a reduction in the weighted frequency of ciliary beating, uniformity of the beating frequency, active area where beating was detected, and power of the beating signal in a concentration-dependent manner. These reductions (of the four different ciliary function readouts) were also observed at the 24, 48, and 72 h post-exposure time points, with the exception of a slight recovery of the power of the beating signal 72 h post-exposure (although the signal remained significantly reduced than that of the air-exposed controls).

Fig. 4 Impact of exposure on ciliary beating functionality. Ciliary beating functionality in the small airway cultures was assessed longitudinally before, immediately after (0 h), and 24 h, 48 h, and 72 h after exposure. (A) Weighted frequency (Hz) is the mean frequency over the pixel, weighted by the fast Fourier transformation (FFT) power at the pixel-dominant frequency. (B) Frequency uniformity (arbitrary unit, AU) is an index expressing the distribution of the detected frequency over the FFT spectrum (0 = blank noise; 1 = unique frequency). (C) The active area (%) is the percentage of detected pixel that differs from the blank noise. (D) The log₁₀(FFT power) is an estimate of the detected power based on the beating signal of the ciliary movement. Nicotine concentrations (mg L⁻¹) in 3R4F smoke or THS2.2 aerosol are indicated for each group.

THS2.2 aerosol exposure at all concentrations tested did not significantly affect the weighted frequency of ciliary beating, uniformity of the beating frequency, active area where beating was detected, or power of the beating signal, almost at all post-exposure time points tested (Fig. 4). An exception was seen at the 24 h post-exposure time point, where the ciliary function readouts were slightly reduced in a concentration-dependent manner compared with the readouts in the air-exposed controls.

Inflammatory responses following exposure

Various publications have reported that organotypic airway cultures can secrete multiple cytokines, chemokines, and growth factors.^18,19,55,56 In this study, we measured the concentrations of secreted pro-inflammatory mediators from the basolateral media of the cultures collected at various post-exposure time points. In general, the concentrations of the secreted mediators increased with post-exposure duration. Fig. 5 shows, the mediators measured at the 72 h post-exposure time point. The concentrations of mediators collected at the 24 h and 48 h post-exposure time points are given in ESI Table 1.† Independent of the post-exposure time points, the effects of exposure followed a similar trend: generally, we typically observed increased mediator levels following 3R4F smoke exposure, and smaller changes following THS2.2 aerosol exposure, relative to the levels seen in the air-exposed samples.

Fig. 5 Profiles of secreted pro-inflammatory mediators following exposure. Mean concentrations of pro-inflammatory mediators measured in the basolateral media of the cultures 72 h after exposure. ↔ indicates p ≤ 0.05 compared with the corresponding air controls. # indicates a p ≤ 0.05 difference from 3R4F smoke exposure at a similar nicotine concentration. In the case of THS2.2 (0.45), 3R4F (0.26) was used as the comparison group.

Fig. 5 shows that when compared with the mediator levels in the air-exposed controls, 3R4F smoke exposure was linked to a concentration-dependent increase in the levels of the majority of mediators measured 72 h post exposure. The concentrations of TIMP-1 and VEGFA were greater in the cultures exposed to 3R4F (0.15) smoke than in the air-exposed controls; however, their levels following 3R4F (0.26) smoke exposure did not differ from the levels following air exposure. In contrast, CSF-2, IL-6, TNFα, and CSF-3 protein levels in the media were higher in cultures exposed to 3R4F (0.26) smoke than in cultures exposed to air. THS2.2 aerosol exposure (at all concentrations tested) was linked to fewer changes in mediator levels (relative to the air exposure) for the majority of mediators, compared with the changes associated with 3R4F smoke exposure. In general, we did not observe a concentration-dependent alteration in mediator concentrations following THS2.2 aerosol exposure—except for VEGFA, for which a significant difference between THS2.2 (0.45) aerosol exposure- and the air-exposure was observed.

Global mRNA and miRNA alterations following exposure

We aimed to deduce possible mechanisms associated with the exposure impact using a systems toxicology approach. We leveraged omics technologies (mRNA and miRNA microarrays) to complement the functional readouts reported in the previous sections (i.e., cytotoxicity, histological assessment, ciliary beating function, and secreted pro-inflammatory mediators). As mentioned in the introduction, the mechanistic understanding of the toxicological response was evaluated following exposure at subtoxic concentrations at a no-observed-adverse-effect level, thus avoiding supraphysiological concentrations of exposures.¹³ Further, cellular and molecular alterations in severely damaged samples would depict merely the cell death response (adverse effects),^16,17 which we demonstrated before.¹⁸ For this reason, the samples exposed to the high concentration of 3R4F smoke, i.e., the 3R4F (0.26) group, were not subjected to the mRNA and miRNA microarrays because overt visible morphological changes had occurred (i.e., the tissue damage shown in Fig. 2).

The global mRNA expression from samples exposed to 3R4F smoke or THS2.2 aerosol was compared with that of their respective air-exposed samples (i.e., systems response profiles⁵⁷). Fig. 6A shows that the greatest number of DEGs was found 24 h post-exposure to 3R4F (0.15) smoke (4343 + 3581 = 7924 genes). The number of DEGs following THS2.2 aerosol exposure at a given concentration peaked at 4 h post-exposure; these numbers also increased in a concentration-dependent manner.

Fig. 6 Alterations in mRNA and miRNAs expression in small airway epithelial cultures following exposure. (A) Systems response profile of the small airway cultures following exposure, referring to the number of differentially expressed genes that were up-regulated or down-regulated following exposure. (B) Significantly altered miRNAs following exposure (listed on the left side of the heatmap). The color gradient reflects the alterations in miRNA expression compared with the respective air controls. The numbering on the right side of the heatmap indicates the cluster of miRNAs in which similar patterns were observed. Nicotine concentrations in the smoke or aerosol are indicated for each group at the bottom of the heatmap (mg L⁻¹). FDR, false discovery rate.

In addition, miRNA profiles were generated from the same culture samples used for the mRNA analysis. Fig. 6B shows that more miRNAs were significantly altered following 3R4F (0.15) smoke exposure than following THS2.2 aerosol exposure at all concentrations tested. The largest number of miRNAs altered following 3R4F (0.15) smoke exposure was found at the 72 h post-exposure time point. Following THS2.2 aerosol exposure at any tested concentration, only a handful of miRNAs were significantly altered: miR-4521, miR-92a, miR-27a, and miR-29b (4 h post-exposure) and miR-138 (48 h post-exposure).

Our approach to quantifying systems perturbation leverages a computational methodology³⁷ and a collection of causal networks depicting biological pathways/processes,³⁸ described in the Materials and methods section. Using the systems response profiles (i.e., the profile of gene alterations following exposure), we derived and quantified the impact of the different exposure conditions on the perturbation of causal biological networks.³⁷ The degree of perturbation is termed the NPA score. Accordingly, the NPA score for a given network reflects a quantitative measure of the exposure-induced impact on the biological processes/pathways modeled in the network.

A heatmap representing the NPA scores for each exposure condition and at each post-exposure time point is shown in Fig. 7A. The network models can be grouped into four major categories: Cell Fate, Cell Proliferation, Cell Stress, and Inflammatory Process Network. Exposure to 3R4F (0.15) smoke was linked to greater NPA scores 4 h post-exposure, among the post-exposure time points tested, for the majority of networks. Exceptions were observed for mTOR, Hox, Hedgehog, Cell Interaction, and Endoplasmic Reticulum Stress networks, for which significant NPA scores were not observed following 3R4F (0.15) smoke exposure at any of the post-exposure time points tested.

Fig. 7 Biological Impact Factor (BIF) derived from cumulated network perturbations following 3R4F and THS2.2 exposure compared with the air controls. (A) Heatmap of network perturbation amplitude (NPA) scores of biological networks impacted by 3R4F smoke and THS2.2 aerosol exposure. The network names are listed on the left side of the heatmap, with the corresponding network family on the right side of the heatmap. The color gradient represents the NPA scores, which were normalized to the maximum NPA score per network. The star symbols (*) in the heatmap indicate that the network is significantly impacted by exposure (i.e., the three values—the confidence interval, *O, and K* statistics—are below 0.05, as described in the Materials and methods section). (B) Percentage relative BIF is plotted on the y-axis. The highest BIF value is taken as 100% as the reference (REF). Nicotine concentrations in the smoke or aerosol are indicated for each group (mg L⁻¹).

At similar nicotine concentrations, THS2.2 (0.14) aerosol exposure was linked to lower NPA scores for most networks, compared with the scores observed for 3R4F (0.15) smoke exposure. For all three THS2.2 aerosol concentrations tested—THS2.2 (0.14), THS2.2 (0.30), and THS2.2 (0.45)—similar patterns of altered NPA scores were observed: almost all networks were impacted most at the 4 h post-exposure time point than at the later time points. For the Hox, Hedgehog, and Cell Cycle networks, higher perturbations were observed at 48 h post-exposure than at the other time points following THS2.2 aerosol exposure at any tested concentration. The perturbation scores for the highest concentration of THS2.2 aerosol exposure, THS2.2 (0.45), were still lower than those observed for 3R4F (0.15) smoke exposure.

A BIF, which is an agglomeration of the NPA scores of all networks, can be computed to deduce the overall impact of exposure on the systems, as described previously.^39,40,57Fig. 7B shows the overall BIF scores for each of the exposure conditions in the context of the biological processes/pathways represented in the networks. The highest BIF score was observed for 3R4F (0.15) smoke exposure 4 h post-exposure (thus, the BIF is considered to be 100% assigned as the reference group marked with “REF” in Fig. 7B), and declined with post-exposure duration. A similar observation was recorded for THS2.2 aerosol exposure at all concentrations tested. The highest BIF scores for THS2.2 (0.14), THS2.2 (0.30), and THS2.2 (0.45) groups were observed consistently at the 4 h post-exposure time point: 14.60%, 40.53%, and 65.19%, respectively.

Discussion

Increasing concerns about the toxicity of airborne contaminants, such as nanoparticles, aerosols, and pollutants, have underscored the need for relevant testing approaches to assess the toxicity risk in humans. Similar importance has been placed on developing new tobacco products, including electronic cigarettes and candidate MRTPs.¹⁰ This study used a recently developed small airway culture model grown at the air–liquid interface (SmallAir™ ⁹), whereby a systems toxicology approach¹² was applied. The study aimed to assess the biological impact of an aerosol from THS2.2, a candidate MRTP, compared with that of smoke from 3R4F cigarettes. The human small airway cultures were reconstituted from primary small airway epithelial cells of a single donor: a 55 years-old female non-smoker, apparently healthy. We acknowledge that the use of a single donor captures only the donor-specific response; however, data from a single-source culture would reduce the influence of a donor-to-donor variability, thus increasing the odds of differentiating the exposure effects. Predicting reliably human population responses to toxic compounds from in vitro data is challenging. Inherent human variations attributed to genetic and other host factors cannot be easily modeled in vitro, and should be addressed before a successful quantitative in vitro-to-in vivo extrapolation.⁵⁸ Nonetheless, in our previous in vitro work,²⁰ we determined that the exposure-related transcriptome changes were largely similar between two donors provided by two different providers, despite minor differences.

In the present study, we used a high concentration 3R4F smoke exposure to show a concentration–response relationship (of progressive toxicity until adverse effects are evident), as observed in other studies.^18–20 Here we showed that marked overt effects in the cultures could be seen following exposure to 3R4F (0.26) smoke—corresponding to a 13% 3R4F smoke dilution. One indication of an adverse effect is a disruption of the system's function.⁵⁹ This was particularly evident from the impaired ciliary beating function observed in this exposure. The levels of the beating frequency, uniformity of the frequency, active area, and power of the beating signals were lower in samples exposed to 3R4F (0.26) smoke than in those exposed to air, immediately after exposure (0 h) and 24, 48, and 72 h post-exposure. This finding suggested that mucociliary clearance following exposure to 3R4F (0.26) smoke was compromised. Mucociliary clearance is a mechanism controlled by a coordinated movement of cilia (the ciliary beating) to transport mucus along the respiratory tract.⁵¹ The reduced uniformity of the beating frequency likely illustrated a disturbed propulsion of the mucus layer. We further postulated that the reduced area where active beating was detected could be caused by loss of cilia; we observed pronounced damaged to the epithelial layer in this group (both 48 and 72 h post-exposure). The adverse effects following 3R4F smoke exposure at this concentration were similarly observed in our previous work.¹⁸ In that study, non-transient and overall perturbations of the network models were detected, confirming that molecular changes following a toxic exposure concentration (one that induces adverse effects) simply reflect the already-visible cellular injury.^16,17 The non-transient perturbations further indicate that the effect is likely to be adverse.⁵⁹

In the context of systems toxicology assessment, the mechanisms of action by which the exposure elicited biological impact were investigated at subtoxic concentrations (at the no observed-adverse-effect level). Such concentrations are likely to be more relevant to revealing the mechanism of the exposure-induced effects: they are sufficient to induce cellular and molecular alterations without resulting in adverse effects.^13,15 Thus, in the present study, we assessed the mechanisms of the exposure-induced effects (evaluated from the global mRNA and miRNA changes) in comparison with those observed in the benchmark 3R4F (0.15) smoke-exposed samples. The network-based analysis of transcriptomes showed that 3R4F smoke exposure at 0.15 mg nicotine per L (corresponding to approximately 13% smoke dilution) had exerted the highest impact on the majority of networks (except for MTor, Hox, and Cell Interaction networks) under the conditions tested. The severity of the exposure impact (reflected from the NPA scores) decreased with post-exposure duration. This pattern of a decreasing impact was observed similarly following THS2.2 aerosol exposure (at all concentrations tested), although the NPA scores were lower and decrease more sharply than NPA scores following 3R4F (0.15) smoke exposure. Notable concentration-dependent perturbations of the Hedgehog, Hox, and Cell Cycle networks 48 h post-exposure to THS2.2 aerosol exposure were observed. The nodes in the Hedgehog network that were markedly impacted included those representing decreased expression of Gli proteins and inhibited catalytic activity of Smoothened (Smo) proteins (ESI Fig. 2†). The Hedgehog signaling pathway has been linked to repair functions in airway epithelia.^60–62 We postulated that the Hedgehog network impacted following THS2.2 aerosol exposure, in conjunction with the perturbed Hox and Cell Cycle networks, simply indicates an adaptive homeostasis response.⁶³ Cells can cope following exposure to multiple forms of stressors; tiny variations in oxygen, oxidants, pH, and even exercise can trigger biochemical, post-translational and gene expression changes. If the toxicity or damage is severe, irreparable injury would occur; such circumstances could trigger genetic and metabolic reactions, including the induction of DNA repair capacities.⁶³ We observed that the perturbations of the response to DNA damage, necroptosis, Wnt, and Mapk networks were significant at all post-exposure time points following 3R4F (0.15) smoke exposure, but only significant at the early time points after exposure to THS2.2 aerosol. It could be further hypothesized that THS2.2 aerosol and 3R4F (0.15) smoke exposures elicited different cellular responses.

The study allowed us to reasonably explore the regulatory relationship between mRNA and miRNA changes following exposure, because the mRNA and miRNA were isolated from the same culture. The number of mRNAs altered following exposures decreased with post-exposure duration; in contrast, the number of altered miRNAs increased with post-exposure duration. This inverse relationship was more perceptible in the samples exposed to 3R4F (0.15) smoke than those exposed to THS2.2 aerosol—mainly because fewer miRNA alterations were detected in the samples exposed to THS2.2 aerosol. We hypothesized that the later changes in miRNAs were a consequence of the perturbation of the mRNA profiles observed at the earlier post-exposure time points. This delay in miRNA alteration may indicate a cellular response to resolve an earlier perturbation in the mRNA profiles.⁶⁴ A similar observation in bronchial cultures following CS exposure was previously reported.²⁰

Concentrations of the mediators secreted from the cells into the basolateral media increased with post-exposure duration. Nonetheless, mediator levels following THS2.2 aerosol exposure were only slightly changed (Fig. 5 and ESI Table 1†), with the exception of VEGFA. A concentration-dependent increase in secreted VEGFA proteins was observed in the culture media following THS2.2 aerosol exposure. VEGFA has been shown to promote proliferation of lung parenchymal cells⁶⁵ and reduce alveolar epithelial apoptosis in the context of wound repair.⁶⁶ These changes may be associated with the recovery phase following exposure (adaptive homeostasis response) postulated before. We also detected greater variability in the concentrations of CCL2, IL-1A, and IL-1B that hinders a meaningful biological interpretation. Nonetheless, the data suggested that the concentrations of secreted CCL2, IL-1A, and IL-1B were not influenced by the exposures.

The concentrations of nicotine in the smoke or aerosol cannot directly infer the actual nicotine dose delivered to the cell cultures because smoke/aerosol is evolving (partitioning of compounds is highly variable throughout the in vitro piping system and, similarly, throughout the human respiratory tract). Nonetheless, the subtoxic and toxic effects observed here were unlikely to be attributed to the nicotine by itself, but the entire gas phase and particulate phase delivered to the cultures. Moreover, particle size largely influences the deposition of aerosols.^67,68 Although not reported here, separate work has been done to evaluate the particle size in 3R4F smoke and THS2.2 aerosol: Schaller and colleagues²² reported that the average mass median aerodynamic diameter of THS2.2 aerosol (0.7 μm, with a mean geometric standard deviation of 1.5) is similar to that of 3R4F smoke (0.8 μm, with a mean geometric standard deviation of 1.3). Furthermore, the efficiency of aerosol deposition in the Cultivation Base Module of the Vitrocell® 24/48 exposure system was investigated in a separate study.⁶⁹ Steiner and colleagues reported a linear correlation between the expected deposited concentrations and the detected values of an aerosol with a mean aerodynamic particle sizes of 0.8 μm. In that study,⁶⁹ the test aerosol was prepared from a mixture of glycerin and disodium fluorescein. Disodium fluorescein was added as a tracer; thus, measurements can be relatively easy, precise, and low-cost. Evaluating particle deposition in the exposure system remains an ongoing activity in our group. Such evaluation can be done using biological test systems (e.g., cell cultures) or a surrogate matrix (e.g., PBS, culture medium). Although cell cultures would be preferable because they reflect the characteristics of the biological test systems used in the respective toxicity studies, 3D organotypic cell culture models are expensive for large-scale production. We further recognize that the dose translatability from an in vitro exposure setup (e.g., in this study) to in vivo conditions (e.g., in smoker) cannot be extrapolated directly. First, the structure of exposure chambers is different from that of the human respiratory tract. Therefore, dynamic changes of aerosol (often referred to as aerosol evolution) along the flow path of an exposure chamber are different from those along the human's respiratory tract. Second, puffing topography differs among smokers, increasing the difficulty in simulating the smoke deposition along the respiratory tract. Puffing parameters used in experimental studies follow a smoking regimen defined by regulatory bodies that allows for standardization across laboratories and studies, but does not represent all human smoking behaviors.⁷⁰ Such data from machine measurements are limited for product hazard assessment but “not intended to be nor valid as a measure of human exposure risk”.⁷⁰ Our group has been working on computational modeling approaches to better characterize aerosol deposition and to estimate the total and regional aerosol deposition in human lungs, similar to the work of Pichelstorfer and colleagues.⁷¹ Nonetheless, despite these limitations, within the experimental setting used in the present study, the biological impact of THS2.2 aerosol exposure on small airway cultures could be compared with that of 3R4F smoke exposure.

Conclusion

With new technologies, investigations at the systems level have become possible. We perceived biological processes/pathways involved in the exposures from the computational network analysis (the NPA methodology), although they were not evident from the more traditional toxicity readouts, i.e., culture morphology, cytotoxicity, ciliary beating function, and secretion of pro-inflammatory mediators. The NPA methodology revealed not only a qualitative measure, but also a quantitative measure of the exposure impact in the context of the biological networks assessed here: the highest biological impact was observed 4 h post-exposure to 3R4F smoke at 0.15 mg nicotine per L (100% impact). In contrast, at the same post-exposure time point (4 h), THS2.2 aerosol at a comparable nicotine concentration (0.14 mg nicotine per L), elicited only a 15% relative biological impact. Overall, compared with the biological impact of 3R4F smoke, the results showed that the aerosol from the candidate MRTP THS2.2 elicited lower impact in all measured endpoints in the human small airway cultures. The effects of THS2.2, if observed, were mostly transient and diminished more rapidly after exposure. Despite the challenges and complexity of utilizing large data sets, this study demonstrated that the systems toxicology approach could provide an additional layer of toxicity-testing data at the pre-clinical stage.

Author contributions

A.R. Iskandar, Y. Martinez, F. Martin, P. Leroy, M.C. Peitsch, and J. Hoeng contributed to the conception and/or design of the work. A.R. Iskandar, Y. Martinez, L.O. Torres, S. Majeed, C. Merg, K. Trivedi, E. Guedj, S. Frentzel, and N.V. Ivanov contributed to the data acquisition. A.R. Iskandar, Y. Martinez, F. Martin, W.K. Schlage, P. Leroy, A. Sewer, C. Mathis, and J. Hoeng contributed to the data analysis and interpretation. A.R. Iskandar, Y. Martinez, F. Martin, W.K. Schlage, C. Mathis, and J. Hoeng drafted the article. All authors critically reviewed and approved the article.

Conflicts of interest

All authors are employees of Philip Morris International R&D, Philip Morris Products S.A. (part of Philip Morris International group of companies). W.K. Schlage is contracted and paid by Philip Morris International R&D. Philip Morris International is the sole source of funding and sponsor of this project.

Acknowledgements

The authors would like to acknowledge Dr Marja Talikka for the scientific review of the manuscript. The authors would like thank the following individuals for contributing their technical expertise: Abdelkader Benyagoub and Camille Schilt for histology; Maica Corciulo for Luminex-based measurements; Karine Baumer and Dariusz Peric for RNA purification; Remi Dulize for microarray experiments.

References

1. P. Burgel, The role of small airways in obstructive airway diseases, Eur. Respir. Rev., 2011, 20, 023–033.
2. P. J. Barnes, Small Airways in COPD, N. Engl. J. Med., 2004, 350, 2635–2637.
3. M. Contoli, J. Bousquet, L. Fabbri, H. Magnussen, K. Rabe, N. Siafakas, Q. Hamid and M. Kraft, The small airways and distal lung compartment in asthma and COPD: a time for reappraisal, Allergy, 2010, 65, 141–151.
4. C. S. Stevenson, S. Sridhar and J. E. Phillips, Predicting drug efficacy using integrative models for chronic respiratory diseases, Inflamm. Allergy Drug Targets, 2013, 12, 124–131.
5. B. J. Canning, Modeling asthma and COPD in animals: a pointless exercise?, Curr. Opin. Pharmacol., 2003, 3, 244–250.
6. P. A. Martorana, E. Cavarra, M. Lucattelli and G. Lungarella, Models for COPD involving cigarette smoke, Drug Discovery Today: Dis. Models, 2006, 3, 225–230.
7. S. Breslin and L. O'Driscoll, Three-dimensional cell culture: the missing link in drug discovery, Drug Discovery Today, 2013, 18, 240–249.
8. K. BeruBe, M. Aufderheide, D. Breheny, R. Clothier, R. Combes, R. Duffin, B. Forbes, M. Gaça, A. Gray and I. Hall, In vitro models of inhalation toxicity and disease, Altern. Lab. Anim., 2009, 37, 89–141.
9. S. Huang, B. Boda, J. Vernaz, E. Ferreira, L. Wiszniewski and S. Constant, Establishment and characterization of an in vitro human small airway model (SmallAir™), Eur. J. Pharm. Biopharm., 2017, 118(Supplement C), 68–72.
10. FSPTCA 2009. Public Law No. 111-131 (June 2009).
11. M. R. Smith, B. Clark, F. Lüdicke, J.-P. Schaller, P. Vanscheeuwijck, J. Hoeng and M. C. Peitsch, Evaluation of the Tobacco Heating System 2.2. Part 1: description of the system and the scientific assessment program, Regul. Toxicol. Pharmacol., 2016, 81, S17–S26.
12. S. J. Sturla, A. R. Boobis, R. E. FitzGerald, J. Hoeng, R. J. Kavlock, K. Schirmer, M. Whelan, M. F. Wilks and M. C. Peitsch, Systems toxicology: from basic research to risk assessment, Chem. Res. Toxicol., 2014, 27, 314–329.
13. T. Hartung, E. van Vliet, J. Jaworska, L. Bonilla, N. Skinner and R. Thomas, Systems toxicology, ALTEX, 2012, 29, 119–128.
14. D. Krewski, D. Acosta Jr., M. Andersen, H. Anderson, J. C. Bailar 3rd, K. Boekelheide, R. Brent, G. Charnley, V. G. Cheung, S. Green Jr., K. T. Kelsey, N. I. Kerkvliet, A. A. Li, L. McCray, O. Meyer, R. D. Patterson, W. Pennie, R. A. Scala, G. M. Solomon, M. Stephens, J. Yager and L. Zeise, Toxicity testing in the 21st century: a vision and a strategy, J. Toxicol. Environ. Health, Part B, 2010, 13, 51–138.
15. M. D. Waters and J. M. Fostel, Toxicogenomics and systems toxicology: aims and prospects, Nat. Rev. Genet., 2004, 5, 936–948.
16. M. A. Davis, S. Eldridge and C. Louden, in Haschek and Rousseaux's Handbook of Toxicologic Pathology, ed. W. M. H. G. R. A. Wallig, Academic Press, Boston, 3rd edn, 2013, pp. 317–352,  DOI:10.1016/B978-0-12-415759-0.00010-8.
17. K. Boekelheide and M. E. Andersen, A mechanistic redefinition of adverse effects - a key step in the toxicity testing paradigm shift, ALTEX, 2010, 27, 243–252.
18. A. R. Iskandar, Y. Xiang, S. Frentzel, M. Talikka, P. Leroy, D. Kuehn, E. Guedj, F. Martin, C. Mathis and N. Ivanov, Impact assessment of cigarette smoke exposure on organotypic bronchial epithelial tissue cultures: a comparison of mono-culture and co-culture model containing fibroblasts, Toxicol. Sci., 2015, kfv122.
19. A. R. Iskandar, C. Mathis, F. Martin, P. Leroy, A. Sewer, S. Majeed, D. Kuehn, K. Trivedi, D. Grandolfo, M. Cabanski, E. Guedj, C. Merg, S. Frentzel, N. V. Ivanov, M. C. Peitsch and J. Hoeng, 3-D nasal cultures: Systems toxicological assessment of a candidate modified-risk tobacco product, Altex, 2017, 34(2), 23–49.
20. A. R. Iskandar, C. Mathis, W. K. Schlage, S. Frentzel, P. Leroy, Y. Xiang, A. Sewer, S. Majeed, L. Ortega-Torres, S. Johne, E. Guedj, K. Trivedi, G. Kratzer, C. Merg, A. Elamin, F. Martin, N. V. Ivanov, M. C. Peitsch and J. Hoeng, A systems toxicology approach for comparative assessment: Biological impact of an aerosol from a candidate modified-risk tobacco product and cigarette smoke on human organotypic bronchial epithelial cultures, Toxicol. in Vitro, 2017, 39, 29–51.
21. Kentucky Tobacco Research & Development Center, 3R4F Preliminary Analysis. Available at: http://ctrp.uky.edu/resources/pdf/webdocs/3R4F20Preliminary20Analysis.pdf.
22. J.-P. Schaller, D. Keller, L. Poget, P. Pratte, E. Kaelin, D. McHugh, G. Cudazzo, D. Smart, A. R. Tricker, L. Gautier, M. Yerly, R. Reis Pires, S. Le Bouhellec, D. Ghosh, I. Hofer, E. Garcia, P. Vanscheeuwijck and S. Maeder, Evaluation of the Tobacco Heating System 2.2. Part 2: Chemical composition, genotoxicity, cytotoxicity, and physical properties of the aerosol, Regul. Toxicol. Pharmacol., 2016, 81, S27–S47.
23. J. P. Schaller, J. P. Pijnenburg, A. Ajithkumar and A. R. Tricker, Evaluation of the Tobacco Heating System 2.2. Part 3: Influence of the tobacco blend on the formation of harmful and potentially harmful constituents of the Tobacco Heating System 2.2 aerosol, Regul. Toxicol. Pharmacol., 2016, 81(Suppl 2), S48–S58.
24. International Organization for Standardization, 1999, ISO 3402: 1999 Tobacco and tobacco products – Atmosphere for conditioning and testing, Available at: //iso.org/iso/home/store/catalogue_tc/catalogue_detail.htm?csnumber=28324.
25. Health Canada, Determination of Tar, Water, Nicotine and Carbon Monoxide in Mainstream Tobacco Smoke, Health Canada Test Method T-11, 1999.
26. S. Majeed, S. Frentzel, S. Wagner, D. Kuehn, P. Leroy, P. A. Guy, A. Knorr, J. Hoeng and M. C. Peitsch, Characterization of the Vitrocell(R) 24/48 in vitro aerosol exposure system using mainstream cigarette smoke, Chem. Cent. J., 2014, 8, 62.
27. C. Mathis, C. Poussin, D. Weisensee, S. Gebel, A. Hengstermann, A. Sewer, V. Belcastro, Y. Xiang, S. Ansari, S. Wagner, J. Hoeng and M. C. Peitsch, Human bronchial epithelial cells exposed in vitro to cigarette smoke at the air-liquid interface resemble bronchial epithelium from human smokers, Am. J. Physiol.: Lung Cell. Mol. Physiol., 2013, 304, L489–L503.
28. W. Fields, A. Maione, B. Keyser and B. Bombick, Characterization and Application of the VITROCELL VC1 Smoke Exposure System and 3D EpiAirway Models for Toxicological and e-Cigarette Evaluations, Appl. In Vitro Toxicol., 2017, 3, 68–83.
29. T. Muller and A. Hengstermann, Nrf2: friend and foe in preventing cigarette smoking-dependent lung disease, Chem. Res. Toxicol., 2012, 25, 1805–1824.
30. N. P. Castro, C. A. Osório, C. Torres, E. P. Bastos, M. Mourão-Neto, F. A. Soares, H. P. Brentani and D. M. Carraro, Evidence that molecular changes in cells occur before morphological alterations during the progression of breast ductal carcinoma, Breast Cancer Res., 2008, 10, R87.
31. F. Zanetti, A. Sewer, C. Mathis, A. R. Iskandar, R. Kostadinova, W. K. Schlage, P. Leroy, S. Majeed, E. Guedj and K. Trivedi, Systems toxicology assessment of the biological impact of a candidate modified risk tobacco product on human organotypic oral epithelial cultures, Chem. Res. Toxicol., 2016, 29(8), 1252–1269.
32. F. Zanetti, B. Titz, A. Sewer, G. Lo Sasso, E. Scotti, W. K. Schlage, C. Mathis, P. Leroy, S. Majeed, L. O. Torres, B. R. Keppler, A. Elamin, K. Trivedi, E. Guedj, F. Martin, S. Frentzel, N. V. Ivanov, M. C. Peitsch and J. Hoeng, Comparative systems toxicology analysis of cigarette smoke and aerosol from a candidate modified risk tobacco product in organotypic human gingival epithelial cultures: A 3-day repeated exposure study, Food Chem. Toxicol., 2017, 101, 15–35.
33. M. N. McCall, B. M. Bolstad and R. A. Irizarry, Frozen robust multiarray analysis (fRMA), Biostatistics, 2010, 11, 242–253.
34. M. Dai, P. Wang, A. D. Boyd, G. Kostov, B. Athey, E. G. Jones, W. E. Bunney, R. M. Myers, T. P. Speed, H. Akil, S. J. Watson and F. Meng, Evolving gene/transcript definitions significantly alter the interpretation of GeneChip data, Nucleic Acids Res., 2005, 33, e175.
35. G. K. Smyth, Linear models and empirical bayes methods for assessing differential expression in microarray experiments, Stat. Appl. Genet. Mol. Biol., 2004, 3, 3.
36. Y. Benjamini and Y. Hochberg, Controlling the false discovery rate: a practical and powerful approach to multiple testing, J. R. Stat. Soc. Series B: Methodol., 1995, 289–300.
37. F. Martin, A. Sewer, M. Talikka, Y. Xiang, J. Hoeng and M. Peitsch, Quantification of biological network perturbations for mechanistic insight and diagnostics using two-layer causal models, BMC Bioinformatics, 2014, 15, 238.
38. S. Boué, M. Talikka, J. W. Westra, W. Hayes, A. Di Fabio, J. Park, W. K. Schlage, A. Sewer, B. Fields, S. Ansari, F. Martin, E. Veljkovic, R. Kenney, M. C. Peitsch and J. Hoeng, Causal biological network database: a comprehensive platform of causal biological network models focused on the pulmonary and vascular systems, Database, 2015, bav030.
39. J. Hoeng, R. Deehan, D. Pratt, F. Martin, A. Sewer, T. M. Thomson, D. A. Drubin, C. A. Waters, D. de Graaf and M. C. Peitsch, A network-based approach to quantifying the impact of biologically active substances, Drug Discovery Today, 2012, 17, 413–418.
40. T. M. Thomson, A. Sewer, F. Martin, V. Belcastro, B. P. Frushour, S. Gebel, J. Park, W. K. Schlage, M. Talikka, D. M. Vasilyev, J. W. Westra, J. Hoeng and M. C. Peitsch, Quantitative assessment of biological impact using transcriptomic data and mechanistic network models, Toxicol. Appl. Pharmacol., 2013, 272(3), 863–878.
41. B. S. Carvalho and R. A. Irizarry, A framework for oligonucleotide microarray preprocessing, Bioinformatics, 2010, 26, 2363–2367.
42. W. Huber, V. J. Carey, R. Gentleman, S. Anders, M. Carlson, B. S. Carvalho, H. C. Bravo, S. Davis, L. Gatto, T. Girke, R. Gottardo, F. Hahne, K. D. Hansen, R. A. Irizarry, M. Lawrence, M. I. Love, J. MacDonald, V. Obenchain, A. K. Oles, H. Pages, A. Reyes, P. Shannon, G. K. Smyth, D. Tenenbaum, L. Waldron and M. Morgan, Orchestrating high-throughput genomic analysis with Bioconductor, Nat. Methods, 2015, 12, 115–121.
43. R Core Team, R: A language and environment for statistical computing, R Foundation for Statistical Computing, Vienna, Austria. Available at: http://www.R-project.org/, 2013.
44. A. Kauffmann, R. Gentleman and W. Huber, arrayQualityMetrics–a bioconductor package for quality assessment of microarray data, Bioinformatics, 2009, 25, 415–416.
45. B. M. Bolstad, R. A. Irizarry, M. Astrand and T. P. Speed, A comparison of normalization methods for high density oligonucleotide array data based on variance and bias, Bioinformatics, 2003, 19, 185–193.
46. A. Kozomara and S. Griffiths-Jones, miRBase: annotating high confidence microRNAs using deep sequencing data, Nucleic Acids Res., 2013, 42(D1), D68–D73.
47. Affymetrix miRNA QCT Tool, Affymetrix miRNA QCTool: User's Guide. Available at: http://media.affymetrix.com/support/downloads/manuals/mirna_qctool_user_manual.pdf, 2011.
48. M. E. Ritchie, B. Phipson, D. Wu, Y. Hu, C. W. Law, W. Shi and G. K. Smyth, Limma powers differential expression analyses for RNA-sequencing and microarray studies, Nucleic Acids Res., 2015, 43(7), e47.
49. U. Bodenhofer, A. Kothmeier and S. Hochreiter, APCluster: an R package for affinity propagation clustering, Bioinformatics, 2011, 27, 2463–2464.
50. W. Stannard and C. O'Callaghan, Ciliary function and the role of cilia in clearance, J. Aerosol Med., 2006, 19, 110–115.
51. A. E. Tilley, M. S. Walters, R. Shaykhiev and R. G. Crystal, Cilia dysfunction in lung disease, Annu. Rev. Physiol., 2015, 77, 379–406.
52. P. J. Stanley, R. Wilson, M. A. Greenstone, L. MacWilliam and P. J. Cole, Effect of cigarette smoking on nasal mucociliary clearance and ciliary beat frequency, Thorax, 1986, 41, 519–523.
53. N. A. Cohen, S. Zhang, D. B. Sharp, E. Tamashiro, B. Chen, E. J. Sorscher and B. A. Woodworth, Cigarette smoke condensate inhibits transepithelial chloride transport and ciliary beat frequency, Laryngoscope, 2009, 119, 2269–2274.
54. D. Kuehn, S. Majeed, E. Guedj, R. Dulize, K. Baumer, A. Iskandar, S. Boue, F. Martin, R. Kostadinova and C. Mathis, Impact assessment of repeated exposure of organotypic 3D bronchial and nasal tissue culture models to whole cigarette smoke, J. Visualized Exp., 2015 , JoVE, (96).
55. M. Talikka, R. Kostadinova, Y. Xiang, C. Mathis, A. Sewer, S. Majeed, D. Kuehn, S. Frentzel, C. Merg and M. Geertz, The response of human nasal and bronchial organotypic tissue cultures to repeated whole cigarette smoke exposure, Int. J. Toxicol., 2014, 33(6), 506–517.
56. X. Li, In vitro toxicity testing of cigarette smoke based on the air-liquid interface exposure: A review, Toxicol. in Vitro, 2016, 36, 105–113.
57. J. Hoeng, M. Talikka, F. Martin, A. Sewer, X. Yang, A. Iskandar, W. K. Schlage and M. C. Peitsch, Case study: the role of mechanistic network models in systems toxicology, Drug Discovery Today, 2014, 19, 183–192.
58. I. Wilk-Zasadna, C. Bernasconi, O. Pelkonen and S. Coecke, Biotransformation in vitro: An essential consideration in the quantitative in vitro-to-in vivo extrapolation (QIVIVE) of toxicity data, Toxicology, 2015, 332, 8–19.
59. R. W. Lewis, R. Billington, E. Debryune, A. Gamer, B. Lang and F. Carpanini, Recognition of adverse and nonadverse effects in toxicity studies, Toxicol. Pathol., 2002, 30, 66–74.
60. D. N. Watkins, D. M. Berman, S. G. Burkholder, B. Wang, P. A. Beachy and S. B. Baylin, Hedgehog signalling within airway epithelial progenitors and in small-cell lung cancer, Nature, 2003, 422, 313–317.
61. A. Sriperumbudur, M. Breitzig, R. Lockey and N. Kolliputi, Hedgehog: the key to maintaining adult lung repair and regeneration, J. Cell Commun. Signal., 2017, 11, 95–96.
62. T. Peng, D. B. Frank, R. S. Kadzik, M. P. Morley, K. S. Rathi, T. Wang, S. Zhou, L. Cheng, M. M. Lu and E. E. Morrisey, Hedgehog actively maintains adult lung quiescence and regulates repair and regeneration, Nature, 2015, 526, 578–582.
63. K. J. A. Davies, Adaptive Homeostasis, Mol. Aspects Med., 2016, 49, 1–7.
64. P. V. Nazarov, S. E. Reinsbach, A. Muller, N. Nicot, D. Philippidou, L. Vallar and S. Kreis, Interplay of microRNAs, transcription factors and target genes: linking dynamic expression changes to function, Nucleic Acids Res., 2013, 41, 2817–2831.
65. M. Yamada, N. Fujino and M. Ichinose, Inflammatory responses in the initiation of lung repair and regeneration: their role in stimulating lung resident stem cells, Inflammation Regener., 2016, 36, 15.
66. J. R. Roberts, G. D. Perkins, T. Fujisawa, K. A. Pettigrew, F. Gao, A. Ahmed and D. R. Thickett, Vascular endothelial growth factor promotes physical wound repair and is anti-apoptotic in primary distal lung epithelial and A549 cells, Crit. Care Med., 2007, 35, 2164–2170.
67. R. B. Schlesinger, Biological Disposition of Airborne Particles: Basic Principles and Application to Vehicular Emissions in Air Pollution, the Automobile, and Public Health, ed. A. Y. Watson, R. R. Bates and D. Kennedy, National Academies Press (US), Washington (DC), 1998, Available from: https://www.ncbi.nlm.nih.gov/books/NBK218161/.
68. N. R. Labiris and M. B. Dolovich, Pulmonary drug delivery. Part I: Physiological factors affecting therapeutic effectiveness of aerosolized medications, Br. J. Clin. Pharmacol., 2003, 56, 588–599.
69. S. Steiner, S. Majeed, G. Kratzer, G. Vuillaume, J. Hoeng and S. Frentzel, Characterization of the Vitrocell(R) 24/48 aerosol exposure system for its use in exposures to liquid aerosols, Toxicol. in Vitro, 2017, 42, 263–272.
70. ISO (2012), Routine analytical cigarette-smoking machine – Definitions and standard conditions, International Organization for Standardization, Available at: https://www.iso.org/obp/ui/#iso:std:iso:3308:ed-5:v1:en.
71. L. Pichelstorfer, W. Hofmann, R. Winkler-Heil, C. U. Yurteri and J. McAughey, Simulation of aerosol dynamics and deposition of combustible and electronic cigarette aerosols in the human respiratory tract, J. Aerosol Sci., 2016, 99, 125–132.

---

Footnotes

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c7tx00152e

‡ These authors contributed equally to this work.

---