Computational investigation of the interaction mechanism between the estrogen related receptor α and its agonists

Abstract

Estrogen related receptor alpha (ERRα), the first identified orphan nuclear receptor, participates in the metabolism of body and adjusts some diseases such as breast cancer, osteoporosis and diabetes. Recent experiments have identified several ligands such as 7a, 7b, DK1 and DK3 toward ERRα as agonists, to have a therapeutic ability in reducing blood glucose, and impact the activity of ERRα receptor. Nevertheless, the detailed interaction modes between the agonists and the ERRα have not been fully understood, which will hinder the design and development of more potent agonists targeting ERRα for the treatment of diabetes. To address this issue, molecular docking and molecular dynamic simulation (MD) were performed to study how the four agonists regulate the behavior of ERRα. It was found that a proved residue Phe232 and some others such as Arg276, Phe286, His398 and Phe399 were the key residues in the ligand binding pocket. Hydrogen bonds between agonists and residues Arg276 and His398 and the pi–pi interactions between agonists and Phe232/Phe286/Phe399 both play important roles for agonist–ERRα binding. Interestingly, we observed that the agonists kept the binding of ERRα and its coactivator PGC-1α by stabilizing the site of helix12 (H12). All these findings are not only beneficial to understand the interaction mechanism between ERRα and its agonists, but also provide clues for designing novel and potent agonists of ERRα for the treatment of diabetes.

---

1 Introduction

Estrogen-related receptors (ERRs), members of nuclear receptor subfamily, contain three isoforms (ERRα, ERRβ and ERR γ).^1,2 Though significant homology is shared between ERRs and estrogen receptors (ERs), ERRs don't bind with endogenous estrogen or any known natural ligands of ERs. Even more, recent functional genomic experiments have clearly demonstrated that ERRs and ERs display strict binding site specificities.³ ERRα is the first orphan nuclear receptor identified by a low stringency hybridization screen in 1988, which is expressed predominately in metabolically active tissues such as intestine, brown adipose tissue and skeletal muscles, and functions as a key regulator of energy metabolism.⁴

Genetic and functional analysis have demonstrated that ERRα plays a critical role in regulating the energy of metabolism by controlling the transcriptional functions or protein expression of key components involved in tricarboxylic acid (TCA) metabolism,⁵ oxidative phosphorylation (OXPHOS) and glycolysis.⁶ Considering the simulation of energy and metabolism, ERRα can be studied as a potential therapeutic target in the treatment of such metabolic disorders as diabetes and obesity.^7–10 Previous studies have shown that ERRα stimulates downstream biological function by recruiting coactivator such as peroxisome proliferator-activated receptor γ coactivator 1α (PGC-1α).^11,12 And PGC-1α is an important coactivator in signal pathway of diabetes, so finding novel small molecules which can influence the activity of ERRα to solve metabolic diseases such as type 2 diabetes is intriguing.

Structural information of ERRα can provide us atomic details about the receptor and advance our understanding to the function of the receptor. Four crystal structures of ERRα have been resolved, which exhibit a canonical three layered “α-helical sandwich” fold composed of 12 α-helices (H1–H12) and two β-sheets.¹³ Two conformations of ERRα exist, the active one complexed with PGC-1α and the inactive one complexed with inverse agonists; both revealed a conserved structure in the binding site region. Nevertheless, significant divergences in both H12 and H3 were observed (Fig. 1). H12 is quite flexible, upon activation H12 binds with coactivator normally (the right figure of Fig. 1A), but in inverse agonist, a major conformational change takes place in the ligand binding domain (LBD) where H12 is repositioned to cap the ligand binding site, so the PGC-1α can't be combined for the position occupied by H12 (Fig. 1B and the left figure of Fig. 1A). The crystal structure¹⁴ of LBD of ERRα complexed with a coactivator peptide PGC-1α reveals a transcriptionally active conformation in the absence of a ligand. PGC-1α thereby plays a role in diverse processes including energy metabolism, mitochondrial biogenesis, and resistance to insulin in type 2 diabetes.^15,16 Binding of activating ligand (agonist) to the ligand binding pocket can also trigger conformational changes and then result in coactivator recruitment. But the ligand binding pocket is highly occupied by residues Phe232, Phe286 and so on, which make it difficult to accommodate synthetic ligands without disrupting its constitutive active conformation (Fig. 2).

Fig. 1 (A) Two crystal structures of ERRα complex in PDB bank, namely, ERRα complexed with the inverse agonist XTC-790 (PDB entry code: 2PJL) and the active structure of ERRα complexed with PGC-1α (PDB entry code: 3D24). (B) States of H3 and H12 showed the most obvious differences between the inactive conformation (2PJL) and active conformation (3D24) of ERRα.

Fig. 2 Chemical structures of ligands 7a, 7b, DK1 and DK3 and the ligand binding pocket of ERRα (PDB entry code: 3D24). The binding pocket is very small and locates in the deep site of ERRα protein.

