1 Introduction

Hypertension is defined as a progressive cardiovascular syndrome arising from complex and interrelated etiologies. Its progression may lead to functional and structural cardiac and vascular abnormalities that damage the heart, kidneys, brain, vasculature, and other organs and lead to premature morbidity and death (Giles et al., 2005, 2009). According to the World Health Organization (World Health Organization, 2009) raised blood pressure is estimated to cause 7.5 million deaths, about 12.8% of the total of all deaths.

The Renin-Angiotensin-Aldosterone system is a major system that regulates blood pressure (BP) and fluid balance in the body. Renin, the first enzyme in the system, cleaves the inactive angiotensinogen which is secreted mainly by the liver and converts it to the decapeptide angiotensin I (Ang I). In a subsequent step, Ang I is converted to the octapeptide hormone angiotensin II (Ang II) by angiotensin converting enzyme (ACE). Ang II mediates vasoconstriction and BP is increased by activation of AT₁R in pathological states. Sartans [AngII receptor blockers - ARBs] constitute a major class of antihypertensive drugs by acting as Trojan horses to angiotensin II type 1 receptor (AT₁R) (Muszalska et al., 2014). Eight ARBs are marketed as commercial drugs: losartan, valsartan, irbesartan, candesartan, telmisartan, eprosartan, olmesartan and azilsartan. Olmesartan [IUPAC name: 5-(2-hydroxypropan-2-yl)-2-propyl-3-[[4-[2-(2H-tetrazol-5-yl) phenyl] phenyl] methyl] imidazole-4-carboxylic acid] is the active metabolite of the prodrug olmesartan medoxomil (Fig. 1).

Close

Figure 1 Chemical structures of olmesartan medoxomil, olmesartan and olmesartan methyl ether as well as their critical dihedral angles and the derivatives of olmesartan R239470, R781253, and R794847.

Export to PPT

Olmesartan exhibits potent inverse agonist activity and is the latest drug before azilsartan. It possesses a unique structure, bearing a hydroxyl group at the imidazole ring and a carboxyl group. For this reason it has received a considerable interest for its pharmacological and biophysical properties. In particular, studies have been already conducted in order to clarify olmesartan’s thermal, dynamic and structural properties in lipid bilayers and comprehend its distinct pharmacological profile (Ntountaniotis et al., 2014b, 2011). In addition, olmesartan’s impurities and degradation products were identified (Babu et al., 2010; Jain et al., 2012; Murakami et al., 2008; Venkanna et al., 2013). Recently, the crystal structure of the human AT₁R in complex with olmesartan was reported by Zhang et al. (2015a).

The in situ formation of olmesartan methyl ether from olmesartan, when the latter is dissolved in methanol was identified using various spectroscopic methodologies and its mechanism of formation was explored using theoretical studies (Ntountaniotis et al., 2014a). In the same study (Ntountaniotis et al., 2014a), olmesartan methyl ether demonstrated in vitro equal AT₁R inverse agonism with that of the prototype drug of this class of molecules, losartan.

It is reported that conformations based on NMR spectroscopy explain better the coherent biological activities of the drug molecules (Kellici et al., 2016; Musafia and Senderowitz, 2010). This statement is examined in the current study for the cases of olmesartan and its methyl ether analogue. The Nuclear Overhauser Effect (NOE) was utilized in this article to acquire distance constraints. The recent literature on the role of methyl groups in medicinal chemistry (Barreiro et al., 2011; Schönherr and Cernak, 2013) urged us to compare the conformational properties of olmesartan and its methyl ether analogue.

Miura et al. (2006) have reported that mutations in the AT₁R such as Y113F, K199Q, H256A and Q257A lower the binding affinity for olmesartan by 25-, 15-, 16- and 103-folds, respectively. Moreover the Y113A mutation did not exert any impact on the binding affinity of olmesartan. In marked contrast these mutations did not affect the binding of olmesartan to AT₁R in a recent study in which different experimental conditions and different expression systems of receptor have been used (Zhang et al., 2015a). These residues did not appear to develop direct interactions with the ligand in the cryo-cooled crystal structure of the AT₁R-olmesartan complex (Zhang et al., 2015a). In this article, docking and molecular dynamics (MD) studies for olmesartan have been performed in an effort to project a putative explanation of its biological effects and to interpret the impact of the mutations in its binding affinity. Crystallographic data of co-crystallized selective AT₁R antagonist ZD7155 with AT₁R (Zhang et al., 2015b) and the more recent crystal structure of AT₁R with olmesartan (Zhang et al., 2015a) were used for the in silico calculations. Induced fit docking calculations were also performed in two homology models of the AT₁R in order to evaluate whether these models provide a better explanation of the biological data. New theoretical approaches to study the molecular docking of drug molecules in the active site of the receptor (Avgy-David and Senderowitz, 2015; Sturlese et al., 2015) allow further accuracy in the results. In addition, comparisons have been conducted with derivative molecules of olmesartan with different biological function: olmesartan methyl ether, R239470 (a neutral antagonist (Miura et al., 2012, 2006)), R781253 (a new inverse agonist (Miura et al., 2012)) and R794847 (a weak partial agonist (Miura et al., 2012; Zhang et al., 2015a)). All these compounds were docked in the AT₁R and then MD simulations were performed in dipalmitoylphosphatidylcholine (DPPC) bilayers.

Finally, statistical treatment of experimental ¹³C NMR chemical shifts of olmesartan/olmesartan methyl ether and the calculated isotropic shielding constants of their respective conformers has been carried out. This was performed in order to assess whether calculated chemical shifts can predict which conformers are the major ones in the conformational equilibrium (Belostotskii, 2008; Cotos-Yáñez et al., 2015).

2 Materials and methods

