Toxicity remission of PAEs on multireceptors after molecular modification through a 3D-QSAR pharmacophore model coupled with a gray interconnect degree method

In the proposed model, the gray interconnect degree method was employed to process the acute toxicity values of phthalate acid esters (PAEs) to green algae, daphnia, mysid, and fish (predicted by EPI Suite software) and to obtain the comprehensive characterization value of the multireceptor toxicity effect (MTE) of PAEs. The 3D-QSAR pharmacophore model indicated that hydrophobic groups significantly affected the MTE of PAEs. Based on this, 16 PAEs derivative molecules with significantly decreased comprehensive characterization value (more than 10%) of the toxic effects of multireceptors were designed. Among them, 13 PAEs derivative molecules reduced the toxicity values (predicted by the EPI Suite software) of four receptor organisms to varying degrees. Finally, two derivative molecules from PAEs were screened and could exist stably in the environment. The derivative molecule’s reduced toxicity to the receptor was obtained through molecular docking methods and simulated the PAEs’ primary metabolic response pathways. The above research results break through the pharmacophore model’s limitation of only being suitable for the single effect of pollutants. Its application provides a new theoretical verification basis for expanding the multieffect pharmacophore model.

organelles and antioxidant systems, resulting in cell deformities and inhibiting algae growth [21]. Long-term exposure to DEHP has a certain inhibitory effect on total reproductive mass, average reproductive mass, and population growth of the large salamander F3 generation [22]. PAEs also have harmful effects on the reproductive and endocrine systems of fish and crustaceans [23], and Patyna found that continuous exposure to low DBP concentrations seriously affect the fertility of Japanese sturgeon offspring [24]. The above literature mainly focuses on PAEs' acute toxicity in a certain organism, and only focuses on some PAE molecules in terms of their exposure pathway, toxicity, and performance of a single biological receptor. The research on the toxicity's molecular mechanism is insufficient.In view of PAEs' increasingly widespread application from an environmental pollution control perspective, it is important to carry out multireceptor low-toxicity activity PAE molecule joint regulation. Therefore, algae, invertebrates, and fish must be included in toxicological data, this article selects four common aquatic organisms (green algae, daphnia, mysid, and fish) that represent different nutritional levels in the water environment to study PAEs' comprehensive toxicity effects on four aquatic organisms and modify multireceptor low-toxicity PAE molecule.
QSAR, as a technology to quantitatively reveal compounds' toxicity and biological activity, can use the validated pharmacophore model [25]. Song et al. [26] proved a pharmacophore model to study the acute toxicity of six naphthoquinone compounds to daphnia magna. The results showed that the compounds' hydrophobicity had a great effect on receptor toxicity. Wang et al. [27] used hydrophobic groups to establish a pharmacophore model of the toxicity of perfluoro carboxylic acids to photobacterium, and the model regression coefficient was high. Qiu et al. [28] used a pharmacophore model to perform hydrophobic group substitution reactions on nine common PAE molecules and selected derivatives, with significantly enhanced Raman characteristic vibration spectra of PAEs. Jiang [29] proved the pharmacophore model could construct a POPs characteristic regulation scheme for PBDEs, and carried out modification designs of representative homologs, which confirmed the pharmacophore model's feasibility in molecular modification. Therefore, in this paper, studying the regulation scheme of the MTE of PAEs to multireceptors can be based on above pharmacophore model design method. In view of the limitation of the pharmacophore model's dependent variable as the pollutant's single pollution effect, this paper uses a grey interconnect degree [30] to deal with the aquatic receptors' toxicity values and calculate the toxicity comprehensive characterization values of PAEs to multireceptors. It is applied to the construction of the pharmacophore model of the MTE of PAEs and the modification design of multireceptor low-toxicity PAEs' derivative molecules, which provide a theoretical basis for constructing a multireceptor comprehensive toxicity effect model of PAEs.

