Human Health Risk Regulation of Reproductive Toxicity, Neurotoxicity, and Endocrine Disruption in Special Populations Exposed to Organophosphorus Flame Retardants

In this study, the joint toxicological characteristics of reproductive toxicity, neurotoxicity, and endocrine disruption (ED) by organophosphorus flame retardants (OPFRs) were regulated by process control. Molecular docking technology, molecular dynamics (MD), 2D-QSAR model, and density functional theory (DFT) were used to develop a health risk regulation scheme for special population such as pregnant women exposed to OPFRs. It was found that MD simulations confirmed the effectiveness of the recommended complementary food scheme (CFS) for the pregnant women with low health risk. When β-lactoglobulin, α-lactoalbumin, milk fat globule membrane (MFGM) protein, ovalbumin (OVA), ovotransferrin (OVT), vitamin, plant pigment, apple polyphenols, and malic acid were present in the CFS, the joint toxicity of OPFRs in pregnant women were significantly decreased by 91.18%. The reproductive toxicity played a dominant role in the joint toxicity and could be reduced by 82.48% under the recommended CFS. There was a competitive relationship between the nutrients in the recommended CFS and OPFRs binding to the joint toxic receptor (JTR). The former could easily occupy the target binding spot of the JTR protein due to the changes of force field of OPFRs docked with JTR, which reduced or prevented the binding of OPFRs to the JTR. In addition, simulation of OPFRs molecular metabolic pathways in pregnant women under the recommended CFS showed that the binding affinity between OPFRs and six metabolic kinases in pregnant women was significantly decreased (–28.85– –87.54%), indicating that the inhibition effect of OPFRs on normal biochemical reactions in the human body was significantly reduced, which to a certain extent verified the effectiveness of the recommended CFS.

As a substitute for POPs such as PDBEs, in terms of harmful effects, OPFRs are not completely separated from POPs. OPFRs still have a variety of biological toxicities, including developmental neurotoxicity, which is associated with adverse reproductive/developmental nervous system effects in animals and humans (Sun et al. 2016;Baldwin et al. 2017;Gu et al. 2018). It was found that the lethal concentration 50% (LC 50 ) and concentration for 50% of maximal effect (EC 50 ) values of OPFRs were similar to the toxicity values of BFRs such as PBDEs, even one or two orders of magnitude less than the LC 50 and EC 50 values of the various flame retardants such as OPFRs, PBDEs, PCBs, HBCD, tetrabromobisphenol A, and DBDPE (Gu et al. 2018;Chu and Li 2019). As the values of LC 50 and EC 50 decreases, the toxicity increases. Therefore, the toxic effect of OPFRs cannot be ignored. TDCPP and TCIPP have been known as carcinogenic and were removed from children's pajamas in the late 1970s, but in recent decades, the production of OPFRs as a flame retardant for other applications has increased significantly (Percy et al. 2020). TCEP and TPHP are neurotoxic to some animals (Gu et al. 2018;Wang et al. 2015). Besides, TPHP had been used as a phthalate substitute in some personal care products (e.g., nail polish), providing another potential exposure route of OPFRs hazards to the environment and human health (Young et al. 2018).
As additive flame retardants, OPFRs exist in products through physical combination. During the manufacturing, processing, and use of products, OPFRs easily enter into the environment and human food chain through leaching, deposition, inhalation, ingestion, and dermal absorption and accumulate in human milk and urine (Wei et al. 2015;Zhang et al. 2016;Ospina et al. 2018;Schreder et al. 2015). Hand-to-mouth contact of the human body (dust inhalation and skin absorption) to OPFRs is one of the most common OPFR absorption pathways for people living in indoor environments (Hou et al. 2016). Dietary intake is also an important way for the human body to absorb OPFRs (Fan et al. 2014). The widespread existence of OPFRs indoors implies that pregnant women are potentially exposed to them, and thus the risk of OPFRs molecular toxicity to pregnant women has aroused widespread concern in late years. Recent studies have shown that TDCPP, TCP, and TPP with a median concentration of about 1-6 µg/g household dust were inadvertently inhaled by pregnant women or growing infants (Meeker et al. 2013;Carignan et al. 2013). TDCPP, TCP, and TPP can be detected in urine, breast milk, and blood samples of pregnant women at concentrations ranging from 1-100 ng/mL or 1-100 ng/g lipids. Metabolites of OPFRs were detected in pregnant women during the perinatal period (Hoffman et al. 2015), which were associated with an increased risk of high body mass index (BMI) and adverse neurological outcomes in children (Boyle et al. 2019;Doherty et al. 2019). Bis(1, 3-dichloroisopropyl) phosphate (BDCIPP), the metabolite of tris(1,3-dichloro-2-propyl) phosphate (TDCPP), was detected in the urine samples of the children across the entire United States at a median concentration of ~ 2.8 ng/mL (Butt et al. 2016). BDCIPP was detected in all the urine samples of the Australian children, with a median concentration of 7.8 µg/g (He et al. 2018). The presence of OPFRs metabolites such as diphenyl phosphate (DPHP), bis(1,3-dichloro-2-propyl) phosphate (BDCIPP), and bis(2-chloroethyl) phosphate (BCEP) in the urine samples of the pregnant women (n = 357) at 16 and 26 weeks of gestation and at birth confirmed that the pregnant women were indeed exposed to OPFRs during fetal development (Percy et al. 2020). This raised a major concern in the scientific community about the continuous OPFRs exposure and their associated health risks to pregnant women. Given this, Zhang et al. (2016) compared the dietary intake (EDI) of OPFRs with their daily reference dose (RfD) (ng/kg bw/ day). It was found that the estimated exposure levels of most of the compounds were at least two orders of magnitude lower than their RfD values. The human health risk assessment for selected OPFRs was probably underestimated because other possible OPFRs already present in the food were not included. Thus, it is necessary to reduce the human health risks of OPFRs in special populations.
Pregnancy is a critical period for the formation and development of fetal tissues and organs and is potentially vulnerable to environmental toxins for pregnant women and their fetuses. Past studies have shown that mammals are particularly sensitive to organic flame retardants (such as reproductive neurotoxicants and endocrine disruptors) during uterine and neonatal development (Schneider et al. 2014). OPFRs in pregnant women affect the hormone secretion sensitivity, lead to reproductive and immune dysfunction, and affect the reproductive process. OPFRs also have adverse effects on the long-term physical development and neurocognitive development of the fetus (Percy et al. 2020;Gu et al. 2018). Therefore, dietary interventions during pregnancy are very important to ensure the health of the pregnant women and fetus, reduce the human health risks of OPFRs molecular reproductive toxicity, neurotoxicity, and endocrine disruption (ED), and minimize the occurrence of adverse pregnancy outcomes . Hence, to reduce the OPFRs molecular joint toxicity risk in pregnant women, it is recommended to formulate a recommended complementary food scheme (CFS) for pregnant women to create in vivo obstruction in the binding process of OPFRs and joint toxicity receptor (JTR) protein. In this study, molecular docking and 2D-QSAR model were used to screen the molecular parameters that significantly affected the joint toxicity of OPFRs and to reveal the joint toxicity mechanism of OPFRs in pregnant women firstly. Second, MD methods based on the Taguchi experiment design were used to develop recommended CFS that significantly improved the human health risks in pregnant women exposed to OPFRs. Finally, MD and DFT methods were used to quantify the weight of neurotoxicity, reproductive toxicity, and ED toxicity of OPFRs in pregnant women. To the best of our knowledge, no study on the CFS to improve or reduce human health risks of OPFRs in pregnant women has been reported in the literature.

