Multivariate Optimization of Tenax TA-Thermal Extraction for Determining Gaseous Phase Organophosphate Esters in Air Samples

Suitable conditions for thermal extraction of semi-volatile organic compounds have largely been arrived at by univariate optimization or based on the recommendations provided by the manufacturers of the extraction equipment. Herein, we demonstrated the multivariate optimization of Tenax TA–thermal extraction for determining organophosphate esters in the gas phase fraction of air samples. Screening and refining experiments were performed using the eighth fraction factorial and Box-Behnken designs, respectively, and satisfactory models were obtained. Subsequently, the process was optimized by Derringer’s desirability function and the global desirability was 0.7299. Following optimization, the analytes were desorbed at 290 °C for 10 minutes at a helium flow of 95 mL min−1, with the transfer line set at 290 °C. The analytes were then cryofocused at 20 °C and then cryodesorbed into the chromatographic column at 295 °C for 6 minutes. Method validation exhibited high linearity coefficients (>0.99), good precision (CV < 14%) and low detection limits (0.1–0.5 ng m−3). The method was tested by pumping 0.024 m3 of real indoor environment air through Tenax TA sorbent tubes. Furthermore, with multivariate optimization, analysis time and other resources were significantly reduced, and information about experimental factor interaction effects was investigated, as compared to the univariate optimization and other traditional methods.