Sources of data
The ECOSAR toxicity prediction module in EPI Suite software was used to predict the toxicity of 14 PAEs molecules to four organisms (green algae, daphnia, mysid, and fish),expressed as the concentration for a 50% maximal effect (EC50) or 50% lethal concentration (LC50), as shown in Table 2.

Calculation of comprehensive characteristic values of the MTE of PAEs using the gray interconnect degree method
A gray relation analysis (GRA) is a multifactor statistical analysis method, based on the similarity or dissimilarity of development trends between factors, which is used to measure the degree of correlation between factors [31]. The GRA results were obtained by the correlation between an indicator and factors that affect the indicator because this method involves longitudinal averaging of gray interconnect coefficients [32]. However, the PAEs molecules (indicators)' comprehensive toxicity effect on four aquatic organisms (factors) was studied in this paper. It was not necessary to obtain the order of the degree of influence between each factor and the indicator, which requires horizontal averaging of gray interconnect coefficients. Because the dimensions of the acute toxicity prediction values of PAEs are the same, no dimensionless processing was required. The acute toxicity classification standard (LC50/EC50 < 1.0 mg/L) was used as the reference sequence, X 0 , and the toxicity values of four organisms to green algae, daphnia, mysid, and fish were used as the comparison sequence, X i (i = 1,2,3,4), the weight of the four groups of comparison sequences was set to 25%. After obtaining the absolute difference between the corresponding points of the reference sequence, X 0 , and the comparison sequence, X i (i = 1,2,3,4, k = 1,2,3, … , 14), and substituting each column's maximum and minimum values of the absolute difference into Eq. (1) to calculate gray interconnect coefficients (ξ 0i (k)) of the four comparison sequences (X i ) and the reference sequence (X 0 ), where ρ is the resolution coefficient, ρ∈(0,1), and generally takes a value of ρ = 0.5.
Eq. (2) was used to calculate the average value of the gray interconnect coefficients horizontally to obtain the gray interconnect degree, y ok , of PAEs and four aquatic organisms. This was used as a comprehensive characterization of the MTE of PAEs, where n = 4. (2)

Construction method of the pharmacophore model of the multireceptor low-toxicity comprehensive effect of PAEs
The structural formulas of 14 PAE molecules were drawn by SYBYL-X2.0 software, entering the molecular construction mode from "sketch" in the toolbar, then optimizing the PAE molecules' force field after drawing, selecting the molecules' lowest energy conformation as the dominant stable conformation, optimizing each molecule's energy in the "Tripos" force field with the molecular program "minimize" and selecting "Gasteiger-Huckel" from the "charges" option menu. Using Powell's energy gradient method, "minimize details" was clicked, selecting the maximum number of repetitions (max. iterations) to 10,000, reducing the energy convergence limit (gradient) to 0.005 [33]. The "gradient" value is a termination criterion and, if the gradient difference calculates twice consecutively below this value, the calculation is terminated and the molecular structure optimization is completed.
The "3D-QSAR pharmacophore model generation" module in Discovery Studio 4.0 software was used to build the pharmacophore model [34]. The selected model's pharmacophore characteristics include: hydrogen bond donor (HBD), hydrogen bond acceptor (HBA), hydrophobic group (H), hydrophobic ring (HA), and ring aromatic (RA). The parameters for generating all molecular conformations were set as follows: "conformation generation" selected the best mode (best), the default energy cutoff of "energy threshold" was 20 kcal/mol, the "maximum conformations" was 255, the "minimum interfeature distance" was 1.5, the number of pharmacophore features was 0-5, and the energy threshold of each homolog to generate a similar conformation was 10, while other parameters adopted default values.
The "Hypo Gen" module in Discovery Studio 4.0 software was selected to evaluate the constructed pharmacophore model. "Cost function", one of the model's evaluation indicators, was used to express and evaluate the model's complexity and chemical characteristics as well as errors between each model's predicted values and experimental data. Each pharmacophore model had its own total consumption (total cost). According to Occam's Razor [35], the lower the "total cost" value, the closer it is to the "fixed cost" value, so the pharmacophore model is more reliable. "Configuration cost", another important parameter, is determined by the model's spatial complexity. The "configuration cost" value of a significant pharmacophore model should not be greater than 17 [36]. The larger the model correlation coefficient "R 2 " (> 0.7), the more predictive the pharmacophore model, and the more likely it is to meet the analytical needs [37]. In addition, "root mean square, " "fit value, " and "error" can be used as the pharmacophore model's evaluation indices.