Materials.Olmesartan was kindly provided by Daiichi Sankyo Pro Pharma. CD₃CN 99.80% was purchased from Euriso-Top, Gif-sur-Yvette, Essonne, France, and CD₃OD (99%+) from Sigma-Aldrich, St. Louis, MO.

NMR spectroscopy. The complete NMR data (¹H, ¹³C, COSY, HSQC, HMBC) have been assigned for both olmesartan (in CD₃OH) and olmesartan methyl ether (in CD₃CN). It must be pointed out that in case of olmesartan methyl ether, dissolved in CD₃OD, the signal of the methoxy group protons overlaps with the residual solvent peak. This necessitated the recording of 1D ¹H and 2D NOESY NMR experiments for olmesartan and olmesartan methyl ether in CD₃OD and CD₃CN solvents (see Supplementary Fig. S1).

2 mg of olmesartan was dissolved in CD₃OD and 2 mg of olmesartan methyl ether was dissolved in either CD₃OD or CD₃CN. The high-resolution NMR spectra were obtained on a Bruker Avance III spectrometer (Bruker Biospin GmbH, Rheinstetten, Germany) operating at 600.11 MHz for 1H and at 150.11 MHz for 13C, with a 5 mm inverse detection probe. The residual 1H signals (3.33 ppm for CD₃OD and 1.94 ppm for CD₃CN) were used as internal standards. 1D and 2D NMR spectra were acquired using the pulse sequences included in the Bruker library of pulse programs.

Experimental data were processed using MestReC software.

The NOEs signals observed (for olmesartan in CD₃OD and for olmesartan methyl ether in CD₃CN and CD₃OD) were quantified. The proton distances (Table 1) were calculated via the mathematical formula (1):

Software. All the computational studies were performed in the Schrödinger 2015.2 suite. The starting structures were sketched in Maestro 10.2 and subsequently ionized at a pH of 7.4 using Epik 3.2, as implemented in Ligprep 3.4 (“Schrödinger Release 2015-2: LigPrep, version 3.4, Schrödinger, LLC, New York, NY, 2015.”). Calculations were performed for both neutral and ionized forms. A schematic representation of the computational workflow is shown in Scheme 1.

Close

Scheme 1 Computational workflow.

Export to PPT

Conformational analysis. The Macromodel 10.8 module (“Schrödinger Release 2015-2: MacroModel, version 10.8, Schrödinger, LLC, New York, NY, 2015.”), was used for the conformational analysis of olmesartan and olmesartan ether. The dielectric constant (ɛ) was set at 35.85 for the simulations in acetonitrile and 33.62 for methanol (at 298 K). The conformational search was performed with distance constraints. The ten distances derived from the NOE (shown in Table 1) were used as constraints during the simulations for olmesartan in methanol. In addition, thirteen distances shown in Table 1 were used as constraints during the simulations for olmesartan ether in methanol. Fourteen distances were used as constraints for the simulation of olmesartan methyl ether in acetonitrile. To account for NOE-derived distance flexibility, a margin of ±0.5 Å was used for every constraint and a force constant of 10 kJ mol⁻¹ Å⁻². All the conformational searches were carried out using the torsional sampling method (a Monte Carlo multiple minimum method). A 21 kJ mol⁻¹ energy cutoff was used to remove the highest energy conformers. Every conformer was minimized using the Polak–Ribiere conjugate gradient (PRCG) method for a maximum of 2000 steps, to a convergence threshold for the gradient of 0.001 kJ mol⁻¹ Å⁻¹ and the OPLS3 force field (Harder et al., 2016; Shivakumar et al., 2010), since it provides accurate energy minimization potential functions for small molecules. In order to work with a manageable set of representative conformers, the resulting structures were clustered using the conformer_cluster.py script. The atomic RMSD matrix was calculated using all the heavy atoms, while the clustering was calculated using an average linkage method.

Geometry optimization. The lowest energy conformers of olmesartan in methanol and olmesartan methyl ether in methanol and acetonitrile were then subjected to further geometry optimization. The optimization was carried out using the Jaguar 8.8 module (Bochevarov et al., 2013; “Schrödinger Release 2015-2: Jaguar, version 8.8, Schrödinger, LLC, New York, NY, 2015.”). The DFT/B3LYP level of theory and the 6-31G^∗∗ basis set were used for the geometry optimization, while the Poisson-Boltzmann model (PBF) was used as solvent (Tannor et al., 1994). The SCF accuracy level was set to ultrafine. NMR chemical shieldings were calculated for all the representative conformers of olmesartan and olmesartan ether and the crystal structure of olmesartan using the GIAO method (Bochevarov et al., 2013; Cao et al., 2005). The chemical shifts were then calculated using the following formula:

Dihedral scan. In order to estimate the flexibility of the τ8 dihedral angle a relaxed conformational scan was performed using both molecular mechanics and ab initio methods. The coordinate scan was achieved in MacroModel 10.8 using the OPLS3 force field and setting the dielectric constant to 35.85 for acetonitrile and 33.62 for methanol. The structures were minimized using the Polak–Ribiere Conjugated Gradient method with a convergence threshold of 0.05. The initial coordinate was set at 0° and the final at 360° with an increment of 1°.

The relaxed potential energy surface (PES) scan was performed using the Jaguar 8.8 (Bochevarov et al., 2013; “Schrödinger Release 2015-2: Jaguar, version 8.8, Schrödinger, LLC, New York, NY, 2015.”) module of Schrödinger for the τ8 dihedral angle for both olmesartan and olmesartan methyl ether. The dihedral was incremented by 5° in each step. The DFT level of theory and the B3LYP/6-31G^∗∗ hybrid functional were used for the optimization of the structures during the coordinate scan for all the intermediate steps. The calculations were performed in vacuo and in methanol or acetonitrile solvents for olmesartan and olmesartan methyl ether respectively using the Poisson-Boltzmann finite element solvent model. Vibrational frequencies and IR intensities were calculated for all the intermediate conformations.