Several classes of inverse agonists have been identified by high-throughput screening and chemical synthesis.^17–22 An X-ray crystallographic study demonstrated that the binding of inverse agonist XTC-790 with the ERRα LBD domain by interacting with Glu331, Val369 and one water in ERRα.²³ All these demonstrated the inverse activation mechanism, but small molecular agonists of ERRα are lacking likely because of the tiny binding pocket of the ERRα LBD.^14,24 Although the failure in identifying good ERRα agonists disappointed many scientists, several agonists, such as 7a, 7b, DK1, and DK3 ^25,26 have been found by the method of chemical synthesis finally and fortunately. All the four agonists (7a, 7b, DK1 and DK3) have a carbonyl group. 7a was an agonist that can effectively enhance ERRα-driven luciferase activity at 10.0 μM after a 24 h incubation, and 7b elevated the transcriptional function of ERRα by 2.7-fold at 10.0 μM because of the little difference between the two agonists. DK1 and DK3 were two agonists with similar scaffold to 7a and 7b. The EC₅₀ values of DK1 and DK3 are not clear, but DK3 should be emphasized because of its specific selectivity to ERRα. The other three can also form complexes with ERRγ and ERs. Further studies showed these agonists could enhance the transcription of ERRα downstream target genes and improve the glucose and fatty acid uptake in C2C12 muscle cells, but it is not yet clear on the exact mechanism about how these agonists activate the receptor.

In this study, in order to identify the activation mechanism, we concentrated on interactions between ERRα and four agonists (7a, 7b, DK1 and DK3). Molecular docking and molecular dynamics (MD) simulation were performed to investigate the binding interactions between the agonists and ERRα as well as the ligand-induced conformational changes of ERRα. With the initial conformation generated from docking, ten systems were constructed for subsequent MD simulations. In addition, we compared the differences between the specific selective agonist DK3 and the other three (7a, 7b and DK1). The results may avail us of understanding the activation mechanism of ERRα by agonists and provide clues for finding novel agonists of ERRα for the treatment of type 2 diabetes.

2 Materials and methods

2.1 Structure preparation

For ligands, the initial structures of the agonists 7a, 7b, DK1 and DK3 were geometrically optimized with the Macro Model module of the Schrödinger software package.²⁷ The optimized conformer was thus obtained based on energy minimization, and was saved in PDB format for input into the docking protocol.

For protein, the initial structure of ERRα was taken from Protein Data Bank (PDB). To date, there are four crystal structures available for ERRα. Two crystal structures complexed with inverse agonists were ignored, and one (PDB entry code: 3D24) of the two apo structures was chosen due to its high resolution (2.11 Å) finally.²¹ The crystal structure was then pretreated with the “Protein Preparation Wizard” workflow in Maestro version 9.3.²⁸ In order to get complete structure, two missing loops were built using Prime.²⁹

2.2 Molecular docking

We used the Autodock4.2 software package to predict the possible binding modes of agonists 7a, 7b, DK1 and DK3 to ERRα (3D24).³⁰ The chemical structures of ligands and the binding pocket of ERRα were shown in Fig. 2. The center of the grid was placed on the ligand binding domain (LBD) and the residues within 20 Å around the center were defined as the binding pocket. Before starting the docking protocol, the target receptor and ligands were prepared using standard docking protocol and were saved as “PDBQT” format and Lamarckian Genetic Algorithm was used to generate 30 outputs for each docking run.^31–33 In docking calculations, the target-ligand poses obtained were ranked using an energy based scoring function. After all outputs were clustered based on the root mean squared deviation (RMSD) values, the top pose conformation of docked ligand with the lowest energy in the biggest cluster was saved.

2.3 MD simulation

2.3.1 MD preparation. The starting structures for MD simulations were obtained from the docking poses and protonation states of the ionizable residues and histidines were determined based on the microenvironment and pK_a values calculated by PROPKA.³⁴ For the purpose of generating parameter files of ligands for MD simulations, geometric optimization and electrostatic potential calculation of the four agonists were conducted at the B3LYP/6-31G(d)³⁵ level using Gaussian 03 software package. Restrained electrostatic potential (RESP) fitting procedure was applied to obtain the partial charges for the atoms in the agonists.³⁶ Both the presence and absence of coactivator PGC-1α were considered in this work, which resulted in ten systems as listed in Table 1.

Table 1 Ten systems used for MD simulations to investigate the interaction mechanism of the four agonists (7a, 7b, DK1 and DK3) with ERRα

Simulation number	System name	Simulation time	Ligand	Coactivator
1	7a_ERRα_PGC-1α	200 ns	7a	PGC-1α
2	7a_ERRα	200 ns	7a	—
3	7b_ERRα_PGC-1α	200 ns	7b	PGC-1α
4	7b_ERRα	200 ns	7b	—
5	DK1_ERRα_PGC-1α	200 ns	DK1	PGC-1α
6	DK1_ERRα	200 ns	DK1	—
7	DK3_ERRα_PGC-1α	200 ns	DK3	PGC-1α
8	DK3_ERRα	200 ns	DK3	—
9	ERRα_PGC-1α	200 ns	—	PGC-1α
10	ERRα	200 ns	—	—

