Computational Analysis of Residue-Specific Binding Free Energies of Androgen Receptor to Ligands

Androgen receptor (AR) is an important therapeutic target for the treatment of diseases such as prostate cancer, hypogonadism, muscle wasting, etc. In this study, the complex structures of the AR ligand-binding domain (LBD) with fifteen ligands were analyzed by molecular dynamics simulations combined with the alanine-scanning-interaction-entropy method (ASIE). The quantitative free energy contributions of the pocket residues were obtained and hotspot residues are quantitatively identified. Our calculation shows that that these hotspot residues are predominantly hydrophobic and their interactions with binding ligands are mainly van der Waals interactions. The total binding free energies obtained by summing over binding contributions by individual residues are in good correlation with the experimental binding data. The current quantitative analysis of binding mechanism of AR to ligands provides important insight on the design of future inhibitors.

Understanding of protein-ligand interaction and the quantitative characterization of binding affinity are very important for the discovery, design, and development of drugs (Cheng et al., 2007;Gilson and Zhou, 2007;Jorgensen and Thomas, 2008;Sun et al., 2013;Hu et al., 2020). Although experiments can study the thermodynamic properties of protein-ligand binding, determination of binging affinity is time-consuming, laborious and expensive. Computationally, the molecular mechanics generalized born surface area (MM/GBSA) method is often used to calculate the free energy of protein ligand binding (Massova and Kollman, 1999;Kollman et al., 2000;Moreira et al., 2007;Genheden and Ryde, 2015). Recently we developed a method called interaction entropy (IE) for practical and efficient calculation of entropy in protein-ligand and protein-protein binding (Duan et al., 2016). This method has been used in combination with alanine scanning (Massova and Kollman, 1999) (AS) and MM/GBSA to obtain the residue-specific contribution of each pocket residue (ASE method) (Yan et al., 2017;Liu et al., 2018;Qiu et al., 2018;Zhou et al., 2018;He et al., 2019;Yang et al., 2019;Wang et al., 2020).
In this study, fifteen AR ligands were analyzed with the ASIE method to quantitively characterize the detailed protein-ligand interactions, and the contribution of key binding residues on AR were identified. In addition, there is a strong correlation between the sum of the contributions of residues and the experimental binding free energy.

Molecular Dynamics Simulations
In this study, 15 androgen receptor systems with experimental ki/kd values and complex structures in the Protein Data Bank (Berman et al., 2000) (PDB) were used for MD simulation (Hamann et al., 1998;He et al., 2004;Bohl et al., 2005;Salvati et al., 2005;Sun et al., 2006;van Oeveren et al., 2006;Ostrowski et al., 2007;Bohl et al., 2008;Nirschl et al., 2009;Saeed et al., 2016). The systems were simulated using pmemd.cuda (Case et al., 2005) in AMBER18 (Pearlman et al., 1995;Case et al., 2019) with the ff14SB force field (Maier et al., 2015). For each system, TIP3P (Jorgensen et al., 1983) water model and periodic boundary conditions were used to solvate the complex, and the minimum distance between solute atoms and periodic boundary was set to 12 Å. Furthermore, Sodium and chloride ions were added to neutralize the system. A two-step minimization process was carried out, where only hydrogen atoms were optimized in the first step, and all atoms were optimized in the second step. The system was then slowly heated to 300 K with Langevin dynamics temperature regulation, followed by an equilibration of 500 ps Finally, 10-ns MD simulations were carried out in an NPT ensemble and 25,000 snapshots were saved for further analysis. Five independent replicates were simulated for each system, and the final results were obtained by averaging the calculated values on the five trajectories.
where the gas-phase: and solvation components: In the IE (interaction entropy) approach (Duan et al., 2016), the gas-phase component is computed by Where E x int and E a int contain the electrostatic and van der Waals energies between the ligand and residues x and Ala. The exponential average was evaluated by discrete time averaging.  After that, Eq. 2 becomes: The solvation free energy was calculated by the MM/GBSA method: The polarization part ΔG gb is obtained by the generalized Born (GB) model. The non-polar term ΔG np can be obtained by using the solvent accessible surface area (SASA) formula: Finally, the free binding energy of protein ligands can be expressed as Zhou et al., 2018) For each system, 25,000 frames were extracted from the entire 10-ns trajectory at an interval of 400 fs for IE calculation. Hundred frames were uniformly extracted from the 25,000 frames for MM/GBSA calculation with "igb" set to 8 (Ryckaert et al., 1977;Nguyen et al., 2013) because igb 8 is the latest GB model, and Nguyen et al. proved that the GB-Neck2 model (igb 8) shows significant improvement in solvation energy and effective radii calculation as compared to GB-OBC (igb 2, igb 5) and GB-Neck (igb 7) (Jorgensen et al., 1983). Following  our previous protocol, the dielectric constant for nonpolar, polar, and charged residues are 1, 3, 5, respectively (Hou et al., 2011;Petukh et al., 2015). We also calculated binding energy with the conventional MM/GBSA method for comparison, where the dielectric constant was set to one for all residues. The free energy and its standard deviation were obtained from five free energy values calculated from five independent trajectories.