Molecular docking and quantum chemical calculation methods
Molecular docking supposes that the binding between the ligand and the receptor conforms to the "lock and key principle", which satisfies the matching of spatial shape and energy, and finally obtains the optimal binding mode and stable composite conformation. Herein, the Lib-Dock quick docking method in Discovery Studio 4.0 softwarewas used [38]. The Poling algorithm performs a conformation search on the ligand molecule and then analyzes the binding site of the receptor and uses the grid-like algorithm to generate a series of polar and nonpolar hot spots. Finally the conformation and hot spots are matched with the energy and geometry to obtain the docking result. Considering that the crystalline water molecule at the binding site may affect the binding of ligand receptor, the water molecule at the protein binding site is eliminated when docking. "Find sites from receptor cavities" under the "Define" module determines the possible binding sites for ligand receptors, followed by selecting "user specified" in "Docking preferences", setting the maximum saved conformation to "10", and the rest of the parameters are the default values. The docking result is expressed using the "Lib-Dock score. " The magnitude of the value represents the strength of the binding ability.
The quantum chemistry calculations herein are based on the Gaussian 09 software package, the computer operating system used is Linux, and the Gaussian calculation results are displayed using the Gauss view 5.0 program. The DFT method is used to calculate the reactants, products, and transition state (TS) of the primary metabolic reaction at the B3LYP/6-31G(d) basis set level, and the reaction energy barrier (ΔE) of the primary metabolic pathway of PAEs molecules is calculated using Eq. (3). The TS has only one imaginary frequency, and the reaction path is verified through the intrinsic reaction coordinate [39].

Calculation of comprehensive characterization of the MTE of PAEs
Using the gray interconnect degree to process the original data (predicted by EPI Suite software), the absolute difference between the corresponding points of X 0 and X i , should be found, and the minimum and maximum values of each column's absolute difference should be obtained (Table 3, for calculation validity, retaining six significant digits). We substitute into Eq. (1) to calculate the gray interconnect coefficient ξ 0i (k) of each corresponding point and then obtain the gray interconnect degree y 0k of PAEs and four aquatic organisms from Eq. (2) ( Table 4).
Comprehensive characterization values of the MTE of PAEs are seen in Table 5.

Construction and evaluation of pharmacophore model of the multireceptor low-toxicity comprehensive effect of PAEs
Herein, 14 PAE molecules were divided into 10 training set molecules used for the construction of pharmacophore models, and four test set molecules were used to validate the pharmacophore models.The pharmacophore models with good performance parameters are listed in Table 6.
It was demonstrated that Hypo 1 had the best evaluation score among the five models. "Total cost" value (51.61) and "RMS" value (0.055) were the smallest, "total cost" was closest to "fixed cost" (33.636), "configuration" value was 16.834 (<17) and "correlation" (0.85) was closest to 1, the absolute value of "error" value of the molecules was less than 2 [40], and was within the tolerance range. Therefore, the model is significant and meets the requirements, so Hypo 1 was selected as the optimal pharmacophore model. The model has a hydrogen bond acceptor, a hydrophobic group, and a hydrophobic ring. The test set was used to verify this pharmacophore model, the "fit values" of the PAEs molecules were high, and the absolute "error" values were less than 2 (shown in Table 7), which shows that Hypo 1 had a stable prediction ability for PAE molecules other than the training set. Table 3. Absolute difference between X 0 (k) and X i (k).

