QSAR模型预测石化废水中芳香族物质对厌氧菌群的综合毒性

吴少奇，王黎,*，夏正海，孙义，张春雨，鲁逸飞

1. 武汉科技大学资源与环境工程学院，环境污染绿色控制与修复技术研究中心，武汉 430081 2. 中韩(武汉)石油化工有限公司发展技术部，武汉 430082

石化、焦化和钢铁制造等行业产生的工业废水成分复杂，且普遍存在着芳香化合物[1]。芳香化合物是一类具有持久性、生物累积性和毒性的有机污染物[2]，对自然环境和人体健康有潜在威胁。含芳香化合物的有机废水处理方式主要包括物理化学法(如超声波处理、光化学氧化和电化学氧化等)[3]和生物法(如好氧降解和厌氧消化等)。其中，厌氧消化很好地克服了物理化学法和好氧降解所面临的高成本、污染物降解不充分以及二次污染等问题[4]。但厌氧微生物对环境较为敏感，高浓度的芳香化合物废水会对其活性产生明显的抑制作用，从而影响出水水质。因此，确定芳香族化合物对厌氧微生物的抑制阈值很有现实意义。

通过传统的逐一毒性试验获取有机物毒性数据，不仅工作量大而且费时费力，已无法满足当下芳香化合物废水处理的实际需求。定量构效关系(quantitative structure activity relationship, QSAR)模型是一种通过统计学手段建立物质结构参数与有机体活性定量关系的分析方法，可实现对有机物性质的有效预测[5]。同时，芳香化合物的毒性大小在很大程度上受其结构参数(如物理化学参数和量子化学参数)的影响[6]。因此，通过建立芳香化合物结构与其毒性关系的QSAR模型，可探究影响芳香化合物毒性大小的主要结构参数，从而合理预测芳香化合物的毒性。

为预测芳香化合物对厌氧微生物的毒性作用，笔者进行了厌氧微生物毒理试验，以19种芳香化合物的48 h半数效应浓度(48 h-EC50)值为样本，基于20个理化和量子化学结构描述符，利用多元逐步线性回归(multiple stepwise linear regression, MSLR)法建立了芳香化合物毒性与结构参数相关联的QSAR模型。

1 材料与方法(Materials and methods)

1.1 实验材料

1.1.1 模拟废水

采用模拟废水，其成分根据石化废水处理某工段进水水质，化学需氧量(COD)为4 000 mg·L-1，用葡萄糖、NH4Cl和KH2PO4按m(COD)∶m(N)∶m(P)=200∶5∶l的比例配制而成，同时加入微量元素Zn2+、Cu2+、Mn2+、Co2+和Ni2+。

1.1.2 接种污泥

试验所用的接种污泥来源于武汉某污水处理厂的厌氧污泥，接种污泥的挥发性固体(VS)、总固体(TS)和pH值分别为2.29%、4.32%和7.1。采集的厌氧污泥先在1.1.1所述的模拟废水中驯化20 d，之后用于厌氧微生物毒理试验。

1.1.3 芳香化合物

用于厌氧微生物毒理试验的芳香化合物(分析纯)共19种，可分为3类：苯类、苯酚类以及苯胺类。苯类包括苯、邻二氯苯、氯苯和对二甲苯；苯酚类包括苯酚、邻甲基苯酚、间甲基苯酚、对甲基苯酚、邻硝基苯酚、对硝基苯酚、间硝基苯酚、3，5-二甲基苯酚、2，4-二氯苯酚、2,6-二氯苯酚、邻(对)苯二酚和五氯苯酚；苯胺类包括苯胺和2,6-二乙基苯胺。

1.2 厌氧微生物毒理试验

在250 mL的消化瓶中加入200 mL模拟废水和接种污泥所构成的混合液，使初始的挥发性悬浮固体(VSS)浓度为15 g·L-1，同时调节pH为7.1±0.1。在消化瓶中分别添加1.1.3中所述的芳香化合物，使瓶中芳香化合物的浓度分别为100、500、1 000、2 000和4 000 mg·L-1。用氮气在液面上部吹扫2 min，然后用胶塞密封瓶口，保持瓶内处于厌氧状态。将消化瓶置于(35±1) ℃的恒温震荡器(ZQTY-70N，上海知楚，中国)中厌氧消化48 h，测定其最终累积产气量。每组实验设置2个平行样。

1.3 分析方法

1.3.1 测定方法

TS、VS、VSS以及pH按照美国公共卫生协会(APHA)的标准方法进行测定[7]；累积产气量采用排水法进行测定[8]。

