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-12-05 12:45:59,431 - numexpr.utils - INFO - Note: NumExpr detected 20 cores but "NUMEXPR_MAX_THREADS" not set, so enforcing safe limit of 8.
2024-12-05 12:45:59,431 - numexpr.utils - INFO - NumExpr defaulting to 8 threads.
/home/sdoerr/miniforge3/envs/htmd/lib/python3.10/site-packages/pandas/core/computation/expressions.py:21: UserWarning: Pandas requires version '2.8.4' or newer of 'numexpr' (version '2.7.3' currently installed).
  from pandas.core.computation.check import NUMEXPR_INSTALLED
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.4.2+20.g74f7135a7.dirty).
%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-12-05 12:46:02,233 - moleculekit.molecule - INFO - Removed 5120 atoms. 2262 atoms remaining in the molecule.
2024-12-05 12:46:02,268 - moleculekit.tools.autosegment - INFO - Created segment P0 between resid 65 and 263.
2024-12-05 12:46:02,268 - moleculekit.tools.autosegment - INFO - Created segment P1 between resid 270 and 352.
2024-12-05 12:46:02,308 - 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-12-05 12:46:04,219 - moleculekit.tools.preparation - INFO - Modified residue ASP   114 B to ASH
2024-12-05 12:46:04,220 - moleculekit.tools.preparation - INFO - Modified residue CYS   140 B to CYX
2024-12-05 12:46:04,220 - moleculekit.tools.preparation - INFO - Modified residue HIS   171 B to HID
2024-12-05 12:46:04,221 - moleculekit.tools.preparation - INFO - Modified residue CYS   217 B to CYX
2024-12-05 12:46:04,221 - moleculekit.tools.preparation - INFO - Modified residue HIS   223 B to HID
2024-12-05 12:46:04,221 - moleculekit.tools.preparation - INFO - Modified residue HIS   297 B to HID
2024-12-05 12:46:04,222 - moleculekit.tools.preparation - INFO - Modified residue HIS   319 B to HIE
2024-12-05 12:46:04,223 - moleculekit.tools.preparation - WARNING - Dubious protonation state: the pKa of 1 residues is within 1.0 units of pH 7.4.
2024-12-05 12:46:04,224 - moleculekit.tools.preparation - WARNING - Dubious protonation state:    ASP   114 B (pKa= 7.85)
prot.view()
/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)
#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-12-05 12:46:07,287 - moleculekit.molecule - INFO - Removed 327 residues from appended Molecule due to collisions.
2024-12-05 12:46:07,493 - htmd.builder.solvate - INFO - Using water pdb file at: /home/sdoerr/Work/htmd/htmd/share/solvate/wat.pdb
2024-12-05 12:46:07,945 - htmd.builder.solvate - INFO - Replicating 8 water segments, 2 by 2 by 2

2024-12-05 12:46:16,364 - htmd.builder.solvate - INFO - 9644 water molecules were added to the system.
2024-12-05 12:46:18,858 - htmd.builder.amber - INFO - Detecting disulfide bonds.
2024-12-05 12:46:18,875 - 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-12-05 12:46:19,910 - htmd.builder.amber - INFO - Starting the build.
2024-12-05 12:46:22,046 - htmd.builder.amber - INFO - Finished building.
2024-12-05 12:46:23,017 - moleculekit.writers - WARNING - Field "resid" of PDB overflows. Your data will be truncated to 4 characters.
2024-12-05 12:46:28,148 - 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-12-05 12:46:28,452 - htmd.builder.ionize - INFO - Adding 15 anions + 0 cations for neutralizing and 68 ions for the given salt concentration 0.15 M.
2024-12-05 12:46:37,889 - htmd.builder.amber - INFO - Starting the build.
2024-12-05 12:46:39,968 - htmd.builder.amber - INFO - Finished building.
2024-12-05 12:46:40,942 - moleculekit.writers - WARNING - Field "resid" of PDB overflows. Your data will be truncated to 4 characters.
2024-12-05 12:46:45,901 - 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-12-05 12:46:45,936 - 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-12-05 12:46:48,396 - moleculekit.tools.sequencestructuralalignment - INFO - Alignment #0 was done on 282 residues: 2-285

Equilibrate#

from acemd.protocols import setup_equilibration

