Molecular Dynamics on GPCRs

Surveying the field of molecular dynamics (MD) in biomolecular systems at any point in time, the strongest conclusion you can draw is that the last generation of simulations underestimated the amount of sampling needed to get accurate results. Fortunately, the march of Moore’s Law eventually solves that problem. Unfortunately, in the mean time, the sophistication and size of the systems people want to simulate keeps getting bigger. But if you want to get a feel for what the state of the art is in terms of how big a system can be simulated for how long, you could do a lot worse than check in on the work being done by David E. Shaw Research.

In a recent paper in PNAS, David E. Shaw and the Dror group at Columbia looked at the binding of antagonists of β2-adrenergic receptor: (S)-propranolol, (S)-alprenol and (S)-dihydroalprenolol. Even more interesting than the idea  of doing MD on a GPCR (which has been difficult until recently due to a lack of X-ray structures), is that they effectively carry out ‘blind’ docking (they reported a similar study with Src kinase a few months ago): they drop 10 ligands at least 30 Å from the binding pocket and 12 Å from the receptor surface, start a few simulations lasting between 1-19 µs (most were 1 µs), and look to see if any of the ligands enter the binding site during that time period. This approach certainly has the advantage of conceptual simplicity, but it’s exactly how you don’t carry out docking, as you’d think it would take far too long for anything to happen.

Remarkably, though, out of 71 simulations, 15 binding events were observed for the antagonists, and of those, seven of the final binding poses showed an RMSD < 1 Å to the corresponding heavy ligand atoms in the alprenolol-β2AR crystal structure (PDB code 3NYA). Admittedly, that leaves eight other results where the ligand ended up in one of two wrong poses, which the authors suggest are metastable poses which convert to the correct pose too slowly for the timescale of the simulation. But that this worked at all is pretty impressive.

The authors were also able to extract a putative binding pathway from the simulation, which had a very interesting feature in that before the antagonists enter the binding pocket, the antagonists spend time in what the authors call the “extracellular vestibule”, defined by a region near the Extracellular Loops 2 and 3, and helices 5-7.  Moreover, entry into the vestibule appears to be associated with the highest energetic barrier in the binding event, rather than entry into the binding pocket. The authors suggest the source of the energy barrier is due to substantial dehydration that has to occur at both the ligand and the vestibule for the ligand to enter the vestibule.

Additionally, the authors calculated binding kinetics and energies for alprenolol that were close to the experimental values (for the binding energy, -13.4 kcal/mol versus the experimental value of -12.2 kcal/mol).

What did it take to achieve all this? The MD simulations took up around 100 µs, and the kinetics calculations required around about the same amount of extra simulation time. The binding energy calculations used double-annihilation free energy perturbation (FEP) with replica exchange, using 50 ns or 10 ns per replica, depending on whether the ligand was simulated in the receptor or solvent respectively. Unfortunately, we’re not told how many of the intermediate stages (where the non-bonded parameters are scaled down) were used in the FEP method.

When I hear about MD simulations, I always want to know how big the system was, what force field was used, and how cut-offs were dealt with. So here goes (warning: if you’re not like me, you may find this a very dry paragraph). All of this was done using a version of the CHARMM 27 force field parameters with modified electrostatics, with explicit solvent and lipid molecules present, in all around 50-60,000 atoms. Equilibration took place for what seems by comparison to the production to be a positively miserly 10 ns on a commodity cluster (again, unfortunately we don’t know how big the cluster was) using the Desmond software, with van der Waals and short range cutoff at 9 Å and using the Particle Mesh Ewald (PME) method for long-range electrostatics. For the actual production step, DES Research’s fearsome Anton MD-specialised supercomputer was used, and the cutoff was increased to 13.5 Å. Instead of PME, the long-range electrostatics were handled using the k-space Gaussian Split Ewald method.

It’s interesting to note that a single simulation was also carried using the agonist (R)-isoproterenol, which entered the binding site, but had much more mobility than the antagonists, probably because it takes milliseconds for the receptor to enter the active conformation. But for now it looks like if you have processes you’re interested in taking place over a few hundred of microseconds, it’s becoming feasible to use MD to explore these, even with GPCRs.