1.3.2 数据分析方法

通过式(1)计算48 h芳香化合物对厌氧微生物的抑制率，然后将得到的48 h各组不同浓度的抑制率代入Logistic方程进行模拟，获得48 h-EC50值[9]。抑制率的计算式如下[10]。

(1)

式中：I为芳香化合物对厌氧微生物活性的抑制率(%)；GA为实验组t时刻厌氧微生物累积产气量(mL)；GB为对照组t时刻厌氧微生物累积产气量(mL)。

1.4 QSAR模型构建

1.4.1 描述符

为了鉴定影响芳香化合物对厌氧微生物毒性的结构参数，选择了20个具有代表性的描述符，涉及理化参数包括正辛醇/水分配系数(logKow)、折射率(n)、水溶解度(logS)、生成焓(heat of form, HOF)、摩尔折射率(molar refractivity, MR)[11]，同时涉及量子化学参数，包括最高占据轨道能(EHOMO)、最低空轨道能(ELUMO)、次最高占据轨道能(ENHOMO)、次最低空轨道能(ENLUMO)、分子偶极矩(dipole moment,μ)、分子总能量(total energy, TE)、分子摩尔体积(molecular volume, MV)、极化率(polarizability,α)、自旋极化率(spin polarization, SP)、零点(zero point energy, ZPE)、超极化率(hyper polarizability,β)、热力学能(thermodynamic energy, ThE)、热容量(heat capacity,HC)、熵(entropy,En)、自洽场能(self-consistent field energy, SCFE)[12]。理化参数通过Chemoffice查询获得，量子化学参数通过Gaussian 09W软件的密度泛函理论杂化泛函(DFT-B3LYP)方法，在6-311G(d, p)水平下对芳香化合物结构优化所得[13]。

1.4.2 QSAR模型构建与分析

在上述3类芳香化合物中随机各抽取一种组成预测集，用于检测模型的预测能力，其余的芳香化合物作为训练集。为了减少数据的跨度，将EC50值转换为10为底的对数值-lgEC50。采用MSLR建立-lgEC50与结构参数的芳香化合物毒性预测模型。

(2)

(3)

模型的适用范围由Williams图(即标准化残差与杠杆值)来表征[15]。计算标准化残差(δ)、杠杆值(h)以及警告杠杆(h*)的公式如下。

(4)

(5)

h*=3(p+1)/n

(6)

式中：p为模型中包含的描述符数量；xi为化合物i的描述符行向量；X为训练集的描述符矩阵。

标准化残差的绝对值>3的化合物被视为异常值。若预测的芳香化合物残差在此范围内，则说明模型的预测是可靠的。采用Origin 2017和IBM SPSS Statistics 25进行上述统计分析。

2 结果与讨论(Results and discussion)

2.1 芳香化合物毒性比较

芳香化合物对厌氧微生物的毒性大小可由48 h-EC50值表示，上述芳香化合物的48 h-EC50值范围为486～2 224 mg·L-1(图1)。EC50值反映了厌氧微生物所能忍受的毒性大小，EC50值越低，则芳香化合物的毒性越高。因此可知，五氯苯酚对厌氧微生物的毒性作用最强，对苯二酚则最弱。

图1 芳香化合物的48 h半数效应浓度(48 h-EC50)值Fig. 1 48 h median effect concentration (48 h-EC50) value of aromatic compounds

从毒理试验结果来看，苯类物质对厌氧微生物的毒性高于苯酚类物质，这与前人所得结论相近[16]。一方面，苯类物质比苯酚类物质具有更强的疏水性，导致前者能更轻易地与微生物细胞相结合，对其产生毒性作用。另一方面，有机化合物的生物降解过程可以间接影响其毒性的作用效果。苯酚类物质的结构比苯类物质更为复杂，导致其难以被厌氧微生物降解。因此，最先与微生物作用的是苯类物质，其毒性作用在反应初期就显现出来，苯酚类物质则与之相反。

2.2 芳香化合物毒性预测的QSAR模型

2.2.1 描述符的筛选和QSAR模型的构建

为了避免忽略具有低相关性但高重要性的描述符，将1.4.1中的20个描述符与训练集中16种芳香化合物进行MSLR分析，回归结果如表1所示。

表1 多元逐步线性回归的芳香化合物定量构效关系(QSAR)模型Table 1 The quantitative structure activity relationship (QSAR) model of aromatic compounds for multiple stepwise linear regression

-lgEC50=1.065+0.196logKow+0.154MR+0.028μ(7)

