PAEs Derivatives’ Design for Insulation: Integrated In-Silico Methods, Functional Assessment and Environmentally Friendly Molecular Modification

As a common substance in production and life, phthalic acid esters (PAEs), the main component of plastics, have brought more and more serious problems to the environment. This study normalized the insulation, toxicity, and bioconcentration data of 13 PAEs to eliminate the dimensional coefficients of each index, and then used the comprehensive index method to calculate the comprehensive effect value of PAEs with three properties. The comprehensive effect value was used as the data source to construct the 3D-QSAR model of PAE molecular comprehensive effect. The DAP was selected as the target molecule, the distribution of each force field in the three-dimensional equipotential map was analyzed, and 30 molecular modification schemes were created. The constructed single-effect models of insulation, toxicity, and bioconcentration of PAEs and the scoring function module of DS software were used to evaluate the stability and environmental friendliness of PAE derivative molecules. Four PAE derivatives were screened for increased comprehensive effects, enhanced insulation, and reduced toxicity and bioconcentration. By calculating the binding energy of the target molecule and the derivative molecule with the degrading enzyme under different applied electric fields, it was found that the binding energy of DAP-1-NO2-2-CH2C6H5 decreases more than DAP does when there is an applied electric field, indicating that the degradation ability of degrading enzymes on PAE derivative molecules is reduced, which indirectly proves that the insulation is enhanced. The innovation of this paper lies in the insulation, toxicity, and bioenrichment data of PAEs being processed by mathematical method for the first time, and PAEs with high insulation, low toxicity, and low bioconcentration were designed by building a comprehensive model.


Introduction
Plastics are polymeric compounds formed through polyaddition or polycondensation and present many advantageous properties such as light weight, stability, high insulation, impact resistance, and abrasion resistance [1]. With the rapid development of society and substantial improvement in the living standards of people, plastic products have been widely used and have gradually become an indispensable part of daily life, work, and travel [2,3]. For different occasions and environments, plastic products should meet certain performance indicators [4,5]. For example, in cable materials, in addition to meeting the requirements of mechanical properties, heat resistance, and abrasion resistance, plastic products should show good electrical insulation [6]. Studies have shown that adding plasticizers to plastic products can greatly reduce the insulation of plastics [7]. Nagy [8] et al. reported the effect of DOP on PVC dielectric properties, indicating that the increase of plasticizer is proportional to the insulation of PVC.
Phthalic acid esters (PAEs) are also called phthalates and are usually used as plasticizers in plastic products to enhance their ductility and flexibility [9]. PAEs are difficult to dissolve in water and are easily soluble in organic solvents [10]. Thus, PAEs are difficult to combine with high polymer chains, but they can be combined with these chains through van der Waals forces or hydrogen bonding, which facilitates them to break from plastic products and enter atmosphere, water, soil, and other environmental media [11,12]. Then, these PAEs can accumulate in organisms through various means such as absorption, ingestion, respiration, and contact [13]. Most PAEs in the air emerge from the production and incineration of plastic products, spraying of coatings, and volatilization of plasticizers in agricultural plastic films; these PAEs have been detected in the atmosphere worldwide [14]. Although the natural complete degradation of plastic is difficult, it can be decomposed into microplastics through processes including sun exposure, weathering, and abrasion, which leads to several adverse effects, such as PAE accumulation and enrichment in organisms through food chains, which endangers organism health and the natural environment [15]. Dubaish et al. [16] reported the frequent appearance of suspended microplastics and black carbon particles in the coastal waters of the southern North Sea. Zhou [17] found that in the Bohai Sea, many benthic organisms, such as larval cockles, portunus snails, and red snails, exhibited microplastic accumulation. Zhao et al. [18] detected 39 volatile and semivolatile organic compounds in the water samples of Taihu Lake, and the content of PAEs was the highest. Wang et al. [19] detected a high concentration of PAEs in greenhouse soil covered with plastic film, and PAEs could accumulate in vegetables grown in this soil.
In recent years, with the increase in awareness regarding the pollution from environmental endocrine disruptors, environmental and health problems resulting from phthalate substances have become highly prominent. Countries such as China, the European Union, Japan, and the United States have all listed PAEs as environmentally superior pollutants [20]. Studies have shown that PAEs can enter organisms through various means, such as air, food, and body surface contact. Moreover, through the food chain, PAEs accumulate in the body, thereby producing multiple toxic effects, such as estrogenic/antiestrogenic effects and effects on androids/antiandrogens and thyroid hormone/antithyroid hormone [21]. Saillenfait [22] reported that exposure to diisooctyl phthalate during pregnancy can cause various reproductive system abnormalities, such as hypospadias and insufficient testis production, and can show an antiandrogen activity, thus disrupting normal male reproductive development. Shi et al. [23] investigated the urine samples of children aged 7-14 years and found a significant correlation between phthalate metabolite concentrations and puberty. Wang et al. [24] found that phthalate compounds can be detected in food containers, and their food simulant migration test proved that when PAEs migrate towards food, their content exceeds the standard limit, which causes toxic effects on the human body.
In view of the consistent performance of the insulation, bioconcentration, and toxicity of the PAEs in natural environments and organisms, the exploration of environmentally friendly phthalate derivatives with high insulation, low bioconcentration, and low toxicity is critical. This paper first normalizes the activity data of the insulation, bioconcentration, and toxicity of PAE molecules, uses the improved comprehensive index method to obtain the comprehensive effect value of PAEs with three properties, and uses this effect value as the data source to construct the 3D-QSAR comprehensive effect model of PAE molecules. Then, using DAP as the target molecule for molecular modification, 30 PAE derivative molecules with enhanced insulation and decreased toxicity, and bioconcentration are designed and screened. Finally, molecular dynamics is used to analyze the differences in the activity and mechanism of PAEs and their derivatives.

