Molecular Modeling of Adenosine Receptors

CHAPTER THREE

Molecular Modeling of Adenosine Receptors Adriano Martinelli1, Gabriella Ortore Dipartimento di Farmacia, University of Pisa, Pisa, Italy 1 Corresponding author: e-mail address: [email protected]

Contents 1. Introduction 2. Adenosine A2a Receptor 2.1 X-ray structures 2.2 Receptor activation 3. Receptor Modeling 3.1 Building adenosine receptor models 3.2 Homology modeling 3.3 MM/PBSA 3.4 Validation methods 3.5 Receptor oligomerization 4. Summary References

38 38 38 40 41 41 42 53 55 56 57 57

Abstract Adenosine is a ubiquitous neuromodulator for which four subtypes are known (A1, A2a, A2b, and A3). Adenosine receptors (ARs) are members of the superfamily of G proteincoupled receptors. Knowledge of the 3D structure of ARs is a fundamental tool for studying how they work and for designing newer, more potent, and selective ligands. Here, we describe the procedures for modeling ARs. The starting point is represented by the resolved X-ray crystallographic structures of the A2a adenosine receptor (AA2AR) while, in the case of the other subtypes (A1, A2b, and A3), homology modeling procedures using AA2AR as a template are necessary. Minimizations and MD simulations are the computational methods used for structural refinement. The ligand–receptor interaction is studied through molecular docking simulations. The necessity of validating the obtained computational models through all of the available experimental data is outlined. Finally, the emerging importance of the homo- and heterodimerization of ARs is taken into account and a procedure for modeling a dimer is reported.

Methods in Enzymology, Volume 522 ISSN 0076-6879 http://dx.doi.org/10.1016/B978-0-12-407865-9.00003-0

#

2013 Elsevier Inc. All rights reserved.

37

38

Adriano Martinelli and Gabriella Ortore

1. INTRODUCTION Four adenosine receptor (AR) subtypes belonging to the rhodopsinlike family of G protein-coupled receptors (GPCRs) have been cloned and characterized in several species including rat, mouse, and human. They are encoded by distinct genes. The A1 (AA1R) and A3 (AA3R) subtypes mainly signal via Gi proteins to mediate the inhibition of adenylyl cyclase and modulate the activity of several Kb and Ca2b channels, whereas the A2a (AA2AR) and A2b (AA2BR) subtypes mainly signal via the Gs proteins, causing the activation of adenylyl cyclase and thus stimulating the formation of cAMP. In recent years, several experiments have provided insights into the physiology and pathophysiology of these four subtypes (Fredholm, Ijzerman, Jacobson, Linden, & Muller, 2011). ARs, like all other GPCRs, are made up of seven transmembrane a-helices (7TM) of approximately 25 residues in length that are connected by intra- and extracellular loops. The helices are placed in a lipidic environment, while the loop regions are surrounded by an aqueous medium. The 3D structure of a protein is usually obtained mainly from X-ray crystallography studies, but for most GPCRs, crystallization is still quite difficult. After the first high-resolution structure of bovine rhodopsin (bRh), published 12 years ago, human AA2AR, turkey b1-adrenoceptor, human b2-adrenoceptor, human D3 dopamine, chemokine receptor CXCR4, and human histamine H1 receptor were also crystallized, providing better templates for the theoretical construction of human GPCRs. The first requirement for an accurate and consistent homology modeling (HM) result is a high homology between the target protein and the template, which is decidedly increased by using human GPCRs as templates rather than bRh. In particular, for ARs, AA2AR is fortunately one of the most experimentally studied GPCRs, and the percentage of homology grows from an average of less than 20%, considering the whole sequence of adenosine subtypes against bRh, to 40% using AA2AR as a template.

2. ADENOSINE A2a RECEPTOR 2.1. X-ray structures In 2008, the first AA2AR was engineered with the introduction of a T4L segment in the IL3 region and crystallized in a complex with the antagonist ZM241385 (Jaakola et al., 2008). The ZM241385 ligand lies perpendicular

Molecular Modeling of Adenosine Receptors

39

to the membrane plane, in a unique ligand-binding pocket, interacting with Glu169 (EL2) and Asn253 (6.55). Furthermore, a p-stacking with Phe168 (EL2) provides an aromatic key interaction. Some cocrystallized water molecules give an H-bond network that connects the ligand scaffold and the binding-site residues. In 2011, the antagonists ZM241385, XAC, and caffeine were successfully crystallized into the receptor pocket of a stabilizing construct (A2A-StaR2) by Dore et al. (2011). The receptor was crystallized in the inactive state conformation in the presence of the ionic lock between the side chains of Glu228 (6.30) and Arg102 (3.50). The ligand-binding cavities of ZM241385, XAC, and caffeine deviate from each other to a detectable extent, suggesting a degree of ligand–receptor-induced fit. The antagonist ZM241385 was recently cocrystallized with a Fab fragment by Hino et al. (2012). The antibody fragment Fab2838 induces an inactive conformation, provoking significant displacements of TM3, TM6, and TM7, without blocking the ligand-binding site. The crystal structures of thermostabilized human AA2AR (A2ARGL31) bound to its endogenous agonist adenosine and the synthetic agonist NECA were reported by Lebon et al. (2011). In the authors’ opinion, this structure represented an intermediate conformation between the active and inactive states because of both structural and pharmacological data. The structure of A2a-T4L complexed with the agonist UK432097 (Xu et al., 2011) is very similar to the structures of A2AR-GL31-NECA in the transmembrane regions, but there are differences in the extracellular surface due to the bulky extensions of UK432097 interacting with ECL3. In Fig. 3.1, the structure of AA2AR complexed with adenosine is shown. In summary, at present, nine structures are present in the Protein Data Bank; their codes are 2YDO, 2YDV (Lebon et al., 2011), 3EML ( Jaakola et al., 2008), 3PWH, 3RFM, 3REY (Dore et al., 2011), 3QAK (Xu et al., 2011), 3VGA, and 3VG9 (Hino et al., 2012). All of these structures are complexes between the receptor and an agonist (2YDO, 2YDV, 3RFM, 3REY, 3QAK), an antagonist (3EML, 3PWH), or an allosteric inverse agonist (3VG9, 3VGA). In order to crystallize these complexes, some stabilizations through mutations were necessary; in the case of 3EML and 3QAK, the stabilization was obtained through the substitution of the intracellular loop 3 with T4/lysozyme, for 3PWH, 3REY, 3RFM, 2YDO, and 2YDV, through the introduction of a small number of point mutations which led to a conformational thermostabilization.

