Mechanisms for allosteric activation of protease DegS by ligand binding and oligomerization as revealed from molecular dynamics simulations

A local perturbation of a protein may lead to functional changes at some distal site, a phenomenon denoted as allostery. Here, we study the allosteric control of a protease using molecular dynamics simulations. The system considered is the bacterial protein DegS which includes a protease domain activated on ligand binding to an adjacent PDZ domain. Starting from crystallographic structures of DegS homo‐trimers, we perform simulations of the ligand‐free and ‐bound state of DegS at equilibrium. Considering a single protomer only, the trimeric state was mimicked by applying restraints on the residues in contact with other protomers in the DegS trimer. In addition, the bound state was also simulated without any restraints to mimic the monomer. Our results suggest that not only ligand release but also disassembly of a DegS trimer inhibits proteolytic activity. Considering various observables for structural changes, we infer allosteric pathways from the interface with other protomers to the active site. Moreover, we study how ligand release leads to (i) catalytically relevant changes involving residues 199–201 and (ii) a transition from a stretched to a bent conformation for residues 217–219 (which prohibits proper substrate binding). Finally, based on ligand‐induced C_α shifts we identify residues in contact with other protomers in the DegS trimer that likely transduce the perturbation from ligand release from a given protomer to adjacent protomers. These residues likely play a key role in the experimentally known effect of ligand release from a protomer on the proteolytic activity of the other protomers. Proteins 2016; 84:1690–1705. © 2016 Wiley Periodicals, Inc.     

How an Inhibitor Bound to Subunit Interface Alters Triosephosphate Isomerase Dynamics

The tunnel region at triosephosphate isomerase (TIM)’s dimer interface, distant from its catalytic site, is a target site for certain benzothiazole derivatives that inhibit TIM’s catalytic activity in Trypanosoma cruzi, the parasite that causes Chagas disease. We performed multiple 100-ns molecular-dynamics (MD) simulations and elastic network modeling (ENM) on both apo and complex structures to shed light on the still unclear inhibitory mechanism of one such inhibitor, named bt10. Within the time frame of our MD simulations, we observed stabilization of aromatic clusters at the dimer interface and enhancement of intersubunit hydrogen bonds in the presence of bt10, which point to an allosteric effect rather than destabilization of the dimeric structure. The collective dynamics dictated by the topology of TIM is known to facilitate the closure of its catalytic loop over the active site that is critical for substrate entrance and product release. We incorporated the ligand’s effect on vibrational dynamics by applying mixed coarse-grained ENM to each one of 54,000 MD snapshots. Using this computationally efficient technique, we observed altered collective modes and positive shifts in eigenvalues due to the constraining effect of bt10 binding. Accordingly, we observed allosteric changes in the catalytic loop’s dynamics, flexibility, and correlations, as well as the solvent exposure of catalytic residues. A newly (to our knowledge) introduced technique that performs residue-based ENM scanning of TIM revealed the tunnel region as a key binding site that can alter global dynamics of the enzyme.     

Investigation of allosteric coupling in human β<Subscript>2</Subscript>-adrenergic receptor in the presence of intracellular loop 3

Background

This study investigates the allosteric coupling that exists between the intra- and extracellular parts of human β₂-adrenergic receptor (β₂-AR), in the presence of the intracellular loop 3 (ICL3), which is missing in all crystallographic experiments and most of the simulation studies reported so far. Our recent 1 μs long MD run has revealed a transition to the so-called very inactive state of the receptor, in which ICL3 packed under the G protein’s binding cavity and completely blocked its accessibility to G protein. Simultaneously, an outward tilt of transmembrane helix 5 (TM5) caused an expansion of the extracellular ligand-binding site. In the current study, we performed independent runs with a total duration of 4 μs to further investigate the very inactive state with packed ICL3 and the allosteric coupling event (three unrestrained runs and five runs with bond restraints at the ligand-binding site).

Results

In all three independent unrestrained runs (each 500 ns long), ICL3 preserved its initially packed/closed conformation within the studied time frame, suggesting an inhibition of the receptor’s activity. Specific bond restraints were later imposed between some key residues at the ligand-binding site, which have been experimentally determined to interact with the ligand. Restraining the binding site region to an open state facilitated ICL3 closure, whereas a relatively constrained/closed binding site hindered ICL3 packing. However, the reverse operation, i.e. opening of the packed ICL3, could not be realized by restraining the binding site region to a closed state. Thus, any attempt failed to free the ICL3 from its locked state due to the presence of persistent hydrogen bonds.