---

2.3.2 Molecular dynamics simulation. Molecular dynamics simulation can provide valuable insights into the structure and thermodynamics of complex biological systems. GROMACS4.6.5 package in a periodic boundary condition was used for our simulations. The similar protocol to that in our previous researches^37,38 was used. The AMBER99SB-ildn force field was applied for the proteins, and TIP3P model was employed to solvate the protein. Periodic boundary condition in all directions was applied in our simulations. Gaff force files were used to parameterize the agonists. Prior to MD simulations; each system was subjected to 50 000-step energy minimization to relieve internal repulsions with the steepest descent algorithm with a maximum force of 10.0 kJ mol⁻¹ nm⁻¹ as the threshold. After energy minimization, NVT equilibration was conducted for 50 ps at 310 K. After stabilizing the temperature, NPT equilibration was performed for another 1 ns at 1 bar using the Parrinello–Rahman pressure coupling. The long-range electrostatics was set to the particle mesh Ewald algorithm,³⁹ and constraints for all bonds with hydrogen atoms were applied using the LINCS algorithm. Following NVT and NPT equilibration, 200 ns production runs were carried out for ten systems using an NPT ensemble with a time step of 2 fs.

2.3.3 Principal component analysis. Principal component analysis (PCA) is another helpful tool for studying the differences in the protein dynamics, particularly those of the sub-structural elements of interest.^40–42 In PCA analysis, these structures were extracted as snapshots from the aligned MD trajectory, and then the covariance was calculated and diagonalized, providing an orthogonal set of eigenvectors representing linearly independent modes of conformational changes. These eigenvectors are the principal components (PCs). The eigenvalue associated with each eigenvectors is a measure of the variance in the original data set associated with the component.⁴³ In order to identify the conformational changes provoked by the LDB region and H12 position changes, the entire MD trajectories were included in these calculations. In this work, the PCAs were performed on the Ca coordinates of key residues in LDB region and PGC-1α, respectively, and the results were presented as the scatter diagram of the first two PCs (PC1 and PC2).

2.4 MM-PBSA binding free energy calculation

The binding free energy of five PGC-1α/ERRα complexes was calculated using g_mmpbsa, which was developed to apply the MM-PBSA method to the GROMACS package. G_mmpbsa implemented the MM-PBSA method in the C programming language, removing any runtime dependency on external software such as DelPhi and APBS.⁴⁴ In general terms, the binding free energy was conceptually summarized using the following equation:

ΔG_binding = G_complex − (G_protein + G_ligand) (1)

where G_complex, G_protein and G_ligand represent free energy of complex, protein and ligand, respectively.

As described in previous publications,⁴⁵ the free energy for each snapshot can be estimated from the following equation:

ΔG_binding = ΔE_MM + ΔG_solv − TΔS (2)

where ΔE_MM represents molecular mechanics potential energy including bonded (ΔE_bonded, always taken as zero) and nonbonded (ΔE_nonbonded) interactions (eqn (3)); ΔG_solv is the solvation free energy, including electrostatic solvation energy (ΔG_polar) and non-electrostatic solvation energy (ΔG_nonpolar) (eqn (5)); TΔS refers to entropic contribution to the free energy, but g_mmpbsa does not include the calculation of entropic terms.

ΔE_MM = ΔE_bonded + ΔE_nonbonded (3)

ΔE_nonbonded = ΔE_vdw + ΔE_elec (4)

In eqn (4), ΔE_vdw and ΔE_elec represent van der Waals and electrostatic energy, respectively.

ΔG_solv = ΔG_polar + ΔG_nonpolar (5)

In current work, we adopted one of the most widely used nonpolar models, SASA-only nonpolar model (eqn (6)), for non-electrostatic solvation energy calculation.

ΔG_nonpolar = γ(SASA) + β (6)

where solvent SASA represents accessible surface area, γ and β are related parameters. To further study the contribution of the agonists, 400 snapshots of each simulated system were used to calculate the binding free energy between coactivator PGC-1α and ERRα protein.

3 Results and discussion

3.1 Docking analysis of agonists in human ERRα

In this study, molecular docking approach was used to determine the binding modes of the four agonists 7a, 7b, DK1 and DK3 in ERRα. Since no prior information is available about agonist–ERRα complexes in Protein Data Bank. Then we selected the conformations of the four agonists in estrogen related receptors with the most favorable energy score. The binding poses of the four ligands were shown in Fig. 3. The four small molecules share very similar scaffold, the main difference between 7a and 7b is the existence or absence of methyl or ethyl in the side chain. The difference between DK1 and DK3 lies in the chlorine atom on the side chain. However, obvious distinctions were observed for the binding poses between 7a and 7b, and between DK1 and DK3 (Fig. 3). The docked conformations of the ligands can be divided into two clusters. One cluster includes 7a and DK3; carbonyl oxygens of them located near to Phe232, while in the other cluster including 7b and DK1, the carbonyl oxygens located near to Phe286. From these two kinds of docking poses, we suspected that agonists with this scaffold have two low energy poses in protein pockets. In addition, benzene rings of agonists 7b and DK1 were surrounded by Phe232, Phe399 and Phe414 in the pocket domain. Two pi–pi interactions were found between the benzene rings of agonists 7b and DK1 and the side chains of residues Phe232 and Phe399. On the contrary, no pi–pi interaction was found in 7a and DK3 docking systems due to the long distance between the benzene ring on the ligand and the residues mentioned above. Based on the docking results, it can be concluded that the Phe232 and Phe399 of ERRα may play important roles for the binding of agonists.

