Network Pharmacology-Based Prediction and Verification of the Potential Mechanisms of He's Yangchao Formula against Diminished Ovarian Reserve

Background He's Yangchao formula (HSYC) has been clinically proven to be effective in treating diminished ovarian reserve (DOR). However, the underlying molecular mechanisms of HSYC in DOR are unclear. Objective This study aims to predict the underlying mechanisms of He's Yangchao formula (HSYC) against DOR through network pharmacology strategies and verify in vivo. Methods Systematic network pharmacology was used to speculate the bioactive components, potential targets, and the underlying mechanism of HSYC in the treatment of DOR. Then, the CTX-induced DOR mouse model was established to verify the effect of HSYC against DOR and the possible molecular mechanisms as predicted in the network pharmacology approach. Results A total of 44 active components and 423 potential targets were obtained in HSYC. In addition, 91 targets of DOR were also screened. The identified hub genes were AKT1, ESR1, IL6, and P53. Further molecular docking showed that the four hub targets were well-bound with their corresponding compounds. In vivo experiments showed that HSYC could promote the recovery of the estrous cycle and increase the number of primordial, growing follicles and corpora lutea. Besides, The results of qRT-PCR showed HSYC could regulate the expression of AKT1, ESR1, P53, and IL6 in DOR mice. Conclusion It was demonstrated that HSYC could increase ovarian reserves, and AKT1, ESR1, IL6, and P53 may play an essential role in this effect, which provided a new reference for the current lack of active interventions of DOR.