Docking procedure. The crystal structures and two homology models of the AT₁R were used for the docking procedure. The first model was kindly provided by Prof. T. Tuccinardi and was generated using the bovine rhodopsin crystal structure as template (PDB ID 1U19) (Tuccinardi et al., 2006). The second model, kindly provided by Dr. M. Matsoukas, was based on the CXCR4 crystal structure with PDB ID 3ODU (Matsoukas et al., 2013a, 2013b). The crystal structures (PDB IDs 4YAY (X-ray free electron laser (XFEL) structure with resolution 2.9 Å) (Zhang et al., 2015b) and 4ZUD (resolution 2.8 Å) (Zhang et al., 2015a)) and the homology models were prepared using the Protein Preparation Wizard (Madhavi Sastry et al., 2013), were ionized at a pH of 7.4 using PROPKA (Olsson et al., 2011), and were minimized using the OPLS3 force field (Shivakumar et al., 2010). Receptor grids for the homology models were generated centroid of the residues Tyr113, Lys199, Tyr253 and His256. The receptor grid for the crystal structure was generated centroid to the co-crystallized ligand ZD7155. The four important mutations according to Miura et al. (2006), for each 3D structure of the protein (Y113F, K199Q, H256A, and Q257A) were generated using the “Residue and Loop mutation” facility in Bioluminate 1.9. All residues within 5 Å of the mutation were minimized. Firstly, Schrödinger’s docking algorithm Glide was used to predict the binding geometries for the ligand-receptor complexes. The docking calculations were performed with GlideXP (extra precision mode) (Friesner et al., 2006; “Small-Molecule Drug Discovery Suite 2015-2: Glide, version 6.7, Schrödinger, LLC, New York, NY, 2015.”), using standard van der Waals scaling of 0.8. All the resulting poses were used as starting conformations for the QM-Polarized Ligand Docking approach (QPLD) (Cho et al., 2005; “Small-Molecule Drug Discovery Suite 2015-2: Schrödinger Suite 2015-2 QM-Polarized Ligand Docking protocol; Glide version 6.7, Schrödinger, LLC, New York, NY, 2015; Jaguar version 8.8, Schrödinger, LLC, New York, NY, 2015; QSite version 6.7, Schrödinger, LLC, New York, NY, 2015.”). In this workflow, accurate quantum mechanical charges are generated from the electrostatic potential energy surface of the ligand using the density functional theory (DFT), the B3LYP/6-31G^∗ basis set and the ultrafine SCF accuracy level within the Jaguar module. Ligands were redocked using GlideXP. The final selection was based on the Emodel score. Since Emodel combines the GlideScore and the internal energy of the ligand conformation, it is a more suitable score to rank the best poses of the same ligand. The five QPLD conformers with best Emodel scores were submitted as starting geometries for the Induced Fit Docking calculations with implicit membrane (Sherman et al., 2006a, 2006b; “Small-Molecule Drug Discovery Suite 2015-2: Schrödinger Suite 2015-2 Induced Fit Docking protocol; Glide version 6.7, Schrödinger, LLC, New York, NY, 2015; Prime version 4.0, Schrödinger, LLC, New York, NY, 2015.”). At the initial stage of docking side chains of residues that are within 5 Å of the ligand were trimmed. Three residues that are within 5 Å of the ligand and have the highest B-factors (above 40) were refined using Prime, version 4.0, Schrödinger, LLC, New York, NY, 2015. The ligands were redocked using the extra precision mode. The implicit membrane is a low-dielectric slab-shaped region, which is treated in the same way as the high-dielectric implicit solvent region. Hydrophobic groups, which normally pay a solvation penalty for creating their hydrophobic pocket in the high dielectric region, do not pay the same penalty while in the membrane slab. Conversely, hydrophilic groups lose any short ranged solvation energy from the high dielectric region when moving into the low dielectric region. The implicit membrane model is intended for use with proteins that span the membrane Prime, version 4.0, Schrödinger, LLC, New York, NY, 2015.

Molecular dynamics (MD) simulations. The poses in the two crystal structures with best Induced Fit Score were subjected to molecular dynamics simulations using Desmond 4.2 (“Schrödinger Release 2015-2: Desmond Molecular Dynamics System, version 4.2, D. E. Shaw Research, New York, NY, 2015. Maestro-Desmond Interoperability Tools, version 4.2, Schrödinger, New York, NY, 2015.”). The receptor was immersed in the DPPC bilayer. The SPC solvent model was used with an orthorhombic box shape. The OPLS3 force field was used during the building of the system (Shivakumar et al., 2010). The system was minimized using steepest descent for 2500 maximum iterations until a gradient threshold of 25 kcal mol⁻¹ Å⁻¹ was reached. The simulation time for the molecular dynamics was set at 300 ns in the NPT ensemble class. The temperature and the pressure were set at 325 K and 1.01325 bar correspondingly. The RESPA integrator was used with a time step of 2.0 fs (Tuckerman et al., 1992). The thermostat method used was the Nose-Hoover chain (Martyna et al., 1992), while the Martyna-Tobias-Klein method was used as barostat (Martyna et al., 1994). The cutoff radius was set at 9.0 Å. MD simulations were run for olmesartan in both crystal structures (4YAY and 4ZUD) and the four mutated structures, while for olmesartan methyl ether and the compounds R239470, R781253 and R794847 MD simulations were run in both crystal structures and repeated twice.

MM-GBSA calculations. The MM-GBSA calculations were performed using Prime. The water and DPPC molecules were removed from the trajectories using the delete_atoms.py utility. The binding energy was calculated for a total of 4151 frames of the MD trajectory starting from the 100th ns until the end of the trajectory at the 300th ns. The binding free energy was calculated by Eq. (3).

