Next Article in Journal
Ameliorative Effect of Ocimum forskolei Benth on Diabetic, Apoptotic, and Adipogenic Biomarkers of Diabetic Rats and 3T3-L1 Fibroblasts Assisted by In Silico Approach
Next Article in Special Issue
Influence of the Active Site Flexibility on the Efficiency of Substrate Activation in the Active Sites of Bi-Zinc Metallo-β-Lactamases
Previous Article in Journal
Enhanced Extraction Efficiency of Flavonoids from Pyrus ussuriensis Leaves with Deep Eutectic Solvents
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Molecules
Volume 27
Issue 9
10.3390/molecules27092799
molecules-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Computational Insight into Biotransformation Profiles of Organophosphorus Flame Retardants to Their Diester Metabolites by Cytochrome P450

by
Yue Jia
1,†,
Tingji Yao
1,2,†,
Guangcai Ma
1,*,
Qi Xu
1,
Xianglong Zhao
1,
Hui Ding
1,
Xiaoxuan Wei
1,
Haiying Yu
1,* and
Zhiguo Wang
3
1
College of Geography and Environmental Sciences, Zhejiang Normal University, Jinhua 321004, China
2
School of Environment, Hangzhou Institute for Advanced Study, University of Chinese Academy of Sciences, Hangzhou 310024, China
3
Institute of Ageing Research, School of Medicine, Hangzhou Normal University, Hangzhou 311121, China
*
Authors to whom correspondence should be addressed.
†
These authors contributed equally to this work.
Molecules 2022, 27(9), 2799; https://0-doi-org.brum.beds.ac.uk/10.3390/molecules27092799
Submission received: 13 April 2022 / Revised: 24 April 2022 / Accepted: 26 April 2022 / Published: 28 April 2022
(This article belongs to the Special Issue Molecular Modeling: Insights into the Enzymatic Reactions and Photochemical Processes in Biomacromolecules)
Browse Figures
Versions Notes

Abstract

:
Biotransformation of organophosphorus flame retardants (OPFRs) mediated by cytochrome P450 enzymes (CYPs) has a potential correlation with their toxicological effects on humans. In this work, we employed five typical OPFRs including tris(1,3-dichloro-2-propyl) phosphate (TDCIPP), tris(1-chloro-2-propyl) phosphate (TCIPP), tri(2-chloroethyl) phosphate (TCEP), triethyl phosphate (TEP), and 2-ethylhexyl diphenyl phosphate (EHDPHP), and performed density functional theory (DFT) calculations to clarify the CYP-catalyzed biotransformation of five OPFRs to their diester metabolites. The DFT results show that the reaction mechanism consists of Cα-hydroxylation and O-dealkylation steps, and the biotransformation activities of five OPFRs may follow the order of TCEP ≈ TEP ≈ EHDPHP > TCIPP > TDCIPP. We further performed molecular dynamics (MD) simulations to unravel the binding interactions of five OPFRs in the CYP3A4 isoform. Binding mode analyses demonstrate that CYP3A4-mediated metabolism of TDCIPP, TCIPP, TCEP, and TEP can produce the diester metabolites, while EHDPHP metabolism may generate para-hydroxyEHDPHP as the primary metabolite. Moreover, the EHDPHP and TDCIPP have higher binding potential to CYP3A4 than TCIPP, TCEP, and TEP. This work reports the biotransformation profiles and binding features of five OPFRs in CYP, which can provide meaningful clues for the further studies of the metabolic fates of OPFRs and toxicological effects associated with the relevant metabolites.

Graphical Abstract

1. Introduction

To improve the fire resistance of household products and fulfill the increasingly rigorous flammability standards, flame retardants (FRs) are often used as additives in various consumer supplies, such as furniture foam, infant products, textiles, and electronics [1,2,3]. The polybrominated diphenyl ethers (PBDEs) have been widely used in polyurethane foam for many years as the mainstream FRs [4]. Since the mid-2000s, PBDEs have been phased out in many countries due to their bioaccumulation, persistence, and toxicity to humans, and thereafter the interest in the alternatives of PBDEs such as organophosphate flame retardants (OPFRs) has been growing [5,6,7]. With the rapid increase in OPFRs production in recent years, global consumption reached 500,000 tons in 2011 and was estimated to be 680,000 tons in 2015 [8]. Because OPFRs are usually water-soluble and are not covalently bound to materials, they thus can easily discharge into the environment through wear, volatilization, and dissolution during use [9,10].
OPFRs have been widely detected in various environmental media, such as indoor and outdoor air, water, soil, dust, and sediments [3,11,12,13,14,15,16,17], and the levels of OPFRs in the indoor environment are significantly higher than those of brominated flame retardants (BFRs) [18]. Humans may be exposed to OPFRs through breathing, skin contact, and dietary intake, which has potential adverse effects on health [19]. Actually, OPFRs have been detected in human hair, breast milk, and urine samples [20,21,22,23]. Toxicological studies have shown that the long-term exposure of OPFRs to animals may lead to adverse effects on their reproductive and endocrine system [24,25,26,27]. In the exposure experiment using a rat brain sphere in vitro model, the enhancive expressions of cytokine gene and receptor demonstrate that OPFRs may induce an inflammatory response [28]. Moreover, OPFRs exposure may reduce the proliferation and growth of human neural stem cells, rat neuronal growth, and network activity [29].
Cytochrome P450 enzymes (CYPs) are phase-I metabolic enzymes involving primarily in the biotransformation of xenobiotic compounds in diverse organisms, which may substantially change the toxicological and physicochemical properties of OPFRs through metabolism [30]. The OPFRs can be metabolized by CYPs to produce hydroxylated metabolite and diester metabolite via the competitive C-hydroxylation and O-dealkylation reactions [31,32]. Although the structural differences among OPFRs can result in various mono-/di-hydroxylated metabolites, diester may be produced as the common metabolite. For instance, both human and Brevibacillus brevis CYPs can metabolize triphenyl phosphate (TPHP) into diphenyl phosphate (DPHP) as the important metabolite, and CYP1A2 and CYP2E1 isoforms are mainly involved in the metabolism in human liver microsomes (HLMs) [6,33]. Incubation of tris(1,3-dichloro-2-propyl) phosphate (TDCIPP) with HLMs results in the formation of bis(1,3-dichloro-2-propyl) phosphate (BDCIPP) and other metabolites of oxidative dehalogenation [32], and CYP3A4 shows the highest activity toward TDCIPP metabolism [34]. Indeed, 2-ethylhexyl diphenyl phosphate (EHDPHP) can be transformed by human CYPs into mono- and di-hydroxylated metabolites, keto metabolites, and diphenyl phosphate (DPHP) as major phase-I outputs [35]. Chen et al. further reported CYP3A4 to be the major activating enzyme involved in EHDPHP metabolism [36]. Hou et al. reported the in vitro metabolism kinetics of tris(2-butoxyethyl) phosphate (TBOEP) and tris(n-butyl) phosphate (TNBP) and identified CYP3A4 and CYP1A as the major CYP isoforms catalyzing the metabolism in fish liver microsomes [37]. However, to date there are few studies have explored the CYP-mediated biotransformation mechanism and reactivity of OPFRs in-depth, and binding interactions between OPFRs and specific CYP isoforms remain unclear, and deserve further research.
In the past two decades, molecular modeling techniques, such as molecular dynamics (MD) simulations and quantum chemical calculations, have been extensively used to shed light on the binding and metabolism of xenobiotics in various biomacromolecules including CYPs [38,39,40,41,42,43]. In this work, five OPFRs detected widely in the environment, including TDCIPP, tris(1-chloro-2-propyl) phosphate (TCIPP), tris(2-chloroethyl) phosphate (TCEP), triethyl phosphate (TEP), and EHDPHP (Figure S1 in the Supplementary Materials), were employed to investigate the reaction mechanism leading to their diester metabolites through the density functional theory (DFT) calculations. Further MD simulations and binding free energy calculations were performed to predict the binding features and affinities of five OPFRs in CYP3A4, and to confirm the DFT results.