Comprehensive Evaluation Method for OPFRs Molecular Joint Toxicity-Molecular Docking and MD Method
The receptor proteins of the human neurotoxicity (glucocorticoid receptor, GR) Frank et al. 2018), reproductive toxicity (peroxisome proliferator-activated receptor, PPAR) (Kojima et al. 2013;Kuwabara et al. 2012), and ED toxicity (estrogen receptor, ER) (Souza et al. 2017) were obtained from Protein Data Bank (PDB, http:// www. rcsb. org/ pdb) and their PDB IDs were 6CFN, 3VI8, and 5TOA, respectively. ZDOCK module in Discovery Studio 4.0 (DS, Accelrys, Inc., San Diego, American) software was used to dock 20 OPFRs molecules with three toxic receptor protein docking complexes in turn. Possible binding sites in the toxic receptor protein docking complexes were found using the "Find Sites from Receptor Cavities" function in the Define module. Binding sites were subsequently modified and defined based on the incorporation of ligand compounds into the protein binding cavity and docking with the ligand protein. The change in binding capacity was determined from the LibDock scores. The binding free energy ( ΔG bind ) characterized joint toxicity of OPFRs reproductive toxicity, neurotoxicity, and ED toxicity in special populations. MD method was used to calculate the ΔG bind of 20 OPFRs molecules and the complex of JTR protein. The calculation was mainly based on the MD simulation module of Gromacs software in the Dell PowerEdge R7425 server. The complex of OPFR molecules and JTR protein was placed in a cube with a side length of 8.3 nm. The GROMOS96 43A1 force field was used for molecular constraint, and Na + was added to neutralize the system charge (Daura et al. 1998;Hansson et al. 1997). The above composite system was set as a group, and the energy minimization simulations based on the steepest gradient method were performed with the simulated steps set to 100,000. The heat bath and pressure simulation time of the composite system was set to 100 ps with a constant standard atmospheric pressure of 1 bar and the dynamic simulation calculation time of each level was equal to 200 ps. It was reported that as the ΔG bind value of the OPFRs molecules and the JTR proteins complex decreased, the binding affinity of the complex was increased, and as the binding affinity of the complex decreased, the joint effects of reproductive toxicity, neurotoxicity, and ED toxicity of OPFRs were decreased. (Yang et al. 2020a). The calculated values of ΔG bind of OPFRs molecule and JTR protein complex are shown in Table 1. The absolute value of ΔG bind was used to characterize the joint toxicity of OPFRs molecules in the special population. With the increase in the absolute value of ΔG bind , the joint effects of reproductive toxicity, neurotoxicity, and ED toxicity of OPFRs were increased (Yang et al. 2020a).

Analysis of the Joint Toxicity Mechanism of OPFRs Molecules by 2D-QASR Model Assisted by Density Functional Theory (DFT)
We used Gaussian 09 to optimize the structures based on the ground state of DFT at B3LYP level and 6-311 G (d,p) basic set and Chemdraw12.0 software to calculate the seven geometrical parameters, five electronic parameters, eight physicochemical parameters, and seven spectral parameters of the OPFRs molecules. The calculated parameters of 20 OPFRs molecules are shown in Table S1 (Chen et al. 2020;Shi et al. 2015). These parameters were selected as the independent variables, while the joint toxic evaluation value of reproductive toxicity, neurotoxicity, and ED toxicity (the ΔG bind of complex of OPFRs molecules and JTP protein) were the dependent variables (Du et al. 2019;Li et al. 2020). Based on these data, we screened out molecular descriptors using the stepwise regression method in SPSS Statistics 2.0 software. A 2D-QASR model of OPFRs molecular reproductive toxicity, neurotoxicity, and ED toxicity was constructed with 15 (≥ 5) molecules (out of 20 molecules) as the training set and five molecules (≥ 5) as the test set (Golbraikh and Tropsha 2002). The sensitivity analysis method was used to analyze the parameters that significantly affected the joint toxicity of OPFRs.
The 2D-QSAR model should have a clear range of applications. The prediction of compounds using the 2D-QSAR model is reliable when the model is within its application domain. On the contrary, it is less reliable when the model is outside its application domain (Du et al. 2019). The Williams diagram was drawn based on standardized residual δ values and the lever h i values (Yang et al. 2020b). This can visually reflect the scope of model application, judge the reliability and robustness of the 2D-QSAR model, and identify the outliers and influential compounds. The formula for the standardized residual δ and lever value h i used is as follows: where Y exp and Y pre represent the experimental and the predicted values of the target molecules, respectively, n represents the number of target molecules, p is the number of model molecular descriptors, X i represents the matrix of i th target compound descriptor, X is the matrix of all target compound descriptors, X T and X T i represent the transpose matrices of X and X i , respectively. The warning lever h * is the limit value of X and the target molecules with h i < h * are considered as the normal value of the structure. On the contrary, it is regarded as the abnormal value of the structure when h i < h * , where h * is defined as: where n represents the number of compounds, and p represents the number of the model molecular descriptors.
The parameters that significantly affected the joint toxicity of OPFR molecules in the special population were determined by calculating the sensitivity parameters (Du et al. 2019). The sensitivity coefficient (SC i ) is the ratio of the relative change of the predicted value to the relative change of the input parameter and was calculated according to the following Eq. (4), where SC i is the sensitivity coefficients of the i th parameter, ΔX i X i represents the rate of change of the i th parameter, and ΔY i Y i is the rate of change of the joint toxicity of OPFR molecules.

Recommended CFS to Alleviate the Joint Toxicity of OPFRs in Pregnant Women-MD Method Based on Taguchi Experiment Design
In this study, we selected 15 regular high protein-rich food items (milk, eggs, and soybean milk), fresh fruits and vegetables (orange, carrot, broccoli, spinach, grapes, jujube, apple, and kiwi fruit), grains (black rice and oat), and drinks (honey and cubilose) that are recommended to pregnant women. The main nutrients present in each food item were considered as external conditions, and their influence on the binding affinity of OPFRs molecules (TCEP, TPHP, and TDCPP) to JTR protein was measured by MD simulations. The MD method assisted by L 32 (2 8 ) Taguchi experiment design was used to further screen the CFS that would be beneficial in reducing the human health risks of OPFRs in pregnant women. The L 32 (2 8 ) Taguchi experimental design was conducted by selecting eight complementary food items (milk (A), egg (B), orange (C), spinach (D), grape (E), jujube (F), apple (G), and kiwi fruit (H)) that effectively reduced the OPFRs molecular joint toxicity (the reduction rate was between 10.74-44.20%) as the variables to generate the orthogonal experimental method. Taguchi experimental design is a special orthogonal experimental method in which fewer experiments are required to analyze a large number of variables (Castorena-Cortés et al. 2009). OPFR molecules and JTR protein complex were placed in a water cube box with a side length of 15 nm, and the nutrients of the complementary food items were added to the constructed water cube box. The MD simulations of the joint toxicity of OPFRs molecules under different CFS were carried out according to the 32 orthogonal experimental groups generated (Table S2). Factorial analysis was then used to determine the CFS that could minimize the joint toxicity of OPFRs in the special population according to the binding free energy.

Joint Toxicity Mechanism of OPFRs Molecules Based on 2D-QSAR Model
The various parameters of OPFRs molecules used to construct the 2D-QSAR model are shown in Table S1. The joint toxicity of OPFRs molecular reproductive toxicity, neurotoxicity and ED toxicity were taken as the dependent variable, while the molecular structure parameters were considered to be the independent variables. The independent variables that significantly affected the joint toxicity of OPFRs were screened by stepwise regression method in SPSS, and the 2D-QASR model between the OPFRs molecular joint toxicity and its parameters was constructed. According to the 2D-QASR model, we studied the influence mechanism of the parameters on the OPFRs molecular joint toxicity. Taking Considering the value T (the absolute value of binding free energy between OPFRs molecules and the JTR complex) of reproductive toxicity, neurotoxicity, and ED joint toxicity as the dependent variable, the 2D-QSAR model equation for the joint toxicity of OPFR molecules was as follows: The R-value of the 2D-QSAR model was 0.85 (> 0.8), and the Sig value was 0 (< 0.05). Sig represents the significance at P = 0.05. The 2D-QSAR model passed the significance test, indicating that the parameters selected in the equation of the 2D-QSAR model were related to the joint toxicity of the OPFR molecules. The R adj value of the 2D-QSAR model was 0.76 (> 0.6), R adj represents the determination coefficient of the 2D-QSAR model after the calibration of the internal validation parameters and the RMSE TR representing the root mean square error value met all the desired requirements. Thus, the 2D-QSAR model showed a good fitting degree (Golbraikh and Tropsha 2002). The Q 2 LOO value of the 2D-QSAR model was 0.85 (> 0.6) and represented the crossvalidation coefficient, indicating that the model had good stability (Qin et al. 2013). The calculated R 2 test value of the 2D-QSAR model was 0.72 (0.6), suggesting that the model had a strong external prediction ability (Roy et al. 2009).
The linear fitting diagrams of the OPFRs molecular predicted values and experimental values of the training set, the test set, and William's diagram within the application domain of the 2D-QSAR model equation are shown in Fig. 1. The predicted and the experimental values of the training and test sets in the 2D-QSAR model were in good agreement with each other. The scope of application was determined by the lever alert values h * and δ. The OPFRs molecules (15 molecules as training set and five molecules as the test set) used for constructing the 2D-QSAR model were introduced into Eq. (1) to calculate h * and δ. The h * value was 0.75 (X-coordinate) and the δ value was ± 3 (Y-coordinate) (Yang et al. 2020b). All the results of the 2D-QSAR model were within the acceptable range, indicating that the 2D-QSAR model was robust and had a wide range of applications. The model was also reliable in predicting the joint toxicity of similar molecules.
The geometric parameter-quadrupole moment Q YZ , electronic parameter-minimum orbital energy E LUMO , energy gap (EG), and spectrogram parameter-Raman C-O stretching vibration frequency Raman − C − O SVF of OPFRs molecules played important roles in the joint toxic effects. To further illustrate the influence of these parameters on the joint toxicity of OPFRs molecules, a sensitivity analysis was conducted for the above four parameters by using Eq. (5). Each parameter was increased by 10%, 20%, 30%, 40%, and 50%. The influence trend and significance of each parameter on the joint toxicity were expressed by relative sensitivity, and the sensitivity coefficient value of each parameter calculated by the change of parameters is shown in Table 2. By comparing the absolute values of the sensitivity coefficients, it was found that the effect of the variable parameters on the joint toxicity of OPFRs molecules acted in the following trend: E LUMO appeared as the important parameter that significantly influenced the joint toxicity of OPFRs molecules and had a positive coefficient value (in Eq. 5), indicating that it positively affected the toxicity of OPFRs. This implied that as the value of E LUMO was increased, the absolute value of the binding free energy of OPFRs molecules and JTR complex also increased subsequently, thereby increasing the joint toxicity of OPFRs molecules. Moreover, E LUMO represented the electrophilicity of compounds , in other words, the electron affinity of OPFRs molecules is considerably affected the joint toxic effect of OPFRs on the JTR proteins in the special population. OPFRs are electrophilic in nature and thus have strong binding affinities for the biomolecules. During ligand binding to the receptor, E LUMO represents the electric field properties around the binding sites of OPFRs molecules and toxic receptor proteins. Therefore, a change in the electric field around the binding sites can be considered to alter the binding affinities of OPFRs and toxic receptor proteins; thus, reducing the OPFRs molecular toxicity. Previous studies have shown that target receptor proteins involved in various biological activities have high specificity and selectivity for ligands with specific structures or molecular force fields (Jones et al. 2021;Gong et al. 2019). Therefore, we considered adding other ligands into the compound system of OPFRs molecules and JTR complex, and by controlling the field properties around the binding sites of OPFRs and toxic receptor proteins, the competitive binding between multiple ligands was regulated; thus, reducing the binding affinity between OPFRs molecules and JTR proteins, hence the joint toxicity of OPFRs.

Screening of CFSs Based on MD Simulations to Regulate Human Health Risks in Pregnant Women
Dietary intake is a potential pathway of human exposure to OPFRs. Zhang et al. investigated the concentrations of typical OPFRs in 50 varieties of rice and 75 varieties of common food in China and found that the concentration of OPFRs was the highest in rice (Zhang et al. 2016). The intake of meat and fish may be related to higher DPP and BDCPP levels (Thomas et al. 2017). Therefore, in this study, 15 complementary food items with low OPFRs content and regular diet suggested during pregnancy were selected to develop a recommended CFS. The main nutrients in complementary food items were added into the complex system of OPFRs and JTR protein and the competitive combination of nutrients and OPFRs molecules was used to restrict or reduce the human health risks of OPFRs. Various supplements have been used in the past to improve the toxicity risks of biological receptors exposed to contaminants (Zhang et al. 2018;Rajabiesterabadi et al. 2020;Yilmaz 2020). The binding of OPFRs molecules to JTR in the special population was a key step in its joint neurotoxicity, reproductive toxicity, and ED toxicity Sun et al. 2016;Baldwin et al. 2017;Gu et al. 2018). TCEP, TPHP, and TDCPP were the most important OPFRs contaminants detected in the environment and pregnant women. MD simulations were used to simulate the influence of the binding ability of major nutrients in each complementary food items and OPFRs to JTR protein. The recommended CFS beneficial to reduce the health risks in the special population exposed to OPFRs were screened out.
The main nutrients present in milk are β-lactoglobulin, α-lactoalbumin, and milk fat globule membrane (MFGM) protein; in the egg are ovalbumin (OVA) and ovotransferrin (OVT); in the soybean milk are plant protein, vitamin B1, vitamin B3, and niacin; in the orange are vitamin C and lycopene; in carrots are vitamin A and folic acid; in the broccoli are vitamins A, B, C, and K; in the spinach are vitamin A and B; in grapes are vitamin A and B and anthocyanin; in the jujube are vitamin A, B, C, and E; in the apple are vitamin C, apple polyphenols, malic acid, anthocyanin, lycopene, and niacin; in the kiwi fruit are vitamin A, C, and E; in black rice in the cereal category are vitamin A and B and folic acid, in oats are vitamin B1, B2, and E and folic acid. The main nutrients found in honey as a drink are vitamin B and C, and in the cubilose are sialic acid and epidermal growth factor. The binding free energies of OPFRs molecules and JTR protein complexes after the addition of main nutrients found in complementary food items were as follows:  of OPFRs molecules and the JTR protein complexes in the presence of 11 complementary food items (milk, egg, soybean milk, orange, spinach, broccoli, grapes, jujube, apple, kiwi fruit, and cubilose) (Fig. 2). However, the opposite was true in the case of milk (A), egg (B), orange (C), spinach (D), grape (E), jujube (F), apple (G), and kiwi fruit (H) where the absolute values of binding energies of OPFRs molecules and JTR protein complex decreased significantly within the range of 10.74-44.20%. Therefore, it was suggested that the above eight complementary food items could significantly reduce the joint toxicity of OPFRs molecules (TCEP, TPHP, and TDCPP) in pregnant women.
In the above eight complementary food items, except milk and eggs, the rest of them were (fruits and vegetables) edible plant supplements. Previous studies showed that medicinal plants and their extracts had antioxidant effects when varieties of fish were exposed to toxic substances such as heavy metals Cu and Zn, bisphenol A, etc.   helped biological receptors exposed to a variety of toxic substances to resist their biotoxicity. The results of these studies supported our observations and were consistent with the efficacy of plant supplements used in this study.

