CN118800357A - A predictive model for biomagnification of organic chemicals in low trophic levels of food chains - Google Patents

Abstract

The present invention discloses a low-trophic-level food chain biomagnification prediction model for organic chemicals, which is the first low-trophic-level food chain biomagnification prediction model for PBDEs and HBCDs based on a QSAR model. The model is obtained through steps such as sample collection and screening, molecular descriptor calculation, model construction, and model verification. The low-trophic-level food chain biomagnification prediction model for organic chemicals constructed using this method can accurately predict the biomagnification factors of organic pollutants in PBDEs and HBCDs chemicals, improve the accuracy of the prediction results, save manpower, material resources and time, is simple, fast and effective, and strictly follows the QSAR model usage rules stipulated by the OECD, explaining the key factors affecting the biomagnification factor from the molecular descriptor structure, which is of great significance to the risk control and environmental safety of toxic chemicals such as PBDEs and HBCDs.

Description

A predictive model for biomagnification of organic chemicals in low trophic levels of food chains

Technical Field

The invention relates to the field of environmental risk assessment of organic chemicals, and in particular to a low trophic level food chain biomagnification prediction model for organic chemicals.

Background Art

Toxic and hazardous chemicals are lipophilic, difficult to degrade, and accumulate and enrich in organisms. In an ecosystem, they exhibit a biomagnification effect as the trophic level increases, causing toxicity to higher organisms and humans. Therefore, studying the content characteristics of organic pollutants in aquatic organisms is of great significance to the health of local residents and food safety.

The degree of enrichment of a compound in high-trophic-level organisms and the human body through the food chain is of great significance for evaluating the ecological and environmental toxicity of the compound. For many years, scientists have been engaged in related work. Through long-term research, they have discovered and established compound accumulation models between various media. There are generally two standards for evaluating whether an organic pollutant has a bioaccumulation effect. The first is the K _OW (n-octanol-water partition coefficient) of the compound. Generally, when log _{K OW} >4-5, the compound may have a bioaccumulation effect, and when log _{K OW} is 5-7, the compound has the greatest bioaccumulation effect; the second is the bioaccumulation factor (BAF). The bioaccumulation factor (BAF) can characterize the relative bioaccumulation capacity of a compound.

The calculation formula of bioaccumulation factor (BAF) in aquatic food chains is as follows:

_Corganisms represents the concentration of pollutants in aquatic organisms, in pg/kg lw; _{Cdissolved phase in water} represents the concentration of pollutants dissolved in water, in pg/L. The concentration in organisms is the fat-normalized concentration.

In the BAF model, when the BAF value of a pollutant is higher than 5000 (or LogBAF>3.7), it can be considered that the pollutant has a bioaccumulation effect in the food chain; when the BAF value is between 2000 and 5000 (or LogBAF>3.7), it is considered to have a potential bioaccumulation effect.

PBDEs and HBCDs are a new type of persistent organic pollutants newly listed in the Stockholm Convention, with typical bioaccumulation. Due to their high lipophilicity and poor metabolic capacity, PBDEs and HBCDs have the potential to be amplified along the food chain. Recent studies have also confirmed that PBDEs and HBCDs can be transferred along the trophic level in the food chain. However, obtaining the biomagnification factor (BMF) of chemicals only through experimental methods is costly, time-consuming and labor-intensive, and it is difficult to meet the needs of ecological risk assessment of chemical substances. At present, there is still a lack of models for the biomagnification effect of chemicals in the food chain. Therefore, it is urgent to develop a scientific, fast and effective theoretical calculation method for food chain transfer models. The Organization for Economic Cooperation and Development (OECD) issued guidelines for the construction and verification of QSAR models in 2007, proposing the standards that QSAR models should meet: ① have clearly defined environmental indicators; ② have clear and clear mathematical algorithms; ③ define the application domain of the model; ④ the model has appropriate goodness of fit, robustness and predictive ability; ⑤ explain the model mechanism as much as possible.

Summary of the invention

