Using systems biology and drug repositioning approaches to discover FDA-approved drugs candidates for endometriosis treatment

Endometriosis impacts around 10% (190 million) of women and girls of reproductive age worldwide [ 1 ]. Endometriosis is a multifaceted and systemic clinical disease that can adversely affect women’s reproductive health and quality of life; it may commence with an individual’s first menstrual period and persist until menopause [ 1 ]. This condition’s pathological hallmark is abnormal endometrial tissue growth in the uterine outer regions [ 2 ]. The ovaries and pelvic peritoneum are the predominant locations for the development of endometriotic lesions. Endometriotic lesions may also occur in other locations, including the fallopian tubes, abdominal wall, intestines, cervix, bladder, and vagina [ 3 , 4 ]. A comprehensive understanding of the pathogenesis of endometriosis may have significant clinical and therapeutic ramifications. Endometriosis has a complex pathophysiology that involves several pathways, including genetics, abnormal endocrine signaling, dysregulated cell growth and death, ectopic endometrial tissue, and altered immunity [ 5 ]. Endometriosis is a chronic condition that induces significant pain during menstruation, sexual intercourse, urination, or defecation. Endometriosis symptoms encompass chronic pelvic pain, bloating, nausea, fatigue, and occasionally depression, anxiety, and infertility. Symptoms of endometriosis encompass abdominal distension accompanied by infection, irregular uterine hemorrhage, vaginal discharge, pelvic pain, and discomfort in the lower abdomen. Additional symptoms of endometriosis encompass constipation accompanied by stomach or intestinal pain [ 1 , 6 , 7 ]. The sole definitive method for diagnosing endometriosis requires surgical intervention or laparoscopy, leading to a latency period of 7–11 years from the onset of symptoms to a conclusive diagnosis [ 8 ], and entails an increased risk of disease progression during that time [ 9 ]. In recent years, noninvasive imaging techniques, particularly transvaginal ultrasonography and magnetic resonance imaging (MRI), have improved diagnostic accuracy for endometriosis, enabling the staging and classification of different types of the endometriosis without the need for surgical procedures [ 10 – 12 ]. Drug repurposing denotes the identification of novel therapeutic applications for established medications to improve their efficacy and optimize their use [ 13 ]. The advantages of Food and Drug Administration (FDA)-approved drug repositioning are numerous, and since many preclinical and clinical trials have already been carried out on these drugs, their use saves time and cost as compared to developing de novo drugs [ 14 ]. The primary hurdles involve identifying drug target proteins (receptors) associated with diseases and discovering pharmacological agents (small molecules) that can mitigate these disorders through contact with the target proteins. The identification of hub genes via bioinformatics has become an effective method for discovering prospective biomarkers and therapeutic targets [ 15 ]. Transcriptomics analysis is a prominent strategy for discovering genetic biomarkers. Proteins produced by genomic biomarkers are regarded as key receptors [ 16 – 18 ]. Endometriosis treatment options currently include drugs, surgical procedures, surgery combined with drugs, and assisted reproductive technology (ART). While surgery remains an appropriate method for endometriosis-related pelvic pain, its efficacy is undermined by limitations, including a 40−50% recurrence rate after surgery [ 19 ]. Administering drugs to reduce related symptoms is preferable for patients without surgical indications. First-line pharmacological treatments comprise non-steroidal anti-inflammatory drugs (NSAIDs), progestins, and oral contraceptives (OCs). Second-line drugs comprise gonadotropin-releasing hormone agonists (GnRH-a) and the levonorgestrel intrauterine delivery method. Nonetheless, the administration of the mentioned drugs, including GnRH-a, is associated with significant adverse effects, such as symptoms of estrogen insufficiency and inhibition of ovulation. Consequently, it is inappropriate for prolonged usage, particularly for individuals with reproductive requirements [ 20 ]. Given these limitations, there is an urgent need for identification of new biological targets and drugs. Zhang et al .‘s study indicated that the primary biological targets of tanshinone IIA, the active component of Chinese medicine Danshen ( Salvia miltiorrhiza Bge.), for endometriosis treatment are vascular endothelial growth factor A (VEGFA), matrix metalloproteinase 9 (MMP-9), estrogen receptor-1 (ESR1), intercellular adhesion molecule (ICAM)‑1, and interleukin-2 (IL‑2), which exert pathological effects on the adhesion, invasion, and angiogenesis of ectopic endometrial tissue within the pelvic and abdominal cavities in endometriosis [ 21 ]. The current study also aimed to find hub genes in endometriosis, identify and evaluate FDA-approved drugs that target the proteins encoded by these genes. In the present study, we first downloaded GSE120103 dataset from the Gene Expression Omnibus (GEO). After pre-processing and normalizing the data by the “Limma” package in R Studio, we identified DEGs. Next, Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed by the ShinyGO 0.81 and DAVID tools, respectively. PPI network was constructed using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING), and Cytoscape software was used to identify cluster modules and hub genes. Then, target transcription factors (TFs) and microRNAs (miRNAs) of selected hub genes were predicted by NetworkAnalyst tool. Subsequently, the druggability of hub genes and their associated TFs was assessed using the Drug-Gene Interaction Database (DGIdb 5.0). The DrugBank database was searched to find FDA-approved drugs that were ligands for the target protein, and then their virtual screening was performed using PyRx 0.8 software. Finally, molecular dynamics (MD) simulation of the target protein and docked complexes were performed by the Assisted Model Building with Energy Refinement (AMBER) 18 program.