Conclusions

Overall, our simulations indicated that starting with very inactive states, the receptor stayed almost irreversibly inhibited, which in turn decreased the overall mobility of the receptor. Bond restraints which represented the geometric restrictions caused by ligands of various sizes when bound at the ligand-binding site, induced the expected conformational changes in TM5, TM6 and consequently, ICL3. Still, once ICL3 was packed, the allosteric coupling became ineffective due to strong hydrogen bonds connecting ICL3 to the core of the receptor.

     

A molecular dynamics study of the C-terminal fragment of the L7/L12 ribosomal protein

The crystallographic dimer of the C-terminal fragment (CTF) of the L7/L12 ribosomal protein has been subjected to molecular dynamics (MD) simulations. A 90 picosecond (ps) trajectory for the protein dimer, 19 water molecules and two counter ions has been calculated at constant temperature. Effects of intermolecular interactions on the structure and dynamics have been studied. The exact crystallographic symmetry is lost and the atomic fluctuations differ from one monomer to the other. The average MD structure is more stable than the X-ray one, as judged by accessible surface area and energy calculations. Crystal (non-dimeric) interactions have been simulated in another 40 ps trajectory by using harmonic restraints to represent intermolecular hydrogen bonds. The conformational changes with respect ot the X-ray structure are then virtually suppressed.The unrestrained dimer trajectory has been scanned for cooperative motions involving secondary structure elements. The intrinsic collective motions of the monomer are transmitted via intermolecular contacts to the dimer structure.The existence of a stable dimeric form of CTF, resembling the crystallographic one, has been documented. At the cost of fairly small energy expenditure the dimer has considerable conformational flexibility. This flexibility may endow the dimer with some functional potential as an energy transducer.     

Mechanism of μ-Opioid Receptor-Magnesium Interaction and Positive Allosteric Modulation

In the era of opioid abuse epidemics, there is an increased demand for understanding how opioid receptors can be allosterically modulated to guide the development of more effective and safer opioid therapies. Among the modulators of the μ-opioid (MOP) receptor, which is the pharmacological target for the majority of clinically used opioid drugs, are monovalent and divalent cations. Specifically, the monovalent sodium cation (Na⁺) has been known for decades to affect MOP receptor signaling by reducing agonist binding, whereas the divalent magnesium cation (Mg²⁺) has been shown to have the opposite effect, notwithstanding the presence of sodium chloride. Although ultra-high-resolution opioid receptor crystal structures have revealed a specific Na⁺ binding site and molecular dynamics (MD) simulation studies have supported the idea that this monovalent ion reduces agonist binding by stabilizing the receptor inactive state, the putative binding site of Mg²⁺ on the MOP receptor, as well as the molecular determinants responsible for its positive allosteric modulation of the receptor, are unknown. In this work, we carried out tens of microseconds of all-atom MD simulations to investigate the simultaneous binding of Mg²⁺ and Na⁺ cations to inactive and active crystal structures of the MOP receptor embedded in an explicit lipid-water environment and confirmed adequate sampling of Mg²⁺ ion binding with a grand canonical Monte Carlo MD method. Analyses of these simulations shed light on 1) the preferred binding sites of Mg²⁺ on the MOP receptor, 2) details of the competition between Mg²⁺ and Na⁺ cations for specific sites, 3) estimates of binding affinities, and 4) testable hypotheses of the molecular mechanism underlying the positive allosteric modulation of the MOP receptor by the Mg²⁺ cation.     

Metropolis Monte Carlo calculations of DNA structure using internal coordinates and NMR distance restraints: An alternative method for generating a high-resolution solution structure

