This is an open access article published under an ACS AuthorChoice License, which permits copying and redistribution of the article or any adaptations for non-commercial purposes.
Article pubs.acs.org/JPCB
Unraveling Binding Interactions between Human RANKL and Its Decoy Receptor Osteoprotegerin Aravinda Munasinghe,† Ping Lin,† and Coray M. Colina*,†,‡ †
Department of Chemistry, University of Florida, Gainesville, Florida 32611, United States Department of Materials Science and Engineering, University of Florida, Gainesville, Florida 32611, United States
‡
S Supporting Information *
ABSTRACT: Recent studies have revealed the importance and the active contribution of the RANKL/OPG/RANK pathway in many bone diseases including different forms of common osteoporosis. In this study, we present an extensive atomistic molecular dynamic study of the OPG/RANKL system. Within the molecular models, we varied the number of OPG molecules bound to the RANKL trimer and carried out a study to determine how the binding affinity of the OPG/RANKL system changes as a function of OPG concentration. The molecular mechanics Poisson−Boltzmann surface area method was used to analyze binding free energies. It is shown that the binding affinity decreases with increasing numbers of OPG molecules. Additionally, conformational changes of RANKL, interactions between the N-terminus outlier module of OPG with RANKL, and residues that play an important role in the binding of OPG to RANKL trimer were investigated. A probable cause for unfavorable binding for a third OPG molecule was found. Along with the currently available experimental studies, this computational study will be valuable for the comprehensive understanding of OPG/ RANKL at the atomistic level.
■
INTRODUCTION New frontiers in the field of bone remodeling science and pharmaceutical research have opened up due to the discovery of important cytokines for osteoclast biology, such as the receptor activator of nuclear factor (NF)-κB (RANK), the receptor activator of nuclear factor (NF)-κB ligand (RANKL), and decoy receptor osteoprotegerin (OPG). Bone remodeling is an important process in the homeostasis of phosphorus and calcium. The remodeling process is maintained by bone matrix formation mediated through osteoblasts and the bone resorption process is facilitated by osteoclasts.1 The equilibrium between osteoclasts and osteoblasts is necessary as it allows the skeletal structure to adapt according to mechanical limitations. However, the equilibrium between osteoblasts and osteoclasts is dependent on many factors, such as hormones, growth factors, and cytokines.2 Any activity which results in disrupting this equilibrium may result in loss of bone density which is mainly associated with development in osteopenic disorders, including postmenopausal osteoporosis, rheumatoid arthritis, and lytic bone metastases.3,4 Also, increased bone mineral density may develop in hepatitis C-associated osteosclerosis.5 Osteoclast cells are multinucleated cells which originate from hematopoiesis, while osteoblast cells originate from bone marrow stem cells. Osteoclast cells mainly differ from their monocyte, myeloid, and macrophage precursors because they have a unique requirement in activating their cell surface receptor activator. RANK is activated by its ligand RANKL.6−8 Lacey et al.9 found RANKL to be both an osteoclast differentiation factor and an activation factor. Their work © XXXX American Chemical Society
shows that the binding of RANKL to its hematopoietic descended cell generates genes to epitomize osteoclast development and also that direct administration of RANKL in vitro activates the osteoclast process. RANKL is recognized as a type II transmembrane protein because it has its carboxy terminus outside the cell.10 Its monomeric molecular weight is approximately 35 kDa, and it contains a short intracellular N-terminal tail with an extracellular region composed of a connecting stalk and a domain for receptor binding. The functional form of RANKL is found to exist as a trimer with a 3-fold axis of symmetry and binds specifically to the RANK receptor to promote the bone resorption process.9,11 The outer surface of the trimerized RANKL has three binding clefts distributed evenly across the cytokine. In the process of maintaining the bone formation and resorption equilibrium, OPG plays a major role as the decoy receptor to RANKL. This decoy receptor, OPG, prevents RANK from interacting with RANKL by binding itself to RANKL. It has been shown, that the mass of normal bones or disease state bones is mainly governed by the RANKL/OPG ratio in the system.12 Hence, manipulation of the OPG/ RANKL ratio within the system has become the basis for developing better therapeutics. OPG lacks a transmembrane region as it is generated as a soluble protein.13 OPG is an approximately 55 kDa monomer Received: July 7, 2017 Revised: September 6, 2017
A
DOI: 10.1021/acs.jpcb.7b06687 J. Phys. Chem. B XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry B
■
that homodimerizes prior to binding with RANKL.14 OPG consists of cysteine rich domains (CRD) in the extracellular region bearing similarities of most members of the tumor necrosis factor receptor superfamily. The extracellular region of OPG consists of four CRDs, and each domain contains topologically distinct modules. Specifically, the CRD1 domain contains the outlier and B2 modules, the CRD2 domain comprises the A1 and B2 modules, the CRD3 domain includes the A1 and B1 modules, and finally, the CRD4 contains the A1 and B1 modules. This classification of modules was first introduced by Naismith and Sprang.15 In the C-terminal region of the molecule, a dimerization cassette consists of two death domains, a basic motif, and a cysteine residue for dimer formation. OPG, RANKL, and RANK molecules are significant in remodeling and maintaining a healthy bone structure. Any disruption to the equilibrium between these molecules will have negative effects on the bone matrix. Therefore, due to the significance of these three molecules, extensive experimental studies have been conducted to explore the underlying chemistry of these systems, to understand their functionality and to gain structural information.16−21 Liu et al.22 have demonstrated the applicability of human RANKL as an immunogen with its rodent counterpart receptor RANK to induce an anticytokine antigen. The authors emphasize two factors that are important in creating an interspecies cytokine immunogen. First, a highly similar structural configuration between the immunogen and the interspecies cytokine needs to be achieved. Second, the immunogen should not invoke any biological response. Hence, they emphasize the importance in identifying a related species which can bind to the respective receptor and the importance of knowing the major residues within the binding cleft of that protein which can be mutated to incur an immune response. Recently, the work done by Tucker et al.23 has shown promising insight toward the applicability of bioconjugates as an effective drug delivery method, along with the prolonged blood circulation time of OPG derivatives. In their work, they have shown that among the three types of molecules they had synthesized, OPG conjugated with a loosely branched polymer has higher activity toward inhibition of osteoclast compared to wild type OPG while being nontoxic. Uncovering associated mechanistic patterns underlying the allosteric effect of branched polymers will assist not only in understanding OPG function, but also in how the polymeric materials affect biomolecular activity. Hur et al.24 recently reported the synthesis of a peptide inhibitor with enhanced RANKL affinity for osteoclastogenesis. In their work, they present experimental data along with a computational molecular dynamic and a docking study. They found that the peptide which shows highest binding affinity toward RANKL forms salt bridges with negatively charged residues in RANKL, which proves crucial for stable binding. In this work, we present an atomistic study on the OPG/ RANKL system to understand the dynamics of the system, as well as a structural analysis on RANKL with the binding of OPG. We hope this study will be helpful for current pharmaceutical research in the process of developing novel therapeutics, as well as in gaining a greater understanding of the RANKL interaction mechanism with OPG.
Article
METHODOLOGY Atomistic MD simulations were performed with the CUDAenabled Amber 14 software package25 along with the amber force field ff14SB.26 The starting structure of the protein was taken from the crystallographic data at Protein Data Bank (PDB), entry 3URF.13 The initial asymmetric crystal structure was processed by removing the N-acetyl-D-glucosamine (NAG). The OPG/RANKL complexes were constructed according to molecular symmetry using the chimera software package27 with the stoichiometric ratios of OPG:RANKL as 1:0 (free OPG), 0:3 (free RANKL trimer), 1:3 (1OPG:3RANKL), 2:3 (2OPG:3RANKL), and 3:3 (3OPG:3RANKL), as shown in Figure 1. Amber’s reduce code28 was used to add hydrogens to
Figure 1. Top view of each complex (a) 3R, (b) 1O3R, (c) 2O3R, and (d) 3O3R. Blue colored molecules represent RANKL and red colored molecules represent OPG.
each complex, and pdb 4amb code was used to define disulfide bonds. The tleap program was used to add an appropriate amount of counterions (sodium or chloride) in order to neutralize the charge of each complex. Protein complexes were solvated using the TIP3P water model29 with a truncated octahedral simulation box. The box size was defined with an initial buffering distance of 10 Å between the protein and each edge of the box. Due to the shape and size of the 3OPG:3RANKL complex, in order to avoid close contacts with periodic images at vertexes, a cubic box was used. Energy minimizations were performed using the sander module in two stages with 5000 steps of the steepest descent method and another 5000 steps of conjugate gradient. During the minimization, a nonbonded cutoff value of 10 Å was employed. In the first stage, the protein was restrained in order to preserve the crystallographic structure, while the counterions and water molecules were optimized. In the second stage, the whole system was optimized without any restraints. Followed by the minimization, each complex was heated by 50 K increments from 0 to 300 K with a 20 ps interval under the NVT ensemble. Each system was further equilibrated for 500 ps under the same conditions. The complexes were restrained during the heating process and the Langevin thermostat was employed with a 1.0 ps−1 collision frequency. The SHAKE algorithm30 was used during the MD simulation to constrain bonds involving hydrogen, allowing a 2.0 fs time step. Long range interactions beyond the nonbonded cutoff B
DOI: 10.1021/acs.jpcb.7b06687 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B were incorporated via the particle mesh Ewald (PME) summation.31 The GPU accelerated PMEMD program25 was used to obtain 70 ns of production simulation runs under the NPT ensemble at 300 K and 1 bar. Atomic coordinates and system properties were recorded every 1 ps. Trajectories were used to analyze root mean square deviation (RMSD), contact count, radii of gyration (Rg), and atomic fluctuations using CPPTRAJ32 in the AmberTools 15 software suite. Binding free energy calculations and alanine scanning mutagenesis (ASM) studies were carried out using the MMPBSA code33 which is distributed along AmberTools 15. All MMPBSA calculations were carried out using only the trajectory of the complex, in which no separate trajectories were used for the free ligands and receptors. However, for structural analysis purposes, MD simulations of free ligands and receptors were carried out. For calculations, the internal dielectric constant was set to 2.0, as suggested by Hou et al. for a moderate solute, and the external dielectric constant was kept at 80.0.34 They have shown that the solute dielectric constant does affect MMPBSA predictions, as it correlates with the characteristics of the binding interface.
Figure 3. (Top) RMSD values of protein backbone for two OPGs (red, green) and three RANKLs (blue, violet, orange) in 2O3R complex. (Bottom) RMSD values of CRD1 (green), CRD2 (red), CRD3 (blue), and CRD4 (orange) domains in one of the OPGs in 2O3R complex.
toward the ability of OPG to preserve its structural integrity and allow interdomain flexibility in between CRD domains. This is consistent with experimental observations on OPG being a very flexible molecule and, therefore, a good competitor for RANK with a higher binding affinity toward RANKL.13 These differences in RMSD values between OPG as a molecule and individually as sub domains were observed in all other OPG complexes, as shown in Figure S1. Moreover, the flexibility between CRD domains allows OPG to take different conformations. It is found that only the CRD2 and CRD3 domains efficiently contribute toward the binding of OPG to RANKL. This can be attributed to the comparatively high RMSD values for CRD1 and CRD4 domains when compared against CRD2 and CRD3 domains. Within the observed dynamic motions in OPG, one of the main findings is the interaction between RANKL and the outlier module7 (Pro5−Gln17) in the OPG CRD1 domain. Throughout each simulation, in each complex, it was observed that the outlier module tends to form a close contact with the nearby RANKL protein. This motion was captured by counting average contacts between one OPG and the associated RANKL, within 2.5 ns time spans. A contact is defined as whenever the distance of an atom in the outlier module is within 6 Å of the nearby RANKL. According to Figures 4 and S2, it can be found that the number of contacts between the
■
RESULTS AND DISCUSSION Even though the OPG/RANKL system is important in the bone remodeling process, dynamics in the solution are difficult to accurately map due to the high flexibility in OPG. In this study, the OPG/RANKL crystal structure was used as the starting structure even though it does not contain the full sequence of the OPG. The crystal structure only contains the RANKL monomer (Ala162−Asp317), and one human OPGCRD monomer (Pro5−Ser165) with three additional His residues and a glycan.13 After processing the asymmetric unit, each complex was constructed by rotating around the 3-fold rotational axis. Dynamic Analysis of OPG/RANKL Complexes. Root mean square deviation (RMSD) of the protein backbone was calculated for each complex as shown in Figure 2. It can be seen
Figure 2. Backbone RMSD of protein backbone values for 3R (green), 1O3R (red), 2O3R (blue), and 3O3R (orange).
that, with the addition of OPG to each complex, the RMSD of the complex increases when compared against the free RANKL trimer complex. In addition, the RMSD values of complexes are fluctuating around 2−5 Å. This is due to the segmental motion between CRDs in OPG, not due to the distortions within domains, as discussed below. All graphs were plotted taking data points at every 100 ps while the rest of the data is not shown for clarity. Contrary to the observed large protein backbone RMSD values in OPG, CRD subdomains have much lower values, as shown in Figure 3. This observation gives further insight
Figure 4. Each graph represents the number of average contacts made between OPG in 3O3R complex, during 2.5 ns time spans. The number of pairs within a distance of 6 Å is shown by a bar in the graph. C
DOI: 10.1021/acs.jpcb.7b06687 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B OPG outlier module and RANKL varies with time. An increment in the number of contacts signifies that the outlier module has moved toward the vicinity of RANKL, thus increasing the number of contacts. Throughout each simulation, it was observed that the outlier module of OPGs attempts to make contact with RANKL. However, no persistent contact is found during the simulation. For the outlier module to interact with the RANKL, it has to overcome the rigidity of the CRD1 domain, due to the presence of two disulfide bonds within the B2 module. However, it is uncertain how important the contribution from this interaction is beyond the 70 ns time frame studied in this work. Nevertheless, this observation may assist in any attempt to discover the mechanism behind increased binding affinity of OPG with the introduction of a conjugated polymer.23 B-factors were also calculated for the 3O3R complex using the respective 70 ns MD simulation and are presented against the experimental results in Figure 5. B-factors obtained from
Figure 6. Variation in the radius of gyration of RANKL trimer in 3O3R (orange), 2O3R (blue), 1O3R (red), and free RANKL trimer (green) with respect to time.
trimer with the addition of each OPG. A small but clear separation in Rg for the RANKL trimer in 2O3R or 3O3R can be observed when compared against the Rg of the RANKL trimer in free and the 1O3R complex. These changes in Rg for the RANKL trimer upon binding of two or more OPGs guided further investigation. In the next section, the strength of the binding of each additional OPG and their impact on the structure of RANKL trimer is presented. Binding Free Energy Study and Structural Analysis of OPG/RANKL Complex. The binding affinity between OPG and RANKL trimer is critical to their biological functions. It is known that OPG forms homodimers in its active form. According to studies carried out by Schneeweis et al.14 a fulllength, homodimeric OPG binds to RANKL trimer with a binding affinity 3 orders of magnitude higher than a truncate of osteoprotegerin without the dimerization domains. However, the crystal structure of truncated OPG and RANKL13 found a stoichiometry of 3:1 between monomeric OPG and RANKL trimer. Since it is possible for crystal packing to impose restraints in the formation of protein complexes, we carried out binding free energy calculations to access the binding affinity of each additional OPG. Binding free energy calculations were carried out using MMPBSA method without taking into account the entropic contribution (see methodology section). It was shown in several studies that adding entropic contribution using normalmode analysis does not improve the result.33,34 The binding free energy for each system was calculated using 200 frames within a 20 ns simulation window. Due to the fluctuations associated with Gibbs free energy calculations, averages were calculated in 5 ns blocks, and standard deviation among blocks is presented as the error associated with each measurement. The individual binding free energy of each OPG in each complex was calculated and reported in Figure 7. Binding free energy values can be found in Table S1. During the binding free energy calculations, one of the OPG molecules in a chosen complex was considered as the ligand and rest of the complex was considered as the receptor; this was done for all OPGs. For example, in the case of calculating binding free energies in the 3O3R complex, one selected OPG was considered as the ligand and the rest of the molecule was assigned as the receptor. The process was replicated for all three OPGs in the 3O3R complex. According to the calculated binding free energies, it was found that one of the three OPGs in the 3O3R complex has a lower binding affinity toward RANKL trimer. This corresponds to the situation where the binding of the third OPG to the 2O3R complex appears unfavorable. However, it is important to keep
Figure 5. B-factor vs residue number in the 3O3R complex. B-factor from the MD simulation is shown with the dark gray line and light gray line shows the B-factor from crystallographic data. Shaded regions are the C−D and D−E loops in each RANKL. Note that residue numbers are accorded by the PDB.
the MD simulation correlate well with the experimental Bfactors. However, as the MD simulations were carried out in the presence of a solvent at 300 K, flexibility in some of the regions is large and showed higher B-factor value compared to solid state crystallographic data. A key observation in Figure 5 is the region shaded in light gray color. Shaded regions correspond to fluctuations in the C−D (His225−Glu234) and D−E (Ser246−Pro250) loops. Nelson et al.7 earlier reported that the conformations of RANKL in the RANK/RANKL complex are surprisingly similar to the cytokine itself. However, as it was observed in those studies, D−E and C−D loops also rearrange themselves to accommodate the binding of OPG. The mobility of the D−E and C−D loops are observed by the relatively high B-factor values in RANKL in the 3O3R complex. This can be seen from the shaded region in Figure 5. Among these two loops, the D−E loop is positioned to form contacts with the first half of the CRD2 domain in OPG. The radius of gyration, Rg, calculation of the RANKL trimer with respect to the simulation time in each complex provides hints toward the fluctuations that may take place in the RANKL trimer as shown in Figure 6. Small fluctuations suggest that conformational changes might take place within the RANKL D
DOI: 10.1021/acs.jpcb.7b06687 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B
binding site as well as the interface residue information obtained from the PDBePisa web interface tool.36 Based on the initial calculation, residues with ΔΔGBind < −1.0 were filtered out (62 residues). Then, the alanine scanning mutagenesis study was carried out by only mutating those selected residues in all other binding clefts. Out of the 62 residues which were analyzed with ASM, residues which showed significant contribution to the binding of OPG to RANKL (residues with ΔΔGBind < −2.0) in 1O3R and 2O3R complexes were selected (15 residues). ASM data for the residues in OPG are given in Table 1 and ASM data for residues in each RANKL are given in Table 2. According to the crystallography study,13 key residues in RANKL that interact with RANK are Arg223, Tyr241, and Lys257. These residues were observed as key residues in OPG/ RANKL interactions during this study as well. However, in addition to these three residues, we found that residues such as Gln237, His253, and Glu300 within the binding cleft of RANKL are also important. Glu300 was also found to be an important residue in the binding process of rodent RANK to human RANKL.22 In addition to investigating the dynamics of the OPG/ RANKL complex and the corresponding affinity study of OPG in RANKL trimer with the binding of OPG, our purpose was to explore possibilities that might open additional avenues in anticytokine therapeutic antibody synthesis. For instance, mutated human RANKL was used as an interspecies cytokine to induce anti-RANKL immune responses in OVX rats.22 This work emphasized the applicability of a RANKL mutant for efficient immunotherapy of osteoporosis. For example, they stated that the structure of the mutated cytokine-induced antibody might not be the appropriate antibody to neutralize the original cytokine. Hence, with the identification of additional important residues within the binding sites of OPG and RANKL, this study may assist in synthesizing mutated cytokines in order to obtain effective antigens against RANKL and OPG. From the ASM analysis, we found residues in OPG which took part in the binding process of OPG to RANKL (see Table 1). From the ASM study, we found that residues in “loop 90” have higher contributions to the binding process of OPG to RANKL when compared with residues in the “loop 50”. A hotspot distribution was found in the OPG/RANKL complex from the ASM data. According to the residue distribution, we can observe two main binding sites in the active site of RANKL trimer, in accordance with experimental studies.13 Additionally, it has been found experimentally that the 90s loop has a higher contribution to binding of OPG toward RANKL. Similarly, in this computational study, residues with a higher contribution toward binding of OPG (dark blue
Figure 7. Binding free energy of each OPG in each complex.
in mind that the calculated binding free energy value of 4.1 ± 1.5 kcal mol−1 does not suggest that the third OPG will never bind to the 2O3R complex as the method employed does not calculate the exact binding free energy. Also, it is important to mention that this calculation does not include the crystal packing contributions. Therefore, it cannot be used to compare the crystallographic observation of 3:1 ratio between monomeric OPG and RANKL trimer, as it is known that artificial crystal contacts formation is highly probable.35 However, this observation clearly correlates with the experimental observation made in the sedimentation analysis14 in which the authors indicate that “the complex observed has a mass and s value consistent with one RANK-L bound to 2−3 monomers of OPG”. In their work, they found that the OPG without the death domain forms a stable 2:1 complex of monomeric OPG to RANKL trimer, while forming complex multimers at a high concentration of OPG. These observations along with the results obtained in this work support the fact that monomeric OPG and RANKL trimer may favor the ratio of 2:1 in solution rather than the ratio of 3:1 which was observed in the crystallographic study. The finding of unfavorable binding of the third OPG, and the structural changes which take place within RANKL trimer with the addition of each OPG, have prompted further investigation of the chemistry in the binding site of each complex. An alanine scanning mutagenesis (ASM) study has been carried out in this work to explore the contribution of each residue in the binding cleft of OPG and RANKL. Gibbs free energies were calculated in the form of ΔΔGBind as given by eq 1. Wild Mutate ΔΔG Bind = ΔG Bind − ΔG Bind
(1)
Initially, a brute force search was employed by carrying out mutational alanine scanning on all the residues near a selected OPG and RANKL binding cleft. The selection of residues (139 residues in total) was done based on the proximity from the Table 1. ASM Data for Residues in OPG Binding Sitea 1O3R
a
2O3R
3O3R
mutated residue
ΔΔGBind (OPG-1)
ΔΔGBind (OPG-1)
ΔΔGBind (OPG-2)
ΔΔGBind (OPG-1)
ΔΔGBind (OPG-2)
ΔΔGBind (OPG-3)
E58A Y61A K67A I94A E95A F96A
−6.9 (3.2) −2.9 (0.2) −4.8 (0.4) −4.3 (0.5) −13.2 (0.1) −7.5 (0.7)
−1.6 (0.1) −4.9 (1.3) −5.0 (2.8) −4.7 (0.6) −11.2 (1.7) −1.9 (0.4)
−2.3 (0.5) −2.7 (0.5) −4.1 (1.3) −2.4 (0.3) −12.9 (1.1) −4.2 (0.7)
−2.0 (0.9) −2.4 (0.2) −3.1 (1.7) −2.2 (0.7) −11.1 (1.5) −8.5 (0.2)
−1.8 (0.5) −2.6 (0.3) −6.2 (1.7) −3.0 (0.2) −12.7 (0.5) −6.2 (0.3)
−2.8 (0.8) −2.1 (0.7) 3.0 (1.2) −4.3 (0.3) −17.5 (2.1) −3.8 (0.3)
All units are in kcal mol−1. E
DOI: 10.1021/acs.jpcb.7b06687 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B Table 2. ASM Data for Residues in RANKL Binding Clefta 1O3R
a
2O3R
3O3R
mutated residue
ΔΔGBind (OPG-1)
ΔΔGBind (OPG-1)
ΔΔGBind (OPG-2)
R191A R223A Q237A Y241A H253A K257A K282A R284A D300A
−2.0 (0.2) −11.6 (0.4) −2.3 (0.5) −3.1 (0.3) −3.7 (0.5) −8.8 (0.2) −5.9 (1.5) −6.3 (2.2) −10.4 (0.6)
−3.6 (0.6) −11.0 (2.3) −3.4 (0.5) −2.4 (0.4) −3.2 (0.5) −7.3 (1.3) −3.1 (1.3) −1.5 (0.3) −10.5 (1.0)
−2.9 (3.2) −10.7 (0.2) −2.6 (1.6) −2.5 (0.4) −1.8 (1.8) −8.0 (0.8) −3.0 (0.2) −2.1 (0.2) −10.7 (0.2)
ΔΔGBind (OPG-1) −1.6 −9.5 −1.7 −3.0 −1.7 −6.3 −2.1 −3.2 −5.8
(0.5) (1.5) (1.1) (1.0) (0.1) (1.0) (0.6) (0.6) (4.9)
ΔΔGBind (OPG-2) −1.9 (0.4) −12.2 (0.9) −2.6 (0.4) −3.4 (0.3) −2.0 (0.9) −9.0 (0.6) −2.7 (0.5) −2.6 (0.5) −6.3 (1.8)
ΔΔGBind (OPG-3) −2.1 −6.7 −2.5 −0.5 −2.6 −5.6 −2.7 −4.6 −0.5
(0.7) (0.3) (1.7) (0.5) (0.8) (2.0) (0.6) (3.0) (0.7)
All units are in kcal mol−1.
residues in Figure 8) to RANKL were found to be largely localized within the 90s loop.
Figure 8. (Left) Residue distribution in RANKL binding cleft. (Right) Residue distribution in OPG. Dark blue residues correspond to residues with a higher contribution for binding free energy, and red colored residues correspond to residues with lower affinity. Figures were drawn using the UCSF-Chimera program.27
Finally, a structural analysis was carried out on the 3O3R complex to find a probable cause for unfavorable binding of the third OPG. As a starting point for this analysis, the ASM data was also used. Even though many residues seem to have relatively similar contributions of ΔΔGBind values when compared to the rest of the key residues in the binding cleft, Asp300 in the third binding cleft showed a difference of 5 kcal mol−1 in ΔΔGBind. Further analysis showed that the distance between Asp300 in RANKL and Lys67 in OPG decreased for the third OPG, when the distance between Lys257 in the adjacent RANKL and Asp300 increased, as shown in Figure 9. According to the graph, we can clearly observe that the fraction of states with Asp300−Lys67 interaction in the OPG-3 binding cleft decreases when comparing to the fraction of states with Asp300−Lys67 interaction in the OPG-2 binding cleft. This distance change has prompted the RANKL−RANKL interaction in the third binding pocket and lowered the affinity of Asp300 toward binding of the third OPG.
Figure 9. (Top) The interactions in Asp300 from a snapshot of 3O3R at third OPG binding site. The orientation of two states of the binding site are shown via violet and green color bond representations. (Bottom) Distance variation in Asp300.
CRD1 domain and CRD4 domain was observed while CRD2 and CRD3 domains were comparatively more stable. Additionally, it was observed that the outlier module in CRD1 domain tends to make contacts with the nearby RANKL molecule. However, the interaction between the outlier module and neighboring RANKL seems not to be in a favorable state due to the rigidity of the B2 module within the CRD1 domain. Further analysis has shown that structural changes were observed with the addition of each OPG to RANKL trimer. The binding free
■
CONCLUSION Fully atomistic simulations were carried out on OPG/RANKL systems with varying stoichiometric ratios of OPG monomer: RANKL trimer as 3:1, 2:1, and 1:1. Each simulation was subjected to dynamic analysis, including analysis of the relative motion of CRD domains in OPG. Large relative motion in the F
DOI: 10.1021/acs.jpcb.7b06687 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B
RANK and the Osteoprotegerin Decoy Receptor. Structure 2012, 20, 1971−1982. (8) Boyce, B. F.; Yao, Z.; Xing, L. Osteoclasts have Multiple Roles in Bone in Addition to Bone Resorption. Crit. Rev. Eukaryotic Gene Expression 2009, 19, 171−180. (9) Lacey, D. L.; Timms, E.; Tan, H.-L.; Kelley, M. J.; Dunstan, C. R.; Burgess, T.; Elliott, R.; Colombero, A.; Elliott, G.; Scully, S.; et al. Osteoprotegerin Ligand is a Cytokine that Regulates Osteoclast Differentiation and Activation. Cell 1998, 93, 165−176. (10) Tanaka, S.; Nakamura, K.; Takahasi, N.; Suda, T. Role of RANKL in Physiological and Pathological Bone Resorption and Therapeutics Targeting the RANKL-RANK Signaling System. Immunol. Rev. 2005, 208, 30−49. (11) Yasuda, H.; Shima, N.; Nakagawa, N.; Yamaguchi, K.; Kinosaki, M.; Mochizuki, S.; Tomoyasu, A.; Yano, K.; Goto, M.; Murakami, A.; et al. Osteoclast Differentiation Factor is a Ligand for Osteoprotegerin/Osteoclastogenesis-Inhibitory Factor and is Identical to TRANCE/ RANKL. Proc. Natl. Acad. Sci. U. S. A. 1998, 95, 3597−3602. (12) Boyce, B. F.; Xing, L. The RANKL/RANK/OPG Pathway. Curr. Osteoporos. Rep. 2007, 5, 98−104. (13) Luan, X.; Lu, Q.; Jiang, Y.; Zhang, S.; Wang, Q.; Yuan, H.; Zhao, W.; Wang, J.; Wang, X. Crystal Structure of Human RANKL Complexed with its Decoy Receptor Osteoprotegerin. J. Immunol. 2012, 189, 245−252. (14) Schneeweis, L. A.; Willard, D.; Milla, M. E. Functional Dissection of Osteoprotegerin and its Interaction with Receptor Activator of NF-κB Ligand. J. Biol. Chem. 2005, 280, 41155−41164. (15) Naismith, J. H.; Sprang, S. R. Modularity in the TNF-Receptor Family. Trends Biochem. Sci. 1998, 23, 74−79. (16) Kitazawa, R.; Kitazawa, S.; Maeda, S. Promoter Structure of Mouse RANKL/TRANCE/OPGL/ODF Gene. Biochim. Biophys. Acta, Gene Struct. Expression 1999, 1445, 134−141. (17) Malyankar, U. M.; Scatena, M.; Suchland, K. L.; Yun, T. J.; Clark, E. A.; Giachelli, C. M. Osteoprotegerin is an αvβ3-Induced, NFκB-dependent Survival Factor for Endothelial Cells. J. Biol. Chem. 2000, 275, 20959−20962. (18) Mizuno, A.; Amizuka, N.; Irie, K.; Murakami, A.; Fujise, N.; Kanno, T.; Sato, Y.; Nakagawa, N.; Yasuda, H.; Mochizuki, S.; et al. Severe Osteoporosis in Mice Lacking Osteoclastogenesis Inhibitory Factor/Osteoprotegerin. Biochem. Biophys. Res. Commun. 1998, 247, 610−615. (19) Bucay, N.; Sarosi, I.; Dunstan, C. R.; Morony, S.; Tarpley, J.; Capparelli, C.; Scully, S.; Tan, H. L.; Xu, W.; Lacey, D. L.; et al. Osteoprotegerin-Deficient Mice Develop Early Onset Osteoporosis and Arterial Calcification. Genes Dev. 1998, 12, 1260−1268. (20) Rodan, G. A.; Martin, T. J. Role of Osteoblasts in Hormonal Control of Bone resorptionA Hypothesis. Calcif. Tissue Int. 1981, 33, 349−351. (21) Simonet, W. S.; Lacey, D. L.; Dunstan, C. R.; Kelley, M.; Chang, M.-S.; Lüthy, R.; Nguyen, H. Q.; Wooden, S.; Bennett, L.; Boone, T.; et al. Osteoprotegerin: A Novel Secreted Protein Involved in the Regulation of Bone Density. Cell 1997, 89, 309−319. (22) Liu, C.; Zhao, Y.; He, W.; Wang, W.; Chen, Y.; Zhang, S.; Ma, Y.; Gohda, J.; Ishida, T.; Walter, T. S.; et al. A RANKL Mutant Used as an Inter-Species Vaccine for Efficient Immunotherapy of Osteoporosis. Sci. Rep. 2015, 5, 14150. (23) Tucker, B. S.; Stewart, J. D.; Aguirre, J. I.; Holliday, L. S.; Figg, C. A.; Messer, J. G.; Sumerlin, B. S. Role of Polymer Architecture on the Activity of Polymer-Protein Conjugates for the Treatment of Accelerated Bone Loss Disorders. Biomacromolecules 2015, 16, 2374− 2381. (24) Hur, J.; Ghosh, A.; Kim, K.; Ta, H. M.; Kim, H.; Kim, N.; Hwang, H.-Y.; Kim, K. K. Design of a RANK-Mimetic Peptide Inhibitor of Osteoclastogenesis with Enhanced RANKL-Binding Affinity. Mol Cells 2016, 39, 316−321. (25) Salomon-Ferrer, R.; Götz, A. W.; Poole, D.; Le Grand, S.; Walker, R. C. Routine Microsecond Molecular Dynamics Simulations with AMBER on GPUs. 2. Explicit Solvent Particle Mesh Ewald. J. Chem. Theory Comput. 2013, 9, 3878−3888.
energy study has shown that binding of the third OPG is unfavorable compared to binding of the first and second OPG to the RANKL trimer. An ASM study has revealed several significant residues which contribute significantly toward the binding of OPG to RANKL and these findings are in good agreement with experimental data. Further, we were able to report extra key residues in OPG and RANKL which take a prominent role in binding of OPG to RANKL in addition to what previous studies have revealed. Finally, the reduced avidity of the third OPG to RANKL trimer points toward conformational changes that leads up to enhanced interaction between Asp300 and Lys257 in RANKL and weakened interaction between Asp300 in RANKL and Lys67 in OPG.
■
ASSOCIATED CONTENT
S Supporting Information *
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.7b06687. RMSD of OPGs, their subdomains, and RANKLs in 1O3R complex and 3O3R complex; average contacts OPGs made during the simulations in 1O3R complex and 2O3R complex; and binding free energy values of each OPG in each complex (PDF)
■
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected]fl.edu; Phone: +1 (352)294 3488. ORCID
Ping Lin: 0000-0002-6141-5424 Coray M. Colina: 0000-0003-2367-1352 Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS The authors acknowledge the financial assistance given by the University of Florida Preeminence Initiative and University of Florida Research Computing for providing computational support that have contributed to the research results reported in this publication.
■
REFERENCES
(1) Clarke, B. Normal Bone Anatomy and Physiology. Clin. J. Am. Soc. Nephrol. 2008, 3 (3), S131−S139. (2) Aoki, S.; Honma, M.; Kariya, Y.; Nakamichi, Y.; Ninomiya, T.; Takahashi, N.; Udagawa, N.; Suzuki, H. Function of OPG as a Traffic Regulator for RANKL is Crucial for Controlled Osteoclastogenesis. J. Bone Miner. Res. 2010, 25, 1907−1921. (3) Rodan, G. A.; Martin, T. J. Therapeutic Approaches to Bone Diseases. Science 2000, 289, 1508−1514. (4) Takayanagi, H.; Iizuka, H.; Juji, T.; Nakagawa, T.; Yamamoto, A.; Miyazaki, T.; Koshihara, Y.; Oda, H.; Nakamura, K.; Tanaka, S. Involvement of Receptor Activator of Nuclear Factor κB Ligand/ Osteoclast Differentiation Factor in Osteoclastogenesis from Synoviocytes in Rheumatoid Arthritis. Arthritis Rheum. 2000, 43, 259−269. (5) Gregson, C. L.; Hardcastle, S. A.; Cooper, C.; Tobias, J. H. Friend or Foe: High Bone Mineral Density on Routine Bone Density Scanning, a Review of Causes and Management. Rheumatology (Oxford, U. K.) 2013, 52, 968−985. (6) Leibbrandt, A.; Penninger, J. M. RANK/RANKL: Regulators of Immune Responses and Bone Physiology. Ann. N. Y. Acad. Sci. 2008, 1143, 123−150. (7) Nelson, C. A.; Warren, J. T.; Wang, M. W.-H.; Teitelbaum, S. L.; Fremont, D. H. RANKL Employs Distinct Binding Modes to Engage G
DOI: 10.1021/acs.jpcb.7b06687 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B (26) Maier, J. A.; Martinez, C.; Kasavajhala, K.; Wickstrom, L.; Hauser, K. E.; Simmerling, C. ff14SB: Improving the Accuracy of Protein Side Chain and Backbone Parameters from ff99SB. J. Chem. Theory Comput. 2015, 11, 3696−3713. (27) Pettersen, E. F.; Goddard, T. D.; Huang, C. C.; Couch, G. S.; Greenblatt, D. M.; Meng, E. C.; Ferrin, T. E. UCSF Chimera−a Visualization System for Exploratory Research and Analysis. J. Comput. Chem. 2004, 25, 1605−1612. (28) Word, J. M.; Lovell, S. C.; Richardson, J. S.; Richardson, D. C. Asparagine and Glutamine: Using Hydrogen Atom Contacts in the Choice of Side-Chain Amide Orientation. J. Mol. Biol. 1999, 285, 1735−1747. (29) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926−935. (30) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. Numerical Integration of the Cartesian Equations of Motion of a System with Constraints: Molecular Dynamics of N-Alkanes. J. Comput. Phys. 1977, 23, 327−341. (31) Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald: An N· log(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089−10093. (32) Roe, D. R.; Cheatham, T. E. PTRAJ and CPPTRAJ: Software for Processing and Analysis of Molecular Dynamics Trajectory Data. J. Chem. Theory Comput. 2013, 9, 3084−3095. (33) Genheden, S.; Ryde, U. The MM/PBSA and MM/GBSA Methods to Estimate Ligand-Binding Affinities. Expert Opin. Drug Discovery 2015, 10, 449−461. (34) Hou, T.; Wang, J.; Li, Y.; Wang, W. Assessing the Performance of the MM/PBSA and MM/GBSA Methods. 1. The Accuracy of Binding Free Energy Calculations Based on Molecular Dynamics Simulations. J. Chem. Inf. Model. 2011, 51, 69−82. (35) Levy, E. D.; Teichmann, S. A. Structural, Evolutionary, and Assembly Principles of Protein Oligomerization. In Progress in Molecular Biology and Translational Science; Elsevier Inc., 2013; Vol. 117, pp 25−51. (36) Krissinel, E.; Henrick, K. Inference of Macromolecular Assemblies from Crystalline State. J. Mol. Biol. 2007, 372, 774−797.
H
DOI: 10.1021/acs.jpcb.7b06687 J. Phys. Chem. B XXXX, XXX, XXX−XXX
