Open-Boundary Molecular Mechanics/Coarse-Grained Framework for Simulations of Low-Resolution G-Protein-Coupled Receptor-Ligand Complexes

Abstract

G-protein-coupled receptors (GPCRs) constitute as much as 30% of the overall proteins targeted by FDA-approved drugs. However, paucity of structural experimental information and low sequence identity between members of the family impair the reliability of traditional docking approaches and atomistic molecular dynamics simulations for in silico pharmacological applications. We present here a dual-resolution approach tailored for such low-resolution models. It couples a hybrid molecular mechanics/coarse-grained (MM/CG) scheme, previously developed by us for GPCR-ligand complexes, with a Hamiltonian-based adaptive resolution scheme (H-AdResS) for the solvent. This dual-resolution approach removes potentially inaccurate atomistic details from the model while building a rigorous statistical ensemble-the grand canonical one-in the high-resolution region. We validate the method on a well-studied GPCR-ligand complex, for which the 3D structure is known, against atomistic simulations. This implementation paves the way for future accurate in silico studies of low-resolution ligand/GPCRs models.

Figure 1

Scheme of the OB-MM/CG setup for hGPCR–ligand complexes. Protein residues belonging to the MM, I, and CG regions are represented in blue, orange, and black, respectively. The ligand is represented in red in the binding site and is modeled at the MM level. The two atomistic MM_w hemispheres are coupled to the coarse-grained reservoir CG_w through the hybrid region HY.

Figure 2

Snapshots with the uncorrected (a) and corrected (b) membrane potentials in OB-MM/CG simulations. The regions parallel to the membrane considered for the calculation of the planar radial density are highlighted with a dashed line. (c) Planar radial density computed from the center of the upper MM_w region in a disk of height 2.0 nm just above the membrane planes, before and after the correction to the membrane potential. The yellow-shaded area corresponds to the presence of the protein. The density value lower than 1 g/cm³ computed in the corrected simulation is attributable to the natural depletion layer in proximity of the membrane plane.

Chart 1. Chemical Structure of the Inverse Agonist S-Carazolola

Figure 3

Histograms of the tetrahedral order parameter q_tet for the water molecules above the binding site, calculated using all-atoms, MM/CG, and OB-MM/CG approaches.

Figure 4

Histograms of the distances of hydrogen bonds between the receptor and the ligand. (a) Distance Asp113(3.32) C_COO–S-car H_OH. (b) Distance Ser203(5.42) O_CO–S-car H_NH. (c) Distance Asn312(7.39) O_CO–S-car N_NH2⁺.

Figure 5

Reorientational tcfs of the ligand in the binding pocket. When the tcfs are fitted with a “model-free” function of the form C(t) = S² + (1 – S²)e^–t/τ, the generalized order parameter S² takes the values 0.963, 0.972, and 0.963 for the all-atom MD, MM/CG, and OB-MM/CG case, respectively.

Figure 6

RMSFs of the C^α atoms in the MM region. The overall protein backbone flexibility is in line with that of previously published atomistic and MM/CG models of β2-AR, with differences typically smaller than 0.3 Å. The shaded areas correspond to the residues stabilizing the ligand through hydrophobic interactions or hydrogen bonds.

Figure 7

Histograms of the number of water molecules at a maximum distance of 5 Å from the ligand in the binding cavity.