Determination of substitution groups and substitution sites of target molecules, DINP, and DEHP, based on the
Hypo 1 optimal pharmacophore model DEHP, which has priority control of pollutants, and DINP, which has the largest comprehensive characterization of the MTE of PAEs in the training set, were selected as target molecules to determine the molecular modification site. Derivative molecules were designed based on this. Figure 1 shows   carboxyl oxygen atom of the DEHP molecule (shown in Figure 2). Therefore, introducing a hydrophobic group at branch positions can affect the PAEs' toxic activity. The positions of the substitution groups, introduced by DINP and DEHP, are shown in Figure 2; that is, molecular modification of this site was determined, which provides a basis for further screening derivative molecules for the MTE of PAEs.

Molecular modification of PAE derivatives based on the multireceptor low-toxicity pharmacophore model
Eleven common hydrophobic group were selected as substituent groups: methyl (-CH 3 ), ethyl (-CH 2 CH 3 ), propyl (-CH 2 CH 2 CH 3 ), vinyl (-CH=CH 2 ), phenyl (-C 6 H 5 ), methoxyl (-OCH3), hypochlorite (-Cl), fluoride (-F), bromo (-Br), sulfydryl (-SH), and nitro (-NO 2 ), to monosubstituted modification for DINP and DEHP, obtaining 22 modified derivative molecules. The constructed optimal pharmacophore model, Hypo 1, was used to predict the comprehensive characterization value of the MTE of PAE derivative molecules and compared with the toxicity comprehensive characterization values of corresponding target molecules,as shown in Table 8. The results showed that 16 PAE derivative molecules with toxicity comprehensive characterization values increased by more than 10%, including nine DINP derivative molecules (an increase of 11.95%-208.12%), and seven DEHP derivative molecules (an increase of 13.02%-48.07%), indicating that the toxicity of 16 derivative molecules was significantly lower than the target molecule.

Evaluation and verification of multireceptor comprehensive toxicity of PAE derivatives 3.5.1. Evaluation and verification of the MTE of PAE derivatives based on the EPI database
The ECOSAR module in the EPI Suite software was used to predict the above 16 PAE derivative molecules' toxicity values to multireceptor model (green algae, daphnia, mysid, and fish), taking the negative logarithmic values, as shown in , and the decline of multireceptors was close to 1:1:1:1. Therefore, a total of 13 PAE derivative molecules were screened, with a significant reduction in toxic activity.

Evaluation and verification of the MTE of PAE derivativesbased on the single receptor pharmacophore model
Based on the negative logarithm of toxicity values (predicted by EPI Suite software) of 14 PAEs on multireceptors as data sources, the abovementioned PAEs' toxicity comprehensive effect pharmacophore model construction method was used to construct green algae, daphnia, mysid, and fish' single receptor optimal pharmacophore models, as shown in Table 10.
The "configuration" values of the four pharmacophore models were 16.674 (<17), while the "correlation" values were all greater than 0.7. This shows that the established pharmacophore model exhibits stable prediction ability.  The above PAEs single receptor pharmacophore model was used to predict the PAEs derivative molecules' toxic activity on the corresponding receptors (negative logarithmic values, Table 11). Among these, DINP-C 6 H 5 and DEHP-F derivatives showed a consistent decrease in toxic activity on the multireceptor model, and this was consistent with the trend in the predicted value of the comprehensive effect pharmacophore model, which further verifies the reliability of the PAEs' multireceptor toxicity comprehensive effect pharmacophore model.

