Membrane Modeling

Written by Stefan Gahbauer, 2019/11/03

In order to account for ligand desolvation and electrostatic interactions in the low-dielectric environment of the hydrophobic membrane core, a lipid-bilayer is generated around the target receptor and included in the docking score grid generation. Aiming at a fast, robust and computationally effective equilibration of the lipid bilayer around the embedded transmembrane receptor, coarse-grained (CG) molecular dynamics (MD) simulations and (if needed) subsequent atomistic simulations are employed.

Required software and datasets

Gromacs (v5 or newer) - Molecular Dynamics software package (http://manual.gromacs.org/)

CHARMM36m force field (http://mackerell.umaryland.edu/charmm_ff.shtml)

MARTINI Coarse-grained force field parameters(http://cgmartini.nl/)

DSSP - Secondary Structure assignment (https://swift.cmbi.umcn.nl/gv/dssp/ , https://anaconda.org/salilab/dssp)

martinize.py - Coarse-graining atomistic protein structures (http://cgmartini.nl/index.php/tools2/proteins-and-bilayers)

insane.py - INSerting proteins in coarse-grained MembrANE (http://cgmartini.nl/index.php/tools2/proteins-and-bilayers)

initram.sh and backward.py - Conversion of coarse-grained system to atomistic resolution (http://cgmartini.nl/index.php/tools2/resolution-transformation)

1) Setting up the coarse-grained system

This is automated in the script:

/mnt/nfs/home/stefan/zzz.scripts/INSERT-MEMBRANE/FILES/0001-prepare-protein-CG-membrane.sh

1.1) Prepare your files

Copy the script above to your working directory.

Copy your rec.pdb to your working directory.

If your rec.pdb has gaps, e.g. unresolved loops between transmembrane helices in case of GPCRs, try to model missing residues.

One way is to use MODELLER following https://salilab.org/modeller/wiki/Missing%20residues.

Corresponding input scripts for modeller can be found in:

/mnt/nfs/home/stefan/zzz.scripts/INSERT-MEMBRANE/FILES/modeller

1.2) Run the script

Login to gimel2.

./0001-prepare-protein-CG-membrane.sh

The script reads rec.pdb and copies all other required files from

/mnt/nfs/home/stefan/zzz.scripts/INSERT-MEMBRANE/FILES/gromacs

Workflow

a. Generate CHARMM36m force field parameters of your protein in a Gromacs-readable format.

Used tool: gmx pdb2gmx

Output files are stored in the generated pdb2gmx directory

-conf.gro / conf.pdb - Gromacs coordinate file

-topol.top / Protein-atomistic.itp - Gromacs topology file, i.e. force field description of your input structure

-posre.itp - Position restraints for heavy atoms of atomistic protein strucutre.

b. Build coarse-grained structure

Used tool: martinize.py

Output files are stored in the generated martini directory.

-chain_.ssd - Output from the DSSP program that is called by martinize.py

-prot-cg.pdb - Coarse-grained protein structure

-prot-cg.top - Coarse-grained Martini topology of system

-Protein.itp - Coarse-grained Martini description of Protein structure

-prot-rot.pdb - Coarse-grained protein structure aligned along z-axis of the simulation box according to the proteins first principal component axis. This ensures the correct placement of the protein during membrane preparation. You may have to adjust the orientation of your input structure prior to membrane modeling.

c. Build coarse-grained membrane

Used tool: insane.py

Here, a lipid bilayer will be created around the protein structure (in the x/y-plane) and water will be added to the system. The default box shape is rectangular and the size is set to x,y=10nm, z=11nm. This can be changed in the ./insane.py command line. The default lipid type is POPC, you can change that to arbitrary lipid compositions using the -l and -u options of of insane.py.

Output files are stored in the martini directory.

-out.top / topol-cg.top - Topology of coarse-grained system

-cg-membrane.gro/.pdb - Coarse-grained system coordinates. Carefully inspect and visualize the cg-membrane.pdb.

Use PyMOL to check if you're protein is embedded correctly in the lipid bilayer.

2) Simulating the coarse-grained system

This is automated in the script:

/mnt/nfs/home/stefan/zzz.scripts/INSERT-MEMBRANE/FILES/0002-run-CG-Minimization-and-MD.sh

2.1) Run the script

Copy the script above to your working directory.

./0002-run-CG-Minimization-and-MD.sh

Workflow

a. Minimize coarse-grained system

Used tools: gmx grompp , gmx mdrun

gmx grompp generates a single .tpr file that contains all information necessary for running a MD simulation or minimization using gmx mdrun.

Minimization parameters are provided in martini_new-rf_min.mdp. The system will be minimized in 500 steps using steepest descent. The protein structure will be frozen during minimization.

Output files are stored in the martini directory.

-min.tpr - MD run input file

-min.log - Output log file from minimization

-min.trr - Minimization trajectory

-min.gro - Minimized system coordinates

b. Simulate coarse-grained system

MD simulation parameters are provided in martini_v2.x_new-rf.mdp. Strong position restraints are applied on the protein structure during the simulation. The system will be simulated for 50ns.

Output files are stored in the martini directory.

-md.tpr - MD run input file

-md.log - Output log file from simulation

-md.trr - lossless trajectory of simulation

-md.xtc - coordinates of simulation trajectory

-md.gro - coordinates of final simulation snapshot

The simulation will run for roughly 3 hours.

3) Converting coarse-grained system to atomistic resolution and select lipid atoms for grid generation

This is automated in the script:

/mnt/nfs/home/stefan/zzz.scripts/INSERT-MEMBRANE/FILES/0003-backmap-and-lpd-selection.sh

3.1) Run the script