Data Source
In this paper, the comprehensive model and single model for PAEs were constructed by using the insulativity, toxicity, and bioconcentration parameters of 13 PAEs as indexes. Among them, the insulativity data of PAEs, which was expressed by their permittivity, was obtained by molecular dynamics. The higher the permittivity, the stronger the insulation [25]. The toxicity data of PAEs were expressed by median lethal concentration of 50% (LC 50 ) of PAEs to fish, cited from Estimation Program Interface (EPI) Suite database [26].
The lower the LC 50 value, the stronger the PAE toxicity. The bioconcentration data (bioconcentration factor, log BCF) of PAEs were predicted using the EPI database [27]. The greater the log BCF value, the stronger the PAE enrichment. Table 1 shows each parameter. The comprehensive index method is a comprehensive evaluation method used to transform multiple variable indicators into a comprehensive index to reflect the overall situation of evaluation subjects [28]. This method can be used not only to reflect the direction and degree of the overall change in complex economic phenomena but also to quantitatively explain actual economic effects caused by the change in phenomena. Especially in a complex model of multiple indicators, this method can consider not only the impact of indicators but also the interrelationship among them [29]. The relevant operation steps are as follows: first, normalization is performed to eliminate the dimensional coefficient between each index [30]; then, the mean method is used to determine the weight for each index; finally, the comprehensive index method is employed to integrate the given weights of indicators into the comprehensive index.
In this study, three property parameters of PAEs were arranged according to the numerical size to separately screen the maximum and minimum values of each indicator. According to the properties of indicator parameters, the positive or negative index formula was selected to normalize the parameters for dimension elimination. Then, the weight of each index was determined as 25:23:52 by using the mean method. Finally, the comprehensive index method was used to process the three weighted indicator parameters to integrate these parameters into the comprehensive effect values of 13 PAEs, which can effectively contain the PAE insulativity, toxicity, and bioconcentration information according to a certain transformation relationship. The specific calculation steps and formulas are as follows [30]: (1) Normalization

For the positive indicator of PAEs
For the negative indicator of PAEs : (2) Comprehensive index method