Evaluation of functional properties of PAE derivatives
The functional characteristics of PAE molecules include stability and insulation. The "total energy, " "energy gap" (which is the difference between the highest occupied orbital energy-E HOMO , and the lowest empty orbital energy-E LUMO of the molecule [41]), "frequency" of the target molecule, and its derivative molecules were calculated using Gaussian software. The "total energy" value represents stability, the "energy gap" value represents insulation, and the larger the energy gap value, the stronger the insulation. At the same time, the "frequency" value (>0) was used to evaluate whether the derivative molecules can exist stably in the environment [42]. As can be seen from Table 12, in the designed PAEs derivative molecules, the "total energy" value of DINP-C 6 H 5 decreased 6.34%, while the DEHP-F increased, but the increase was small (<5%), indicating the PAEs derivative molecules' stability had improved or remained unchanged compared with the target molecule. The "energy gap" value demonstrated a smaller change, indicating that the PAEs derivative molecules' stability had increased, while insulation was less affected. Both PAEs derivatives had "frequency" values greater than zero, indicating that their structures could exist stably in the environment.

Evaluation of persistent organic pollutants (POPs) properties of PAE derivatives
EPI Suite software was used to predict the bioaccumulation, long distance migration, and persistence of two PAE derivative molecules [43]. Based on POPs characteristic parameter values of DINP and DEHP, from the analysis in Table 13, the "LOGK OW " values of the two PAE derivative molecules had decreased by 6.56%-22.84%. The "LOGK OA " values had also decreased by 8.88%-12.56%, indicating that PAE derivative molecules had significant bioaccumulation and long distance migration in the environment. Because the PAE molecule itself is not a persistent organic pollutant (half-life in air>two days), the increase in the "half-life" value of the PAE derivative molecules had no significant effect on its persistence in the environment.

Analysis of toxicity mechanism of PAE derivatives based on molecular docking
Under external stress, a large number of oxygen radicals were generated in the algae cells. At this time, the antioxidant system was activated, and the peroxidase catalysis promptly removed a large amount of reactive oxygen species. PAEs caused oxidative damage to algae cells by acting on mitochondria (Mn-SOD) and cytoplasm (Cu/Zn-SOD) [44]. Glutathione (GSH) is an important antioxidant in living organisms. Copepods are rich in unsaturated fatty acids and,  when environmental stress exceeds the capacity of copepods, unsaturated fatty acids are degraded and lipid peroxidation occurs [45]. Chitinase is closely related to shrimp growth, food digestion, and disease defense [46], while PAEs may have toxic effects on shrimp by affecting chitinase gene expression [47]. Peroxisome proliferator-activated receptors (PPARs) control many intracellular metabolic processes. Among these, PPAR-α receptors are abundantly expressed in liver cells, and activation is a necessary condition for phthalate compounds to cause toxic hepatic reactionsin fish [48]. This article downloaded the MN-SOD enzyme crystal structure (1BA9), glutathione peroxidase crystal structure (3DWV), chitinase crystal structure (3ZXX), and PPAR-α protein crystal structure (3KDT) from the PDB protein database, which respectively represent green algae, daphnia, mysid, and fish receptors, docked with PAE molecules before and after modification, and expressed the recipient organism's toxic activity by the molecular docking ability [49]. The target molecules (DEHP, DINP) and derivative molecules (DEHP-F, DINP-C 6 H 5 ) were molecularly docked with the four enzyme proteins by using Discovery Studio 4.0 software and the corresponding scoring function values were calculated. The lower the scoring function value, the weaker the binding ability between molecules and enzyme protein, and the lower the toxic effect on the recipient organism. Table 14 shows that the scoring function values of the two PAE derivative molecules docking with four enzyme proteins were lower than the target molecules (a decrease of 3.9%-19.8%), indicating that the designed PAE derivative molecules had a weaker receptor binding ability, reducing toxicity to the receptor organism.
After the molecule binds to the receptor protein, it falls into the pocket formed by the amino acid residues around the receptor protein and reacts with the receptor mainly through hydrogen bonding, charge, or polar interaction, followed by the main chain and side chain of amino acids, generally interact with acceptor molecules in the form of hydrogen bonds. The docking results show that the main forces when PAEs and their derivatives bind to 3KDT and 1BA9 proteins are electrostatic force and van der Waals force, and the main forces, when they bind to 3DWV and 3ZZX proteins, include the electrostatic force, van der Waals force, and hydrophobic interaction. Compared with the target molecule, when the derivative molecule binds to the receptor protein, the number of surrounding amino acid residues that interact with it decreases and the binding ability becomes weaker. Therefore, it is possible to explain the decrease in the scoring function values for the binding of the derivative molecule to the receptor protein. When DINP-C6H5 binds to 3DWV, the number of surrounding amino acid residues that interact with 3DWV is larger than that of the target molecule; however, the scoring function value is lower. The possible cause for this is that when the target molecule is combined with 3DWV, it forms hydrogen bonds with amino acid residues TRPB137 and HOHB3105 and forms a π-bond interaction with TRPB137, whereas DINP-C6H5 only forms a π-bond interaction with TRPB137 when combined with 3DWV.