The present invention, for the first time, constructs a low trophic level food chain biomagnification prediction model for PBDEs and HBCDs based on the QSAR model, as shown below:

PEC _{oral,predator} = PEC _water *BCF _fish *BMF (1)

BMF＝-5.04472+0.8374*GGI3-35.46426*Mor21v (2)

The meaning of each parameter is as follows:

The present invention provides a method for constructing the above-mentioned low trophic level food chain biomagnification model, which is specifically as follows:

⑴Sample collection and screening

The laboratory data obtained contains the biomagnification factor data of 9 organic compounds, which cover organic substances such as PBDEs and HBCDs. In order to establish an effective QSAR model, the data set is first divided into a training set and a test set. To ensure the representativeness of the training set compounds, the grouping method used in this work is the Kennard & Stone method. This method can avoid the uneven distribution of training set samples to a certain extent, and can well divide the data set into training set and test set. This data set is divided into 6 training sets and 3 validation sets.

⑵ Calculation of molecular descriptors

This method first constructs the molecular structures of 9 organic compounds in ChemDraw software, and then imports them into the HypeChem program to optimize the molecules. The optimization is divided into two steps: first, the MM+ molecular force field method is used for preliminary energy optimization, and then the semi-empirical quantum mechanics AM1 method is used to optimize the structure more accurately. The optimized structure is imported into DRAGON5.4 software to calculate 1664 different types of theoretical molecular descriptors. These descriptors were preprocessed before modeling, that is, constant terms, terms close to constants, and molecular descriptors with high correlation (the correlation coefficient of two molecular descriptors with a smaller correlation coefficient with the target value is greater than 0.96) were deleted. Finally, the remaining 1169 descriptors were used for the subsequent variable selection process.

⑶Model construction

A genetic algorithm was used to select a set of descriptors that were highly correlated with bioaccumulation, and this process was implemented in MobyDigs. After genetic algorithm variable selection, a linear QSAR model, namely the GA-MLR model, was established using the multivariate linear regression (MLR) method. The model evaluation function selected the leave-one-out cross validation, that is, when the performance of the model did not change significantly after adding a descriptor (the increase in Q2 by adding a descriptor was less than 0.02), the optimal number of descriptors was reached. In this method, the optimal number of descriptors was 7. The relevant parameters in the modeling were set as follows: population size was 100, the maximum allowed variables allowed in the initial model was 7, the mutation trade-off (T) was 0.5, and the crossover and mutation probabilities were based on the T parameter.

(4) Model verification

After genetic algorithm variable selection, the linear QSAR model, namely the MLR model, was established using the multivariate linear regression method. The linear MLR equation is as follows:

BMF＝-5.04472+0.8374*GGI3-35.46426*Mor21v

n _tr ＝6Q ² _LOO ＝0.9981R ² _fitting ＝0.9996R ² _adj ＝0.9994RMSE _tr ＝0.0665R ² _boot ＝0.7175

n _ext ＝3R ² _ext ＝0.836, Q ² _ext ＝0.8662R ² _adj ＝0.9994RMSE _ext ＝0.0289

Among them, GGI3 represents the third-order topological charge index, Mor21v represents 3D-MoRSE-weighted atomic van der Waals volume, GGI3 is positively correlated with the biomagnification factor, Mor21v is negatively correlated with the biomagnification factor, the RMSE of the training set and the validation set are 0.0665 and 0.0289, respectively, and the model prediction effect is good.

Table 1 Experimental and predicted values of BMF for PBDEs and HBCDs

The advantages of the low-trophic-level food chain biomagnification prediction model for organic chemicals established by this method are: the concentration of pollutants in aquatic organisms, the content level of PBDEs and HBCDs, the concentration of pollutants dissolved in water bodies, and then the calculated biomagnification time of chemicals in the food chain is long and the cost is high. The low-trophic-level food chain biomagnification prediction model for organic chemicals constructed by this method can accurately predict the biomagnification factors of organic pollutants such as PBDEs and HBCDs chemicals, improve the accuracy of the prediction results, save manpower, material resources and time, is simple, fast and effective, and strictly follows the QSAR model usage rules stipulated by the OECD, explaining the key factors affecting the biomagnification factor from the molecular descriptor structure, which is of great significance to the risk control and environmental safety of toxic chemicals such as PBDEs and HBCDs.

