Computational Methods

Summary. 
The long-term goal of the Sakmar Laboratory is to provide a comprehensive understanding of ligand recognition in GPCRs.  Since the determination of the first high-resolution structure of a GPCR, the visual photoreceptor rhodopsin, the Sakmar Laboratory has used computer homology modeling and molecular dynamics (MD) simulations to build receptor models in native-like phospholipid bilayers.  These all-atom models typically contain more than 50,000 atoms and are simulated on a sub-microsecond timescale using massive parallel processing (MPP) architectures on the TeraGrid, such as the PSC BigBen Cray XT3.  When combined with biochemical and biophysical experiments to map the thermodynamics of ligand binding in rhodopsin, this computational approach provides significant insight toward understanding ligand binding and recognition in GPCRs.  Recently, the Sakmar Laboratory collaborated with Dr. Xavier Periole and Prof. Siewert-Jan Marrink, University of Groningen, The Netherlands, [http://md.chem.rug.nl/~marrink/science.html] to develop coarse-grain molecular dynamics methods to model the behavior of multiple GPCRs in bilayers, including the spontaneous self-assembly of receptor dimers and oligomers. 


Homology Modeling and All Atom Molecular Dynamics Simulations. 
GPCRs constitute a superfamily of heptahelical receptors that respond to a variety of external stimuli such as amines, peptides, enzymes, hormones, odors, and light.  Rhodopsin is the defining member of the largest family group, the rhodopsin-like class A GPCRs, with 569 possible candidates identified in the human genome.  Class A GPCRs comprise the opsins and protein, peptide, and cationic amine receptor subfamilies, constituting ~90% of all known GPCRs.  The sequences of related GPCRs can be used to build homology models based on rhodopsin and other GPCRs for which high-resolution structural data are available, rendering in silico studies of class A GPCRs feasible.
The Sakmar Laboratory recently performed MD simulations of the β2-adrenergic receptor and compared the results to what had been done earlier with rhodopsin.  In order to investigate the binding mode of catecholamines in the context of the new structures, an MD simulation was devised to place the receptor bound to carazolol, an inverse agonist ligand, in a typical bilayer membrane environment.  Next, to address the question of whether or not adrenaline, the receptor agonist, would be stable in fundamentally the same receptor conformations used to bind carazolol additional simulations were performed.  The nanosecond timescale of these MD simulations is certainly not sufficient to capture the slow transition to the active receptor conformation.  However, a principal component analysis of the movements of the backbone Cα atoms demonstrated a global change of the receptor structure in its transition from the carazolol to the adrenaline bound form.  The change is consistent with the elongated structure of the antagonist pushing the extracellular ends of transmembrane helices 2 and 6 outward into the bilayer.


Molecular Dynamics Simulations of Ligand Binding Pathways. 
Many GPCRs have hydrophobic ligands that are likely to partition into the membrane bilayer.  In extended MD simulations of rhodopsin in the membrane, spontaneous structural fluctuations were observed leading to a temporary opening of the protein molecular surface to provide access to the ionone ring moiety of the 11-cis-retinal ligand.  In contrast, the homologous region of a homology model of the chicken red cone pigment, iodopsin, which is known for its very fast retinal binding kinetics, revealed a continuous opening with direct view on the ionone ring of the ligand.
Careful analysis of MD models of rhodopsin in the membrane suggested that this transient pore in the membranous region of the receptor between transmembrane helices 5 and 6 might represent the primary ligand entry site.  To test this hypothesis, Dr. Huber performed preliminary multi-nanosecond free energy calculations of ligand binding in a membrane model of rhodopsin based on reversible MD simulations using umbrella potentials.  The weighted histogram analysis method (WHAM) was employed to eliminate the effects of the biasing umbrella potentials and to recover potential of mean force (PMF) free energy profiles.  The results were compared with irreversible steering MD simulations applying the Jarzynski nonequilibrium equality, and with (multi-configurational) thermodynamic integration.
The ligand binding and un-binding simulations supported the notion that the transient pore structure might act as a channel for the ligand to enter the receptor active site.  Site-directed mutagenesis experiments confirmed that substitutions of amino acid side-chains in the putative binding pocket affected the kinetics of ligand binding and release. 


Coarse-Grain Molecular Dynamics (CGMD) Simulations of Complex Systems. 
Coarse-grain molecular dynamics (CGMD) simulations hold considerable promise to elucidate the behavior of membrane receptors in native-like membrane environments.  By reducing the degree of complexity of the system from atoms to chemical groups formed by 3–6 non-hydrogen atoms, this technique allows the required increase in system size and time scale as compared with a full atomistic approach.  CGMD has been validated to assure that the physicochemical properties of the model system are conserved.  The Sakmar Laboratory recently extended the CGMD approach to address the question of how the physicochemical properties of the membrane lipids might affect self-assembly of GPCRs into higher-order structures.  First, simulations of rhodopsin as a monomer were used to characterize the response of both the protein and the bilayer to the presence of hydrophobic mismatch.  Next, multi-copy rhodopsin simulations (16 proteins per unit cell) at a protein-to-lipid ratio 1:100 in the same bilayer environments were carried out.  Rhodopsin progressively self-assembled into aggregates and ordered linear arrays.  Localized adaptation of the membrane bilayer to the protein and surface complementarity were key factors that defined the rate, extent, and orientational preference of protein–protein association.


 

figure 1
figure 2
figure 3
figure 4
figure 5
figure 6
figure 7