Prediction of the chemical shielding constants of olmesartan and olmesartan methyl ether in their binding sites. From the MD trajectories of olmesartan and olmesartan methyl ether in the receptor, a snapshot closest to the calculated average structure for each trajectory of the protein-ligand complex was created. Chemical shielding constants were calculated using QSite (Murphy et al., 2000). The QM region consisted of the ligands which were treated using the DFT/B3LYP/6-31G^∗∗ basis set. The rest of the protein was treated using molecular mechanics and the OPLS-2005 force field with the maximum number of iterations set to 100 cycles. All calculations were converged before this limit using the conjugated gradient algorithm. The NMR shifts were calculated from the shielding constants using the formula (2).

3 Results and discussion

3.1

3.1 Conformational analysis of olmesartan and olmesartan methyl ether in solution using NMR and theoretical calculations

¹H NMR spectra of olmesartan in CD₃OD and olmesartan methyl ether in CD₃CN have been previously reported (Ntountaniotis et al., 2014a). The ¹H NMR spectrum of olmesartan methyl ether in CD₃OD is illustrated in Supplementary Fig. S1. Their structural identification is shown in Supplementary Table S3. Distance NOE constraints are reported in Table 1 and were obtained after integration of the contour cross-peaks from the three 2D NOESY spectra.

The representative conformations of olmesartan in methanol, derived using NOE constraints described in Table 1 are shown in Fig. 2a. In the article by Musafia and Senderowitz (2010) it is stated that NMR derived structures were more appropriate for ligand based drug design. We have evaluated this hypothesis by applying conformational analysis on olmesartan and its methyl ether by considering the NOE data and comparing the resulted conformations with those obtained without implementing NOE constraints. We are aware that the use of NOE constraints for structure determination has limitations on the basis that flexible molecules adopt a plethora of conformations that are projected on the finally recorded NOEs. However, if all NOEs are counted together, this dynamic process of the various conformations is lost. As it can be observed by comparing the conformations in Fig. 2 and the ones in Figs. S2–S4, conformational analysis with and without NOE constraints provided almost identical results. Five of the eight representative conformations are syn; thus, the tetrazole and imidazole moieties adopt the same orientation with respect to the ring plane of phenyl A. Half of the representative conformations of olmesartan methyl ether, in both solvents, are anti; thus, the tetrazole and imidazole moieties adopt the opposite orientation with respect to the A phenyl ring plane (Fig. 2b and c). In order to understand which conformation (syn or anti) possesses a lower energy a relaxed PES scan was performed around the τ8 dihedral using both force field and ab initio methods (Figs. S5 and S6). Both techniques showed that for olmesartan in methanol the syn conformation is more stable (by 0.3 kcal/mol) while for olmesartan methyl ether the anti conformation is more stable (0.18 kcal/mol).

Close

Figure 2 Representative low energy conformations of (a) olmesartan in methanol and olmesartan methyl ether in: (b) methanol and (c) acetonitrile as derived from MacroModel using NOE derived constraints.

Export to PPT

The optimized low energy structure of olmesartan using NOE constraints is compared with the crystallized structure of olmesartan from ethanol by Yanagisawa et al. (1996) (CSD reference code: ZOGSOD). The RMSD between the crystal structure and the lowest energy conformer in olmesartan was found to be 0.2 Å (Fig. 3a and b). Another low energy structure was found for olmesartan which is the conformational enantiomer of the lowest energy conformer. The optimized lowest energy structure of olmesartan methyl ether in methanol as shown in Fig. 3c adopts an anti conformation as opposed to that of olmesartan that adopts a syn conformation in the same solvent. Olmesartan methyl ether in acetonitrile adopts a syn conformation (Fig. 3d). These results indicate that the syn-anti interconversion in the two solvents is achieved without a significant energy cost.

Close

Figure 3 (a) Crystal structure of olmesartan (CSD reference code: ZOGSOD); (b) low energy conformers of olmesartan in methanol optimized using the HF/6.31G^** basis set; (c) optimized lowest energy conformer of olmesartan methyl ether in methanol; (d) optimized lowest energy conformer of olmesartan methyl ether in acetonitrile.

Export to PPT

3.2

3.2 Induced fit docking calculations

Olmesartan docked in the crystal structure 4YAY: As it can be predicted, most interactions of olmesartan are developed with Arg167. Both, the tetrazole and the carboxyl moiety of the molecule develop hydrogen bonds with Arg167. Furthermore, cation-π interactions are developed between Arg167 and the A phenyl group of olmesartan. The hydroxyl group and the imidazole group of olmesartan develop hydrogen bonds with Tyr35 (Fig. 4a left). The imidazole group is also implicated in π-π interactions with Trp84. Polar interactions are also observed between the tetrazole group and Ser105, Ser109 and Lys199. The alkyl group develops hydrophobic interactions with Ile288, Phe77 and Tyr292. The RMSD between the docked structure of olmesartan in the crystal structure 4YAY and the crystallized structure of olmesartan in the AT₁R receptor (PDB ID 4ZUD) was found to be 0.3 Å. This pinpoints that the crystal structure 4YAY produces reliable poses and that the selected virtual workflow is suitable. As can be observed in Fig. S7a–d and Supplementary Table S1 the mutations Y113F, K199Q, H256A and Q257A do not modify the binding mode of olmesartan at the AT₁R receptor. Olmesartan in all these cases adopts an anti configuration. Occasionally, the tetrazole moiety develops hydrogen bonds with Ser105 (Y113F, K199Q and Q257A). Also the hydroxyl group occasionally develops hydrogen bonds with Thr88 (Y113F, H256A).

Close

Figure 4 Most energetically favorable poses of olmesartan (left) and olmesartan methyl ether (right) in the AT₁R after applying induced fit docking in the (a) crystal structure (PDB ID 4YAY), (b) CXCR4 homology model, (c) rhodopsin homology. Blue lines represent hydrogen bonds, and cyan lines represent π-cation interactions.