BRIEF DESCRIPTION OF THE DRAWINGS

Fig. 1BMF prediction model fitting diagram.

Fig. 2 Representation diagram of the application domain of the BMF prediction model.

DETAILED DESCRIPTION

The present invention is further described below in conjunction with specific embodiments.

Example 1

The specific steps for constructing a low-trophic-level food chain biomagnification prediction model for organic chemicals are as follows:

(1) Data collection, setting training set and validation set sample compounds

The laboratory obtained the biomagnification factors (BMFs) of 9 organic chemicals. A total of 6 sample compounds were selected for the training set and 3 sample compounds were selected for the validation set.

(2) Calculate descriptors

The compound structures were pre-optimized using MM+ molecular mechanics in Hyperchem 7.0 software, and the compound structures were optimized using the semi-empirical AM1 method. Based on the optimized structures, descriptors were calculated using Dragon 5.4 software, and a preliminary screening of the 1664 calculated descriptors was performed.

(3) Model construction

Genetic algorithm (GA) in MobyDigs software was used for variable selection. Based on the selected variables, a prediction model was established using the multivariate linear regression (MLR) method, namely the GA-MLR model:

BMF＝-5.04472+0.8374*GGI3-35.46426*Mor21v

Among them, GGI3 represents the third-order topological charge index, Mor21v represents 3D-MoRSE-weighted atomic van der Waals volume, GGI3 is positively correlated with the biomagnification factor, Mor21v is negatively correlated with the biomagnification factor, and the RMSE of the training set and validation set are 0.0665 and 0.0289, respectively.

(4) Model verification

According to the OECD guidelines for QSAR models, the constructed model needs to be internally validated (goodness of fit and robustness assessment) and externally validated (predictive ability assessment). The square of the correlation coefficient (R ² _adj ) and the root mean square error (RMSE) between the corrected experimental value and the fitted value are used to characterize the goodness of fit of the model:

Where n represents the number of compounds, m is the number of predictor variables, _yi and represent the experimental value and predicted value of the activity index of the i-th compound respectively; It is the average value of the experimental values of the compound activity index.

The cross validation coefficient (Q ² _LOO ) and Bootstrapping method (Q ² _BOOT ) were used to characterize the stability of the model:

in, It represents the average experimental value of the activity index of the training set compound. The bootstrapping method uses 1/5 cross validation and is repeated 5000 times.

External validation correlation coefficient (Q ² _EXT ), R ² _EXT , RMSE _EXT were used to characterize the prediction ability of the model:

Where n _EXT represents the number of compounds in the validation set, Represents the average value of the experimental value and the predicted value of the activity index of the validation set compound. The characterization and evaluation parameters of the model are obtained:

n _tr ＝6Q ² _LOO ＝0.9981R ² _fitting ＝0.9996R ² _adj ＝0.9994RMSE _tr ＝0.0665R ² _boot ＝0.7175

n _ext ＝3R ² _ext ＝0.836, Q ² _ext ＝0.8662R ² _adj ＝0.9994RMSE _ext ＝0.0289

Where n _tr and n _ext are the number of compounds in the training set and validation set, respectively, and p is the significance level. ^{R 2} _adj is the determination coefficient corrected for degrees of freedom; RMSE is the root mean square error; Q ² _LOO is the leave-one-out cross-validation coefficient; Q ² _BOOT is the validation coefficient of the Bootstrapping method; R ² _ext is the correlation coefficient between the experimental value and the predicted value, Q ² _ext is the external validation determination coefficient, and RMSE _ext is the root mean square error of the validation set.

The results show that the model has good predictive ability and robustness.

Example 2

This embodiment characterizes the application domain of the above prediction model.