Fig. 3 The overall situation of the small molecule agonists docked into the protein crystal structure. Agonist 7a and DK3 share similar docking poses, and agonist 7b share similar docking pose with DK1.

3.2 Stability of the systems and activation mechanism of the receptor

3.2.1 Stability of each simulated system. In order to get the more accurate interaction mode and dynamic characteristics of agonists with ERRα, a 200 ns MD simulation was conducted for each of ten systems in this work. Firstly, the RMSD values of protein Cα atom and the ligand atoms relative to their initial structures were calculated to examine the structural stability of the ten systems during MD simulations. The RMSD values of agonists 7a, 7b, DK1 and DK3 in the first eight systems (Table 1) exhibited a fluctuation from 0 ns to 80 ns during the MD simulations and then became equilibrated. The atoms in the protein did not exhibit significant deviations from their initial structures. Based on the equilibrium trajectories of MD simulations, a variety of analyses were performed to find out the key residues in the binding pocket and to study the mechanism of action of ERRα.

The RMSF values of the ten systems were also calculated which were shown in Fig. 4. The H2–H3 loop (loop between helix 2 and helix 3, residues 214–219) and H10–H11 loop (loop between helix 10 and helix 11, residues 367–373) in protein, which were broken in the crystal structure and repaired by the method of prime, had extremely larger fluctuation than any other parts including other loops. When agonists were complexed, the conformation changed slightly, for only some residues around the ligand binding pocket and H12 performed conformational changes. Ignoring the two unstable repaired loops, other minor fluctuations were caused by the movement of several key residues such as Phe232, Phe286, His398 and Phe399. Fluctuations of the ten systems are similar no matter the agonists and coactivator PGC-1α was complexed or not. But the fluctuation value in the two DK3 systems seems a little lower than the others, which implicates that systems complexed with DK3 are more stable.

Fig. 4 RMSF of ten simulated systems. General trend of ten simulated systems' fluctuations is similar, but the value of DK3 systems seems a little lower than the others.

3.2.2 Interactions between agonists and ERRα. The interactions between agonists and ERRα of each simulated system were investigated in detail with the equilibrium trajectories of MD simulations. Although no significant fluctuation was observed after 80 ns in all ten systems, some slight fluctuations were still observed in systems DK3_ERRα and ERRα_PGC-1α, and these two systems were more stable in the last 20 ns. Hence, the following analysis was based on the last 20 ns equilibrium trajectory. We clustered the 400 snapshots which were saved from the last 20 ns equilibrium trajectory with an interval of 50 ps, and selected the center structure of the biggest cluster as the representational conformation of each simulated systems. All the 10 snapshots, including 8 complex systems and two apo structures, have been selected to analyze the interaction modes of agonists with ERRα.

From the results of interaction analysis (Fig. 5), we found that the 8 complex systems had different interaction modes. For the representational structure of system 7a_ERRα_PGC-1α, a hydrogen bond formed between the carbonyl oxygen of agonist 7a and the hydrogen atom of the imino group of the side chain of residue His398, and a pi–pi interaction was also found between the benzene ring of 7a and the residue Phe232. But for the structure of system 7a_ERRα, no hydrogen bond or pi–pi interaction were found. For agonist 7b, the snapshot of system 7b_ERRα_PGC-1α showed a hydrogen bond between the carbonyl oxygen of 7b and the hydrogen atom of the imino group of the side chain of residue His221, and a pi–pi interaction between the benzene ring of 7b and residue Phe232, too. But no hydrogen bond or pi–pi interaction was found in the snapshot of system 7b_ERRα. For agonist DK1, no hydrogen bond was found in systems DK1_ERRα_PGC-1α and DK1_ERRα, but two pi–pi interactions between DK1 and residues Phe286 and His398 were found in both systems with DK1. Finally, for agonist DK3, two hydrogen bonds formed between DK3 and ERRα. One hydrogen bond was between the carbonyl oxygen of DK3 and the hydrogen atom binding to the nitrogen atom of the guanidine group of Arg276, and the other formed between the nitrogen atom connecting with methyl in DK3 and the hydrogen atom binding to the other nitrogen atom of the guanidine group of Arg276. Besides the two hydrogen bonds, there was also a pi–pi interaction between the benzene ring of DK3 and Phe232 in system DK3_ERRα_PGC-1α. But no hydrogen bond was found in system DK3_ERRα, only the pi–pi interaction between DK3 and Phe286 existed. From above comparison, we can speculate that more hydrogen bonds and pi–pi interactions will form when PGC-1α is bound to ERRα simultaneously.