2. Computational Procedures

2.1. Model System

The reactive iron(IV)−oxo porphyrin species Compound I (Cpd I) of CYP was mimicked by a simplified model, Fe4+O2−(C20N4H12)−1(HS)−1, to unravel the mechanistic details of biotransformation of 5 OPFRs. Numerous previous studies have proven the high reliability and effectiveness of this simple model to accurately reproduce the electronic structures and reactivity of Cpd I [44,45,46,47].

2.2. DFT Methodology

All DFT calculations were performed by Gaussian 09 program [48]. The geometry optimizations and frequency analyses were carried out using the spin-unrestricted UB3LYP functional [49,50,51], integrating the LANL2DZ pseudopotential basis set [52] for the Fe atom and 6-31G(d,p) basis set [53] for other atoms (denoted BS1). Frequency analyses were used to obtain the thermal corrections to the Gibbs free energy, and to confirm the nature of all optimized geometries, among which all reactant/intermediate structures showed only real frequencies and the transition states had a single imaginary frequency. Based on the optimized geometries, single-point energy (SPE) calculations were performed using a larger def2-TZVP basis set [54] for all atoms (BS2) to obtain more accurate energies. Dispersion corrections were introduced into SPE calculations using the DFT-D3(BJ) method [55] to further improve the accuracies of B3LYP energies. Furthermore, the SMD solvation model [56] with nonpolar chlorobenzene and polar water was used to implicitly simulate CYP active site and aqueous-phase surroundings at the UB3LYP-D3/BS2 level, respectively. The relative energies reported in this work are thus the single-point energies with the inclusion of solvation and dispersion effects and Gibbs free energy corrections. It should be noted that Cpd I shows two closely lying spin states, including high-spin (HS) quartet and low-spin (LS) doublet states [57,58], and thus the energy profiles of biotransformation routes of 5 OPFRs were evaluated in both HS and LS states.

2.3. Molecular Docking