Construction of Comprehensive Model for PAEs' Insulativity, Toxicity, and Bioconcentration
Using the Minimize module in Sybyl-X 2.0 software, the three-dimensional structures of PAEs were optimized using the Powell conjugate gradient method to obtain the lowenergy conformation of each molecule, with selected Gasterger-Hückle charge and Tripos molecular force field. The energy convergence standard and iteration number were set to 0.005 kcal mol −1 and 10,000, respectively [31].
The PAE molecule with a high activity was selected as the template, the Align Database module was used to superpose the skeletons of multiple PAE compounds, and the molecular parameters of CoMFA and CoMSIA fields were calculated. Model construction involved the following steps: First, the molecular structures of PAEs available in the training set were imported into the training table, and CoMSIA field parameters were calculated using the Autofill module. Among the CoMSIA fields, three-dimensional hydrophobic, electrostatic, hydrogen bond donor, and hydrogen bond acceptor fields were selected as the molecular field. Relative permittivity was selected by referring to the distance, and the threshold was set to 125.4 kJ/mol. The remaining parameters followed the default system settings. Subsequently, when partial least squares method was used to analyze the conformational relationship between molecular structure and biological activity [32], leave-one-out method was selected for cross-validation of training set compounds [33] to calculate the cross-validation coefficient (q 2 > 0.5) and the optimal principal component number (n). No validation regression analysis was performed to obtain the non-cross validation coefficients (r 2 > 0.6), standard error of estimate (SEE), and test value (F) [34]. Finally, the test set constituting the remaining molecules was used to verify the constructed model. The constructed model exhibited high stability and universality when the calculated cross-validation coefficients (r 2 pred > 0.6) and standard error of predict (SEP) met the conditions [35,36].
In this paper, 10 PAE molecules were randomly selected from 13 PAEs (including template molecules) with comprehensive effect values as training sets in a ratio of about 3:1, and a 3D-QSAR comprehensive model with PAEs' insulation, toxicity, and bioconcentration was constructed. The remaining 4 PAEs (including template molecules) were used as the test set for the internal verification of the constructed comprehensive model.

Calculation of the Comprehensive Effect Value of PAEs with Insulation, Toxicity, and Bioconcentration
The insulation, toxicity, and bioconcentration data of PAE molecules processed using the normalization method were weighted (25:23:52), and then, the comprehensive effect value of 13 PAEs was calculated by employing the synthetic index method ( Table 2). The data presented in Table 2 were used to construct the CoMSIA model of PAEs' comprehensive effect with insulation, toxicity, and bioconcentration. Then, the threedimensional contour maps were used to analyze the proportion of the target molecule (DAP) in each force field. The model evaluation parameters are presented in Table 3. The optimal principal component (n) of the comprehensive effect CoMSIA model is 6, and the cross-validation coefficient (q 2 ) is 0.747 (>0.5) ( Table 3), indicating that the prediction ability of the developed model has a certain degree of credibility. The non-cross validation coefficient (R 2 ) is 0.929 (>0.9), the standard deviation (SEE) is 0.100, and the F value is 26.209, indicating that the model exhibits good fitting ability and robustness [37,38]. The model has certain reliability and universality in the prediction of the PAEs' comprehensive effect value with insulation, toxicity, and bioconcentration. Simultaneously, the CoMSIA model was used to observe the contribution rate of each force field. The contribution rates of DAP to the steric, electrostatic, hydrophobic, hydrogen bond donor, and hydrogen bond acceptor fields are 31.5%, 13.3%, 30.9%, 0.00%, and 12.3%, respectively. This finding showed that the comprehensive effect CoMSIA model can have a relatively large effect on the steric space, electrical distribution, hydrophobicity, and hydrogen bond acceptor of PAE molecules, while the influence of the hydrogen bond donor field on PAE molecules can be ignored.