Export to PPT

Olmesartan in two homology models: Although the mutations reported by Miura et al. (2006) cause significant effects on the biological activity of olmesartan and were explained based on the bovine rhodopsin crystal structure, they do not result in important differences in the Induced Fit Docking scores. In their reported findings, olmesartan was manually docked in the model and then the system was minimized using the CHARMM force field. The tetrazole moiety was suggested to interact with Gln257, the carboxyl group to interact with Lys199 and the hydroxyl group to interact with Tyr113 (Miura et al., 2006). Miura et al. in a more recent model based on the CXCR4 structure (PDB ID 3OE0) (S. Miura et al., 2013) found that both His256 and Gln257 develop hydrogen bonds with the tetrazole moiety, Lys199 develops hydrogen bonds with both the tetrazole and carboxyl groups and Tyr113 develops a hydrogen bond with the hydroxyl group. In order to evaluate whether AT₁R homology models based on rhodopsin and CXCR4 explain in a more optimized way the impact of these mutations, IFD calculations were performed in two reliable models already published (Matsoukas et al., 2013a; Tuccinardi et al., 2006).

1. Olmesartan in the CXCR4 homology model

The most remarkable feature observed is the fact that the predicted binding pose of olmesartan in the homology model is very similar to the one that olmesartan adopts in the crystal structure. Olmesartan adopts the anti configuration. The tetrazole moiety develops hydrogen bonds with Arg167, Ser105, Ser109 and Lys199 and π-π interactions with Phe182. One of the differences with the crystal model is the crucial role that Lys199 plays in this model. Lys199 develops hydrogen bonding with the carboxyl group and cation-π interactions with the A phenyl group. Hydrogen bonds are also observed between the hydroxyl group of olmesartan and the carbonyl group of Pro192 and the imidazole group and Asp263 (Fig. 4b left). The mutations here also do not cause modification in the binding pose of olmesartan in the homology model (Fig. S8a–d). Even the mutation K199Q, that was expected to have a great impact on the binding mode due to the importance of Lys199 in this model, does not affect drastically the binding mode. A hydrogen bond is developed between Gln199 and the tetrazole moiety. The application of unconstrained IFD calculations did not reveal the pose reported by Miura et al. (2013). Since in this model Arg167 has been located in the ligand binding site, olmesartan consistently forms hydrogen bonds with it.

• Olmesartan in the rhodopsin homology model

In this homology model the residues Gln257, His256, Phe182, Tyr184, Thr260 and Asn200 play a crucial role while the residues Lys199 and Tyr113 are in the vicinity of the binding site. More specifically, the tetrazole moiety develops hydrogen bonds with the residues Gln257, Thr260 and Asn200. The carboxyl group develops hydrogen bonds with Tyr184 and Phe182. The imidazole and hydroxyl groups develop hydrogen bonds with His256 (Fig. 4c left). The mutations Y113F and K199Q do not modify the binding mode of olmesartan or the Induced Fit Scores (Fig. S6a and b). The mutation H256A seems to cause some changes in the interactions developed by olmesartan. The hydrogen bond between Gln257 and Ser109 is interrupted. The two hydrogen bonds developed by His256 (Fig. 4c) are substituted by one hydrogen bond with Ala256 (Fig. S9c). Changes are also observed with the mutation Q257A. The tetrazole moiety develops hydrogen bonds only with Asn200. In addition the carboxyl group of olmesartan develops hydrogen bonds with Phe182, Tyr184, Ser109 and Tyr113. Although the interactions seem to change for olmesartan in these two mutations, the binding conformation is always the anti one.

Olmesartan methyl ether in the crystal structure and the two homology models. Olmesartan’s methyl ether adopts a very similar conformation to that of olmesartan inside the crystal structure. The anti conformation is favored. Both the tetrazole and carboxyl moieties of the molecule form hydrogen bonds with Arg167. π-π interactions are developed between tetrazole and Phe182 and between the imidazole group and Trp84. The imidazole group also forms a hydrogen bond with Tyr35 (Fig. 4a right). The methyl ether moiety of olmesartan develops hydrophobic contacts with the residues Pro285, Tyr35, Ile31 and Ile288 (Supplementary Table S2).

The homology model based on CXCR4 predicts the anti conformation of olmesartan methyl ether. The residue Arg167 forms hydrogen bonds with both the tetrazole and the carboxyl groups. The tetrazole moiety develops multiple hydrogen bonds with Tyr87. The carboxyl group develops a hydrogen bond with Lys199. The imidazole group forms a hydrogen bond with His256 and the oxygen of the methyl ether group develops a hydrogen bond with Tyr292 (Fig. 4b right).

The homology model based on the rhodopsin structure predicts the formation of hydrogen bonds between tetrazole and Asn200 and Gln257 (Fig. 4c right). The carboxyl group is predicted to develop hydrogen bond with Tyr184. The imidazole moiety forms π-π interactions with His256 and Phe182.

3.3

3.3 Molecular dynamics simulations

MD simulations have been performed as it is well known that calculations are affected by the environment surrounding the studied system. The following reasons favoured the use of DPPC bilayers exclusively in this computational study: (a) DPPC’s partition coefficient with respect to its aqueous environment, especially in the fluid state, resembles that of natural plasma membranes of the vasculature (Netticadan et al., 1997; Oliveira et al., 2009); and (b) Phosphatidylcholines (PCs) are the most abundant lipid species in the plasma membranes of the vascular smooth muscle cells (Oliveira et al., 2009) and sarcolemma cardiac membranes (Netticadan et al., 1997). The RMSDs for the Cα of the receptor (in blue) and for the ligand (in red) are shown for the following: (i) olmesartan in the crystal structure (Fig. S10a), (ii) the four mutations (Fig. S10b–e) and (iii) olmesartan methyl ether in the crystal structure (Fig. S10f). As can be seen after 100 ns the protein’s structure has been stabilized.