Summary A new method, a restrained Monte Carlo (rMC) calculation, is demonstrated for generating high-resolution structures of DNA oligonucleotides in solution from interproton distance restraints and bounds derived from complete relaxation matrix analysis of two-dimensional nuclear Overhauser effect (NOE) spectral peak intensities. As in the case of restrained molecular dynamics (rMD) refinement of structures, the experimental distance restraints and bounds are incorporated as a pseudo-energy term (or penalty function) into the mathematical expression for the molecular energy. However, the use of generalized helical parameters, rather than Cartesian coordinates, to define DNA conformation increases efficiency by decreasing by an order of magnitude the number of parameters needed to describe a conformation and by simplifying the potential energy profile. The Metropolis Monte Carlo method is employed to simulate an annealing process. The rMC method was applied to experimental 2D NOE data from the octamer duplex d(GTA-TAATG)·d(CATTATAC). Using starting structures from different locations in conformational space (e.g. A-DNA and B-DNA), the rMC calculations readily converged, with a root-mean-square deviation (RMSD) of <0.3 Å between structures generated using different protocols and starting structures. Theoretical 2D NOE peak intensities were calculated for the rMC-generated structures using the complete relaxation matrix program CORMA, enabling a comparison with experimental intensities via residual indices. Simulation of the vicinal proton coupling constants was carried out for the structures generated, enabling a comparison with the experimental deoxyribose ring coupling constants, which were not utilized in the structure determination in the case of the rMC simulations. Agreement with experimental 2D NOE and scalar coupling data was good in all cases. The rMC structures are quite similar to that refined by a traditional restrained MD approach (RMSD<0.5 Å) despite the different force fields used and despite the fact that MD refinement was conducted with additional restraints imposed on the endocyclic torsion angles of deoxyriboses. The computational time required for the rMC and rMD calculations is about the same. A comparison of structural parameters is made and some limitations of both methods are discussed with regard to the average nature of the experimental restraints used in the refinement.Abbreviations MC Monte Carlo - rMC restrained Monte Carlo - MD molecular dynamics - rMD restrained molecular dynamics - DG distance geometry - EM energy minimization - 2D NOE two-dimensional nuclear Overhauser effect - DQF-COSY double-quantum-filtered correlation spectroscopy - RMSD root-mean-square deviation To whom correspondence should be addressed.     

Structure and dynamics of the full-length M1 muscarinic acetylcholine receptor studied by molecular dynamics simulations

The three-dimensional structure of full-length structure of the M₁ muscarinic receptor was obtained through the fragmental homology modeling procedure. A 10-ns molecular dynamics (MD) simulation of the protein imbedded in a lipid slab and surrounded by water molecules was further used to relax the model. It was found that the homology model corresponded to the conformation in the ground state, since no significant motions of the backbone of transmembrane domains were observed. Furthermore, the reliability of the model was validated by analyzing key inter-helical contacts, sidechain-sidechain interactions, the formation of stable aromatic microdomains (clusters) and the docking of acetylcholine to its binding site. Moreover, a few conserved interactions observed in the X-ray structure of rhodopsin, such as inter-helical sidechain-sidechain hydrogen bonds were accurately reproduced in the MD simulation. The coupling of ACh to its binding site was found to be dominated by π-cation and salt bridge interactions, while its conformational space was restrained through van der Waals and hydrogen bond interactions. In general, such features were in very good agreement with the available experimental as well as with theoretical data. Considering the above, the structural information obtained in this study can be used a starting point to investigate the activation mechanism of the receptor and the ability to develop selective agonists and allosteric modulators which could be used for the treatment of Alzheimer’s disease.     

Negative allosteric modulators of cannabinoid receptor 2: protein modeling,binding site identification and molecular dynamics simulations in the presence of an orthosteric agonist

Abstract

Selective activation of the cannabinoid receptor subtype 2 (CB2) shows promise for treating pain, inflammation, multiple sclerosis, cancer, ischemic/reperfusion injury and osteoporosis. Target selectivity and off-target side effects are two major limiting factors for orthosteric ligands, and therefore, the search for allosteric modulators (AMs) is a widely used drug discovery approach. To date, only a limited number of negative CB2 AMs have been identified, possessing only micromolar activity at best, and the CB2 receptor’s allosteric site(s) are not well characterized. Herein, we used computational approaches including receptor modeling, site mapping, docking, molecular dynamics (MD) simulations and binding free energy calculations to predict, characterize and validate allosteric sites within the complex of the CB2 receptor with bound orthosteric agonist CP55,940. After docking of known negative CB2 allosteric modulators (NAMs), dihydro-gambogic acid (DHGA) and trans-β-caryophyllene (TBC) (note that TBC also shows agonist activity), at the predicted allosteric sites, the best total complex with CB2, CP55,940 and NAM was embedded into a hydrated lipid bilayer and subjected to a 200 ns MD simulation. The presence of an AM affected the CB2–CP55,940 complex, altering the relative positioning of the toggle switch residues and promoting a strong π–π interaction between Phe117^3.36 and Trp258^6.48. Binding of either TBC or DHGA to a putative allosteric pocket directly adjacent to the orthosteric ligand reduced the binding free energy of CP55,940, which is consistent with the expected effect of a negative AM. The identified allosteric sites present immense scope for the discovery of novel classes of CB2 AMs.     