Recommendation of CFS Based on Taguchi Experimental Design to Improve Human Health Risks in Special Population
MD results showed that milk (A), egg (B), orange (C), spinach (D), grape (E), jujube (F), apple (G), and kiwi fruit (H) could significantly reduce the joint toxicity of OPFRs molecules (TCEP, TPHP, and TDCPP) in special population. The above eight supplemental foods were considered as the experimental variables, and the binding free energies of TCEP, TPHP, and TDCPP molecules and the JTR protein complex were calculated by MD simulations according to the complementary food combinations generated by the L 32 (2 8 ) Taguchi orthogonal experimental design (Fig. 3).
In Fig. 3, the dark blue color represented the presence of complementary food whereas, the light blue color represented the absence. It was observed that the heat map of scheme 33 (blank group, without complementary food items) had the darkest color. The absolute values of binding free energies of 32 groups of CFS were lower than that of scheme 33, indicating that all the CFS were able to reduce the binding ability of OPFRs and the JTR complex protein. The absolute values of the binding free energies of 32 groups of CFS were changed in the range of −14.16% to −85.19%, in comparison to the blank group. When milk (A), orange (C), grape (E), jujube (F), and kiwi fruit (H) were all added at the same time (the main nutrition: β-lactoglobulin, α-lactoalbumin, MFGM protein, vitamin A, B, C, D, and E, anthocyanin, and lycopene), the binding ability of TCEP, TPHP, and TDCPP to JTR protein in special population was significantly reduced (85.19% lower than that in the blank group); consequently, the health risks of OPFRs to humans was significantly reduced.
The absolute values of binding energies of TCEP, TPHP, and TDCPP with the JTR protein in the special populations were taken as the response value and the factorial analysis of Taguchi experimental design was performed subsequently. As the absolute value of average binding free energy was decreased, consequently, the binding affinity between OPFRs molecules and the JTR protein and the toxicity of the OPFRs in special population was also decreased. Hence, according to Table 3, the recommended CFS was: milk (A), orange (C), spinach (D), and jujube (F), their main nutrients were β-lactoglobulin, α-lactoalbumin, MFGM protein, Vitamin A, B, C, D, E, lycopene). The binding free energy of TCEP, TPHP, and TDCPP with the JTR protein (−8.563 kJ/mol) was calculated in the presence of recommended CFS. When compared with the binding free energy of the blank group, the change rate of the binding free energy under recommended CFS was −91.18%. It was the lowest absolute value of binding energy among all CFS and was consistent with the results of the Taguchi experimental design factor analysis. The main nutrition of CFS could reduce the joint toxicity of OPFRs to the maximum extent (85.19%) in 32 groups of Taguchi experimental design schemes [comprised of milk (A), orange (C), grape (E), jujube (F), and kiwi fruit (H)] and the recommended CFS containing milk (A), orange (C), spinach (D) and jujube (F) (sources of main milk protein, vitamins, and plant pigments). These results indicated that β-lactoglobulin, α-lactoalbumin, MFGM protein, vitamins A, B, C, D, and E, and plant pigments were the main factors that significantly reduced the binding ability of OPFRs to the JTR protein in the special population.
Yousef et al. evaluated the immune stimulation, and anti-inflammatory effects of Roselle added to the diet under normal conditions and ammonia exposure on rainbow trout (Oncorhynchus mykiss). Roselle contains a lot of natural pigments such as multi-vitamins and anthocyanins. The results showed that Roselle supplementation significantly increased the white blood cells, plasma total immunoglobulins, alternative complement pathway (ACH50), bactericidal activity, and skin mucus/plasma lysozyme activity. It was demonstrated that Roselle was capable of augmenting immune response and mitigate inflammation in rainbow trout, leading to better health following ammonia toxicity (Yousefi et al. 2021). Supplementation of vitamin C in feed could improve iron toxicity in aquatic Ictalurus punctatus (Yadav et al. 2020). The antioxidant protective effects of multi-vitamins have been shown to counteract the toxicity of hydrophilic metals such as lead (Pb) in many species of fish (Nourian et al. 2019;Shahsavani et al. 2017). Besides, Harsij et al. (2020) found that fish fed with diets containing nanoselenium, vitamins C and E, and antioxidants and exposed to sublethal ammonia showed significantly better growth performance, immune and antioxidant responses than those in the control group. All of these studies have shown that the dietary plants rich in vitamins and natural pigments can reduce the toxicity of the recipient organisms exposed to a variety of toxic materials. The nutrients screened in this study that can significantly reduce the joint toxicity of OPFRs in special population were also vitamins and plant pigments, which were consistent with the above Fig. 3 Heat map of binding free energies of OPFRs molecules-JTR protein complexes with different CFS research results. Therefore, it can be inferred that vitamin A, B, C, D, E and plant pigments (anthocyanin, lycopene) screened in this study were likely to be the main factors that significantly reduce the binding ability of OPFRs to the JTR protein in the special population.
In accordance with the response value (absolute value of binding free energy), the CFS for reducing the binding affinity of OPFRs molecules and JTR followed the sequence: egg (B) > orange (C) > grape (E) > apple (G) > milk (A) > fresh jujube (F) > spinach (D) > kiwi fruit (H). The average binding free energy analysis showed that out of the top three complementary food items in the above sequence, i.e., egg (B), orange (C), and grape (E), the addition of orange only could reduce the binding affinity between OPFRs molecules and the JTR. In 32 groups of the Taguchi experiment scheme, the schemes 10, 12, 26, and 28 met the above-mentioned requirements, where CFS contained only orange (C) and no egg (B) and grape (E). The absolute value of binding free energy of OPFRs molecules and JTR in the presence of the above-mentioned complementary foods was in the range of 58.84%-78.58%, lower than that of the blank group. On the contrary, the schemes 5, 7, 21, and 23 [with egg (B) and grape (E) and without oranges (C)] had a lower absolute value of binding free energy (18.52-34.72%) of OPFRs molecules and JTR than that of the blank group. The reduction degree of the absolute value of binding free energy was much less than that of the CFS determined according to the average response values and the rank of the factorial analysis results, indicating that the factorial analysis results of the Taguchi experimental design were representative and reliable. Therefore, according to the results of factorial analysis of Taguchi experimental design, the recommended CFS with the milk (A), orange (C), spinach (D), and jujube (F) (main nutrient: β-lactoglobulin, α-lactoalbumin, MFGM protein, and vitamins) significantly reduced the binding free ability of OPFRs to the toxic receptor protein in the special population. The main proteins present in the milk and vitamins were the recommended complementary nutrient combinations for the prevention of health risks in special populations exposed to OPFRs. However, with the presence of milk (A), eggs (B), spinach (D), and apple (G) (primary nutrients: β-lactoglobulin, α-lactoalbumin, MFGM protein, OVA, OVT, vitamin, plant pigments, apple polyphenols, and malic acid) at the same time, the absolute value of binding free energy of TCEP, TPHP, and TDCPP with JTR protein was minimally decreased (only 14.16%). Therefore, it was not recommended to consume OVA, OVT, β-lactoglobulin, α-lactoglobulin, and MFGM protein together with the complementary foods rich in vitamins, plant pigments, apple polyphenols, and malic acid. In addition, although no dietary factors were significantly associated with OPFRs and their metabolite concentrations, some studies showed that increased dairy product and fresh intake could be associated with the lower levels of OPFRs metabolites DPP, BDCPP, and ip-DPP (M.B. Thomas et al. 2017) and supported the introduction of the CFS (with milk, orange, spinach and jujube present at the same time) recommended in the present work capable of significantly reducing the binding affinity of OPFRs with the JTR protein.