Analysis of Three-Dimensional Contour Maps of PAEs' Comprehensive Effect Model
The DAP molecule was selected as the target molecule, the 3D-QSAR model of the combined effect with insulation, toxicity, and bioconcentration was imported, and the three-dimensional contour maps were analyzed by observing and comparing the color blocks of the DAP molecule in different force fields to obtain the factors affecting the DAP molecule activity under different force fields [39]. The three-dimensional contour map can visually and effectively show the interaction between small molecules and their receptors. It provides a certain basis for molecular modification, and the closer the molecule is to the color patch area of the contour map, the more beneficial the enhancement of molecular activity. There are four force fields in the CoMSIA model that can have different effects on the molecule activity in the field (Figure 1). In the steric field, the green region indicates that the introduction of bulky groups in this area can enhance the overall molecule activity. In the electrostatic field, the blue and red regions indicate that the introduction of positively and negatively charged groups, respectively, in this area improves the molecule activity. In the hydrophobic field, the yellow region indicates that the introduction of hydrophobic groups enhances the compound activity. In the hydrogen bond acceptor field, the purple-red region indicates that the increase in the hydrogen bond acceptor is conducive to molecule activity enhancement. Combined with the contribution rate of molecular field and the three-dimensional contour maps, this study ignored the influence of electrostatic field, hydrogen bond, and hydrogen receptor field when there was molecular modification of the target molecule, and only considered the influence of steric field and hy-drophobic field with the highest contribution rate in the molecular force field. In Figure 1A, it can be observed that the H atom attached to C-2 and the vicinity of the C-2 atom of the DAP molecule are all green, and the yellow area in Figure 1C basically covers the side chain on the benzene ring of the DAP molecule. This indicates that the introduction of bulky or hydrophobic groups at the sites covered by the side chains of DAP molecules can improve the PAEs' comprehensive effect values of insulation, toxicity, and bioconcentration.

Modification and Evaluation of DAP Derivatives Based on the CoMSIA Model
Comprehensively considering the three-dimensional equipotential map information and the hydrophobic fields, the modification sites shown in Figure 2 were selected to replace and modify the DAP molecule. According to the molecular modification site map, 13 groups with large volume and hydrophobic characteristics, including -F, -CL, Br, -SH, -CH 3 , -CH 2 CH 3 , -CH 2 CH 2 CH 3 , -C 6 H 5 , -CH 2 C 6 H 5 , -NO 2 , -CH 2 NO 2 , and -OCH 3 , were selected to modify the DAP at single site and dual site. A total of 56 DAP derivative molecules were designed, among which 27 exhibited single-site modification and 29 showed dual-site modification. The comprehensive effect and single effect of 56 designed PAE derivatives were predicted using the constructed PAEs' comprehensive effect model with insulation, toxicity, and bioconcentration and the single-effect models of insulation, toxicity, and bioconcentration, respectively. Consequently, a total of 30 DAP derivative molecules with higher comprehensive effect of insulation, toxicity, and bioconcentration than that of the target molecule (DAP), no reduction in insulation, and decreased toxicity and bioconcentration were screened (Table 4). This result indicated that the designed comprehensive effect model of PAE derivatives can effectively combine the molecular information of insulation, toxicity, and bioconcentration and provide a reliable basis for the modification of PAEs. The specific prediction results are presented in Table 4.  Discovery Studio (DS) software is used for molecular modelling and simulation and presents a very wide range of applications in molecular design. The scoring function value for DS can be used to characterize the binding strength between the ligand and receptor [40]. Therefore, insulative, toxic, and bioenriching receptors were selected, and the pharmacophore module of DS was used for the molecular docking of the target molecule and 30 PAE derivatives. The insulative receptor was derived from the carboxylic acid esters of thermophilicbacterium bacillus [41] (PDB ID: 5UOH), the toxicity receptor was deoxytocin synthase (PDB ID: 6XXK), and the bioconcentration receptor was obtained from cytochrome P-450 [42] (PDB ID: 4R20). The scoring function value shows that among the 30 PAE derivatives, 4 PAE derivatives have lower toxicity and biological concentration than that of the target molecule without reducing the insulation. Compared with the DAP molecule, the insulation, toxicity, and bioconcentration of DAP-1-NO 2 -2-C 6 H 5 increased by 21.23%, decreased by 3.73%, and decreased by 10.72%, respectively. This finding showed that the designed derivative molecules have good parameter performance (high insulation, low toxicity, and low bioconcentration), and the PAEs' comprehensive effect model can simultaneously consider the information of insulation, toxicity, and bioconcentration. The specific analysis results are shown in Table 5.

