Combined QM/MM Study of Thyroid and Steroid Hormone Analogue Interactions with αvβ3 Integrin

Recent biochemical studies have identified a cell surface receptor for thyroid and steroid hormones that bind near the arginine-glycine-aspartate (RGD) recognition site on the heterodimeric αvβ3 integrin. To further characterize the intermolecular interactions for a series of hormone analogues, combined quantum mechanical and molecular mechanical (QM/MM) methods were used to calculate their interaction energies. All calculations were performed in the presence of either calcium (Ca2+) or magnesium (Mg2+) ions. These data reveal that 3,5′-triiodothyronine (T3) and 3,5,3′,5′-tetraiodothyroacetic acid (T4ac) bound in two different modes, occupying two alternate sites, one of which is along the Arg side chain of the RGD cyclic peptide site. These orientations differ from those of the other ligands whose alternate binding modes placed the ligands deeper within the RGD binding pocket. These observations are consistent with biological data that indicate the presence of two discrete binding sites that control distinct downstream signal transduction pathways for T3.


Introduction
Integrins are plasma membrane proteins that generate complex and intracellular signals during morphogenesis, tissue remodeling, and repair [1]. Interactions of integrins with specific extracellular matrix proteins and growth factors or with certain small molecules generate ligand-specific signals. The integrin superfamily of heterodimeric glycoproteins consists of αand β-monomers that associate in defined combinations that include more than twenty different mammalian subtypes, of which at least eight are characterized by the presence of an Arg-Gly-Asp (RGD) site used in the process of recognition of different protein ligands [1][2][3]. The combination of αvβ3 subunits is the most abundant integrin expressed in mammalian cells [1].
Structural data have revealed that the extracellular segment of integrin αvβ3 (i.e., lacking the transmembrane and cytoplasmic tails) is V shaped with the 4-domain αv subunit and the 8-domain β3 subunit bent by 135 • [4][5][6]. Data for the RGD cyclic pentapeptide complex reveals that the RGD peptide occupies a shallow crevice between the propeller and βA domains in the integrin head.
The X-ray structure of apo αvβ3 has been determined in the presence of Ca 2+ ions [6] and in the presence of Mn 2+ as a complex with the cyclic pentapeptide Arg-Gly-Asp-{D-Phe}-{N-methyl-Val} ( Figure 1) [4]. These data showed that replacement of Ca 2+ with Mn 2+ in the αvβ3 integrin did not result in any important structural changes to the protein.
There are also relatively small changes in the αvβ3 structure with and without the cyclic pentapeptide bound in the RGD site of integrin.
Ligand binding to integrins is dependent upon bivalentcation interactions that are generally stimulated by Mg 2+ or Mn 2+ and inhibited by Ca 2+ . Competitive binding data 2 Journal of Biomedicine and Biotechnology revealed preferential binding sites for one or more cations that could also be displaced by other cations and that their specific occupancy can lead to enhanced ligand binding in an allosteric or synergistic manner [1,7]. The central core of the β-propeller is one of the primary sites for cation binding. A common structural motif observed for integrins is the presence of a metal binding site, termed metal ion-dependent adhesion site (MIDAS), that is unique for Mg 2+ /Mn 2+ and is located in the integrin α subunit I domain. Mutagenesis data have shown that the MIDAS motif and the nearby solvent exposed side chains are required for ligand binding [3].
The integrin protein family has been studied both experimentally and theoretically [8][9][10][11][12][13][14][15]. Studies using NMR and molecular dynamics methods to map the interaction of the cyclic β-peptide isoDGR, a competitive antagonist of RGD-containing ligands of αvβ3 integrin that inhibits endothelial cell adhesion, proliferation, and tumor growth, revealed a new recognition motif of the RGD binding pocket of αvβ3 integrin [9]. Docking studies for a series of cyclic peptides further defined the structure-function relationships in binding to the RGD recognition site of αvβ3 integrin.
Studies using a combination of experimental and computational approaches defined the role of divalent cations in mediating the binding of integrins to RGD peptides. Structural data revealed that the divalent cation at the MIDAS site in the βA domain of integrin contacts the Asp of the RGD peptide. This site is occupied by the divalent ion only in the RGD-liganded state [15]. The results of these molecular dynamics studies support the stability of the orientation of Asp residue of the RGD recognition site that interacts with the MIDAS which is essential for binding, whereas the orientation of Glu220 carboxyl oxygens of the β3 domain for the ligand associated metal binding sites (LIMBSs) was maintained throughout the molecular dynamics simulations [10,15]. The steered molecular dynamics studies of the unbinding of RGD cyclic peptides from αvβ3 integrin revealed the importance of single water molecules tightly coordinated to the MIDAS ion. These data suggest that this process is a general strategy for the stabilization of protein-protein adhesion against cell-derived forces through divalent cations [15].
Recently, a cell surface receptor for thyroid hormone has been identified on the extracellular domain of the αvβ3 integrin [16] that is relevant to the proliferative [16][17][18][19] and proangiogenic activity of thyroid hormone [16,20,21]. Transduction of the hormone signal into proliferation and angiogenesis at the integrin is achieved by the mitogenactivated protein kinase (MAPK; extracellular regulated kinase 1/2, ERK1/2) pathway [16,19,22,23] or by the phosphatidylinositol 3-kinase (P3K)/Akt pathway [24][25][26]. L-Thyroxine (T 4 ) and 3,5,3 -triiodo-L-thyronine (T 3 ) are agonist thyroid hormone analogues at the integrin receptor, whereas a deaminated derivative of T 4 , tetraiodothyroacetic acid (T 4 ac), is a nonagonist inhibitor of the binding of T 4 and T 3 to integrin. Competition data revealed that RGD peptides also block hormone binding by integrin, suggesting that the hormone-binding site is near the RGD recognition site on αvβ3. Further analysis revealed that RGD peptides could inhibit the binding of other biologically active small molecules with estrogen-like structures, such as cis-and trans-resveratrol and stilbene ( Figure 2) [18,[27][28][29]. These data reveal that T 4 promotes cell proliferation in cancer cells whereas resveratrol is proapoptotic and that further shows that T 4 inhibits resveratrol-induced apoptosis [18]. This action is prevented by T 4 ac, shown to be an inhibitor of T 4 binding to integrin. In the presence of T 4 , T 4 ac does not affect the induction of apoptosis by resveratrol [18,30]. These observations suggest that there is a discrete binding site for resveratrol and thyroid hormone on integrin.
The noniodinated thyroid hormone analogue GC-1 lacks the amino group of the alanine side chain, lacks the ether linkage between phenyl rings, and has methyl groups at the 3,5-positions and an isopropyl group at the 3 -position, compared with the thyroid hormone T 3 ( Figure 2). This thyroid hormone analogue is a selective agonist for nuclear thyroid hormone receptor TRβ1 and has the desirable metabolic effect of lowering circulating cholesterol without inducing tachycardia [31,32]. GC-1 was also demonstrated to be proangiogenic at the plasma membrane and require interaction with αvβ3 integrin [20].
To understand how these hormones and their analogues interact with integrin, modeling studies were carried out using the coordinates of the extracellular fragment of αvβ3 integrin (1L5G.pdb) [4,5] as a starting model for their interactions near the RGD recognition site [27]. These preliminary models revealed that the hormone analogues shown in Figure 2 can bind in the RGD recognition site and that there may be two closely related binding pockets that can bind these hormone analogues. The smaller of these sites may be preferential for the more planar compounds such as estradiol, stilbene, or resveratrol. Trans-Stilbene

Material and Methods
To provide more definitive details of the intermolecular interactions of these hormone analogues with the RGD recognition site in αvβ3 integrin, quantitative molecular modeling including molecular dynamics (MD) simulations and combined quantum mechanical-molecular mechanical (QM/MM) methods were employed. These calculations also incorporated the effects of solvent and divalent cations on ligand-binding interactions. In order to calculate the desolvation energies of these ligands, an explicit water environment based on molecular dynamics (MD) simulations was used in the QM/MM calculations.
In the present study we utilized our own method of combined QM/MM calculations which are based on MD. In this approach, standard MD of the protein solvated by a sphere of explicit water molecules was carried out. Then a randomly chosen set of snapshots was selected from the MD, and for each snapshot, full geometry optimization of the ligand (QM) in a fixed protein matrix (MM) was performed. In order to calculate the desolvation energy of the ligand, the same QM/MM was used to calculate the ligand surrounded by a sphere of explicit water molecules. The ligand-binding energy to the protein was calculated as a difference between the average ligand energy in the protein and the water environment. These computations were performed using the quantum-chemical program Q-Chem [33], which includes a modified version of an interface between QM and MM atoms [34]. Since the combined QM/MM calculations for each snapshot are independent of the other snapshots, this method referred to as the large-scale QM/MM method has the advantage of perfect parallel scalability [35]. All calculations were performed using a large 1.000 node Linus cluster where each node consists of two 3.2 GHz Xeon processors and 2 GB of shared memory.

Computational Details
Coordinates from the X-ray structure of the αvβ3 receptor (1L5G PDB) were used as a starting model [4]. The cyclic RGD pentapeptide was removed from the file and original Mn 2+ ions were replaced by Mg 2+ or Ca 2+ ions. Parameters of the NAG residues observed in this structure were obtained from the ANTECHAMBER program of the AMBER programming suite [36]. All hydrogen atom positions were populated based on the AMBER force field. A total of 25 Na + cations were used to neutralize the entire protein. The positions of these cations were obtained from grid calculations of the protein electrostatic potential [37,38]. The ligand structures ( Figure 2) were obtained using quantum mechanical (QM) calculations at the B3LYP/6-31+G * level of theory. The QM calculations consisted of full geometry optimization of the ligands performed in the gas phase. These calculations, including iodine atoms of the thyroid hormones (T 4 , T 3 , T 4 ac), were obtained using the lanl2dz effective core potential with a corresponding basis set [38].
Docking calculations were carried out using the Auto Dock 3 program [39]. Because of the lack of molecular mechanical parameters for iodine in the Auto Dock program, this atom was temporary replaced in the docking process by a carbon atom. The docking procedure was performed using a fixed internal geometry of these ligands. After docking, one snapshot of these compounds was randomly selected for subsequent study. The precise selection of a docked structure is not important because geometry and orientation of the molecule inside the protein active site were recalculated in the following MD simulations and in the QM/MM calculations. Once these ligands were docked to the protein near the RGD recognition site of αvβ3 integrin, another round of MD simulations was carried out. Parameters of the ligands were obtained using the ANTECHAMBER program in the AMBER suite.
The molecular system, including the protein, the ligand, and the metal ions, was placed into a center of TIP3P water sphere of radius 80Å. The MD consisted of 100 ps heating dynamics from 0 to 300 K, followed by equilibration dynamics performed for 5 ns. The MD simulation was also performed at constant temperature and volume, with application of constrained harmonic potentials outside the water sphere and for the metal ions. The MD has been performed separately for each ligand ( Figure 2) and separately for the protein system with Ca 2+ and Mg 2+ ions. Following the MD simulations, combined QM/MM calculations of the interaction between the ligand and the protein environment were carried out. First, a randomly selected protein snapshot from MD was chosen that was then divided into QM part and MM contributions. The QM part included the ligand and the MM part included the rest of the protein system with water molecules and the metal ions. Finally, full geometry optimization of the QM part inside the fixed positions of all atoms of the MM part was performed using the B3LYP/6-31+G * /AMBER method [34], as implemented in the Q-Chem program [33]. During geometry optimization, the ligands were allowed to change their internal geometry and position relative to the fixed protein matrix. The calculations of the ligands, including those with iodine atoms (T 4 , T 3 , T 4 ac), were performed using the lanl2dz effective core potential with a corresponding basis set.
The interaction energy between the ligand and the protein was calculated as an energy difference between the energy of the molecule inside the protein and the energy of the molecule in the gas phase which was obtained from the optimal geometry of the molecule in the protein and in the gas phase, respectively. The same procedure was then applied to 600 randomly selected protein snapshots from the last 200 ps of MD, and the results of the calculations are reported statistically based on the selected protein ensemble (Tables 1  and 2).
Calculations of an induced dipole moment of the ligands ( Figure 2) bound to the protein and in water were also carried out ( Table 2). The induced dipole moment is reported as a difference between the permanent dipole moment of the ligand in the protein (water) and the permanent dipole moment of the ligand in the gas phase, obtained from the optimal geometry of this molecule in the protein (water) and the gas phase, respectively. The calculated values of the interaction energies and induced dipole moments are reported as average values based on an ensemble of 600 protein (solvent) samples. Use of a larger number of samples in the ensemble would represent the protein (solvent) system much more realistically; however it requires more timeconsuming calculations. Therefore the number of samples (snapshots) chosen in the calculations was a compromise between computational accuracy and efficiency.

Computational Models of Ligand Interactions.
Structural data for the extracellular portion of αvβ3 integrin reveals Table 1: The interaction energy between the ligands (Figure 2) and the protein (or water). Numbers in parenthesis indicate the interaction energy of the ligand bound in the protein relative to the interaction energy of this molecule in water. The interaction energy is calculated as a difference between the energy of the ligand bound in the protein (or water) and the energy of this molecule in the gas phase, for the optimal geometry of this molecule bound in the protein (water) and the gas phase, respectively.   (Figure 2) bound in the protein (water). Numbers in parenthesis indicate the induced dipole moment of the ligand in the protein relative to the induced dipole moment of this molecule in water. The induced dipole moment is calculated as a difference between the permanent dipole moment of the ligand in the protein (water) and the permanent dipole moment of this molecule in the gas phase, for the optimal geometry of this molecule in the protein (water) and the gas phase, respectively. a heterodimer with 12 domains assembled into an ovoid "head" and two "tails." The αv and βA domains form a sevenbladed propeller and the β3 portion of the integrin "head" is composed of the βA and hybrid domains ( Figure 1) [4,5].

Ligands
The RGD cyclic peptide containing the RGD sequence binds in a crevice between the propeller and the βA domains of the integrin head (Figure 1) [5].  Figure 3 is the binding orientation of T 4 ac (red) and estradiol (yellow) as docked in the RGD cyclic peptide binding pocket.
Only small conformational changes were noted near the βA domain between the apo integrin structure and the RGD cyclic peptide complex [5]. These changes are also involved the metal (Ca 2+ or Mn 2+ ) binding site in this region.
To model thyroid or steroid hormone analogue interactions near the RGD cyclic peptide binding site, molecular docking experiments using quantum chemical calculations (QM/MM) were carried out in the presence of Ca 2+ or Mg 2+ for each of the ligands shown in Figure 2. Furthermore, these calculations were carried out with the ligand in two different starting positions near the RGD binding site (indicated as Ca1 and Ca2) to remove modeling bias from the calculations. Similar calculations were made with the ligand starting position about 15Å away from the RGD binding site. The results of the distant starting site calculations showed much smaller stabilization energies and suggest that the interaction between the ligand and metal cations is an important part of ligand binding (data not shown). The results of these quantitative calculations are compared with previously reported qualitative docking studies [27].

Thyroid Hormone Binding.
The results of the QM/MM modeling for the thyroid hormones, T 4 , T 3 , and T 4 ac, based on two starting orientations near the RGD binding site at the interface of the αv and β3 domains, show that the thyroid hormones can bind in a wide range of orientations. In the case of T 4 , the different starting positions result in binding modes in which the amino and hydroxyl ends of the molecule are flipped in the binding pocket. In one orientation, the 4 -hydroxyphenyl ring occupies a binding pocket deeper within the crevice between the αvβ3 domains, while in the other binding mode, the phenolic ring occupies the interface between the integrin αvβ3 domains (Figure 4). Consistent with pharmacokinetic/pharmacodynamic modeling reported elsewhere [26], the QM/MM modeling results for T 3 differ significantly from those of T 4 ( Figure 5). In one mode, the tyrosyl moiety lies along the Arg side chain of the cyclic RGD peptide ( Figure 6). The second binding mode for T 3 places the 4 -phenolic ring in a different binding pocket than observed for T 4 ( Figure 5). The results for the two starting positions of T 4 ac reveal yet another variation in binding modes (Figure 7). In this case, one of the T 4 ac orientations aligns along the side chain of Arg in the RGD peptide binding site, as was observed for one of the orientations of T 3 ( Figure 6). Comparison of the other starting position for T 3 and T 4 ac reveals a binding mode in which the phenolic ring of T 4 ac binds deeper within the RGD binding pocket, similar to that of T 4 . The phenolic ring of T 3 occupies an alternate pocket not occupied by either T 4 or T 4 ac (Figure 8). The one starting position calculated for the noniodinated thyroid analogue GC-1 binds in a similar manner to that of T 3 (Ca2) (Figure 9). The differences in binding modes possible for the thyroid hormones are in agreement with competition data that indicate that RGD peptides block hormone binding by αvβ3 integrin. These studies further show that T 4 ac is a nonagonist inhibitor of T 4 and T 3 binding to αvβ3 integrin.

Steroid-Like Ligands.
The results of the QM/MM modeling of ligand interactions with cis-and trans-reservatrol based upon two starting positions show that these ligands have similar binding orientations near the RGD binding site ( Figure 10). In the case of the cis-resveratrol, both starting positions placed the ligand in similar orientations such that the dihydroxyphenyl ring makes interactions with Asn215 of the β3 domain. In the case of trans-resveratrol, the two starting positions result in the molecule being flipped such that the dihydroxyphenyl ring is bound the opposite direction. The 4 -hydroxyphenyl ring still maintains contact to Asn215 of the β3 domain ( Figure 10). These results are in agreement with biochemical data that indicate only the β3 monomer was essential for activation and induction of apoptosis by resveratrol and that the occlusion of the RGD binding site blocks the cellular actions of resveratrol [30].
The modeling results for the QM/MM binding modes of cis-and trans-stilbene, the nonhydroxylated parent compound of resveratrol, revealed that the energy minimum conformation for the two starting positions of cis-stilbene were perpendicular to each other with one conformer bound deeper in the binding pocket than the other cis-conformer ( Figure 11). The results for the two starting positions of trans-stilbene show one molecule bound deeper in the RGD pocket than the other which is outside the binding pocket. The more hydrophobic nature of cis-and transstilbene compared with cis-and trans-resveratrol results in the greater variability of its binding modes as these molecules do not optimize binding involving specific hydrogen bonding interactions with residues in the RGD binding site.
The binding data for estradiol show that the 3-hydroxyl can interact with Asn215, as observed for resveratrol. An alternate orientation for estradiol binding shows the molecule perpendicular to that of the first and filling an alternate binding site (Figure 12).

Computational Models.
The interaction energies between thyroid hormone analogues and the RGD recognition site of αvβ3 integrin are presented in Table 1 that also shows the interaction energies between these ligands and water. As illustrated in Table 1, all calculated interaction energies in the protein are more negative than the corresponding interaction energies of those molecules in the water environment, which indicates that these molecules are more strongly stabilized by the protein environment than by the water environment. Among these ligands, estradiol shows the strongest interaction energy with the protein containing Mg 2+ ions, namely, −72 kcal/mol, while GC-1 shows the strongest interaction energy with the protein containing Ca 2+ , namely, −33 kcal/mol. Both values of the interaction energy include the calculated desolvation energies obtained in this study. T 3 shows a relatively strong interaction with the αvβ3 integrin and the type of metal binding present; namely, the interaction energy with the protein containing the Mg 2+ ions has a value of −60 kcal/mol, and the interaction energy with the protein containing the Ca 2+ ions has a value of   It is interesting to note that the interaction energy of T 3 with the protein containing both divalent metal ions is larger than the interaction energy of T 4 , even though both compounds have similar interaction energies with water. This observation indicates that an additional iodine atom in T 4 actually decreases its interaction energy with the protein. Moreover, the interaction energy of T 4 is larger than the interaction energy of T 4 ac indicating that the NH 2 -group which is present in T 4 and is absent in T 4 ac increases the interaction energy between the ligand and protein. Generally, the protein environment with the Ca 2+ ions stabilizes ligands such as resveratrol and stilbene, stronger than the protein environment with the Mg 2+ ions. However, for the other ligands studied the protein environment with Mg 2+ stabilizes the ligand more strongly than the protein environment with   greater opportunity to penetrate deeper into the protein active site. Deeper penetration of the protein active site by estradiol, GC-1, T 3 , T 4 , and T 4 ac, increases the interaction between the ligand and the divalent cations, particular for Mg 2+ , and this effect is observed in these calculations. Table 2 shows the values of the induced dipole moment of the ligands, obtained in the protein and water environments. The largest induced dipole moment is observed for estradiol bound in the presence of Mg 2+ , which has a value of 3.8 debye relative to the induced dipole moment of this ligand in water. The similar increased induced dipole moment in the protein with Mg 2+ is also observed in these calculations for T 4 .
Although it is not systematically observed for all ligands, there is a correlation between a value of the induced dipole moment of the ligand and its interaction energy with the surrounding environment that is caused by two effects. The first is related to the change of the molecular geometry by steric interactions between atoms of the ligand and atoms of the active site pocket, or atoms of a first solvation shell of water. The second effect is related to polarization of the electronic density of the molecule by an electrostatic field of the protein or water surroundings. These calculations show the second effect is more important. The change of the electronic density of the ligand inside the protein or water, observed as an induced dipole moment, generates a stronger attraction between polarized electrons of the ligand and the electrostatic field of the protein or water. As a consequence, the stronger attraction increases the total interaction energy between the ligand and the protein or water. The interaction energies show a strong stabilization effect of the bound ligand to integrin. However, these QM/MM calculations do not include polarization effects of the protein or water by the ligand that is inside these two environments. The QM/MM calculations have been performed using an approach in which all atoms of the protein or water are represented by molecular mechanical atoms, and therefore they cannot be polarized. Consequently, if polarization of the protein (or water) is included in the QM/MM calculation, the calculated interaction energy between the ligand and the surrounding environments is expected to be even larger. These studies reveal that there are two important effects responsible for the stabilization of the small molecules inside the active site of integrin binding domain. The first effect is related to internal flexibility of the molecule which yields better penetration of the molecule into the protein active site, and consequently it increases its stabilization. The second effect is related to polarization of electronic density of the molecule, by the electrostatic field of the protein (or water). Greater molecular polarization is responsible for a larger induced dipole moment, which consequently interacts more strongly with the surrounding protein (or water) environments. Therefore, it can be finally concluded that both molecular flexibility and molecular polarizability of the ligands are the main forces responsible for strength of the binding process of the small molecules to αvβ3 integrin.

Conclusions
Recent data have shown that the hormone receptor on αvβ3 integrin enables the thyroid hormones T 4 and T 3   to stimulate cancer cell proliferation and angiogenesis and to regulate the activity of certain membrane pumps [28]. The results of these QM/MM modeling studies for a series of thyroid or steroid-like ligands (Figure 2) bound to the αvβ3 integrin reveal that these ligands have multiple modes of binding near the RGD peptide binding site that was shown to occupy a shallow crevice between the propeller and βA domains of the integrin head [4] (Figure 1). The most notable observation from these studies is that the thyroid hormone analogues probe significantly different conformational space near the RGD recognition site and that these modes can be correlated with their biologic response to signaling via the integrin. Furthermore, the binding modes for T 3 and T 4 ac are consistent with pharmacodynamic data that indicate the presence of two discrete binding sites for T 3 that control distinct downstream signal transduction pathways [28,40]. Comparison of the low-energy binding modes for thyroid hormones suggest that the positions of T 4 and T 4 ac illustrated in Figure 8 are consistent with data that shows T 4 ac blocks T 4 binding to integrin. The different binding positions of T 3 ( Figure 5) are also consistent with data that suggest two discrete, but close, binding sites that activate different downstream signaling pathways [28].
The binding of resveratrol has been shown to interact primarily with the β3 monomer of integrin, activating the ERK1/2 pathway [18,30]. These QM/MM modeling data for cis-and trans-resveratrol indicate that both conformers have intermolecular interactions with the side chain residues of the β3 domain, although the cis-conformer which is flipped end to end with those of the trans conformers can potentially be involved in more hydrogen bond interactions ( Figure 10). Furthermore, comparison of the binding of resveratrol with T 4 and T 4 ac ( Figure 13) suggests a mechanism of ligand recognition that is compatible with the observed downstream biological responses.
The QM/MM modeling data presented here show that these ligands have flexible binding modes and that some ligands can bind in either a forward or reverse direction near the RGD recognition site. This flexibility in their low energy conformations suggests that each ligand can selectively take advantage of specific interactions within the local environment and with the metal binding sites. The observation of multiple binding pockets near the RGD recognition site also indicates that there may be subtle changes in the integrin binding site that can dictate binding to alternate sites that influence downstream signaling pathways.