# Use a 10A flat bottom potential on the ligand coordinates to prevent the ligand from diffusing
# from its original position during equilibration.
# You can refer to https://software.acellera.com/acemd/manual.html#extforces-options
# for more information on the restraint options
restr = {
    "type": "positionalRestraint",
    "sel": "resname MOL and noh",
    "fbwidth": [10, 10, 10],
    "setpoints": ["5@0ns"],
}

setup_equilibration(
    os.path.join(path,'build'),
    os.path.join(path,'equil'),
    run="40ns",
    temperature=300,
    barostatconstratio=True,
    extforces=[restr]
)
2024-12-05 12:46:48,911 - rdkit - INFO - Enabling RDKit 2024.03.5 jupyter extensions
2024-12-05 12:46:48,993 - acemd - INFO - # You are on the latest ACEMD version (4.0.1).
2024-12-05 12:46:48,994 - acemd - INFO - Copying ./01_prepare/build/structure.pdb to ./01_prepare/equil/structure.pdb
2024-12-05 12:46:48,997 - acemd - INFO - Copying ./01_prepare/build/structure.prmtop to ./01_prepare/equil/structure.prmtop
2024-12-05 12:46:50,432 - acemd - WARNING - Found cis peptide bond with dihedral angle -3.32 deg in the omega diheral (HID 160 CA  ) (HID 160 C  ) (PRO 161 N  ) (PRO 161 CA  ) with indexes [2528 2541 2543 2553]
mdx = LocalGPUQueue()
mdx.submit(os.path.join(path, 'equil'))
mdx.wait()
2024-12-05 12:46:51,491 - jobqueues.util - INFO - Trying to determine all GPU devices
2024-12-05 12:46:51,547 - jobqueues.localqueue - INFO - Using GPU devices 0
2024-12-05 12:46:51,548 - jobqueues.util - INFO - Trying to determine all GPU devices
2024-12-05 12:46:51,585 - jobqueues.localqueue - INFO - Queueing /tmp/testmor/01_prepare/equil
2024-12-05 12:46:51,586 - jobqueues.localqueue - INFO - Running /tmp/testmor/01_prepare/equil on device 0
2024-12-05 12:48:19,471 - jobqueues.localqueue - INFO - Completed /tmp/testmor/01_prepare/equil

Production#

from acemd.protocols import setup_production

# Apply a flat bottom potential to prevent the ligand from changing from the
# extra- to the intra-cellular side of the membrane from the periodic image.
# Link the box to the center of mass of the membrane but offset the box
# center in the positive z direction by 30 Angstrom so that it's situated
# above the membrane. The box will have 60A width in the z direction.

restr = {
    "type": "positionalRestraint",
    "sel": "resname MOL and noh",
    "fbwidth": [70, 70, 60],
    "fbcenter": "lipids",
    "fbcenteroffset": [0, 0, 30],
    "setpoints": ["5@0ns"],
}

setup_production(
    os.path.join(path,'equil'),
    os.path.join(path,'prod'),
    run="50ns",
    temperature=300,
    extforces=[restr]
)
2024-12-05 12:50:16,763 - acemd - INFO - Copying ./01_prepare/equil/output.coor to ./01_prepare/prod/input.coor
2024-12-05 12:50:16,766 - acemd - INFO - Copying ./01_prepare/equil/output.xsc to ./01_prepare/prod/input.xsc
2024-12-05 12:50:16,767 - acemd - INFO - Copying ./01_prepare/equil/structure.prmtop to ./01_prepare/prod/structure.prmtop
2024-12-05 12:50:17,710 - acemd - WARNING - Found cis peptide bond with dihedral angle -10.44 deg in the omega diheral (HID 160 CA  ) (HID 160 C  ) (PRO 161 N  ) (PRO 161 CA  ) with indexes [2528 2541 2543 2553]
# Read in the last frame of the equilibration
mol = Molecule(os.path.join(path,'equil','structure.prmtop'))
mol.read(os.path.join(path,'equil','output.coor'))
mol.read(os.path.join(path,'equil','output.xsc'))
mol.view('not water')

# Visualize the FB box
width = np.array([70, 70, 60])
fbcenter = mol.getCenter("lipids", com=True)
fbcenteroffset = np.array([0, 0, 30])
b = VMDBox(np.vstack((fbcentre + fbcenteroffset - width/2, fbcentre + fbcenteroffset + width/2)).T.flatten())
2024-12-05 12:55:59,602 - moleculekit.writers - WARNING - Field "resid" of PDB overflows. Your data will be truncated to 4 characters.