宋韬，黄青,*，徐琪依，柴源

1. 暨南大学环境学院， 广州510632 2. 广东省环境污染与健康重点实验室， 广州 510632

有机污染物暴露水平与健康风险问题一直是环境与健康领域的主要研究方向之一，有重要的现实意义和研究价值。多环芳烃(polycyclic aromatic hydrocarbons，PAHs)作为环境中最常见的一类有机污染物，其来源很大一部分为人为源，主要体现为工业源、交通源和生活源等[1]。其中，生活源由于浓度低、频率不固定等特点，其暴露风险容易被人们忽视，烧烤源就是其中一个典型代表。随着人们生活方式的改变，更多的人们喜欢烧烤类的食物，甚至形成一种独特的饮食方式。在烧烤活动中，生物质燃料的不完全燃烧或肉类食品脂肪的高温热解会产生大量PAHs类化合物，对近距离接触烧烤源的厨师、服务员等人群以及食客产生较大的健康影响[2-3]。通常，PAHs会通过呼吸作用、皮肤吸附、膳食暴露等途径进入人体[4-6]，在细胞色素酶的代谢作用下PAHs在体内降解，并通过排泄物排出体外[5-8]。

生理药代动力学(physiologically-based pharmacokinetic，PBPK)模型是广泛适用于各类暴露场景[9-14]的环境模型之一，其特点在于可以定量预测目标化合物在暴露对象体内的吸附、分布、代谢动力学过程，并能够根据实际或假设的环境条件进行暴露模拟。PBPK模型已经被研究人员大量地应用于多种化合物暴露剂量的模拟，如苯并[a]芘 (B[a]P)[9-11]、溴代二氯甲烷[12]、双酚A[13]、PAEs[14]、镉[15]、DDTs[16]等。我们尝试使用基于中国人群参数建立的PBPK模型来解决本地化烧烤场景中国居民PAHs的暴露剂量相关问题。

目前，已经有大量文献报道不同人群尿液中羟基多环芳烃(OH-PAHs)浓度水平[17-19]。然而，人体组织内暴露标志物浓度的浓度确因样品难以获得与实验困难等原因，数据则相对较少。本文利用蒙特卡洛法和本地化的PBPK模型分别模拟烧烤源暴露人群PAHs外暴露剂量与各人群体内组织中PAHs暴露标志物浓度变化过程，计算了人群终身致癌风险和暴露标志物富集水平。探讨了暴露途经、性别等因素对暴露剂量、标志物浓度的影响，为环境健康模型模拟提供可参考案例。

1 材料与方法(Materials and methods)

1.1 实验方案

为了分析不同人群暴露于烧烤源后的PAHs外暴露水平与模拟人群各组织中PAHs暴露标志物浓度，本研究建立了如下相应的蒙特卡洛模型与PBPK模型。年龄、性别、职业差异是环境暴露水平与人体健康评估中3种重要的影响因素。本研究中，我们选取了4组不同人群作为暴露对象，分别有：成年男性，成年女性，儿童(6～12岁，不区分性别)，职业人群。职业人群是指近距离接触烧烤源，日均暴露时间和频率均高于普通居民的厨师、服务员等相关从业人群。为了简化模型，假定职业人群的暴露参数除呼吸速率(IRinh)外与成年男性相同。选择2013年一次野外采样[20]的烧烤源环境数据(距离烤炉2 m 和10 m 处)作为PAHs暴露条件，进行内、外暴露剂量评估。

外暴露剂量使用B[a]P等效浓度(B[a]P-eq)计算得到。暴露条件根据Wu 等[20]的研究确定为：气相B[a]P-eq(3.7±5.5) ng·m-3、颗粒相(10±8.7) ng·m-3(2 m处)和气相(2.5±3.8) ng·m-3、颗粒相(2.5±0.5) ng·m-3(10 m处)。而对于内暴露，则采用PAHs的一种单体——芘在人体内的“组织-时间序列浓度”(tissues time-course concentration)模拟描述，暴露条件根据Wu等[20]和Li等[21]的研究确定为空气环境中芘浓度为(155±5) ng·m-3，烧烤类食物中芘含量平均值为99 μg·kg-1。把实验对象中的成年男性、成年女性、儿童作为非职业组，暴露途经有呼吸、皮肤接触和膳食暴露，假定其平均进食时间作为暴露时长，约为1 h·d-1。而职业人群则通过呼吸、皮肤接触暴露于PAHs，并选择法定每日工作时间上限8 h(8 h·d-1)作为日均暴露时长。

