Ctrl+K

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
from moleculekit.config import config

config(viewer='webgl')
datadir = home(dataDir='mor')

2024-06-11 15:51:03,138 - numexpr.utils - INFO - Note: NumExpr detected 20 cores but "NUMEXPR_MAX_THREADS" not set, so enforcing safe limit of 8.
2024-06-11 15:51:03,138 - numexpr.utils - INFO - NumExpr defaulting to 8 threads.
2024-06-11 15:51:03,219 - rdkit - INFO - Enabling RDKit 2022.09.1 jupyter extensions

Please cite HTMD: Doerr et al.(2016)JCTC,12,1845. https://dx.doi.org/10.1021/acs.jctc.6b00049
HTMD Documentation at: https://software.acellera.com/htmd/

You are on the latest HTMD version (2.3.28+0.g1f64e666d.dirty).

Prepare the protein#

View the file as it comes from the OPM database:

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

NGLWidget()

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

from moleculekit.opm import get_opm_pdb
prot, thickness = get_opm_pdb('4dkl', keep=False)
thickness

32.0

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

prot.filter('(protein and noh and chain B) or (water and within 5 of (chain B and protein))')

2024-06-11 15:51:05,557 - moleculekit.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, sel="protein")
prot.set("segid", "W", sel="water")

2024-06-11 15:51:05,587 - moleculekit.tools.autosegment - INFO - Created segment P0 between resid 65 and 263.
2024-06-11 15:51:05,588 - moleculekit.tools.autosegment - INFO - Created segment P1 between resid 270 and 352.

prot = systemPrepare(prot)

2024-06-11 15:51:05,684 - moleculekit.tools.preparation - WARNING - Both chains and segments are defined in Molecule.chain / Molecule.segid, however they are inconsistent. Protein preparation will use the chain information.

---- Molecule chain report ----
Chain A:
    First residue: HOH   717
    Final residue: HOH   717
Chain B:
    First residue: MET    65
    Final residue: HOH   735
---- End of chain report ----

2024-06-11 15:51:07,626 - moleculekit.tools.preparation - INFO - Modified residue ASP   114 B to ASH
2024-06-11 15:51:07,627 - moleculekit.tools.preparation - INFO - Modified residue CYS   140 B to CYX
2024-06-11 15:51:07,627 - moleculekit.tools.preparation - INFO - Modified residue HIS   171 B to HID
2024-06-11 15:51:07,627 - moleculekit.tools.preparation - INFO - Modified residue CYS   217 B to CYX
2024-06-11 15:51:07,628 - moleculekit.tools.preparation - INFO - Modified residue HIS   223 B to HID
2024-06-11 15:51:07,628 - moleculekit.tools.preparation - INFO - Modified residue HIS   297 B to HID
2024-06-11 15:51:07,628 - moleculekit.tools.preparation - INFO - Modified residue HIS   319 B to HIE
2024-06-11 15:51:07,629 - moleculekit.tools.preparation - WARNING - Dubious protonation state: the pKa of 1 residues is within 1.0 units of pH 7.4.
2024-06-11 15:51:07,630 - moleculekit.tools.preparation - WARNING - Dubious protonation state:    ASP   114 B (pKa= 7.85)

prot.reps.add(sel='segid P0', style='NewCartoon', color=1)
prot.reps.add(sel='segid P1', style='NewCartoon', color=2)
prot.view()

NGLWidget()

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()

NGLWidget()

Embed the protein into a membrane#

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

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)

2024-06-11 15:51:08,522 - moleculekit.molecule - INFO - Removed 319 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()

NGLWidget()

Add a ligand#

import random
lig = Molecule(join(datadir, 'MOL.cif'))
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)))

2024-06-11 15:51:10,118 - htmd.builder.solvate - INFO - Using water pdb file at: /home/sdoerr/Work/htmd/htmd/share/solvate/wat.pdb
2024-06-11 15:51:10,551 - htmd.builder.solvate - INFO - Replicating 8 water segments, 2 by 2 by 2
Solvating: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 8/8 [00:10<00:00,  1.27s/it]
2024-06-11 15:51:21,290 - htmd.builder.solvate - INFO - 8853 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()

NGLWidget()

Build with AMBER#

molbuilt = amber.build(smol, param=[join(datadir, 'MOL.frcmod')], topo=[join(datadir, 'MOL.cif')], outdir='./final-build', saltconc=0.15)

2024-06-11 15:51:26,607 - htmd.builder.amber - INFO - Detecting disulfide bonds.
2024-06-11 15:51:26,614 - htmd.builder.builder - INFO - One disulfide bond was added

Disulfide Bond between: UniqueResidueID<resname: 'CYX', chain: 'B', resid: 78, insertion: '', segid: 'P0'>
                   and: UniqueResidueID<resname: 'CYX', chain: 'B', resid: 155, insertion: '', segid: 'P0'>

2024-06-11 15:51:27,682 - htmd.builder.amber - INFO - Starting the build.
2024-06-11 15:51:29,609 - htmd.builder.amber - INFO - Finished building.
2024-06-11 15:51:30,594 - moleculekit.writers - WARNING - Field "resid" of PDB overflows. Your data will be truncated to 4 characters.
2024-06-11 15:51:35,176 - htmd.builder.builder - WARNING - Found cis peptide bond in 1 frames: [0] in the omega diheral "Angle of (HID 160 CA  ) (HID 160 C  ) (PRO 161 N  ) (PRO 161 CA  ) " with indexes [2528, 2541, 2543, 2553]
2024-06-11 15:51:35,516 - htmd.builder.ionize - INFO - Adding 15 anions + 0 cations for neutralizing and 64 ions for the given salt concentration 0.15 M.
2024-06-11 15:51:44,191 - htmd.builder.amber - INFO - Starting the build.
2024-06-11 15:51:46,213 - htmd.builder.amber - INFO - Finished building.
2024-06-11 15:51:47,141 - moleculekit.writers - WARNING - Field "resid" of PDB overflows. Your data will be truncated to 4 characters.
2024-06-11 15:51:51,546 - htmd.builder.builder - WARNING - Found cis peptide bond in 1 frames: [0] in the omega diheral "Angle of (HID 160 CA  ) (HID 160 C  ) (PRO 161 N  ) (PRO 161 CA  ) " with indexes [2528, 2541, 2543, 2553]
2024-06-11 15:51:51,569 - py.warnings - WARNING - /home/sdoerr/miniforge3/envs/htmd/lib/python3.10/site-packages/Bio/pairwise2.py:278: BiopythonDeprecationWarning: Bio.pairwise2 has been deprecated, and we intend to remove it in a future release of Biopython. As an alternative, please consider using Bio.Align.PairwiseAligner as a replacement, and contact the Biopython developers if you still need the Bio.pairwise2 module.
  warnings.warn(

2024-06-11 15:51:54,048 - moleculekit.tools.sequencestructuralalignment - INFO - Alignment #0 was done on 282 residues: 2-285

Visualize built system#

molbuilt.view()

2024-06-11 15:51:54,436 - moleculekit.writers - WARNING - Field "resid" of PDB overflows. Your data will be truncated to 4 characters.

NGLWidget()

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

On this page