Hybrid molecular mechanics/coarse-grained simulation for G-protein coupled receptors

G protein-coupled receptors (GPCRs) are the most frequent protein drug target with about 34% of all FDA approved drugs. Unfortunately, no experimental structures are available for most human GPCRS (about 90%) and the sequence identity shared with the available templates for homology modeling is low, resulting in low-resolution models. Thus, the molecular dockings that can be performed in these models are inaccurate. The hybrid molecular mechanics/coarse-grained (MM/CG) simulation was developed to refine molecular docking poses based on low-quality homology models. This hybrid scheme divides the system in three regions with different resolution (for more inforamtion see: Schneider et al. 2018 (preprint), Leguèbe et al. 2012, Marchiori et al. 2013, Fierro et al. 2019):

• (MM) The region around the ligand binding site and the extracellular part of the receptor, as well as the water droplet, are represented in all-atom resolution. Amber 14SB, GAFF or GAFF2 and TIP3P are the force fields used in the MM region for the receptor, the ligand and water, respectively.
• (CG) Residues far from the ligand binding site towards the inside of the cell are simulated in coarse-grained representation. Each residue is reduced to one bead at the position of the alpha carbon atom. A Gō-like potential is used for the CG region.
• (I) The region between the MM and CG regions is an interface. Residues in this region interact with the MM and CG residues with the respective interactions. The height of the interface layer is chosen in order to prevent direct interactions of MM and CG residues.

The receptor is embedded in five potential walls:

• (1+2) Two planes represent the upper and lower limit of the membrane and are placed at the height of the membrane lipid heads. The planes prevent the water from penetrating the membrane volume. The planes are purely repulsive.
• (3+4) Two hemispheres prevent the water from evaporation (3) and cap the intracellular loops of the protein. The hemispheres are repulsive.
• (5) The fifth wall embraces the transmembrane part of the protein mimicking the interactions with the membrane lipids. Specific residue types have a Lennard-Jones-like interaction with the fifth wall.

What is this server doing?

The set-up of MM/CG simulations involves many consecutive steps using different third party software and many scripts. In order to prevent users from downloading and compiling multiple tools, using different webservers and figuring out different formats of in- and output thereof, this webserver was built. Indeed, starting from a GPCR/ligand complex structure (in PDB format), the server builds the input files for the hybrid MM/CG simulation. The server can work in a step-by-step approach, giving the user the option to change and check important parameters and intermediate results, or in a fully automatic way, using default values. Moreover, at the end of the protocol, the server will run an energy minimization and some equilibration steps.

Example input and output

Example input PDB file: A2A receptor in complex with it‘s antagonist caffeine.
Example output
Example preparation: Click "1-click preparation" or "interactive preparation" for an example preparation and simulation. The whole preparation takes about 5 minutes excluding equilibration and MD. The same result can also be seen in "Example output" above.

Preparation steps

Upload

Choose a name for the simulation and upload a PDB file fulfilling the requirements given below. The preparation can be done in a fully automatic way using default parameters or interactively. The interactive option allows the user to change parameters and analyze the intermediate results.
The modeling and docking of a GPCR/ligand complex can be performed by our GOMoDo webserver. The GOMoDo webserver is directly linked to this server which allows the user to follow the entire procedure of modeling, docking and refinement with a few clicks.

Receptor protonation state

In order to impose protonation states of individual receptor residues that are different from the default states produced by pdb2gmx, the residue names in the input PDB file can be adapted accordingly, e.g. from "HIS" to "HIP" for a doubly protonated histidine.

Example

To run the example system with the A2A receptor in complex with it‘s antagonist caffeine (CFF), do not select a file to upload.

These criteria can help to spot errors in the input PDB file:

• The first residue in the input file should have residue number "1".
• Only the receptor and ONE ligand are present in the file. No water, no lipids, no cryo-agents.
• Residue names and atom names have to match the ones used in the AMBER force field, e.g. "HIS" and not "HSE", "HSP" or "HSD" for histidine.
• Capping group names have to match the ones used in the AMBER force field, e.g. "ACE" or "NME", but not "NAC"
• Capping groups, if present, have to have their own residue number.
• The residue number of the ligand as well as the atom numbers of the ligand should continue the numbers of the receptor, e.g. not ligand residue number "1" after receptor residue number "330" and not ligand atom numbers "1", "2"... after receptor atom number "3500".
• The residue ID of the ligand (e.g. "CFF" for caffeine) cannot be a residue name or capping group name from the AMBER force field.
• Receptor atom lines start with "ATOM" and ligand lines with "ATOM" or "HETATM".

Receptor positioning into the membrane planes

The PPM webserver (positioning of proteins in membranes) is used to position the uploaded protein inside the membrane planes.

The PPM server calculates rotational and translational positions of transmembrane and peripheral proteins in membranes using their 3D structure (PDB coordinate file) as input. It can be applied to newly determined experimental protein structures or theoretical models.

Reference:
Lomize M.A., Pogozheva I,D, Joo H., Mosberg H.I., Lomize A.L. OPM database and PPM web server: resources for positioning of proteins in membranes. Nucleic Acids Res., 2012, 40(Database issue):D370-6

Receptor topology construction

The topology of the receptor (without the ligand) is built using pdb2gmx. Protonation states are assigned by pdb2gmx (default). The Amber 14SB force field is used. The TIP3P force field is defined for the water.

Ligand topology construction

The protonation state at pH=7 is assigned to the ligand by openbabel. The topology of the ligand is created using ACPYPE. The GAFF or GAFF2 force field for organic and pharmaceutical molecules and AM1-BCC charges are used for the ligand parametrization. The ligand can be neutral or charged.

Solvation of the system

The receptor and ligand topologies are combined, the simulation box is created and water is added to the box. The position of the receptor in the box can be corrected by changing the parameter below, use the default if no problems occur.

Shaping of the water hemisphere and partial coarse-graining of the Receptor

Water far away from the binding site will be deleted in order to create a hemispherical water droplet around the extracellular site of the protein. A region of the receptor towards the intracellular site will be coarse-grained based on the chosen cutoff (use the default if unsure).

Minimization

Minimize the system using the steepest descent method. A plot of the potential energy can be seen in the Results section after finishing the set up.

Surface construction

The five surfaces are built around the receptor. Two hemispheres capping the intracellular part of the protein as well as the extracellular part and the water hemisphere. One surface embracing the receptor and mimicking the interactions with membrane lipids and two planes at the height of the lipids polar heads preventing the water from penetrating the membrane volume. Surface parts close to the ligand binding site are deleted.

Short equilibration and production hybrid MM/CG MD simulation

Scripts for the equilibration and MD run are prepared and sent to the queue. As soon as the equilibration starts, estimated completion times will be shown in the Results section. If not logged in, you can retrieve your results by saving the URL of the Results section. If logged in, all your setups will be shown in a comprehensive table in the Results section. Registration does not require any personal information, just a username and a password.

Analysis

The minimization, the short equilibration and the MD can be analyzed in the Results section as soon as the individual stage is completed. 3D visualizations, potential energy and temperature plots are provided. All files to continue and further analyze the simulation can be downloaded.

Run times

The preparation takes about 6 min for the example system while the optional equilibration (2ns) and first production (2ns) take 20 h together. The extensions take 10 h each. To get 10 ns production on the server it takes about 2.5 days runtime for the example system. Downloading the system after the preparation and running the equilibration and production on more processors will speed up the whole simulation by up to one order of magnitude (10x).

Continue simulations

In order to continue the minimization, equilibration or MD simulation, files from the preapration can be downloaded using the respective button in the Results section. A script to download the GROMACS code together with the MM/CG patch can be obtained in the Download Code section or here.