Metformin-Driven Activation of Polymorphic Follicle-Stimulating Hormone Receptors for Polycystic Ovary Syndrome Treatment: A Computational Study

Introduction

Polycystic ovarian syndrome (PCOS) is the most prevalent endocrinological disorder, affecting approximately 4–12% of women. It is recognized as a leading cause of infertility among women. Both genetic and environmental factors are hypothesized to contribute to PCOS development. Notably, PCOS is characterized by alterations in luteinizing hormone (LH), estrogen, testosterone, and, particularly, follicle-stimulating hormone (FSH) [1]. PCOS is often associated with an elevated LH: FSH ratio, accompanied by hyperandrogenism, menstrual irregularities, anovulatory infertility, and insulin resistance. FSH plays a crucial role in female reproductive health, contributing to follicle development, oocyte maturation, and regulation of steroidogenesis. FSH exerts its biological effects (gametogenesis and steroidogenesis) through the follicle-stimulating hormone receptor (FSHR), which is a member of the G-protein-coupled receptor family. The FSHR gene is located on chromosome 2p21, comprising 10 exons, 9 introns, and a promoter region [2]. Any dysfunction in FSHR could potentially disrupt ovarian function and folliculogenesis in PCOS. Such dysfunction includes increased aromatase expression upon FSH binding, as documented by Rice et al (2013) through the activation of cAMP-responsive transcriptional regulatory proteins in women with PCOS [3].

Metformin (N, N-dimethyl biguanide) has been widely used in individuals with PCOS and insulin-resistant conditions that modulate FSHR activity. The use of metformin offers several benefits, including improved insulin sensitivity, ovulation restoration, lower androgen levels, weight reduction, and a decreased risk of gestational diabetes mellitus [3,4]. However, the exact mechanism by which metformin targets FSHR and contributes to its therapeutic effects remains unidentified. FSHR polymorphisms were observed in women with PCOS, which may influence metformin-based treatment outcomes. Several studies emphasized the role of FSHR genetic variations in ovarian dysfunction, infertility, and pregnancy-related complications [5–9]. For instance, Zilaitiene et al reported a link between FSHR gene variants and conception [10], while other studies have shown that FSHR genotypes influence pregnancy outcomes for in vitro fertilization (IVF) [11]. These findings suggest that FSHR polymorphisms play a critical role in determining PCOS treatment outcomes. Given the impact of FSHR genetic variations on receptor function and treatment response, personalized therapeutic approaches based on FSHR genotyping may enhance treatment efficacy. While studies suggest that metformin indirectly influences FSHR expression, the precise molecular mechanisms through which these alterations affect ovulation and infertility remain unclear. The present study was performed to bridge this knowledge gap by using computational approaches to evaluate the agonistic (activation) effect of metformin on FSHR, as well as the structural and functional impact of FSHR variants in individuals with PCOS. Therefore, a series of computational approaches were used to assess the binding affinity of metformin with wild-type and variant FSHR structures. The workflow consisted of: (1) generating and evaluating the quality of wild-type and variant FSHR structures; (2) assessing protein stability influenced by these variants; and (3) performing molecular docking and dynamic studies with metformin.

Material and Methods

FSHR SEQUENCE AND PCOS VARIANTS:

The wild-type human FSHR protein sequence (accession ID: P23945) was obtained from the UniProt database, which offers extensive annotations, including protein sequences, single-nucleotide variants, structural characteristics, and functional insights. A literature review was performed to identify missense mutations in FSHR. Notably, polymorphisms at codons 307 and 680, particularly the A307T and N680S variants, have been frequently associated with PCOS [12–15]. Accordingly, the A307T and N680S variant sequences were generated by substituting threonine and asparagine residues at positions 307 and 680, respectively, from the wild-type structure.

HOMOLOGY MODELLING AND STRUCTURE VALIDATION:

The wild-type and variant FSHR sequences were analyzed using MODELLER (version 10.6) for subsequent structural and functional analysis. MODELLER efficiently generates 3D structural models referencing the RCSB Protein Data Bank (PDB) (https://www.rcsb.org/) to identify suitable and highly similar templates for structural modelling (http://www.salilab.org/modeller). Among the multiple generated models, the model with the lowest DOPE score was selected as the optimal structure, following the methodology of Haneef et al (2019) [16]. The DOPE (Discrete Optimized Protein Energy) function in MODELLER assesses structural quality by referencing a refined state. After optimization, the 3D model of FSHR was validated following the methodologies of Messaoudi et al (2013, Pourhajibagher (2022), and Mani et al (2023) [17–19]. The Structural Analysis and Verification Server (SAVES) version 6.0 (https://saves.mbi.ucla.edu/) was employed to validate and select the best model. The SAVES provides multiple tools, such as ERRAT, VERIFY3D, and PROCHECK, to assess the structure’s quality. The PROCHECK tool assesses the stereochemical quality of the protein structure, while VERIFY 3D evaluates the compatibility between the atomic model (3D) and its corresponding amino acid sequence (1D). Additionally, the ERRAT (http://services.mbi.ucla.edu/ERRAT/) distinguishes correctly and incorrectly modelled regions based on non-bonded atomic interactions.

ANALYSIS OF STRUCTURAL STABILITY:

Both the wild-type and variant (A307T and N680S) structures were thoroughly investigated to comprehend the structural and functional consequences of the variants. Mu-Pro [20] and I-Mutant 2.0 [21] tools were employed to assess the impact of each variant on FSHR structural stability. These tools are crucial for determining whether a mutation stabilizes or destabilizes the protein structure, offering valuable insights into agonist binding, and ultimately affecting therapeutic efficacy. Additionally, the HOPE server was used to investigate structural changes caused by the A307T and N680S mutations, focusing on variations in amino acid size, charge, and hydrophobicity. This analysis is essential for understanding how these mutations influence the formation of salt bridges and hydrogen bonds, which are critical for maintaining the structural integrity and functional stability of FSHR [22]. Furthermore, the RING web server (https://ring.biocomputingup.it/) was used to generate a protein residue interaction network for both wild-type and variant structures [23]. This network-based approach enables the identification of differences in key interactions, including hydrogen bonds, ionic interactions, π-cation and π-π stacking, Van der Waals forces, and disulfide bonds. These interactions highlight the impact of the substituted residue in intra-protein signal transduction and receptor function.

METFORMIN STRUCTURE RETRIEVAL AND DOCKING PREPARATION:

After examining the wild-type and FSHR variants, the potential agonistic binding of metformin to FSHR was explored to assess its impact on receptor activation. The chemical structure of metformin was retrieved in Structure Data File (SDF) format from the PubChem database (NCBI) using CID: 4091. Subsequently, the structure was optimized using the OPLS_2005 force field as mentioned by Roshni et al (2024) in the LigPrep module of Schrödinger, Maestro 11.2 [24]. The ligand size was restricted to approximately 500 atoms, and potential ionization states were determined using Epik at pH 7.0±2.0. The optimization process included desalination and tautomer generation while maintaining specified chiralities. Finally, the chemical structure of metformin was used to generate at least 32 distinct conformers. Similarly, the receptor grid for FSHR was generated using the Receptor Grid Generation module in Schrödinger. The grid box was carefully positioned around the active site of FSHR to ensure accurate docking simulations. The coordinates for the active site were derived from a previous report, which describes key allosteric sites for receptor activation and contributing to the agonist effect [25].

GLIDE MOLECULAR DOCKING PROTOCOL:

After grid generation, both wild-type and variant structures of FSHR and metformin were processed for molecular docking using the Glide docking tool in Schrödinger, Maestro 11.2 [24]. To achieve high-precision docking, the Glide XP (extra precision) refinement protocol was employed in flexible mode, allowing metformin to adopt the lowest energy binding conformation within or near the active site of FSHR. The optimal docking poses, characterized by the least binding energy, were selected for further analysis. The interactions between metformin and FSHR were then visualized using a two-dimensional interaction plot, providing insights into key binding interactions, hydrogen bonding, and spatial orientation within the receptor’s active site.

MOLECULAR DYNAMIC SIMULATION:

Molecular dynamics (MD) simulations were conducted using GROMACS 2021 [26] to assess the interactive stability of 3 complexes (metformin-wild-type, metformin-A307T, and metformin-N680S). The simulation system was built utilizing the CHARMM36 force field, which offers precise modelling of protein-ligand interactions and enhanced correlation with experimental data, as described by Huang et al (2013) [27]. The complex structures were placed in a dodecahedron box, ensuring appropriate distances from each edge. The system was then solvated using the TIP3P water model, and sodium or chloride ions were added to maintain neutrality to depict the biological conditions. To regulate system conditions, the Berendsen NVT ensemble was used to maintain a constant temperature of 300 K, while the NPT ensemble ensured a stable pressure of 1 atm with a fixed number of particles. Electrostatic interactions were handled using the Particle Mesh Ewald method with a 1.0 nm cutoff, and bond constraints were enforced using the LINCS algorithm. The complex structure underwent energy minimization using the steepest descent approach for 50,000 steps. Subsequently, MD simulations were conducted for 100 ns for each complex, following Nasebi et al (2024), which state that 100 ns duration provides reliable and quantitative results when using explicit solvent models like TIP3P [28]. For each complex, the trajectories, such as root mean square deviation (RMSD), root mean square fluctuation (RMSF), radius of gyration (Rg), solvent accessible surface area (SASA), and number of hydrogen bonds (HB), were calculated.

Results

RETRIEVAL OF WILD-TYPE AND VARIANTS OF FSHR SEQUENCE:

The wild-type FSHR sequence obtained from UniProt was modified by introducing amino acid variants, resulting in 2 distinct sequences for analysis: (1) A307T and (2) N680S. The MODELLER prioritizes templates for homology modelling based on sequence coverage, functional residue conservation, and structural resolution. The software aligns the target sequence with known protein structures from PDB databases, performing sequence alignment and calculating sequence identity to identify structurally similar templates. All 3 sequences (wild-type, A307T, and N680S) showed 52.56% similarity to the human thyrotropin receptor (PDB ID: 7UTZ), which was determined via electron microscopy at a resolution of 2.4 Å. Given that sequence identity above 30% is a strong indicator of model accuracy, as stated by Fiser et al (2014) [29], the 52.56% identity with the human thyrotropin receptor confirms its suitability as a template for homology modelling. After identifying the most appropriate template, the software generated comparative models for the wild-type, A307T, and N680S variants by sustaining spatial restraints derived from the template structure.

MODEL SELECTION AND QUALITY CONFIRMATION:

Among the various generated models, the final model for wild-type, A307T, and N680S was chosen based on its low DOPE scores. The selected models achieved a DOPE score of −0.46, reflecting favorable energy characteristics. Further, the models underwent comprehensive quality assessments via SAVES V6.0 as outlined by Ahmad et al (2023) [30]. The SAVES V6.0 incorporates several validation tools to assess structural quality. For instance, PROCHECK generated a Ramachandran plot for FSHR, showing that 88.8% of residues were in the “favoured” regions, 9.5% in the “allowed” regions, 1.1% in the “generously allowed” regions, and only 0.6% in the “prohibited” regions, which underscores the model’s better stereochemistry. In addition, the ERRAT analysis yielded a score of 95.14, far exceeding the acceptable threshold of 50, further confirming the structure’s high quality. The Verify3D evaluation also indicated that 86.37% of the residues were of high quality above the 80% cutoff required for a reliable model. Overall, all 3 models (wild-type, A307T, and N680S) met the structural accuracy criteria, ensuring reliability for further analysis.

STRUCTURE STABILITY PREDICTION AND NETWORK ANALYSIS:

Using stability prediction tools, the structural and physicochemical impacts of the A307T and N680S mutations were analyzed. The A307T substitution decreased protein stability, with IMutant2.0 reporting a DDG of −0.96 and Mu-Pro indicating a ΔΔG of −0.52. Similarly, the N680S mutation was predicted to reduce stability, as shown by a ΔΔG of −0.72 from IMutant2.0 and a ΔΔG of −0.56 from Mu-Pro. The HOPE analysis indicated that both variant residues, A307T and N680S, are located near a highly conserved region. In the wild-type, the residue is smaller and less hydrophobic, a characteristic that influences hydrogen bonding. Such surface residue disrupts interactions with other chemicals or protein components. Alterations at positions 307 and 680 can change the overall size of the protein compared to the wild-type. Specifically, in the A307T mutation, the substituted residue is larger, which might lead to steric clashes (“bumps”) and a loss of hydrophobic contacts due to its differing hydrophobicity. Conversely, in the N680S mutation, the variant residue is smaller than its wild-type counterpart, potentially resulting in fewer interactions. We analyzed the signal residue distribution in both wild-type and variant structures using the RING server. In the wild-type, the protein residue interaction network included 593 hydrogen bonds, 12 ionic interactions, 22 stacking interactions, 6 disulfide bonds, and 378 Van der Waals interactions. In contrast, the A307T and N680S variants showed 63 hydrogen bonds, 11 ionic interactions, 16 stacking interactions, 6 disulfide bonds, and 360 Van der Waals interactions. These differences in interaction counts indicate that the variants may exhibit altered protein stability compared to the wild-type, potentially impacting structural integrity, agonist binding, and receptor function.

MOLECULAR DOCKING:

After analyzing the structural features of the wild-type, A307T, and N680S variants, the agonistic effect of metformin on FSHR and its variants under insulin-resistant PCOS conditions was evaluated. The Glide docking tool was then used to predict the binding position and minimum energy conformation of metformin within the active regions of wild-type, A307T, and N680S FSHR. The active sites reported by Aathi MS et al (2022) [24] suggest that agonist binding can trigger FSHR activation. The Glide docking process involves generating ligand conformers, performing rough scoring, grid energy minimization, and final scoring. The analysis yielded the lowest binding energies of −4.205 kcal/mol for the wild-type, −4.321 kcal/mol for A307T, and −4.294 kcal/mol for N680S. Figure 1 illustrates the binding pose and interactions of metformin with wild and variant forms of FSHR. Notably, metformin interacts with GLU454 residue in all 3 protein structures through hydrogen bonding for its agonistic activity.

MD:

To confirm the interactive stability between wild-type FSHR and its variants with metformin, MD simulations were performed mimicking the dynamic biological environment. Various parameters (eg, RMSD, RMSF, Rg, SASA, and HB) were computed to assess the structural stability, flexibility, compactness, solvent accessibility, and molecular interactions for the metformin-bound wild-type, A307T, and N680S variants (Figures 2–6). As shown in the RMSD plot (Figure 2), the wild-type (black) displayed considerable fluctuations throughout the 0–100 ns simulation, suggesting reduced stability when metformin is bound. In contrast, the A307T (red) and N680S (grey) variants exhibit limited fluctuations during the same period, implying a more stable binding conformation. The RMSF was assessed over the simulation period to determine the average variations of the FSHR residue, with peak intensity indicating the level of each amino acid fluctuation (Figure 3). The amino acids between 293 and 330, as well as at terminal protein regions, showed marginal fluctuation, while the rest of the amino acids showed similar intensities between wild and FSHR variants. However, the RMSF trajectories of these protein-ligand complexes were favorable throughout the simulation time. The Rg analyses in Figure 4 and SASA in Figure 5 show no significant differences in compactness or environmental accessibility between the wild-type and variant-metformin complexes. Further, the time-dependent HB formation (Figure 6) was determined to evaluate the strong binding of metformin to the FSHR. Notably, the A307T-metformin complex (n=540) has a higher hydrogen bond compared to the wild-type (n=534) and N680S (n=533) FSHR-metformin complexes. Overall, the findings demonstrate the more stable binding behavior of metformin with the A307T and N680S variants compared to the wild-type throughout the 100-ns simulation.

Discussion

FSHR is an essential receptor involved in regulating the reproductive system in fertile women. Metformin, widely prescribed to improve insulin sensitivity, has also been found to indirectly regulate FSHR expression in women with PCOS [2,3]. Studies suggest that genetic variations in FSHR are linked to PCOS and can influence the therapeutic effectiveness of metformin. The present study examined the wild-type FSHR along with its reported polymorphic variants (A307T and N680S) to evaluate the agonistic binding of metformin for receptor activation and PCOS treatment in insulin-resistant individuals. The study utilized a systematic workflow integrating homology modelling, pathogenicity screening of variants, molecular docking, and MD simulations. The literature analysis revealed widely reported polymorphisms in the FSHR gene associated with PCOS. Among them, the A307T and N680S variants were the most frequently reported and demonstrated significant associations with the condition. Notably, these polymorphisms, located in exon 10 of FSHR, have been frequently observed in women with PCOS across different populations. Several studies have reported that alterations at residues 307 and 680 in FSHR may have biological significance by influencing ovarian response to hormones, altering the menstrual cycle, regulating ovarian hyperstimulation, and contributing to the development of PCOS. The present structural analysis of FSHR carrying the A307T and N680S polymorphisms revealed a significant impact on protein stability, as assessed by tools such as I-Mutant2.0 and Mu-Pro. The A307T substitution was predicted to decrease stability, with ΔΔG values of −0.96 (I-Mutant2.0) and −0.52 (Mu-Pro), while the N680S variant also exhibited reduced stability, with ΔΔG values of −0.72 and −0.56, respectively. Given the critical role of FSHR in follicular development and ovarian function, these destabilizing effects may impair receptor activity, potentially contributing to PCOS and influencing therapeutic responses. The HOPE analysis highlights significant structural changes and altered intra-protein interactions resulting from the A307T and N680S substitutions compared to the wild-type. Such alterations may impact the stability and function of FSHR.

Furthermore, molecular docking was performed using the Glide module of Schrödinger to assess the binding affinity and potential agonistic activity of metformin with the wild-type FSHR and its A307T and N680S variants. The docking results revealed binding scores of −4.321 and −4.294 kcal/mol for the A307T and N680S variants, respectively, compared to the wild-type (−4.205 kcal/mol) FSHR (Figure 2). Interestingly, metformin exhibited an increased binding affinity across both the variant FSHR (A307T and N680S) structures compared to wild-type, suggesting that these polymorphisms directly impact metformin’s interaction with the receptor. With a significant hydrogen-bonded interaction at glutamic acid 454, the docking poses demonstrated that metformin bound to (A) wild-type, (B) A307T, and (C) N680S. Further, MD simulations were then used to validate the docking results and examine the stability of wild-type, A307T, and N680S variants upon binding to metformin in a dynamic environment (Table 1). Various statistical parameters, including RMSD, RMSF, Rg, SASA, and HB (Figures 2–6), were analyzed based on their average and standard deviation (S.D) values.

RMSD generally measures the structural deviation of a protein over time. The RMSD analysis demonstrated that all studied complexes achieved equilibrium after 30 ns of simulation. Lower RMSD values indicate stable ligand binding. However, the wild-type complex exhibited greater structural fluctuations compared to the A307T and N680S variants. The average RMSD values for the wild-type, A307T, and N680S were 1.12±0.17 nm, 0.91±0.11 nm, and 0.844±0.11 nm, respectively. Notably, both variants showed reduced RMSD values, suggesting a more stable structure with metformin than the wild-type (Figure 2). The RMSF measures the flexibility of individual amino acid residues and quantifies how much each residue deviates from its average position over time (Figure 3). The average RMSF values for the wild-type, A307T, and N680S were 0.41±0.36 nm, 0.38±0.29 nm, and 0.39±0.33 nm, respectively. Both variants had slightly reduced residual fluctuations and exhibited a more rigid conformation, indicating localized stability with metformin compared to wild-type. The Rg reflects the overall compactness and structural stability of a protein. The average Rg values for the wild-type, A307T, and N680S were 4.06±0.08 nm, 4.07±0.06 nm, and 4.13±0.07 nm respectively. As shown in Figure 4, both the A307T and N680S variants, upon metformin binding, maintained consistently favorable Rg values, indicating a stable structure for binding. Likewise, SASA measures the extent of a protein’s surface exposed to the surrounding solvent. The average SASA values for the wild-type, A307T, and N680S were 362.39±2.56 nm2, 362.12±2.36 nm2, and 364.27±2.21 nm2, respectively, demonstrating a more compact structure for metformin binding. These results indicate that the A307T and N680S mutations do not substantially alter the solvent accessibility or structural compactness of the FSHR protein, indicating that these variants likely maintain their functional integrity (Figure 5).

Additionally, during the 100-ns simulation, HB analysis was performed to investigate the interactions within the FSHR-metformin complex. Figure 6 illustrates the number of hydrogen bonds formed in both the wild-type and variant systems. The average HB values for the wild-type, A307T, and N680S were 498.66±10.58, 501.10±11.40, and 486.82±11.15, respectively. Notably, A307T forms slightly more hydrogen bonds with metformin throughout 100 ns, potentially contributing to its stability. Overall, the present simulation results indicate that metformin exhibits stronger binding to the A307T and N680S variants compared to the wild-type, maintaining agonist activity, as supported by docking scores, minimal RMSD values, and low residual fluctuations.

The present computational analyses underscore the critical role of genetic variations in FSHR (A307T and N680S) in the context of metformin treatment for insulin-resistant PCOS conditions. Notably, metformin maintains a consistent ability to regulate FSHR in both A307T and N680S variant forms compared to wild-type individuals. However, this study has certain limitations, including the inability of docking scores to fully capture the dynamic nature of protein-ligand interactions, and the relatively short simulation timescale, and force field constraints that may influence long-term complex behavior, interaction energies, and structural stability. Further experimental validation through site-directed mutagenesis and biochemical binding assays is needed to confirm these findings. Future research should focus on identifying patients with specific FSHR variants to optimize therapeutic strategies, enabling more effective and targeted interventions for PCOS. Additionally, exploring metformin’s efficacy in relation to other FSHR polymorphisms and clinical factors, such as individual metabolic responses contributing to PCOS pathology, and identifying broader pathological targets, may aid in refining treatment strategies.

Conclusions

This study investigated the structural and functional behavior of the FSHR protein and its binding interactions with metformin in the presence of A307T and N680S variants associated with PCOS. Amino acid substitutions influenced the stability of the FSHR protein, leading to increased rigidity upon metformin binding. This analysis suggests that individuals with the A307T and N680S variants may experience therapeutic benefits comparable to those with the wild-type FSHR when treated with metformin.