Olmesartan’s interactions in the crystal structure 4ZUD: As can be observed in Fig. 5a olmesartan develops hydrogen bonds with Arg167 during the trajectory time of the simulation. Hydrogen bonds between Arg167 side chain and the carboxyl group of olmesartan occur 95% of the simulation time, while those with tetrazole occur 35% of the simulation time. Interactions such as π-cation occur for 89% of the time with the second phenyl group and 42% of the time with the tetrazole group. Ionic interactions are also developed with the carboxyl group for 41% of the simulation time. The other key residue in the interaction diagram is Lys199. Hydrogen bonds are formed between Lys199 side chain and two nitrogen atoms of the tetrazole group for 11% of the simulation time. Hydrophobic interactions such as π-cation are also developed between Ly199 and the first phenyl group. Water bridges are formed by Lys199 as well as ionic interactions with the negatively charged nitrogen atom in tetrazole. Another important residue is Trp84. It develops π-cation interactions with the positively charged nitrogen atom of the imidazole aromatic ring for 66% of the simulation time. The residue Tyr35 is mainly involved in water bridges establishing connections with the imidazole ring. The residue Phe182 develops mainly π-π interactions with the tetrazole moiety of the drug as well as water bridges that connect it with the carboxyl moiety. The residues Cys180, Tyr184, Tyr92 and Tyr87 are implicated in water bridges that allow them to interact with the carboxyl group (although these contacts are only observed for 10–15% of the simulation time importantly they occur periodically). The residues Leu81, Val108, Leu112, Tyr113, Met284, Ile288 and Tyr292 develop hydrophobic interactions with olmesartan.

Close

Figure 5 Protein-ligand interactions during the whole time of simulation for (a) olmesartan; and (b) olmesartan methyl ether.

Export to PPT

In order to better understand the preferred conformation of olmesartan in the receptor, the torsional angles of olmesartan were analyzed (Fig. 6). The dihedrals τ6, τ7, and τ8 constrain the molecule to an anti conformation during the whole time of simulation. This is probably due to three important interactions: (i) tetrazole with Lys199 and Arg167 that constrains the dihedral τ8, (ii) carboxyl group with Arg167 and (iii) imidazole group with Trp84 that forces the dihedrals τ6 and τ7 to adopt a narrow range of values.

Close

Figure 6 Modification of olmesartan’s dihedral angles as a function of time during MD simulations in the crystal structure 4ZUD surrounded by DPPC bilayers.

Export to PPT

The results obtained from the MD calculations have been compared with those obtained with the most important experimental NOEs (H9-H4 and H19-H12) (Fig. S13). The average values for the distance H9-H4 were found to be 4.07 Å and for the distance H19-H12 2.70 Å. These values are in accordance with those obtained from the NOE.

Olmesartan’s interactions in the mutated crystal structures Y113F, Y113A, H256A, Q257A: The interactions of olmesartan in the mutants Y113F, H256A and Q257A are shown in Fig. S11a, c and d. These three mutations do not seem to affect the conformation of the drug inside the active site or the interactions it develops with the key residues. More specifically, in these three mutations the most important interactions with residues such as Arg167, Lys199, Tyr35, Trp84, Tyr87, Thr88, Tyr92, Val108, Ile288 and Tyr292 remain intact. The same interactions were observed for the mutated structure Y113A (Fig. S12). For all critical dihedral angles that define olmesartan only τ8 has been useful for further elaboration. All other dihedral angles remained constant during the trajectory in both the wild type and the mutated form. In order to show this fact, the evolution of the critical dihedral τ8 for all the simulations is shown in Fig. 7b (Y113F), 7d (H256A), 7e (Q257A) and 7f (Y113A).

Close

Figure 7 Modification of τ8 during the trajectory for (a) olmesartan; (b) olmesartan in the Y113F mutant; (c) olmesartan in the K199Q mutant; (d) olmesartan in the H256A mutant; (e) olmesartan in the Q257A mutant; (f) olmesartan methyl ether and (g) olmesartan in the Y113A mutant.

Export to PPT

Olmesartan’s interactions in the mutated crystal structure K199Q: While in the previous cases the critical dihedrals remained the same, in the case of K199Q the values for the torsional angle τ8 change drastically (Fig. 7c). It is clear that this mutation causes the greatest effect in the ligand orientation inside the receptor. Seemingly, it stabilizes the position of the tetrazole moiety of the molecule through hydrogen bonding and π-cation interactions. Also the interactions that are lost due to the replacement of Lys with Gln are not balanced by any other residue. This shows that a change to this particular residue causes simultaneously a substantial change in the conformation of olmesartan and the loss of critical interactions that could explain the reduction in biological activity.

Olmesartan methyl ether in the crystal structure: Olmesartan methyl ether develops the same interactions as olmesartan (Fig. 5b). The carboxyl group forms hydrogen bonds with Arg167 for 90% of the simulation time, as well as ionic interactions. This moiety also develops water bridges with residues such as Tyr87, Cys180, and Phe182. The tetrazole group develops: hydrogen bonds with Arg167 (11% of the time) and Lys199 (13%), water bridges with Ser109 (14%), and π-cation interactions with Lys199 (47%) and Arg167 (37% of the time). The imidazole group interacts with Trp84 with π-π interactions (12%) and π-cation interactions (46%) and with Tyr35 mainly forming water bridges (25%). The residue Arg167 also develops π-cation interactions with the A phenyl group. The conformation that olmesartan methyl ether adopts in the binding site is always anti.