1.2 外暴露模拟与健康风险

日均摄入量(estimated daily intake, EDI)根据暴露途径分为呼吸(inh)、皮肤接触(der)和膳食摄入(diet)，分别根据式(1～4)在Matlab平台上进行蒙特卡洛模拟运算，循环次数为1 000次。每组1 000个独立的模拟结果再通过正态拟合函数normfit得出α=0.05条件下的EDI均值与标准差。

EDIinh(ng·d-1) = (Cg+Cp)· IRinh· ET

(1)

EDIder(ng·d-1) =Cg·kp· SA ·fsa· ET

(2)

EDIdiet(ng·d-1) =Cfood· IRdiet

(3)

EDI = EDIinh+ EDIder+ EDIdiet

(4)

其中，Cg,Cp(ng·m-3)分别表示气相、颗粒相B[a]P-eq浓度；IRinh(m3·h-1) 表示呼吸速率，对于男、女、儿童、职业人群取值分别为U(6.7, 11.1)，U(6.3, 8.7)，U(7.0, 9.2)，U(17.2, 28.5)(U为均匀分布)；kp(m·h-1)，B[a]P的皮肤渗透系数，取值为2.5 m·h-1[6]；SA (m2)是皮肤表面积；fsa是皮肤暴露比例，取值25%[20]；ET (h·d-1)表示暴露持续时长，对非职业人群和职业人群分别为1 h·d-1和8 h·d-1；Cfood(ng·g-1)是食物中的B[a]P-eq浓度，数值取0.4～0.6 ng·g-1 [13]。IRdiet(g·d-1) 代表日均食物摄入量，取值分别为成人为98.4 g·d-1, 儿童为U(68.1, 88.9) g·d-1。

人群终身致癌风险按照式(5)计算：

ILCR =∑CSFi· EDI· ED· EF · cf · bw-1· AT-1

(5)

其中，CSF，经呼吸、皮肤接触、膳食暴露途经的致癌斜率因子分别为3.14、37.4、7.3 kg·d·mg-1(US EPA, 2013)；ED，终身暴露年数，成人取52年，儿童取6年；EF，暴露频率，假定职业人群每周工作5 d，可得到260 d·a-1；cf，转化系数，10-6；AT，平均时长，25 550 d。

1.3 内暴露模拟(PBPK模型建立)

芘作为美国环保署(US EPA)推荐的16种优控PAHs之一，芘和它的主要代谢物(1-羟基芘)无论是对非职业或职业人群都被认为是非常有效的PAHs暴露标志物[17-18]。因此，本研究选择模拟人暴露于芘后的体内血液循环过程、肺部气血交换过程、皮肤吸附、肝代谢等生化过程。基于先前的研究[9,11,13-14]，本研究根据物质平衡原理建立了如图 1中所示的PBPK模型结构。房室具体包括肺、脂肪、皮肤、肾、肝、血液，均认为是灌注限制(perfused-limited) 组织。

1.3.1 模型参数

本研究中，所用的PBPK模型生理参数包括不同人群的组织体积(Vi)、组织血流量(Qi)、体重(bw)、皮肤表面积(SA)、呼吸速率(IRinh)(表2)。通常情况下，模型中组织体积与组织血流量可以分别通过对象的体重和心输出量(cardiac output, CO)计算得出(式6，7)[22-23]。式(6，7)中，i代表脂肪、皮肤等组织，比率系数(Vic,Qic)在表1中给出。心输出量同样也可通过暴露对象体重计算得出[13](式7)。同时，体重、皮肤表面积、呼吸速率等暴露参数也是模型生理参数的一部分(表2)，取值参考自《中国人群暴露参数手册》。模型中的物化参数只与化合物本身性质有关，用于表示化合物经皮接触、经口摄入以及组织-血液、空气-血液分配系数，不同的污染物有各自的分配系数和皮肤渗透系数。其中皮肤渗透系数(kp)，气-血分配系数参考Heredia-Ortiz等的结果[9]。组织-血液分配系数参考Brown等[23]的结果。