40

Adriano Martinelli and Gabriella Ortore

Figure 3.1 The X-ray structure of AA2AR complexed with adenosine (PDB ID: 2YDO).

To build AA2AR models, the structural modification of the experimental structure with respect to the wild one and the type of form (active or inactive) of the experimental structures needs to be taken into account. GPCRs exist in various conformationally different forms. At present, there are no complete and reliable studies of the dynamics of GPCRs, and generally, only two conformations are taken into consideration: the active one and the inactive one.

2.2. Receptor activation Agonist binding shifts the equilibrium toward the active conformation, and the transmembrane helices are rearranged on the intracellular side to bind the G protein. Agonists can bind to the binding site with a high affinity only if it is in the active conformation. An antagonist can bind to both conformations of the binding site, but it is unable to shift the equilibrium. Finally, an inverse agonist locks the receptor in the inactive conformation (Hino et al., 2012). Analyzing the structures of the studied complexes in both active and inactive form, a reasonable hypothesis for AR activation appears quite clear: a small conformational change within the ligand-binding cavity induces

41

Molecular Modeling of Adenosine Receptors

Glu(6.30) Arg(3.50)

TM7

TM3

TM6

Figure 3.2 On the left: superimposition of helices TM3, TM6, and TM7 in the inactive (lighter colors, 3VG9) and active (darker colors, 3QAK) states of AA2AR. On the right: ionic lock between Arg102(3.50) and Glu228(6.30) unlocked in the active state and locked in the inactive state.

helical rearrangement over the 7TM domains that spread toward the intracellular side. These changes include an outward tilt of TM6, together with a movement of TM5, and an inward tilt of TM7. As a result, the intracellular interface is reshaped to promote the binding of G proteins. The most evident differences between the active and the inactive conformations are the position and conformation of TM3, TM6, and TM7, and the presence in the inactive form of the so-called ionic lock, an ionic bond between two residues highly conserved in GPCRs, Arg102 (3.50) and Glu228 (6.30), in the AA2AR. In Fig. 3.2, the conformation of TM3, TM6, and TM7 in the active and inactive states are superimposed and the ionic lock in the two states is shown.

3. RECEPTOR MODELING 3.1. Building adenosine receptor models The building of a protein model is a common task for which a great number of computer procedures have been developed (Martinelli & Tuccinardi, 2006). The most favorable case is that in which an experimental structure of the protein to be modeled is available. As regard the ARs, fortunately, A2AR is one of the most studied GPCRs experimentally. However, also in the case of A2AR, the X-ray crystal structure cannot be considered the final one

42

Adriano Martinelli and Gabriella Ortore

and some optimization may be necessary. Particularly, the kind of form (active or inactive) of the experimental structures needs to be taken into account. The procedure for optimizing an X-ray structure is simple and generally consists of a simple minimization and the addition of hydrogen atoms. However, if the X-ray structure needs some modification, and especially if the cocrystallized ligand has to be replaced by a different one, the procedure is more complicated and needs to be more precise. In these cases, an optimization procedure analogous to that described in Section 3.2.4 is used.

3.2. Homology modeling When there is no experimental structure for a GPCR, it can be built through a HM procedure, in which the target protein is built, starting from the experimentally known 3D structure of related proteins. In such a procedure, the sequence alignment of homologue proteins can be converted into 3D models using the assumption that, during their evolution, the 3D structure of the proteins has been better conserved than their individual sequences (Martinelli & Tuccinardi, 2006). HM is a common computational model-building procedure mainly used for GPCRs and schematically consists of three steps: (1) the alignment of the receptor sequence on a suitable template, (2) the modeling of the receptor (the core of the TM domains, loops, and side chains), (3) the optimization of the structure through molecular dynamic (MD) simulations and molecular mechanic (MM) minimizations. The critical factor in this procedure is the degree of homology of the template with respect to the protein to be modeled (Martinelli & Tuccinardi, 2008). 3.2.1 Alignment This first step consists of the alignment between the primary sequence of the HM receptor and that of a template. In the case of AA1R, AA2BR, and AA3R, the most suitable template is surely one of the released structures of AA2AR. For this task, the set of conserved residues on the 7TM domains is fundamental. The alignment of the sequences of the four ARs is shown in Fig. 3.2. The best results have generally been obtained by multialigning all the receptors of the same subclass with the template. In this way, the evolutionary relationships of the subclass are taken into account and alignment accuracy improves.

