The Use of a Physiologically Based Pharmacokinetic Modelling in a “Full-Chain” Exposure Assessment Framework: A Case Study on Urban and Industrial Pollution in Northern Italy

: Background and goals: The estimate of the internal dose provided by physiologically based pharmacokinetic (PBPK) modelling is a big step forward in the frame of human health risk assessment (HRA) from contaminating sources. The PBPK model included in the MERLIN-Expo platform was here tested with data collected in a human biomonitoring (HBM) pilot study to check model e ﬃ cacy in predicting concentrations in human blood and urine of people exposed to a modern solid waste incinerator (SWI). The aim of the study was to investigate if the use of a PBPK model integrated in a computational platform could replace more expensive and invasive pilot studies. Twenty eight subjects living and working within 4 km of the incinerator (exposed) and 21 subjects living and working outside this area (unexposed) were selected among the population recruited in the HBM study. The group of exposed (E) subjects and the group of non-exposed (NE) subjects were comparable for all relevant anthropometric characteristics and exposure parameters except for the exposure to SWI emissions. Three di ﬀ erent scenarios were created: an “only diet-scenario” (DS), a “worst case scenario” (WCS) and a “most likely scenario” (MLS). The platform was tested for blood-lead (B-Pb), urinary-lead (U-Pb), urinary-anthracene (U-Ant) and urinary-ﬂuoranthene (U-Flt). Average estimated U-Pb was statistically equal to the measured one (est. 0.411~0.278; meas. 0.398~0.455 µ g / L) and estimated vs. measured U-Ant di ﬀ er by one order of magnitude only (est. 0.018~0.010; meas. 0.537~0.444 ng / L) while for U-Flt and B-Pb, the error was respectively of two and four orders of magnitude. It is likely that the extremely high accuracy in the Pb concentration input values referring to diet led to the very accurate estimate for this chemical in urine, but the higher error in the B-Pb computed value suggests that PBPK model equations cannot entirely capture the dynamics for blood compartments. MERLIN-Expo seems a very promising tool in saving time, energy and money in the screening step of the HRA framework; however, many software validations are still required.