经口摄入系数(oral absorption coefficient,ko)用于表示食入剂量被人体(肝组织)吸收的速率大小，它直接决定了膳食暴露途经中，组织浓度的动力学过程。本研究中，综合了Haddad 等[25]、Crowell等[11]和Viau等[4]3篇早期关于大鼠/小鼠PAHs经口芘暴露研究的结果，取他们通过实验数据得到的优化ko值(分别为0.035 min-1、0.005～0.01 min-1、0.014 min-1)的平均值(0.023 min-1)作为本研究的ko参数(表3)。吸附比例(absorption fraction，fabs)定义为进入代谢组织的芘最终以代谢物形式通过尿液代谢出的比例[4]。在本研究中，我们综合了Viau等[4]、Chien等[5]和Li等[21]关于芘膳食暴露研究中使用或得到的优化参数fabs，取它们的平均值3.76%为本研究使用参数(表3)。米-门(Michaelis-Menten)机理的一阶线性估计描述肝组织代谢过程是最常用的代谢机理之一(式10)。其关键代谢参数——消除速率常数(intrinsic clearance, CLint)等于最大代谢速率常数和米氏常数的比值(Vmax/Km)[9]。本文采用了Haddad等[26]通过体外模拟实验得到的用于人肝代谢芘的CLint(0.068±0.025 L·min-1·kg-0.74)作为代谢参数(表3)。

图1 PBPK模型结构Fig. 1 PBPK model structure

表1 生理参数表：血流量比例，组织体积比例与组织-血分配系数Table 1 Physiological parameters: blood flow fraction, volume fraction and tissues-blood partition coefficients

表2 中国人群暴露参数Table 2 Chinese population exposure parameters

注： “—”表示本研究中PBPK模型与外暴露量估算不需要用到的参数。Note: “—” means unused parameters in PBPK and external dose modeling.

表3 芘的物化参数和代谢参数Table 3 Physico-chemical and metabolic parameter for pyrene

注：pba=2.04，kp=0.0119 cm·h-1原本在Heredia-Ortiz等[9]的研究中用于B[a]P模型。Note:pba=2.04,kp=0.0119 cm·h-1have been used in B[a]P model in the reference[9].

Vi(L) =Vic· bw

(6)

Qi(L·min-1) =Qic· CO

(7)

CO (L·min-1) = 15.87 · bw0.75

(8)

1.3.2 PBPK模型运行

根据物质守恒原理，建立图1所示的PBPK模型，各组织浓度用微分方程组(式8～12)表示。利用Matlab对建立的PBPK微分方程组进行数值求解，以分钟作为时间步长，初始条件设定为各组织浓度均为零。

脂肪和肾：i=f或k

Vi· dCi/ dt=Qi· (CA -Ci/pci)

(9)

皮肤：

Vs· dCs/ dt=

Qs· (CA -Cs/pcs) +kps· SA ·fsa·Cg

(10)

理论上，肺、皮肤生物转化，肝代谢(hepatic metabolism)作用均会使体内的PAHs含量减少。然而，Heredia-Ortiz等[9]认为因肝代谢是模型中主要PAHs代谢过程，肺、皮肤生物转化过程相对前者可以忽略。又因肝组织是食物经胃肠通道进入人体的目标组织，PAHs在肝组织的浓度可以表示为:

Vl· dCl/ dt=Ql· (CA-Cl/pcl) +fabs·ko·Ad

-CLint·Cl/pcl

(11)

其中，Ad为目标化合物食入量(ng)。

血液：

CO · CA = ΣQi·Ci

(12)

气血交换过程：

CA = (CO · CV + IR ·Cg) / (CO + IR /pba)

(13)

2 结果与讨论(Results and discussion)

2.1 外暴露模拟

2.1.1 B[a]P-eq日均摄入量

根据B[a]P-eq浓度计算得出的4类人群(男、女、儿童、职业人群)EDI分别为(54±5.9)，(54±5.8)，(43±5.4)，(178±98) ng·d-1(图2)。对于非职业人群，EDIinh和EDIder都表现出“成年男性>成年女性>儿童”的规律。并且由于EDIdiet的人群差异主要由日均进食量造成，《手册》中的儿童日均肉类消耗推荐值小于成人推荐值，因此，儿童EDIdiet也小于成年人。职业人群EDIinh和EDIder均显著大于非职业人群。

对于非职业人群，EDI贡献比例的分析显示：第一，膳食暴露对总EDI的贡献度均大于90%，是绝对主导因素。第二，皮肤接触产生的EDI几乎是呼吸产生的2倍，即皮肤接触暴露更应引起防护关注。以上分析结果与段小丽等[27]的计算结论一致。