Fig. 5 Binding modes of the eight simulated systems complexed with agonists and key residues of the two simulated systems without agonists.

The hydrogen bond interaction is crucial for the binding and specificity of ligands to the receptor, so hydrogen bond analysis along the whole simulation time and hydrogen bonding frequency analysis of the 20 ns equilibrium time were performed. During the whole simulation time, the agonist 7a formed only one hydrogen bond with the residue His398; 7b formed only one hydrogen bond with the residue His221. Differently, DK1 did not form any hydrogen bond with either of these two histidine residues. While DK3 formed three hydrogen bonds with other residues (Lys208, Ala211 and Arg276) of ERRα, but the hydrogen bond between DK3 and Ala211, with only 1% hydrogen bond occupancy, was extremely weak. Considering the results of hydrogen bonding frequency analysis in Table 2, and neglecting the systems complexed with DK1, which had no any hydrogen bond, we can further know that when PGC-1α is bound, the chance of forming hydrogen bond is larger no matter the agonist is 7a, 7b or DK3.

Table 2 Hydrogen bonding frequency analysis of eight simulated MD systems

Simulation number	System name	Residues related	Hydrogen bond occupancy
1	7a_ERRα_PGC-1α	His398	40%
2	7a_ERRα	His398	2%
3	7b_ERRα_PGC-1α	His221	80%
4	7b_ERRα	—	—
5	DK1_ERRα_PGC-1α	—	—
6	DK1_ERRα	—	—
7	DK3_ERRα_PGC-1α	Arg276	100%
		Lys208	14%
		Ala211	1%
8	DK3_ERRα	—	—

---

Among the four agonists studied, DK3 is the only agonist which was known to selectively activate the ERRα. We further analyzed the behaviors of DK3 in the binding pocket of ERRα to discover the reason of the selectivity. Two obvious differences between DK3 and other three agonists appeared via above analyses. First, the lower RMSF value (Fig. 4) showed that the systems with DK3 are more stable. Second, we speculate that the hydrogen bond between DK3 and Arg276 (Table 2) contributes to the stability of DK3 simulated systems. What had already been proved was that the ligand would interact with Arg276 of ERRα, which had also been observed between ER/ERRγ and their selective ligands.⁴⁶ It was clearly demonstrated that Arg276 played a real important role for the binding of agonists to ERRs.

From all these analyses, we can conclude that residues Phe232, which had been verified as a key residue via the mutation experiment,⁴⁷ and Phe286 are the dominant residues which can form pi–pi interactions with ligands in most simulated systems, and His221, Arg276 and His398 are other three crucial residues to form hydrogen bonds with ligands in systems complexed with 7a, 7b and DK3, respectively. From this study, we knew those several residues around ligand binding pocket of ERRα are really important to impact the binding of ERRα and coactivator PGC-1α. All these interactions between key residues and ligands, especially the hydrogen bond of selective agonist DK3 can clearly show us the distinction of selective agonists of ERRα and help us find more novel ERRα agonists for the treatment of type 2 diabetes.

3.3 Impact of agonists on the binding of coactivator PGC-1α to ERRα

3.3.1 Key residues for ERRα/PGC-1α binding. To further study how these agonists influence the binding of ERRα and coactivator, we compared binding modes of the five simulated systems complexed with coactivator PGC-1α. We found the phenomenon that several residues (Glu419, Ala420) in H12 changed from helix to loop in ERRα_PGC-1α system, but the uncoiling had not been found in other four simulated systems with agonists (Fig. 6). In the system without agonist, we observed that M422, the C-terminal residue of H12, moved to the middle of H8–H9 loop and PGC-1α (Fig. 7), as a result the interaction between the loop and the coactivator was cut off by M422. But in other four systems with agonists, H12 is quite stable; the residue of M422 has not moved to the middle of H8–H9 loop and PGC-1α, either (Fig. 7). Two kinds of relative positions of H8–H9 loop and PGC-1α were found. One occurred in 7a and DK1 systems; the side chain of Arg205 in PGC-1α faced down and formed a gap between H8–H9 loop and PGC-1α (Fig. 7). While in 7b and DK3 simulated systems, it faced up, and the closer distance make the binding between H8–H9 loop and PGC-1α more easily. But all of the four systems with agonists have the capability to form hydrogen bonds via water molecule without the obstacle of M422. And a mentionable thing which can prove what we found reasonable is that one of the N-terminal flanking residues of PGC-1α, Arg205, contributes to the specific ERRα/PGC-1α interaction. Arg205 of PGC-1α forms a hydrogen bond with Gln262 (H4), and water-mediated hydrogen bonds with Ser337 (H8–H9 loop) and the main chain carbonyl oxygen of Ala420 (H12).⁴⁸ So, we can speculate that the agonists facilitate the binding of coactivator PGC-1α to ERRα via enhancing the stability of residues in H12, such as M422.