Molecular dynamics re-refinement of two different small RNA loop structures using the original NMR data suggest a common structure

Restrained molecular dynamics simulations are a robust, though perhaps underused, tool for the end-stage refinement of biomolecular structures. We demonstrate their utility-using modern simulation protocols, optimized force fields, and inclusion of explicit solvent and mobile counterions-by re-investigating the solution structures of two RNA hairpins that had previously been refined using conventional techniques. The structures, both domain 5 group II intron ribozymes from yeast ai5γ and Pylaiella littoralis, share a nearly identical primary sequence yet the published 3D structures appear quite different. Relatively long restrained MD simulations using the original NMR restraint data identified the presence of a small set of violated distance restraints in one structure and a possibly incorrect trapped bulge nucleotide conformation in the other structure. The removal of problematic distance restraints and the addition of a heating step yielded representative ensembles with very similar 3D structures and much lower pairwise RMSD values. Analysis of ion density during the restrained simulations helped to explain chemical shift perturbation data published previously. These results suggest that restrained MD simulations, with proper caution, can be used to "update" older structures or aid in the refinement of new structures that lack sufficient experimental data to produce a high quality result. Notable cautions include the need for sufficient sampling, awareness of potential force field bias (such as small angle deviations with the current AMBER force fields), and a proper balance between the various restraint weights.     

A molecular dynamics investigation of mono and dimeric states of the outer membrane enzyme OMPLA

OMPLA is a phospholipase found in the outer membranes of many Gram-negative bacteria. Enzyme activation requires calcium-induced dimerisation plus bilayer perturbation. As the conformation of OMPLA in the different crystal forms (monomer versus dimer; with/without bound Ca(2+)) is remarkably similar we have used multi-nanosecond molecular dynamics (MD) simulations to probe possible differences in conformational dynamics that may be related to enzyme activation. Simulations of calcium-free monomeric OMPLA, of the Ca(2+)-bound dimer, and of the Ca(2+)-bound dimer with a substrate analogue covalently linked to the active site serine have been performed, all with the protein embedded in a phospholipid (POPC) bilayer. All simulations were stable, but differences in the dynamic behaviour of the protein between the various states were observed. In particular, the stability of the active site and the hydrophobic substrate-binding cleft varied. Dimeric OMPLA is less flexible than monomeric OMPLA, especially around the active site. In the absence of bound substrate analogue, the hydrophobic substrate-binding cleft of dimeric OMPLA collapses. A model is proposed whereby the increased stability of the active site in dimeric OMPLA is a consequence of the local ordering of water around the nearby calcium ion. The observed collapse of the substrate-binding cleft may explain the experimentally observed occurrence of multiple dimer conformations of OMPLA, one of which is fully active while the other shows significantly reduced activity.     

Change in allosteric network affects binding affinities of PDZ domains: analysis through perturbation response scanning

The allosteric mechanism plays a key role in cellular functions of several PDZ domain proteins (PDZs) and is directly linked to pharmaceutical applications; however, it is a challenge to elaborate the nature and extent of these allosteric interactions. One solution to this problem is to explore the dynamics of PDZs, which may provide insights about how intramolecular communication occurs within a single domain. Here, we develop an advancement of perturbation response scanning (PRS) that couples elastic network models with linear response theory (LRT) to predict key residues in allosteric transitions of the two most studied PDZs (PSD-95 PDZ3 domain and hPTP1E PDZ2 domain). With PRS, we first identify the residues that give the highest mean square fluctuation response upon perturbing the binding sites. Strikingly, we observe that the residues with the highest mean square fluctuation response agree with experimentally determined residues involved in allosteric transitions. Second, we construct the allosteric pathways by linking the residues giving the same directional response upon perturbation of the binding sites. The predicted intramolecular communication pathways reveal that PSD-95 and hPTP1E have different pathways through the dynamic coupling of different residue pairs. Moreover, our analysis provides a molecular understanding of experimentally observed hidden allostery of PSD-95. We show that removing the distal third alpha helix from the binding site alters the allosteric pathway and decreases the binding affinity. Overall, these results indicate that (i) dynamics plays a key role in allosteric regulations of PDZs, (ii) the local changes in the residue interactions can lead to significant changes in the dynamics of allosteric regulations, and (iii) this might be the mechanism that each PDZ uses to tailor their binding specificities regulation.     

