How does stress affect human being-a molecular dynamic simulation study on cortisol and its glucocorticoid receptor.

Stress can be either positive or negative to human beings. Under stressful conditions, the mental and physical conditions of human can be affected. There exists certain relation between stress and illness. The cortisol and other glucocorticoids bind to the same receptor, which is called glucocorticoid receptor. Some evidences indicated that cortisol molecule binding to its glucocorticoid receptor was necessary for the stress response. Up to now, the structure-function relationships between cortisol molecule and its glucocorticoid receptor have not been deliberated from the atomic-level. In order to get a detailed understanding of the structure-function relationships between the cortisol molecule and glucocorticoids receptor, we have carried out molecular dynamic (MD) simulations on glucocorticoid receptor (Apo system) and cortisol with its glucocorticoid receptor complex (HCY system). On the basis of molecular dynamic simulations, a couple of key residues were identified, which were crucial for the binding of cortisol molecule. The results of binding free energy calculations are in good agreement with the experiment data. Our research gives clear insights from atomic-level into the structural-functional aspects of cortisol molecule and its glucocorticoid receptor, and also provides valuable information for the design of drug which can treat stress related illnesses.

Introduction

As the pace of life grows faster and faster and when we feel everything has become too much, we generally use the word “stress”. We are overloaded and wondering whether we really can cope with the pressures. That is how we feel the stress in daily life. From the physiological or biological point of view, the stress can be described in such way as when there is a stressor, like an environmental condition or stimulus, the organisms themselves can respond to it. Using this method, the body can react to stress. Regarding stressful event, the body's response is to activate sympathetic nervous system. As it is impossible for the body to keep this state for a long time, the body will return to a normal physiological condition with the help of parasympathetic system. For human beings, stress can be either a positive or a negative condition. Under such conditions, the mental and physical properties of human can be affected. There exists certain connection between stress and illness. Several studies (Schneiderman et al., 2005) pointed out that both acute and chronic stresses can cause illness. Stress can also make human more susceptible to physical illness such as common cold (Cohen et al., 1997). It is of great importance to develop certain drugs which can relieve stress.

The cortisol and other glucocorticoids bind to the same receptor, which is called glucocorticoid receptor. Like other steroid receptors (Kumar and Thompson, 1999), the glucocorticoid receptor represents as a modular structure (Kumar and Thompson, 2005) and consists of the following domains (marked A–F) including A/B-N-terminal regulatory domain; C-DNA-binding domain; D-hinge region; E-ligand-binding domain and F-C-terminal domain. The ligand binding domain is the region where cortisol binds. It has been reported that cortisol molecule binding to its glucocorticoid receptor is necessary for the stress response (Kolodkin et al., 2013). Unfortunately, the detailed information about the cortisol molecule binding to its glucocorticoid receptor has not been studied yet. To this end, we have investigated the atomic-level structural characterization of glucocorticoid receptor (Apo system) and cortisol with its glucocorticoid receptor complex (HCY system). Molecular dynamics (MD) simulations can be used as an effective way to study the conformational changes on atomic level (Wang et al., 2013). In this study, MD simulations for Apo and HCY systems were carried out. The aims of this work are to figure out the details about cortisol molecule binding to glucocorticoid receptor and to identify the key residues which are responsible for the cortisol binding. Our work provides detailed atomistic insights into the structure–function relationships between cortisol molecule and its glucocorticoid receptor, and also provides valuable information for the design of drug which can treat stress related illnesses.

Computational methods

Initial structures