Molecular Modeling of Adenosine Receptors

43

Many programs are available for alignment, but the most accurate is the CLUSTAL W software (Thompson, Higgins, & Gibson, 1994). CLUSTAL W assigns individual weights to each sequence in partial alignment for downweighting near-duplicate sequences and upweighting the most variable ones. The gaps are treated in such a manner as to favor gaps in potential loop regions and penalize those in potential helix regions. 3.2.2 Modeling Based on the alignment step, the TM domains of the receptor model are generally constructed starting from the X-ray of the selected template. The loop fragments are more variable in conservation and especially in length and thus can require special modeling treatment. In the choice of the X-ray structure to be used as a template, some aspects are to be taken into account: • the kind of state, active or inactive, of the template chosen, on the basis of the pharmacological properties of the ligand to be inserted into the model; • the similarity of shape and chemical features between the ligand to be inserted into the model and the ligand cocrystallized in the template; • the resolution of the template and the presence in it of mutations and engineered portions. The MODELLER program is the most used (Sali & Blundell, 1993) for performing the modeling step; it models the 3D structure of proteins by using spatial restraints (distances and torsion angles) derived from the template structure. The final 3D model is then obtained by optimizing a molecular probability density function. The following step is the adjustment of the side chain conformations. The side chain conformation depends on the residue type and on the conformational state of the closer residues. This step in the modeling procedure is usually performed with simulated annealing methods, which sample all of the possible side chain conformations and evaluate their relative stability. It is always advisable, and in some cases necessary, to build a certain number of model conformations and, after the optimization, to evaluate them and select the most reliable. 3.2.3 Detailed procedure of the alignment and modeling steps The selected crystal structure of human AA2AR bound to a ligand is taken from the Protein Data Bank, while all of the primary sequences of the other ARs are obtained from the SWISS-PROT protein sequence database.

44

Adriano Martinelli and Gabriella Ortore

Figure 3.3 Alignment of the human adenosine receptors. The highly conserved patterns are marked with black. The other identical residues are in bold. The first and second lines of the alignment scheme are reported the numerations of the human A1AR, while the TM domains of AA2AR are reported in gray.

The sequential alignment is performed by means of CLUSTAL W using the Blosum series as a matrix, with a gap open penalty of 10 and a gap extension penalty of 0.05. In Fig. 3.3, the so-obtained alignment of the human ARs is reported. As an example, let us consider the modeling of human AA1R using 3EML as templates (Tuccinardi et al., 2011). The alignment was guided by the highly conserved amino acid residues, including the asparagine residue N1.50, the LA-AD (L2.46, A2.47, A2.49, and D2.50), and the D/ERY-V (D/E3.49, R3.50, Y3.51, and V3.54) motif, the highly conserved tryptophan W4.50, the two prolines P4.59 and P6.50, and the NPXXY motif in TM7 (N7.49, P7.50, and Y7.53). About 45% of identical residues were found, with an alignment score of 48. The human A1R was constructed directly from the coordinates of the corresponding amino acids of the template through the MODELLER program. Then 10 structures were generated using the “very slow MD annealing” refinement method, and the best receptor model was chosen following the discrete optimized protein energy assessment method. The backbone conformation of the model was validated by inspection of the Psi/Phi Ramachandran plot obtained from PROCHECK analysis (Laskowski, Macarthur, Moss, & Thornton, 1993). No or very few residues should be found in a disallowed geometry.

Molecular Modeling of Adenosine Receptors

45

3.2.4 Optimization A combination of MD and minimization methods should be used with the aim of finding and correcting possible alignment errors and making the model more stable. However, this step should be carried forward cautiously, as MD refinements may distort the homology models rather than improving them (Martinelli & Tuccinardi, 2006). The aim of a dynamic simulation is really to stabilize the structure without inducing important modification, which in most cases is comprised of artifacts created by the calculation. A less used alternative to MD simulation is the use of a simulated annealing MD simulation applied to the receptor–ligand complexes. Briefly, this method consists of a series of cooling cycles, starting from a high temperature (400–1000 K), and then cooling the system to a lower temperature (100–200 K). In this way, several low-energy conformations are obtained, and the selection of the best one is made, taking into account the position of the residues suggested by mutagenesis data and by structure–activity relationship (SAR) studies as (a) critical for the ligand affinity and (b) probably shaping the binding site. Mutagenesis data are the second-most useful data, just after the crystallographic ones, for modeling receptors. These data suggest which are the most important residues and what their role might be in determining the biological behavior of the receptor. A list of the most important mutation data for the ARs can be found in Table 3.1. The easiest way to optimize a receptor model requires vacuum MD simulations followed by minimization protocols, thereby skipping the explicit environment, but in this case, a set of restraints is required to replace the natural stabilizing effects of the membrane bilayer on the TM domains. At present, the computational power normally available allows the use of an explicit model of membrane and solvent. These MD simulations are run with different types of lipid molecules, mainly 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphatidylcholine, dimyristoylphosphatidylcholine, dipalmitoyl phosphatidylcholine, and palmitoyl oleylphosphatidyl ethanolamine. The receptor modeling has as its final aim: the study of the interaction of ligands, agonist or antagonist, with the receptor; therefore, the optimization protocol is normally applied to a receptor model complexed with a highaffinity ligand. In this way, the modifications of the binding-site geometry induced by the interaction with the ligand can be taken into account. The ligand insertion inside the receptor model can be obtained through either a manual or automated docking procedure. In both cases, the ligand