Analysis of toxicity mechanism of PAE derivatives based on metabolic response
The primary metabolite of PAEs, phthalate monoesters, has been detected in aquatic environments [50]. Ge Jian et al. [51] studied the metabolism of DEP, DBP, BBP, and DEHP in grass carp organs, and results showed that the main metabolites were corresponding phthalate monoesters. When studying the degradation products of black algae, Chen Bo [52] found that the phthalate monoesters, MBP and MEHP, were distributed in black algae. The primary metabolic pathway of PAEs in aquatic organisms is hydrolysis to the corresponding monoester compounds under the action of enzymes. According to these mimics the primary metabolic processes of PAEs, and their derivative molecules, in aquatic organisms,the products of DINP, DEHP, and derivative molecules DINP-C 6 H 5 , DEHP-F after primary metabolism are MINP, MEHP, MINP-C6H5, and MEHP-F ( Figure 3). Toxicity values of the primary metabolites of PAE target and derivative molecules (phthalate monoesters) to multireceptors were predicted using the EPI database (Table 15). The toxicity of PAE derivative monoester molecules was significantly lower than that of the target monoester molecules (green algae, 146.1%-2683.4%; daphnia, 125.5%-1880.3%; mysid, 188.4%-5070.1%; and fish, 112.6%-1473.4%).
Gaussian software was used to calculate the reaction energy barriers of the primary metabolic pathways of PAE target and derivative molecules and to determine whether the reaction could proceed, and how easy it was, by comparing the activation energy barriers in the transition states of the primary metabolic pathway response before and after molecular modification [53] as shown in Table 16. The reaction energy barrier of DINP was 51.77 KJ/MOL and DINP-C 6 H 5 was 5.96 KJ/MOL; a reduction of 88.5% compared with DINP. The reaction energy barrier of DEHP was 4.08 KJ/MOL and DEHP-F was 0.31 KJ/ MOL; a reduction of 92.4% compared with DEHP, indicating that the energy required for the first-order metabolism of PAE derivative molecules was greatly reduced compared with the target molecules, and the derivative molecules were more easily metabolized in aquatic organisms, causing the toxic activity of the organism to be reduced significantly.

Conclusion
Herein, the gray interconnect degree method assisted the PAE multireceptor low-toxicity effect and the pharmacophore model were established, passing validation of the traditional pharmacophore model and successfully applying it to the PAE multireceptor low-toxicity comprehensive effect of molecular modification. Based on the evaluation of the functional characteristics and POPs characteristics, two PAE derivative molecules were screened, with a reduction in comprehensive toxicity of 13.29% and 21.89%. Molecular docking and simulation methods of primary metabolic mechanisms in aquatic organisms confirmed the reason for the decrease in the multireceptor low-toxicity comprehensive effect of PAE derivative molecules. The method established in this article broke through the limitations of traditional pharmacophore models for single effect modeling of pollutants and provided theoretical support for building pharmacophore models that can simultaneously control multiple effects of pollutants.