| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
* Department of Biosciences and Nutrition, Karolinska Institutet, and Karo Bio AB, Novum, SE-141 57 Huddinge, Sweden
Correspondence: Address reprint requests to Lennart Nilsson, E-mail: lennart.nilsson{at}biosci.ki.se.
 |
ABSTRACT |
|---|
 TOP ABSTRACT INTRODUCTIONMATERIALS AND METHODSRESULTSDISCUSSIONCONCLUSIONSSUPPLEMENTARY MATERIALACKNOWLEDGEMENTSREFERENCES |
|---|
|
INTRODUCTION |
|---|
|
TOPABSTRACTINTRODUCTIONMATERIALS AND METHODSRESULTSDISCUSSIONCONCLUSIONSSUPPLEMENTARY MATERIALACKNOWLEDGEMENTSREFERENCES |
|---|
) to explicit evaluations of the absolute free energy of binding (2). Equally important properties in the drug development process are selectivity and efficacy; and resolution of the free energy of binding into its kinetic components (on- and off-rates) may provide a way of understanding and predicting these phenomena (3). Specifically, experimental studies of binding kinetics based on the surface plasmon resonance technique show distinct separation of on-rates between estrogen receptor agonists and antagonists (4). Thus, detailed knowledge of how ligands enter and exit the receptor adds useful information to the drug development process. Theoretical modeling of kinetic constants requires identification and characterization of the binding pathways to find the most plausible routes for the ligand from solvent to binding pocket.
A crystallographic structure model of receptor with bound ligand offers a well-defined starting point for identification of unbinding pathways. Assuming microscopic reversibility, the results from an unbinding simulation are valid also for the reverse process. Molecular dynamics (MD) simulations describe the evolution over time of a system in atomic detail, but the application of MD on macromolecular systems is severely limited by the short time ranges (tens of nanoseconds) which are routinely accessible today. Biologically relevant events such as rearrangement of tertiary structure in proteins or ligands entering a binding site typically occur on a much longer timescale and will therefore not be observed in ordinary MD simulations. Several techniques have been presented to overcome this limitation, e.g., locally enhanced sampling molecular dynamics (ESMD) (5), steered molecular dynamics (SMD) ((6) and references therein), and random expulsion molecular dynamics (REMD) (7).
In REMD, an artificial force of constant magnitude acting on the bound ligand is added to the standard force field to accelerate the unbinding process. The direction of the force is randomly varied to ensure that the ligand searches for many possible escape routes during a relatively short timespan. REMD was developed on and successfully applied to cytochrome P450 systems (7) but suggested to be generally applicable to ligand-receptor systems. Originating from simulations of atomic force microscopy applications, SMD is an alternative approach which also applies an artificial force to a specific part of the system but in a predetermined direction. SMD was used to study the unbinding of RA from RAR (6), and several plausible binding routes for the ligand were investigated. SMD allows a more detailed mapping of the energy profile along the specified pathway compared to REMD but is limited in ligand-binding studies by the subjective choice of pathway direction. Consequently, Lüdemann et al. used REMD for an unbiased search for pathways (7) and followed up the analysis by applying SMD on the resulting pathways (8).
In this work, we describe the application of REMD to study unbinding of RA from RAR. Our focus is on the identification of pathways and evaluation of a computational method not previously applied to the nuclear hormone receptor family. The RAR system was selected for several reasons. First, conclusions from this model should be applicable to other nuclear receptors, including many important drug targets, due to the highly conserved tertiary structure within this superfamily. Second, this system has been studied with alternative computational techniques (6,9), allowing methodological comparisons. Third, the extended and flexible nature of RA in combination with the presence of a strong electrostatic interaction with the receptor is a challenge for the REMD method and would enhance our understanding of the application of REMD. Last, high quality starting structures are available through the Protein Data Bank (10). The physical characteristics of RA prompted us to undertake a specific study of the REMD methodology as applied to this type of ligand. For this purpose we constructed a minimal two-dimensional (2D) model system in which the influence of parameters such as force application scheme, intrinsic flexibility, and charge state could be controlled. The very small number of particles in this model system enabled extensive testing because of the short simulation times required. Conclusions from this study assisted interpretation of our results from the full RAR model.
|
MATERIALS AND METHODS |
|---|
|
TOPABSTRACTINTRODUCTIONMATERIALS AND METHODSRESULTSDISCUSSIONCONCLUSIONSSUPPLEMENTARY MATERIALACKNOWLEDGEMENTSREFERENCES |
|---|
). Hydrogen atoms were added by the HBUILD (12) command in CHARMM (13). The structure was relaxed to the CHARMM22 (14) all-atom force field by two rounds of adopted-basis Newton-Raphson minimization: first 500 steps with all 
-carbons harmonically restrained (force constant 20.0 kcal/mol/Å2) followed by 1500 steps without restraints. The complex was solvated in a 20-Å-radius sphere of TIP3P (15) water, where all noncrystallographic water molecules having the oxygen atom closer than 2.8 Å of any protein atoms were removed. A very brief, unconstrained steepest descent minimization for 20 steps finalized the model setup. SHAKE was used to constrain bond lengths and angles during the MD simulations.
Two-dimensional model system
Carbon atoms of type CT3 (14) were used to define a ligand and a wall or a cage (Fig. 1, a and b). The uncharged wall and cage atom centers were placed 5.0 Å apart, not connected to each other, and positionally restrained by a harmonic potential (force constant 5.0 kcal/mol/Å2). The ligand was modeled by 4, 7, or 15 carbon atoms arranged as an extended alkane chain, using C-C bond lengths and C-C-C angles from the force field. The ligand atoms had no partial charge unless otherwise stated. A circular boundary potential was used in the wall model to continue the ligand movements in the x, y plane. Verlet dynamics was carried out with all coordinates and velocities in the z-direction zeroed.
|
|
) for the rigid ligand structure using the HF/6-31G* basis set and the electrostatic potential (ESP) (17) charge-fitting algorithm. Force field parameters for RA are listed in the Supplementary Material (Tables S1 and S2). In the initial setup the tolerance for energy change between MD steps was set to 2000 kcal/mol, and a time step of 0.002 ps was used.
In the final setup the constraints were moved to residues 325–336, which correspond to helix 8 in the RAR model. The initial HF/6-31G* RA ligand charges were restrained by the RESP (18) program as provided from the AMBER home page (19) to average charges on methyl group hydrogen atoms and reduce the overall magnitude of partial charges. It was assumed that this procedure yielded a sufficiently accurate electrostatic description of the ligand for use in conjunction with a strong external force. A stronger REMD force was used, up to 2000 pN, and the energy change tolerance level was set to 200 kcal/mol. The time step was decreased to 0.001 ps.
REMD-specific parameters
The magnitude of the REMD force was set to 10–30 pN per atom in the ligand, corresponding to 500–1000 kJ/mol/nm as applied on the center of mass of camphor by Lüdemann et al. (7). In our initial simulations, intervals of N = 10, 20, and 30 MD steps (N) were tested. The velocity of the ligand was periodically evaluated as the average velocity over the last N steps. The variation in key results such as fraction of successful expulsions was very small, and a fixed interval of 20 steps (20 fs) was used throughout the final simulation setup. The initial minimum velocity of 
v
min = 2 x 10–4 Å/fs led to very few changes of directions, and therefore the limit was increased fivefold to 1 x 10–3 Å/fs. The REMD simulation was aborted either when the ligand atom C6 had moved more than 26 Å from the starting position (approximately the radius of the receptor protein) or when the maximum simulation time of 200 ps was reached.
|
RESULTS |
|---|
|
TOPABSTRACTINTRODUCTIONMATERIALS AND METHODSRESULTSDISCUSSIONCONCLUSIONSSUPPLEMENTARY MATERIALACKNOWLEDGEMENTSREFERENCES |
|---|
,21) guided our setup of a RAR model for REMD simulations. RAR was solvated in a spherical water droplet to enable complete solvation of surface residues with a limited number of water molecules. The starting structure with the 12 helices (H1-H12) fold possessed characteristic receptor-ligand interactions as described for the crystallographic structure (see Fig. 4, a and d) (11). The important Arg-278, Ser-289, and Phe-288 contacts with the ligand were all present, and ligand binding pocket (LBP) residues on H1, H3, H5, and the ß-turn were also within van der Waals (vdW) contact distance of the ligand.
|
-carbon atom of the most central residue (Ile-270) were tested to prevent translation but allow the protein to rotate and thereby possibly facilitate expulsion through a curved pathway. However, all of these attempts resulted in failed simulations due to excessive energy changes between MD steps. A contrasting approach using weak restraints (force constant 0.1 or 0.01 kcal/mol/Å2) on all carbons resulted in fewer failures but also very few successful expulsion events. Initially, the most successful scheme was to apply moderate restraints (force constant 1.0) on the carbons of the distal H9, residues 348–370. This helix is far from the binding site as well as far from previously suggested pathways near H12 (6) and H7 (22). In the final setup for the RAR model, the central H8 (residues 325–336) was restrained instead of H9, since recent simulations of the structurally related glucocorticoid receptor (20) show that this is the least fluctuating helix. The optimal magnitude of the REMD force is a compromise between high expulsion efficiency (high force) and minimal structural distortion (low force), which we monitored by the change in total energy of the system in one MD step. Finding parameters which yield a reasonable number of expulsions within a given time limit, without violation of the energy change tolerance criterion, turned out to be nontrivial. To allow a sufficiently large REMD force necessary for the charged ligand with MD parameters close to the default values, the time step was reduced from 0.002 to 0.001 ps, but still a large number of simulations failed due to large energy changes between successive MD steps in the initial setup. A detailed examination of the model revealed large fluctuations in intermolecular electrostatic energies, and therefore the influence of various charge schemes was investigated. These are summarized in Table 1. The deprotonated ligand with 6-31G*/ESP charges, which was our starting point, yielded a large number of aborted simulations due to violation of the energy tolerance. Protonation of the carboxylic moiety of RA reduced the magnitude of the partial charges in the carboxylic end of the ligand. This ligand model yielded no expulsions when a weak force (200–500 pN) was applied to C6 or C15 and 100% violation of the energy change tolerance criterion when a force of 750 pN was used. Increasing the energy change tolerance 1000-fold to 20,000 kcal/mol led to successful expulsions in all cases with protonated RA and 750-pN pull force on C6 or C15. To clarify if the charge-charge interactions between RA and protein residues caused these rapid energy changes, a completely uncharged, protonated RA was tested. With 750-pN force and default energy change threshold of 20 kcal/mol no simulations crashed, and 60–90% yielded ligand expulsion within 200 ps. The same setup for a deprotonated RA using restrained electrostatic potential (RESP) charges based on 6-31G* calculations resulted in few expulsions, even during 400 ps of simulation, but no energy violation. Further methodological discussion about the default MD parameters in a REMD simulation is given in Supplementary Material Part S4.
|
The two-dimensional model system
The details of the wall and cage geometries (Fig. 1, a and b) are given in the Materials and Methods section. Our model ligand was a carbon atom chain of varying length where the parameters for bond lengths, angles, and vdW radius were taken from aliphatic carbons in the force field. The expulsion ratio was measured as the fraction of successful expulsions from 30–50 simulations for every new parameter setup. A rough error estimate for the expulsion ratio based on the standard deviation in repeated simulation series is 5% (Supplementary Material, Table S3).
The simplest model system consisted of a wall with a gap and a ligand restrained by a circular potential (Fig. 1 a). Three ligand sizes were tested (4, 7, and 15 atoms), and the gap was positioned at the center of the wall. For each ligand, three different simulations were carried out: ordinary MD without REMD force, REMD with the force distributed on all atoms, or REMD with the external force applied to one atom at the end of the molecule. The successful expulsion ratios are presented in Table 2. Our results from the centered gap simulations show that the force application mode is of minor importance for small ligands, whereas the expulsion ratio is increased from 70% to nearly 100% for the extended ligands when the force application is changed from all atoms to a single atom. In the centered gap model it is easier to pull out a large ligand when the force is applied to one atom. Small ligands tumble around in the model cavity and therefore exhibit a lower expulsion frequency. Note that in this model the ends of the ligand are initially equally far from the exit gap. Typically the ligands will rotate and exit with the pulled atom first. This is particularly pronounced for the large ligand where all expulsions take place in this manner. In addition, the comparison with ordinary MD shows that REMD is particularly efficient when applied to large ligands.
|
|
|
|
|
|
A second pathway was found between H7, the ß-sheet loop, and the H1-H3 loop (Pathway B, Fig. 4 e). This part of the receptor is denser, but still the ß-sheet loop and the H1-H3 loop provided sufficient receptor flexibility. Pathway B was observed when the REMD force was applied to C10 or C15, but at low frequencies. During expulsion along pathway B, mainly H5 and S1 residues such as Thr-277, Thr-280, and Thr-287 were within 3.0 Å. The ligand also passes by the charged Arg-278 and polar Ser-289, crucial for the initial positioning of the ligand in the LBP through the interaction with the RA carboxyl group. Initially attached water molecules stay in close interaction with the ligand during unbinding through pathway B.
A third and unexpected exit pathway was observed in the sterically crowded region between the ß-sheet, H7, and H8 (Pathway C, Fig. 4 f). Escapes along this pathway often caused major rearrangements of the local receptor structure since the ligand had to make way for its exit (Fig. 5 a). This pathway was observed only when the force was applied to RA atom C6. Expulsion along pathway C involved contacts with a number of residues in H5 (Leu-273 and Cys-276) and H7 (Phe-312, Ala-315, Gly-316, and Leu-319) residues. Also Met-286 in S1 and residues at the loop connecting H5-S1 are contacted (Tyr-279 and Pro-281). The escape through pathway C causes a reduction of the solvent interactions at the ligand carboxy end. The initially attached water molecules stay in the LBP, whereas the ligand is expelled through the narrow pathway.
|
|
DISCUSSION |
|---|
|
TOPABSTRACTINTRODUCTIONMATERIALS AND METHODSRESULTSDISCUSSIONCONCLUSIONSSUPPLEMENTARY MATERIALACKNOWLEDGEMENTSREFERENCES |
|---|
Our pathway A is similar to a possible exit pathway suggested by Kosztin et al. (6). Their SMD characterization of this pathway shows a low energy profile and little receptor deformation compared to a pathway similar to our pathway D as well as a proposed entry site between H11 and H12 which was never observed in our simulations. In the proposed mechanism for ligand escape via pathway A as defined by Kosztin et al. (6), the carboxyl end exits first. No attempt to rotate the ligand molecule before the expulsion was made. In contrast, our results indicate that both the RA head (C6) and tail (C15) might be the first parts to exit. The former case suggests an alternative binding mechanism where the carboxylate end of RA is attracted to the polar "window" (6) on the receptor surface and anchored by electrostatic interactions, followed by an incorporation of the ß-ionone ring into the hydrophobic LBP driven by desolvation of RA. Given this order of events, the forced unbinding by this mechanism observed in our simulations seems less likely than the binding of RA through this pathway, thus suggesting pathway A to be an entry pathway. In any case, this region of the nuclear receptor structure has also been proposed to host a ligand binding/unbinding pathway for 3,5,3'-L-triiodothyronine to the thyroid hormone receptor (22,23).
Pathway B, which is located on the other side of the H1-H3 loop, is similar to pathway A but more sterically hindered, resulting in a narrower and less varied pathway. Exit through pathway B is less frequently observed in our simulations and to our knowledge not previously described for RAR. However, recent modeling studies have proposed the existence of an alternative binding pocket of nuclear receptors (24) related to the fast biological functions observed for some nuclear receptors. The kinetics for the fast biological response is several orders of magnitude faster than the genomic response and may therefore require a more accessible LBP. Specifically, for the vitamin D receptor and the estrogen receptor, the alternative LBP was proposed to be located between the H1-H3 loop and the ß-hairpin loop. Interestingly, this would correspond to the portal of our pathway B, and it is feasible that ligands which are transiently bound to the alternative LBP occasionally enter the classic LBP through this pathway.
Pathway C is found in a sterically dense region of the receptor. The observed exit frequency for this pathway is low and only observed while pulling the ligand at atom C6 in the ß-ionone ring of RA. Due to the limited space, initially attached water molecules are left in the LBP. This desolvation should increase the energy profile for this pathway. The protein is severely deformed in these simulations, especially in the region close to the exit pathway where both the sheet section and the H7 move 
5 Å (Fig. 5 a). Our simulations involving pathway C indicate an unphysical unbinding sequence where the hydrophobic part of RA enters the solvent before the hydrophilic end, or a binding sequence where the carboxylic moiety is buried in the hydrophobic interior of the protein before the entry of the hydrophobic part of RA, and may therefore be an artifactual result of the strong REMD force. In contrast to pathway A, no electrostatic field is present in this region which could guide the association of the ligand-receptor complex. More physically sensible is a binding mechanism where the ß-ionone ring enters the hydrophobic interior through pathway C, and the position of the carboxylate near Arg-278 is subsequently obtained by a rearrangement of the ligand position within the LBP. This rearrangement could in principle be shared with the initial step of the two-state unbinding mechanism along pathways near H11 and H12, as suggested by others (5).
Pathway D was only occasionally observed in the final setup but more frequently in the initial setup. We believe that the reduced energy tolerance levels of the main setup limit the number of expulsion events through this pathway, which requires large structural rearrangements. However, our 2D model simulations indicate that this pathway is disfavored due to the application of the REMD force on an atom far from the exit and suggest that the actual frequency may be twofold higher than the observed frequency (see below). Kosztin et al. (6) also try to unbind RA through this pathway, but RA never leaves the receptor interior during a 750-ps simulation despite several attempts to vary the direction of the applied force and considerable deformation of the receptor structure. The linear pathway implied by SMD might not be appropriate if the exit route has a curved nature. REMD could in principle find curved pathways, but still pathway D is observed at a low frequency and the receptor structure is severely affected (Fig. 5 b).
In this study, our focus has been to identify all possible escape pathways and our REMD simulations do not allow a detailed quantification of the likelihood of each proposed pathway. The large conformational changes observed for pathways C and D show that ligand escape is sterically possible along these pathways, but due to the short simulation time they do not represent equilibrated conformations. Support for the existence of the proposed pathways comes from experimental studies on several nuclear receptors where point mutations of amino acids along our pathways have been shown to affect receptor function in general and binding kinetics in particular. Our pathway D exits between the H11-H12 loop and the N-terminal half of H3, and in this region mutations are studied in the estrogen receptor (25,26), the glucocorticoid receptor (27,28), the androgen receptor (29), and the thyroid hormone receptor (30). Many of these mutations comprise residues which govern the interaction between the H11-H12 loop and H3 and could therefore also modulate the binding kinetics for a pathway such as our pathway D. Further support comes from studies by Gee and Katzenellenbogen (31), who showed that the estrogen receptor may undergo partial unfolding in the LBP region. It was hypothesized that this transient conformational change influences ligand binding, and as such, the partial unfolding may confer larger conformational changes than observed for our pathway D.
Ligand expulsion effect on helix 12
H12 never leaves the canonical agonist position (32) in our REMD simulations, although it is not artificially constrained to this position. Since the REMD methodology does not favor any particular direction, we expected H12 to be pushed away if that would lead to a favorable unbinding pathway for the ligand in addition to the trivial case where H12 is displaced from the agonist position leaving the LBP wide open. Instead, the ligand finds several other ways out of the LBP in our simulations. Reversing the unbinding process, our results suggest that it is possible for an agonist ligand to bind to a nuclear hormone receptor in which H12 is already positioned in the agonist conformation. We believe that this could explain observations of large differences in on-rates between estrogen receptor agonist and antagonist ligands (4). Possibly, the larger antagonist ligands can only enter the binding site when H12 is not bound to the receptor surface in the agonist position as postulated by the mouse trap mechanism (11), whereas agonists may bind virtually any nonliganded conformation of the receptor. This is in line with the conformational selection model (33–35) where ligands continuously sample an ensemble of protein conformations and ligand binding results in selective stabilization of certain conformations leading to a shifted equilibrium of the conformational ensemble. In this model, antagonist ligands would sample a smaller subset of conformations, i.e., face a lower effective concentration of appropriate receptors, thus having a lower association rate.
The similar results from our REMD simulations and the SMD simulations by Kosztin et al. (6) are somewhat contrasted by ESMD results (9) where RA exits from RAR exclusively via pathway D, and the unbinding involves a two-step mechanism involving ligand dissociation from specific interactions in the LBP followed by a protein reorganization phase. Crucial for unbinding through this mechanism is the detachment of H12, which is not observed in our model or mentioned in the work by Kosztin et al. (6). The different results can be method dependent or model dependent. The ESMD method is conceptually similar to REMD in using an increased temperature factor for the ligand, and it is evident that the accelerated motion of the ligand affects the protein structure in both methods, but probably in different ways. The use of a single ligand copy in REMD causes a distinct effect on the protein structure at the collision site, whereas the multiple ligand copies in ESMD should lead to a more wide-spread protein-ligand interaction. Possibly the more broadly distributed but weaker protein-ligand interactions of ESMD simultaneously destabilize many interactions around the LBP, which eventually leads to the detachment of H12 making the exit pathway D more likely. The differences between these investigations can also be sought in the RAR model. The ESMD simulations were based on a homology model of the RAR receptor without addition of structural water molecules. Later, it has been shown that structural water molecules may indeed bridge H12 to the rest of the RAR receptor and thereby slow down the dissociation rate considerably (36). In a more recent ESMD study of the thyroid hormone receptor three ligand unbinding pathways are found, two of which do not require opening of H12 (22). All three pathways, however, deform various parts of the secondary structure.
Somewhat counterintuitive, the strong force used in our REMD simulations may actually be the reason for the limited deformation of the overall protein structure. Whereas individual side chains are mobile and quickly react to sterical changes in their surroundings, the secondary structure elements along the observed pathways may not have time to react to the rapidly moving ligand on the subnanosecond timescale. Indeed, the strong REMD force applied in this study makes most of our escape routes nearly straight (Fig. 3). To fully utilize the inherent capability of the REMD method to explore curved escape routes, one would need to remove strong electrostatic interactions which require a large force to break. A smaller force in a system dominated by weak short-range interactions would result in a greater probability for the ligand of halting and changing directions. In this context, it would be of interest to compare proposed pathways in different systems for ligands of varying size relative to its receptor. In our case, the ligand is approximately equal in size to a single secondary structure element such as a ß-strand or a medium sized -helix. Such large ligands may have predominantly linear routes to and from buried binding sites, whereas small ligands such as diatomic gases in myoglobin may more frequently use more complicated routes possibly involving intermediate binding sites. The overall shape of escape routes would then essentially be determined by the magnitude of the conformational change required for ligand travel.
Residues involved in the unbinding process
Imperative for an informative REMD simulation is a physically relevant starting structure, where protein-ligand interactions are consistent with available experimental data. However, the magnitude of the artificial force in REMD is much larger than the normal intermolecular interactions between ligand and receptor, and therefore we can only qualitatively judge the importance of individual residues located along the unbinding pathway using this method. We noted that the characteristics of the contacted amino acids differ along the four pathways. Along pathway A, there is a mix of polar, charged, and hydrophobic residues. Pathway B exhibits a highly hydrophilic character. Expulsion through this part of the receptor involves several charged and polar side chains, which suggests that solvent may play an important role for ligand expulsion in these directions. Pathways C and D are mainly lined by hydrophobic amino acids. Residues contacted along pathway C have rather small side chains, which may be required for ligand expulsion through this dense region of the receptor.
Effect of ligand size and force application site
The RAR system we are interested in differs from the cytochrome P450 system which Lüdemann et al. (7) studied first by REMD in several ways. Our primary concern was the RA ligand, which has an extended conformation with more flexibility than camphor, the main ligand in the cytochrome P450 study. Lüdemann et al. included the extended ligand palmitoleic acid and noted that it required a significantly larger rupture force due to the strong charge-charge interaction between the acid moiety and the protein, but it is not clear how and to what extent the flexible and extended structure of this ligand influenced the outcome of the REMD simulation. Should the force be applied to the center of mass or to individual atoms? Will application of force at one end of the ligand bias the result? How should the strong charge-charge interaction best be handled? To be able to study these kinds of questions we used a 2D model system, which offered precise control over simulation parameters and very short simulation times compared to the full three-dimensional protein model.
For an extended and somewhat flexible ligand, we expected different results when the force was distributed to all atoms compared to application of the full force on a single atom. The all-atom force distribution is, from a technical viewpoint, the same as application of the force to the center of mass. In terms of expulsion frequency, this is consistent with our observations of the basic 2D system as shown in Table 2. Whereas the small ligands are expelled in 60–70% of the simulations, independent of force application scheme, the expulsion ratio is increased from 70% to nearly 100% for the longest ligand when the force is applied to a single atom. As expected, application of the REMD force on one atom makes this atom the leading atom, and expulsion of the ligand will therefore take place with this atom first. We can interpret this by an analogy: it is easier to pull a trailer forward with your car, than reversing. The smaller ligands (4 and 7 atoms, respectively) are probably too small to have an effective "trailer", whereas the longest ligand has effectively one trailer in each direction when the REMD force is distributed evenly. When the force is distributed to all atoms, the orientation by which the ligand exits the gap is more random and close to 50% for both ends. This shows that our 2D model works as expected and sets a standard for expulsion ratios for further experiments.
REMD was developed to overcome the problem of sampling events such as ligand binding and unbinding that normally take place on timescales much longer than ordinary MD simulations. In our model system, we expected a very small ligand molecule to be able to diffuse out through the gap, essentially unhindered, whereas a large ligand should be trapped. Our results from the 2D simulations show that both the 4-atom ligand and the 7-atom ligand may diffuse out during an ordinary 200-ps MD simulation, whereas the longer ligand needs the assistance from an external force to pass the gap within that time. In this study, the magnitude of the REMD force has simply been scaled to the number of atoms in the ligand, independent of atom types. In the 2D model there is only one atom type, but when including hydrogen atoms as in the full protein model, inverse weighting of the force to the atomic mass may be an alternative. However, limited tests with realistic small molecules in our 2D model setting showed no marked differences such as structural deformations due to exaggerated forces on hydrogen atoms, which is why this may be a relatively minor issue.
The expulsion ratio is slightly decreased for the small ligands when the gap is moved to the far end of the wall (Table 3), whereas the largest ligand increases its ratio to almost 100% expulsion also with a distributed force. A plausible explanation is that the smaller ligands may tumble freely in the relatively large cavity defined by the wall and circular potential, but when the gap is in the middle of the wall, the small ligands have only about half the space for tumbling before they are likely to hit the gap. All expulsions with the largest ligand take place head first, showing that there are no complete rotations, thus the ligand mainly moves back and forth or follows the cavity borders.
Ligand flexibility and electrostatic interactions
Increasing the angular and dihedral force constants to create a more rigid ligand decreases the expulsion ratio (Table 4). A very rigid ligand has a lower expulsion frequency compared to a more flexible one, indicating that a certain degree of flexibility facilitates expulsion. This effect is particularly pronounced when the force is distributed over all ligand atoms. In contrast, for a flexible ligand the highest expulsion frequency is achieved when the REMD force is applied on a single atom. Our results indicate that an alkane chain built on standard force field parameters can be considered to be completely flexible in this context.
The introduction of a charged site in the system strongly decreases the expulsion frequencies, as seen in our 2D "cage" model (Table 5). The lowest expulsion rates are observed when the charge is placed at the interior of the system, opposite to the exit. Our results clearly show that the easiest way to disrupt an electrostatic interaction and get ligand expulsion is to apply the REMD force on the charged atom. Force application to other atoms resulted in lower expulsion frequency at the same level. Consequently, the exit closest to the charge will be favored, and in our 2D model study the expulsion ratio was 1:2 for the inner versus the exit charge, respectively. The actual frequency of escapes along a pathway where the ligand is anchored by an inner charge may therefore be underestimated by a factor of 2.
To summarize, our REMD simulations showed that the force application mode is a critical parameter for the outcome of a REMD simulation. A single atom application is preferred if the ligand is flexible and charged. There is a substantial bias for escape routes near the atom which is being pulled by the REMD force, but in a tight compartment the ligand might exit with an atom far from the force application site first. Given more room, the molecule will most likely rotate and exit with the pulled atom first. The relative orientation of the ligand in the exiting process is therefore highly dependent of the REMD method.
|
CONCLUSIONS |
|---|
|
TOPABSTRACTINTRODUCTIONMATERIALS AND METHODSRESULTSDISCUSSIONCONCLUSIONSSUPPLEMENTARY MATERIALACKNOWLEDGEMENTSREFERENCES |
|---|
REMD is a powerful yet not fully developed method for elucidating unbinding pathways. Here, we find that application of the REMD force to all ligand atoms (or the center of mass) is generally preferred. This application mode gives the most unbiased search for expulsion pathways, and the main advantage of REMD is then fully utilized. For a ligand which is tightly bound in the binding pocket through a strong electrostatic interaction, application of the REMD force on a single atom is preferred. The expulsion frequency will then increase, although the resulting pathways are biased toward the direction where the REMD force is applied.
Furthermore, we suggest that the problem of finding possible dissociation pathways is separated from the problem of elucidating the force required for dissociation by applying the REMD method to completely uncharged ligand models. The magnitude of the REMD force is best determined by systematic increments starting from a relatively low strength. More studies are required to find a general protocol for setting the REMD parameters.
|
SUPPLEMENTARY MATERIAL |
|---|
|
TOPABSTRACTINTRODUCTIONMATERIALS AND METHODSRESULTSDISCUSSIONCONCLUSIONSSUPPLEMENTARY MATERIALACKNOWLEDGEMENTSREFERENCES |
|---|
|
ACKNOWLEDGEMENTS |
|---|
|
TOPABSTRACTINTRODUCTIONMATERIALS AND METHODSRESULTSDISCUSSIONCONCLUSIONSSUPPLEMENTARY MATERIALACKNOWLEDGEMENTSREFERENCES |
|---|
Submitted on February 7, 2006; accepted for publication July 3, 2006.
|
REFERENCES |
|---|
|
TOPABSTRACTINTRODUCTIONMATERIALS AND METHODSRESULTSDISCUSSIONCONCLUSIONSSUPPLEMENTARY MATERIALACKNOWLEDGEMENTSREFERENCES |
|---|
2. Dixit, S. B., and C. Chipot. 2001. Can absolute free energies of association be estimated from molecular mechanical simulations? The biotin-streptavidin system revisited. J. Phys. Chem. A. 105:9795–9799.[CrossRef]
3. Markgren, P. O., W. Schaal, M. Hamalainen, A. Karlen, A. Hallberg, B. Samuelsson, and U. H. Danielson. 2002. Relationships between structure and interaction kinetics for HIV-1 protease inhibitors. J. Med. Chem. 45:5430–5439.[CrossRef][Medline]
4. Rich, R. L., L. R. Hoth, K. F. Geoghegan, T. A. Brown, P. K. LeMotte, S. P. Simons, P. Hensley, and D. G. Myszka. 2002. Kinetic analysis of estrogen receptor/ligand interactions. Proc. Natl. Acad. Sci. USA. 99:8562–8567.
5. Elber, R., and M. Karplus. 1990. Enhanced sampling in molecular-dynamics—use of the time-dependent Hartree approximation for a simulation of carbon-monoxide diffusion through myoglobin. J. Am. Chem. Soc. 112:9161–9175.[CrossRef]
6. Kosztin, D., S. Izrailev, and K. Schulten. 1999. Unbinding of retinoic acid from its receptor studied by steered molecular dynamics. Biophys. J. 76:188–197.
7. Ludemann, S. K., V. Lounnas, and R. C. Wade. 2000. How do substrates enter and products exit the buried active site of cytochrome P450cam? 1. Random expulsion molecular dynamics investigation of ligand access channels and mechanisms. J. Mol. Biol. 303:797–811.[CrossRef][Medline]
8. Ludemann, S. K., V. Lounnas, and R. C. Wade. 2000. How do substrates enter and products exit the buried active site of cytochrome P450cam? 2. Steered molecular dynamics and adiabatic mapping of substrate pathways. J. Mol. Biol. 303:813–830.[CrossRef][Medline]
9. Blondel, A., J. P. Renaud, S. Fischer, D. Moras, and M. Karplus. 1999. Retinoic acid receptor: a simulation analysis of retinoic acid binding and the resulting conformational changes. J. Mol. Biol. 291:101–115.[CrossRef][Medline]
10. Berman, H. M., J. Westbrook, Z. Feng, G. Gilliland, T. N. Bhat, H. Weissig, I. N. Shindyalov, and P. E. Bourne. 2000. The Protein Data Bank. Nucleic Acids Res. 28:235–242.
11. Renaud, J. P., N. Rochel, M. Ruff, V. Vivat, P. Chambon, H. Gronemeyer, and D. Moras. 1995. Crystal structure of the RAR-gamma ligand-binding domain bound to all-trans retinoic acid. Nature. 378:681–689.[CrossRef][Medline]
12. Brunger, A. T., and M. Karplus. 1988. Polar hydrogen positions in proteins: empirical energy placement and neutron diffraction comparison. Proteins. 4:148–156.[CrossRef][Medline]
13. Brooks, B. R., R. E. Bruccoleri, B. D. Olafson, D. J. States, S. Swaminathan, and M. Karplus. 1983. CHARMM: a program for macromolecular energy, minimization, and dynamics calculations. J. Comput. Chem. 4:187–217.[CrossRef]
14. MacKerell, A. D., D. Bashford, M. Bellott, R. L. Dunbrack, J. D. Evanseck, M. J. Field, S. Fischer, J. Gao, H. Guo, S. Ha, D. Joseph-McCarthy, L. Kuchnir, K. Kuczera, F. T. K. Lau, C. Mattos, S. Michnick, T. Ngo, D. T. Nguyen, B. Prodhom, W. E. Reiher, B. Roux, M. Schlenkrich, J. C. Smith, R. Stote, J. Straub, M. Watanabe, J. Wiorkiewicz-Kuczera, D. Yin, and M. Karplus. 1998. All-atom empirical potential for molecular modeling and dynamics studies of proteins. J. Phys. Chem. B. 102:3586–3616.
15. Jorgensen, W. L., J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein. 1983. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 79:926–935.[CrossRef]
16. Frisch, M. J., G. W. Trucks, H. B. Schlegel, P. M. W. Gill, B. G. Johnson, M. A. Robb, J. R. Cheeseman, T. Keith, G. A. Petersson, J. A. Montgomery, K. Raghavachari, M. A. Al-Laham, V. G. Zakrzewski, J. V. Ortiz, J. B. Foresman, J. Cioslowski, B. B. Stefanov, A. Nanayakkara, M. Challacombe, C. Y. Peng, P. Y. Ayala, W. Chen, M. W. Wong, J. L. Andres, E. S. Replogle, R. Gomperts, R. L. Martin, D. J. Fox, J. S. Binkley, D. J. Defrees, J. Baker, J. P. Stewart, M. Head-Gordon, C. Gonzalez, and J. A. Pople. 1995. Gaussian 94, Revision C.3. Gaussian, Inc., Pittsburg, PA.
17. Duan, Y., C. Wu, S. Chowdhury, M. C. Lee, G. Xiong, W. Zhang, R. Yang, P. Cieplak, R. Luo, T. Lee, J. Caldwell, J. Wang, and P. Kollman. 2003. A point-charge force field for molecular mechanics simulations of proteins based on condensed-phase quantum mechanical calculations. J. Comput. Chem. 24:1999–2012.[CrossRef][Medline]
18. Bayly, C. I., P. Cieplak, W. D. Cornell, and P. A. Kollman. 1993. A well-behaved electrostatic potential based method using charge restraints for deriving atomic charges—the RESP model. J. Phys. Chem. 97:10269–10280.[CrossRef]
19. Amber Home Page. http://amber.scripps.edu/.
20. Carlsson, P., K. Koehler, and L. Nilsson. 2005. Glucocorticoid receptor point mutation V571M facilitates coactivator and ligand binding by structural rearrangement and stabilization. Mol. Endocrinol. 19:1960–1977.
21. Kunz, S., R. Sandoval, P. Carlsson, J. Carlstedt-Duke, J. W. Bloom, and R. L. Miesfeld. 2003. Identification of a novel glucocorticoid receptor mutation in budesonide-resistant human bronchial epithelial cells. Mol. Endocrinol. 17:2566–2582.
22. Martínez, L., M. T. Sonoda, P. Webb, J. D. Baxter, M. S. Skaf, and I. Polikarpov. 2005. Molecular dynamics simulations reveal multiple pathways of ligand dissociation from thyroid hormone receptors. Biophys. J. 89:2011–2023.
23. Wagner, R. L., J. W. Apriletti, M. E. McGrath, B. L. West, J. D. Baxter, and R. J. Fletterick. 1995. A structural role for hormone in the thyroid hormone receptor. Nature. 378:690–697.[CrossRef][Medline]
24. Norman, A., M. Mizwicki, and D. Norman. 2004. Steroid-hormone rapid actions, membrane receptors and a conformational ensemble model. Nat. Rev. Drug Discov. 3:27–41.[CrossRef][Medline]
25. Zhang, Z. P., J. M. Hutcheson, H. C. Poynton, J. L. Gabriel, K. J. Soprano, and D. R. Soprano. 2003. Arginine of retinoic acid receptor beta which coordinates with the carboxyl group of retinoic acid functions independent of the amino acid residues responsible for retinoic acid receptor subtype ligand specificity. Arch. Biochem. Biophys. 409:375–384.[CrossRef][Medline]
26. Zhao, C., A. Koide, J. Abrams, S. Deighton-Collins, A. Martinez, J. A. Schwartz, S. Koide, and D. F. Skafar. 2003. Mutation of Leu-536 in human estrogen receptor-alpha alters the coupling between ligand binding, transcription activation, and receptor conformation. J. Biol. Chem. 278:27278–27286.
27. Roux, S., B. Terouanne, P. Balaguer, N. Jausons-Loffreda, M. Pons, P. Chambon, H. Gronemeyer, and J. C. Nicolas. 1996. Mutation of isoleucine 747 by a threonine alters the ligand responsiveness of the human glucocorticoid receptor. Mol. Endocrinol. 10:1214–1226.[Abstract]
28. Vottero, A., T. Kino, H. Combe, P. Lecomte, and G. P. Chrousos. 2002. A novel, C-terminal dominant negative mutation of the GR causes familial glucocorticoid resistance through abnormal interactions with p160 steroid receptor coactivators. J. Clin. Endocrinol. Metab. 87:2658–2667.
29. Langley, E., J. A. Kemppainen, and E. M. Wilson. 1998. Intermolecular NH2-/carboxyl-terminal interactions in androgen receptor dimerization revealed by mutations that cause androgen insensitivity. J. Biol. Chem. 273:92–101.
30. Collingwood, T. N., R. Wagner, C. H. Matthews, R. J. Clifton-Bligh, M. Gurnell, O. Rajanayagam, M. Agostini, R. J. Fletterick, P. Beck-Peccoz, W. Reinhardt, G. Binder, M. B. Ranke, A. Hermus, R. D. Hesch, J. Lazarus, P. Newrick, V. Parfitt, P. Raggatt, F. de Zegher, and V. K. Chatterjee. 1998. A role for helix 3 of the TRbeta ligand-binding domain in coactivator recruitment identified by characterization of a third cluster of mutations in resistance to thyroid hormone. EMBO J. 17:4760–4770.[CrossRef][Medline]
31. Gee, A. C., and J. A. Katzenellenbogen. 2001. Probing conformational changes in the estrogen receptor: evidence for a partially unfolded intermediate facilitating ligand binding and release. Mol. Endocrinol. 15:421–428.
32. Brzozowski, A. M., A. C. W. Pike, Z. Dauter, R. E. Hubbard, T. Bonn, O. Engstrom, L. Ohman, G. L. Greene, J. A. Gustafsson, and M. Carlquist. 1997. Molecular basis of agonism and antagonism in the oestrogen receptor. Nature. 389:753–758.[CrossRef][Medline]
33. Freire, E. 1998. Statistical thermodynamic linkage between conformational and binding equilibria. Adv. Protein Chem. 51:255–279.[Medline]
34. Kumar, S., B. Y. Ma, C. J. Tsai, N. Sinha, and R. Nussinov. 2000. Folding and binding cascades: dynamic landscapes and population shifts. Protein Sci. 9:10–19.[Abstract]
35. Teague, S. J. 2003. Implications of protein flexibility for drug discovery. Nat. Rev. Drug Discov. 2:527–541.[CrossRef][Medline]
36. Sonoda, M. T., N. H. Moreira, L. Martinez, F. W. Favero, S. M. Vechi, L. R. Martins, and M. S. Skaf. 2004. A review on the dynamics of water. Braz. J. Phys. 34:3–16.
This article has been cited by other articles:
 |
|
 |
|
  M. T. Sonoda, L. Martinez, P. Webb, M. S. Skaf, and I. Polikarpov Ligand Dissociation from Estrogen Receptor Is Mediated by Receptor Dimerization: Evidence from Molecular Dynamics Simulations Mol. Endocrinol., July 1, 2008; 22(7): 1565 - 1578. [Abstract] [Full Text] [PDF]
|
|
| |||||||||||||||||||||||||||||