Induced Effects of Sodium Ions on Dopaminergic G-Protein Coupled Receptors

G-protein coupled receptors, the largest family of proteins in the human genome, are involved in many complex signal transduction pathways, typically activated by orthosteric ligand binding and subject to allosteric modulation. Dopaminergic receptors, belonging to the class A family of G-protein coupled receptors, are known to be modulated by sodium ions from an allosteric binding site, although the details of sodium effects on the receptor have not yet been described. In an effort to understand these effects, we performed microsecond scale all-atom molecular dynamics simulations on the dopaminergic D₂ receptor, finding that sodium ions enter the receptor from the extracellular side and bind at a deep allosteric site (Asp2.50). Remarkably, the presence of a sodium ion at this allosteric site induces a conformational change of the rotamer toggle switch Trp6.48 which locks in a conformation identical to the one found in the partially inactive state of the crystallized human β₂ adrenergic receptor. This study provides detailed quantitative information about binding of sodium ions in the D₂ receptor and reports a possibly important sodium-induced conformational change for modulation of D₂ receptor function.     

The conformational feasibility for the formation of reaching dimer in ASV and HIV integrase: a molecular dynamics study

Retroviral integrases are reported to form alternate dimer assemblies like the core–core dimer and reaching dimer. The core–core dimer is stabilized predominantly by an extensive interface between two catalytic core domains. The reaching dimer is stabilized by N-terminal domains that reach to form intermolecular interfaces with the other subunit’s core and C-terminal domains (CTD), as well as CTD–CTD interactions. In this study, molecular dynamics (MD), Brownian dynamics (BD) simulations, and free energy analyses, were performed to elucidate determinants for the stability of the reaching dimer forms of full-length Avian Sarcoma Virus (ASV) and Human Immunodeficiency Virus (HIV) IN, and to examine the role of the C-tails (the last ~16–18 residues at the C-termini) in their structural dynamics. The dynamics of an HIV reaching dimer derived from small angle X-ray scattering and protein crosslinking data, was compared with the dynamics of a core–core dimer model derived from combining the crystal structures of two-domain fragments. The results showed that the core domains in the ASV reaching dimer express free dynamics, whereas those in the HIV reaching dimer are highly stable. BD simulations suggest a higher rate of association for the HIV core–core dimer than the reaching dimer. The predicted stability of these dimers was therefore ranked in the following order: ASV reaching dimer < HIV reaching dimer < composite core–core dimer. Analyses of MD trajectories have suggested residues that are critical for intermolecular contacts in each reaching dimer. Tests of these predictions and insights gained from these analyses could reveal a potential pathway for the association and dissociation of full-length IN multimers.     

Insight derived from molecular dynamics simulation into dynamics and molecular motions of cuticle-degrading serine protease Ver112

Cuticle-degrading serine protease Ver112, which derived from a nematophagous fungus Lecanicillium psalliotae, has been exhibited to have high cuticle-degrading and nematicidal activities. We have performed molecular dynamics (MD) simulation based on the crystal structure of Ver112 to investigate its dynamic properties and large-scale concerted motions. The results indicate that the structural core of Ver112 shows a small fluctuation amplitude, whereas the substrate binding sites, and the regions close to and opposite the substrate binding sites experience significant conformational fluctuations. The large concerted motions obtained from essential dynamics (ED) analysis of MD trajectory can lead to open or close of the substrate binding sites, which are proposed to be linked to the functional properties of Ver112, such as substrate binding, orientation, catalytic, and release. The significant motion in the loop regions that is located opposite the binding sites are considered to play an important role in modulating the dynamics of the substrate binding sites. Furthermore, the bottom of free energy landscape (FEL) of Ver112 are rugged, which is mainly caused by the fluctuations of substrate binding regions and loops located opposite the binding site. In addition, the mechanism underlying the high flexibility and catalytic activity of Ver112 was also discussed. Our simulation study complements the biochemical and structural studies, and provides insight into the dynamics-function relationship of cuticle-degrading serine protease Ver112.     

Anesthetic interaction with ketosteroid isomerase: insights from molecular dynamics simulations