n=16，R2=0.928，F=47.581，P=0.000

2.2.2 模型检验

对于内部验证，式(1)的R2为0.928，逐步回归的显著性P=0.000<0.05，F=47.581>临界值F0.05(3,14)=3.34，说明该模型与训练集拟合优度较好，-lgEC50与logKow、MR和μ所构成的回归方程可靠度很高。从表2可知，logKow、MR和μ三者之间的相关系数分别为0.745、0.116和0.188，均<0.8[17]，且模型的VIF值介于l～10之间，这说明3个自变量之间没有潜在的共线性，故从统计学角度所构模型可以接受。模型的R2(>0.6)，F(>15)，P(<0.05)和VIF(<10)值均满足最小统计稳健性标准[18-19]。

2.2.3 外部预测能力

图2 -lgEC50的实验值与计算值相关性图Fig. 2 Correlation between experimental and calculated values of -lgEC50

此外，利用Williams图分析了QSAR模型对芳香化合物毒性预测的适用范围(图3)。水平虚线表示标准化的剩余异常值(±3)，垂直虚线表示警示杠杆值(h*=0.750)。由图3可知，训练集和测试集中的所有数据点都属于适用范围，表明QSAR模型对于通过-lgEC50与结构参数的线性关系预测芳香化合物对厌氧微生物的毒性是可靠的。

图3 QSAR模型的Williams图Fig. 3 Williams plot of the QSAR model

表2 回归方程中各自变量系数相关矩阵及其方差膨胀因子(VIF)Table 2 Correlation matrix of each variable coefficient and variance inflation factor (VIF) in the regression equation

表3 模型参数以及模型验证结果Table 3 Model parameters and model verification results

2.3 芳香化合物的毒性与描述符的作用机理

QSAR模型表示为式(7)，表明芳香化合物对厌氧微生物的毒性影响与logKow、MR和μ有关，且-lgEC50均与三者成正比。从表2可知，logKow、MR和μ对此模型的贡献程度分别为51.30%、27.63%和21.07%，说明logKow可能是影响芳香化合物对厌氧微生物毒性大小的主要因素。

由图4(a)可知，logKow与-lgEC50值相关性显著(r2=0.8105)。芳香化合物的毒性随logKow增大而增大。logKow反映了芳香化合物的溶脂性大小，而微生物的细胞膜是由磷脂双分子层所构成。芳香化合物的logKow越大，其溶脂能力越强，可能越容易通过细胞膜，从而表现出的毒性作用也就越大。

图4 模型参数与-lgEC50的相关性图注：(a) logKow与-lgEC50的相关性；(b) MR与-lgEC50的相关性；(c) μ与-lgEC50的相关性。Fig. 4 Correlation between model parameters and -lgEC50Note: (a) the correlation between logKow and -lgEC50; (b) the correlation between MR and -lgEC50; (c) the correlation between μ and -lgEC50.

MR是由于在光的照射下分子中电子云相对于分子骨架相对运动的结果，可作为分子中电子极化率的度量[21]。由图4(b)可知，MR与-lgEC50值相关性显著(r2=0.7128)，芳香化合物对厌氧生物的毒性随MR增大而增大。

μ体现了芳香化合物影响微生物细胞极性的能力。由图4(c)可知，μ与-lgEC50值相关性不显著(r2=0.0695)。Su等[19]探究了四胺基大环配体(TAML)/H2O2)氧化芳香族有机污染物的定量构效关系，发现所构多元线性模型的自变量#H:C与因变量lnkobs的相关性并不显著。因此，应将μ与其他因素结合起来考虑。

芳香化合物的浓度会显著影响厌氧消化过程中厌氧微生物活性，从而引起出水水质的波动。本研究通过厌氧微生物毒理试验，以芳香化合物48 h-EC50值为样本，基于理化和量子化学结构描述符，采用MSLR建立了QSAR模型，得出如下结论。

(1)QSAR模型为：-lgEC50=1.065+0.196logKow+0.154MR+0.028μ，R2=0.928，此模型可用于芳香化合物的毒性大小预测(h*=0.750)，有助于提高石化类工业废水的厌氧生物处理效率。

(2)芳香化合物对厌氧微生物的毒性影响与logKow、MR和μ有关，三者对模型的贡献程度分别为51.30%、27.63%和21.07%，logKow可能是影响芳香化合物对厌氧微生物毒性大小的主要因素。

(3)-lgEC50值与logKow和MR相关性显著，r2分别为0.8105和0.7128；-lgEC50值与μ相关性不显著，需要与其他因素综合起来考虑。