Results screening factorial design. The investigated factors and their ranges have been given in the methods section (Screening design). The design matrix together with the response factors' results are presented in Table S1. The runs at the center were repeated twice (six replicates in total) in this case to ensure that the repeatability was within the normal dispersion range, and in terms of the coefficient of variation (CV), repeatability ranged from 4.0% for TCEP to 16.6% for EHDPP. Analysis of variance (ANOVA) results of the screening design experiments are presented in Table 1 in form of the relative significance (with their signs) of the main effects associated to each factor as well as the quadratic interaction of desorption flow. The positive sign means that the change of the experimental factor from the low level to high level increases the magnitude of the response factor, while a negative sign means that the change in the experimental factor from the low level to high level decreases the magnitude of the response factor. Desorption flow rate, TDS transfer line temperature, desorption hold time, desorption temperature and cryodesorption time had non-significant effects on the chromatographic peak areas of the target compounds. Cryofocusing temperature (E), cryodesorption temperature (F) and the quadratic effect of desorption flow rate (AA), had significant effects on the peak areas of at least one of the target compounds and their effects were therefore, evaluated by response surface methodology (RSM). The rest of the main experimental factors were not considered for optimization but were set at the central points when the BBD experiments were being performed.
Response Surface Methodology (RSM). The letter codes for the studied factors in the BBD are A, B and C, respectively, for desorption flow rate, cryofocusing and cryodesorption temperatures. All experiments were performed in a randomized manner to minimize the bias effects of uncontrolled variables 25 . Three central points were incorporated to estimate the experimental error that is independent of the fitted model and to ensure that the variability was normally dispersed. The CVs of the results obtained from the central point experimental runs ranged from 6.8% for TEP to 19.6% for TPHP. The plot and experimental domain matrix for this design are shown in the Fig. S2. Additionally, the design matrix (with actual experimental factors' values) together with the response factors' results are presented in Table S2. The regression coefficients of the second order polynomial models (equation 6) expressing the relationship between the three experimental factors and response factors were computed based on the generated results and presented in Table 2. Analysis of variance (ANOVA) was used to evaluate the appropriateness of the model through the coefficients of determination as well as their corresponding optimization process. For each of the responses, the partial desirability values were computed and are presented in Table 3, together with the global desirability that was obtained by optimizing the process simultaneously by searching from all design points and vertices. Under the aforementioned optimization considerations, the optimum conditions for the factors that were investigated in this design were as follows: desorption flow rate = 95.1069 mL min −1 , cryofocusing temperature = 31.71 °C and cryodesorption temperature = 294.698 °C with a global desirability of 0.7299 (Table 3).
Chromatographic Separation. The compounds were completely separated under the chromatographic conditions that were used in our previous study 8 , and the results are presented in Table 4. Furthermore, the total ion chromatograms (TIC) for the field blank tube, intermediate calibration concentration and one of the real air samples are presented in Fig. S1.
Method performance and validation. For purposes of calibration, five concentration levels were prepared from 50 to 2000 pg corresponding to 2.08 to 83.33 ng m −3 calculated for a sample volume of 0.024 m 3 and experiments at each level, were performed in triplicate. Table 4 shows the results/values of statistical and analytical parameters obtained for each of the target OPEs.
Application to real indoor air samples. We collected air samples from six different indoor environments (office (n = 3), bedroom (n = 2) and car (n = 1)). For all sampling campaigns, the ambient temperature was 23.2 ± 2.3 °C. As mentioned in the sampling section (under methods and materials), a QFF was used to seal the sampling end of the sorbent tube to prevent particles from entering the tube. Although some studies have indicated that some OPEs can get adsorbed on a QFF during sampling 27 , the sampling time of two hours was too short to influence the results significantly due the adsorption artifact 28 . The results are presented in Table 5. TPP, TNBP, TCEP, TCIPP, TPHP and TEHP were quantified in at least one of the sampled environments. Our results are in agreement with those obtained in some previous studies in which at least one of the target OPEs in this study were investigated in the gaseous fraction of air samples 5,29 .

Discussion
Full factorial design (FD) structures have benefits of efficiency, as well as allowing the experimenter to study the main factors and their interaction effects on process response factors. To allow the experimenter get information on the main effects and lower order interactions without having to run the full factorial design (FD), modification to fractional factorial designs (FFD) is performed. In factorial designs, each experimental factor is usually studied at two levels (low (−1) and high (+1)) although experimental studies at the central point (0) are also possible 25 . As a consequence, factorial designs are gaining interest in preliminary studies or in the initial steps of optimizing a given process 30,31 .
The fundamentals, advantages, limitations and applications of the response surface methodology (RSM) based on the BBD in optimizing chemical processes, were comprehensively reviewed by Ferreira et al. 24 . The analysis of variance (ANOVA) results (  Table 1. Experimental domain and the relative significance (with their sign) of the main effects associated to each factor as well as the quadratic interaction of desorption flow. +++ or −−− indicate a statistically significant (p < 0.05) positive or negative effect. ++ or −−indicate a positive or negative effect that was close to statistical significance. + or − indicate a positive or negative effect that was far from reaching statistical significance.
www.nature.com/scientificreports www.nature.com/scientificreports/ determination (r 2 ) values ranging from 0.7661 for TNBP to 0.9724 for TBOEP. The adjusted coefficient of determination (r adj 2 ) values were also evaluated and are presented in Table 2. In general, the values of r 2 and r adj 2 indicated that the modelled response factors cover the point cloud of the experimental results satisfactorily. Besides, for each of the responses, the Durbin-Watson (DW) statistic had a p-value greater than 0.05, showing no indication of serial autocorrelation in the residuals at the 95% confidence level.
As depicted in Fig. 1(a), the peak area of TEP increased when desorption flow increased together with the simultaneous increase in cryofocusing temperature that was apparent below −20 °C. Thereafter, peak area reduced with increasing cryofocusing temperature and desorption flow as the linear interaction of desorption flow and cryofocusing temperature and quadratic interaction of cryofocusing temperature which have negative effects became important. The reason could be that the more volatile TEP is less effectively trapped at higher temperatures and flow rates. In Fig. 1(b), the simultaneous increase in desorption flow and cryofocusing temperature increases the peak area of TPP. However, after 20 °C, the interaction of the two has a negative effect on peak area although not significant. In Fig. 1c to j, the peak areas increase with the simultaneous increase in desorption flow and cryofocusing temperatures. Moreover, the interactions that are significant for example AB for TCEP, TBOEP and EHDPP, have positive effects on the peak areas. AA has a significant negative effect on the peak area of EHDPP as seen clearly in Fig. 1(j).
We observed that for all the compounds, except TEP (the most volatile), the peak areas generally increased with increasing cryofocusing temperature. Leon et al. 19 attributed the same observation to the fact that, at lower cryofocusing temperatures, the transfer of the trapped analytes from the PTV injector to the chromatographic column may be ineffective due to frosting and water condensation or any other physical phenomena. Additionally, the peak areas were observed to increase with increase in desorption flow rate probably due to the adequate transfer of the analytes to the PTV injector. However, the peak areas are observed to decrease at higher desorption flow rates probably due to inadequate trapping of the analytes though for most of the compounds, the peak areas remained almost constant past the optimum desorption flow rate. It is worth noting that in RSM www.nature.com/scientificreports www.nature.com/scientificreports/ experiments, all variables are important, as there has been a preliminary screening to finding important variables, prior to exploring the surface 25 . Therefore, we never excluded any effects in the derivation of the second order polynomial models that give the relationship between the peak areas and the experimental factors.
Regarding optimization, the TE conditions would be considered as optimum if the peak area yields approached the maximum values. As a matter of fact, the goal of the process was to maximize the response factors. As seen in Fig. 1(a-j), the peak areas for different compounds are maximized by different conditions. Therefore, to obtain the optimal conditions (for the three factors evaluated) that are suitable for the extraction of all the compounds, we used the desirability function that was suggested by Derringer and Suich 32 , to simplify the response factor (peak area) for each compound. According to this function, each individual response "Y i " was transformed into a desirable dimensionless value, partial desirability (d i ). The scale of the desirability function  Table 3. Modelled and adopted optimum conditions for the significantly influential experimental factors, partial and global desirability of the peak areas, the predicted areas (and their 95% confidence limits) and the experimentally obtained peak areas.
www.nature.com/scientificreports www.nature.com/scientificreports/ ranged from d = 0 for a completely undesirable response to d = 1 for a fully desired response. The 10 individual partial desirability functions were then weighted into a single composite response known as the global desirability function, D.
However, considering the operational parameters of the equipment, desorption flow and cryodesorption temperature were modified to 95 mL min −1 and 295 °C. Furthermore, by exploring the desirability change with decreasing cryofocusing temperature from the optimum setting of 31.71 °C, keeping the desorption flow at 95 mL/min ( Fig. 1(k)), we adopted 20 °C as a compromised optimum temperature 25 for trapping the analytes without affecting the value of the global desirability. Besides, at these conditions of desorption flow rate and cryofocusing temperature, it is clear from the response surface plots ( Fig. 1(a to j) that the peak areas of the analytes remained within the desirable ranges. Additionally, at this temperature, the consumption of liquid nitrogen (coolant for the PTV) could be minimized 23 since higher cryofocusing temperatures require a longer time for the CIS to equilibrate thereby consuming more liquid nitrogen prior to actual sample extraction. Cryofocusing at 20 °C has also been shown to produce good responses in previous studies 19,23 in which the TE of SVOCs from water samples was investigated.
Under the modified optimal conditions for desorption flow, cryofocusing temperature and cryodesorption temperature, we fixed the insignificantly influential experimental factors (desorption temperature, TDS transfer line temperature, desorption time and cryodesorption time) as elucidated by the results of the screening design, at their central and high points, respectively, and in each case, the experiments were run in triplicate. We observed that the peak areas of the less volatile compounds (TDCIPP, TBOEP, TPHP and EHDPP) reduced at high levels of insignificant parameters. This could be attributed to the partial thermal decomposition that these compounds may undergo at higher desorption and cryodesorption temperatures 7 . By performing a paired samples t -test, it was revealed that the difference between the peak areas obtained at both the central and high points for the insignificant factors was not significant (p > 0.05). We therefore, used the central points for the insignificantly influential factors (to minimize thermal decomposition of the less volatile analytes), together with the modified optimal conditions of the significantly influential factors to test for the thermal extraction efficiency and validate the method.   www.nature.com/scientificreports www.nature.com/scientificreports/ The experimentally obtained peak area results when the insignificant factors were set are the central points are presented in Table 3. By using the predicted mean peak areas, prediction 95% confidence intervals were constructed and presented in Table 3. The experimentally obtained/validated mean peak areas (Table 3) generally fell  with these confidence limits for TEP, TPP, TNBP, TCIPP and TPHP. For the rest of the compounds, the experimentally determined mean peak areas fell slightly below the lower limits of the 95% confidence intervals. From the results of a t -test, there was no significant difference (p > 0.05) between these mean peak values and their corresponding 95% confidence interval lower limits. This implies that the model predictions could be accepted for these compounds. Further, a scatter plot of the validated/experimentally obtained areas against the model predicted areas was generated and is presented in Fig. S3. The high value of the coefficient of determination (r 2 = 0.9906) indicates a very strong agreement between the predicted and experimental results as well as the high accuracy of the proposed second order polynomial models. Moreover, the repeatability (CV%, Table 3) of the optimized procedure under the final modified optimal conditions was satisfactory, ranging from 6.7% for TCIPP to 12.3% for TBOEP.
The thermal desorption efficiency investigations were carried out by analyzing 3 separately spiked sorbent tubes, followed by desorbing the tubes for the second time. Desorption efficiency was calculated as the ratio of the peak area obtained on the first desorption to the total peak area obtained from the two desorptions and expressed as a percentage. The range of the desorption efficiency (%± CV) was from 99.8 ± 0.02% for TPHP to 100.0 ± 0.05% for TPP. Notably, these results showed no signs of carryover for all the target compounds.
Using the optimized conditions, we tested the performance of the method by evaluating appropriate validation parameters. The coefficient of determination (r 2 ) values obtained ranged from 0.9942 for TEHP to 0.9998 for TEP, indicating that the responses were highly linear in the concentration range that was evaluated for all compounds. The values obtained were acceptable and in agreement with those obtained in the previous studies in which TE was investigated for flame retardants 7,8 . In addition, the linearity of the calibration curves was tested using Pearson correlation coefficient at a probability of 99%. In all cases, this test was passed, i.e., all the linearity coefficients were significant with 99% confidence. We also evaluated the relative response factor (RFF) CV values all of which were less than 10% for all calibration data points ( Table 4, column 8).
Regarding the sensitivity of the method, the calculated instrumental LODs and LOQs ranged from 2.7 to 12.0 and 9.0 to 44.0 pg, respectively. These values correspond to LODs of 0.11 to 0.50 ng m −3 and LOQ of 0.38 and 1.83 ng m −3 , respectively basing on the sampling volume of 0.024 m 3 . It is noteworthy to mention that our method is characterized by at least 6 times lower LODs than those obtained by Aragón et al. 5 . However, they are comparable with those obtained by Ramírez et al. 33 but slightly lower than those obtained by Lazarov et al. 7 probably due to differences in sample volumes. Table 4 summarizes the LOD and LOQ results obtained for all the target compounds. The results obtained demonstrated that the optimized method was effective in determining the target OPEs especially at low concentration levels.
The accuracy of the method was investigated at three calibration concentration levels: the lowest, intermediate and highest. Intra-day precision (n = 5) expressed in terms of coefficient of variation (CV) ranged from 1.5 to 9.9%, 2.2 to 13.2% and 0.7 to 9.8% at the lowest, intermediate and highest levels, respectively (Table 4). Generally, our results on intra-day precision of the method were in agreement with those obtained in some previous TE studies on SVOCs in which the CV values were less than 20% 5,7,8,33,34 . Moreover, the obtained CV values of less than 14%, indicated that the optimized procedure was characterized by good repeatability of the results. The recoveries of the target compounds ranged from 82% for EHDPP to 133% for TPHP, 90% for TEP to 115% for TDCIPP, and 96% for TEP to 124% for TBOEP, at the lowest, intermediate and highest calibration concentration levels, respectively. Generally the recovery range of our results fell with the range of 70 to 120% that is acceptable during the validation of analytical procedures depending on matrix complexity 35 . A few exceptions that were observed at the lowest and highest calibration concentration levels, were not significantly different (p > 0.05, student t -test) from the maximum acceptable limit (120%). These values affirm that the optimized method is accurate for the determination of the 10 OPEs in air samples. Furthermore, the obtained values indicated that the TE procedure was highly efficient to isolate gaseous OPEs present in the air sample matrix and hence applicable in the determination of these compounds in air samples collected from different environments. The results of inter-day (intermediate) precision (n = 5) CVs ranged from 4.4 to 11.8%, 2.5 to 12.5% and 3.2 to 12.5% at the lowest, intermediate and highest levels, respectively (Table 4). Our results agree with those obtained in some previous TE studies on SVOCs 8,33 in which CVs less than 20% were obtained. The results indicated that the optimized procedure was characterized by high reproducibility, i.e., the results obtained were not sensitive to changes or variations in experimental times. It is worth noting that the precision values obtained for all the target compounds were less than 20% that is usually the acceptable limit for analytical method validation 36 .
With regard to the selectivity of the method, we compared the chromatograms obtained for the field blank, the intermediate calibration level and the real air sample (car) as shown in Fig. S1. There were no peak interferences observed at the retention times of the analytes in the SIM mode in the field blank except for TCEP. Some other peaks that were observed in the field blank tube are attributed to the internal standards and thermal decomposition which could result into the formation of other fragments 11 , or any other physical phenomena. These findings suggested that the spectrometric conditions that we used 8 , accorded high selectivity to the method for the target analytes. At all levels, all blank detections were less than the LOQs exept for TCEP, of which appropriate blank corrections were performed during the treatment of data obtained from real air sampling. Although some traces of the internal standards could be observed in the post conditioning analyses, perhaps due to some cold spots in the TDS, the percentage was as low as 4.3% and 2.6% for d 27  The aim was to study the effect of changing the cryofocusing temperature by 10% past the modified optimum, on the peak area of EHDPP. Substituting for A = 95 mL min −1 and C = 295 °C, equation 2 becomes; And when x = 20 °C, y = 2861758.54 Δy = 26865.846 which corresponds to a 1% increase in peak area. This result shows that the method is not sensitive to random variations in the experimental conditions. It is therefore, highly rugged.
Concerning the breakthrough tests, after analyzing the two tubes for each test cycle, the concentration, in duplicate of each compound, obtained from the inlet tube was expressed as a percentage of the corresponding concentration obtained by analyzing tubes spiked with 1000 pg of the standard solution. The relative percentage recovery of all the compounds was greater than 95% up to when a volume of 0.024 m 3 was passed through the tubes that were being tested. When 0.048 m 3 was pumped through the tubes, the recovery of the compounds ranged from 50% for TEP to 90% for TNBP. We therefore considered 0.024 m 3 as the appropriate sampling volume in real sampling applications. Moreover, breakthrough is considered to have occurred when the recovery in the inlet tube is less than 95% 9 .
Our results were compared with other studies as presented in Table 6. The optimum value for desorption flow in our study was 95 mL min −1 and is comparable with 100 mL min −1 which was used in the extraction of personal care products 33 but which was used based on the supplier's recommendation. Besides, our result was different from those used by Aragón et al. 5 and Matsiko et al. 8 because in these studies, univariate optimization was used. In Lazarov et al. 7 , no information is available on how the flow of 20 mL min −1 was selected. Our final desorption temperature value of 290 °C was comparable with those that were used in previous studies 5,7,8,33 , although in these studies, these values were selected based on the maximum temperature that is recommended for the sorbents that were used for extraction except for 300 °C which was the optimized value obtained by Matsiko et al. 8 . The final desorption time that we used in this study was comparable with those optimized in previous studies, whether optimized 5,8 or selected 7,33 . Information on the TDS transfer temperature was not provided in some studies 5,7,33 , but the final value that we used in this study is comparable with the one optimized by Matsiko et al. 8 . For cryofocusing temperature, the modified optimum value in our study was different from the values used in the previous studies because the used values were either selected 5,7,33 based on previous studies or arrived at by univariate optimization 8 . The optimum value for cryodesorption temperature in this study was different from those used in other studies due to the fact in some studies, this value was selected based on the maximum allowable value for the trap sorbent 5,33 or univariate optimization was used 7,8 . The same reason applies to the difference in cryodesorption hold time. The detection and quantification limits in ng m −3 were calculated by dividing corresponding values in pg by the sampling volume (0.024 m 3 ) and were comparable with those obtained in the previous studies 5,7,33 on TE for determining SVOCs in air samples (Table 6). It is worth noting that the MQLs obtained in this study are also comparable with the results obtained in some studies in which solid phase extraction (SPE) was used for indoor active sampling of OPEs [37][38][39] . This indicates that the sensitivity of our method is satisfactory for the determination of gaseous phase OPEs in the indoor environments even at low concentration levels.  www.nature.com/scientificreports www.nature.com/scientificreports/ In conclusion, this work demonstrated the applicability of multivariate modelling in optimizing Tenax TAthermal extraction for the determination of 10 OPEs in the gas phase fraction of air samples. With this optimization approach we were able to investigate a wider experimental domain, i.e., effect of the main factors, linear interactions and quadratic effects of the experimental factors on the magnitude of the response factors. Additionally, this optimization approach coupled with thermal extraction have proved to be time saving, and the amount of other resources such as the standard reagents, liquid nitrogen and helium are tremendously reduced. Without doubt, the versatility of multivariate optimization coupled with thermal extraction renders it a fast and cheap alternative to the commonly used univariate optimization and other multi-step extraction procedures such as accelerated solvent extraction. The relevance of the method was tested by collecting air samples from six indoor environments and TCEP was quantified in all the samples. Although this method is characterized by low method detection limits (0.02-0.5 ng m −3 ), the values are higher than those obtained in our previous study 8 , probably due to different sampling volumes. Noteworthy to mention, this method is characterized by low sampling volumes which may limit the detection of some target analytes, especially those that mainly partition in the particulate phase, in cases where the concentrations are at trace levels. Notwithstanding, the applicability multivariate optimization approach followed by TE in multi-residue analysis of SVOCs will be of interest.

Materials and Methods
Reagents and chemical standards. Acetone (pesticide grade) was obtained from Fisher Scientific, Fair Lawn, NJ, USA. Toluene (99.7% pure) was purchased from J.T. Baker, USA. Details of the chemical and internal standard solutions of the target OPEs have been provided in our previous study 8 . Stock solutions were prepared in toluene and appropriate dilutions were made with acetone. Additionally, the physical properties of the target compounds are presented in Table S4. sampling cartridges. Standard glass sorbent cartridges (180 mm, length 6 mm, outer diameter and 4 mm inner diameter, Gerstel), each packed with approximately 60 mm/approximately 180 mg of Tenax TA adsorbent, were used for adsorption/sampling and then thermal extraction, prior to chromatographic analysis. Tenax TA is a polymeric material based on 2, 6-diphenyl-p-phenylene oxide that is used as an adsorbent material in various analytical applications. Its physical properties with regard to topography, crystal structure and thermal stability were evaluated by Alfeeli et al. 40 . Before use, the sorbent packed tubes were thermally cleaned in a Gerstel Tube conditioner for eight hours at a temperature of 320 °C and under nitrogen gas flowing at 75 mL min −1 . Subsequent pre-cleaning was carried out for 2 hours at the same conditions of temperature and gas flow rate. Thereafter, the clean tubes were placed in storage tubes which were then capped with long-term storage plastic caps combined with PTFE ferrules and then stored in airtight sealable plastic jars, which were then kept under a nitrogen atmosphere to prevent contamination from the immediate surroundings. tube loading and desorption. To optimize, calibrate and validate the method for the TE procedure under the current study, it was necessary to spike the sorbent tubes with standard solutions. For all the calibration concentration levels used at each stage, 1 μL of the standard solution in acetone was spiked onto the sorbent tube using a 10 μL GC/MS auto sampler syringe and thereafter, purged with 99.99% nitrogen gas at a flow rate of 75 mL min −1 for 2 minutes. The liquid mixture vapourized in this gas flow constitution and the analytes were adsorbed onto the sorbent in vapour phase. Purging was carried out to ensure that the optimization and validation procedures were consistent with real air sampling. Acetone was chosen as the solvent because of its weak retention characteristics on Tenax TA 40 . Notably, 2 minutes of nitrogen gas flow (150 mL) was adequate to dry the solvent while making sure that the analytes of interest were not expelled out of the tube. The spiked tubes were then placed in storage tubes which were capped with long-term storage plastic caps combined with PTFE ferrules, wrapped in aluminium foil and stored at −20 °C for about 4 hours for the sorbent bed to equilibrate, prior to extraction. A standard two stage desorption strategy was used for analysis. experimental Design/statistical treatment. All experiments in the screening and refining/modelling procedures were designed using Statgraphics Centurion 18 (StatPoint Technologies, Inc) software. The goal of the designs used in this study was to maximize the response factors (compound specific chromatographic peak areas). Other statistical analyses were performed by Microsoft excel, IBM ® SPSS ® Statistics Version 21 and Origin ® Pro Version 9.1. screening. We planned screening experiments using the eighth fractional factorial design (2 K−3 , K = 7) with 4 center points to determine which factors affected the compounds' peak areas significantly. The seven factors that were investigated, together with their ranges are as follows; desorption flow rate (20-100 mL min −1 ), desorption temperature (260-320 °C), thermal desorption system (TDS) transfer line temperature (260-320 °C), desorption hold time (5-15 min), cryofocusing temperature (−100-+60 °C), cryodesorption temperature (260-320 °C) and cryodesorption hold time (2-10 min). These factors and their ranges were considered basing on preliminary experiments, previous studies [5][6][7][8][21][22][23] , boiling points of the analytes (Table S4), and the thermal stability of Tenax TA 40 .
A second order polynomial model was selected in which all the main factor effects and the quadratic effect of desorption flow rate were automatically included in the design. By using the forward algorithm, 19 desired runs were selected, i.e., 16 runs at the factorial points and 3 runs at the center of the design. The central point was repeated twice (six replicates in total) to ensure that the repeatability was within the normal dispersion range. The experiments were ran in a random manner in order to minimize the effect of uncontrolled variables 25 . These experiments were carried out in duplicate by spiking 1 μL of the standard solution containing 1000 pg of each target compound on separate sorbent tubes. (2019) 9:3330 | https://doi.org/10.1038/s41598-019-40119-2 www.nature.com/scientificreports www.nature.com/scientificreports/ Modelling and optimization. The significantly influential factors as revealed by the results of the screening experiments, were used to develop second order polynomial models for evaluating the response factors. The experiments were planned using the Box -Behnken design (BBD) with 15 experimental runs (12 at factorial points and 3 at the center). The significance level for each experimental factor in the response factor mathematical models was evaluated by the analysis of variance (ANOVA). These experiments were carried out in duplicate by spiking 1 μL of 2000 pg standard solution on separate sorbent tubes. The central point was also repeated twice to ensure that repeatability was within the normal dispersion range. In the selected second order models (equation 6), the main factors, linear interaction and quadratic effects were evaluated. 0  1  2  3  1 1  2  12  13  22  2  23  33  2 where, Y is the single response (compound specific peak area) to be modelled, β is the regression coefficient, A, B and C are the experimental factors and ε is the experimental error. Subsequently, the models, together with Derringer's desirability function, were used to determine the optimal conditions for the TE procedure under the current study.

Method validation and Quality Assurance/Control (QA/QC).
To test for the linearity of the method, the tubes were spiked separately in triplicate with 1 μL of each of the calibration concentration solutions (50-2000 pg) containing 100 pg of the internal standards 9 . Calibration was then done by plotting the area ratio, A A ( / ) i sur (where A i is the peak area of the target analyte and A sur is the peak area of the surrogate standard) against the concentration of the target analyte. Due to the unavailability of certified reference materials, the accuracy of the method in terms of recovery, method repeatability and intermediate precision was tested by spiking five tubes separately with the lowest, intermediate and highest concentration levels. Repeatability was tested by performing the experiments on the same day, while intermediate precision was tested by performing the experiments for five consecutive days, one experiment per day for each of the three calibration levels. In order to ascertain whether the target compounds could be present in the air samples, we examined the limits of detection and quantitation (LOD and LOQ). These were calculated using the formulae; 3σ and 10σ where σ is the standard deviation of the analytical results obtained by analyzing 5 desorption tubes spiked with the lowest calibration concentration.
QA/QC was ensured at different levels. Post conditioning analyses were carried out by analysing six pre-cleaned tubes without spiking them with the internal standards. We also investigated post conditioning analyses when the tubes were spiked with the internal standards. Laboratory and field blanks were carried out by connecting six tubes separately to the pump and then exposing them to the laboratory air and an office environment air, respectively, in the passive diffusive mode (without starting the pump) for one minute 34 . They were then capped, wrapped in a clean aluminuim foil and kept at −20 °C until analysis. Breakthrough volume tests were carried out by connecting two pre-cleaned sorbent tubes in series with a silicone tube, with the inlet tube spiked with 1000 pg of the standard solution. Volumes of clean laboratory air ranging from 0.003 to 0.048 m 3 were pumped through the tubes at a rate of 200 mL min −1 . At the end of each experiment, the two tubes were analysed separately. For each volume, the experiments were run in duplicate.
sampling. The air samples were collected based on active air sampling with a variable flow air pump (QC-IC portable membrane pump, 50-500 mL min −1 , Beijing Labour Protection Science Research Institute, CN). Air samples were pumped through the preconditioned sorbent tubes at a flow rate of 200 mL min −1 for 2 hours with a total volume of 0.024 m 3 of air collected. Prior to sampling, we calibrated two pumps using a soap film flowmeter (0101-0113 (1-10-100 mL)) as illustrated in Fig. S4A. The inlet port of the pump was connected to the Tenax TA packed Gerstel tube via a silicone delivery tube. A quartz fiber filter (QFF) was used to seal the sampling end to prevent air borne particles from entering the tube. The values obtained for the two tested pumps (n = 3, for each set rate) were plotted and fitted with exponential decay (order 1) functions as shown in Fig. S4B. By inspection, it can be seen that pump A gave the best calibration results. It is worth noting that, calibration was carried out before and after each sampling session to ensure consistency in the sampled volume. The pump was positioned on a stand, 0.8 m above the floor level in each of the studied indoor environments. Instrumental analysis. Thermal desorption was performed using a commercial desorption unit, TDS-3 (Gerstel) connected to a programmed-temperature vapouriser (PTV) injector/cooled injection system (CIS -3) (Gerstel) by a heated transfer line. Other analysis details are provided in Text S1.

Data Availability
The data associated with this manuscript has been availed as far as possible.