46

Adriano Martinelli and Gabriella Ortore

Table 3.1 AR residues analyzed with mutagenesis studya Residue/mutational results

1.36

AA2BR

V11I: no variation

1.37

AA1R

G14T: increased agonists affinity

AA2BR

A12T: no variation

AA1R

E16A: reduced agonists affinity E16Q: reduced antagonists affinity

AA2AR

E13Q: reduced agonists affinity E13A: reduced CGS21680 affinity, increased LUF5834 affinityb

1.48

AA1R

P25L: reduced agonists affinity

1.50

AA3R

N30A: decreased affinity

1.54

AA1R

I31C: no variation

IL1

AA2BR

N36D: no variation

AA2BR

T42A: decreased agonists affinity

2.41

AA1R

C46A/S: no variation

2.45

AA1R

S50A: no variation

2.50

AA1R

D55A: increased agonists affinity

AA3R

D58N: no variation

2.51

AA2BR

V54L: decreased agonists affinity

2.55

AA2BR

L58V: no variation

2.56

AA2BR

F59L: loss of affinity

2.60

AA1R

L65F: no variation L65T: decreased agonists affinity

2.64

AA1R

I69S/T/V: decreased agonists affinity

EL1

AA2BR

F71V/H/D/C/S: decreased agonists affinityc F71I/E/Y/L/M no variationc F71L/D74G: increased agonists affinityc

3.25

AA1R AA2BR

C80A/S: no detectable binding C78S: loss of activityd

3.27

AA1R

M82F: no variation

3.30

AA1R

C85A: no variation C85S: reduced agonists affinity

AA3R

C88F: decreased agonists affinity

1.39

Molecular Modeling of Adenosine Receptors

47

Table 3.1 AR residues analyzed with mutagenesis study—cont'd Residue/mutational results

3.31

AA1R

P86F: reduced agonists affinity

AA2BR

F84L/S: decreased agonists affinity

AA1R

V87A: reduced antagonists affinity

AA2AR

V84L: marginal variation

AA1R

L88A: reduced affinity

AA2AR

L85A: decreased agonists affinityb L85R: decreased agonists affinityb

AA1R

T91A: reduced affinity

AA2AR

T88A/R/S: reduced affinityb T88D/E: loss of agonists binding

AA3R

T94A: decreased affinity

AA1R

Q92A: reduced affinity

AA2AR

Q89A/D: increased affinity

AA3R

H95A: decreased affinity

AA1R

S93A: no variation

AA2AR

S90A: marginal variation

AA2BR

S91G: no variation

AA1R

S94A: no detectable binding S94T: reduced antagonists affinity

AA2AR

S91A: marginal variation

3.49

AA3R

D107K/N/R: marginal variation

3.50

AA3R

R108E/H/K/N: marginal variation

3.51

AA3R

Y109F: decreased agonists affinity

4.43

AA1R

A125K: no variation

4.49

AA1R

C131A/S: no variation

4.53

AA1R

S135A: no variation

4.56

AA2BR

I136T: increased agonists affinitye

4.59

AA1R

T141A: no variation

4.61

AA2BR

F141L: increased agonists affinitye

3.32

3.33

3.36

3.37

3.38

3.39

Continued

48

Adriano Martinelli and Gabriella Ortore

Table 3.1 AR residues analyzed with mutagenesis study—cont'd Residue/mutational results

4.62

AA1R

F144L: no variation

EL2

AA3R

K152A: decreased antagonists affinity

AA2AR

E151A/D/Q: loss of affinity

AA2BR

C154S: decreased antagonist affinity, no variation on agonists affinityd

AA2AR

E161A: increased antagonists affinity

AA2BR

C166L: increased agonists affinitye C166S: no variationd

AA2BR

C167S: increased agonists affinitye; C167S no variationd

AA2BR

L168P: increased agonists affinitye

AA1R

C169A: no detectable binding

AA2BR

C171S: loss of activityd

AA3R

Q167A/E/R: decreased affinity

AA2AR

F168A: no detectable bindingf; decreased agonists affinityb F168Y: decreased antagonist and agonists affinityb,f F168W: decreased antagonist and marginal variation on agonists affinityb,f

AA2BR

F173: increased agonists affinitye

AA2AR

E169A: loss of antagonists affinity, no variation in agonists affinityb E169Q: no variation

AA2AR

D170K: no variation

AA2AR

P173R: no variation

5.38

AA2AR

Met177A: decreased antagonist affinity; decreased CGS21680 affinity; no variation on NECAf

5.41

AA2AR

F180A: marginal variation

5.42

AA2AR

N181S: reduced agonists affinity

Molecular Modeling of Adenosine Receptors

49

Table 3.1 AR residues analyzed with mutagenesis study—cont'd Residue/mutational results

5.43

AA2AR

F182A: loss of affinity F182Y/W: reduced agonists affinity

AA3R

F182A: decreased antagonists affinity

5.56

AA2BR

V200G: increased agonists affinitye

5.58

AA2BR

Y202S: increased agonists affinitye

6.34

AA3R

A229E: no variation

6.48

AA3R

W243A/F: decreased antagonists affinity

AA2AR

W246A: decreased CGS21680 affinity, increased LUF5834 affinityb W246Y: increased agonists affinityb

6.49

AA3R

L244A: no variation

6.51

AA2AR

L249A no detectable bindingf

6.52