Fig. 6 (A) The change of secondary structure of H12 in ERRα_PGC-1α simulated system. (T: turn. E: extended conformation. B: isolated bridge. H: alpha helix. G: 3–10 helix. I: Pi-helix. C: coil (none of the above)). (B) In the simulated systems without agonists, several residues in H12 changed from helix to loop. Magenta: H12 of ERRα_PGC-1α system. Green: H12 of 7a_ERRα_PGC-1α system. Purple: H12 of 7b_ERRα_PGC-1α system. Blue: H12 of DK1_ERRα_PGC-1α system. Pink: H12 of DK3_ERRα_PGC-1α system.

Fig. 7 Binding mode of coactivator PGC-1α and protein ERRα. The terminal residue M422 of H12 located between H8–H9 loop and PGC-1α and hindered their binding when agonist is not complexed.

3.3.2 Stability of H12. Subsequently, the stability of H12 and PGC-1α were discussed because the state of H12 directly decides whether the PGC-1α can bind to ERRα or not. Representative structures selected with the same protocol as described in Section 3.2.2 were used to study the binding-mode of coactivator PGC-1α and the protein ERRα. It was known that the vertical H12, especially the terminal residue M422, can't promote the formation of the ERRα/PGC-1α complex. So, we speculate that the agonists improve the binding of ERRα and the coactivator via stabilizing the level state of H12.

To test the stability of H12 in our simulations, RMSD of H12 for all the 200 ns simulation trajectory in each system was calculated and shown in Fig. 8. RMSD value of H12 in ERRα_PGC-1α simulated system without agonist reached a value of ∼5 to 8 Å, which is much higher than the other four ERRα_PGC-1α systems with agonists (about 1 Å), meaning H12 in the system of ERRα_PGC-1α without agonist was not stable at all, while H12 in the ERRα_PGC-1α systems with agonists were much more stable. Through comparing the RMSD of H12 in the five systems, we believed that agonists contribute to the stability of H12 in the estrogen related receptor, and then the stable conformation of H12 is helpful to the recruitment of the coactivator. These results support our inference that the agonists improve the binding of ERRα and the activation factor via stabilizing the level state of H12.

Fig. 8 RMSD of H12 in the five simulated systems complexed with PGC-1α.

3.3.3 PCA analysis of key residues. To further study the blocking effect of Met422 of H12 in the system without agonist, the variations of the conformation of three important residues (Ser337, His341 and Met422), which are directly related to the binding between ERRα and coactivator, were studied by PCA analysis. All the 200 ns MD simulation trajectories were used for PCA analysis and the results were shown in Fig. 9. The 2D projection clearly suggested that ERRα_PGC-1α without agonist had greater conformational variations in the three residues, whereas in case of agonist-bound complexes, no such evident variations were observed. This further suggests that when agonists were not complexed, the three residues presented a variety of different conformations, while fewer kinds of conformations were found when agonists were complexed.

Fig. 9 PCA analysis of key residues in H12 and H8–H9 loop.

3.3.4 Binding free energy between ERRα and co-activator PGC-1α. Binding free energies of five PGC-1α-bound simulated systems were calculated to further study the activation mechanism of agonists to ERRα. The calculated energies including the total binding energy and the separated energy components were listed in Table 3. As shown, the calculated binding free energy of the system without agonist (−206.36 ± 2.43 kJ mol⁻¹) is higher than that of the other systems with agonists. The calculated binding free energy of the four systems complexed with agonists (7a, 7b, DK1 and DK3) is −19.47 kJ mol⁻¹, −91.14 kJ mol⁻¹, −52.99 kJ mol⁻¹ and −23.70 kJ mol⁻¹ lower than the system without agonist, respectively. The results demonstrated that the promotion of the binding of ERRα and coactivator PGC-1α is the activation mechanism of agonists to ERRα.

Table 3 The binding free energy of ERRα combined with PGC-1α, calculated by g_mmpbsa (kJ mol⁻¹)

System name	ΔGvdw	ΔGelec	ΔGsolv	ΔGSASA	ΔGbinding
ERRα_PGC-1α	−292.64 ± 0.79	−793.89 ± 2.51	914.46 ± 3.36	−34.33 ± 0.06	−206.36 ± 2.43
7a_ERRα_PGC-1α	−209.40 ± 3.96	−659.30 ± 19.06	668.66 ± 19.39	−26.35 ± 0.35	−225.83 ± 12.70
7b_ERRα_PGC-1α	−302.20 ± 1.97	−816.30 ± 10.45	857.57 ± 13.51	−36.79 ± 0.20	−297.50 ± 6.93
DK1_ERRα_PGC-1α	−275.38 ± 4.24	−636.50 ± 12.57	686.31 ± 15.43	−33.63 ± 0.41	−259.35 ± 10.99
DK3_ERRα_PGC-1α	−302.85 ± 2.21	−660.89 ± 6.47	769.42 ± 9.26	−36.25 ± 0.23	−230.06 ± 8.32

---

From various analyses, we concluded that the small molecular agonists allosterically affect the stability of H12, then contribute to the recruitment of coactivator, which is similar to the verified activation mechanism in other agonists and nuclear receptor interaction systems.^49–51

4 Conclusions