All four groups of the GSE120103 dataset underwent normalization, and a box plot was generated to illustrate the data distribution before and after to normalization ( Fig 1A , B ). A total of 4764 DEGs were identified between the FE and FC groups, with 2618 up-regulated and 2146 down-regulated. On the other hand, 9736 DEGs were found between the IE and IC groups, comprising 4266 up-regulated and 5470 down-regulated DEGs (Supplementary S1 Table ). The DEGs between the FE and FC groups, as well as between the IE and IC groups, were depicted using a volcano plot ( Fig 1C , D ). According to the Venn diagram, 708 up-regulated and 414 down-regulated DEGs were found to be common in both FE and IE ( Fig 1E , F ). (A) Box plot of the expression data before and (B after normalization. (C) Volcano plot of DEG between FE and FC groups. (D) Volcano plot of DEG between IE and IC groups. (E) The Venn diagram of the overlapping up-regulated DEGs in FE and IE. (F) The Venn diagram of the overlapping down-regulated DEGs in FE and IE. The ShinyGO 0.81 tool provided 355 GO items for common up-regulated DEGs between the FE and IE groups, including 269 BP items (Supplementary S2 Table ), 42 MF items (Supplementary S3 Table ), and 44 CC items (Supplementary S4 Table ). Similarly, 405 GO items for common down-regulated DEGs between the FE and IE groups were obtained, comprising 294 BP (Supplementary S5 Table ), 44 MF (Supplementary S6 Table ), and 67 CC (Supplementary S7 Table ) items. Fig 2 illustrates the 10 most significant functional enrichment analyses of BP, CC, and MF catalogs for up- and down-regulated DEGs. The DAVID tool identified up-regulated DEGs in 18 KEGG enrichment pathways, with 10 being significant (P-value < 0.05) (Supplementary S8 Table ). Down-regulated DEGs were associated in 31 KEGG enrichment pathways, 21 of which were significant (Supplementary S9 Table ). (A) BP corresponding to common up-regulated DEGs, (B) MF corresponding to common up-regulated DEGs, (C) CC corresponding to common up-regulated DEGs, (D) BP corresponding to common down-regulated DEGs, (E) MF corresponding to common down-regulated DEGs, (F) CC corresponding to common down-regulated DEGs. The PPI network of shared up-regulated DEGs in the FE and IE groups, consisting of 614 nodes and 236 edges, was generated using STRING ( Fig 3 ). To enhance clarity, we removed isolated genes that lacked interactions with other genes. Supplementary S10 Table and Fig 4 indicate the top 10 genes identified through the five methods provided by Cytoscape’s CytoHubba plugin. The Venn diagram indicated a presence of two common genes, namely IL-6 and Kinase Insert Domain Receptor (KDR), among the genes identified by all five approaches ( Fig 5 ). Consequently, we selected these genes as potential hub gene candidates. The network comprises 614 nodes and 236 edges. Each node signifies a protein, while each edge denotes a protein-protein association. (B) The top 10 genes identified using BottleNeck method. (C) The top 10 genes identified using Closeness method. (D) The top 10 genes identified using Degree method. (E) The top 10 genes identified using Stress method. The colors represent high (red) to low (yellow) scores. The miRNAs/TFs-hub genes network was constructed using the NetworkAnalyst tool. The miRNAs/TFs-IL6 network comprises 7 TFs and 33 miRNAs (Supplementary S11 Table ), while the miRNAs/TFs-KDR network consists of 9 TFs and 22 miRNAs (Supplementary S12 Table ). In the network, hsa-miR-155-5p, hsa-miR-335-5p, hsa-miR-199a-3p, FOXC1, and STAT3 interacted with both hub genes ( Fig 6 ). The red circle, green diamond, and purple circle indicate hub genes, miRNAs, and TFs, respectively. The DGIdb 5.0 revealed that among the hub genes and their corresponding TFs, IL-6, KDR, and STAT3 had druggability potential ( Table 1 ). We chose the KDR (kinase insert domain receptor) gene for further investigation because it belongs to the kinase family, is located on the outside of the plasma membrane, and exists in both druggable genome lists from Hopkins and Groom (2002) and Russ and Lampel (2005). Synonyms for the KDR gene include vascular endothelial growth factor receptor 2 (VEGFR2), VEGFR, fetal liver kinase 1 (FLK1), and CD309, with VEGFR2 utilized in this article henceforth. The PDB format structure of the human VEGFR2 protein (ID: 1YWN, resolution: 1.71 Å) in complex with 4-amino-furo[2,3-d] pyrimidine (also known as leukaemia inhibitory factor (LIF)) was downloaded from RCSB PDB, and the preparation process was conducted on it (Supplementary S13 File). A collection of 16 FDA-approved drugs in SDF format was obtained from the DrugBank database (Supplementary S14 File ). This library included a list of drugs that possessed both an available structure and inhibitory characteristics ( Table 2 ). The molecular docking analysis, which used FDA-approved drugs as ligands targeting human VEGFR2, revealed valuable insights on molecular interactions. Molecular docking was carried out using a grid box including all of the binding residues, including Leu35, Gly36, Arg37, Gly38, Gly41, Val43, Ala61, Lys63, Glu80, Ile83, Leu84, Ile87, Val93, Val94, Val111, Glu112, Phe113, Cys114, Leu164, His171, Leu180, Ile189, Cys190, and Asp191. The results of molecular binding analysis revealed that 16 FDA-approved drugs had binding affinities ranging from −7.6 to −9.6 kcal/mol, with only one drug (ponatinib) having a greater binding affinity than the co-crystal ligand (LIF) and being chosen for further study. Table 2 presents the binding affinities of the FDA-approved drugs targeting human VEGFR2. The binding orientation and interaction of ponatinib and LIF within the binding pocket of human VEGFR2 were predicted using the UCSF Chimera software and ProLIF package, respectively. Fig 7A illustrates that ponatinib overlaps with the LIF and occupies the active site. LIF established different interactions with human VEGFR2, including 6 van der Waals (vdW) contacts (leu35, Val43, Glu112, Lys115, Gly117, Asn118) and 5 hydrophobic interactions (Leu35, Val43, Ala61, Gly117, Leu180) ( Fig 7B ). Likewise, ponatinib formed many interactions with human VEGFR2 including 10 vdW contacts (Leu35, Cys114, Lys115, Gly117, Asn118, Arg177, Leu180, Phe192, Arg196, Asp197) and 6 hydrophobic interactions (Val43, Asn118, Arg177, Phe192, Arg196, Asp197) ( Fig 7C ). (A) 3D interactions of LIF (red) and ponatinib (blue) with human VEGFR2 (tan). (B) Interactions map of LIF with human VEGFR2. (C) Interactions map of ponatinib with human VEGFR2. The RMSD is an fundamental parameter employed to evaluate structural deviations and the overall stability of proteins [ 62 ]. The RMSD plot of free VEGFR2 increased from 1.1 Å to 3.9 Å during the first 60 ns of the simulation, then fluctuated between 2.5 Å and 3.5 Å until the end. The RMSD plot of the VEGFR2-LIF complex indicates a minimum value of 0.9 Å and a maximum value of 2.9 Å, whereas the RMSD plot for the VEGFR2-ponatinib complex reveals a minimum value of 0 Å and a maximum value of 2.9 Å ( Fig 8A ). The fluctuations of residues during the simulation were evaluated using an RMSF analysis. The RMSF plot for the initial 75 residues of the free VEGFR2 exhibited greater fluctuations compared to the VEGFR2 complexed with LIF and ponatinib. Beyond this segment, with the exception of minor fluctuations, the RMSF plots s for VEGFR2, VEGFR2-LIF, and VEGFR2-ponatinib complexes overlapped throughout the simulation ( Fig 8B ). Rg denotes the folding rate and quantifies the protein’s compactness during simulation [ 63 ]; typically, a protein with a greater Rg indicates a less densely packed structure [ 64 ]. The maximal Rg values for VEGFR2, VEGFR2-LIF, and VEGFR2-ponatinib complexes were 21.05 Å, 21.2 Å, and 21.25 Å, respectively ( Fig 8C ). SASA is a measurement of a protein structure’s surface area that solvent molecules can access [ 65 ]. The SASA plots for VEGFR2, VEGFR2-LIF, and VEGFR2-ponatinib were nearly overlapping during the initial 23 ns. Within the range of 23–50 ns, the least SASA value was seen for VEGFR2, whereas the maximum SASA value was recorded for VEGFR2-ponatinib. During the final 20 ns of the simulation, the VEGFR2 plot fluctuated, whereas the VEGFR2-LIF and VEGFR2- ponatinib plots reached almost steady state ( Fig 8D ). (A) RMSD plot. (B) RMSF plot. (C) Rg plot. (D) SASA plot. The MM-GBSA and MM-PBSA methods were used to calculate the binding free energies of VEGFR2-LIF and VEGFR2-ponatinib complexes. The total MM-GBSA binding free energy was −82.51 ± 5.20 kcal/mol and −85.52 ± 5.06 for VEGFR2-LIF and VEGFR2-ponatinib complexes, respectively. The total MM-PBSA binding free energy was −82.38 ± 5.34 kcal/mol and −84.89 ± 5.65 kcal/mol for VEGFR2-LIF and VEGFR2-ponatinib complexes, respectively. The computed values indicate that the contributions of vdW energy are substantial in both methods. The contributions of each energy component are presented in Table 3 .

The computational findings of this study suggest that VEGFR2 is a promising target for endometriosis treatment. This study assessed FDA-approved drugs for their ability to target the VEGFR2 protein in diseased situations. The virtual screening approach commenced with molecular docking of 16 FDA-approved pharmaceuticals, followed by an evaluation of binding affinity. Ponatinib exhibited superior binding affinity compared to the ligand within the structure. Moreover, interaction analysis demonstrated that both occupy the identical binding pocket and interact with essential amino acids. MD simulations lasting 100 ns were used to analyze the docked complexes. MD simulations validated the structural stability of the ponatinib-VEGFR2 complex. Based on these findings, ponatinib has been discovered as a viable endometriosis treatment. However, more experimental and clinical trials studies are needed to corroborate these findings and investigate the therapeutic potential of ponatinib for endometriosis.

The GSE120103 dataset was obtained from the GEO ( https://www.ncbi.nlm.nih.gov/geo/(database (accessed on 3 February 2025) [ 22 ]. This gene expression profile was generated using the GPL6480 (Agilent-014850 Whole Human Genome Microarray 4x44K G4112F) platform and derived from a study conducted by Bhat et al [ 23 ]. The GSE120103 dataset contains 9 human endometrial tissue samples from each of the four groups: fertile women with endometriosis (FE), fertile women without endometriosis (FC), infertile women with endometriosis (IE), and infertile women without endometriosis (IC). Using “Limma” package [ 24 ] in R Studio [ 25 ], we identified differentially expressed genes (DEGs) between the FE and FC groups, as well as between the IE and IC groups, in the GSE120103 dataset. DEGs were selected based on |log2 fold change (FC)| > 1 and adjusted P-value < 0.05. The “ggplot2” package in RStudio was used to create the volcano plot [ 26 ]. The “Venn” package in R Studio was used to screen for common up- and down-regulated DEGs in the FE and IE groups [ 27 ]. GO functional analysis and KEGG pathway assessment of common up- and down-regulated DEGs between the FE and IE groups were conducted using ShinyGO 0.81 ( https://bioinformatics.sdstate.edu/go/ ) (accessed on 5 February 2025) [ 28 ] and DAVID ( https://davidbioinformatics.nih.gov/ ) (accessed on 5 February 2025) [ 29 ], respectively. The GO terms are divided into three categories: biological processe (BP), molecular function (MF), and cellular component (CC). In the analysis conducted in the ShinyGO 0.81 tool, False Discovery Rate (FDR) value of 0.05 was used as the significance cutoff, as well as, the significance level in the analysis conducted in the DAVID tool was set to P- value < 0.05. The ShinyGO 0.81 tool used an FDR value of 0.05 as the significance cutoff, while the DAVID tool used a P-value < 0.05 as the significance level. We used the STRING (version: 12.0) ( https://string-db.org/ ) (accessed on 6 February 2025) to develop a protein-protein interaction (PPI) network of shared up-regulated DEGs in the FE and IE groups. A minimum required interaction score of 0.700 (high confidence) was set for the generation of this interaction network. The database currently comprises 12535 organisms, 59.3 million proteins, and more than 20 billion documented interactions [ 30 ]. We then used the Cytoscape software (version 3.10.3) to further analyze and visualize the PPI network. Cytoscape is an open-source software platform designed for visualizing intricate networks and merging them with many types of attribute data [ 31 ]. We next used the CytoHubba plugin [ 32 ] to identify hub genes based on five topological analysis algorithms: Betweenness [ 33 ], BottleNeck [ 34 ], Closeness [ 35 ], Degree [ 36 ], and Stress [ 37 ]. The top 10 genes were chosen based on each method, and the shared genes among these five algorithms were considered as the final hub genes. The NetworkAnalyst tool ( https://www.networkanalyst.ca/ ) (accessed on 7 February 2025) was used to map interactions among important transcriptional regulatory transcription factors (TFs) and post-transcriptional regulatory microRNAs (miRNAs) with hub genes [ 38 – 40 ]. The JASPAR [ 41 ] and miRTarBase v9.0 [ 42 ] databases were used to generate the TFs-hub genes network and the miRNAs-hub genes network, respectively. The DGIdb 5.0 ( https://dgidb.org/ ) (accessed on 8 February 2025) was used to analyze the druggability of hub genes and their related TFs identified in our study. The DGIdb provides information regarding drug-gene interactions and druggability sourced from publications, databases, and various online resources [ 43 ]. This database categorizes genes as potentially druggable based on their presence in selected pathways, molecular functions, and gene families from the Gene Ontology, the Human Protein Atlas, IDG, “druggable genome” lists from Hopkins and Groom (2002) [ 44 ] and Russ and Lampel (2005) [ 45 ]. The target protein’s X-ray crystal structure was obtained from the RCSB Protein Data Bank (PDB) ( https://www.rcsb.org/ ) (accessed on 9 February 2025) in PDB format (Supplementary S1 File ) [ 46 ]. We removed all co-crystallized ligands, heteroatoms, and associated crystal water molecules from the parent structure using the UCSF Chimera software (version 1.10.2) [ 47 ] to ensure that pre-bound molecules did not influence the binding affinity of the predicted drugs. To fill in the missing residues in the target protein’s crystal structure, homology modeling was performed using SWISS-MODEL software ( https://swissmodel.expasy.org/ ) (accessed on 9 February 2025) and the parent sequence as a template [ 48 ]. The pKa values of the ionizable groups in the protein were ascertained using PROPKA version 3.1, facilitating the optimization of the hydrogen bond network at pH 7.4 [ 49 ]. Energy minimization was done using UCSF Chimera software (version 1.10.2) with default parameters such as steepest descent steps of 1000, steepest descent step size of 0.02 Å, conjugate step gradient of 10, conjugate steps of 0.02 Å, and update interval of 100. The DrugBank database ( https://go.drugbank.com/ ) (accessed on 10 February 2025) [ 50 ] was searched to identify the FDA-approved drugs that formed ligands for the target protein during the docking analysis. Among the identified drugs, those with inhibitory activity and available structures were retrieved from the PubChem database ( https://pubchem.ncbi.nlm.nih.gov/ ) (accessed on 10 February 2025) in structure-data file (SDF) format [ 51 ]. Molecular docking-based virtual screening was conducted to identify drugs with high binding affinity for the target protein. The AutoDock Vina Wizard [ 52 ] in PyRx 0.8 software was utilized for screening [ 53 ]. The appropriate file formats for the receptor and ligand for molecular docking were generated in PDBQT format using PyRx 0.8 software. Active sites were defined as regions within the target protein structure that were less than 5 Å distant from the co-cocrystal ligands. The molecular docking analysis was conducted using dimensions of X = 29.3683 Å, Y = 23.2777 Å, and Z = 17.8267 Å, respectively, and the center coordinates for X, Y, and Z were set as −1.5459, 35.5985, and 14.8893, respectively. The grid box contained all of the residues required for binding, including Leu35, Gly36, Arg37, Gly38, Gly41, Val43, Ala61, Lys63, Glu80, Ile83, Leu84, Ile87, Val93, Val94, Val111, Glu112, Phe113, Cys114, Leu164, His171, Leu180, Ile189, Cys190, and Asp191. To validate the docking method, we also performed redocking of the co-crystal ligand with the target protein. Following docking, hit drugs candidates were chosen based on binding affinity scores. The analysis of protein-ligand interactions is a crucial component of the drug discovery process, since it elucidates the molecular mechanisms that may account for the efficacy or ineffectiveness of specific drugs [ 54 ]. In this study, the ProLIF package was used to visualize docking results and analyze interactions between receptor protein and ligand residues in the docked complexes. ProLIF (Protein-Ligand Interaction Fingerprints) is a Python utility for generating interaction fingerprints for molecular complexes. These fingerprints are vector representations of molecular interactions in three dimensions. This occurs frequently between proteins and ligands [ 55 ]. The dynamic behavior of unbound receptor and docked complexes was investigated using MD simulation. The simulation protocol was executed with the AMBER 18 program [ 56 ]. The primary coordinates were derived from the unbound receptors, docked complexes, and topology files were prepared utilizing the tLEAP interface of AMBER 18. System solvation was performed with three-point transferable intermolecular potential (TIP3P) water, and the force fields used for computations were General Amber Force Field (GAFF) [ 57 ] and ff99SB. To neutralize the system’s overall charge, sodium counter-ions were introduced. To eliminate steric clashes, the docked protein complex underwent minimization through 1500 steepest descent and 1000 conjugate gradient steps. Langevin dynamics were utilized for system heating over 10 picoseconds (ps) [ 58 ], followed by 100 ps of equilibration in the canonical (NVT) ensemble. During the production run, hydrogen bonds were constrained using the SHAKE algorithm [ 59 ]. The temperature of the system was incrementally raised from 0 to 300 K over 200 ps at constant volume, after which the system was equilibrated at constant pressure. The production run lasted for a total of 100 nanoseconds (ns). The final MD trajectories for unbound receptor and docked complexes was examined. The analysis includes generation root mean square deviation (RMSD), root mean square fuctuation (RMSF), radius of gyration (Rg), and solvent accessible surface area (SASA) values for free receptor and docked complexes. The molecular mechanics Poisson–Boltzmann surface area (MM-PBSA) and molecular mechanics generalized Born surface area (MM-GBSA) methods were employed to compute the binding free energies of the docked complexes [ 60 ]. A total of 450 frames each after 0.2 ns were extracted from the complete MD trajectory and analyzed using the MM-PBSA computation via the MMPBSA.py module [ 61 ] of AMBER18. To determine the binding free energies (△G bind ) values, the following equation was used: ΔG complex represents the total free energy of the protein-ligand complex, while ΔG receptor and ΔG ligand indicate the total free energies of the separated protein and ligand in solvent, respectively. T denotes the temperature, while S signifies the entropy contribution to ligand binding, determined by established approximations. The gas phase energy (E gas ) is commonly derived from the force field’s MM. It encompasses contributions from internal energy, electrostatic interactions, and van der Waals interaction energies as follows: The ΔG solv term is computed via an implicit solvent model and divided into electrostatic and non-polar components. The ethical committee of Semnan University of Medical Sciences approved this study with the number: IR.SEMUMS.REC.1403.248.

(DOCX) (DOCX) (DOCX) (DOCX) (DOCX) (DOCX) (DOCX) (DOCX) (DOCX) (DOCX) (DOCX) (DOCX) (PDB). (PDB) (PDB). (PDB) (ZIP)