AA1R

H251L: no variation

AA2AR

H250A: loss of affinity H250F/Y: decreased agonists affinity H250N: increased agonists affinity

AA3R

S247A: no variation

AA2AR

N253A: loss of affinityb

AA3R

N250A: loss of affinity

AA1R

C255A/S: no variation

AA2AR

C254A: marginal variation

6.59

AA2AR

F257A: loss of affinity

EL3

AA1R

C260A/S: no variation

AA1R

C263A/S: no variation

AA2AR

C262G: no variation

AA2AR

H264A: decreased CGS21680 affinity; increased LUF5834 affinityb

7.36

AA2BR

N273Y: no variation

7.39

AA2AR

I274A: loss of affinity

6.55

6.56

Continued

50

Adriano Martinelli and Gabriella Ortore

Table 3.1 AR residues analyzed with mutagenesis study—cont'd Residue/mutational results

7.42

AA1R

T277A/S: reduced affinity

AA2AR

S277A: reduced agonists affinity; increased LUF5834 affinityb S277C: reduced agonists affinity S277N/E/T: no variation

AA1R

H278L: loss of affinity

AA2AR

H278A: loss of affinity H278D/E: no variation H278Y: decreased affinity

AA3R

H272E: decreased affinity

7.45

AA3R

N274A: decreased affinity

7.46

AA2AR

S281A: loss of affinity S281T: increased affinity S281N: increased agonists affinity; decreased antagonists affinity

7.53

AA3R

Y282F: decrease of agonists affinity

7.43

a

Dal Ben, Lambertucci, Marucci, Volpini, and Cristalli (2010). Lane et al. (2012). c Peeters et al. (2011). d Schiedel et al. (2011). e Peeters, Li, van Westen, and Ijzerman (2012). f Jaakola et al. (2010). b

disposition should be in agreement with ligand-based information and mutagenesis data, following the golden rule that molecular modeling needs to be driven and validated by as much experimental data as possible (Martinelli & Tuccinardi, 2006). In the early stages of receptor–ligand MD simulations, the main ligand–receptor interactions can be stabilized through suitable constraints, thereby allowing some adjustment of the receptor-binding site. However, the main part of the simulation is to be performed without any unnatural restraint. As regards the length of the MD simulations, some experiments suggest that too extensive refinement generally tends to bring a homology model farther from rather than closer to the true structure (Martinelli & Tuccinardi, 2006). On average, GPCR modeling should use few nanoseconds of MD simulation length and, in general, the more complex the model (including

Molecular Modeling of Adenosine Receptors

51

explicit environment), the longer the simulation time. In any case, a simulation cannot be stopped until a certain structural stability is reached, and this is obtained when the root mean square deviation (RMSD) of heavy atoms remains quite constant. If this stability is not observed after some (2–4) nanoseconds, it could mean that the starting structure is probably not reliable. A validation of the final model is always necessary to verify that the ligand possesses some of the interactions suggested by the mutagenesis data and/or the experimental SAR. As already pointed out, more conformations are generated in the modeling step and all them are optimized, but then the most reliable one has to be selected. This task can be performed through one of these approaches or a combination of them: 1. a visual evaluation of the ligand–receptor interactions, on the basis of the mutagenesis data, of the SAR analysis of the ligand and of the interaction shown by the cocrystallized ligands in the templates; 2. a rescoring of the complexes through one or more docking programs and scoring functions (Tuccinardi, Ortore, Manera, Saccomanni, & Martinelli, 2006), or through building a receptor-based 3D-QSAR model using the docking poses as an alignment tool (Ortore, Tuccinardi, Bertini, & Martinelli, 2006); 3. a calculation of the value of the interaction energy through MM calculations (Tuccinardi et al., 2008) or, more rigorously, through the Molecular Mechanic/Poisson-Boltzmann Surface Area (MM/PBSA) method (Tuccinardi et al., 2011). 3.2.5 Detailed procedure of the optimization step As an example, a procedure recently used for modeling a complex between AA1R and an antagonist is described below (Tuccinardi et al., 2011). A compound (it can be named, e.g., as C) possessing a high potency for AA1R was selected (actually, C has to be selected before starting the alignment and the modeling steps). C is an antagonist; therefore, a template in the inactive form was chosen (in this case 3EML), whereas if it had been an agonist, a different template would have been chosen. Moreover, if more templates with the right form are available, it may be better to select one that is complexed with a ligand more sterically similar to C. C was built using Maestro and subjected to a conformational search (CS) of 1000 steps, using a water environment model (generalized-Born/surfacearea model) by means of Macromodel (Monte Carlo method, MMFFs force field, and a distance-dependent dielectric constant of 1.0). C was then

52

Adriano Martinelli and Gabriella Ortore