Also in this case, the distances H19-H12, H9-H4, H9-H5, H9-H23 and H10-H23 were monitored (Fig. S13). The average values were found to be 2.7 Å, 4.0 Å, 5.0 Å, 5.6 Å and 6.0 Å in accordance with the NOEs.

3.5

3.5 The interactions of R239470, R794847 and R781253 with the crystal structure 4ZUD

Zhang et al. (2015a) suppose that the carbamoyl moiety of R239470 forms additional hydrogen bond interactions with Tyr87 and Tyr92. In the MD simulation of R239470 a very stable intramolecular hydrogen bond is formed between the —NH₂ moiety of the carbamoyl group and the hydroxyl group of the molecule (Fig. 8a). This interaction seems to “lock” the structure in a slightly different position compared to the one that olmesartan adopts, and results in forming an additional stable interaction between the imidazole moiety and Tyr92 for 42% of the simulation time. The rest of the developed interactions are the same for both molecules of olmesartan and R239470. In the case of R794847 this intramolecular hydrogen bond does not occur. The —NH₂ group of the carbamoyl moiety forms a hydrogen bond with Tyr87 for 72% of the simulation time. The 4-hydroxybenzyl moiety develops interactions mainly with Thr260 (hydrogen bond for 81% of the simulation time), Lys199 (π-cation interactions for 68% of the simulation time), Gln257 (hydrogen bond for 20% of the simulation time) and Trp253 (π-π interactions for 15% of the simulation time) (Fig. 8b). On the other hand, the hydroxyl group of the 4-hydroxybenzyl moiety of R781253 forms stable hydrogen bonds with Gln257 (80% of the simulation time) and Trp253 (17% of the simulation time) and π-cation interactions with His256 (37% of the simulation time). The rest of the interactions that are developed by these molecules at their binding site on the AT₁R follow the same pattern of olmesartan and olmesartan methyl ether and need no further discussion (Fig. 8c).

Close

Figure 8 Interactions fraction and binding modes of (a) R239470; (b) R794847; (c) R781253.

Export to PPT

From the MM-GBSA calculations the most favored binding was exerted by R781253 showing highly favored VdW and Coulombic interactions (Table 2). The favored binding of R794847 is attributed also to Coulombic and VdW interactions. Olmesartan shows similar binding to wild type and mutated Y113F and K199Q receptors. However, it shows less favored binding with the mutated receptors H256A and Q257A. In H256A the coulombic interactions are disfavored and in Q257A the Solv_GB. Coulombic interactions are responsible for the low binding of olmesartan methyl ether.

Table 2 Prime MM-GBSA energies for olmesartan, binding at the wild- type AT₁R and its four mutated structures Y113F, K199Q, H256A and Q257A, olmesartan methyl ether, R781253, R794847 and R239470 binding at the wild type of AT₁R.

	ΔG Binding	Coulomb	Covalent	Hbond	Lipo	Packing	Solv_GB	vdW	Kd (nM)
Olmesartan Wild Type	−94.94	−8.82	11.33	−7.05	−37.72	−4.62	6.58	−54.63	2.2 ± 0.7 (Miura et al., 2006)
Olmesartan Y113F	−96.35	3.93	11.24	−7.22	−39.48	−5.36	−2.01	−57.44	52 ± 7 (Miura et al., 2006)
Olmesartan K199Q	−101.40	−2.38	9.92	−7.15	−40.46	−4.75	−1.55	−55.02	33 ± 6 (Miura et al., 2006)
Olmesartan H256A	−84.23	37.69	9.92	−7.11	−38.03	−5.28	−26.38	−55.02	36 ± 4 (Miura et al., 2006)
Olmesartan Q257A	−83.49	−16.18	8.35	−7.15	−31.52	−3.34	16.87	−50.51	226 ± 4 (Miura et al., 2006)
Olmesartan methyl ether	−82.75	37.72	8.54	−6.40	−40.30	−4.10	−22.14	−56.07	pIC50 = 8.44 ± 0.04a (Ntountaniotis et al., 2014a)
R781253	−121.48	−17.82	8.09	−7.87	−51.07	−3.32	20.57	−70.06	21 ± 10 (Miura et al., 2012)
R794847	−110.63	−51.99	8.91	−3.72	−47.04	−4.03	57.06	−69.80	48 ± 12 (Miura et al., 2012)
R239470	−80.03	9.78	7.22	−2.94	−31.84	−3.62	−8.57	−50.05	0.8 ± 0.3 (Miura et al., 2012)

Coulomb: Coulomb energy; Covalent: Covalent binding energy; vdW: Van der Waals energy; Lipo: Lipophilic energy; Solv_GB: Generalized Born electrostatic solvation energy; Hbond: Hydrogen-bonding energy; Packing: Pi-pi packing energy.

^a In reference Ntountaniotis et al. (2014a) is mentioned only the pIC₅₀ value.

4 Conclusions

In a previous publication we found that olmesartan in a methanolic solution even at low temperatures is converted to its methylated ether analogue. The in vitro experimental results showed that this methylated analogue possesses equipotent activity with AT₁R antagonist olmesartan. This finding triggered our interest to study the conformational properties both of olmesartan and its olmesartan methyl ether in the new co-crystallized AT₁R receptor with ZD7155. Indeed, the molecular docking and dynamics results can provide a plausible explanation for the in vitro results.