Hot-Spots Residues in the Complex Structures
We used the ASIE method to quantitatively analyze the contribution of each pocket residue on AR when binds to different ligands. Figure 2 shows the specific binding free energy values of the pocket residues in the 15 systems. First of all, 704LEU, 707LEU, 745MET, 749MET and 764PHE contributed much more to the binding free energy than other residues in these 15 complexes, which are identified as hotspot residues. Secondly, there are clear differences in the contribution of residues to binding energy in these four types of ligands. For example, in cluster 2, residues 741TRP, 895MET, and 899ILE generally have higher binding energy contribution compared to those in cluster 3. In cluster 3, residues 701LEU, 742MET, 780MET, 787MET, 873LEU, 876PHE, and 880LEU have a more prominent binding energy contribution compared to cluster 2. 3g0w (Nirschl et al., 2009) and 1i37 (Sack et al., 2001) from cluster3 and cluster4 are selected for detailed analysis of hot spot residues, because these two systems have the strongest experimental value of binding free energy (Bohl et al., 2008;Nirschl et al., 2009), and the ligands belong to non-steroidal and steroid respectively. In 1i37, 704LEU, 873LEU, 764PHE, 742MET, 745MET, 780MET, 749MET, 787MET, 877THR, 741TRP, and 895MET contribute more than 1 kcal/mol to the total binding energy and are identified as hot residues ( Table 1).
The crystal structure of 1i37 shows that, LEU704, which is the major contributing residue, forms two Alkyl interactions with the ligand dihydrotestosterone at a distance of 4.8 Å and 5.1 Å ( Figure 3A).
Next, for the different mechanism of those four types of ligands binding to androgen receptor, we analyzed cluster2 and cluster3, among which the most obvious difference was found. As can be seen from Figure 2, residues 741TRP, 895MET and 899ILE in cluster2 has more binding free energy contribution than those in cluster3. Figure 4A,B respectively show the relative positions of these residues and ligands in the crystal structure of 3b67 (Cluster2) and 3g0w (Cluster3). It can be seen that the position of ligands in 3b67 is closer to these three  residues than that in 3g0w. In cluster3, 701LEU, 742MET, 780MET, 787MET, 873LEU, 876PHE and 880LEU generally have a stronger binding free energy contribution than those in cluster2. In addition, Figure 4C,D suggest that these residues are generally closer to the ligand in 3g0w (cluster3), which is in good agreement with our calculation results. The quantitative analysis of these residues specific binding energies provides important clues to design high affinity AR LBD ligands. As shown in Figure 2, the substituents of 2ax6 are shorter than other ligands in Cluster2, AR LBD has fewer residues to interact with it, and its binding free energy is also the lowest in Cluster2. 704LEU has a strong binding free energy contribution in all 15 systems, and in 2hvc, 704LEU has a very high binding energy contribution. This is due to its interaction with the benzene ring on the 2hvc ligand, and the trifluoromethyl on the ligand also interacts with it. Due to its polycyclic structure, cluster4's two steroidal compounds can interact with more residues in AR LBD than the other three types of molecules, which is also of certain reference significance for the design of AR LBD ligand.

The Total AR-Ligand Binding Energy
We next calculate the total binding energy by summing up the free energy contribution of each residue and compare the results with those obtained using the conventional MM/GBSA method. Table 3 shows the contributions of each energy to the binding free energy of the 15 systems, including enthalpy and entropy. Furthermore, enthalpy is decomposed into van der Waals, electrostatic energy, GB, and nonpolar solvation energy, and it is clear that van der Waals interaction provided most of the enthalpy of each systems. Table 4 shows the computational details of the fifteen systems, including the enthalpy, entropy and binding free energy calculated by ASIE method, the   Figure 5A shows the correlation between the binding free energy obtained from experiments and the binding free energy calculated by ASIE method. The calculated binding energy shows good correlation with the experimental values (R 0.85), although the experimental values range is very narrow. Figure 5B shows the correlation between the enthalpy calculated by the MM/GBSA method and the experimental binding free energy (R 0.34). Furthermore, even if the system with the largest error (2HVC) was removed from the MM/GBSA results, the correlation (0.66) between the calculated value and the experimental value was still lower than that of ASIE.

CONCLUSION
Androgen receptor (AR) is an important target for many diseases. In our study, we used the ASIE method to quantitatively analyze the contribution of residues when AR binds to different ligands. The residues that contribute most to each ligand are determined. Furthermore, the sum of the contributions of each residue was relatively consistent with the experimental binding energy values. The results of these calculations will be useful for the design and analysis of AR ligands.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
GS performed MD simulation JB helped with ASIE coding XP helped with system analysis XH participated in design of the project and discussion YQ help wrote the paper JZ supervised the overal project.

FUNDING
This work was supported by National Key R&D Program of China (Grant Nos. 2016YFA0501700 and 2019YFA0905201), the National Natural Science Foundation of China (Grant Nos. 21922301, 21761132022, 21673074, 91753103, and 21933010), and the Natural Science Foundation of Shanghai Municipality (Grant No. 18ZR1412600).