Mechanism Analysis for the Improvement of OPFRs Molecular Neurotoxicity, Reproductive Toxicity, and ED Toxicity Under Recommended CFS
To investigate the toxic effects that played a major role in joint toxicity in pregnant women exposed to OPFRs, three OPFRs molecules (TCEP, TPHP, and TDCPP) were docked with the neurotoxic receptor, reproductive toxic receptor, and endocrine disrupting toxic receptors, respectively. Then, each molecule of TCEP, TPHP and TDCPP was docked with neurotoxic receptor, reproductive toxic receptor and endocrine disrupting toxic receptor in sequence. The binding free energies of 12 groups of docking complexes were calculated by MD simulations, and their absolute values are shown in Fig. 4. The darkest shades of red, green, and blue colors in the outermost part of the annular histogram in Fig. 4 represented the absolute values of binding free energies of the three OPFRs molecules with the JTR complex. The absolute value of binding energy represented the binding affinity of OPFRs and JTR protein. The colors varying from dark to light in the legend represented the binding affinity of different ligands and each toxic receptor in the annular histogram. Fig. 4 Annular histogram of binding affinity of OPFRs molecules (TCEP, TPHP, and TDCPP) simultaneously or in sequence binding with the neurotoxic receptor, reproductive toxic receptor, and endocrine disrupting toxic receptor It was observed from Fig. 4 that the absolute value of binding free energy of three OPFRs molecules with the JTR complex, neurotoxic receptor, reproductive toxic receptor, and endocrine disrupting toxic receptor were 97. 089, 98.352, 116.27, and 89.089, respectively. The toxicity increases with the increase in absolute values of binding free energies. The results showed that the toxicity effect of OPFRs molecules and reproductive toxic receptor complex was most prominent and increased by 19.76% in comparison to the joint toxicity effect of OPFRs molecules and JTR complex. While the toxicity of TCEP, TPHP, and TDCPP with neurotoxic receptors increased slightly (1.30%), it decreased slightly with ED toxic receptor (8.24%), suggesting that the reproductive toxicity of OPFRs molecules was greater than neurotoxic and endocrine disrupting toxicity and played a dominant role in the joint toxicity. The results showed that the absolute values of binding free energies of TCEP and TPHP with reproductive toxicity receptor (67.876-60.612) were higher than those with neurotoxicity receptor (56.238-48.971) and endocrine disrupting toxicity receptor (48.502-42.640). It can be said that the reproductive toxicity of TCEP and TPHP molecules was greater than neurotoxicity and endocrine toxicity. On the contrary, the TDCPP molecule showed significant neurotoxicity in comparison to TCEP and TPHP, which could be related to the strong neurotoxicity of the TDCPP molecule itself (Ali et al. 2016;Wang et al. 2015;Dishaw et al. 2011). The binding free energy results of the three OPFRs molecules in simultaneous or sequential binding with neurotoxicity receptor, reproductive toxicity receptor, and endocrine toxicity receptor showed that the effect of reproductive toxicity was significantly greater than that of the other two toxic effects and was the main influencing factor in the joint toxicity. Therefore, it was necessary to explore the regulatory mechanism of the joint multi-toxicity in terms of reproductive toxicity.
The recommended CFS proposed in this study could significantly reduce the joint toxicity of OPFRs in the special population (decreased by 91.18%). To investigate whether the recommended CFS also had a significant low health risk effect on the main toxicity (reproductive toxicity) of the joint toxic effects, the binding free energies of OPFRs with the reproductive toxic receptor, neurotoxicity receptor, and endocrine toxic receptor were calculated under the recommended CFS. The results showed that the absolute values of binding free energies of three OPFRs with reproductive toxicity receptor, neurotoxicity receptor, and endocrine toxicity receptor were 20.37, 52.79, and 47.49. However, without recommended CFS, the binding affinities of OPFRs to reproductive toxicity receptor (116.27), neurotoxicity receptor (98.53), and endocrine disrupting toxic receptor (89.09) were decreased by 82.48%, 46.55%, and 46.69%, respectively. The ratio of reduction in reproductive toxicity, neurotoxicity, and endocrine toxicity was 4.69:2.65:2.66. The results showed that the recommended CFS could not only significantly reduce the joint toxicity of OPFRs, but also effectively reduce the single effects of reproductive toxicity, neurotoxicity, and endocrine toxicity, and it was worth mentioning that the reproductive toxicity in the joint toxicity decreased the most. Effects of reproductive toxicity on animals mainly included decreased pregnancy rate, miscarriage, stillbirth, teratogenesis, and fetal developmental disorders (Guerby et al. 2020). Therefore, reproductive toxicity could be considered as a major health risk factor in pregnant women during delivery. The presence of the recommended CFS reduced reproductive toxicity by 176.98% and 176.32% when compared with the reduction in neurotoxicity and endocrine toxicity, respectively. This indicated that the recommended CFS had a significant positive effect in the improvement of reproductive toxicity as the major health risk in pregnant women.