The crystal complex structure of cortisol and its glucocorticoid receptor was retrieved from the RCSB Brookhaven Protein Data Bank (PDB entry: 4P6X (He et al., 2014), which served as the starting structure for the following molecular dynamic (MD) simulations. Only the chain A of the crystal complex remained. The protonation states of ionizable residues were determined at pH = 7.0 using H++ server (Gordon et al., 2005), which can predict the pKa value of protein residues at a given pH. The prepared complex structure (HCY system) was used as the starting structure of the subsequent MD simulations. The cortisol molecule was removed from this prepared structure to create the Apo form of glucocorticoid receptor.

Molecular dynamic (MD) simulations

MD simulations for both HCY and Apo systems were implemented in AmberTools15 by using sander.MPI module. The 99SB force field (Hornak et al., 2006) was chosen to be the force field for the protein. The force field parameters for cortisol molecule were supplied by general AMBER force eld (Wang et al., 2004). Two sodium ions (Na) were added to each of the two systems using coulomb potential grid in order to keep the whole system neutral. TIP3P water model (Jorgensen et al., 1983) was selected to solvate both systems using a truncated octahedron box. The size of the water box was set to 10 A distance around the solute molecule. The two systems were first carried out for 2000 steps minimization by employing the decent method and then for 3000 steps conjugate minimization of the entire systems. Then the two systems were heated from 0 to 300 K. The time scale for this process was 1000 ps. The ensemble for the heating process was the canonical ensemble (NVT ensemble).

During this process, a force constant of 10.0 kcal mol− and a harmonic restraint were applied on the protein and small cortisol molecules. The Langvie thermostat was employed to maintain the temperature, and then the two systems were equilibrated for 2000 ps. During this process, the NPT ensemble was adopted and the constant pressure was set to 1.0 bar. The total relaxation time for the barostat bath was set to 2.0 ps. In the end, the Apo and HCY systems were both simulated for 100 ns. The periodic boundary conditions were employed in this research. The long range electrostatics was handled by the particle-mesh Ewald (PME) method (Darden et al., 1993). The cut-off value for short range interactions was set to 10.0 A. Shake algorithm was employed to hold fixed bonds involving hydrogen. The time step for all the simulations was all set to 2.0 f s.

MD trajectories analysis

The MD simulations were carried out for both Apo and HCY systems for 100 ns. To obtain thoughtful insights into the motion behavior during the 100 ns simulation time, the trajectory obtained by MD simulation was analyzed. These trajectories were processed with AmberTools1.5 module. Root-mean-square deviation (RMSD) was employed to quantify the conformational changes of the same protein. This value was an important criterion in judging the structures of protein. In this study, C-RMSD was calculated for all systems with the first frame as reference structure. Root-mean-square fluctuations (RMSF) were used to evaluate the fluctuation of each residue during the simulation time. Hydrogen bonds are of great importance to biological molecules. We employed the following criteria for the hydrogen bonds analysis: The cut-off value of distance between the two heavy atoms was set to 3.0 A; the angle between acceptor and donor atom for hydrogen bonds employed a 120 cut-off value. The cluster analysis was also employed for the trajectories analysis. In order to visualize the trajectory and to present the structures, VMD (Humphrey et al., 1996), Chimera (Pettersen et al., 2004) and PyMOL (DeLano, 2002) softwares were used.

MM-GB/SA calculations

The MM-GB/SA methods (Wang et al., 2013) were applied to estimate the binding free energies between the ligand and its receptor. The binding free energy (Gbind) in MM-GB/SA between a ligand (L) and a receptor (R) to form a complex RL can be calculated as below

In the Eq. (2), the EMM represents the molecular mechanics component, which was determined in gas phase. Gsol is the stabilization energy caused by solvation. TS represents the vibrational entropy term. The EMM term is a sum of three terms as shown in the Eq. (3). Eint, Eele, and Evdw are internal, Coulomb and van der Waals interaction terms, respectively. The Eq. (4) gives the solvation energy, Gsol. As shown in this equation, the solvation energy can be divided into two terms including the electrostatic solvation free energy (GGB) and the nonpolar solvation free energy (GSA). The GGB term is obtained by the Generalized Born (GB) method. The GSA term can be proportional to the molecular solvent accessible surface area (SASA) method (Hou et al., 2008). The binding free energies were obtained by averaging over the values calculated for 3000 snapshots from the last 30 ns of the trajectories at 5 ps intervals for the complex structure.

Energy decomposition was carried out to anchor the pivotal residues which are crucial and responsible for the binding process of the cortisol molecule from the energetic point of view. We only performed the pre-residue decomposition calculation with the aid of MMGBSA module in AMBER11, wherein the decomposition energy was calculated by the following equation:

Results and discussion

The overall structural features or Apo and HCY systems

As a crucial and fundamental criterion for MD simulation, Root-mean-square deviation (RMSD) can provide an overall assessment of the simulated structures on the conformational changes. By means of contrasting the variations of positions of all atoms of the system with a selected or reference structure, we can obtain the RMSD value of the system. In this study, the reference structure during the RMSD calculations was set as the initial frame of each system. The RMSD plots for Apo and HCY systems are manifested in Fig. 1, from which we can see the structures of Apo and HCY systems reach equilibrium during the 100 ns simulation time. However, there are some details to which we should pay attention. The structure of Apo system changes significantly referenced to initial structure during the time of the first 20 ns while the structure of HCY system keeps fairly stabilized during the same period of time. From 20 to 70 ns, the RMSD value of Apo system increases rapidly, with the highest RMSD value almost reaching 3A. The structure of HCY system during this time still remains stable. Both systems reach equilibrium in the last 30 ns.

Fig. 1

Calculated root-mean-square deviations (RMSD) of the backbone atoms referenced to the corresponding starting structure. Black lines represent Apo system and red lines represent HCY system.

Root-mean-square fluctuation (RMSF) is also an important criterion for simulated structures. This criterion can afford specifics on the fluctuation of every single residue of the simulation system. Residues with high RMSF values are of large flexibility, which are marked in Fig. 2. As illustrated in Fig. 2, the overall structural flexibility is quite similar. All residues with high RMSF values locate in the loop region of the structure, staying far away from the cortisol binding region. Therefore, these flexible residues may not affect the function of the Apo and HCY systems.

Fig. 2

Root-mean-square fluctuations (RMSF) of the Apo and HCY system during 100 ns MD simulations. Flexible residues are labeled as Ser92, Leu111 Ser184 and Asn253.

Hydrogen bonds analysis

As we all know, hydrogen bonds make an important contribution to the structural and functional fields of biological macromolecules. To this end, hydrogen bonds between cortisol molecule and surrounding residues in HCY system are supervised during the 100 ns. The results of the analysis are shown in Fig. 3 and Table 1.

Fig. 3

Root-mean-square fluctuations (RMSF) of the Apo and HCY systems during 100 ns MD simulations. Flexible residues are labeled as Ser92, Leu111 Ser184 and Asn253.

Table 1

The properties of hydrogen bonds between cortisol molecule and its surrounding residues.

Hydrogen bonds	Occupancies (%)	Distances (A)	Angle ()
N39-HCY265	97.27	2.87		15.18
T214-HCY265	96.03	2.93		22.51
HCY265-N39	42.15	3.19		28.01
Q117-HCY265	37.03	2.79		15.69
HCY265-Q45	6.06	3.08		23.64
HCY265-T214	6.05	3.01		44.46
HCY265-Q45	5.17	3.05		31.16

As can be seen in Fig. 3 and Table 1, cortisol molecule in HCY system totally forms seven hydrogen bonds with its adjacent residues. The residues which are responsible for the binding of cortisol molecule are Asn39, Gln45, Gln117 and Thr214. From Fig. 3 we can see that the hydrogen bond Asn39@OD1-HCY265@HO2 and Thr214@OG1-HCY265@H5 exist during almost whole simulation time, with occupancies of 97.27% and 96.03%, respectively. These two hydrogen bonds are mainly in charge of the binding of cortisol molecule. The oxygen atom of amide group of Asn39 and Thr214 form hydrogen bonds with cortisol molecule. With the help of these two hydrogen bonds, the cortisol molecule stays stable in the binding pocket.

For the sake of getting a better comprehension of the electrostatic properties of cortisol molecule, electrostatic potential surface analysis was carried out. The multiwfn software (Lu and Chen, 2012, Ha, 2016) was used for electrostatic potential analysis and the VMD (Humphrey et al., 1996) was employed for depiction. The results of electrostatic potential analysis are shown in Fig. 4, from which we can see that the cortisol molecule itself is a negative and positive clearly molecule. That is, cortisol possesses a positive maximal and negative minimal point. It has been reported that cortisol molecule binding to its glucocorticoid receptor is necessary for the stress response (Kolodkin et al., 2013). Therefore getting a better understanding of the properties of cortisol molecule is important for further design of drug which can treat stress. Based on the electrostatic potential analysis, we can design the inhibitor of glucocorticoid receptor from the electrostatic point.

Fig. 4

The electrostatic potential surface of the cortisol molecule. The positive maximal point and negative minimal point are labeled.

Cluster analysis

Clustering analysis is a powerful technique for analyzing the trajectory produced by MD simulation. This technique can provide in-depth insight into the structural specifics of Apo and HCY systems. By means of the average linkage algorithm (Shao et al., 2007, Liu, 2014), the two trajectories were separated into five clusters. Among these five clusters, five snapshots were selected as the representative structures of each cluster. These five structures are named as follows: C0, C1, C2, C3 and C4, wherein C0 is the initial structure. The clustering analysis results are illustrated in Fig. 5 and Table 2.

Fig. 5

The results of cluster analysis. The representative structure of Apo system is colored in turquoise and that of HCY system in medium purple. A close view for the cortisol binding pocket is on the right. Key residues are labeled.

Table 2

Populations of representative structures and their RMSD values to crystal structure.

System	Structure	Population (%)	RMSD value to crystal structure (A)
Apo	C4	48.3	1.33
HCY	C0	59.2	0.96

From Table 2 we can see that the most popular cluster for Apo system is C4, with a population of 48.3%, while the most popular cluster for HCY system is C0, with a population of 59.2%. The C4 and C0 were superimposed on the crystal structure (PDB ID: 4P6X) for the purpose of verifying the structural differentiation between each representative structure and crystal structure. As shown in Table 2, Apo system possesses a higher RMSD value (1.33 A) as compared with HCY system, which indicates that the overall structure of glucocorticoid receptor changed when cortisol molecule was absent. A close view was made in order to find out the detailed structural differentiation between the two systems. It has been discussed in hydrogen bond analysis section that Gln45 formed hydrogen bond with cortisol molecule, while inn Apo system, the spatial position of this residue changed quite a lot. This information may be of great importance in the design of glucocorticoid receptor inhibitor to treat stress.

Relative motions of Apo and HCY systems

Principal component analysis (PCA) was designed to separate a conformational space of proteins into essential subspaces. The essential subspace containing different degrees of freedom can give a description of protein motions (Amadei et al., 1993). In this study, PCA was carried out with the aid of the ProDy software package (Bakan et al., 2011). VMD (Humphrey et al., 1996) and its plugin NWWiz (Bakan et al., 2011) were employed for the visualization of 3D structural snapshots. The results of the PCA are shown in Fig. 6, from which it is easy to find the remarkable motions in the two systems points on different parts of the structures. For HCY system, the most distinctive motion was the upper region of the structure. While for Apo system, no significant motions were observed. The relative motion information for HCY system is also helpful to the design of effective inhibitor of glucocorticoid receptor.

Fig. 6

A porcupine plot in stereo showing the Apo and HCY structures with cones signifying the first eigenvectors movements.

Binding free energy Aspect for HCY system

To describe the binding ability of cortisol molecule, MM-GB/SA calculations was performed. The 3000 frames of last 30 ns trajectory were used for binding free energy calculations. The results of binding free energy calculations are illustrated in Table 3, from which we can see the calculated binding free energy (−12.48 kcal mol−) is in good agreement with the experiment data (−7.096 kcal mol−). It was worthy to note that the electrostatic interaction (Eele) played an important role in binding cortisol molecule of HCY system. As mentioned above, cortisol molecule possessed a positive maximal and negative minimal point, which enabled the Eele term to be the driving force in binding cortisol to glucocorticoid receptor.

Table 3

Binding free energy (kcal mol−1) and its components of HCY system.

System	Eele	Evdw	GGB	GSASA	GMM GB = SAa	−T S	GT OTb	Gexpc
HCY	−23.72 5.39	−53.83 2.59	46.89 4.67	−3.34 0.072	−34.00 4.20	−21.51 3.77	−12.48 5.64	−7.096 (Kolodkin et al., 2013)

GMM GB = SA = Eele + Evdw + GGB + GSASA.

GT OT = GMM GB = SA.

Gexp = RTlnKD.

The residues which are responsible for the cortisol molecule binding were explored by the per-residue binding free energy decomposition analysis. The results are manifested in Fig. 7 and Table 4. As illustrated in Fig. 7, the residues with a lower than −1.5 kcal mol− contribution of binding free energy are labeled. The results of hydrogen bond analysis indicate that Asn39 and Thr214 form high occupancies of hydrogen bonds with the cortisol molecule. The energy decomposition results show that these two residues possess the two highest decomposition energies among all the residues, which are −2.71 kcal mol− for Asn37 and −2.63 kcal mol− for Thr214. As listed in Table 4, the interactions between the cortisol molecule and these two residues are mainly composed of electrostatic interaction. This result is consistent with the previous electrostatic potential surface and binding free energy analysis. Therefore, the electrostatic interaction is the driving force in the cortisol binding process. The crucial residues identified by energy decomposition analysis are of great importance for further design of inhibitor of glucocorticoid receptor.

Fig. 7

The decomposition binding free energy of HCY system. The residues that contribute over −1.5 kcal mol− are labeled.

Table 4

The decomposition binding free energy (kcal mol−1) and its components of HCY system.

Residue	van der Waals	Electrostatic	Polar solvation	Non-polar solvation	Total
Met35	−1.67	−0.58	0.45	−0.093	−1.90
Leu38	−2.51	0.26	0.15	−0.23	−2.33
Asn39	−1.50	−4.47	3.39	−0.14	−2.71
Met76	−1.66	−0.34	0.59	−0.20	−1.61
Met79	−1.85	0.51	−0.17	−0.20	−1.70
Phe98	−1.37	0.38	−0.52	−0.15	−1.67
Gln117	−0.85	−3.73	3.18	−0.10	−1.51
Thr214	−0.56	−2.45	0.42	−0.036	−2.63

Conclusions

MD simulation and binding free energy calculation have become an important and powerful tool. This tool is extremely effective in computational area applying to the study of receptor-ligand interactions (Wang et al., 2013). For the purpose of investigating the structural details on binding cortisol molecule to glucocorticoid receptor, molecular dynamic simulation was carried out for Apo and HCY systems. According to the MD simulation results, the structure of HCY system was more stable than that of the Apo system. Hydrogen bond analysis identified a couple of residues which were responsible for the binding of cortisol molecule. Gln45 was an important residue, of which the spatial position changed a lot from Apo system to HCY system. It is believed that such position changes may be related to the stress conduction. The results of binding free energy calculation were in good agreement to experiment data (Kolodkin et al., 2013). The energy decomposition analysis identified certain crucial residues which were responsible for the binding of the cortisol molecule. Based on the results of electrostatic potential surface and MM-PB/SA calculations, it is of great significance to take the electrostatic properties of cortisol molecule into consideration. Our research gives clear insights from atomic-level into the structural–functional aspects of cortisol molecule and its glucocorticoid receptor, and also provides valuable information for the design of drug which can treat stress related illnesses.