Enhanced Biodegradation of Phthalic Acid Esters’ Derivatives by Plasticizer-Degrading Bacteria (Burkholderia cepacia, Archaeoglobus fulgidus, Pseudomonas aeruginosa) Using a Correction 3D-QSAR Model

A phthalic acid ester’s (PAEs) comprehensive biodegradability three-dimensional structure-activity relationship (3D-QSAR) model was established, to design environmentally friendly PAE derivatives, which could be simultaneously degraded by plasticizer-degrading bacteria, such as Burkholderia cepacia, Archaeoglobus fulgidus, and Pseudomonas aeruginosa. Only three derivatives of diethyl phthalate (DEP (DEP-27, DEP-28 and DEP-29)) were suited for their functionality and environmental friendliness, which had an improved stability in the environment and improved the characteristics (bio-toxicity, bioaccumulation, persistence, and long-range migration) of the persistent organic pollutants (POPs). The simulation inference of the microbial degradation path before and after DEP modification and the calculation of the reaction energy barrier exhibited the energy barrier for degradation being reduced after DEP modification and was consistent with the increased ratio of comprehensive biodegradability. This confirmed the effectiveness of the comparative molecular similarity index analysis (CoMSIA) model of the PAE’s comprehensive biodegradability. In addition, a molecular dynamics simulation revealed that the binding of the DEP-29 derivative with the three plasticizer-degradation enzymes increased significantly. DEP-29 could be used as a methyl phthalate derivative that synergistically degrades with microplastics, providing directional selection and theoretical designing for plasticizer replacement.


Introduction
In industrial production, a plasticizer is an indispensable part of microplastics and increases the flexibility and durability of plastic products. The phthalic acid esters (PAEs) are the most widely used plasticizers [1]. About 8.4 million tons of PAE plasticizers are used annually, accounting for 70% of its total use worldwide. Among these, diethyl phthalate (DEP) accounts for a high proportion [2]. The extensive use of phthalate plasticizers creates great commercial value, though it also poses environmental health risks that we cannot ignore. PAE plasticizers are not directly connected to plastic polymers, and they easily release into the environment during use [3]. Agricultural plastic film is one of the primary sources of microplastics in the soil. The PAE content in a vegetable greenhouse base could reach up to 9.68 mg/kg [4]. The PAEs in microplastics belong to refractory organic compounds, which are detected in soil [5], rivers [6], drinking water [7], food [8], household garbage [9], sewage sludge [10], 2 of 16 industrial wastewater [11], marine sediment [12] and landfill leach [13]. Toxicological studies proved that PAE plasticizers remaining in the environment could enter human and animal bodies through inhalation, diet and skin contact, causing great harm to human health and environmental safety [14]. Although PAEs are not acutely toxic to organisms, exposure to large doses could lead to teratogenic, carcinogenic and mutagenic effects in animals [15]. In addition, animal experiments confirmed that some phthalates cause liver and kidney damage in animals [16]. The US Environmental Protection Agency, the European Union, and the China National Environmental Monitoring Center have listed PAEs as the priority pollutants [17]. The degradation and transformation of phthalate plasticizers in environmental residual has become a topic of interest in recent years.
Microbial degradation is a key process that mineralizes the organic pollutants in the environment. Hoellein et al. [18] studied the biological effects of microplastics in the water surrounding sewage treatment plants. They observed that only specific microorganisms could affect migration behavior and the biodegradability of plasticizer. During microbial degradation, the PAEs could simultaneously be degraded by a variety of plasticizer-degrading bacteria. A study reported that 50% of the PAEs degraded after an inoculation of aerobic microorganisms in sewage sludge for 28 days [19]. Sugatt et al. [20] studied the biodegradability of 14 PAE plasticizers commercially used in microplastics using a shake flask test. They detected the PAE plasticizers were sensitive to the mixed microbial community in the natural environment, and they gradually degraded and mineralized. To improve the microbiological degradation efficiency of plasticizers, several scientists attempted to isolate a single species of bacteria with a high-efficiency of degradation from nature. Wu et al. [21] screened and isolated Ochrobactrum lupini and Agrobacterium tumefaciens that used dibutyl phthalate (DBP) as their only carbon source and energy source from polluted river sludge and completely mineralized it. The degrading enzymes (esterase, carboxylesterase and hydratase) play a crucial role in the synergistic degradation of PAE plasticizers in microplastics [22]. Under different environmental conditions, microorganisms degrade the PAE plasticizers, but the degradation cycle is relatively long [23]. Therefore, it is of immense significance to improve the biodegradability of the PAE plasticizers to inhibit their adverse effects on the environment and human health.
In this paper, we employed the range normalization method combined with the entropy weight method to characterize the comprehensive biodegradability of PAEs by enzymes from three plasticizer-degrading bacteria [24]. We constructed a three-dimensional structure-activity relationship (3D-QSAR) model of the comprehensive biodegradability of the PAEs. DEP was taken as the target molecule for molecular modification, and we screened for DEP derivatives suitable for degradation by plasticizer-degrading bacteria. Further, we evaluated the DEP derivatives for their functionality and environmental friendliness, providing directional selection for plasticizer replacement.