而对于职业人群，由于接触污染源的时间更长，其EDI总量大于非职业人群，呼吸贡献约84%，皮肤贡献约16%。因此，对于烧烤从业人员防护重点应在于佩戴口罩，加强厨房烟气排出等措施。

2.1.2 健康风险评价

应用式(5)计算得出，通过3种主要暴露途径对男性、女性、儿童和职业人群造成终身累计超额致癌风险(ILCR)总和为3.87×10-6、4.81×10-6、7.57×10-7、1.28×10-5。通常认为ILCR值10-6～10-4为致癌风险的最大可接受水平，因此，4类人群致癌风险均在可接受范围内，职业人群风险高于非职业人群风险一个数量级，儿童致癌风险可以忽略。

图2 不同人群苯并[a]芘等效剂量日均摄入量估算结果比较Fig. 2 Comparisons among estimated daily intake of B[a]P-eq dose of different population

2.2 内暴露模拟

尿液中的代谢标志物——羟基多环芳烃(OH-PAHs) 在我国很早就被用于评价焦化工人、交警等高暴露职业人群的PAHs暴露水平，也被用于儿童、学生等敏感人群[17-19,28]。然而，通过尿液中检测出的羟基多环芳烃浓度代表的PAHs内暴露量，有可能低估了真实暴露量，特别是对于5～6环结构的PAHs单体。首先，PAHs代谢物可以通过其他途经排出(粪便)或被其他组织(大脑)富集[8]；其次，在体内酶生物转化作用下，体内PAHs代谢物会以非羟基结构存在，并被储存或排出体外。因此，使用模型模拟的方法得到组织中暴露标志物的浓度，相对于最常用的尿液检测方法具有一定的优势。本文中，我们分别根据职业与非职业人群的暴露特点，给出了他们相应的组织内暴露标志物浓度模拟结果。图3为成年男性、成年女性、儿童和职业人群的各组织中暴露标志物(芘)浓度—时间曲线。

2.2.1 非职业人群

图3中，成年男性、成年女性、儿童各组织最大浓度分别是脂肪0.36、0.56、0.43 ng · L-1；皮肤0.64、0.75、0.71 ng·L-1；肾0.97、1.12、1.05 ng·L-1；肝6.52、7.29、8.67 ng·L-1和静脉血0.71、0.80、0.94 ng·L-1。3类非职业人群组织浓度排序均为：肝>肾脏>静脉血≥皮肤>脂肪。统计分析显示，成年女性组织浓度水平高于男性(P<0.01)，性别差异特点与Brown 等[22]报导的相同，儿童组织中只有脂肪浓度显著高于成人水平(P<0.01)，其余组织与成人无显著差别(P>0.05)。本文中，较大的Vfc和较小的CLint导致女性人群摄入量较大而代谢量较少，是导致模拟组织浓度大于男性的主要原因。综上，儿童和女性作为敏感人群在剂量评估过程中是值得特别关注的。

图3 PBPK模拟人群PAHs暴露后体内组织芘浓度时间变化Fig. 3 PBPK modelling of population tissues concentrations of pyrene

表4 人群健康风险评价Table 4 Health risk assessment of four subjects

2.2.2 职业人群

图3中，职业人群各组织最大浓度排序为：脂肪(2.97 ng·L-1)>皮肤(1.14 ng·L-1)≥肾脏(1.14 ng·L-1)>肝(0.57 ng·L-1)>静脉血(0.17 ng·L-1)。职业人群经连续8 h呼吸及皮肤接触暴露后，脂肪、皮肤和肾模拟浓度为非职业人群的1～2倍，静脉血浓度约为非职业人群的1/5，肝组织浓度约为非职业人群的1/13。

对于职业人群，脂肪是富集芘浓度最大的组织，这与赵肖等[16]报道的分布情况一致，可能的原因是血-脂肪组织分配系数高(pcf= 11.84)，造成脂肪组织是芘在人体内的“汇”。同时，由于芘的皮肤渗透系数(kp)较小，皮肤组织对空气中气相芘的吸附作用有限，皮肤血流量较小，导致皮肤接触途经对体内组织浓度的贡献相当小。

2.2.3 人体总富集量