The nature and the sites of interactions between anesthetic halothane and homodimeric Delta5-3-ketosteroid isomerase (KSI) are characterized by flexible ligand docking and confirmed by 1H-15N NMR. The dynamics consequence of halothane interaction and the implication of the dynamic changes to KSI function are studied by multiple 5-ns molecular dynamics simulations in the presence and absence of halothane. Both docking and MD simulations show that halothane prefer the amphiphilic dimeric interface to the hydrophobic active site of KSI. Halothane occupancy at the dimer interface disrupted the intersubunit hydrogen bonding formed either directly through side chains of polar residues or indirectly through the mediation of the interfacial water molecules. Moreover, in the presence of halothane, the exchange rate of the bound waters with bulk water was increased. Halothane perturbation to the dimer interface affected the overall flexibility of the active site. This action is likely to contribute to the halothane-induced reduction of the KSI activity. The allosteric halothane modulation of the dynamics-function relationship of KSI without direct competition at the enzymatic active sites may be generalized to offer a unifying explanation of anesthetic action on a diverse range of multidomain neuronal proteins that are potentially relevant to clinical general anesthesia.     

Control of cation permeation through the nicotinic receptor channel

We used molecular dynamics (MD) simulations to explore the transport of single cations through the channel of the muscle nicotinic acetylcholine receptor (nAChR). Four MD simulations of 16 ns were performed at physiological and hyperpolarized membrane potentials, with and without restraints of the structure, but all without bound agonist. With the structure unrestrained and a potential of −100 mV, one cation traversed the channel during a transient period of channel hydration; at −200 mV, the channel was continuously hydrated and two cations traversed the channel. With the structure restrained, however, cations did not traverse the channel at either membrane potential, even though the channel was continuously hydrated. The overall results show that cation selective transport through the nAChR channel is governed by electrostatic interactions to achieve charge selectivity, but ion translocation relies on channel hydration, facilitated by a trans-membrane field, coupled with dynamic fluctuations of the channel structure.     

Dynamic properties of the N-terminal swapped dimer of ribonuclease A

Bovine pancreatic ribonuclease (RNase A) forms two 3-dimensional domain-swapped dimers with different quaternary structures. One dimer is characterized by the swapping of the C-terminal region (C-Dimer) and presents a rather loose structure. The other dimer (N-Dimer) exhibits a very compact structure with exchange of the N-terminal helix. Here we report the results of a molecular dynamics/essential dynamics (MD/ED) study carried out on the N-Dimer. This investigation, which represents the first MD/ED analysis on a three-dimensional domain-swapped enzyme, provides information on the dynamic properties of the active site residues as well as on the global motions of the dimer subunits. In particular, the analysis of the flexibility of the active site residues agrees well with recent crystallographic and site-directed mutagenesis studies on monomeric RNase A, thus indicating that domain swapping does not affect the dynamics of the active sites. A slight but significant rearrangement of N-Dimer quaternary structure, favored by the formation of additional hydrogen bonds at subunit interface, has been observed during the MD simulation. The analysis of collective movements reveals that each subunit of the dimer retains the functional breathing motion observed for RNase A. Interestingly, the breathing motion of the two subunits is dynamically coupled, as they open and close in phase. These correlated motions indicate the presence of active site intercommunications in this dimer. On these bases, we propose a speculative mechanism that may explain negative cooperativity in systems preserving structural symmetry during the allosteric transitions.     

Can molecular dynamics simulations improve the structural accuracy and virtual screening performance of GPCR models?

The determination of G protein-coupled receptor (GPCR) structures at atomic resolution has improved understanding of cellular signaling and will accelerate the development of new drug candidates. However, experimental structures still remain unavailable for a majority of the GPCR family. GPCR structures and their interactions with ligands can also be modelled computationally, but such predictions have limited accuracy. In this work, we explored if molecular dynamics (MD) simulations could be used to refine the accuracy of in silico models of receptor-ligand complexes that were submitted to a community-wide assessment of GPCR structure prediction (GPCR Dock). Two simulation protocols were used to refine 30 models of the D₃ dopamine receptor (D₃R) in complex with an antagonist. Close to 60 μs of simulation time was generated and the resulting MD refined models were compared to a D₃R crystal structure. In the MD simulations, the receptor models generally drifted further away from the crystal structure conformation. However, MD refinement was able to improve the accuracy of the ligand binding mode. The best refinement protocol improved agreement with the experimentally observed ligand binding mode for a majority of the models. Receptor structures with improved virtual screening performance, which was assessed by molecular docking of ligands and decoys, could also be identified among the MD refined models. Application of weak restraints to the transmembrane helixes in the MD simulations further improved predictions of the ligand binding mode and second extracellular loop. These results provide guidelines for application of MD refinement in prediction of GPCR-ligand complexes and directions for further method development.