The newly-found four agonists of ERRα (7a, 7b, DK1 and DK3) were docked into the LBD of ERRα, which provided insights into the probable binding modes of agonists with ERRα. Subsequently, long-time MD simulations were performed for both apo-ERRα and ERRα/agonists complex systems, and binding mode of agonists to ERRα was depicted in detail. The docking results and interaction mode analysis of the MD trajectory displayed that Phe232, Arg276, Phe286, His398 and Phe399 around the binding pocket of ERRα are the key residues for the binding of agonists to ERRα. The results of MD simulations also reflected that the agonists may promote the binding of coactivator PGC-1α to ERRα by stabilizing the conformation and the site of H12. When agonists were absent, the site of H12 changed and hence hindered the binding of coactivator factor and protein, although the allosteric mechanism is not clear enough and more research is still needed, our model and findings may not only contribute to understand the binding mode of agonists to ERRα, but also supply useful information for searching more agonists of ERRα to help anti-diabetic drug design.

Acknowledgements

We gratefully acknowledged the financial supports from the National Natural Science Foundation of China (Grant 81273438) and the Fundamental Research Funds for the Central Universities (Grant WY1113007).

Notes and references

1. S. Sanyal, J.-Y. Kim, H.-J. Kim, J. Takeda, Y.-K. Lee, D. D. Moore and H.-S. Choi, J. Biol. Chem., 2002, 277, 1739–1748.
2. S. N. Schreiber, D. Knutti, K. Brogli, T. Uhlmann and A. Kralli, J. Biol. Chem., 2003, 278, 9013–9018.
3. G. Deblois, J. A. Hall, M.-C. Perry, J. Laganière, M. Ghahremani, M. Park, M. Hallett and V. Giguère, Cancer Res., 2009, 69, 6149–6157.
4. D. J. Mangelsdorf, C. Thummel, M. Beato, P. Herrlich, G. Schütz, K. Umesono, B. Blumberg, P. Kastner, M. Mark, P. Chambon and R. M. Evans, Cell, 1995, 83, 835–839.
5. A. Pedram, M. Razandi, B. Blumberg and E. R. Levin, FASEB J., 2016, 30, 230–240.
6. V. Giguère, Endocr. Rev., 2008, 29, 677–696.
7. C. Handschin and V. K. Mootha, Drug Discovery Today, 2005, 2, 151–156.
8. F. X. Soriano, M. Liesa, D. Bach, D. C. Chan, M. Palacín and A. Zorzano, Diabetes, 2006, 55, 1783–1791.
9. R. J. Patch, L. L. Searle, A. J. Kim, D. De, X. Zhu, H. B. Askari, J. C. O'Neill, M. C. Abad, D. Rentzeperis, J. Liu, M. Kemmerer, L. Lin, J. Kasturi, J. G. Geisler, J. M. Lenhard, M. R. Player and M. D. Gaul, J. Med. Chem., 2011, 54, 788–808.
10. V. Mootha, C. Handschin, D. Arlow, X. Xie, J. St Pierre, S. Sihag, W. Yang, D. Altshuler, P. Puigserver, N. Patterson, P. Willy, I. Schulman, R. Heyman, E. Lander and B. Spiegelman, Proc. Natl. Acad. Sci. U. S. A., 2004, 101, 6570–6575.
11. Y. Zhang, L. W. Castellani, C. J. Sinal, F. J. Gonzalez and P. A. Edwards, Genes Dev., 2004, 18, 157–169.
12. R. Bakshi, S. Mittal, Z. Liao and C. R. Scherzer, J. Parkinson's Dis., 2016, 5, 1–9.
13. R. S. Savkur and T. P. Burris, J. Pept. Res., 2004, 63, 207–212.
14. J. Kallen, J.-M. Schlaeppi, F. Bitsch, I. Filipuzzi, A. Schilb, V. Riou, A. Graham, A. Strauss, M. Geiser and B. Fournier, J. Biol. Chem., 2004, 279, 49330–49337.
15. J. Lin, C. Handschin and B. M. Spiegelman, Cell Metab., 2005, 1, 361–370.
16. B. N. Finck and D. P. Kelly, J. Clin. Invest., 2006, 116, 615–622.
17. O. Lanvin, S. Bianco, N. Kersual, D. Chalbos and J.-M. Vanacker, J. Biol. Chem., 2007, 282, 28328–28334.
18. H. Greschik, J.-M. Wurtz, S. Sanglier, W. Bourguet, A. van Dorsselaer, D. Moras and J.-P. Renaud, Mol. Cell, 2002, 9, 303–313.
19. B. B. Busch, W. C. Stevens, R. Martin, P. Ordentlich, S. Zhou, D. W. Sapp, R. A. Horlick and R. Mohan, J. Med. Chem., 2004, 47, 5593–5596.
20. F. Wu, J. Wang, Y. Wang, T.-T. Kwok, S.-K. Kong and C. Wong, Chem.-Biol. Interact., 2009, 181, 236–242.
21. J. Wang, F. Fang, Z. Huang, Y. Wang and C. Wong, FEBS Lett., 2009, 583, 643–647.
22. S. Xu, X. Zhuang, X. Pan, Z. Zhang, L. Duan, Y. Liu, L. Zhang, X. Ren and K. Ding, J. Med. Chem., 2013, 56, 4631–4640.
23. J. Kallen, R. Lattmann, R. Beerli, A. Blechschmidt, M. J. J. Blommers, M. Geiser, J. Ottl, J.-M. Schlaeppi, A. Strauss and B. Fournier, J. Biol. Chem., 2007, 282, 23231–23239.
24. S. M. Hyatt, E. L. Lockamy, R. A. Stein, D. P. McDonnell, A. B. Miller, L. A. Orband-Miller, T. M. Willson and W. J. Zuercher, J. Med. Chem., 2007, 50, 6722–6724.
25. L. Peng, X. Gao, L. Duan, X. Ren, D. Wu and K. Ding, J. Med. Chem., 2011, 54, 7729–7733.
26. K. Ding, C. Wong, Z. Kang and X. Zhou, CN Pat., ZL200810026782.0, 2011.
27. C. A. Romero, T. Grkovic, J. Han, L. Zhang, J. R. J. French, D. I. Kurtboke and R. J. Quinn, RSC Adv., 2015, 5, 104524–104534.
28. S. Saxena, P. B. Devi, V. Soni, P. Yogeeswari and D. Sriram, J. Mol. Graphics Modell., 2014, 47, 37–43.
29. Y. Tsutsumi, Funkcialaj Ekvacioj, 1987, 30, 115–125.
30. X. Hou, J. Du, J. Zhang, L. Du, H. Fang and M. Li, J. Chem. Inf. Model., 2013, 53, 188–200.
31. G. M. Morris, D. S. Goodsell, R. S. Halliday, R. Huey, W. E. Hart, R. K. Belew and A. J. Olson, J. Comput. Chem., 1998, 19, 1639–1662.
32. M. F. Sanner, J. Mol. Graphics Modell., 1999, 17, 57–61.
33. G. M. Morris, R. Huey, W. Lindstrom, M. F. Sanner, R. K. Belew, D. S. Goodsell and A. J. Olson, J. Comput. Chem., 2009, 30, 2785–2791.
34. C. R. Søndergaard, M. H. M. Olsson, M. Rostkowski and J. H. Jensen, J. Chem. Theory Comput., 2011, 7, 2284–2295.
35. J. E. Del Bene, W. B. Person and K. Szczepaniak, J. Phys. Chem., 1995, 99, 10705–10707.
36. M. K. Gilson, K. A. Sharp and B. H. Honig, J. Comput. Chem., 1988, 9, 327–335.
37. X. Sun, H. Ågren and Y. Tu, J. Phys. Chem. B, 2014, 118, 10863–10873.
38. X. Sun, J. Cheng, X. Wang, Y. Tang, H. Ågren and Y. Tu, Sci. Rep., 2015, 5, 8066.
39. T. Darden, D. York and L. Pedersen, J. Chem. Phys., 1993, 98, 10089–10092.
40. S. Wold, K. Esbensen and P. Geladi, Chemom. Intell. Lab. Syst., 1987, 2, 37–52.
41. M. A. Balsera, W. Wriggers, Y. Oono and K. Schulten, J. Phys. Chem., 1996, 100, 2567–2572.
42. M. L. Teodoro, G. N. Phillips and L. E. Kavraki, J. Comput. Biol., 2003, 10, 617–634.
43. F. Fratev, Phys. Chem. Chem. Phys., 2015, 17, 13403–13420.
44. R. Kumari, R. Kumar and A. Lynn, J. Chem. Inf. Model., 2014, 54, 1951–1962.
45. L. Zeng, M. Guan, H. Jin, Z. Liu and L. Zhang, Chem. Biol. Drug Des., 2015, 86, 1438–1450.
46. L. Wang, W. J. Zuercher, T. G. Consler, M. H. Lambert, A. B. Miller, L. A. Orband-Miller, D. D. McKee, T. M. Willson and R. T. Nolte, J. Biol. Chem., 2006, 281, 37773–37781.
47. W. Wei, A. G. Schwaid, X. Wang, X. Wang, S. Chen, Q. Chu, A. Saghatelian and Y. Wan, Cell Metab., 2016, 23, 479–491.
48. H. Greschik, M. Althage, R. Flaig, Y. Sato, V. Chavant, C. Peluso-Iltis, L. Choulier, P. Cronet, N. Rochel, R. Schüle, P.-E. Strömstedt and D. Moras, J. Biol. Chem., 2008, 283, 20220–20230.
49. C. Helsen and F. Claessens, Mol. Cell. Endocrinol., 2014, 382, 97–106.
50. L. Celik, J. D. D. Lund and B. Schiøtt, Biochemistry, 2007, 46, 1743–1758.
51. G. Costantino, A. Entrena-Guadix, A. Macchiarulo, A. Gioiello and R. Pellicciari, J. Med. Chem., 2005, 48, 3251–3259.

---