Introduction
Finding a precise, reliable and standardized Health Risk Assessment (HRA) framework to investigate the risk of a population exposed to pollutants in the environment is one of the main issues Atmosphere 2020, 11 in Risk Assessment science nowadays [1]. At the same time, the use of a fast method giving a hint of the entity of the risk for the target population is of equal importance for Local Authorities and environmental protection agencies to design and carry out timely, efficient and cost-effective site remediation when the situation is not so critical. To give priority to the most important contaminations, a fast but reliable screening analysis is needed. In the last few years, some new online or offline computational platforms have been developed in this area, integrating toxicological and epidemiological approaches. Among these, MERLIN-Expo software (https://merlin-expo.eu/, [2][3][4]) is a free tool developed within the 4FUN Project (funded by the EU 7th Framework Programme) created for Environmental Exposure Assessment (EEA) and Human Exposure Assessment (HEA). It models complex environmental phenomena and pollution pathways from the sources to their deposition and accumulation in human organs and tissues [5]. Environmental epidemiology has, for many years, adopted a "black box" approach, where human health outcomes were associated to pollutant concentrations in environmental matrixes, considered as proxies of the effective exposure. Thanks to computers' computational capacity increase, Geographic Information Systems (GIS) development and pollutants dispersion modelling reliability, in the last two decades the estimation of the exposure at the individual level was significantly improved. Today, a further step in the exposure assessment is the estimation of the uptake dose of each organ and scientific and technologic improvements of the last few years allowing the study of the toxicokinetics and toxicodynamics of chemicals in the human body using mainly in vitro and in silico approaches rather than in vivo tests. The consequences of this is a change in perspective also in the epidemiologic field in assessing human health risk to chemicals and a growing interest in molecular epidemiology. Computational platforms like MERLIN-Expo allows us to perform a fast "full-chain" exposure assessment analysis that starts from contaminating sources and ends with an estimation of pollutant concentrations in the organs and tissues of the human body. In this HRA framework, PBPK models represent a promising tool aimed at describing the adsorption, distribution, metabolism and excretion (ADME) processes of chemicals in the human body quickly and reliably [6]. The MERLIN-Expo tool was used in this study in the exposure assessment to chemicals (and related risks for human health) related to a specific industrial point source. More specifically, the aim of this study is the comparison of pollutant concentrations measured in samples obtained through a human biomonitoring (HBM) study with the MERLIN-Expo software output, in order to test this tool, and in particular its PBPK model, with a statistic analysis. Testing software like MERLIN and PBPK models [7][8][9] is of extreme importance if Local Authorities want to avoid a huge waste of time, energy and money in HEA and HRA screening steps. Indeed, if these software are sufficiently accurate and precise, they could replace long, invasive and expensive pilot studies in finding essential information to develop a more complex study or to plan public health urgent actions [10,11].

HBM Study
In 2010 a human bio-monitoring (HBM) cross-sectional pilot study was conducted by Modena Local Health Unit in collaboration with the Regional Agency for Prevention, Environment and Energy of Emilia-Romagna to test the use of specific biomarkers in assessing local population exposure to municipal solid waste incinerator (SWI) plant emissions [12]. More specifically, the aim of the study was to investigate the potential role of some specific pollutants as tracers of the exposure to SWI emissions of subjects living near the plant. There were many contaminants, among which there were heavy metals and polycyclic aromatic hydrocarbons (PAHs), investigated which were found in the blood and urine of the 168 subjects enrolled. Inhabitants having residence and workplaces in an area of a 4 km radius centered on the SWI were considered exposed to emissions and those living and working outside that area were considered unexposed. The subjects were chosen with similar job tasks and no significant occupational exposures. A self-administered questionnaire was distributed to investigate diet, life-style, anthropometric characteristics and health status. The ADMS dispersion model was used to produce monthly fall-out maps that showed the pathway of total suspended particles (TSP) emitted from the SWI three emission lines that manage waste for a total of 180,000 tons per year. The exposure time window was assumed as being 30 days before the sampling date, which occurred for all the subjects in the period between May and June 2010. The daily exposure was assumed as 16 h at residence and 8 h at workplace. The associations between the chosen biomarkers and the exposure indexes was investigated through a multivariate linear regression analysis that showed no differences between the exposed and unexposed groups for the majority of the contaminants, among which there were blood-lead (B-Pb), urinary-lead (U-Pb) and urinary-fluoranthene (U-Flt), while few other contaminants, including urinary-anthracene (U-Ant), were found to be higher in exposed subjects.
For the present study, 49 subjects (21 exposed and 28 non exposed) between 36 and 58 years of age have been selected among the 168 people recruited in the HBM study in order to reduce particulate matter (PM) confounding. Two exclusion criteria were adopted: the use of a fireplace in the house to exclude the main source of exposure inside the house, and cotinine levels in serum higher than 3 µg/L, to exclude ex-and current smokers. Information was obtained from the questionnaire answers and the blood analysis data were collected in the HBM study. The platform was tested for B-Pb, U-Pb, U-Ant and U-Flt.
The group of exposed (E) subjects and the group of non-exposed (NE) subjects were chosen to be statistically homogeneous (α = 0.01) for all relevant anthropometric characteristics, i.e., age, height, weight and Body Mass Index (BMI) and exposure parameters (background PM 2.5 , in addition to cigarette smoke [13] and the presence of the fireplace) except for the exposure to SWI emissions.
The Chi-squared test and the analysis of the Skewness and Kurtosis were performed to test for normal distribution of data. The two-tailed homoskedastic Student test was chosen for normally distributed variables, while the K-S test was used for non-Normal distributions. The two groups turned out to be different in the number of males and females only, with 20 out of 28 women in the E group (71%) and 7 out of 21 (33%) in the NE. The results of the statistical analysis are summarized in Tables 1 and 2. A statistical analysis was performed with STATA software (IC 13 StataCorp LP college station, TX, USA) to compare MERLIN vs. HBM concentration distributions. A chi squared test for small samples and a kurtosis and skewness analysis were used to test normality. When the data series were both normally distributed, a t test was applied to compare them, otherwise the non-parametric 2 samples Kolmogorov-Smirnov (K-S) test was used. The K-S test was considered the most suitable test given the 49 values in the distributions while the Wilcoxon or Mann-Whitney tests were used for samples with more than 50 data. Pearson correlation index was used for data normally distributed, while Spearman was chosen for non-normal distributions. Table 1. Chi-square test for normality (χ 2 ) results and a comparison between the E and NE group using t-test (T) and Kolmogorov-Smirnov test (K-S). E = exposed, NE = non exposed. Moreover, two analyses were carried out to investigate a potential correlation between U-Pb values and the distance of the subjects' residences from SWI, using both U-Pb estimated and measured.

Exposure Assessment
Data collected during the exposure assessment process were background particulate matter 2.5 (background PM 2.5 ) and SWI PM 10 concentrations in atmosphere, background PM 2.5 and SWI PM 10 chemical speciation, time spent indoor (home and workplace), diet habits, pollutant average concentration in food in the area of Modena, sex, age, physical parameters and pharmacokinetic parameters.
Subjects' residence and workplace were geo-referenced in the HBM pilot study Using ArcMap 10.1 (REF. ESRI) and all the life-style data were collected by the self-administered questionnaire [12].

Air Route
PM 2.5 daily average concentration inhaled by a subject at his residence was obtained by combining measured and modelled data on particulate matter. PM 2.5 daily average concentration data for the period 2007-2014 were collected from the "Parco Ferrari" urban background station (Modena) and from "via Tagliati" station (Modena), located in the same industrial area of the SWI. A linear regression analysis was used to estimate the missing data in the time series produced by the "Parco Ferrari" station ( Figure 1). The validity of the regression analysis procedure was assured by the small number of the missing data in the time series (<25%) [14]. PM 2.5 annual average concentration at the "Parco Ferrari" station was calculated and the PM 2.5 annual average concentration data at residence were obtained by ARPAE simulations with data from the NINFA-PESCO dispersion model (cell resolution 1 km × 1 km) collected for each subject. The ratio between the two averages was used to estimate PM 2.5 daily average concentration inhaled at residence. PM 2.5 daily average concentration inhaled by a subject in the workplace was obtained following the same procedure. PM 10 monthly average concentrations from SWI to residence and to workplace for the period of interest (April-June 2010) were taken from monthly fall-out maps created by ARPAE with Air Dispersion Modelling System (ADMS) Urban dispersion model [12]. PM 2.5 daily average concentration for the Modena municipality was taken by annual (2010) average concentration in Modena estimated with the NINFA-PESCO model (cell resolution 1 km × 1 km) and was used to assess outdoor exposure. Outdoor activities included in the questionnaire questions were physical activities, car driving and other activities. Given the nature of these activities and of the spaces where they likely took place, the Modena average concentration of PM was considered to be a reliable proxy of the outdoor exposure. The great volume of cleaner air inhaled during a run in a park for few hours per week were in fact balanced by the lower volume of dirtier air inhaled while driving in traffic for a longer time. PM2.5 daily average concentration inhaled by a subject in the workplace was obtained following the same procedure.
PM10 monthly average concentrations from SWI to residence and to workplace for the period of interest (April-June 2010) were taken from monthly fall-out maps created by ARPAE with Air Dispersion Modelling System (ADMS) Urban dispersion model [12].
PM2.5 daily average concentration for the Modena municipality was taken by annual (2010) average concentration in Modena estimated with the NINFA-PESCO model (cell resolution 1 km × 1 km) and was used to assess outdoor exposure. Outdoor activities included in the questionnaire questions were physical activities, car driving and other activities. Given the nature of these activities and of the spaces where they likely took place, the Modena average concentration of PM was considered to be a reliable proxy of the outdoor exposure. The great volume of cleaner air inhaled during a run in a park for few hours per week were in fact balanced by the lower volume of dirtier air inhaled while driving in traffic for a longer time.
PM10 monthly average concentration from the SWI for Modena municipality was taken by PM fall-out maps realized with ADMS simulations. Assuming that a subject spent more outdoor time in proximity of his residence and workplace, for the exposed group the parameter was estimated to be 1 ng/m 3 in April and 1.5 ng/m 3 in May and June, while for the unexposed, it was 0.5 ng/m 3 in April and 1 ng/m 3 in May and June.
Considering the height of the houses where the subjects lived, PM2.5 concentration does not vary significantly with the altitude (less than 5.7% in the cold season and less than 4.1% in the warm season; [15]), and the floor of the house was not taken into account. Vehicle traffic volume was not relevant for this study, being a great source of coarse particles only [16]. In Modena city, the difference in PM2.5 concentrations between heavy and low traffic areas is less than 15% only [15].

Indoor/Outdoor Coefficient
Starting from PM2.5 background concentrations measured by ARPA stations and PM10 concentrations from SWI estimated with ADMS, PM indoor concentrations contributed from outdoors was obtained using an indoor/outdoor coefficient. In fact the coefficient defines how much PM2.5 coming from outside manages to enter indoor spaces. The choice of the right coefficient to assess the real concentration inhaled by subjects in indoor places was an extremely interesting issue [17], as the questionnaire answers showed that the average time spent indoors by the study population was 91% of the total (min. 82%, max. 99%). A two-tailed homoscedastic Student test assured the homogeneity of E and NE for the number of hours spent outdoors per week (p = 0.037; E: 17.3~7.5, NE: 13.0~7.4 h). PM 10 monthly average concentration from the SWI for Modena municipality was taken by PM fall-out maps realized with ADMS simulations. Assuming that a subject spent more outdoor time in proximity of his residence and workplace, for the exposed group the parameter was estimated to be 1 ng/m 3 in April and 1.5 ng/m 3 in May and June, while for the unexposed, it was 0.5 ng/m 3 in April and 1 ng/m 3 in May and June.
Considering the height of the houses where the subjects lived, PM 2.5 concentration does not vary significantly with the altitude (less than 5.7% in the cold season and less than 4.1% in the warm season; [15]), and the floor of the house was not taken into account. Vehicle traffic volume was not relevant for this study, being a great source of coarse particles only [16]. In Modena city, the difference in PM 2.5 concentrations between heavy and low traffic areas is less than 15% only [15].

Indoor/Outdoor Coefficient
Starting from PM 2.5 background concentrations measured by ARPA stations and PM 10 concentrations from SWI estimated with ADMS, PM indoor concentrations contributed from outdoors was obtained using an indoor/outdoor coefficient. In fact the coefficient defines how much PM 2.5 coming from outside manages to enter indoor spaces. The choice of the right coefficient to assess the real concentration inhaled by subjects in indoor places was an extremely interesting issue [17], as the questionnaire answers showed that the average time spent indoors by the study population was 91% of the total (min. 82%, max. 99%). A two-tailed homoscedastic Student test assured the homogeneity of E and NE for the number of hours spent outdoors per week (p = 0.037; E: 17.3~7.5, NE: 13.0~7.4 h).
The physical-chemical transformation processes (deposition, aggregation and volatilization) that the particulate undergoes when it enters buildings changes its chemical composition [18]. During these processes the proportion of some chemicals (e.g., K, nitrates, sulphates) in the PM varies significantly, but this is not true for the substances considered in the present study, e.g., lead (Pb), anthracene (Ant) and fluoranthene (Flt). Due to differences in house ventilation that occurs during the year, a value of 0.91 for the "warm season" (from May to August), of 0.51 for the "cold season" (from November to February) and of 0.71 for the "mild seasons" was chosen [15,19].
A two sample K-S test was used to compare the E and NE groups for the PM concentrations in inhaled air at residence (Figures 2 and 3). At the significance level of 0.01, the null hypothesis (equality of concentrations samples) of the K-S test was refused for residential PM 10 from SWI, while it was accepted for the residential background PM 2.5 . Background PM 2.5 concentration at the workplace was higher for the NE group, while PM 10 from SWI at the workplace was higher for the E group.
A two sample K-S test was used to compare the E and NE groups for the PM concentrations in inhaled air at residence (Figures 2 and 3). At the significance level of 0.01, the null hypothesis (equality of concentrations samples) of the K-S test was refused for residential PM10 from SWI, while it was accepted for the residential background PM2.5. Background PM2.5 concentration at the workplace was higher for the NE group, while PM10 from SWI at the workplace was higher for the E group.  A two sample K-S test was used to compare the E and NE groups for the PM concentrations in inhaled air at residence (Figures 2 and 3). At the significance level of 0.01, the null hypothesis (equality of concentrations samples) of the K-S test was refused for residential PM10 from SWI, while it was accepted for the residential background PM2.5. Background PM2.5 concentration at the workplace was higher for the NE group, while PM10 from SWI at the workplace was higher for the E group.

Chemical Speciation
A background PM 2.5 chemical speciation was conducted by ARPAE Emilia-Romagna in Bologna city. Since the Modena municipality is near Bologna and urban and territorial features, weather characteristics and main pollutant sources are very similar, chemical concentration data generated by this analysis were considered as valid for our study. In 2016, incinerator PM 10 chemical speciation data were collected for the Parma incinerator plant [20]. This new generation plant is very similar to the one in Modena, so those data were used for this study. A chemical composition comparison between background PM 2.5 and incinerator PM 10 is possible since PM from an incinerator plant is mainly composed of ultrafine particles [20]. Incinerator PM 2.5 mass is 94% of the PM 10 and 87% of the total emitted particles (Figure 4), so the operation of adding pollutant incinerator PM 10 and background PM 2.5 into the estimate of the total concentration inhaled produced an error of 6% only.
weather characteristics and main pollutant sources are very similar, chemical concentration data generated by this analysis were considered as valid for our study. In 2016, incinerator PM10 chemical speciation data were collected for the Parma incinerator plant [20]. This new generation plant is very similar to the one in Modena, so those data were used for this study. A chemical composition comparison between background PM2.5 and incinerator PM10 is possible since PM from an incinerator plant is mainly composed of ultrafine particles [20]. Incinerator PM2.5 mass is 94% of the PM10 and 87% of the total emitted particles (Figure 4), so the operation of adding pollutant incinerator PM10 and background PM2.5 into the estimate of the total concentration inhaled produced an error of 6% only. The chemical speciation showed the presence of similar chemicals but in different percentages between incinerator and background PM (  The chemical speciation showed the presence of similar chemicals but in different percentages between incinerator and background PM (Table 3, Figures 5 and 6).     Subjects' exposure to atmospheric contaminants was described adopting the following scheme. Daily average concentration of PM 10 from SWI in inhaled air was multiplied by indoor/outdoor coefficient, by the fraction of contaminant in the PM and by the fraction of time spent in a place to obtain the daily average concentration of contaminants in inhaled air at residence, at the workplace and in outdoor places. The same procedure was repeated for background PM 2.5 . Time spent at the workplace was set at 40 h per week and the time at home was derived by the difference between total time and the sum of the outdoor time and the time at workplace. The time series created were given in input to the software in the "Time series" window. MERLIN "Interpolation-use end values" mode was chosen among the others to pass from the discrete to continuous data distribution. That mode was judged to be the most reliable one to model contaminant concentration changes in the atmosphere because it does not provide very different estimated values in infinitesimal time steps but uses a linear interpolation between input data. The time of exposure considered went from the 31st to the 1st day before the day of the urine and blood samples. A value of 0.75 L was set in the software for the inhaled volume of air during one breath. This is the value for an adult awake in a sitting position and represents a good average of the values assumed by the volume throughout the daily routine. An example of the daily amount of inhaled Ant, Flt and Pb is shown in Figure 7.
atmosphere because it does not provide very different estimated values in infinitesimal time steps but uses a linear interpolation between input data. The time of exposure considered went from the 31st to the 1st day before the day of the urine and blood samples. A value of 0.75 liters was set in the software for the inhaled volume of air during one breath. This is the value for an adult awake in a sitting position and represents a good average of the values assumed by the volume throughout the daily routine. An example of the daily amount of inhaled Ant, Flt and Pb is shown in Figure 7.

Diet Route
Data on weekly intake of portions of food were collected in the HBM questionnaire. The "standard portion" for each type of food was taken from the Italian society of human nutrition (SINU) that defined the quantitative standards for portions in Italy in accordance with consumer expectations [21]. The number of portions of food per week taken for each subject from the questionnaire (summarized in Table 4) was multiplied by the weight of the "standard portion" for each type of food. After some operations, the daily intake of food was estimated and, by multiplying it by the contaminant concentrations in food, the daily amount of Ant, Flt and Pb ingested was found. The "weekly resolution" of these data does not allow us to know the rate of consumption of a certain food at the end of the period of exposure. This can lead to an accuracy loss in the final result. So, some different scenarios were created and some simulations were carried out to test the entity of this error. Since the order of magnitude of the final output did not change for the most reliable scenarios, the entity of the error was eventually considered almost irrelevant. Pb concentration in each type of food and beverage was taken from the IZSLER database and was authorized by the Modena Local Health Unit (USL). Pb time series referred to the period 2011-2017 and food random samples were taken from companies and farms that served markets located in Modena or near municipalities. One concentration value for each sample and for each type of food was collected and an average value was calculated. Concerning data lower than the LOD (1 μg/kg), the value of 0.9

Diet Route
Data on weekly intake of portions of food were collected in the HBM questionnaire. The "standard portion" for each type of food was taken from the Italian society of human nutrition (SINU) that defined the quantitative standards for portions in Italy in accordance with consumer expectations [21]. The number of portions of food per week taken for each subject from the questionnaire (summarized in Table 4) was multiplied by the weight of the "standard portion" for each type of food. After some operations, the daily intake of food was estimated and, by multiplying it by the contaminant concentrations in food, the daily amount of Ant, Flt and Pb ingested was found. The "weekly resolution" of these data does not allow us to know the rate of consumption of a certain food at the end of the period of exposure. This can lead to an accuracy loss in the final result. So, some different scenarios were created and some simulations were carried out to test the entity of this error. Since the order of magnitude of the final output did not change for the most reliable scenarios, the entity of the error was eventually considered almost irrelevant. Pb concentration in each type of food and beverage was taken from the IZSLER database and was authorized by the Modena Local Health Unit (USL). Pb time series referred to the period 2011-2017 and food random samples were taken from companies and farms that served markets located in Modena or near municipalities. One concentration value for each sample and for each type of food was collected and an average value was calculated. Concerning data lower than the LOD (1 µg/kg), the value of 0.9 µg/kg was set to produce an overestimation and to balance the lower estimation due to the missing foods and beverages in the questionnaire questions. The Pb average concentration in meat had taken into account 76 pieces of data collected between 2011 and 2012 that referred to cow, chicken, pig and rabbit muscles. Moreover, the data for milk, eggs and wine were collected in that period while for the remaining foods, the samples were taken between 2013 and 2017. Concentration average values for Flt were taken from ISS reports [22] and are estimates at the national level derived from a European and extra-European literature analysis [23][24][25][26][27][28]. Since the reliability and accuracy of Flt data were significantly lower than those of the Pb data, two different scenarios were created for each exposed subject: the worst case scenario (WCS) and the most likely scenario (MLS) ( Table 5). In the first scenario, following the precautionary principle, the highest values found in Italy were chosen (to balance the underestimation due to missing foods and beverages in the questionnaire questions) and in the second, the average values were chosen. More precisely, in the WCS, in case a class of food included more maximum data, one for each type of food in that class and a maximum value for the class was produced by computing the average of all the maximum values. Data referred to food from highly contaminated areas were not taken into account. Only the WCS was chosen for the unexposed subjects. A ratio between Ant and Flt average concentrations in food from the European Union (ISS data) was computed to adopt a likable value for Ant concentration in food. The value of this ratio was 0.34. The Ant concentration for the food in Modena was obtained by multiplying this value for Flt concentrations. Chemical concentrations found in each type of food was multiplied by MERLIN for the daily amount of food ingested (data from questionnaire) to find the daily average intake of chemicals due to that type of food. The MERLIN "exposure" module computed the total intake of chemicals and gave it as an input to the PBPK model. An "only diet scenario" (DS) was created by considering the contribution of the diet route only (diet data assumed like in the WCS) and then, by comparing the WCS and the DS, the atmospheric contribution was isolated.

Use of MERLIN-Expo
The MERLIN-Expo tool is composed of many environmental and transfer models that follow the contaminants of interest from their sources to the human/population target and a pharmacokinetic-based model (PBPK model) that computes the specific tissue concentration for each element and compound of the mixture of substances entered in the organism via ingestion or inhalation. The models have different structures and tasks and can be linked to create a specific pathway, allowing us to model chemical transfers from pollutant sources to environmental matrixes, to animals and humans. The software was used in this study mainly to gather all the exposure information and estimate B-Pb, U-Pb, U-Ant and U-Flt. Diet info and atmospheric exposure data were given as inputs to the MERLIN PBPK model to obtain chemical doses at organs. The software does not take into account dermal route, but only ingestion and inhalation. Since the people analyzed in this study were all adults working in offices, the dermal route seemed to be negligible for the pollutants considered. In order to achieve this goal, in the MERLIN "Model" window a conceptual model was created linking the sub-models "Human Intake", where all the exposure data were collected together, and "Man", was in the PBPK model. The contaminant concentration in air and food and assumption rates were given as input in the "Human Intake" sub-model and "Total concentration inhaled" and "total quantity Ingested" were the requested outputs. These outputs were given as inputs to the "Man" sub-model, which provided the contaminant concentration in tissues and the "Transfer of contaminants from kidneys to urine" as outputs.
To consider the three different scenarios created (WCS, MLS and DS), 119 simulations were performed with the MERLIN-Expo tool. The simulation output was the concentration of Pb, Ant and Flt in the 21 compartments outlined in the software and the "transfer of contaminants from kidneys to urine", that is the mass flow of pollutants transferred to urine daily. The chemical concentrations in urine was calculated by dividing this last output by the volume of urine per day. The daily flow rate of urine was obtained from the ratio between daily mass of creatinine excreted from the human body (men: 1.642 g/day, women: 1.041 g/day; [29]) and creatinine concentration in urine. The use of generic values from the literature for men and women for the daily mass of creatinine excreted seemed an effective way to obtain a reliable estimate of the urine flow rate when measured values were not available, even if the daily amount of creatinine produced usually differed consistently from one person to another. The estimated daily average urine volume and the standard deviations turned out to be very similar for the two groups, but the standard deviations were quite wide due to differences in creatinine mass excreted among subjects of different sexes and ages (E: 1.17~0.48 L/day; NE: 1.22~0.84 L/day).

PBPK Model Parametrization
In MERLIN "Parameters" window chemical-physical and pharmacokinetic properties of pollutants are investigated. The value of a parameter is required in input and is fixed during each simulation but it can be modified from one simulation to another without varying other parameters to perform sensitivity analysis and calibrate parameters themselves. A database with the values of principal pharmacokinetic parameters for some chemicals is available in MERLIN-Expo software. Since the list of parameterized substances includes Pb but not Ant and Flt, the missing values for the pharmacokinetic parameters required by the MERLIN PBPK model were taken from the INTEGRA platform (http://www.integra-lri.eu/) [30].
The Michaelis-Menten Equation is used in MERLIN PBPK model to describe chemical metabolism [31] and the Michaelis-Menten constant is an equation parameter that depends on pollutant concentration and solubility [32,33]. In the INTEGRA database, the value of the constant for Ant is equal to 42 (µmol/L). This value was converted into 7.49 (mg/L)using Ant molecular weight and was given as input in MERLIN with this unit of measurement. The same procedure was used for Flt, obtaining a value of 8.50 (mg/L). The value of this constant for Pb was kept at 1 (mg/L), as fixed in MERLIN. The maximum velocity per kg of bodyweight is another parameter of the Michaelis-Menten Equation that is required in input. The Ant INTEGRA software provides a value of 480 (µmol/h) for this parameter. The value was converted to 1.43 (mg/min) using the molecular weight of the chemical to give it as input in MERLIN and individual values were derived using body weight data from the questionnaire. The same procedure was adopted for Flt and a value of 1.62 (mg/min) was obtained while the Pb value of the parameter proposed by MERLIN was not changed. The "Main cytochrome P450 involved in metabolism" parameter was selected as "others" for all three contaminants, since the main cytochrome for Ant, Flt and Pb is not selectable among the MERLIN options. The value proposed by MERLIN for the "biliary excretion rate" parameter was not modified for the three pollutants. The excretion rate in kidneys was assumed to be 0.248 (1/h) for Flt [29] and Ant. This value was then divided for the body weight of the subject to find individual values for the parameter "excretion rate per kg of bodyweight". The values of this parameter for kidneys and liver referred to Pb are given by MERLIN and they were not modified. The "blood:air partition coefficient" parameter is set as equal to 10 99 for non-volatile compounds and Pb is part of this class of substances. Since INTEGRA suggests a value of 3,563,640 for PAHs, this was adopted for Ant and Flt. The "tissue-blood partition coefficients" for Ant and Flt were provided by INTEGRA, while for Pb the values adopted by MERLIN were not changed. The "saturation constant for partitioning between blood and plasma" was not relevant for Ant and Flt, since their concentration in blood was surely too low to saturate the plasma, and so the default value indicated by MERLIN (1 mg/L) was not changed. The "Ingestion via the liver" mode was selected and the parameter 'fraction of contaminant absorbed via ingestion' was set equal to 0.9 according to INTEGRA developers, for the two PAHs and equal to 0.11 for Pb, as chosen in MERLIN. Moreover, for the parameter "fraction of contaminant absorbed via inhalation", INTEGRA proposes the value of 0.9, while MERLIN assigns a value of 0.5 to this fraction for Pb. The MERLIN database provided the values for the "relative blood flows in adults", the "density of organs" and the "relative organ weight in adults" parameters [34].

Comparison between the MERLIN Platform and the HBM Study Results
The average estimated U-Pb was statistically equal to the measured one (est. 0.411~0.278; meas. 0.398~0.455 µg/L) and the average estimated U-Ant differed from the measured U-Ant from one order of magnitude only (est. 0.018~0.010; meas. 0.537~0.444 ng/L), while for U-Flt and S-Pb, the error was respectively of two and four orders of magnitude ( Table 6). The chi-squared test showed low p-values (<0.01) for all the distributions and the null hypothesis of the K-S test was not rejected for U-Pb only (p-value: 0.106) at significance level of α = 0.01. A correlation analysis was performed between urinary lead (U-Pb) concentration values estimated by MERLIN software and those measured reported in the HBM study ( Figure 8 MERLIN software slightly overestimated U-Pb measured concentrations for values less than 0.25 µ g/L and slightly underestimated it for values over 1 µ g/l, while maximum precision and accuracy was in the range of 0.25-0.75 µ g/L (Figure 9).   MERLIN software slightly overestimated U-Pb measured concentrations for values less than 0.25 µ g/L and slightly underestimated it for values over 1 µ g/l, while maximum precision and accuracy was in the range of 0.25-0.75 µ g/L (Figure 9).  The DS showed that diet was the main route of exposure for all the contaminants investigated (Pb: 98%, Ant: 95%, Flt: 89%), highlighting that the contribution of air pollution contribution to uptake (background and from SWI) was almost negligible. This result is perfectly in line with the scientific literature on this topic.

Exposed vs. Non-Exposed
Furthermore, the comparison between MERLIN and HBM data was analyzed by taking into account the E group and the NE group separately ( Table 7). The analysis of the results for the E and NE groups considered separately was very similar to the previous analysis for all the subjects. U-Pb concentration was correctly estimated, while U-Ant, U-Flt and B-Pb estimated averages were considerably lower than the measured ones, as shown in Tables 7 and 8. A K-S test showed no differences between the E and NE groups for all the contaminant average concentrations estimated with MERLIN at a significance level of 0.05 (p-values: U-Pb: 0.896, B-Pb: 0.134, U-Ant: 0.340, U-Flt: 0.400), while in the HBM study U-Ant was higher for the E group (Tables 7 and 8).

Concentration in Urine
The shape of the urinary concentration curve of the chemical varies consistently from one pollutant to another but not so much between subjects. After one week from the beginning of the exposure, the U-Pb transferred from kidneys to urine-directly related to its concentration in urine-is about 79% of the value found after one month and the value after 24 days is 99% of the value after 30 days for all the subjects recruited in this study ( Figure 10). Conversely, in the first days after the exposure, the concentration increase is fast, showing a near-logarithmic trend of the curve. U-Ant and U-Flt concentration curves are very similar: after one day the concentration is nearly the same as after one month, even if from day 2 to day 30 values fluctuate around an average value ( Figure 11). For all the pollutants, the concentration value after one month is almost the same as the value after one year. Those trends support the choice to start the exposure one month before the sampling, in particular for U-Pb.
Atmosphere 2020, 11, x FOR PEER REVIEW 15 of 22 For all the pollutants, the concentration value after one month is almost the same as the value after one year. Those trends support the choice to start the exposure one month before the sampling, in particular for U-Pb.

Concentration in Organs and Tissues
The MERLIN software was tested here only for blood and urine, but the next studies should focus attention also on the concentration estimates that the software provides for the 22 organs and tissues included in the PBPK model. The tool showed for all the subjects that lung tissue is the compartment with the highest Pb values (e.g., subject 1: 5.5 × 10 −5 mg/L, after an exposure time of one month). The second compartment most subject to Pb accumulation was bone (one order of magnitude less than lung), the third blood, then brain, urinary tract and kidneys ( Figure 12). PAH's estimated concentrations in tissues are much lower than the Pb ones ( Figure 13). After a period of exposure equal to one month, for both Flt and Ant, the compartment with the highest values was the liver (Flt: 7.83 × 10 −7 mg/L), followed by blood (Flt: 3.62 × 10 −9 mg/L), lungs (Flt: 2.49 × 10 −9 mg/L),

Concentration in Organs and Tissues
The MERLIN software was tested here only for blood and urine, but the next studies should focus attention also on the concentration estimates that the software provides for the 22 organs and tissues included in the PBPK model. The tool showed for all the subjects that lung tissue is the compartment with the highest Pb values (e.g., subject 1: 5.5 × 10 −5 mg/L, after an exposure time of one month). The second compartment most subject to Pb accumulation was bone (one order of magnitude less than lung), the third blood, then brain, urinary tract and kidneys ( Figure 12). PAH's estimated concentrations in tissues are much lower than the Pb ones ( Figure 13). After a period of exposure equal to one month, for both Flt and Ant, the compartment with the highest values was the liver (Flt: 7.83 × 10 −7 mg/L), followed by blood (Flt: 3.62 × 10 −9 mg/L), lungs (Flt: 2.49 × 10 −9 mg/L), muscles (Flt: 3.91 × 10 −12 mg/L), gut (Flt: 3.27 × 10 −12 mg/L) and kidneys (Flt: 3.01 × 10 −12 mg/L). Especially for Pb, the shape of the curve was found to be quite similar between the different organs for the same pollutant. However, the slope of the curves is different, so the steady state is reached first in some compartments (e.g., adipose tissue and kidneys) and then in others (e.g. bones). Especially for Pb, the shape of the curve was found to be quite similar between the different organs for the same pollutant. However, the slope of the curves is different, so the steady state is reached first in some compartments (e.g., adipose tissue and kidneys) and then in others (e.g. bones).

Distance Analysis
No correlation (Pearson coefficient = −0.06) was found between the U-Pb estimated concentrations and the distance of the residence from the SWI plant ( Figure 14). A similar result was found by carrying out the same analysis using HBM data (Pearson coeff. = −0.04).
Atmosphere 2020, 11, x FOR PEER REVIEW 17 of 22 No correlation (Pearson coefficient = −0.06) was found between the U-Pb estimated concentrations and the distance of the residence from the SWI plant ( Figure 14). A similar result was found by carrying out the same analysis using HBM data (Pearson coeff. = −0.04). No correlation was found between chemicals concentration values in urine and blood and the presence of a boiler inside the house.

Discussion
We tested a PBPK model inside a computational platform for the human "full-chain" exposure assessment of chemicals. We considered three environmental contaminants (Pb, Ant, Flt) measured in the blood (Pb only) and urine of people exposed to urban and industrial pollution, and we compared the measured concentrations with MERLIN output. We also focused attention on the No correlation was found between chemicals concentration values in urine and blood and the presence of a boiler inside the house.

Discussion
We tested a PBPK model inside a computational platform for the human "full-chain" exposure assessment of chemicals. We considered three environmental contaminants (Pb, Ant, Flt) measured in the blood (Pb only) and urine of people exposed to urban and industrial pollution, and we compared the measured concentrations with MERLIN output. We also focused attention on the contribution of the near SWI to the total exposure.
The results found in this and other studies [20,35] highlight the necessity to continue the monitoring of all incinerator plants, even if they are modern and new. In fact, even if they emit a very low amount of particles, they are richer in pollutants than those forming the background PM. In 2010, the background PM 2.5 annual average concentration at subjects' residences was 20.98 (~1.44) µg/m 3 for exposed subjects and 20.92 (~1.11) µg/m 3 for unexposed subjects, while during the period of April-June 2010, the incinerator PM 10 average concentrations were 1.30 (~1.16) ng/m 3 and 0.67 (~0.23) ng/m 3 , respectively. Taking into account the results of the chemical speciation, the amount of chemicals inhaled by subjects from an incinerator is two (Pb) or three (Ant and Flt) orders of magnitude less than that inhaled from all other sources.
According to the analysis described in this study, MERLIN-Expo software produced encouraging results, but several open issues remain to be solved.

Comparison between MERLIN Software and HBM Study Results
Likely, the reliability and the accuracy in the Pb concentration input data led to a very accurate estimate for this chemical in urine. In particular, the concentration values in food and beverages taken from the IZSLER database led to a very precise description of the amount of Pb ingested with the diet, which seemed to be the main route of exposure for all the pollutants studied. The major error in estimated B-Pb suggests that some mistakes may exist in PBPK model equations and/or in the pharmacokinetic parameters for the blood compartment. This is the most likely conclusion, since the procedure, the input data and the tools used were the same for the U-Pb estimation.
The need to take values for some PBPK model parameters from other software probably affected U-Ant and U-Flt outputs; however the values of concentration in food taken from the literature might be considered responsible for the errors in estimated concentrations. However, the difference between MERLIN outputs and HBM study results for Flt and Ant were probably not due to poor data accuracy, since the adoption of the WCS should have shown the opposite result (a higher MERLIN average value than HBM one). A comparison between MERLIN and INTEGRA pharmacokinetic parameters for Pb would have been interesting in order to understand the entity of error for Flt and Ant. However no data were available in the INTEGRA database for Pb or other heavy metals with similar physical-chemical properties. Anyway, a significant difference between the two pieces of software in parameter values for benzo(a)pyrene suggests that the hypothesis of a different parameterization for PAH should really be taken into account. U-Ant and U-Flt concentration values in WCS and MLS were of the same order of magnitude. This fact confirms the hypothesis of an error in the MERLIN-Expo software in modelling the ADME of these pollutants or in INTEGRA pharmacokinetic parameters and assures that the final result is not affected by the inaccurate input data in the diet route.
The medium/low Spearman correlation coefficients between MERLIN and HBM ordered values for all the chemicals mainly exist for two reasons. First of all, the adoption of only two standard (and not individual) values (one for men and one for women) for daily creatinine production by subjects led to an error in urine daily volume estimation. Secondly, the assumption of a constant daily intake of chemicals through diet route did not allow us to gain an accurate estimation of pollutant concentrations in short periods of time, since from the questionnaire it was impossible to know the amount of hours passed from the intake of that specific food to the blood and urine sampling. These facts led to the lowering of the correlation coefficient for all the pollutants but, likely, did not affect the population average values consistently.

Exposed vs. Non-Exposed
Since K-S tests showed a significantly higher PM concentration from the SWI inhaled by the exposed subjects and a similar exposure to background PM for the E and the NE groups, the differences in chemicals concentrations between both groups were due to SWI contributions only.
By using the MERLIN-Expo tool, the U-Ant was not different between the exposed and non-exposed groups while, by using data from the HBM study, a difference between the two existed. As a consequence, the exclusive use of the software in the exposure assessment would have led to an error in considering U-Ant concentration equal for both groups, while U-Ant was found higher in the E group. As concerns B-Pb, U-Pb and U-Flt, final considerations of the comparison between the E and the NE groups shown with MERLIN were found to be the same in an analysis of the two groups using measured data.
The comparison between WCS and DS showed that the contribution of atmospheric exposure was higher in the E group than in the NE group only to U-Flt estimated concentrations, but when the contribution was referred to the measured concentrations, the difference between the two groups disappeared.

Diet Scenario
Since Flt and Ant concentrations in foods consumed in Modena were surely less than those adopted in the WCS (maximum concentrations found in Italian food), the error in the estimated concentrations in urine (lower than measured concentrations) is mainly attributable to the tool and not to the low accuracy of the input data.
Without the DS, even if MERLIN had shown different values between E and NE, the increase in pollutant concentrations in urine and blood of the E group would have not been clearly and exclusively associated with the incinerator, because the little differences in diet between the two groups could have played an important role.
Moreover, concerning the E group, the U-Ant and U-Flt estimated average concentrations in the WCS and in the MLS were of the same order of magnitude (Ant WCS: 0.016~0.008 (ng/L), Ant MLS: 0.004~0.002 (ng/L), Flt WCS: 0.051~0.026 (ng/L), Flt MLS: 0.015~0.007 (ng/L)), despite the large difference in PAH concentration in food. Once again, this fact supports the thesis that the error in concentration estimates in urine is not due mainly to inaccurate data in the diet route.
Since more accurate and precise data were used to describe the contributions of the air route than the concentrations used in the diet route, a comparison between the values of the two groups was carried out by considering the difference between the WCS and the DS to analyze the contribution of the SWI more precisely. The atmosphere contribution to the exposure was isolated (Table 9) thanks to this procedure and, the two populations being homogeneous for the exposure to background PM 2.5 , the differences between the two groups could be attributed to the SWI only. A two sample K-S test (Table 10) showed that the average concentration was not significantly different between the E group and the NE group (null hypothesis never rejected). The percentage of pollutant concentrations in urine and blood due to atmospheric contribution only on the total value estimated with MERLIN (Table 11, M) and on the total value measured in the HBM study (Table 11, HBM) showed a very low contribution of the air route to the total exposure. The estimates on the measured values are not affected by errors in the diet route assessment but by those on the air route only, and these are low, but the differences between measured and estimated total concentrations are too large to give credibility to these values, except for U-Pb. According to MERLIN simulations, for both groups the atmospheric contribution to the total exposure seems very small compared with the diet one (Tables 12 and 13). Table 9. Contaminant concentrations due to the air route only. E = exposed, NE = non exposed.   Table 11. Ratio (%) between the contaminant concentration in blood or urine due to the air route only and total concentration estimated by MERLIN (M) or total concentration measured in the Bio-monitoring study (HBM). E = exposed, NE = non exposed. A.V. = average value. Pb = Lead, Ant = Anthracene, Flt = Fluoranthene.

Statistical indices
Pb   Table 12. Ratio (%) between the estimated contaminant concentration in blood or urine due to the diet route only and total concentration estimated by MERLIN. E = exposed, NE = non exposed.  Table 13. Ratio (%) between the estimated contaminant concentration in blood or urine due to the diet route only and total concentration measured in the HBM study. E = exposed, NE = non exposed.

Mass Balance
An equation that aims to estimate the contaminant mass losses at a fixed time was included in MERLIN to control the mass balance of the system of differential equations. The smaller the mass losses are, the lower the error in predicting concentration in organs is, and the higher the precision of the results is. For all the subjects, the order of magnitude of this mass loss was equal to 10 −6 mg. Since this error is distributed over all the compartments, its order of magnitude becomes 10 −7 mg when it is referred to the estimated concentration in urine. As a consequence, the error is negligible when it is referred to the estimated mass of Pb transferred from kidneys to urine daily (10 −4 mg/day), but it is not so when it is compared to PAH's values (10 −8 mg/day) and Pb concentrations in blood (10 −7 mg/day). This confirms the high precision of the software in Pb estimate and more doubts are put on equations for PAHs and B-Pb. By using a sub-group of the wider sampled population chosen in the HBM study, the significant difference in U-Flt concentrations between the two groups disappears in this study.

Conclusions
The aim of this study was to test the PBPK model in MERLIN-Expo software through a comparison between software outputs and data collected as part of the HBM study. Following this purpose, B-Pb, U-Pb, U-Ant and U-Flt concentration data collected as part of the HBM study were compared with estimated concentration values shown in the software outputs through a statistic analysis. The extreme accuracy of MERLIN software in estimating U-Pb concentrations encourages the use of PBPK models and HEA tools in HRA screening studies. However, the low accuracy in B-Pb, U-Ant and U-Flt estimated concentrations, even if tested on 49 subjects only, highlights the need to improve the software. As a consequence, first of all, the provision of a value for the missing parameters should increase output accuracy and reduce time. Secondly, the corrections to the equations for HPAs and blood compartments seem to be necessary. Moreover, testing the other compartments of the MERLIN-Expo PBPK model with measured concentrations from biopsies should be seen as a very important issue in order to validate the model completely. In conclusion, MERLIN-Expo seems a very promising tool in saving time, energy and money in the screening step of HRA frameworks. However, a more comprehensive evaluation of other similar tools is recommended to provide important information on how to obtain a sufficiently accurate, fast and cheap tool to plan public health protection interventions.