minimized using the conjugated gradient method until a convergence value of 0.05 kcal/(A˚ mol) (same forcefield and parameters used for the CS) was obtained. Automated docking was performed through the GOLD program (Verdonk, Cole, Hartshorn, Murray, & Taylor, 2003). The binding site for C, defined by taking into account the interaction of the cocrystallized ligand in the X-ray structure of AA2AR, was individuated in the region between TM2, TM3, and TM5–TM7. The region of interest identified by GOLD was defined as that containing the residues within 10 A˚ from the original position of the cocrystallized ligand in the AA2AR X-ray structure. The “allow early termination” option was deactivated, while the possibility of clustering the results was activated (the RMSD distance clus˚ ). All other parameters were left at their default value. tering was set to 1.5 A C was submitted to 30 genetic algorithm runs by applying the GoldScore scoring function. The clusters showing a scoring value less than 10 units lower than the best docked conformation were considered; in this way, four docking conformations were found. All four A1–ligand complexes were used in turn as starting points for 5 ns of MD simulations. A fully hydrated phospholipid bilayer environment made up of palmitoyloleoylphosphatidylcholine (POPC) was used for the MD simulations. An explicit solvent model TIP3P water was used, and the system was solvated on the “extracellular” and “intracellular” side with ˚ water cap. a 12 A All simulations were performed using AMBER 10 and the modified parm94 as a forcefield. The creation of the phospholipid bilayer and the insertion of the receptor–ligand complexes were carried out through the facility of the VMD program (Humphrey, Dalke, & Schulten, 1996). Chlorine ions were added as counterions to neutralize the system. The general Amber force field parameters were assigned to POPC molecules. The partial charges were calculated using the AM1-BCC method, as implemented in the Antechamber suite of AMBER 10. The system therefore contained 230 POPC molecules, 23,415 water molecules, 12 chlorine atoms, the AA1R, and compound C, for a total of 106,121 atoms. Three minimization steps were carried out before MD: 1. the protein and phospholipids were fixed with a position restraint of 100 kcal/mol/A˚2, and only the positions of the water molecules were minimized; 2. the phospholipid–water system was minimized, while a position restraint of 100 kcal/mol/A˚2 was applied to the protein; 3. a restraint of 30 kcal/mol/A˚2 was applied only on the a carbons of the receptor.

Molecular Modeling of Adenosine Receptors

53

The three minimization stages consisted of 10,000 steps; the first 1000 steps were Steepest Descent, and the last 9000 were conjugate gradient. MD trajectories were run, starting from the minimized structure. Particle mesh Ewald electrostatics and periodic boundary conditions were used in the simulation. The time step of the simulations was 2.0 fs, with a cutoff ˚ for the nonbonded interaction. The SHAKE option was employed of 12 A to keep all bonds involving hydrogen atoms rigid. A first constant-volume MD equilibration was carried out for 200 ps in order to raise the temperature from 0 to 300 K (the Langevin dynamics method was used). Then, 4800 ps of constant-pressure MD were carried out at 300 K. During the first 400 ps, all of the a carbons of the receptor remain blocked with a harmonic force constant, which decreased during these 400 ps from 10 to ˚ . In the last 4400 ps, there were no constraints. 1 kcal/mol
A Each of the four final structures was obtained as the average of the last 4400 ps of MD minimized using the CG method until a convergence of 0.05 kcal/A˚ mol was obtained. The backbone conformation of the resulting receptor structure was evaluated by inspection of the Ramachandran plot obtained from PROCHECK analysis. In Fig. 3.4, the monitoring of a MD runs through the analysis of the trend of the total energy and the displacement of the structure from the starting one during 5 ns of simulation.

3.3. MM/PBSA The MM/PBSA approach is a method to evaluate the free energies of binding. This approach averages the contributions of gas-phase energy, solvation free energy, and solute entropy from snapshots (calculated through MD simulations) of the complex molecule, as well as the unbound components. A snapshot of every 50 ps is extracted from the last 4400 ps of MD for a total of 88 snapshots (at time intervals of 50 ps) for each species (complex, receptor, and ligand). The total MM-PBSA free energy of binding is then computed as: G ¼ Gpolar þ Gnonpolar þ EMM  TS The EMM term, which includes electrostatic, van der Waals, and internal energies, is calculated using AMBER 10. The polar energy term (Gpolar) is calculated from the PBSA module of the AMBER 10 program through the Poisson–Boltzman method using dielectric constants of 1 and 80 to

54

Adriano Martinelli and Gabriella Ortore

Total energy (kcal/mol)

A

-150,000

-160,000

-170,000

-180,000

0

1000

2000

3000

4000

5000

4000

5000

Time (ps) B

3.5 3.0

RMSD (Å)

2.5 2.0 1.5 1.0 0.5 0.0 0

1000

2000

3000

Time (ps)

Figure 3.4 Analysis of the MD simulation of the ligand C complexed with AA1R. In the upper plot, the total energy of the system versus the time; in the lower plot, the RMSD of the a carbons of the TM helices and the heavy atoms of the receptor with respect to the starting structure versus the time.

represent the gas and water phases, respectively. The nonpolar energy term (Gnonpolar) is calculated through the MOLSURF program. In order to compare the binding free energy of the various conformers with the receptor, only the first three terms are taken into account, while the entropic term TS is considered approximately constant.

Molecular Modeling of Adenosine Receptors

55