The Williams plot is a model application domain defined by the standard residual (δ) and the leverage value (represented by h _i , where i represents different compounds). δ is calculated using the following formula:

The leverage value (leverage, h _i ) of the training set compounds can be calculated using the following formula:

h _i ＝ x _i ^T (X ^T X) ^–1 x _i (8)

Where _xi is the row vector of the molecular structure descriptor of the i-th compound. The warning value (h ^* ) is defined as:

h ^* = 3(k + 1)/n (9)

Among them, k is the number of descriptors and n is the number of training sets.

The results of the model application domain characterization are shown in Figures 1 and 2. In Figure 1, h ^* = 3(k+1)/n = 3(2+1)/6 = 1.5. The vertical axis of the Williams graph uses the standard residuals of the experimental values and the predicted values to characterize the degree of dispersion of the experimental values. When the absolute value of the standard residual δ of the compound is greater than 3.0, it is considered an outlier. The horizontal axis represents the h _i value of the compound in the training set. When h _i is greater than the warning value (h* = 1.5), it means that the substructure of the substance appears less frequently in the training set, which will have a significant impact on the model prediction results.

As can be seen from the figure, the leverage values h of all compounds are within the warning leverage value h ^* , indicating that the structure of this compound has a certain similarity with the structure of the training set compound. The standard residuals are all in the range of (-3, +3), indicating that this model is suitable for the prediction of BDE153 (CAS: 68631-49-2), BDE 183 (CAS: 207122-16-5) and BDE-99 (CAS: 60348-60-9). The standard residuals all fall within the range of (-3, +3), indicating that this model is suitable for the prediction of these three substances and can be well predicted, as shown in Figure 2.

Example 3

The model constructed in Example 1 was used to predict the food chain amplification factors of 9 persistent organic pollutants. The results are shown in Table 2. The R ² _adj of the model = 0.999, indicating that the model has a strong fitting ability. Q ² _LOO = 0.9981, Q ² _BOOT = 0.7175, indicating that the model is relatively robust. R ² _ext = 0.9994, Q ² _ext = 0.8662. Golbraikh et al. believed that the acceptable standards for QSAR models are Q ² > 0.50 and R ² > 0.60. The results show that the model has good predictive ability and can be successfully applied to compounds outside the training set. Figure 1 is a fitting diagram of the predicted value and experimental value of the food chain magnification factor (BMF) of persistent organic chemicals. It can be seen from Figure 1 that the predicted value and experimental value of most substances fit well, and BDE153 (CAS: 68631-49-2), BDE183 (CAS: 207122-16-5) and BDE-99 (CAS: 60348-60-9) can be well predicted.

Example 4

The biomagnification factor (BMF) of BDE 203 (SMILES: BrC1=C(OC2=CC(Br)=C(Br)C(Br)=C2Br)C(Br)=CC(Br)=C1Br) was predicted using the model constructed in Example 1. First, based on the molecular structure of the chemical substance, two descriptors, GGI3 and Mor21v, were calculated using the Dragon software, which were 1.937 and -0.06, respectively, and Hat was 0.435, which was within the application domain of the model.

BMF＝-5.04472+0.8374*GGI3-35.46426*Mor21v

BMF＝-5.04472+0.8374*(1.937)-35.46426*(-0.06)＝3.25

The predicted value of BDE 203 (BMF) is 3.25, which is close to the experimental result (1.73).

Example 5

The biomagnification factor (BMF) of BDE 196 (SMILES: BrC1=CC=C(OC2=C(Br)C(Br)=C(Br)C(Br)=C2Br)C(Br)=C1Br) was predicted using the model constructed in Example 1. First, based on the molecular structure of the chemical substance, two descriptors, GGI3 and Mor21v, were calculated using the Dragon software, which were 2.375 and -0.016, respectively, and Hat was 0.651, which was within the application domain of the model.

BMF＝-5.04472+0.8374*GGI3-35.46426*Mor21vBMF＝-5.04472+0.8374*(2.375)-35.46426*(-0.016)＝2.05, so the predicted value of BDE 196 (BMF) is 2.05, which is close to the experimental measurement result (1.43).