Introduction
Diminished ovarian reserve (DOR) defines as a decrease in the quantity and quantity of oocytes, which is a major cause of female infertility and assisted reproductive technology (ART) failure [1][2][3]. Its occurrence rate has shown rapid growth in recent years. e prevalence has risen from 19% to 26% over the past decades, representing a great challenge in modern infertility management [4,5]. As a complex disease, DOR seriously affects the mental and physical health of women. Currently, the main pathogen of DOR is unclear. It has been reported that factors such as genetic predisposition, autoimmunity, oncologic treatments, cigarette smoking, ovarian surgery, and psychological stress are all related to the onset of DOR [6][7][8].
e decline of ovarian reserve is a continuous and gradual process. e premature depletion of ovarian reserve eventually results in premature ovarian insufficiency, which is a more severe condition that makes it much harder to treat. erefore, rapid treatment is crucial for DOR women. Various treatment strategies are being used to fight against DOR, including estrogen supplementation, dehydroepiandrosterone (DHEA), ovulation induction therapy, ART, and ovum donation. DHEA is an endogenous steroid that originates from the zona reticularis of the adrenal cortex and the ovarian theca cells in women. It is produced by the conversion of cholesterol, acting an important role in the formation of estradiol in peripheral tissues [9]. DHEA is often used in the treatment of DOR. However, long-term use of DHEA can lead to some side effects [10]. erefore, it is an urgent issue to develop a more economical, safe, and comprehensive treatment for DOR.
Traditional Medicine has been included in the influential global medical compendium by the World Health Organization. Traditional Chinese Medicine (TCM) is a branch of Traditional Medicine and is an essential element of China's national health policies. It exerts beneficial therapeutic efficacy and reduces side effects through various herbal combinations, thereby preventing and treating various diseases [11]. As shown in Table 1, He's Yangchao formula (HSYC) comprises eight-herbal medicines. It is a particular herbal prescription for DOR that originates from the ancestral experience of "He's TCM Gynecology," which is the famous school doctrine of Gynecology of TCM in China and has lasted for more than 170 years [12]. HSYC effectively improved the clinical symptoms, increased the AMH level and the antral follicle count, and reduced serum follicle-stimulating hormone (FSH) in DOR patients [13]. However, the mechanism of HSYC in the treatment of DOR remains unclear.
Compared with western medicine, TCM is characterized by multicomponent and multitarget and exerts its effect on multiple biological processes to treat diseases via its diverse bioactive components, which act on multiple targets [14]. Network pharmacology's emergence provides us with new ideas and methods to explore the complex network relationship between multicomponent and multitarget [15]. Network pharmacology is a branch of pharmacology, which takes the theory of systems biology and multiple pharmacology as a foundation and biomolecular network as the main means of identifying the active compounds in herbal medicines and potential therapeutic targets [16].
is study identified the bioactive compounds, potential targets, and the underlying mechanism of HSYC against DOR. In addition, in vivo experiments were performed to validate the molecular mechanism. Figure 1 shows the workflow of this study. is study would provide novel insights into the mechanisms underlying TCM-mediated efficacy in DOR and provide the theoretical basis for further experimental investigations. ). TCMSP is a complete library of Chinese herbal ingredients containing more than 36,000 chemical molecules, ingredients-related targets, and their relationships with diseases [11]. TCM is mostly oral preparation and exerts various biological effects by reaching organs and tissue after absorption, distribution, metabolism, and excretion (ADME) in vivo, which is called the pharmacokinetics of TCM. Oral bioavailability (OB) and drug-like (DL) are the important parameters of ADME [17,18]. e components with OB ≥ 30% and DL ≥ 0.18 were selected as candidate main bioactive components. TCMSP, the Encyclopedia of Traditional Chinese Medicine (ETCM, https://www.tcmip.cn/ ETCM/index.php/Home/Index), SymMap (https://www. symmap.org/ (accessed on 15 March 2020)), and Pub-Chem (https://pubchem.ncbi.nlm.nih.gov (accessed on 15 March 2020)) were used to retrieve and predict the potential targets of bioactive components in HSYC. e obtained potential targets, including the name, gene ID, and organism, were normalized based on the UniProt protein database (https://www.uniprot.org (accessed on 18 March 2020)).

Identification of DOR-Related Targets.
e DOR-related targets were collected from three databases. GeneCards (https://www.genecards.org) is a database for collecting a comprehensive compendium of human genes. e keyword "diminished ovarian reserve" was used to search GeneCards database (accessed on 18 March 2020). Next, "diminished ovarian reserve" was used as a keyword to search DisGeNET (https://www.disgenet.org (accessed on 18 March 2020)). In addition, we searched the National Centre for Biotechnology Information Gene (NCBI Gene, https://www.ncbi.nlm.nih. gov/gene (accessed on 18 March 2020)).

Network Construction.
Based on the above analyses, a visual interaction network based on the herbs-compoundstargets interactions was established by Cytoscape software. Meanwhile, to show the mechanism of action of HSYC, the protein-protein interaction (PPI) network of potential targets of the bioactive compounds in HSYC was constructed by the STRING database (https://string-db.org (accessed on 18 March 2020)) with the cut-off criterion of combined score > 0.9 and the species "Homo sapiens." In addition, the targets of the active compounds in HSYC were matched with the DOR targets to obtain the intersected genes. en, the intersected genes of compounds-DOR were submitted to the STRING database to structure another PPI network.  Evidence-Based Complementary and Alternative Medicine was established by Cytoscape software to holistically explain the mechanism of HSYC in the treatment of DOR.

Screening of Hub Targets and Molecular Docking
Simulation. e intersected genes of the three subnetworks based on three topology parameters, including degree, betweenness, and closeness, were intersected with the 22 potential therapeutic targets to get hub genes. e crystal structures of protein receptors of hub targets were obtained from RCSB Protein Data Bank (https://www.rcsb.org/ (accessed on 29 March 2020)). PyMol 2.4.0 was used to remove the original ligand, solvent molecules, and redundant protein chains and add polar hydrogen for the crystal structures of protein receptors. en AutoDock Tools 1.5.6 was used to compute Gasteiger charges and determine the center and size of the docking box for further docking. e 3D structure of the bioactive compounds was obtained from the PubChem database. e structure was assigned MMFF94 force field and charges, made energy minimization and added polar hydrogen. AutoDock Vina, a freely available open-source program for molecular docking, was used to evaluate the binding of the bioactive compounds and their corresponding targets. Before molecular docking, all the protein structures and the bioactive compounds were converted to PDBQT format for further docking using AutoDock Tools 1.5.6. en the bioactive compounds were docked into their corresponding proteins using AutoDock Vina. Finally, the binding affinity was calculated by Auto-Dock Vina and visualized the docking result by the opensource version PyMol 2.4.0 software.

2.7.
e Construction of DOR Models. irty-six female C57BL/6 mice (8 weeks old) were purchased from Shanghai SIPPR-Bk Lab Animal Company, with the license number SCXK (Shanghai) 2018-0006. All animal experimental procedures were approved by the Institutional Animal Care and Use Committee of the Zhejiang Chinese Medical University (Approval number: IACUC-20180402-02). e mice were maintained under standard conditions (20°C-25°C, natural light cycle, 40%-60% relative humidity) with ad libitum access to rodent chow and water. After 1week adaptation, mice were randomly divided into six groups (n � 6/group): control, DOR, low (Low-HSYC), moderate (Mod-HSYC), high (High-HSYC), and DHEA. A DOR mice model induced by cyclophosphamide (CTX) was established based on the previous study [19]. All mice in the DOR, HSYC, and DHEA groups received an intraperitoneal (ip) injection once of 90 mg/kg CTX (Sigma-Aldrich, Germany), while mice in the control group received an equal volume of saline. e DOR mice in HSYC groups received 5.88 mg/kg, 11.75 mg/kg, and 23.5 mg/kg body weight HSYC via intragastric administration daily. e DOR mice in the DHEA group were treated with 9.75 mg/kg body weight DHEA (Sigma-Aldrich, Germany) via intragastric administration daily. e mouse administration dose was calculated by referring to Experimental Methodology of Pharmacology compiled by Wei et al. [20]. e mice in the control and DOR groups received an equal volume of saline.
It takes approximately three weeks for primordial follicles in mice to pass through the primary, secondary, preantral, and preovulatory stages [21]. erefore, the administration should take more than three weeks to evaluate the short effects. On the 28th day, all animals were sacrificed by cervical dislocation by a skilled expert to minimize suffering. Cervical dislocation is a standard method for euthanasia of experimental mice (weight < 200 g) [22]. e blood samples and ovarian tissues were collected from each animal.

Evaluation of Estrous Cycle by Vaginal Cytology.
Vaginal smears were obtained daily at 09 : 00 am to observe the estrus cycle. Vaginal smear samples were stained with toluidine blue (Qiagen) and analyzed under a light microscope (Olympus, Tokyo, Japan). e estrous cycle was determined according to the cytological appearance as previously reported [23].

Follicular Counts.
e ovarian tissues were fixed in paraformaldehyde, embedded in paraffin, and then cut into 5 μm thick sections. e number of follicles was observed by hematoxylin and eosin (HE) staining according to the methods previously described [24]. All HE sections were examined by two researchers (L.Q.Y. and Y.Z.).

Quantitative Real-Time PCR (qRT-PCR) Analysis.
Total RNA was extracted from ovarian tissues using TRIzol reagent (Invitrogen, Shanghai, China). e first-strand cDNA was reversely transcribed using a Bestar qPCR RT Kit (DBI Bioscience, Ludwigshafen, Germany). 1 µg of DNA sample was used for the PCR amplification of related genes according to the reaction conditions as follows: 95°C for 4 min, 38 cycles of 94°C 30 s, 60°C 30 s, 72°C 30 s. Primers were synthesized by Sangon (Shanghai, China) and listed in Table S1. A DBI Bestar ® SybrGreen qPCR Master Mix kit (DBI Bioscience) and Agilent Stratagene Mx3000P Realtime PCR instrument (DBI Bioscience) were used for the amplification. e relative expression level of each gene was calculated using the 2 −ΔΔCt method. β-Actin was used as the internal reference gene.

Statistical Analysis.
Animal data were expressed as the mean ± standard deviation (SD) from 6 mice. Statistical analysis was conducted in GraphPad Prism 6.0. Differences between groups were analyzed using an unpaired t-test, and differences across six groups were analyzed using one-way ANOVA. A significant level was set at P < 0.05.

Herb-Compound-Target Network.
A total of 44 kinds of bioactive components in HSYC with good ADME properties were obtained. e detailed information on bioactive compounds is described in Table S2. Moreover, there were 423 targets of bioactive compounds obtained from several databases (Table S3). en, the herb-compound-target network was constructed with 411nodes and 1429 edges (see 4 Evidence-Based Complementary and Alternative Medicine Figure 2). In the visual interaction network, the size of the node represents the size of the degree value. e color of the nodes is illustrated from red to cyan in descending order of degree values. It can be seen from Figure 2 that quercetin, 3′methoxydaidzein, and stigmasterol compounds, as well as ESR1, PTGS2, AR, PTGS1, and PGR genes, play a crucial role in the network.

PPI Network of the HSYC's Targets.
e PPI network of targets in HSYC with the highest confidence was constructed using the STRING database, including 288 nodes and 1395 edges (see Figure 3). e color of the nodes is illustrated from red to cyan in descending order of degree values. So MAPK1, P53, AKTI, MAPK3, JUN RELA, HSP90AA1, TNF, and IL6 may be the main targets of biological networks regulated by HSYC.

PPI Network of the Potential erapeutic Targets.
A total of 91 DOR-related targets were obtained from the human genomic database. e number of DOR-related targets in GeneCards, DisGeNET, and NCBI Gene was 44, 43, and 28, respectively (Table S4). Twenty-two common HSYC and DOR targets (potential therapeutic targets) were obtained using the Venn Diagram tool (Figures 4(a) and 4(b)), showing the PPI network of the potential therapeutic targets.

KEGG Pathway Analysis and GO Functional
Enrichment.
e intersected targets of HSYC and DOR were analyzed by the R package. KEGG pathway analysis revealed that the targets of the PPI network were assigned to 89 KEGG pathways based on P < 0.01. GO annotation showed that the targets of the PPI network were classified into 1139 biological process (BP) terms, 20 cellular component (CC) terms, and 43 molecular function (MF) terms according to P < 0.01 (see Figure 5).
As shown in Figure 5(a), the targets of the PPI network in the BP ontology are primarily associated with epithelial cell proliferation, regulation of epithelial cell proliferation, positive regulation of MAPK cascade, protein kinase B signaling, reproductive structure development, reproductive system development, positive regulation of DNA-binding transcription factor activity, gliogenesis, positive regulation of reactive oxygen species metabolic process, and regulation of reactive oxygen species metabolic process. Based on GO annotation of MF, it could be seen that the targets were mainly involved in receptor regulator activity, cytokine receptor binding, growth factor receptor binding, protein phosphatase binding, cytokine activity, steroid hormone receptor activity, antioxidant activity, steroid binding, and nuclear receptor activity. For the CC ontology, the targets were mainly located in the vesicle lumen, endoplasmic reticulum lumen, nuclear chromatin, cytoplasmic vesicle lumen, platelet alpha lumen, platelet alpha, PML body, postsynaptic cytosol, region of cytosol, and the myelin sheath.
As shown in Figure 5(b), KEGG pathway enrichment analysis suggested that the targets were mainly related to signal pathways, such as the PI3K-AKT signaling pathway, MAPK signaling pathway, HIF-1 signaling pathway, ErbB signaling pathway, JAK-STAT signaling pathway, FoxO signaling pathway, estrogen signaling pathway, mTOR signaling pathway, VEGF signaling pathway, and NOD-like receptor signaling pathway. e results indicated that the targets of HSYC were distributed in different metabolic pathways.
e "multicomponent, multitarget, and multipathway" mutual regulation is the possible mechanism for the treatment of DOR.

Screening of Hub Targets and Docking Results Analysis.
Network analyzer, a plug-in of Cytoscape, was used to calculate the node topology parameters of the PPI network of the bioactive compounds' potential targets in HSYC. e top 20 genes screened by degree, betweenness, and closeness were used to construct the subnetworks (Figures 6(a)-6(c)). In our study, four hub targets, including AKT1, ESR1, P53, and IL6, were identified from the Venn map ( Figure 6(d)).
Text, the four hub targets were selected to perform molecular docking with their corresponding compounds. AutoDock Vina evaluated the small molecules' binding strength with proteins mainly through the binding energy (affinity). If the affinity is less than zero, it means that the ligand can spontaneously bind to the receptor. e smaller the value is, the easier the bioactive component is to bind to the receptor. As shown in Table 2, a total of 36 pairs of docking results were obtained. Among the docking results, most binding complexes possessed high binding affinity, which means these compounds may bind to the four hub targets to improve ovarian reserve. Figure 7 shows the binging mode of hub targets with their corresponding compounds with high binding affinity.

Herb-Component-Target-Pathway Interaction Network.
en, a herb-compound-target-pathway network (see Figure 8) was constructed. In the visual interaction network, the squares represent herbs, components, targets, and pathways; the edges represent the interaction of each other. In the network, the four hub genes (AKT1, ESR1, IL6, and P53) regulate the top 10 significant pathways, targeted by an average of 8.5 compounds and all the herbs. Figure 9(a), mice in the DOR groups exhibited significantly lower body weight than mice in the other groups on the 7th, 21st, and 28th day (P < 0.05). In addition, e mice's estrous cycle was continuously measured to evaluate the protective effects of HSYC on ovarian function. e diestrus stage of the DOR group was significantly prolonged during the 28 days of experiments, as presented in Figure 9(b). is abnormality did not completely recover to normal lengths, but it was attenuated in the DHEA and HSYC groups, with a better effect in the High-HSYC group. e normal estrous cycle has four stages, including proestrus, estrus, metestrus, and diestrus stages, as shown in Figure 9(c).

HSYC Improves the Ovarian Structure and Follicle Development in DOR Mice.
e histological analysis of the ovaries is presented in Figures 10(a)-10(d) . Compared with the control group, except for the High-HSYC group, the numbers of primordial follicles, growing follicles (including primary follicles, secondary follicles, and anal follicles), and corpus luteum were significantly reduced in the other groups ( Figure 10(a)-10(d), P < 0.05). is result shows that highdose HSYC can prevent the loss of different types of follicles caused by CTX. Both HSYC and DHEA can significantly increase the number of growing follicles and the number of corpus luteum (P < 0.05), and the effect of HSYC is better than that of DHEA (Figures 10(a), 10(c), and 10(d), P < 0.05).

HSYC Regulates the mRNA Expressions of the Hub
Targets. In contrast to the control group, the mRNA expressions of AKT1 and ESR1 declined in the DOR group (P < 0.05); the expressions then elevated after HSYC and DHEA treatment (Figures 11(a) and 11(b), P < 0.05). In addition, the expressions of IL6 and P53 have the opposite tendency.

Discussion
Nowadays, delaying pregnancy and parenthood has become a global phenomenon, especially in high-resource countries. However, females' ovarian reserve falls sharply after the age of 35. In addition to causing reproductive dysfunction, DOR also increases cardiovascular disease and depression [25].
Clinical findings have confirmed that HSYC is effective in treating DOR, but the underlying mechanism remains unknown. e network pharmacology approach, which integrates systems biology and in silico technologies, may offer a direction for the mechanistic study of complicated TCM [16]. In the current study, we used this approach to clarify the pharmacological mechanism of HSYC on DOR.
In the system of TCM, compounds lacking proper pharmacokinetic properties could not reach the target organs to deliver the biological activities [26]. In the current study, the compounds in HSYC with OB ≥ 30% and DL ≥ 0.18 were considered pharmacokinetically active as they are possibly absorbed and distributed in the human body. In the herb-compound-target-disease network, compounds with high-degree may account for the major Evidence-Based Complementary and Alternative Medicine therapeutic effects of HSYC on DOR. In our study, there were 44 bioactive compounds in HSYC. Quercetin and β-sitosterol have a high degree, which suggested that both of them may have major therapeutic effects. Quercetin has extensive biological properties, including antioxidant, antiinflammatory, and antiapoptosis [27,28]. It has been shown that quercetin significantly increased the primordial follicle number and AMH level and decreased atretic follicle number in CTX-induced DOR mice [29,30]. β-Sitosterol is another significant bioactive compound. e pharmaceutical properties of β-sitosterol, like anti-infection, anti-inflammatory activity, antioxidant, and estrogenic functions, have been reported [31,32]. Mohammadrezaei et al. reported that the treatment with β-sitosterol increased testosterone and ethoxyresorufin-O-deethylase levels in Kutum's fertilized egg [33]. ese findings might indicate the potential therapeutic efficacy of our present eight-herbal medicine in DOR. e KEGG pathways of HSYC against DOR include the PI3K-AKT signaling pathway, MAPK signaling pathway, and HIF-1 signaling pathway. Existing studies reported that PI3K/AKT had been shown to have an important role within DOR. Studies have shown that the PI3K/AKT pathway's activation plays a crucial role in folliculogenesis, including follicle recruitment, development, and maturation [34][35][36]. Besides, suppressing MAPK signaling could inhibit DOR and preserve ovarian function and structure [37,38]. HIF-1 is a master regulator of cellular adaptation to hypoxia and has been suggested as a potential therapeutic target in DOR. HIF-1 contributed to granulosa cell proliferation, which is crucial for ovarian follicle growth by regulating cell proliferation and follicle-stimulating hormone-mediated autophagy [39]. e four hub targets, AKT1, ESR1, P53, and IL6, were created as 3D graphs of docking and analysis. Further molecular docking showed that the four hub targets were well-bound with their corresponding compounds. AKT1 is the important gene in the PI3K-AKT signaling pathway. Also, AKT1 is a multifunctional protein regulating cell growth, survival, and proliferation [40]. It plays a critical role in both the growth and maturation of the oocyte. AKT1 could inhibit follicular cell apoptosis by reducing BAX and FOXO3a in vitro cultures [41]. AKT1−/− females mice display a reduced number of growing antral follicles and abnormal estrous cyclicity [42]. As the primary female sex hormones, estrogens exert their functions by binding to  Evidence-Based Complementary and Alternative Medicine estrogen receptor α (ERα), encoded by the gene ESR1 [43,44]. ESR1 knockout mice showed reproductive impairment as anovulation and complete infertility [45]. Several studies have reported that ESR1 polymorphism is associated with primary ovarian failure (POF) [46,47]. Professor Zi-Jiang Chen's team conducted a genetic analysis on 371 unrelated idiopathic POF patients and 800 normal women (all of Han nationality) and found that ESR1 gene mutations are significantly related to POF [48]. Many studies have also revealed that CTX-induced loss of ovarian reserve involves the activation of P53 [49]. P53 can directly activate the transcription of various proapoptotic genes in CTXinduced ovarian damage mice [49]. Similarly, Stefania et al.
used ChIP experiments to confirm that P53 plays a direct role in CTX-induced apoptosis in mouse oocytes [50]. Besides, IL6 is one of the inflammatory cytokines and plays an irreplaceable role in the inflammatory response. Recent studies have shown that inflammation is a crucial marker of the aging ovarian stroma and is considered a new DOR mechanism [51,52]. Clinical studies in humans have shown that inflammatory marker levels are associated with the risk of DOR. Yue et al. found that compared to the control group, the serum level of IL6 in the primary ovarian failure group was significantly higher [53]. By molecular docking, we found that AKT1, ESR1, P53, and IL6 have the good binding ability with bioactive compounds in HSYC. erefore, based on the molecular docking results and combining these hub targets' physiological roles in the ovaries, we speculate that HSYC could play a beneficial role in DOR mice via AKT1, ESR1 P53, and IL6.
Next, we studied the efficiency and mechanism of HSYC on DOR mice through animal experiments. A regular estrous cycle that lasts for 4 to 6 days in mice can be a key index of the integrity of the HPO reproductive axis [54]. In the present study, HE staining results revealed that HSYC could significantly improve the estrous cycle. ese results indicate a protective role of HSYC on the DOR mice HPO axis. In addition, the primordial follicle pool is nonrenewable. eir development ends with either ovulation or atresia. erefore, it represents the most important follicular [54]. is    P53 and IL6 was upregulated in DOR mice. HSYC could promote the expression of AKT1 and ESR1 and inhibit the expression of P53 and IL6. ese are consistent with previous reports. ese results were consistent with the results predicted by network pharmacology and molecular docking. As with any study, this study did have some limitations. First, due to public databases being updated in real-time, this study can only partially elucidate the molecular mechanism by which HSYC improves ovarian reserve. Second, only four hub genes were verified by molecular docking, and in vitro experiments, other possible pathways should be performed through in vivo and in vitro experiments in the future.

Conclusions
Taken together, HSYC may ameliorate ovarian reserve against CTX-induced DOR via the regulation of multiple targets. Our study further suggested that combined network pharmacology prediction and experimental validation study may offer a valuable tool to characterize the action mechanism of TCM.

Data Availability
e data used to support the findings of this study are included within the article and the supplementary materials.  Figure 11: e mRNA expression of four hub targets. e relative mRNA expression levels of AKT1, ESR1, IL6, and P53 genes, respectively. Data are expressed as mean ± SD. * P < 0.05, compared with the control group; # P < 0.05, compared with the DOR group.  Table S1. e primer sequences used in this present study. Table S2. Information of bioactive compounds in HSYC with good ADME properties. Table S3. Targets of bioactive compounds obtained from databases. Table S4. DOR-related targets. (Supplementary Materials)