3.4. Validation methods A general rule of molecular modeling states that each new model should be “validated” in order to evaluate its reliability (Martinelli & Tuccinardi, 2006). It is necessary to verify whether or not the model is able to explain the biopharmacological properties of known ligands different from those used for building the model or to predict new active ligands. The computational procedure is to dock a series of ligands into the model (Tuccinardi et al., 2006). This docking might be accomplished through automated procedures using suitable docking software, or through a manual disposition of the ligands inside the receptor, eventually followed by some optimization. The first approach can be used when the structure, and especially the steric hindrance of the considered ligands, is not too different from that of the ligand used for modeling the receptor. On the contrary, the other approach is used when the ligands to be docked are structurally different and therefore are unable to interact with the receptor without some adjustments of its binding site. In this case, a remodeling of the receptor through MD simulation may be necessary (Cosimelli et al., 2008). The procedure to be used in this case is similar to that already described. In addition, often the length of such MD simulation can be shorter (1–2 ns). In any case, the ligand–receptor interactions shown by the docking analysis need to be able to explain ligand activities. It is important that potent ligands show interactions involving residues suggested to be important by mutagenesis studies and that these interactions remain stable during an MD simulation. From a quantitative point of view, the ligands docked into a receptor model can also be evaluated through a scoring function so that a quantitative correlation between the experimental and calculated activity can be used as a substantial validation term. However, the nature of all scoring functions is strongly empirical, thus determining a low reliability of the results so obtained. An improvement can be sought by combining more docking programs and scoring functions. Another possibility for improving the quantitative evaluation of ligand–receptor interactions is to make use of the ligand-docking process for a receptor-based 3D-QSAR analysis as an alternative scoring method able to predict the activity of ligands. The successful design of new ligands is to be considered a strong validation method. Virtual screening is a technique that is used for filtering a molecular database and extracting compounds that promise to have a certain

56

Adriano Martinelli and Gabriella Ortore

biological activity and, therefore, become candidates for experimental biological tests. Structure-based virtual screenings use a docking analysis to filter the database, but due to the great number of compounds in the database, it is advisable to apply in advance some other filters able to discharge compounds not possessing drug-like characteristics.