Claims (6)

1. A low nutrient grade food chain bio-scale predictive model of an organic chemical, characterized in that the low nutrient grade food chain bio-scale predictive model is as follows:

PEC_{oral,predator} ＝ PEC_water*BCF_fish*BMF (1)

BMF＝-5.04472+0.8374*GGI3-35.46426*Mor21v (2)

Wherein PEC _{oral,predator} refers to the in vivo concentration of low nutrient level predators, mg.kg _wetfish ^-1;PEC_water refers to the predicted concentration in water, mg/L; BCF _fish refers to the biological enrichment factor of fish, L.kg _wetfish ^-1; BMF refers to biological amplification factor; GGI3 refers to the 3rd order topological charge index; mor21v refers to 3D-MoRSE-weighted atomic Van der Waals volumes.

2. The method for constructing a low nutrient grade food chain bio-amplification predictive model of an organic chemical according to claim 1, comprising the steps of:

(1) Sample collection and screening;

(2) Calculating a molecular descriptor;

(3) Constructing a model;

(4) And (5) model verification.

3. The method for constructing a low nutrient level food chain bio-scale predictive model of an organic chemical according to claim 2, wherein step (1) is specifically: laboratory data obtained biological amplification factor data containing 9 organic compounds, which cover the organics of PBDEs and HBCDs, were first divided into training and test sets using the grouping method Kennard & Stone method, which was divided into 6 training and 3 validation sets.

4. The method for constructing a low nutrient level food chain bio-scale up prediction model of an organic chemical according to claim 2, wherein the step (2) is specifically: firstly, constructing molecular structures of 9 organic compounds in ChemDraw software, and then introducing HypeChem programs to optimize the molecules; the optimization is divided into two steps: firstly, performing preliminary energy optimization by using an MM+molecular force field method, then performing more accurate configuration optimization on a structure by using a semi-empirical quantum mechanical AM1 method, and introducing the optimized structure into DRAGON5.4 software to calculate 1664 theoretical molecular descriptors of different types; these descriptors are preprocessed prior to modeling, i.e., constant terms, near-constant terms, and molecular descriptors with high correlation are deleted, and 1169 descriptors are finally left for the following variable selection process.

5. The method for constructing a low nutrient level food chain bio-scale up prediction model of an organic chemical according to claim 2, wherein the step (3) is specifically: using a genetic algorithm to select a set of descriptors that have a high correlation with biological enrichment, this process is implemented in MobyDigs; after genetic algorithm variable selection, a linear QSAR model is established by a multiple linear regression MLR method, and model evaluation function selection is subjected to one-way interactive inspection, namely when the performance of the model is not obviously changed after one descriptor is added, the optimal descriptor number is reached; in the method, the number of the optimal descriptors is 7; the relevant parameters in modeling are set as follows: the population size is 100, the allowed number of variables to be made large by the initial model is maximum allowed variables, the variation equilibrium value is mutation trade-off, T is 0.5, and the cross-over and variation mutation probabilities are based on the T parameter.

6. The method for constructing a low nutrient level food chain bio-scale predictive model of an organic chemical according to claim 2, wherein step (4) is specifically: after genetic algorithm variable selection, a linear QSAR model, namely an MLR model, is established by a multiple linear regression method, and a linear MLR equation is as follows:

BMF＝-5.04472+0.8374*GGI3-35.46426*Mor21v

n_tr＝6Q² _LOO＝0.9981R² _fitting＝0.9996R² _adj＝0.9994RMSE_tr＝0.0665R² _boot＝0.7175

n_ext＝3R² _ext＝0.836,Q² _ext＝0.8662R² _adj＝0.9994RMSE_ext＝0.0289

Wherein GGI3 represents 3-order topological charge index, mor21v represents 3D-MoRSE-weighted atomic Van der Waals volume, GGI3 has positive correlation with biological amplification factor, mor21v has negative correlation with biological amplification factor, training set and validation set RMSE are 0.0665 and 0.0289 respectively.