Mechanism Analysis of Insulation Improvement of PAE Derivatives Based on Molecular Dynamics
In order to further explore the related mechanism and effect of the improvement of the performance of PAE derivative molecules compared with target molecules, in this paper, by means of molecular dynamics simulation, carboxylesterase [43] (PDB ID: 5UOH) was used as PAEs' degradation enzyme, and DAP and its derivatives were respectively docked using SYBYL-2.0 X software to obtain the corresponding protein ligand complex. Then, GROMACs molecular dynamics software and simulation method of g-MMPBSA were used to simulate the binding free energy (∆G) of DAP and its derivatives with degradation enzyme under different applied electric field conditions (Table 6) [44]. Among them, the intermolecular binding free energy can quantitatively indicate the strength of the intermolecular interaction to a certain extent [45]. The negative value indicates that the receptor and ligand can effectively bind and produce corresponding effects, and the smaller the ∆G of PAEs and degrading enzymes, the stronger the degradation of PAEs by the enzymes. Studies have found that the applied electric field can promote cycloalkane degradation [46]. Tiehm et al. [47] reported that the applied electric field can stimulate the ability of noncharged compounds to degrade pollutants.
Under the conditions of no applied electric field (0 V) and applied electric field (5 V), the binding free energies of DAP, its three derivatives, and degradation enzyme receptor (5UOH) are <0, which indicates that this enzyme has a certain degradation effect on PAEs and their derivatives. In addition to DAP-1-NO 2 -2-CH 2 C 6 H 5 , the binding free energy values of the remaining two derivatives (DAP-1-NO 2 -2-CH 2 CH 2 CH 3 and DAP-1-NO 2 -2-OCH 3 ) decreased to different degrees (∆G 0V = −14.71%, 70.10%; ∆G 5V = −110.79%, −290.82%) compared with the target molecule (DAP), indicating that the designed deriva-tives have better stability and environmental friendliness. Moreover, the comparison of the free binding energy changes of the target molecule and the derivative molecules of PAEs in the absence and presence of the applied electric field showed that the free binding energy increases with the applied electric field. The binding free energy of DAP with applied electric field was 44.932 kJ/mol higher than that without applied electric field. The binding free energy of DAP-1-NO 2 -2-CH 2 C 6 H 5 with the applied electric field increased by 69.630 kJ/mol compared to without the applied electric field. The results show that the degradation ability of 5UOH to the derivative molecule is reduced, that is, the insulation of DAP-1-NO 2 -2-CH 2 C 6 H 5 molecule is enhanced, and this molecule can be used as the preferred derivative molecule with enhanced PAE insulation.

Conclusions
The 3D-QSAR model of PAEs' comprehensive effects with insulation, toxicity, and bioconcentration was constructed by using the comprehensive index method for the first time. A total of 30 PAE derivatives with enhanced insulation and reduced toxicity and bioconcentration were designed and screened by using the three-dimensional equipotential map of PAEs' multi-effect model and the models of comprehensive effect and single effect of insulation, toxicity, and bioconcentration. Four PAE derivative molecules with high insulation, low toxicity, and low bioconcentration were screened by comparing the scoring function values of insulation, toxicity, and bioconcentration of PAEs and their derivatives before and after modification. Analysis of molecular dynamics simulations showed that the increase in the free binding energy of DAP-1-NO 2 -2-CH 2 C 6 H 5 molecule was higher than that in the DAP molecule when the electric field was applied, which can explain the mechanism of enhanced insulation of the designed PAE derivative molecule. This paper provides a theoretical basis for the design of novel PAE plasticizers with high insulation, low toxicity, and low bioaccumulation.