Mechanism Analysis of the Competitive Binding of Nutrients and OPFRs Molecules with JTR Complex
Since main nutrients such as β-lactoglobulin, α-lactoalbumin, MFGM protein, vitamin A, B, C, D, and E, and plant pigments (lycopene and anthocyanin) in CFS significantly reduced the binding ability of OPFRs to JTR, hence, the competitive binding mechanism of main nutrients and OPFRs molecules with the JTR complex was determined by molecular docking and DFT. Varieties of hydrophobic organic contaminants (HOCs) were concentrated in the zucchini plants and contaminated the agricultural soil through the roots of the plants. Major latex-like proteins (MLPs) in the zucchini family played an important role in binding with HOCs. Fujita et al. inhibited the binding of HOCs and MLPs by introducing the compounds with indoles and quinazolines-like structures to facilitate easy binding to MLP. When zucchini plants were cultivated in the contaminated soil with 1.25 mmol/kg pyrene and 12.5 mmol/kg dieldrin, the concentration of pyrene and dieldrin in xylem sap was significantly reduced by 30% and 15%, respectively. It was demonstrated that the pesticides bonded to MLPs, competitively inhibited the binding of MLPs to pyrene and dieldrin in roots, resulting in the reduction of agricultural soil pollution caused by HOCs (Fujita et al. 2020). PUM2 could facilitate the stemness of breast cancer cells by competitively binding to neuropilin-1 3' UTR with miR-376a . The above studies have shown that the competitive binding between different ligands could achieve the optimal regulation of ligand-receptor binding and accomplish the established goals. Therefore, it can be speculated that when ligands were added to the OPFRs-JTR complex system, nutrients modified the nature of the molecular force field around the target binding spot of ligand molecules (including nutrients and OPFRs molecules) and JTR. Consequently, nutrients competitively inhibited the binding of OPFRs to JTR in special populations and changed the binding affinity of OPFRs and JTR, resulting in a subsequent reduction in the human health risks in special populations exposed to OPFRs.
Firstly, nutrients that significantly reduced the binding affinity between OPFRs and JTR protein and the corresponding OPFRs molecules were sequentially docked with JTR by the molecular docking method. The total score (TS) of ligand molecules docked with receptor protein represented their binding affinities. As the TS increases, the binding affinity of ligands and receptors increases. Difference from the molecular docking used in comprehensive evaluation, the purpose of the molecular docking methods in competitive binding was to reflect the binding affinity between OPFRs and toxic proteins, and the binding affinity between nutrients and toxic proteins. In Fig. 5, the TS of vitamin A, B, C, D, and E, lycopene, and anthocyanin docked with JTR complex was the Z value of spherical coordinate with a radius equal to one. (The TS of β-lactoglobulin, α-lactoalbumin, and MFGM protein docked with JTR was not considered, because β-lactoglobulin, α-lactoalbumin, and MFGM protein were directly docked with the JTR protein as a protein complex; thus TS could not be calculated.) The TS of three OPFRs molecules (TCEP, TPHP, and TDCPP) docked with JTR complex was the Z value of spherical coordinate, and the radius was equal to 1.5. The TS of the nutrients docked with JTR slightly fluctuated in the range of -25.63-16.48%, in comparison to the TS score of the three OPFRs molecules. The TS of β-lactoglobulin, α-lactoalbumin, MFGM protein, vitamin A, B, C, D, and E, lycopene, anthocyanin, and the three OPFRs molecules docked with JTR was 81.01% higher than that of the only three OPFRs molecules docked with the JTR, indicating that when β-lactoglobulin, α-lactoalbumin, MFGM protein, vitamin A, B, C, D, and E, lycopene, anthocyanin were present along with the three OPFRs, the binding affinity of nutrients to the JTR protein was better than that of only OPFRs. It can be speculated that OPFRs molecules and the nutrients had competitive binding with the JTR protein. Furthermore, the coexistence of these ten nutrients and OPFRs molecules facilitated the easier binding of nutrients to JTR and occupancy of more target binding spots of JTR by nutrients and reduced or prevented the binding of OPFRs molecules to JTR protein, thereby reducing the joint toxicity of OPFRs in special population.
In the docking system of OPFRs molecules and JTR protein, OPFRs molecules were taken as the object of study. According to the force field diagram of OPFRs and nutrients docked with JTR ( Figure S1), the main non-bonded interactions acted between OPFRs molecules and amino acid residues of the JTR protein were hydrogen bonding and charged or polar interactions, which played an important role in the stability of the ligand-receptor protein complex. While, in the presence of nutrients of the CFS when the nutrients and OPFRs were simultaneously docked with the JTR protein, van der Waals forces with low binding ability was the major nonbonded interaction between OPFRs molecules and the amino acid residues of the JTR protein. This indicated that the presence of nutrients changed the nature of the force field of the target binding spot of the OPFRs-JTR docking system (from hydrogen bond and electrostatic force with the strong binding ability to van der Waals force with small binding ability) (Jiang  and Li 2016;Qu et al. 2012). From the perspective of the total number of forces in the OPFRs-JTR docking system, TCEP, TPHP, and TDCPP generated 12, 24, and 18 non-bonded interactions with amino acid residues of JTR, respectively. However, when OPFRs molecules and multiple nutrients of CFS were simultaneously docked with the JTR protein, the number of non-bonded forces generated by TCEP remained the same but were decreased to 12 and 11 in the case of TPHP and TDCPP, respectively. In conclusion, the nutrients of the recommended CFS presented in the OPFRs-JTR complex altered the properties of the amino acid field of the OPFRs molecules.
According to the results of 2D-QSAR model, E LUMO was the main factor that significantly affected the joint toxicity of OPFRs. In other words, the electron affinity of OPFRs was the main factor affecting the toxic effects of OPFRs on JTR in the special population. Gaussian 09 was used to calculate the electronic parameter E LUMO of TCEP, TPHP, and TDCPP molecules and nutrients that could minimize the binding affinity of OPFRs molecules to the JTR. The E LUMO values of TCEP, TPHP, TDCPP, vitamin A, B, C, D, and E, lycopene, and anthocyanin were 1.32, 1.69, 1.61, −0.07, −0.75, −1.14, −0.19, −1.78, −2.45, and −2.09, respectively. The E LUMO of all the nutrients was significantly lower than that of OPFRs molecule. According to DFT, the E LUMO of a compound was related to its ionization potential. Moreover, it also characterized the binding sensitivity of a compound toward an organic nucleophile (JTR protein). JTR proteins preferred to bind to the molecules with low E LUMO (Karelson et al. 1996;Mumit et al. 2020). This explains the easy binding of the nutrients found in dairy products and fresh food to the receptor protein than OPFRs molecules due to the significantly lower value of E LUMO than OPFRs molecules, thus, altering the force field of amino acid around the OPFRs-JTR complex system to a greater extent and decreasing the binding effect of OPFRs to the target binding spot of the JTR. It is a common phenomenon that multiple ligands compete to bind to a receptor (Fujita et al. 2020;Jones et al. 2020). Wang et al. (2017) showed that among several strong nonbonded interactions of amino acids surrounding a ligandreceptor protein docking complex, the electrostatic force was the main force that kept the HSA-Ceviprex complex stable and targeted the receptor protein selectively bound to the ligand. These results were in complete agreement with our results, validating the effectiveness of the 2D-QSAR model, and thus, it was suggested that the joint toxicity of OPFRs in the special population could be controlled by regulating the intake of nutrients in complementary food.

