Preparation of the \(\mu\) opioid receptor with ligand#
This is a complex build system as it has several components, the protein, a sodium ion, the ligand and of course the membrane.
from htmd.ui import *
from htmd.home import home
#get the files
shutil.copytree(home(dataDir="mor"),'/tmp/testmor/pdb')
os.chdir('/tmp/testmor')
path='./01_prepare/'
2024-06-12 13:36:44,454 - numexpr.utils - INFO - Note: NumExpr detected 20 cores but "NUMEXPR_MAX_THREADS" not set, so enforcing safe limit of 8.
2024-06-12 13:36:44,454 - numexpr.utils - INFO - NumExpr defaulting to 8 threads.
2024-06-12 13:36:44,569 - 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+20.g8b8902599).
%ls /tmp/testmor/pdb
4dkl.pdb ff.rtf MOL.cif QM-min.pdb
ff.prm membrane80by80C36.pdb MOL.frcmod sod.pdb
Build#
#Protein 4dkl is taken from opm
prot = Molecule('pdb/4dkl.pdb')
prot.filter('(protein and noh and chain B) or (water and within 5 of (chain B and protein))')
pcenter = np.mean(prot.get('coords','protein'), axis=0)
prot = autoSegment(prot, sel='protein')
prot = systemPrepare(prot)
2024-06-12 13:38:49,041 - moleculekit.molecule - INFO - Removed 5120 atoms. 2262 atoms remaining in the molecule.
2024-06-12 13:38:49,076 - moleculekit.tools.autosegment - INFO - Created segment P0 between resid 65 and 263.
2024-06-12 13:38:49,077 - moleculekit.tools.autosegment - INFO - Created segment P1 between resid 270 and 352.
2024-06-12 13:38:49,184 - 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-12 13:38:51,090 - moleculekit.tools.preparation - INFO - Modified residue ASP 114 B to ASH
2024-06-12 13:38:51,091 - moleculekit.tools.preparation - INFO - Modified residue CYS 140 B to CYX
2024-06-12 13:38:51,091 - moleculekit.tools.preparation - INFO - Modified residue HIS 171 B to HID
2024-06-12 13:38:51,092 - moleculekit.tools.preparation - INFO - Modified residue CYS 217 B to CYX
2024-06-12 13:38:51,092 - moleculekit.tools.preparation - INFO - Modified residue HIS 223 B to HID
2024-06-12 13:38:51,092 - moleculekit.tools.preparation - INFO - Modified residue HIS 297 B to HID
2024-06-12 13:38:51,093 - moleculekit.tools.preparation - INFO - Modified residue HIS 319 B to HIE
2024-06-12 13:38:51,094 - moleculekit.tools.preparation - WARNING - Dubious protonation state: the pKa of 1 residues is within 1.0 units of pH 7.4.
2024-06-12 13:38:51,094 - moleculekit.tools.preparation - WARNING - Dubious protonation state: ASP 114 B (pKa= 7.85)
prot.view()
#Add sodium in the receptor
sod = Molecule('pdb/sod.pdb')
sod.set('segid','S1')
prot.append(sod)
#Use a POPC membrane created with vmd and C36
memb = Molecule('pdb/membrane80by80C36.pdb')
mcenter = np.mean(memb.get('coords'),axis=0)
memb.moveBy(pcenter-mcenter)
mol = prot.copy()
mol.append(memb, collisions=True) # Append membrane and remove colliding atoms
#Add ligand, previously parametrized using gaamp
lig = Molecule('pdb/MOL.cif')
lig.set('segid','L')
lcenter = np.mean(lig.get('coords'),axis=0)
newlcenter = [np.random.uniform(-10, 10), np.random.uniform(-10, 10), 43]
lig.rotateBy(uniformRandomRotation(), lcenter)
lig.moveBy(newlcenter - lcenter)
mol.append(lig)
#Add water
coo = mol.get('coords','lipids or protein')
m = np.min(coo,axis=0) + [0,0,-5]
M = np.max(coo,axis=0) + [0,0,20]
mol = solvate(mol, minmax=np.vstack((m,M)))
#Build
topos = amber.defaultTopo() + ['pdb/MOL.cif']
params = amber.defaultParam() + ['pdb/MOL.frcmod']
mol = amber.build(mol, topo=topos, param=params, outdir=os.path.join(path,'build'), saltconc=0.15)
2024-06-12 13:39:40,456 - moleculekit.molecule - INFO - Removed 327 residues from appended Molecule due to collisions.
2024-06-12 13:39:40,768 - htmd.builder.solvate - INFO - Using water pdb file at: /home/sdoerr/Work/htmd/htmd/share/solvate/wat.pdb
2024-06-12 13:39:41,223 - htmd.builder.solvate - INFO - Replicating 8 water segments, 2 by 2 by 2
Solvating: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 8/8 [00:10<00:00, 1.26s/it]
2024-06-12 13:39:51,927 - htmd.builder.solvate - INFO - 9644 water molecules were added to the system.
2024-06-12 13:39:54,236 - htmd.builder.amber - INFO - Detecting disulfide bonds.
2024-06-12 13:39:54,245 - 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-12 13:39:55,271 - htmd.builder.amber - INFO - Starting the build.
2024-06-12 13:39:57,377 - htmd.builder.amber - INFO - Finished building.
2024-06-12 13:39:58,428 - moleculekit.writers - WARNING - Field "resid" of PDB overflows. Your data will be truncated to 4 characters.
2024-06-12 13:40:03,262 - 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-12 13:40:03,614 - htmd.builder.ionize - INFO - Adding 15 anions + 0 cations for neutralizing and 68 ions for the given salt concentration 0.15 M.
2024-06-12 13:40:13,210 - htmd.builder.amber - INFO - Starting the build.
2024-06-12 13:40:15,339 - htmd.builder.amber - INFO - Finished building.
2024-06-12 13:40:16,352 - moleculekit.writers - WARNING - Field "resid" of PDB overflows. Your data will be truncated to 4 characters.
2024-06-12 13:40:20,922 - 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-12 13:40:23,463 - moleculekit.tools.sequencestructuralalignment - INFO - Alignment #0 was done on 282 residues: 2-285
Equilibrate#
from htmd.protocols.equilibration_v3 import Equilibration
from htmd.mdengine.acemd.acemd import GroupRestraint
# Use a 10A flat bottom potential to prevent the ligand from diffusing from original position during equilibration
width = np.array([10, 10, 10])
flatbot = GroupRestraint('segname L and noh', width, [(5, '0ns')])
md = Equilibration()
md.runtime = 40
md.timeunits = 'ns'
md.temperature = 300
md.restraints = [flatbot] + md.defaultEquilRestraints('20ns')
md.useconstantratio = True
md.write(os.path.join(path,'build'), os.path.join(path,'equil'))
2024-06-12 13:40:35,913 - 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]
# Visualize the flat bottom potential box
mol.view('not water')
fbcentre = mol.get('coords', sel='segid L').mean(axis=0).squeeze()
b = VMDBox(np.vstack((fbcentre - width/2, fbcentre + width/2)).T.flatten())
2024-06-12 13:40:37,971 - py.warnings - WARNING - /tmp/ipykernel_388679/442516882.py:3: RuntimeWarning: Mean of empty slice.
fbcentre = mol.get('coords', sel='segid L').mean(axis=0).squeeze()
2024-06-12 13:40:37,990 - py.warnings - WARNING - /home/sdoerr/miniforge3/envs/htmd/lib/python3.10/site-packages/numpy/core/_methods.py:184: RuntimeWarning: invalid value encountered in divide
ret = um.true_divide(
/usr/local/lib/vmd/vmd_LINUXAMD64: /lib/x86_64-linux-gnu/libGL.so.1: no version information available (required by /usr/local/lib/vmd/vmd_LINUXAMD64)
mdx = LocalGPUQueue()
mdx.submit(os.path.join(path, 'equil'))
mdx.wait()
Production#
# Read in the last frame of the equilibration
mol = Molecule(os.path.join(path,'equil','structure.psf'))
mol.read(os.path.join(path,'equil','output.xtc'))
mol.dropFrames(keep=mol.numFrames-1)
from htmd.protocols.production_v6 import Production
from htmd.mdengine.acemd.acemd import GroupRestraint
# Apply a flat bottom potential to prevent the ligand from entering from periodic image of the protein
width = np.array([70, 70, 60])
# Center the box at residue 218 which is on the upper side of the protein
fbcentre = mol.get('coords', sel='protein and resid 218').mean(axis=0).squeeze()
flatbot = GroupRestraint('segname L and noh', width, [(5, '0ns')], fbcentre=fbcentre)
md = Production()
md.runtime = 50
md.timeunits = 'ns'
md.temperature = 300
md.restraints = flatbot
md.write(os.path.join(path,'equil'), os.path.join(path,'prod'))
mol.view('not water')
b = VMDBox(np.vstack((fbcentre - width/2, fbcentre + width/2)).T.flatten())