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 moleculekit.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).