Average geometrical values (Å) over the last 50 ns of 100-ns-long MD simulations.
The selective α1-adrenergic receptor antagonist doxazosin is used for the treatment of hypertension. More recently, an experimental report demonstrated that this compound exhibits antiproliferative activity in breast cancer cell lines with similar inhibitory activity to gefitinib, a selective inhibitor of EGFR in the active state (EGFRAC). This experimental study provided evidence that doxazosin can be employed as an anticancer compound, however, the structural basis for its inhibitory properties is poorly understood at the atomic level. To gain insight about this molecule, molecular dynamics (MD) simulation with the molecular mechanics generalized Born surface area (MMGBSA) approach was employed to explore the structural and energetic features that guide the inhibitory properties of doxazosin and gefitinib in overexpressing EGFR/HER2 cell lines. Our result suggest that doxazosin exerts its inhibitory properties in breast cancer cell lines by targeting EGFR/HER2 but mainly HER2 in the inactive state (HER2IN), whereas gefitinib by targeting mainly EGFRAC, in line with previous literature. Decomposition of the binding affinity into individual contributions of HER2IN-doxazosin and EGFRAC-gefitinib systems detected hot spot residues but also showed polar interactions of Met801/Met793 with the quinazoline ring of both compounds. Principal component (PC) analysis revealed that the molecular recognition of the HER2IN-doxazosin system was linked to conformational changes but EGFRAC-gefitinib was not.
- MD simulations
Human epidermal growth factor receptor 1 (EGFR) and 2 (HER2) form part of a family of human epidermal grown factor receptors (EGFRs), and whose phosphorylation impacts cell proliferation, differentiation, and migration . The cytoplasmic tyrosine kinase domain (TKD) is considered one of the most studied receptors for developing new anticancer drugs . Activation of EGFR starts the molecular recognition of endogenous growth factors at the extracellular domain that, at the time, promotes the formation of homo- and heterodimers among the different members of EGFR, with HER2 the preferred member of EGFRs to form heterodimers [3, 4, 5]. The transition from a monomeric to dimeric state in EGFR is coupled to a conformational change in the TKD from an inactive to active state [6, 7, 8], whereas that for HER2 transitions from inactive, intermediate, and inactive states [9, 10, 11]. Generally, the signaling activity regulated by EGFR/HER2 is under control, however, mutations in TKD give place to constitutive activation of these receptors, which results in the development of different types of cancer, such as lung  and breast cancer . In addition, overexpression of EGFR/HER2 also happens with radiotherapy and chemotherapy resistance [14, 15, 16].
Based on the ability of tyrosine kinase inhibitors (TKIs) to inhibit EGFR, they can be divided in two types, those targeting the active state, such as Iressa, and those targeting the inactive state, such as erlotinib and lapatinib [17, 18, 19, 20]. Lapatinib showed dual activity on EGFR/HER2 [21, 22, 23, 24, 25, 26]. Despite the benefits of using these TKIs, the employment of them has been linked to severe side effects and drug resistance [27, 28, 29, 30]. Therefore, it is necessary to identify new compounds, either through drug design or drug repurposing, that target EGFR and/or HER2 receptors and are effective for cancer therapy. In this context, the combination of docking and molecular dynamics (MD) simulations has been widely exploited to generate new information about the binding properties between natural or synthetic TKIs and EGFR/HER2 [10, 11, 19, 20, 31, 32, 33, 34, 35, 36]. In a previous study, Hui et al. explored the inhibitory properties of doxazosin, an α-1 antagonist used for the treatment of hypertension, in two human breast cancer cell lines: BCC MDA-MB-231 and MCF-7 cells . MDA-MB-231 and MCF-7 cells are estrogen receptor (ER) positive and ER negative, respectively , and both cell lines also expressed EGFR and HER2; however, MDA-MB-231 expressed both receptors in higher concentrations than MCF-7 . Although EGFR and HER2 are important regulators for normal cellular processes, their dysregulation has been associated to protein overexpression that leads to the development of different types of cancer [1, 5]. They demonstrated that doxazosin induces apoptosis in breast cancer cell lines similar to Iressa (Gefitinib), reducing phosphorylated EGFR by a mechanism that does not involve the α1-adrenergic receptor, however, the structural and energetic basis for its inhibitory properties is poorly understood. In addition, Sharkawi et al. identified similar experimental antiproliferative activity of doxazosin in an MCF-7 cell line through the inhibition of EGFR . Thus, more robust structural and energetic studies are required to provide structural insight into the affinity of doxazosin for EGFR/HER2 compared with gefitinib. Structural data, docking, and molecular dynamics (MD) simulations combined with the MMGBSA approach were used to elucidate the molecular mechanism through which doxazosin and gefitinib inhibit EGFR/HER2.
2.1 Structural modeling
The free forms of EGFR in the inactive (EGFRIN) and active (EGFRAC) states were taken from the crystallographic structures of EGFRIN (PDB entry 1XKK) and EGFRAC (PDB entry 1 M17) conformations. The free forms of HER2 in the inactive (HER2IN) and active (HER2AC) states were taken from previous MD simulation studies; HER2IN  and HER2AC  conformations. Amino acid residues missing in the electron density map of EGFR structures were built with MODELER Version 9.14 .
2.2 Docking studies
Docking calculations were carried out using AutoDock 4.2 and AutoDock Tools 1.5.6 software . The ligand structures were built and optimized with the Gaussian package . The initial geometries of ligands were optimized at the AM1 level. Hydrogen atoms were added to ligands and receptors, and Kollman and Gasteiger partial charges were assigned for ligand and proteins, respectively. The affinity grid maps were constructed on the receptor using a grid size of 70 × 70 × 70 Å and 0.370 Å of spacing. Due to the stochastic nature of the Lamarckian algorithm, 20 runs were performed for each compound, and 30 conformations of the ligand (binding poses) were observed between ligand and protein. The best binding poses were selected using the criteria of having the lowest energetic conformations at the receptor binding site.
2.3 Molecular dynamics simulations
The protein-ligand results obtained by docking were checked through MD simulation studies. MD simulations were carried out using the AMBER16 package , in conjunction with the ff14SB force field . The systems simulated were put into a space-filling dodecahedric box of 12 Å, solvated with TIP3P water model , and neutralized with sodium and chloride ions (0.10 M) to create a physiological concentration. The parameterizations of the ligands were performed assigning AM1-BCC atomic charges and matching the atoms with the general Amber force field (GAFF) . Once the systems were constructed, they were minimized using steepest descent with position restraint of the ligands, followed by steepest descent without position restraint and conjugate gradients. The minimized systems were then submitted to 100 ns-long MD simulations using an NPT ensemble with the velocity rescaling arrangement to simulate a constant temperature at 310 K. A constant temperature and pressure (1 atm) were maintained using the weak-coupling algorithm , with coupling constants τT and τP of 1.0 and 0.2 ps, respectively. The electrostatic term was described by the PME method , and a 10 Å cut-off was selected for the van der Waals interactions. The time step for the MD simulations was set to 2.0 fs. The SHAKE algorithm  was employed to reset bonds to their right lengths after an unconstrained update. The conformations obtained from MD simulations at intervals of 20 picoseconds (ps) were analyzed using the cpptraj tool in Amber16. Plots of variation of root mean squared deviation (RMSD) and radius of gyration (RG) were generated to evaluate convergence. Clustering analysis using a cutoff of 2.5 Å was performed to identify the most populated conformation in the simulation. Principal components (PC) analysis along the most essential eigenvectors was carried out to evaluate total flexibility. A map of interactions was generated using Maestro Version 10.1, 2015–1 .
2.4 Affinity prediction and per-residue decomposition
The binding free energy (Δ
3. Results and discussion
3.1 Convergence and equilibrium
The stability of the evaluated systems was observed by measuring two geometrical parameters. The root mean squared deviation (RMSD) and the radius of gyration (RG) were determined to identify the time at which the systems reached convergence (Table 1). RMSD analysis showed that free and bound EGFR/HER2 systems reached stability between 20 to 50 ns with RMSD values, which oscillated between 1.40 and 4.20 Å. RG examination revealed that free and bound EGFR/HER2 systems exhibited stability from 20 to 50 ns with values oscillating between 18.8 and 20.2 Å. Based on this result, further analysis was carried out discarding the first 50 ns.
|EGFRAC||2.1 ± 0.20||19.0 ± 0.12|
|EGFRAC-doxazosin||1.7 ± 0.10||19.1 ± 0.10|
|EGFRAC-gefitinib||2.2 ± 0.20||19.2 ± 0.10|
|EGFRIN||1.8 ± 0.20||18.8 ± 0.10|
|EGFRIN-doxazosin||2.2 ± 0.20||19.0 ± 0.10|
|EGFRIN-gefitinib||2.7 ± 0.21||19.0 ± 0.10|
|HER2AC||3.6 ± 0.17||20.0 ± 0.10|
|HER2AC-doxazosin||1.7 ± 0.20||20.0 ± 0.14|
|HER2AC-gefitinib||1.4 ± 0.20||19.9 ± 0.10|
|HER2IN||3.9 ± 0.40||20.0 ± 0.01|
|HER2IN-doxazosin||3.4 ± 0.40||20.2 ± 0.10|
|HER2IN-gefitinib||4.2 ± 0.22||19.6 ± 0.13|
3.2 Structural analysis of complexes between doxazosin and EGFRAC/EGFRIN
To explore the structural differences between doxazosin and gefitinib on EGFR/HER2, the most populated receptor-ligand conformations were retrieved over the equilibrated simulation time (the last 50 ns) through clustering analysis. Analysis of the complex between doxazosin and EGFRAC showed that the ligand was stabilized through van der Waals interactions with Phe723, Val726, Leu792, Met793, Pro794, Phe795, Cys797, Leu799, and Leu844, and polar interactions with Gly719, Lys745, Gly796, Asp800, Asp837, Arg841, Asn842, and Asp855 (Figure 1A). In contrast, Val718 formed both van der Waals interactions and a hydrogen bond with the quinazoline ring of doxazosin. In the complex with EGFRIN, dozaxosin was bound through van der Waals interactions with Phe723, Val718, Val726, Ala743, Leu792, Phe795, Tyr801, and Leu844. The polar contact took place with Lys745, Gln791, Gly796, Glu804, and His805. Met793 formed both van der Waals interactions and one hydrogen bond with the benzodioxin moiety of doxazosin. Pro794 established both van der Waals contacts and one hydrogen bond with the quinazoline ring of doxazosin, whereas Glu804 formed hydrogen bonds with the quinazoline ring of doxazosin (Figure 1B). Stabilization of doxazosin did not establish interactions with Thr790 and Met766, two residues whose mutations have been linked to EGFR drug resistance [57, 58]. In addition, the characteristic interactions between Met793 and the quinazoline moiety were not observed, which has previously been observed for other TKIs of EGFR .
3.3 Structural analysis of complexes between doxazosin and HER2AC/HER2IN
Dozaxosin in complex with HER2AC was bound through van der Waals interactions with Leu726, Val734, Ala751, Phe864, Leu852, and Leu800 and polar interactions with Gly727, Ser783, Thr798, Gln799, Arg849, Thr862, Asp863, and Lys921 (Figure 1C). For the complex between doxazosin and HER2IN, the ligand was stabilized by Val754, Leu755, Met774, Leu785, Leu796, Pro802, Cys805, Leu852, and Phe864 and polar contacts with Lys753, Arg756, Gln799, Thr862, and Asp863. Tyr803 and Met801 formed van der Waals and hydrogen bonds with polar groups of the quinazoline ring, whereas Ser783 and Thr798 formed polar contacts with the linker between piperazine and the benzodioxin moiety (Figure 1D). Structural comparison of the complexes of doxazosin with HER2AC/HER2IN showed that doxazosin was better coordinated on HER2IN than HER2AC through more well-adjusted types of van der Waals and hydrogen bonds. In addition, the characteristic hinge hydrogen bond between Met801 and the polar atoms of the quinazoline moiety of several TKIs [31, 59] was present only in for the complex with doxazosin and HER2IN.
3.4 Structural analysis of complexes between gefitinib and EGFRAC/EGFRIN
Analysis of complexes between gefitinib and EGFRAC illustrated that the ligand was bound through van der Waals interactions with Leu718, Val726, Ala743, Met766, Leu788, Ile789, Leu792, Pro794, Phe795, Cys797, Leu844, and Phe856. The polar interactions were through contacts with Gly719, Ser720, Gly721, Lys745, Glu762, Thr790, Gln791, Gly796, Arg803, Thr854, and Asp855. Met793 formed van der Waals interactions and one polar interaction with the quinazoline ring of gefitinib, whereas Asp800 formed a hydrogen bond with one of the substituents at the quinazoline ring (Figure 2A). In complex with EGFRIN, gefitinib was stabilized by van der Waals contacts with Phe723, Val718, Val726, Ala743, Cys775, Leu792, and Cys797. The polar contacts were through Gly719, Ser720, Gly721, Gly724, Thr790, Gln791, Arg841, Asn842, Thr854, and Asp 855. Met793 established both van der Waals interactions and one hydrogen bond with the quinazoline ring, whereas Lys745 made a hydrogen bond with one of the substituents of the quinazoline ring (Figure 2B). Comparison of the map of interactions of both complexes showed that gefitinib was better coordinated on EGFRIN than EGFRAC. In both complexes, gefitinib established interactions with Thr790, a residue whose mutation is linked to EGFR drug resistance [57, 58]. In addition, in both complexes, the characteristic interactions between Met793 and the quinazoline moiety of ligand were observed, which has been reported elsewhere .
3.5 Structural analysis of complexes between gefitinib and HER2AC/HER2IN
Gefitinib in complex with HER2AC was bound through van der Waals interactions by Leu726, Val734, Leu800, Cys805, Leu807, Val851, and Leu852 (Figure 2C). Polar interactions were stabilized by Gly727, Ser728, Gly729, Asp808, Asp845, Arg849, Asn850, Thr862, and Gly865 residues, whereas Asp863 formed a hydrogen bond with one of the substituents of the quinazoline ring (Figure 2C). Gefitinib formed a complex with HER2IN, coordinated by Leu726, Val734, Val731, Cys805, Leu807, Leu866, Ala867, Leu869, Tyr923, Ala920, and Pro922 residues through van der Waals interactions. Polar interactions took place by Gly727, Ser728, Lys753, Arg811, Arg849, Asp863, and Lys921 residues, whereas Asp808 formed a hydrogen bond with one of the quinazoline ring substituents (Figure 2D). Structural comparison of both systems depicted that gefitinib was better stabilized on HER2IN than HER2AC. In addition, in both complexes, the characteristic polar interaction between Met801 and polar atoms of the quinazoline moiety of ligand was not observed [31, 59], as observed for the complex between doxasozin and HER2IN (Figure 1C).
3.6 Binding free energy
Determination of the Δ
|EGFRAC-doxazosin||−43.52 (5.13)||22.33 (4.91)||−2.69 (0.80)||−5.46 (0.50)||−29.33 (6.14)|
|EGFRIN-doxazosin||−38.78 (3.90)||3.25 (0.91)||9.68 (2.80)||−4.62 (0.37)||−30.46 (4.20)|
|HER2AC-doxazosin||−36.77 (6.0)||40.15 (12.0)||−21.06 (11.0)||−4.36 (0.55)||−22.05 (4.60)|
|HER2IN-doxazosin||−48.93 (4.0)||23.99 (9.0)||−7.39 (1.50)||−5.99 (0.30)||−38.32 (4.0)|
|EGFRAC-gefitinib||−55.01 (0.16)||−11.11 (0.65)||−25.40 (0.60)||−7.16 (0.01)||−47.88 (0.17)|
|EGFRIN-gefitinib||−44.10 (0.15)||27.00 (0.69)||−16.48 (0.65)||−5.74 (0.01)||−39.32 (0.23)|
|HER2AC-gefitinib||−34.76 (4.0)||38.46 (16.0)||−22.98 (5.0)||−4.90 (0.70)||−24.18 (5.0)|
|HER2IN-gefitinib||−38.78 (5.7)||−31.92 (13.0)||49.70 (1.50)||−4.90 (0.50)||−25.91 (6.0)|
3.7 Decomposition of the per-residue free energy
This analysis identified the residues that contributed the most to the Δ
Table 4 shows that Leu726, Gly727, Ser728, Val734, Asp863, and Phe864 contributed the most to the Δ
3.8 Principal component analysis
We evaluated the differences in mobility for the free and bound EGFR/HER2 systems via PC analysis. Evaluation of the quantification of the diagonalized covariance matrix based on covariance showed the following values: HER2AC, 15.0 nm2; HER2AC-doxazosin, 10.26 nm2; HER2AC-gefitinib, 8.83 nm2; HER2IN, 28.9 nm2; HER2IN-doxazosin, 18.21 nm2; HER2IN-gefitinib, 20.41 nm2; EGFRAC, 10.79 nm2; EGFRACAC-doxazosin, 9.46 nm2; EGFRAC-gefitinib, 7.69 nm2; EGFRIN, 8.17 nm2; EGFRIN-doxazosin, 11.55 nm2; and EGFRAC-gefitinib, 10.0 nm2. This result indicates that the molecular recognition of doxazosin or gefitinib on HER2AC, HER2IN, and EGFRAC decreased the number of conformational states compared to that of free HER2AC, HER2IN, and EGFRAC states. However, this conformational reduction was more significant for the free and bound HER2IN system. In contrast, a small increase in the conformational mobility was experienced upon the coupling of gefitinib by EGFRIN. This indicated that doxazosin and gefitinib binding to HER2IN was linked with reduced heterogeneity, which suggests that this molecular recognition was associated with an unfavorable entropy contribution that could contribute to decrease the favorable Δ
In this chapter, we explored the structural and energetic features that guide the similar inhibitory properties of doxazosin with gefitinib in overexpressing EGFR/HER2 cell lines combining docking and MD simulation with the MMGBSA approach. Based on these studies, we identified that doxazosin was able to target the active and inactive states of EGFR and HER2, however, its inhibitory activity against breast cancer cell lines was mainly by targeting HER2IN. Similarly, although gefitinib was able to target the inactive and inactive states of EGFR and HER2, its activity mainly targeted EGFRAC, in line with previous reports. Per-residue free energy analysis identified the key residues stabilizing HER2IN-DOX and EGFRAC-GEF systems, showing that in the stabilization of both systems, Met793 and Met801 were involved for EGFR and HER2, respectively. These residues stabilized HER2IN-DOX and EGFRAC-GEF systems through the formation of hydrogen bonds with the quinazoline ring, as reported for other TKIs. This study provides structural and energetic information that can be used to design new inhibitors for HER2IN or EGFRAC using doxazosin or gefitinib, respectively, as a pharmacophoric model.
This is an unpublished original article. The work was supported by grants from CONACYT (CB-A1-S-21278) and SIP/IPN (20201015).