We selected CYP3A4 as the potential receptor of OPFRs because this isoform has a higher abundance and larger active pocket than the other CYPs in human liver tissue [59,60,61]. Moreover, in vitro experiment has identified CYP3A4 as the major CYP isoform catalyzing the metabolism of TDCIPP and EHDPHP [34,36]. CYP3A4 may also be involved in the metabolism of TCIPP, TCEP, and TEP due to the structural similarity of these 3 OPFRs with TDCIPP.
On the basis of the crystal structure of human 3A4 (PDB code: 2v0m) obtained from the Protein Data Bank (https://www.rcsb.org/, accessed on 15 April 2021) [59,62], molecular docking simulations were performed using Autodock Vina program [63] to construct the initial 3A4-OPFR binding complexes. Before docking, the bound ligand and water molecules in 3A4 were manually removed, and then the missing residues were complemented using the Chimera program [64]. During the docking, CYP3A4 was kept as the rigid receptor, while OPFRs were set as the flexible ligand. The grid box with suitable three dimensions was set to accommodate the ligand and CYP active site. The obtained 5 CYP3A4-OPFR complexes with the lowest binding affinities were shown in Figure S2, which were selected as the initial conformations for the following MD simulations.

2.4. Molecular Dynamics Simulations

All MD simulations were carried out using Amber12 program [65] with the ff14SB force field [66]. In CYP3A4, all Asp and Glu residues were set to be deprotonated, while all Lys and Arg were set to be protonated. The protonation states of all His residues were determined based on the predicted pKa values by the PDB2PQR Server [67] and the visual inspections of surrounding hydrogen-bonding networks. Specially, His30, His54, His65, His287, and His402 were singly Nε-protonated, His324 was singly Nδ-protonated, and His267 was doubly Nδ- and Nε-protonated. The force field parameters of 5 OPFRs were generated using the Antechamber module of AmberTools. All CYP3A4-OPFR complexes were solvated using the truncated octahedral TIP3P water box [68], and then the generated solvation models were neutralized by adding counterions. Before the production simulation, each solvated system was pre-equilibrated through 4000 ps energy minimization, 500 ps heating from 0 K to 300 K, 500 ps density equilibration at 300 K, and then 1000 ps constant pressure equilibration at 300 K. The final production simulations with different time scales were performed to equilibrate the solvated systems.

2.5. Binding Free Energy Calculations

After the production simulations, 200 snapshots were extracted from the last 20 ns MD trajectories to calculate the binding free energies (ΔGbind) between CYP3A4 and OPFRs. Based on the molecular mechanics/generalized Born surface area (MM/GBSA) method [69], ΔGbind was calculated using the following Equations (1)–(4):
ΔGbind = Gcomplex − (GCYP3A4 + GOPFRs)
Gbind = EMM + Gsolv − TS
EMM = Eele + EvdW
Gsolv = GGB + GSA
Gcomplex, GCYP3A4, and GOPFRs refer to the free energies of the binding complex, CYP3A4, and OPFRs, respectively. EMM represents the molecular mechanics energy in the gas phase, consisting of electrostatic energy (Eele) and van der Waals interaction energy (EvdW). Gsolv refers to the solvation free energy, which can be divided into a polar term (GGB) and a nonpolar term (GSA). TS is the entropy contribution. The calculated ΔGbind can be used to compare the binding affinities of different OPFRs in CYP3A4. Furthermore, the active site residues contributing significantly to the ligand-binding were also identified by binding energy decomposition calculations (energy contribution < −1.0 kcal/mol), which is essential for understanding the molecular mechanism of ligand-receptor binding [70].

3. Results and Discussion

3.1. Biotransformation Profiles of OPFRs by CYP

Biotransformation of TDCIPP to BDCIPP. Previous studies have detected BDCIPP as the main product in the phase-I metabolism of TDCIPP [1,71,72], and CYPs may be involved in the metabolic process [32]. Our calculations confirm that the biotransformation of TDCIPP to BDCIPP proceeds through the Cpd I-mediated C-hydroxylation and intramolecular O-dealkylation processes. C-Hydroxylation begins with the hydrogen atom transfer (HAT) from Cα atom (denoted C1) of TDCIPP to FeIV=O unit of Cpd I, generating an unstable C1-radical intermediate (IM1-H) with reduced FeIV-OH species. C1-radical further captures OH from FeIV-OH to form a hydroxylated intermediate (IM1-OH) with resting FeIII-porphyrin species. Subsequently, IM1-OH undergoes O-dealkylation to produce the final diester metabolite BDCIPP (PBDCIPP) with 1,3-dichloroacetone.
The free energy profiles for biotransformation of TDCIPP and the optimized geometries of relevant reaction species in both HS and LS states are shown in Figure 1, and the structural characteristics of transition states of HAT and OH rebound are also summarized in Table S1. HAT proceeds through a transition state TS1-H with an LS/HS energy barrier of 27.89/27.52 kcal/mol relative to the reactant complex (RC). The formation of intermediate IM1-H is endothermic by 5.73 kcal/mol in the LS state, while it is exothermic by 3.05 kcal/mol in the HS state. The subsequent OH rebound to C1-radical is essentially a barrier-free process, although it undergoes a transition state 4TS1-OH with a tiny barrier of 0.36 kcal/mol in the HS state. The collapse of IM1-H to IM1-OH is strongly exothermic, which facilitates the follow-up O-dealkylation of IM1-OH. As shown in Figure 2b, the intramolecular O-dealkylation proceeds through H2O-assisted proton transfer from OH to P=O group, undergoing a transition state TSBDCIPP with a barrier of 8.26 kcal/mol to generate the exothermic PBDCIPP.
Biotransformation of TCIPP to BCIPP. TCIPP shows highly structural similarity with TDCIPP, and therefore they may have similar metabolic routes. Figure 2 shows the free energy profiles for biotransformation of TDCIPP to bis(1-chloro-2-propyl) phosphate (BCIPP). The LS and HS barriers of HAT are calculated to be 24.27 and 21.17 kcal/mol, respectively, which are lower than those in TDCIPP metabolism (27.89 and 27.52 kcal/mol). The barrier difference may be attributed to the chlorine substituent. Compared with TDCIPP, TCIPP has only one chlorine substituent in each alkyl side chain, which may reduce the steric hindrance effect and further enhance the reactivity of C1-hydroxylation. The barrier-free collapse of C1-radical intermediate IM1-H results in the hydroxylated intermediate IM1-OH with a dramatically exothermic effect (51.81 to 58.20 kcal/mol relative to the respective RCs). The H2O-assisted O-dealkylation of IM1-OH undergoes transition state TSBCIPP with an energy barrier of 7.64 kcal/mol to yield the exothermic product PBCIPP with chloroacetone (Figure 2b).
Biotransformation of TCEP to BCEP. TCEP can undergo C1-hydroxylation and O-dealkylation to produce bis(2-chloroethyl) phosphate (BCEP) and chloroacetaldehyde as metabolites. The free energy profiles are shown in Figure 3. The LS/HS barrier of HAT is 18.04/16.78 kcal/mol, much smaller than those in TDCIPP (27.89/27.52 kcal/mol) and TCIPP (24.27/21.17 kcal/mol) biotransformation, suggesting that TCEP has higher metabolic activity to form the diester metabolite. After OH rebound, the resultant formation of IM1-OH is strongly exothermic by 55.89/55.59 kcal/mol. Moreover, the O-dealkylation of IM1-OH proceeds through transition state TSBCEP with barrier of 7.45 kcal/mol to form the exothermic PBCEP.
Biotransformation of TEP to DEP. Compared with TCEP, TEP has no chlorine substituent in its ethyl groups. The HAT barrier of TEP is 17.34/17.58 kcal/mol in the LS/HS state (Figure 4), which is comparable with that of TCEP (18.04/16.78 kcal/mol), suggesting that chlorine substituent is incapable of differentiating the metabolic activities of TCEP and TEP. OH rebound in HS surface cross the transition state 4TS1-OH with a minor barrier of 2.25 kcal/mol to yield exothermic C1-hydroxyTEP, followed by the H2O-triggered O-dealkylation of IM1-OH to generate the exothermic diethyl phosphate (DEP) with acetaldehyde.
Biotransformation of EHDPHP to DPHP. The free energy profiles of biotransformation of EHDPHP to DPHP are shown in Figure 5. The corresponding HAT barrier is 17.48/16.89 kcal/mol in the LS/HS surface, which has a slight difference with those of TCEP and TEP (18.04/16.78 and 17.34/17.58 kcal/mol). The results state clearly that these three OPFRs have comparable metabolic activities leading to their diester metabolites. The formation of remarkably exothermic IM1-OH facilitates the ensuing O-dealkylation to produce exothermic DPHP with 2-ethylhexanal. Moreover, the C1-hydroxyEHDPHP has much smaller O-dealkylation barrier than the other 4 C1-hydroxyOPFRs (0.27 vs. 6.80 to 8.26 kcal/mol). It is notable that prior studies have confirmed that TPHP containing three phenyl groups can be transformed to diester metabolite DPHP by CYPs [6,32]. Given the structural similarity of EHDPHP with TPHP, it is also possible that EHDPHP metabolism produce 2-ethylhexyl phenyl phosphate (EHPHP) as another diester metabolite.
Overall, these results demonstrate that the HAT barriers are much higher than the O-dealkylation barriers for these 5 OPFRs, suggesting HAT to be the rate-determining step of biotransformation routes. The O-dealkylation of hydroxyOPFRs is a moderately exothermic process with a relatively low energy barrier, facilitating the formation of diester metabolites. In addition, a comparison of HAT barriers of 5 OPFRs reveals that the metabolic activity towards diester metabolites may follow the order of TCEP ≈ TEP ≈ EHDPHP > TCIPP > TDCIPP.

3.2. Binding Modes of 5 OPFRs in CYP3A4

We further selected CYP3A4 isoform as a potential receptor of OPFRs and performed MD simulations to unravel the binding interactions between OPFRs and CYP3A4 and confirm the DFT results discussed above. We calculated the root-mean-square deviations (RMSDs) of CYP3A4 and the bound OPFRs to evaluate the conformational stability of the binding complexes. The smooth RMSD curves indicate the dynamic equilibrium of binding conformations along the simulation time (Figure S3). Figure 6 shows the averaged binding conformations extracted from the last 20 ns MD trajectories. Structural analyses of the binding conformations reveal that CYP3A4 can trap the flexible OPFRs into the active pocket, and thus it may have the potential to metabolize OPFRs.
Among the 5 binding complexes, OPFRs show different binding features. Except for EHDPHP, the remaining 4 OPFRs can orientate one of their Cα atoms (denoted C1 above) toward the Fe=O unit. As shown in Figure 6, the averaged distances between the H atom linked to Cα atoms and the O atom of the Fe=O unit are calculated to be 3.05, 2.75, 2.90, and 3.17 Å for TDCIPP, TCIPP, TCEP, and TEP, respectively, which state clearly that the Cα atoms of these 4 OPFRs are the preferred sites of metabolism by CYP3A4 to yield their Cα-hydroxylated metabolites, followed by O-dealkylation to generate the diester metabolites. These results accord well with the DFT results discussed above. Interestingly, for the CYP3A4-EHDPHP complex, the Cα atom of the 2-ethylhexyl group of EHDPHP stays away from the Fe=O unit, rejecting undoubtedly the formation of DPHP metabolite. However, the averaged distance between the para-H atom of phenyl and the O atom of Fe=O is 3.35 Å, facilitating the preferential formation of para-hydroxyEHDPHP metabolite. Therefore, we conclude that DPHP is not the candidate metabolite involved in the CYP3A4-mediated metabolism of EHDPHP.

3.3. Binding Affinities of 5 OPFRs in CYP3A4

Based on the last 20 ns MD snapshots, we calculated the binding free energies (ΔGbind) to estimate the binding affinities of 5 OPFRs in CYP3A4. The calculated ΔGbind and its energy components for 5 CYP3A4-OPFR complexes are shown in Table 1. Among the 5 OPFRs, EHDPHP has the most negative value of ΔGbind (−26.71 kcal/mol), followed by TDCIPP (−25.02 kcal/mol), TCIPP (−20.70 kcal/mol), TCEP (−13.32 kcal/mol), and TEP (−12.59 kcal/mol). Consequently, the binding affinities of OPFRs in CYP3A4 follow the order of EHDPHP > TDCIPP > TCIPP > TCEP > TEP. The further energy decomposition of ΔGbind indicates that the van der Waals interaction energies (ΔEvdW: −25.28 to −43.85 kcal/mol) possess the most significant contribution to OPFRs binding, followed by the nonpolar solvation free energies (ΔGSA: −4.02 to −6.88 kcal/mol) and electrostatic interaction energies (ΔEele: −0.56 to −4.55 kcal/mol) that have only minor contributions. Instead, entropy component (TΔS: −15.64 to −21.72 kcal/mol) and polar solvation free energies (ΔGGB: 2.92 to 6.30 kcal/mol) show negative contributions to OPFRs binding.

3.4. Effects of CYP3A4 Active Site Residues on OPFRs Binding

Based on the binding energy decomposition calculations, we further identified the key CYP3A4 residues contributing significantly to OPFRs binding (Figure 6). Residues with an energy contribution value of less than −1 kcal/mol are regarded as the key residues [73]. The binding conformations of 5 OPFRs have certain differences in the CYP3A4 active site, and only hydrophobic Ile301 and Phe304 are the shared residues contributing to the binding of these ligands. In particular, Ile301 residue gives the major driving force for the binding of TDCIPP, TEP, and EHDPHP with energy contributions of −2.21, −1.53, and −2.54 kcal/mol (Figure 6a,d,e). Besides, Arg105 residue plays important role in promoting the binding of TDCIPP, TCIPP, and EHDPHP, while it has no contribution to binding TCEP and TEP. Compared with TECP and TEP, TDCIPP and TCIPP have more complex Cl-substituted alkyl groups, and thus their binding involves more key residues. Seven residues provide important contributions to TDCIPP binding, including Arg105, Ser119, Leu211, Phe241, Ile301, Phe304, and Leu482 (Figure 6a). TCIPP has only one Cl substituent in each propyl group compared with TDCIPP, and six residues are considered to be important to its binding, including Arg105, Ser119, Ile301, Phe304, Ala370, and Leu482 (Figure 6b). Further analysis shows that 5 shared residues (Arg105, Ser119, Ile301, Phe304, and Leu482) contribute significantly to the binding of both TDCIPP and TCIPP, which suggests that the different Cl substitution effects can moderately differentiate their binding features. For CYP3A4-EHDPHP complex, a total of 8 key residues are identified to contribute EHDPHP binding, including Arg105 (−1.50 kcal/mol), Phe108 (−2.23 kcal/mol), Ile120 (−1.17 kcal/mol), Leu211 (−2.00 kcal/mol), Phe241 (−1.44 kcal/mol), Ile301 (−2.54 kcal/mol), Phe304 (−1.30 kcal/mol), and Leu482 (−1.39 kcal/mol) (Figure 6e). The number of key residues with significant contributions in the CYP3A4-EHDPHP complex is higher than that in other complexes, and thus EHDPHP shows a stronger binding affinity to CYP3A4. Moreover, EHDPHP orients one of its phenyl groups toward Cpd I and moves the other two side chains into a cavity formed by the hydrophobic residues. Actually, except for Arg105 and Ser119, the other key residues are hydrophobic in CYP3A4-OPFR complexes, which can provide a hydrophobic cavity to stabilize the ligands binding.

4. Conclusions

CYP-catalyzed biotransformation of OPFRs may alert their environmental behavior and toxicological effects, and thus the identification of biotransformation profiles and relevant metabolites can provide meaningful guidance for environmental risk assessment of OPFRs exposure. In this work, we first performed quantum chemical calculations to unravel the biotransformation mechanism of OPFRs to their diester metabolites. The results show that Cpd I-mediated biotransformation of OPFRs undergoes Cα-hydroxylation comprising of rate-determining HAT and barrier-free OH rebound to generate Cα-hydroxyOPFRs, followed by H2O-triggered O-dealkylation to yield the respective diester metabolites. Metabolic activities of OPFRs follow the order of TCEP ≈ TEP ≈ EHDPHP > TCIPP > TDCIPP. Compared with the energy profile of TCIPP biotransformation, the extra Cl substituent in TDCIPP enhances the HAT barrier and thus reduces the metabolic activity toward BDCIPP. However, the Cl substitution cannot differentiate the reactivity of TCEP and TEP.
We further performed MD simulations to explore the binding interactions of OPFRs in CYP3A4 at the molecular level. The MD results demonstrate that CYP3A4 can accommodate these 5 OPFRs with varying binding features. Binding mode analyses indicate that the Cα atoms of TDCIPP, TCIPP, TCEP, and TEP are the preferred sites of metabolism leading to their diester metabolites, while EHDPHP may have more preference for para-hydroxyEHDPHP metabolite rather than DPHP. Based on the binding free energy calculations, we estimate the binding affinities in the order of EHDPHP > TDCIPP > TCIPP > TCEP > TEP and van der Waals interactions contribute to the major driving forces for OPFRs binding. Furthermore, the hydrophobic residues in the CYP3A4 active site play a crucial role in stabilizing the OPFRs binding, among which Ile301 and Phe304 are the common key residues contributing to the binding. Overall, the present work reports the CYP-mediated biotransformation profiles of OPFRs to their diester metabolites, which may provide an essential understanding of the metabolic fates of OPFRs.

Supplementary Materials

The following supporting information can be downloaded at: https://0-www-mdpi-com.brum.beds.ac.uk/article/10.3390/molecules27092799/s1, Figure S1: Structural diagram of 5 OPFRs used in this work; Figure S2: Optimal docking conformations of OPFRs in 3A4 active site; Figure S3: Root mean square deviations (RMSDs) of 5 CYP3A4-OPFR complexes; Table S1: Structural characteristics and activation energy barriers of the transition states of hydrogen atom transfer and OH rebound; Tables S2–S6: Absolute energies and relative energies for the reaction species involved in Cpd I-catalyzed biotransformation of OPFRs to diesters; Tables S7–S11: Calculated spin densities for the molecular species involved in OPFRs C1-hydroxylation; Cartesian coordinates of all DFT-optimized structures.

Author Contributions

Conceptualization, T.Y. and G.M.; methodology, Y.J., T.Y. and G.M.; software, Q.X. and X.Z.; validation, H.D. and X.W.; formal analysis, Q.X., X.Z. and Z.W.; investigation, Y.J. and T.Y.; resources, G.M. and Z.W.; data curation, H.D. and X.W.; writing—original draft preparation, Y.J. and T.Y.; writing—review and editing, G.M. and H.Y.; visualization, X.W. and Z.W.; supervision, G.M.; project administration, G.M. and H.Y.; funding acquisition, G.M. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Natural Science Foundation of China (22176177), Natural Science Foundation of Zhejiang Province (LY20B070005), and National College Students Innovation and Entrepreneurship Training Program (201910345038).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Sample Availability

Not available.

References

  1. Hoffman, K.; Daniels, J.L.; Stapleton, H.M. Urinary metabolites of organophosphate flame retardants and their variability in pregnant women. Environ. Int. 2014, 63, 169–172. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  2. Stapleton, H.M.; Klosterhaus, S.; Eagle, S.; Fuh, J.; Meeker, J.D.; Blum, A.; Webster, T.F. Detection of Organophosphate Flame Retardants in Furniture Foam and US House Dust. Environ. Sci. Technol. 2009, 43, 7490–7495. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  3. Kim, J.-W.; Isobe, T.; Chang, K.-H.; Amano, A.; Maneja, R.H.; Zamora, P.B.; Siringan, F.P.; Tanabe, S. Levels and distribution of organophosphorus flame retardants and plasticizers in fishes from Manila Bay, the Philippines. Environ. Pollut. 2011, 159, 3653–3659. [Google Scholar] [CrossRef] [PubMed]
  4. Stapleton, H.M.; Sharma, S.; Getzinger, G.; Ferguson, P.L.; Gabriel, M.; Webster, T.F.; Blum, A. Novel and High Volume Use Flame Retardants in US Couches Reflective of the 2005 PentaBDE Phase Out. Environ. Sci. Technol. 2012, 46, 13432–13439. [Google Scholar] [CrossRef]
  5. Butt, C.M.; Congleton, J.; Hoffman, K.; Fang, M.; Stapleton, H.M. Metabolites of Organophosphate Flame Retardants and 2-Ethylhexyl Tetrabromobenzoate in Urine from Paired Mothers and Toddlers. Environ. Sci. Technol. 2014, 48, 10432–10438. [Google Scholar] [CrossRef]
  6. Zhang, Q.; Ji, S.; Chai, L.; Yang, F.; Zhao, M.; Liu, W.; Schueuermann, G.; Ji, L. Metabolic Mechanism of Aryl Phosphorus Flame Retardants by Cytochromes P450: A Combined Experimental and Computational Study on Triphenyl Phosphate. Environ. Sci. Technol. 2018, 52, 14411–14421. [Google Scholar] [CrossRef]
  7. Zhang, Q.; Lu, M.; Dong, X.; Wang, C.; Zhang, C.; Liu, W.; Zhao, M. Potential Estrogenic Effects of Phosphorus-Containing Flame Retardants. Environ. Sci. Technol. 2014, 48, 6995–7001. [Google Scholar] [CrossRef]
  8. Wang, R.; Tang, J.; Xie, Z.; Mi, W.; Chen, Y.; Wolschke, H.; Tian, C.; Pan, X.; Luo, Y.; Ebinghaus, R. Occurrence and spatial distribution of organophosphate ester flame retardants and plasticizers in 40 rivers draining into the Bohai Sea, north China. Environ. Pollut. 2015, 198, 172–178. [Google Scholar] [CrossRef] [Green Version]
  9. Lai, S.; Xie, Z.; Song, T.; Tang, J.; Zhang, Y.; Mi, W.; Peng, J.; Zhao, Y.; Zou, S.; Ebinghaus, R. Occurrence and dry deposition of organophosphate esters in atmospheric particles over the northern South China Sea. Chemosphere 2015, 127, 195–200. [Google Scholar] [CrossRef]
  10. Marklund, A.; Andersson, B.; Haglund, P. Traffic as a source of organophosphorus flame retardants and plasticizers in snow. Environ. Sci. Technol. 2005, 39, 3555–3562. [Google Scholar] [CrossRef]
  11. Cequier, E.; Marcé, R.M.; Becher, G.; Thomsen, C. A high-throughput method for determination of metabolites of organophosphate flame retardants in urine by ultra performance liquid chromatography–high resolution mass spectrometry. Anal. Chim. Acta 2014, 845, 98–104. [Google Scholar] [CrossRef] [PubMed]
  12. Chung, H.-W.; Ding, W.-H. Determination of organophosphate flame retardants in sediments by microwave-assisted extraction and gas chromatography–mass spectrometry with electron impact and chemical ionization. Anal. Bioanal. Chem. 2009, 395, 2325–2334. [Google Scholar] [CrossRef] [PubMed]
  13. Gao, Z.; Deng, Y.; Yuan, W.; He, H.; Yang, S.; Sun, C. Determination of organophosphorus flame retardants in fish by pressurized liquid extraction using aqueous solutions and solid-phase microextraction coupled with gas chromatography-flame photometric detector. J. Chromatogr. A 2014, 1366, 31–37. [Google Scholar] [CrossRef]
  14. Staaf, T.; Östman, C. Organophosphate triesters in indoor environments. J. Environ. Monit. 2005, 7, 883–887. [Google Scholar] [CrossRef] [PubMed]
  15. Jiao, E.; Hu, X.; Li, L.; Zhang, H.; Zhu, Z.; Yin, D.; Qiu, Y. Occurrence and risk evaluation of organophosphorus flame retardants in two urban rivers in Yangtze River Delta. Environ. Monit. Assess. 2021, 193, 146. [Google Scholar] [CrossRef]
  16. Zhang, L.; Lu, L.; Zhu, W.; Yang, B.; Lu, D.; Dan, S.F.; Zhang, S. Organophosphorus flame retardants (OPFRs) in the seawater and sediments of the Qinzhou Bay, Northern Beibu Gulf: Occurrence, distribution, and ecological risks. Mar. Pollut. Bull. 2021, 168, 112368. [Google Scholar] [CrossRef]
  17. Chupeau, Z.; Bonvallot, N.; Mercier, F.; Le Bot, B.; Chevrier, C.; Glorennec, P. Organophosphorus Flame Retardants: A Global Review of Indoor Contamination and Human Exposure in Europe and Epidemiological Evidence. Int. J. Environ. Res. Public Health 2020, 17, 6713. [Google Scholar] [CrossRef]
  18. Ali, N.; Dirtu, A.C.; Van den Eede, N.; Goosey, E.; Harrad, S.; Neels, H.; Coakley, J.; Douwes, J.; Covaci, A. Occurrence of alternative flame retardants in indoor dust from New Zealand: Indoor sources and human exposure assessment. Chemosphere 2012, 88, 1276–1282. [Google Scholar] [CrossRef]
  19. Or, P. TSCA Work Plan Chemical Problem Formulation and Initial Assessment Chlorinated Phosphate Ester Cluster Flame Retardants; United States Environmental Protection Agency: Washington, DC, USA, 2015.
  20. Zhang, X.; Zou, W.; Mu, L.; Chen, Y.; Ren, C.; Hu, X.; Zhou, Q. Rice ingestion is a major pathway for human exposure to organophosphate flame retardants (OPFRs) in China. J. Hazard. Mater. 2016, 318, 686–693. [Google Scholar] [CrossRef]
  21. Kucharska, A.; Cequier, E.; Thomsen, C.; Becher, G.; Covaci, A.; Voorspoels, S. Assessment of human hair as an indicator of exposure to organophosphate flame retardants. Case study on a Norwegian mother–child cohort. Environ. Int. 2015, 83, 50–57. [Google Scholar] [CrossRef]
  22. Kucharska, A.; Covaci, A.; Vanermen, G.; Voorspoels, S. Development of a broad spectrum method for measuring flame retardants-overcoming the challenges of non-invasive human biomonitoring studies. Anal. Bioanal. Chem. 2014, 406, 6665–6675. [Google Scholar] [CrossRef] [PubMed]
  23. Chen, X.; Zhao, X.; Shi, Z. Organophosphorus flame retardants in breast milk from Beijing, China: Occurrence, nursing infant’s exposure and risk assessment. Sci. Total Environ. 2021, 771, 145404. [Google Scholar] [CrossRef] [PubMed]
  24. Dishaw, L.V.; Powers, C.M.; Ryde, I.T.; Roberts, S.C.; Seidler, F.J.; Slotkin, T.A.; Stapleton, H.M. Is the PentaBDE replacement, tris (1,3-dichloro-2-propyl) phosphate (TDCPP), a developmental neurotoxicant? Studies in PC12 cells. Toxicol. Appl. Pharmacol. 2011, 256, 281–289. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  25. Liu, X.; Ji, K.; Choi, K. Endocrine disruption potentials of organophosphate flame retardants and related mechanisms in H295R and MVLN cell lines and in zebrafish. Aquat. Toxicol. 2012, 114–115, 173–181. [Google Scholar] [CrossRef] [PubMed]
  26. Porter, E.; Crump, D.; Egloff, C.; Chiu, S.; Kennedy, S.W. Use of an avian hepatocyte assay and the avian toxchip polymerse chain reaction array for testing prioritization of 16 organic flame retardants. Environ. Toxicol. Chem. 2014, 33, 573–582. [Google Scholar] [CrossRef]
  27. van der Veen, I.; de Boer, J. Phosphorus flame retardants: Properties, production, environmental occurrence, toxicity and analysis. Chemosphere 2012, 88, 1119–1153. [Google Scholar] [CrossRef]
  28. Hogberg, H.T.; da Silveira E Sa, R.d.C.; Kleensang, A.; Bouhifd, M.; Cemiloglu Ulker, O.; Smirnova, L.; Behl, M.; Maertens, A.; Zhao, L.; Hartung, T. Organophosphorus flame retardants are developmental neurotoxicants in a rat primary brainsphere in vitro model. Arch. Toxicol. 2021, 95, 207–228. [Google Scholar] [CrossRef]
  29. Behl, M.; Hsieh, J.-H.; Shafer, T.J.; Mundy, W.R.; Rice, J.R.; Boyd, W.A.; Freedman, J.H.; Hunter, E.S.; Jarema, K.A.; Padilla, S.; et al. Use of alternative assays to identify and prioritize organophosphorus flame retardants for potential developmental and neurotoxicity. Neurotoxicol. Teratol. 2015, 52, 181–193. [Google Scholar] [CrossRef]
  30. Kirchmair, J.; Göller, A.H.; Lang, D.; Kunze, J.; Testa, B.; Wilson, I.D.; Glen, R.C.; Schneider, G. Predicting drug metabolism: Experiment and/or computation? Nat. Rev. Drug Discov. 2015, 14, 387–404. [Google Scholar] [CrossRef] [Green Version]
  31. Hou, R.; Xu, Y.; Wang, Z. Review of OPFRs in animals and humans: Absorption, bioaccumulation, metabolism, and internal exposure research. Chemosphere 2016, 153, 78–90. [Google Scholar] [CrossRef]
  32. Van den Eede, N.; Maho, W.; Erratico, C.; Neels, H.; Covaci, A. First insights in the metabolism of phosphate flame retardants and plasticizers using human liver fractions. Toxicol. Lett. 2013, 223, 9–15. [Google Scholar] [CrossRef] [PubMed]
  33. Wei, K.; Yin, H.; Peng, H.; Lu, G.; Dang, Z. Bioremediation of triphenyl phosphate by Brevibacillus brevis: Degradation characteristics and role of cytochrome P450 monooxygenase. Sci. Total Environ. 2018, 627, 1389–1395. [Google Scholar] [CrossRef] [PubMed]
  34. Chai, L.; Zhang, H.; Song, R.; Yang, H.; Yu, H.; Paneth, P.; Kepp, K.P.; Akamatsu, M.; Ji, L. Precision Biotransformation of Emerging Pollutants by Human Cytochrome P450 Using Computational–Experimental Synergy: A Case Study of Tris(1,3-dichloro-2-propyl) Phosphate. Environ. Sci. Technol. 2021, 55, 14037–14050. [Google Scholar] [CrossRef] [PubMed]
  35. Ballesteros-Gómez, A.; Erratico, C.A.; Eede, N.V.d.; Ionas, A.C.; Leonards, P.E.G.; Covaci, A. In vitro metabolism of 2-ethylhexyldiphenyl phosphate (EHDPHP) by human liver microsomes. Toxicol. Lett. 2015, 232, 203–212. [Google Scholar] [CrossRef] [PubMed]
  36. Chen, Z.; Xie, J.; Li, Q.; Hu, K.; Yang, Z.; Yu, H.; Liu, Y. Human CYP enzyme-activated clastogenicity of 2-ethylhexyl diphenyl phosphate (a flame retardant) in mammalian cells. Environ. Pollut. 2021, 285, 117527. [Google Scholar] [CrossRef] [PubMed]
  37. Hou, R.; Huang, C.; Rao, K.; Xu, Y.; Wang, Z. Characterized in Vitro Metabolism Kinetics of Alkyl Organophosphate Esters in Fish Liver and Intestinal Microsomes. Environ. Sci. Technol. 2018, 52, 3202–3210. [Google Scholar] [CrossRef]
  38. Cui, S.; Yu, Y.; Zhan, T.; Gao, Y.; Zhang, J.; Zhang, L.; Ge, Z.; Liu, W.; Zhang, C.; Zhuang, S. Carcinogenic Risk of 2,6-Di-tert-Butylphenol and Its Quinone Metabolite 2,6-DTBQ Through Their Interruption of RARβ: In Vivo, In Vitro, and In Silico Investigations. Environ. Sci. Technol. 2022, 56, 480–490. [Google Scholar] [CrossRef]
  39. Han, C.; Zhu, W.; Ma, G.; Chen, Y.; Li, X.; Wei, X.; Yu, H. Computational insight into biotransformation of halophenols by cytochrome P450: Mechanism and reactivity for epoxidation. Chemosphere 2022, 286, 131708. [Google Scholar] [CrossRef]
  40. Ma, G.; Geng, L.; Lu, Y.; Wei, X.; Yu, H. Investigating the molecular mechanism of hydroxylated bromdiphenyl ethers to inhibit the thyroid hormone sulfotransferase SULT1A1. Chemosphere 2021, 263, 128353. [Google Scholar] [CrossRef]
  41. Ma, G.; Yu, H.; Han, C.; Jia, Y.; Wei, X.; Wang, Z. Binding and Metabolism of Brominated Flame Retardant β-1,2-Dibromo-4-(1,2-dibromoethyl)cyclohexane in Human Microsomal P450 Enzymes: Insights from Computational Studies. Chem. Res. Toxicol. 2020, 33, 1487–1496. [Google Scholar] [CrossRef]
  42. Zhang, R.; Li, P.; Shi, X.; Zhang, R.; Wang, J.; Li, Y.; Zhang, Q.; Wang, W. Insights into the metabolic mechanism of PBDEs catalyzed by cytochrome P450 enzyme 3A4: A QM/MM study. Chemosphere 2021, 278, 130430. [Google Scholar] [CrossRef] [PubMed]
  43. Mirzaei, M.S.; Ivanov, M.V.; Taherpour, A.A.; Mirzaei, S. Mechanism-Based Inactivation of Cytochrome P450 Enzymes: Computational Insights. Chem. Res. Toxicol. 2021, 34, 959–987. [Google Scholar] [CrossRef] [PubMed]
  44. Fu, Z.; Wang, Y.; Chen, J.; Wang, Z.; Wang, X. How PBDEs Are Transformed into Dihydroxylated and Dioxin Metabolites Catalyzed by the Active Center of Cytochrome P450s: A DFT Study. Environ. Sci. Technol. 2016, 50, 8155–8163. [Google Scholar] [CrossRef] [PubMed]
  45. Ji, L.; Ji, S.; Wang, C.; Kepp, K.P. Molecular Mechanism of Alternative P450-Catalyzed Metabolism of Environmental Phenolic Endocrine-Disrupting Chemicals. Environ. Sci. Technol. 2018, 52, 4422–4431. [Google Scholar] [CrossRef] [Green Version]
  46. Ji, L.; Schueuermann, G. Model and Mechanism: N-Hydroxylation of Primary Aromatic Amines by Cytochrome P450. Angew. Chem. Int. Ed. 2013, 52, 744–748. [Google Scholar] [CrossRef]
  47. Wang, X.; Wang, Y.; Chen, J.; Ma, Y.; Zhou, J.; Fu, Z. Computational Toxicological Investigation on the Mechanism and Pathways of Xenobiotics Metabolized by Cytochrome P450: A Case of BDE-47. Environ. Sci. Technol. 2012, 46, 5126–5133. [Google Scholar] [CrossRef]
  48. Frisch, M.J.; Trucks, G.W.; Schlegel, H.B.; Scuseria, G.E.; Robb, M.A.; Cheeseman, J.R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G.A.; et al. Gaussian 09, Revision D.01; Gaussian, Inc.: Wallingford, CT, USA, 2013. [Google Scholar]
  49. Becke, A.D. Density-functional exchange-energy approximation with correct asymptotic behavior. Phys. Rev. A Gen. Phys. 1988, 38, 3098–3100. [Google Scholar] [CrossRef]
  50. Becke, A.D. Density-functional thermochemistry. III. The role of exact exchange. J. Chem. Phys. 1993, 98, 5648–5652. [Google Scholar] [CrossRef] [Green Version]
  51. Lee, C.; Yang, W.T.; Parr, R.G. Development of the ColleSalvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B Condens. Matter 1988, 37, 785–789. [Google Scholar] [CrossRef] [Green Version]
  52. Hay, P.J.; Wadt, W.R. Abinitio effective core potentials for molecular calculations-potentials for the transition-metal atoms sc to hg. J. Chem. Phys. 1985, 82, 270–283. [Google Scholar] [CrossRef]
  53. Hehre, W.J.; Ditchfield, R.; Pople, J.A. Self-consistent molecular orbital methods. XII. Further extensions of Gaussian-type basis sets for use in molecular orbital studies of organic molecules. J. Chem. Phys. 1972, 56, 2257–2261. [Google Scholar] [CrossRef]
  54. Weigend, F.; Ahlrichs, R. Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy. Phys. Chem. Chem. Phys. 2005, 7, 3297–3305. [Google Scholar] [CrossRef] [PubMed]
  55. Grimme, S. Semiempirical GGA-type density functional constructed with a long-range dispersion correction. J. Comput. Chem. 2006, 27, 1787–1799. [Google Scholar] [CrossRef] [PubMed]
  56. Marenich, A.V.; Cramer, C.J.; Truhlar, D.G. Universal Solvation Model Based on Solute Electron Density and on a Continuum Model of the Solvent Defined by the Bulk Dielectric Constant and Atomic Surface Tensions. J. Phys. Chem. B 2009, 113, 6378–6396. [Google Scholar] [CrossRef] [PubMed]
  57. Shaik, S.; Cohen, S.; Wang, Y.; Chen, H.; Kumar, D.; Thiel, W. P450 Enzymes: Their Structure, Reactivity, and Selectivity-Modeled by QM/MM Calculations. Chem. Rev. 2010, 110, 949–1017. [Google Scholar] [CrossRef] [PubMed]
  58. Shaik, S.; Kumar, D.; de Visser, S.P.; Altun, A.; Thiel, W. Theoretical perspective on the structure and mechanism of cytochrome P450 enzymes. Chem. Rev. 2005, 105, 2279–2328. [Google Scholar] [CrossRef]
  59. Ekroos, M.; Sjögren, T. Structural basis for ligand promiscuity in cytochrome P450 3A4. Proc. Natl. Acad. Sci. USA 2006, 103, 13682–13687. [Google Scholar] [CrossRef] [Green Version]
  60. Bren, U.; Oostenbrink, C. Cytochrome P450 3A4 Inhibition by Ketoconazole: Tackling the Problem of Ligand Cooperativity Using Molecular Dynamics Simulations and Free-Energy Calculations. J. Chem. Inf. Model. 2012, 52, 1573–1582. [Google Scholar] [CrossRef]
  61. Guengerich, F.P. Cytochrome P450 and Chemical Toxicology. Chem. Res. Toxicol. 2008, 21, 70–83. [Google Scholar] [CrossRef]
  62. Burley, S.K.; Bhikadiya, C.; Bi, C.; Bittrich, S.; Chen, L.; Crichlow, G.V.; Christie, C.H.; Dalenberg, K.; Di Costanzo, L.; Duarte, J.M.; et al. RCSB Protein Data Bank: Powerful new tools for exploring 3D structures of biological macromolecules for basic and applied research and education in fundamental biology, biomedicine, biotechnology, bioengineering and energy sciences. Nucleic Acids Res. 2020, 49, D437–D451. [Google Scholar] [CrossRef]
  63. Trott, O.; Olson, A.J. AutoDock Vina: Improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading. J. Comput. Chem. 2010, 31, 455–461. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  64. Pettersen, E.F.; Goddard, T.D.; Huang, C.C.; Couch, G.S.; Greenblatt, D.M.; Meng, E.C.; Ferrin, T.E. UCSF Chimera—A visualization system for exploratory research and analysis. J. Comput. Chem. 2004, 25, 1605–1612. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  65. Case, D.A.; Cheatham Iii, T.E.; Darden, T.; Gohlke, H.; Luo, R.; Merz, K.M., Jr.; Onufriev, A.; Simmerling, C.; Wang, B.; Woods, R.J. The Amber biomolecular simulation programs. J. Comput. Chem. 2005, 26, 1668–1688. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  66. Maier, J.A.; Martinez, C.; Kasavajhala, K.; Wickstrom, L.; Hauser, K.E.; Simmerling, C. ff14SB: Improving the accuracy of protein side chain and backbone parameters from ff99SB. J. Chem. Theory Comput. 2015, 11, 3696–3713. [Google Scholar] [CrossRef] [Green Version]
  67. Dolinsky, T.J.; Nielsen, J.E.; McCammon, J.A.; Baker, N.A. PDB2PQR: An automated pipeline for the setup of Poisson–Boltzmann electrostatics calculations. Nucleic Acids Res. 2004, 32, W665–W667. [Google Scholar] [CrossRef]
  68. Jorgensen, W.L.; Chandrasekhar, J.; Madura, J.D.; Impey, R.W.; Klein, M.L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79, 926–935. [Google Scholar] [CrossRef]
  69. Kollman, P.A.; Massova, I.; Reyes, C.; Kuhn, B.; Huo, S.; Chong, L.; Lee, M.; Lee, T.; Duan, Y.; Wang, W.; et al. Calculating Structures and Free Energies of Complex Molecules:  Combining Molecular Mechanics and Continuum Models. Acc. Chem. Res. 2000, 33, 889–897. [Google Scholar] [CrossRef]
  70. Gohlke, H.; Kiel, C.; Case, D.A. Insights into Protein–Protein Binding by Binding Free Energy Calculation and Free Energy Decomposition for the Ras–Raf and Ras–RalGDS Complexes. J. Mol. Biol. 2003, 330, 891–913. [Google Scholar] [CrossRef]
  71. Carignan, C.C.; McClean, M.D.; Cooper, E.M.; Watkins, D.J.; Fraser, A.J.; Heiger-Bernays, W.; Stapleton, H.M.; Webster, T.F. Predictors of tris(1,3-dichloro-2-propyl) phosphate metabolite in the urine of office workers. Environ. Int. 2013, 55, 56–61. [Google Scholar] [CrossRef] [Green Version]
  72. Cooper, E.M.; Covaci, A.; van Nuijs, A.L.N.; Webster, T.F.; Stapleton, H.M. Analysis of the flame retardant metabolites bis(1,3-dichloro-2-propyl) phosphate (BDCPP) and diphenyl phosphate (DPP) in urine using liquid chromatography-tandem mass spectrometry. Anal. Bioanal. Chem. 2011, 401, 2123–2132. [Google Scholar] [CrossRef] [Green Version]
  73. Ji, B.; Liu, S.; He, X.; Man, V.H.; Xie, X.-Q.; Wang, J. Prediction of the Binding Affinities and Selectivity for CB1 and CB2 Ligands Using Homology Modeling, Molecular Docking, Molecular Dynamics Simulations, and MM-PBSA Binding Free Energy Calculations. ACS Chem. Neurosci. 2020, 11, 1139–1158. [Google Scholar] [CrossRef] [PubMed]
Molecules 27 02799 g001 550
Figure 1. Free energy profiles for (a) C1-hydroxylation of TDCIPP to IM1-OH by Cpd I and (b) O-dealkylation of IM1-OH to BDCIPP, along with the optimized geometries of relevant reaction species in HS and LS states. Bond lengths are given in angstroms, while relative energies are given in kcal/mol.
Figure 1. Free energy profiles for (a) C1-hydroxylation of TDCIPP to IM1-OH by Cpd I and (b) O-dealkylation of IM1-OH to BDCIPP, along with the optimized geometries of relevant reaction species in HS and LS states. Bond lengths are given in angstroms, while relative energies are given in kcal/mol.
Molecules 27 02799 g001
Molecules 27 02799 g002 550
Figure 2. Free energy profiles for (a) C1-hydroxylation of TCIPP to IM1-OH by Cpd I and (b) O-dealkylation of IM1-OH to BCIPP, along with the optimized geometries of relevant reaction species in HS and LS states. Bond lengths are given in angstroms, while relative energies are given in kcal/mol.
Figure 2. Free energy profiles for (a) C1-hydroxylation of TCIPP to IM1-OH by Cpd I and (b) O-dealkylation of IM1-OH to BCIPP, along with the optimized geometries of relevant reaction species in HS and LS states. Bond lengths are given in angstroms, while relative energies are given in kcal/mol.
Molecules 27 02799 g002
Molecules 27 02799 g003 550
Figure 3. Free energy profiles for (a) C1-hydroxylation of TCEP to IM1-OH by Cpd I and (b) O-dealkylation of IM1-OH to BCEP, along with the optimized geometries of relevant reaction species in HS and LS states. Bond lengths are given in angstroms, while relative energies are given in kcal/mol.
Figure 3. Free energy profiles for (a) C1-hydroxylation of TCEP to IM1-OH by Cpd I and (b) O-dealkylation of IM1-OH to BCEP, along with the optimized geometries of relevant reaction species in HS and LS states. Bond lengths are given in angstroms, while relative energies are given in kcal/mol.
Molecules 27 02799 g003
Molecules 27 02799 g004 550
Figure 4. Free energy profiles for (a) C1-hydroxylation of TEP to IM1-OH by Cpd I and (b) O-dealkylation of IM1-OH to DEP, along with the optimized geometries of relevant reaction species in HS and LS states. Bond lengths are given in angstroms, while relative energies are given in kcal/mol.
Figure 4. Free energy profiles for (a) C1-hydroxylation of TEP to IM1-OH by Cpd I and (b) O-dealkylation of IM1-OH to DEP, along with the optimized geometries of relevant reaction species in HS and LS states. Bond lengths are given in angstroms, while relative energies are given in kcal/mol.
Molecules 27 02799 g004
Molecules 27 02799 g005 550
Figure 5. Free energy profiles for (a) C1-hydroxylation of EHDPHP to IM1-OH by Cpd I and (b) O-dealkylation of IM1-OH to DPHP, along with the optimized geometries of relevant reaction species in HS and LS states. Bond lengths are given in angstroms, while relative energies are given in kcal/mol.
Figure 5. Free energy profiles for (a) C1-hydroxylation of EHDPHP to IM1-OH by Cpd I and (b) O-dealkylation of IM1-OH to DPHP, along with the optimized geometries of relevant reaction species in HS and LS states. Bond lengths are given in angstroms, while relative energies are given in kcal/mol.
Molecules 27 02799 g005
Molecules 27 02799 g006 550
Figure 6. Averaged binding conformations of (a) TDCIPP, (b) TCIPP, (c) TCEP, (d) TEP, and (e) EHDPHP in CYP3A4 extracted from the last 20 ns MD trajectories. The key residues contributing to OPFRs binding are shown in the images, and corresponding binding energy contributions are given in kcal/mol and shown in the brackets. Distances are given in angstroms.
Figure 6. Averaged binding conformations of (a) TDCIPP, (b) TCIPP, (c) TCEP, (d) TEP, and (e) EHDPHP in CYP3A4 extracted from the last 20 ns MD trajectories. The key residues contributing to OPFRs binding are shown in the images, and corresponding binding energy contributions are given in kcal/mol and shown in the brackets. Distances are given in angstroms.
Molecules 27 02799 g006
Table
Table 1. Binding free energies between CYP3A4 and OPFRs.
Table 1. Binding free energies between CYP3A4 and OPFRs.
CYPOPFRsEnergy Components (kcal/mol)
ΔEeleΔEVDWΔGGBΔGSATΔSΔGbind
3A4TDCIPP−4.55−40.926.30−5.67−19.82−25.02
TCIPP−0.56−35.352.92−4.93−17.22−20.70
TCEP−3.15−27.704.59−4.23−17.17−13.32
TEP−2.44−25.283.51−4.02−15.64−12.59
EHDPHP−3.31−43.855.61−6.88−21.72−26.71
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Jia, Y.; Yao, T.; Ma, G.; Xu, Q.; Zhao, X.; Ding, H.; Wei, X.; Yu, H.; Wang, Z. Computational Insight into Biotransformation Profiles of Organophosphorus Flame Retardants to Their Diester Metabolites by Cytochrome P450. Molecules 2022, 27, 2799. https://0-doi-org.brum.beds.ac.uk/10.3390/molecules27092799

AMA Style

Jia Y, Yao T, Ma G, Xu Q, Zhao X, Ding H, Wei X, Yu H, Wang Z. Computational Insight into Biotransformation Profiles of Organophosphorus Flame Retardants to Their Diester Metabolites by Cytochrome P450. Molecules. 2022; 27(9):2799. https://0-doi-org.brum.beds.ac.uk/10.3390/molecules27092799

Chicago/Turabian Style

Jia, Yue, Tingji Yao, Guangcai Ma, Qi Xu, Xianglong Zhao, Hui Ding, Xiaoxuan Wei, Haiying Yu, and Zhiguo Wang. 2022. "Computational Insight into Biotransformation Profiles of Organophosphorus Flame Retardants to Their Diester Metabolites by Cytochrome P450" Molecules 27, no. 9: 2799. https://0-doi-org.brum.beds.ac.uk/10.3390/molecules27092799

Article Metrics

Back to TopTop