Miura et al. (2006) have reported that mutations in the AT₁R such as Y113F, K199Q, H256A and Q257A lower the binding affinity of olmesartan. It is possible that these residues do not contact the ligand and their mutation affects ligand binding indirectly, by altering the overall conformation of AT₁R. This could be supported by the fact that mutation of Y113F reduced olmesartan affinity, whereas the mutation Y113A substitution for this residue, did not affect ligand binding. Alanine mutation is a well-tolerated substitution that does not significantly affect receptor conformation (Cunningham and Wells, 1989). Alternatively the inconsistency between the results from mutagenesis studies and the crystal structure of AT₁R-olmesartan complex could be due to a possible conformational heterogeneity of Y113, K199, H256 and/or Q257 and to the fact that the crystal structure of a receptor represents only one single “snapshot” conformation. In specific, the crystal structure of the inactive β2 adrenergic receptor revealed the absence of an electrostatic interaction between the intracellular ends of TM3 and TM6, which has been shown in mutagenesis and molecular modeling studies to stabilize the inactive state of this receptor and in the inactive crystal structure of rhodopsin (Rasmussen et al., 2011). Dror et al. have suggested that the inactive β2 adrenergic receptor exists in equilibrium between conformations with the electrostatic interaction formed and broken, and that the crystal structure of this receptor represents one single conformation (Dror et al., 2011). Interestingly, a previous study has shown that even different copies of the same crystal structure of the D3 dopamine receptor had different structures, thus highlighting particular areas of conformational flexibility. In specific the second intracellular loop of D3 dopamine receptor formed a 2.5 turn alpha helix in one copy whereas in the other it was unstructured (Chien et al., 2010). The presence of multiple conformations of a single receptor is further supported by the fact that even a single agonist could stabilize different conformations, and thus activating different signaling pathways. For example, the corticotropin releasing factor-related agonist, sauvagine, stimulated the adenylate cyclase activity in a bell-shaped manner, which could be partially due to that at subnanomolar concentrations, the cAMP production was raised after activation of Gs, whereas at higher concentrations, the cAMP production was inhibited through the stimulation of Gi (Beyermann et al., 2007). The concept of the presence of different receptor conformations that bind differentially the same agonist and activate differentially various signaling pathways has also been introduced by Maggi and Schwartz, based on their experimental findings (Maggi and Schwartz, 1997). Different AT₁R conformations stabilized by arrestin types 1 and 2 could be responsible for the 9-fold difference in the binding affinities of SII-Ang II to these entities (Sanni et al., 2010). The complexity in the formation of different conformations of a receptor bound with the same ligand is further increased by the ability of ligands to change their conformations in their unbound and bound states. A previous study has shown that hydrophobic ligands unfold during binding (Perola and Charifson, 2004). More interestingly, Stockwell et al. have shown that three of the most ubiquitous biological ligands in nature, ATP, NAD and FAD bind to proteins in widely varying conformations (Stockwell and Thornton, 2006). This is further supported by the fact that the retinal bound to rhodopsin upon light absorption is transformed from the 11-cis form into the all-trans retinal. During this process the receptor is transformed to its active state metarhodopsin, passing through several intermediate conformations, such bathorhodopsin and lumirhodopsin.

The above discussion brings to the point that the crystal structure and homology models of AT₁R do not simulate adequately the real environment of AT₁R. AT₁R is incorporated in the lipid bilayer which may interfere in the binding of AT₁R antagonists. Mutation of AT₁R may break the communication of the lipid bilayer core and AT₁R and thus affects significantly the binding of the AT₁R. The absence of a real lipid bilayer environment in docking and MD calculations make these inaccurate to simulate the biological results. An effort is under progress to increase the complexity of the study systems in order to examine whether we can reconcile the discrepancy between the experimental and in silico results.

Prime MM-GBSA calculations for the five molecules studied (olmesartan, olmesartan methyl ether, R781253, R794847 and R239470) revealed similarities and differences in their binding energies and its components. Olmesartan methyl ether has less favoured ΔG_binding. This is attributed to its disfavored Coulombic interactions. The interactions of olmesartan in the crystal structure of AT₁R (PDB ID 4ZUD) did not provide a certain trend. For example the interactions of olmesartan at the wild type and mutated type K199Q do not differ significantly. R239470 shows higher ΔG due to disfavored Coulombic interactions. In contrast R781253 and R794847 have more favored Coulombic interactions than olmesartan and this results in a better binding energy in respect to olmesartan. Thus, their favored Coulombic interactions compensate successfully their Generalized Born electrostatic solvation energy.

The program Jaguar was used to compare the simulated chemical shifts obtained from different conformations with experimental NMR data. Such comparisons are useful to trigger the interest for using more accurate quantum mechanical calculations.

The ideal situation would be to study the interactions of drugs in a real environment. NMR methodologies applied also to AT₁R even in its crystal structure would be also very beneficial for the advancement of the field (Unione et al., 2014). There are currently in the literature several discrepancies, especially in the field of membrane proteins, on defining their real structure (as was the case of the EmrE protein). Herein, we illustrated that a ligand-based approach can be used and is able to lead to logical results on the interaction of a ligand with its GPCR. Indeed, methods such as paramagnetic relaxation enhancement (PRE) and residual dipolar couplings (RDCs), could be utilized; however, as has been capitalized in a recent (2016) perspective in Nature Structural and Molecular Biology (Liang and Tamm, 2016) “although membrane proteins constitute more than half of all drug targets 50, saturation transfer difference (STD) so far has seen only limited use for membrane proteins, probably because of the complicating presence of the membrane mimetic.” In the very same paper it was pinpointed that “nonphysiological constraints imposed by the crystal lattice rather than the micelle environment are problematic for these proteins.” The associated problems due to the presence of the membrane mimetic are not only due to the huge difficulties in the sample preparation (every batch could be different) but also that the recorded effects based on the technique suggested by the reviewer could be mediated due to complicated interactions of the ligand with the membrane mimetic, or the membrane proteins or both. Furthermore, STD or saturation transfer double difference (STDD) NMR, works best for weakly bound ligands, but not for tightly bound ligands. Although membrane proteins constitute a very important field of drug discovery, since this class of proteins bears more than 50% of current drug targets, current challenges in determining the ligand-binding interaction profile exist. Herein, we suggest a workflow, following appropriate validations, that could assist on the determination of the ligand binding interaction poses of membrane proteins.