Copy the script above to your working directory.

./0003-backmap-and-lpd-selection.sh

Workflow

a. Backmapping from coarse-grained to atomistic

Used tool: gmx trjconv, initram.sh

gmx trjconv can perform a variety of conversions of MD trajectory, e.g. making molecules broken over the periodic boundary conditions whole again.

initram.sh calls the backward.py program which performs the backmapping of input coarse-grained to atomistic systems, and performs a small series of short minimizations and simulations to relax the backmapped system.

Output files are stored in the generated backmap directory.

-0-backmapped.gro / projected.gro - initial backmapped coordinates

-backmapped.top - Topology of atomistic system

-1-EM*/2-EM* - Output from minimizations

-3-mdpr*/4-mdpr*/5-mdpr*/6-mdpr* - Output from simulations

-backmapped.gro - Coordinates of final backmapped and relaxed system

This may need a few attempts to work all the way through. There is a while-loop that only stops until all relaxation steps have finished.

b. Replacing backmapped protein with initial atomistic protein structure

Used tool: PyMOL script align.pml

The pymol script will align the initial "Gromacs"-protein structure (conf.pdb) onto the backmapped structure and combine the fitted protein coordinates with the coordinates of the lipid and solvent environment.

Output files are stored in the generated backmap/prepare_AA_system directory.

-conf-fitted.pdb - Fitted initial protein structure

-backmapped-environment.pdb - All membrane and water coordinates

-fitted_system.pdb - Complete system containing fitted protein and environment coordinates

Be sure that fitted_system.pdb has the same number of coordinates as backmapped-mol.pdb. If there is a discrepancy there might be an issue with the PyMOL version you're using to run align.pml. Using PyMOL v2 or newer seems to avoid any issues. You can also generate backmapped-environment.pdb manually by taking all POPC and Water cooridnates from backmapped-mol.pdb.

c. Run minimizations of atomistic system

Used tools: gmx grompp, gmx mdrun

Output files are stored in the backmap/prepare_AA_system directory.

Tow minimizations will be calculated.

1) Minimization with frozen protein coordinates: 1,500 steps steepest descent (min_freeze.mdp).

-min_freeze* - Output files of first minimization

2) Minimization of full system: 500 steps (min.mdp).

-min* - Output files of second minimization

d. Select lipid atoms for DOCK grid generation

Used tools: PyMOL script prepare.pml

This will select carbon and hydrogen atoms of the hydrophobic lipid tail segments in a radius of 1.7 nm around the protein and assign them to the atom type "LPD"

You need to provide the rec.pdb to want to use for docking (potentially with missing loops) at that step as xtal-prot.pdb

Output files are stored in the generated backmap/prepare_AA_system/prepare_min directory.

-shell-LPD.pdb - all LPD atoms selected for grid generation

Add these coordinates to your docking protein structure and provide a amb.crg.oxt file adding

C lpd 0.000 LIPID SPHERE

Now you can run blastermaster.