Data Source
The Burkholderia, Archaeoglobus [22] and Pseudomonas [25] species are typical plasticizer-degrading bacteria that effectively degrade PAE plasticizers. We obtained the structures of phthalate dioxygenase reductase from Burkholderia cepacia (PDB ID: 2PIA), esterase from Archaeoglobus fulgidus (PDB ID: 2ZYI) and carboxylesterase from Pseudomonas aeruginosa (PDB ID: 3CN7) from the Protein Data Bank (PDB (http://www.rcsb.org/pdb)) [26]. Using the molecular docking module in Sybyl-x2.0 software, we docked 17 different PAEs with the three degrading enzymes to obtain the corresponding scoring function value [27], which formed the source of the degradability data of the PAEs.

Construction of the Comparative Molecular Similarity Index Analysis (CoMSIA) Model of the Comprehensive Biodegradability of PAEs
We employed the range normalization method to standardize the degradation data of the three plasticizer-degrading bacteria and the entropy weight method to objectively weigh the standardized data [28]. Finally, we obtained the value of the comprehensive biodegradability of the PAE molecules, degraded by the three bacteria. The following are the range normalization formulas of the scoring function values of 17 PAE molecules after docking with the three plasticizer-degrading enzymes: Equation (1) is the processing formula for the positive index, where a higher index is desirable. Equation (2) is the processing formula for the negative index, where a lower index is desirable. Among them, i represents the different PAEs (i = 1, 2, 3, . . . , 17); j represents the degradability index of the different degradative enzymes (j = 1, 2, 3; corresponding to the enzymes from B. cepacia, P. aeruginosa and A. fulgidus, respectively); Y ij is the standardized value of each scoring function; X ij is the scoring function value of the j-th biodegradability index of the i-th PAE; X maxj and X minj are the maximum and minimum values of the corresponding indexes, a ∈ (0,1), which are generally 0.9. The scoring function value characterizes the index of microbial degradability with a higher score being desirable. Therefore, we used the processing formula for the positive index to standardize the microbial degradability data.
The entropy weight method is an objective method for weight processing [29], which uses information entropy for processing the weight. The calculation steps are as follows: Step 1: calculate the information entropy (E j ) of the j-th biodegradability index of the i-th PAE using the following formula: where k = 1/lnm, m is the total number of PAE molecules (i.e., 17), and E j is the information entropy of the j-th biodegradability index.
Step 2: calculate the difference coefficient (H j ) of the j-th biodegradability index, reflecting the difference in the degree of each index using the following formula: where H j is the difference coefficient of the j-th biodegradability index.
Step 3: normalize the index difference coefficient (H j ) to estimate the weight (W j ) of each biodegradability index using the following formula: where W j is the weight of the j-th biodegradability index, and n is the total number of biodegradability indicators (i.e., 3).
Step 4: calculate the comprehensive value of the biodegradability of the three plasticizer-degrading bacteria for each PAE using the following formula: where Z in is the comprehensive value of the biodegradability of the three plasticizer-degrading bacteria of the i-th PAE, and n is the total number of biodegradability indicators (i.e., 3). The molecular structure of the PAEs and the comprehensive biodegradability CoMSIA model of the three plasticizer-degrading bacteria were mainly constructed by using the 3D-QSAR module in The molecular structure of the PAEs and the comprehensive biodegradability CoMSIA model of the three plasticizer-degrading bacteria were mainly constructed by using the 3D-QSAR module in Sybyl-x 2.0 software (Tripos, Princeton, NJ, USA). The minimize module in the software optimized the molecular structure of the PAEs [30].To obtain a stable conformation of the lowest energy state of the molecule, the Powell's conjugate gradient method was employed and combined with the Tripos molecular force field. The latter sets the charge of the molecule to the Gasterger-Hückle charge [31]. The energy convergence standard was adjusted to 0.005 kcal mol −1 , and the number of iterations was set to 10,000 times [32]. In this paper, the diundecyl phthalate (DUP) with the highest comprehensive value of the biodegradability index of the three plasticizer-degrading bacteria was used as the template molecule ( Figure 1). The Align Database module of the software was selected to select the common part of the optimized PAE molecular structure as the common skeleton for superposition [33]. Sixteen molecules were randomly selected from the comprehensive biodegradation values of 17 PAEs molecules and three plasticizer degrading bacteria. . According to the ratio of 4:1, 13 PAE molecules were used as the training set, and the remaining three PAE molecules were used as the test set. The template molecule DUP was used in both the test set and the training set. The least squares method (PLS) was used to analyze the training set compounds. The training set was cross-validated under the extraction method module to obtain the cross-validation coefficient (q 2 ) and the optimal number of principal components (n). The training set compounds were analyzed using the no validation regression to obtain the cross-validation coefficient (R 2 ), standard deviation (SEE), test value (F) and the contribution rate of the force field (stereo field, electrostatic field). The scrambling stability test assessed the robustness of the model by evaluating the scrambling stability test parameters (Q 2 ), cross-validated standard error of prediction (cSDEP), and dQ 2 /dr 2 yy. The cross-validation method assessed the external prediction ability of the built model to the test set compounds with the external verification parameter (r 2 pred). Based on the information prompted by the contour maps of the constructed CoMSIA model, the substitution sites and substitution groups that had a greater impact on the comprehensive score value were screened, forming the theoretical basis for the subsequent molecular modification. In this paper, DEP was used as the target molecule for modification, which is the most widely used and frequently detected in the environment. Figure  2 displays the molecular structure of the target molecule DEP.  Sixteen molecules were randomly selected from the comprehensive biodegradation values of 17 PAEs molecules and three plasticizer degrading bacteria. According to the ratio of 4:1, 13 PAE molecules were used as the training set, and the remaining three PAE molecules were used as the test set. The template molecule DUP was used in both the test set and the training set. The least squares method (PLS) was used to analyze the training set compounds. The training set was cross-validated under the extraction method module to obtain the cross-validation coefficient (q 2 ) and the optimal number of principal components (n). The training set compounds were analyzed using the no validation regression to obtain the cross-validation coefficient (R 2 ), standard deviation (SEE), test value (F) and the contribution rate of the force field (stereo field, electrostatic field). The scrambling stability test assessed the robustness of the model by evaluating the scrambling stability test parameters (Q 2 ), cross-validated standard error of prediction (cSDEP), and dQ 2 /dr 2 yy. The cross-validation method assessed the external prediction ability of the built model to the test set compounds with the external verification parameter (r 2 pred ). Based on the information prompted by the contour maps of the constructed CoMSIA model, the substitution sites and substitution groups that had a greater impact on the comprehensive score value were screened, forming the theoretical basis for the subsequent molecular modification. In this paper, DEP was used as the target molecule for modification, which is the most widely used and frequently detected in the environment. Figure 2 displays the molecular structure of the target molecule DEP. The molecular structure of the PAEs and the comprehensive biodegradability CoMSIA model of the three plasticizer-degrading bacteria were mainly constructed by using the 3D-QSAR module in Sybyl-x 2.0 software (Tripos, Princeton, NJ, USA). The minimize module in the software optimized the molecular structure of the PAEs [30].To obtain a stable conformation of the lowest energy state of the molecule, the Powell's conjugate gradient method was employed and combined with the Tripos molecular force field. The latter sets the charge of the molecule to the Gasterger-Hückle charge [31]. The energy convergence standard was adjusted to 0.005 kcal mol −1 , and the number of iterations was set to 10,000 times [32]. In this paper, the diundecyl phthalate (DUP) with the highest comprehensive value of the biodegradability index of the three plasticizer-degrading bacteria was used as the template molecule ( Figure 1). The Align Database module of the software was selected to select the common part of the optimized PAE molecular structure as the common skeleton for superposition [33]. Sixteen molecules were randomly selected from the comprehensive biodegradation values of 17 PAEs molecules and three plasticizer degrading bacteria. . According to the ratio of 4:1, 13 PAE molecules were used as the training set, and the remaining three PAE molecules were used as the test set. The template molecule DUP was used in both the test set and the training set. The least squares method (PLS) was used to analyze the training set compounds. The training set was cross-validated under the extraction method module to obtain the cross-validation coefficient (q 2 ) and the optimal number of principal components (n). The training set compounds were analyzed using the no validation regression to obtain the cross-validation coefficient (R 2 ), standard deviation (SEE), test value (F) and the contribution rate of the force field (stereo field, electrostatic field). The scrambling stability test assessed the robustness of the model by evaluating the scrambling stability test parameters (Q 2 ), cross-validated standard error of prediction (cSDEP), and dQ 2 /dr 2 yy. The cross-validation method assessed the external prediction ability of the built model to the test set compounds with the external verification parameter (r 2 pred). Based on the information prompted by the contour maps of the constructed CoMSIA model, the substitution sites and substitution groups that had a greater impact on the comprehensive score value were screened, forming the theoretical basis for the subsequent molecular modification. In this paper, DEP was used as the target molecule for modification, which is the most widely used and frequently detected in the environment. Figure  2 displays the molecular structure of the target molecule DEP.

Evaluation of Functionality and Environmental Friendliness of DEP Derivatives Based on Density Functional Theory (DFT)
Stability was used as an evaluation index for the molecular functionality of the DEP derivatives [34] with the frequency, total energy and energy gap as the parameters. The above parameters were calculated by Gaussian 2009 software (Gaussian Inc., Wallingford, CT, USA) based on density functional theory (DFT) at the unit level of b3pw91/6-31G* [35]; environmental friendliness was evaluated by the four persistent organic pollutants' (POPs) characteristics of toxicity (LC 50 ), persistence (logt 1/2 ), bioconcentration (logBCF) and migration (−logPL) [36].  Table 1. DUP exhibited the highest comprehensive biodegradability among all the PAEs. It was the template molecule for the common skeleton superpositioning of the training set compounds to construct a CoMSIA model of the comprehensive biodegradability of the PAEs by the three plasticizer-degrading bacteria. In Table 1, a and b represent the training set and the test set of the model, respectively. Table 2 lists the evaluation parameters of the CoMSIA model. The cross-validation coefficient (q 2 ), calculated by the cross-validation analysis, was 0.731 (>0.5), and the principal component number (n) of the model was 8. Both of these indicate that the CoMSIA model had good internal prediction capabilities [37]. The standard deviation (SEE) calculated by the non-cross-validation analysis was 0.013 (<0.95), the test value (F) was 578.738, and the non-cross-validation coefficient (R 2 ) was 0.999 (>0.9). These indicate that the model had a high fitting ability [38]. The perturbation stability test parameter (Q 2 ) was 0.85, cSDEP was 0.251, and dq 2 /dr 2 yy was 1.488, indicating that the model had a high stability [39]. The external validation coefficient (r 2 pred ), estimated by the external validation analysis of the test set compound, was 0.761 (>0. 6), indicating that the model exhibited a good external prediction ability [40]. In addition, the contribution rate of each force field was 34.7% in the steric field (S), 11.1% in the electrostatic field (E), 46.8% in the hydrophobic field (H), 7.4% in the hydrogen bond acceptor field (D) and 0.0% in the hydrogen bond donor field (A). The contribution rates of the steric, electrostatic field and hydrophobic field were the highest. Therefore, the force field information of the electrostatic and hydrophobic fields in the contour maps of the comprehensive biodegradability of the DEP molecule was analyzed ( Figure 3). The steric field information map (Figure 3a) displays that introducing small-volume groups into the yellow region could improve the comprehensive biodegradation value of the DEP molecules by the three plasticizer-degrading bacteria. In other words, the groups with volumes less than that of -CH 2 CH 3 were introduced at the C 1 position. In the electrostatic field (Figure 3b), introducing positively charged groups into the blue region or negatively charged groups into the red region could effectively improve the comprehensive biodegradability of the DEP molecules [41]. However, the red region was located in the common skeleton and was difficult to replace and modify. Therefore, a group with a positive charge greater than -H introduced into the C 1 position could achieve the purpose of increasing the comprehensive biodegradability. In the hydrophobic field (Figure 3c), introducing strong hydrophilic substituents in the white area at the C 2 position was conducive to increase the comprehensive biodegradability of the DEP molecules. In summary, to improve the comprehensive biodegradability of DEP, it was necessary to modify its structure by introducing single and double substitutions: (i) at the C1 position, a group with a volume smaller than -CH 2 CH 3 and a group with an electropositivity greater than -H, and (ii) at the C 2 position, a hydrophilic group. had good internal prediction capabilities [37]. The standard deviation (SEE) calculated by the non-cross-validation analysis was 0.013 (<0.95), the test value (F) was 578.738, and the non-cross-validation coefficient (R 2 ) was 0.999 (>0.9). These indicate that the model had a high fitting ability [38]. The perturbation stability test parameter (Q 2 ) was 0.85, cSDEP was 0.251, and dq 2 /dr 2 yy was 1.488, indicating that the model had a high stability [39]. The external validation coefficient (r 2 pred), estimated by the external validation analysis of the test set compound, was 0.761 (>0. 6), indicating that the model exhibited a good external prediction ability [40]. In addition, the contribution rate of each force field was 34.7% in the steric field (S), 11.1% in the electrostatic field (E), 46.8% in the hydrophobic field (H), 7.4% in the hydrogen bond acceptor field (D) and 0.0% in the hydrogen bond donor field (A).

Contour Map Analysis of the CoMSIA Model
The contribution rates of the steric, electrostatic field and hydrophobic field were the highest. Therefore, the force field information of the electrostatic and hydrophobic fields in the contour maps of the comprehensive biodegradability of the DEP molecule was analyzed ( Figure 3). The steric field information map (Figure 3a) displays that introducing small-volume groups into the yellow region could improve the comprehensive biodegradation value of the DEP molecules by the three plasticizer-degrading bacteria. In other words, the groups with volumes less than that of -CH2CH3 were introduced at the C1 position. In the electrostatic field (Figure 3b), introducing positively charged groups into the blue region or negatively charged groups into the red region could effectively improve the comprehensive biodegradability of the DEP molecules [41]. However, the red region was located in the common skeleton and was difficult to replace and modify. Therefore, a group with a positive charge greater than -H introduced into the C1 position could achieve the purpose of increasing the comprehensive biodegradability. In the hydrophobic field (Figure 3c), introducing strong hydrophilic substituents in the white area at the C2 position was conducive to increase the comprehensive biodegradability of the DEP molecules. In summary, to improve the comprehensive biodegradability of DEP, it was necessary to modify its structure by introducing single and double substitutions: (i) at the C1 position, a group with a volume smaller than -CH2CH3 and a group with an electropositivity greater than -H, and (ii) at the C2 position, a hydrophilic group.

Molecular Modification and Prediction of Comprehensive Biodegradability of DEP
According to the force field information from the contour maps of the CoMSIA model, the groups was substituted at the C 1 and C 2 positions to conduct single and double substitutions on the DEP molecule. The groups with volumes smaller than -CH 2 CH 3 (-CH 3 , -OH, -H, -CN, -NH 2 , -CHO) and groups with an electropositivity greater than -H (-CH 3 , -CH 2 CH 3 , -CH(CH 3 ) 2 , -C(CH 3 ) 3 ) were introduced at the C 1 position. At the C 2 position, the groups with higher hydrophilicities (-OH, -CHO, -COOH, -NH 2 , -COCH 3 , -CONH 2 ) were introduced. These substitutions would improve the comprehensive biodegradability of derivative molecules. Thus, a total of 30 kinds of DEP derivatives were designed accordingly. Using the constructed CoMSIA model to predict the comprehensive biological properties of the modified derivatives, it was found that the comprehensive biodegradability of most of the DEP derivatives enhanced significantly (Table 3). Among them, the comprehensive biodegradability of 12 DEP derivatives increased by more than 15%, and the comprehensive biodegradation value of DEP-26 molecules had the highest growth rate, reaching 50.37%. DEP-23 molecules had the same substitution groups at site C 2 as in the DEP-26 molecules. However, the H atom at site C 1 (H 1 atom) in DEP-26 resulted in a higher comprehensive biodegradation value than the corresponding substituent in DEP-23. DEP-21 was introduced to more hydrophilic groups, increasing its comprehensive biodegradability value to more than that of DEP-15. In summary, the prediction of the comprehensive biodegradability of the DEP derivatives by the three plasticizer-degrading bacteria was consistent with the information presented in the contour maps of the CoMSIA model.

Verification of the CoMSIA Model for the Comprehensive Biodegradability of PAEs Molecules
In this paper, Sybyl-x2.0 software was also used to construct the 3D-QSAR model of the single biodegradability of the PAEs by the three plasticizer-degrading bacteria ( Table 4). The DEP derivative molecules were screened through the biodegradable 3D-QSAR model of the PAEs' molecular plasticizer-degrading bacteria, whose degradation efficiency ratio was close to the weight ratios of the entropy weight method in the comprehensive biodegradability model ( Table 3). The docking score of the DEP derivative with phthalate dioxygenase reductase was 4.939-6.725, with a growth rate of up to 20.74%. Compared with the growth rate before modification, the docking score values of the DEP derivative binding with esterase ranged from −2.46% to 27.71%. The scoring functions for carboxylesterase and DEP derivatives exhibited a significant increase, with the growth rate ranging from 15.86% to 24.96%. The ratio of the increase in the biodegradation effect value of DEP-23-DEP-30 by three plasticizer-degrading bacteria was the closest to that calculated and weighted by the entropy weight method. These results verified that the CoMSIA model effectively provides the biodegradability information of DEP derivatives when degraded by three plasticizer-degrading bacteria. Additionally, this result demonstrated that the CoMSIA model modified by the entropy weight method is reliable with a good predictive ability, and can be applied to the design modifications in PAE molecules.

Evaluation of the Functionality and Environmental Friendliness of DEP Derivatives
Molecular stability was the primary index for evaluating the functionality of DEP derivatives, characterized by the energy value, energy gap value, and positive frequency value. The persistent organic pollutant (POP) characteristics determined the environmental friendliness of the DEP derivatives, before and after modification, by evaluating their bio-toxicity (logLC 50 ), bioaccumulation (logBCF), persistence (logt 1/2 ) and long-range migration (logK OA ) [36]. The energy value was inversely proportional to the stability of the molecules in the environment [42]. A higher positive frequency value indicated that the DEP derivative molecules could exist in the environment [43]. The bio-toxicity and long-range migration of the DEP derivatives were inversely proportional to the predicted logLC 50 and logK OA values, respectively. On the contrary, the bioaccumulation and persistence of the DEP derivatives were directly proportional to the predicted logBCF and logt 1/2 values, respectively. Table 5 summarizes the predicted values of the parameters evaluating the functionality and environmental friendliness of the DEP derivatives, before and after modification. Table 5 displays that the energy values of the DEP derivatives were significantly lower than that of the DEP molecule, with a reduction of 19.91-37.39%. The energy gap values did not change significantly before and after modification, indicating that the modified DEP derivatives were highly stable. In addition, the positive frequency values of the DEP derivatives were greater than zero, indicating that the modified molecules could exist in the environment. Among the eight DEP derivatives, only DEP-27, DEP-28, and DEP-29 exhibited no significant changes in their bio-toxicity compared with the unmodified molecule. The bio-toxicity of the remaining five derivatives significantly increased, with the highest increase-rate going up to 38.91%. An organic compound with a logBCF value less than 100 does not easily accumulate in organisms and has a small impact on the environment [44]. Except for the DEP-25 and DEP-26 molecules, the logBCF values of the DEP derivatives were much lower than 100. Among them, the logBCF values of DEP-27, DEP-28, and DEP-29 after modification did not change significantly. The predicted half-lives of the DEP derivatives were slightly higher than that of DEP, indicating that the degradability of these derivatives improved marginally. The predicted logK OA values of the eight DEP derivatives were significantly higher than that of DEP, suggesting that their mobilities were decreasing. In conclusion, among the eight DEP derivatives, only DEP-27, DEP-28 and DEP-29 exhibited significantly higher biodegradation values, and their functionalities and environmental friendliness were also better than the unmodified DEP molecule.

Simulation of Microbial Degradation of DEP and Its Derivative Molecules
Plasticizer residue poses a serious threat to the environment. Among the possible ways of degradation, microbial degradation is the primary elimination process of PAE plasticizers. It achieves the purpose of recycling elements and balancing the ecosystem [45]. Gram-negative and Gram-positive bacteria degrade PAEs via different pathways, but both eventually form protocatechate [46]. Protocatechate can convert into pyruvate, succinate, and oxaloacetate, entering the tricarboxylic acid cycle and finally mineralizing into CO 2 and H 2 O [47]. Ren et al. [48] found that esterase plays a crucial role in the microbial degradation of PAEs while studying their microbial degradation mechanism. Ester bond hydrolysis forms the key initial step in the microbial degradation of PAEs [49]. Based on the microbial degradation path of PAEs, the screened DEP derivative molecules (DEP-27, DEP-28 and DEP-29) were taken as examples to simulate and derive the microbial degradation path before and after DEP molecular modification (Figure 4). A Gaussian calculation was carried out for the energy barrier of the reaction to compare the difficulty of the biodegradation process before and after the DEP molecular modification (Table 6). degradation mechanism. Ester bond hydrolysis forms the key initial step in the microbial degradation of PAEs [49]. Based on the microbial degradation path of PAEs, the screened DEP derivative molecules (DEP-27, DEP-28 and DEP-29) were taken as examples to simulate and derive the microbial degradation path before and after DEP molecular modification (Figure 4). A Gaussian calculation was carried out for the energy barrier of the reaction to compare the difficulty of the biodegradation process before and after the DEP molecular modification (Table 6).   As shown in Figure 4, the microbial degradation transformation paths of DEP and its derivatives are divided into the transformation paths of Gram-negative bacteria and Gram-positive bacteria. Under the action of the plasticizer-degrading enzymes, DEP first hydrolyzed to phthalate monoesters (M0-1) and then hydrolyzed to form phthalate (M0-2). Phthalate 4,5-dioxygenase in Gram-negative bacteria oxidized the hydroxyl phthalate to produce 4,5-dihydroxyphthalic acid (M0-3), which decarboxylated to protocatechate (M0-5) [50]. Gram-positive bacteria hydrolyzed phthalic acid at C 3 and C 4 to produce 3, 4-dihydroxyphthalic acid (M0-4), which decarboxylated to protocatechate (M0-5) [51]. Finally, protocatechate mineralized into CO 2 and H 2 O through the tricarboxylic acid cycle, culminating in the complete degradation of PAEs [47]. When the C 2 site of DEP was replaced with -CONH 2 , the resulting derivative was difficult to hydrolyze into phthalic acid [52]. So the degradation process was slightly different from that of DEP.

Calculation of the Reaction Energy Barrier for Microbial Degradation Transformation Paths of DEP and Its Derivatives
The reaction energy barrier required for the microbial degradation of the three DEP derivatives was lower than that required for the degradation of DEP (Table 6). A smaller reaction barrier is desirable as it improves the likelihood of the reaction path to occur. This indicates that the microorganisms could easily degrade the three DEP derivatives compared with DEP. The change rates of the reaction energy barrier for DEP-27, DEP-28, and DEP-29 were −20.15%, −23.42%, and −30.26%, respectively. These were consistent with the enhanced comprehensive biodegradability predicted by the comprehensive biodegradability model (DEP-27: 23.33%, DEP-28: 27.04% and DEP-29: 31.85%). This validates the reliability of the comprehensive biodegradation model, and indicates that the biodegradability of DEP can improve by molecular modification.

Simulation and Verification of the Molecular Dynamics of the Microbial Degradation of DEP and Its Derivatives
In this paper, GROMACS 4.6.5 software was used to simulate the molecular dynamics (MD) of the ligand complex structures of DEP and its derivatives in Sybyl-x2.0, docked with phthalate dioxygenase reductase, esterase and acid esterase. The Poisson-Boltzmann surface area (MM-PBSA) method was used to calculate the binding free energy in the molecular dynamics [53]. Compared with the scoring function of molecular docking, the degree of binding energy characterized the biodegradability of DEP molecules before and after modification (Table 7). Higher docking scores, indicate a stronger binding of the molecules to the degradative enzyme [54]. Additionally, the binding energy is directly proportional to the affinity between the molecule and the degrading enzyme [55], indicating a higher biodegradability of the molecule. Compared with the DEP molecules, DEP-27 and DEP-28 had lower scoring functions and higher free energy of binding to the degrading enzymes ( Table 7). The scoring functions of only DEP-29 binding to the three degrading enzymes increased significantly, and its binding free energy also decreased significantly. The results demonstrated that DEP-29 was the only derivative that could easily bind to the three degrading enzymes at the same time and that the degree of binding was significant. This not only proves the rationality of the comprehensive biodegradation model of PAEs to modify the derivatives, but also forms the means to screen the PAE derivatives that synergistically degrade with microplastics.

Conclusions
In this paper, a 3D-QSAR model of the comprehensive biodegradability of the phthalic acid esters (PAEs) was constructed by combining the range normalization and entropy weight methods. In combination with molecular modifications, this model was successfully applied to design environmentally friendly PAEs that can co-degrade with microplastics. These PAE derivatives significantly improve the biodegradability of PAE plasticizers in microplastics, and can reduce the residual PAEs in the natural environment, relieve the adverse effects of plasticizer residues in the human body and environment, and provide the directional selection and theoretical support for the replacement of plasticizers. Further research work will be focused on whether the designed PAE derivatives could be degraded synergistically with the plasticizer-degrading bacteria or not, and the synthesis and biodegradation of these derivatives would provide an effective verification for the molecular modification method available in this study.