System Building: μ-opioid Receptor in Membrane with morphinan antagonist

In this tutorial, we will showcase how to build a protein system embeded in a membrane with a ligand in one of the compartments for simulating binding. The sample system is a mu-opioid receptor (the protein) in a membrane and a morphinan antagonist (the ligand).

Let’s start by doing some imports and definitions:

from htmd.ui import *
from htmd.home import home
from os.path import join
config(viewer='webgl')
datadir = home(dataDir='mor')
Please cite HTMD: Doerr et al.(2016)JCTC,12,1845.
https://dx.doi.org/10.1021/acs.jctc.6b00049
Documentation: http://software.acellera.com/
To update: conda update htmd -c acellera -c psi4

You are on the latest HTMD version (unpackaged : /home/joao/maindisk/software/repos/Acellera/htmd/htmd).

Prepare the protein

View the file as it comes from the OPM database:

Molecule(join(datadir, '4dkl.pdb')).view()
A Jupyter Widget

Retrieve the structure from OPM, do not keep the DUM atoms, and print the thickness as calculated by OPM:

from htmd.util import opm
prot, thickness = opm('4dkl', keep=False)
thickness
2018-03-17 01:54:49,428 - htmd.molecule.molecule - INFO - Removed 2546 atoms. 4836 atoms remaining in the molecule.
32.0

Remove non-protein atoms, and keep only a monomer.

prot.filter('protein and noh and chain B or water within 5 of (chain B and protein)')
2018-03-17 01:54:51,462 - htmd.molecule.molecule - INFO - Removed 2574 atoms. 2262 atoms remaining in the molecule.
array([   0,    1,    2, ..., 4808, 4809, 4810], dtype=int32)

Automatically detecting segments and assigning names to them.

prot = autoSegment(prot,'protein')
2018-03-17 01:54:53,807 - htmd.builder.builder - INFO - Created segment P0 between resid 65 and 263.
2018-03-17 01:54:53,808 - htmd.builder.builder - INFO - Created segment P1 between resid 270 and 352.

Build the protein

topos  = charmm.defaultTopo() + [join(datadir, 'ff.rtf')]
params = charmm.defaultParam() + [join(datadir, 'ff.prm')]

prot = charmm.build(prot, topo=topos, param=params, outdir='./tmp-build', ionize=False)
2018-03-17 01:54:57,172 - htmd.builder.charmm - INFO - Writing out segments.
2018-03-17 01:54:57,371 - htmd.builder.builder - INFO - One disulfide bond was added
Bond between A: [serial 3005 resid 140 resname CYS chain B segid P0]
             B: [serial 3615 resid 217 resname CYS chain B segid P0]
2018-03-17 01:54:57,912 - htmd.builder.charmm - INFO - Starting the build.
2018-03-17 01:54:58,939 - htmd.builder.charmm - WARNING - Failed to guess coordinates for 1 atoms due to bad angles.
2018-03-17 01:54:58,941 - htmd.builder.charmm - WARNING - Poorly guessed coordinates for 78 atoms.
2018-03-17 01:54:58,941 - htmd.builder.charmm - WARNING - Please check /data/joao/maindisk/software/repos/Acellera/htmd/tutorials/tmp-build/log.txt for further information.
2018-03-17 01:54:58,943 - htmd.builder.charmm - INFO - Finished building.
prot.reps.add(sel='segid P0', style='NewCartoon', color=1)
prot.reps.add(sel='segid P1', style='NewCartoon', color=2)
prot.view()
A Jupyter Widget

Add a sodium atom in the receptor

sod = Molecule(join(datadir, 'sod.pdb'))
sod.set('segid','S1')
prot.append(sod)
prot.reps.add(sel='ions', style='VDW', color='green')
prot.view()
A Jupyter Widget

Embed the protein into a membrane

memb = Molecule(join(datadir, 'membrane80by80C36.pdb'))

Center the membrane onto the protein center

pcenter = np.mean(prot.get('coords','protein'),axis=0)
mcenter = np.mean(memb.get('coords'),axis=0)
memb.moveBy(pcenter-mcenter)

And now embed.

The two are equivalent - append with collisions=True only adds atoms if they do not clash sterically.

mol = prot.copy()
mol.append(memb, collisions=True)
# mol = embed(prot,memb)
2018-03-17 01:55:22,309 - htmd.molecule.molecule - INFO - Removed 339 residues from appended Molecule due to collisions.

Visualize the embedded system

mol.reps.add(sel='protein', style='NewCartoon', color='Secondary Structure')
mol.reps.add(sel='ions', style='VDW', color='green')
mol.reps.add(sel='lipids', style='Lines')
mol.view()
A Jupyter Widget

Add a ligand

import random
lig = Molecule(join(datadir, 'QM-min.pdb'))
lig.set('segid','L');
lcenter = np.mean(lig.get('coords'),axis=0)
newlcenter=[random.uniform(-10, 10), random.uniform(-10, 10),  43 ]
lig.rotateBy(uniformRandomRotation(), lcenter)
lig.moveBy(newlcenter-lcenter)
mol.append(lig)

Put it in a water box

coo = mol.get('coords','noh and (lipids or protein)')
m = np.min(coo, axis=0) + [0, 0, -5]
M = np.max(coo, axis=0) + [0, 0, 20]
smol = solvate(mol, minmax=np.vstack((m,M)))
2018-03-17 01:55:41,667 - htmd.builder.solvate - INFO - Using water pdb file at: /data/joao/maindisk/software/repos/Acellera/htmd/htmd/builder/wat.pdb
2018-03-17 01:55:43,163 - htmd.builder.solvate - INFO - Replicating 8 water segments, 2 by 2 by 2
Solvating: 100%|██████████| 8/8 [00:17<00:00,  2.20s/it]
2018-03-17 01:56:03,643 - htmd.builder.solvate - INFO - 9117 water molecules were added to the system.

Visualize

smol.reps.add(sel='segid L', style='Licorice')
smol.reps.add(sel='water', style='Lines')
smol.view()
A Jupyter Widget

Build with CHARMM

molbuilt = charmm.build(smol, topo=topos, param=params, outdir='./final-build', saltconc=0.15)
2018-03-17 01:56:35,041 - htmd.builder.charmm - INFO - Writing out segments.
2018-03-17 01:56:53,448 - htmd.builder.builder - INFO - One disulfide bond was added
Bond between A: [serial 1212 resid 140 resname CYS chain B segid P0]
             B: [serial 2448 resid 217 resname CYS chain B segid P0]
2018-03-17 01:56:53,971 - htmd.builder.charmm - INFO - Starting the build.
2018-03-17 01:56:54,837 - htmd.builder.charmm - INFO - Finished building.
2018-03-17 01:56:58,567 - htmd.builder.ionize - INFO - Adding 14 anions + 0 cations for neutralizing and 66 ions for the given salt concentration.
2018-03-17 01:56:59,001 - htmd.builder.ionize - INFO - Min distance of ions from molecule: 5A
2018-03-17 01:56:59,002 - htmd.builder.ionize - INFO - Min distance between ions: 5A
2018-03-17 01:56:59,003 - htmd.builder.ionize - INFO - Placing 80 ions.
2018-03-17 01:57:36,506 - htmd.builder.charmm - INFO - Writing out segments.
2018-03-17 01:57:54,859 - htmd.builder.charmm - INFO - Starting the build.
2018-03-17 01:57:55,740 - htmd.builder.charmm - INFO - Finished building.

Visualize built system

molbuilt.view()
A Jupyter Widget

The molbuilt is a Molecule object that contains the built system, but the full contents to run a simulation are located in the outdir (./final-build in this case).