3.5. Receptor oligomerization 3.5.1 Adenosine receptor homo- and heterooligomerization Some experimental evidence indicates that many GPCRs can form homoand/or heteroassemblies of two or more monomers (Guidolin et al., 2011). Among them, AA2AR shows a similar behavior. In particular, the existence of A2AR homodimers was demonstrated by time-resolved fluorescence resonance energy transfer. At present, it seems that the homodimer may be a functional form of this receptor. 3.5.2 Procedure for modeling an AA2AR dimer The modeling starts with rigid-body docking simulations through, for example, ZDOCK software (Chen & Weng, 2003) using two identical copies of the AA2AR monomer (Fanelli & Felline, 2011). For this purpose, one of the available X-ray structures is used. One of the monomers is maintained fixed, while the other one is rotated in 6 steps and moved around in order to find favorable interactions. In this way, several thousand possible dimer conformations are generated and ranked according to a scoring function (e.g., a ZDOCK score). All of the dimer arrangements characterized by a too-great tilt angle (>23 ) or too˚ ) are discarded, and a clusterization of the great vertical displacement (>6.0 A remaining ones leads to the prediction of some dimers that are manually adjusted to relieve bad contacts. The predicted dimers are then subjected to energy minimizations and MD simulations following a protocol similar to that described above for modeling a single monomer. In these calculations, a cocrystallized ligand is retained at the binding site of both monomers. If, due to the large number of atoms present in the dimer, it is impossible to use an explicit membrane and solvent, the structure of the two monomers has to be stabilized by applying suitable distance restraints between the residues of the binding site and the ligand, as well as between the helices.

Molecular Modeling of Adenosine Receptors

57

4. SUMMARY In the research fields targeting drug discovery or understanding biological mechanisms, molecular modeling of protein targets is a necessary tool. A protein model can be considered reliable only if it is supported by experimental data, and the most favorable case is that in which experimental X-ray or NMR structures are available. Mutagenesis data are also very important for evidencing the role of residues fundamental for interaction with ligands. At present, high-resolution experimental structures of AA2AR in both active and inactive states are available and allow the building of reliable molecular models. Moreover, the present available computational power allows MD simulation on complex systems with explicit membranes and solvents, thereby avoiding the introduction of unnatural restraints. However, we need to remember the empirical nature of molecular modeling. The possibility of using templates structurally very close to the target protein is not sufficient to assure the reliability of the model. Additionally, a model of AA2AR in a complex with a ligand different with respect to the cocrystallized ones needs to be validated using all of the available experimental data.

REFERENCES Chen, R., & Weng, Z. (2003). A novel shape complementarity scoring function for protein–protein docking. Proteins, 51, 397–408. Cosimelli, B., Greco, G., Ehlardo, M., Novellino, E., Da Settimo, F., Taliani, S., et al. (2008). Derivatives of 4-amino-6-hydroxy-2-mercaptopyrimidine as novel, potent, and selective A3 adenosine receptor antagonists. Journal of Medicinal Chemistry, 51, 1764–1770. Dal Ben, D., Lambertucci, C., Marucci, G., Volpini, R., & Cristalli, G. (2010). Adenosine receptor modeling: What does the A2A crystal structure tell us? Current Topics in Medicinal Chemistry, 10, 993–1018. Dore, A. S., Robertson, N., Errey, J. C., Ng, I., Hollenstein, K., Tehan, B., et al. (2011). Structure of the adenosine A2A receptor in complex with ZM241385 and the xanthenes XAC and caffeine. Structure, 19, 1283–1293. Fanelli, F., & Felline, A. (2011). Dimerization and ligand binding affect the structure network of A2A adenosine receptor. Biochimica et Biophysica Acta, 1808, 1256–1266. Fredholm, B. B., Ijzerman, A. P., Jacobson, K. A., Linden, J., & Muller, C. E. (2011). International Union of Basic and Clinical Pharmacology I. XXXI. Nomenclature and classification of adenosine receptors—An update. Pharmacological Reviews, 63, 1–34. Guidolin, D., Ciruela, F., Genedani, S., Guescini, M., Tortorella, C., Albertin, G., et al. (2011). Bioinformatics and mathematical modelling in the study of receptor–receptor interactions and receptor oligomerization. Focus on adenosine receptors. Biochimica et Biophysica Acta, 1808, 1267–1283.

58

Adriano Martinelli and Gabriella Ortore

Hino, T., Arakawa, T., Iwanari, H., Yurugi-Kobayashi, T., Ikeda-Suno, C., NakadaNakura, Y., et al. (2012). G-protein-coupled receptor inactivation by an allosteric inverse-agonist antibody. Nature, 482, 237–240. Humphrey, W., Dalke, A., & Schulten, K. (1996). VMD—Visual molecular dynamics. Journal of Molecular Graphics, 14, 33–38. Jaakola, V. P., Griffith, M. T., Hanson, M. A., Cherezov, V., Chien, E. Y., Lane, J. R., et al. (2008). The 2.6 angstrom crystal structure of a human A2A adenosine receptor bound to an antagonist. Science, 322, 1211–1217. Jaakola, V. P., Lane, J. R., Lin, J. Y., Katritch, V., IJzerman, A. P., & Stevens, R. C. (2010). Ligand binding and subtype selectivity of the human A2A adenosine receptor: Identification and characterization of essential amino acid residues. The Journal of Biological Chemistry, 285, 13032–13044. Lane, J. R., Herenbrink, C. K., van Westen, G. J. P., Spoorendonk, J. A., Hoffmann, C., & IJzerman, A. P. (2012). A novel nonribose agonist, LUF5834, engages residues that are distinct from those of adenosine-like ligands to activate the adenosine A2A receptor. Molecular Pharmacology, 81, 475–487. Laskowski, R. A., Macarthur, M. W., Moss, D. S., & Thornton, J. M. (1993). Procheck— A program to check the stereochemical quality of protein structures. Journal of Applied Crystallography, 26, 283–291. Lebon, G., Warne, T., Edwards, P. C., Bennett, K., Langmead, C. J., Leslie, A. G. W., et al. (2011). Agonist-bound adenosine A(2A) receptor structures reveal common features of GPCR activation. Nature, 464, 521–525. Martinelli, A., & Tuccinardi, T. (2006). An overview of recent developments in GPCR modelling: Methods and validation. Expert Opinion on Drug Discovery, 1, 459–476. Martinelli, A., & Tuccinardi, T. (2008). Molecular modeling of adenosine receptors: New results and trends. Medicinal Research Reviews, 28, 247–277. Ortore, G., Tuccinardi, T., Bertini, S., & Martinelli, A. (2006). A theoretical study to investigate D2DAR/D4DAR selectivity: Receptor modeling and molecular docking of dopaminergic ligands. Journal of Medicinal Chemistry, 49, 1397–1407. Peeters, M. C., Li, Q., van Westen, G. J., & Ijzerman, A. P. (2012). Three “hotspots” important for adenosine A(2B) receptor activation: A mutational analysis of transmembrane domains 4 and 5 and the second extracellular loop. Purinergic Signalling, 8, 23–38. Peeters, M. C., van Westen, G. J. P., Guo, D., Wisse, L. E., Mu¨ller, C. E., Beukers, M. W., et al. (2011). GPCR structure and activation: An essential role for the first extracellular loop in activating the adenosine A2B receptor. The FASEB Journal, 25, 632–643. Sali, A., & Blundell, T. L. (1993). Comparative protein modelling by satisfaction of spatial restraints. Journal of Molecular Biology, 234, 779–815. Schiedel, A. C., Hinz, S., Thimm, D., Sherbiny, F., Borrmann, T., Maass, A., et al. (2011). The four cysteine residues in the second extracellular loop of the human adenosine A2B receptor: Role in ligand binding and receptor function. Biochemical Pharmacology, 82, 389–399. Thompson, J. D., Higgins, D. G., & Gibson, T. J. (1994). CLUSTAL W: Improving the sensitivity of progressive multiple sequence alignment through sequence weighting, position-specific gap penalties and weight matrix choice. Nucleic Acids Research, 22, 4673–4680. Tuccinardi, T., Ortore, G., Manera, C., Saccomanni, G., & Martinelli, A. (2006). Adenosine receptor modelling. A1/A2a selectivity. European Journal of Medicinal Chemistry, 41, 321–329. Tuccinardi, T., Schenone, S., Bondavalli, F., Brullo, C., Bruno, O., Ranise, A., et al. (2008). Substituted pyrazolo[3,4-b]pyridines as potent and selective A1 adenosine antagonists: Synthesis, biological evaluation and development of an A1 bovine receptor model. ChemMedChem, 3, 898–913.

Molecular Modeling of Adenosine Receptors

59

Tuccinardi, T., Zizzari, A. T., Brullo, C., Daniele, S., Musumeci, F., Schenone, S., et al. (2011). Substituted pyrazolo3,4-bpyridines as human A(1) adenosine antagonists: Developments in understanding the receptor stereoselectivity. Organic & Biomolecular Chemistry, 9, 4448–4455. Verdonk, M. L., Cole, J. C., Hartshorn, M. J., Murray, C. W., & Taylor, R. D. (2003). Improved protein-ligand docking using GOLD. Proteins, 52, 609–623. Xu, F., Wu, H., Katritch, V., Han, G. W., Jacobson, K. A., Gao, Z.-G., et al. (2011). Agonist bound structure of the human adenosine A2a receptor. Science, 332, 322–327.