MD Simulations of the OPFRs Molecular Metabolic Process in the Human Body Under the Recommended CFS
In humans and animals, once the activities of certain metabolic enzymes are induced and activated, the metabolism of toxic components is accelerated; thus, reducing the time for toxic components to accumulate in vivo and the toxicity of the pollutants. For example, Aconitum alkaloids were metabolized by the metabolic enzyme CYP450 into single-ester alkaloids with low toxicity (Park et al. 2016). Some metabolic enzymes belong to kinase, which can activate the toxicity of pollutants in vivo. For example, CYP3A65 is an important metabolic enzyme in zebrafish that activates BDE47 to cause developmental toxicity, resulting in delayed embryo hatching and abnormal embryonic neurodevelopment. This could be due to the decreased thyroxine levels and the subsequent destruction of thyroid hormone homeostasis through the up-regulation of Dio3A and Dio3B gene expressions by CYP3A65. However, it was found that CYP3A65 knockout could significantly improve the thyroid hormone homeostasis and reduce the toxic effect of BDE47 on zebrafish embryos (Yang et al. 2017).
Six metabolic kinases [thyroid hormone receptor (TRα 1 ), retinoic acid receptor α (RARα), retinoid X receptor (RXRα), pregnane X receptor (PXR), and liver X receptor (LXR)] in pregnant women were selected to simulate the metabolic activation processes with OPFRs molecules by MD method. The activation processes were simulated with and without the recommended CFS. The OPFRs molecules bonded to the above six metabolic kinases could activate the toxicity of OPFRs in vivo, which in turn could inhibit or block normal biochemical reactions in the human body; thus, inducing hepatomegaly, thyroid hormone secretion disorder, and other toxic effects (Kojima et al. 2013;Belcher et al. 2014;Jiang et al. 2020). Therefore, it can be deduced that as the binding affinity between OPFRs and metabolic kinases decreases, the activation degree of OPFRs molecular toxicity also decreases thus, reducing the health risks of OPFRs in special populations.
MD simulations were carried out to calculate the binding free energy of OPFRs and six metabolic kinases receptors in the presence of the recommended CFS comprised of milk (A), orange (C), spinach (D), and jujube (F), as well as without them (Fig. 6). The absolute value of binding free energies of OPFRs molecules and TRα1, TRβ1, RARα, RXRα, PXR, and LXR in the presence of recommended CFS were decreased by 41.32%, 60.70%, 87.54%, 39.77%, 28.85%, 68.67%, respectively, considering the binding energy values without the recommended CFS as the reference. These results verified that the recommended CFS could not only significantly reduce the joint toxicity of OPFRs in the special population but also reduce the binding affinity between OPFRs and metabolic kinases receptors, inhibiting the toxic effects of OPFRs. Hence, it can be anticipated that the recommended CFS had certain theoretical guiding significance to reduce the human health risk of OPFRs in pregnant women.

Conclusions
In summary, the joint toxicity mechanism of OPFRs molecules in the special population was studied via integrated molecular docking and MD simulations based on the 2D-QSAR model, DFT, and Taguchi experiment design. The nutrients in the recommended CFS could easily occupy the target binding spot of the JTR protein, which reduced or prevented the binding of OPFRs to the JTR. Besides, the health risk regulation scheme of human exposure to OPFRs for the special population was presented, where the recommended CFS could theoretically reduce the joint effects of neurotoxicity, reproductive toxicity, and ED toxicity in the special population. When β-lactoglobulin, α-lactoalbumin, milk fat globule membrane (MFGM) protein, ovalbumin (OVA), ovotransferrin (OVT), vitamin, plant pigment, apple polyphenols, and malic acid were present in the CFS, the joint toxicity of OPFRs in pregnant women were significantly decreased. The methods constructed in this study provided theoretical support for mitigating the potential health risks in the special population exposed to OPFRs in the environment.