对于一般暴露实验，组织样品难以获得，因此相对缺少组织的实测数据。然而，根据内暴露模拟结果，我们可以计算出暴露过程中人体组织对PAHs的总富集量(式14)。

A=Cf·Vf+Cs·Vs+Ck·Vk

+Cl·Vl+ 0.75 · CV ·Vb

(14)

其中，人体血液总体积(Vb)可用体重的7.9%换算得到，且静脉血占血液总体积的75%[12]。本研究中，一次暴露过程中人体组织芘总富集量(A)最大值分别为非职业人群6～11 ng·d-1，职业人群48 ng·d-1。同时我们注意到暴露对象在离开暴露环境或进食一段时间过后，A会较峰值时降低约80%。

2.2.4 小结

“膳食暴露”和“呼吸与皮肤接触暴露”过程体现出完全不同的代谢动力学特性：以膳食暴露为主的非职业人群组织浓度动力学主要体现为先上升再下降的“单峰”特性，表现出与Keys等[14]、Draper等[29]相似的趋势，以呼吸与皮肤接触暴露为主的职业人群组织浓度则在暴露期间很快达到最大浓度并保持稳定，结束暴露后迅速下降，动力学过程表现出与Kim等[30]相似的趋势(图4)。

2.3 不确定性分析

模拟结果的变异系数CV越大，则表明不确定性越大。图5为本研究模型的不确定性分析结果。本研究中，选取了预测的各组织暴露标志物浓度为模型不确定性的评价标准，非职业人群标志物浓度的平均CV为34%～36%，职业人群的平均CV为12%。与其他研究相比，我们的不确定性较低，这说明单独的参数不确定性对预测结果的不确定性贡献有限，应该注意到建模过程中对暴露源的理想化处理及人体组织、生化过程的简化会产生更大的不确定性。

图4 模拟得到人群暴露于航空燃油后呼出气中标志物萘浓度(Kim 等[30], 2007)Fig. 4 Modelling of exhaled air naphthalene concentration for occupational people exposed to jet fuel (JP-8) (Kim et al., 2007)[30]

图5 模型模拟浓度的平均变异系数Fig. 5 The average coefficient of variation (CV) of pyrene concentration in different models

我们还基于模型参数的概率分布表述(表4)，采用蒙特卡洛法评估了模型参数的不确定性贡献。基于芘总富集量的百分位分布(表5)可知，“膳食暴露”模型的最大的不确定参数来源是fabs。“呼吸及皮肤暴露”模型的参数不确定性源是IRinh，而当标志物的kp较小时，SA的参数不确定性作用可忽略。

综上所述：本文利用蒙特卡洛法与建立的本地化的PBPK模型分别定量模拟了不同人群暴露于烧烤源PAHs后内、外暴露剂量，终身致癌风险，标志物总富集量。结果表明，非职业人群B[a]P-eq日均摄入剂量为“男性≥女性>儿童”，膳食暴露贡献最大，各组织中芘最大浓度为0.36～8.67 ng·L-1，肝组织浓度最大；职业人群B[a]P-eq日均摄入剂量大于非职业人群，呼吸暴露贡献大于皮肤接触暴露，组织中芘最大浓度范围是0.17～2.97 ng·L-1，脂肪组织浓度最大。计算得到的致癌风险7.57×10-7～1.28×10-5，均在可接受范围内(10-6～10-4)。模型不确定性较低，fabs对“膳食模型”的不确定性贡献最大，IRinh对“呼吸及皮肤暴露模型”的不确定性贡献最大。

本文重点介绍了用于模拟组织内暴露标志物的模型方法，采用了大量本地化的暴露参数，比较分析了年龄、性别、职业等因素对PAHs暴露结果的影响。未来应在此基础上开展更深入的模拟研究，建立各种污染物的模型参数数据库，丰富模型结构，加强验证性实验等。

表4 模型参数概率分布表：每个参数独立用于模型的不确定性估算Table 4 Probabilistic distribution of model parameters: each parameter was used to assess model uncertainty respectively

注：U, 均匀分布；N(μ,σ2)，正态分布。Note: U, uniform distribution; N(μ, σ2), normal distribution.

表5 不确定性分析：基于模型参数和暴露参数概率分布计算的芘总富集量(ng·d-1)百分位值Table 5 Uncertainty analysis: the calculated percentiles of pyrene accumulated amounts (ng·d-1) based on probabilistic distribution of model and exposure parameters