About Issues Tutorials Documentation
</span>
This notebook prepares a co-crystallized protein / small molecule ligand structure from the PDB database and prepares it for molecular dynamics simulation.
In [ ]:
import moldesign as mdt
import moldesign.units as u
First, we'll download and investigate the 3AID crystal structure.
In [ ]:
protease = mdt.from_pdb('3AID')
protease
In [ ]:
protease.draw()
This structure is not ready for MD - this command will raise a ParameterizationError
Exception. After running this calculation, click on the Errors/Warnings tab to see why.
In [ ]:
amber_ff = mdt.forcefields.DefaultAmber()
newmol = amber_ff.create_prepped_molecule(protease)
You should see 3 errors:
ARQ
not recognizedHD1
in residue HIS69
, chain A
was not recognizedHD1
in residue HIS69
, chain B
was not recognized(There's also a warning about bond distances, but these can be generally be fixed with an energy minimization before running dynamics)
We'll start by tackling the small molecule "ARQ".
In [ ]:
sel = mdt.widgets.ResidueSelector(protease)
sel
In [ ]:
drugres = mdt.Molecule(sel.selected_residues[0])
drugres.draw2d(width=700, show_hydrogens=True)
In [ ]:
drugmol = mdt.tools.set_hybridization_and_saturate(drugres)
drugmol.draw(width=500)
In [ ]:
drugmol.draw2d(width=700, show_hydrogens=True)
In [ ]:
drug_parameters = mdt.create_ff_parameters(drugmol, charges='gasteiger')
Section II. dealt with getting forcefield parameters for an unknown small molecule. Next, we'll prep the other part of the structure.
Waters in crystal structures are usually stripped from a simulation as artifacts of the crystallization process. Here, we'll remove the waters from the protein structure.
In [ ]:
dehydrated = mdt.Molecule([atom for atom in protease.atoms if atom.residue.type != 'water'])
Histidine is notoriously tricky, because it exists in no less than three different protonation states at biological pH (7.4) - the "delta-protonated" form, referred to with residue name HID
; the "epsilon-protonated" form aka HIE
; and the doubly-protonated form HIP
, which has a +1 charge. Unfortunately, crystallography isn't usually able to resolve the difference between these three.
Luckily, these histidines are pretty far from the ligand binding site, so their protonation is unlikely to affect the dynamics. We'll therefore use the guess_histidine_states
function to assign a reasonable starting guess.
In [ ]:
mdt.guess_histidine_states(dehydrated)
In [ ]:
amber_ff = mdt.forcefields.DefaultAmber()
amber_ff.add_ff(drug_parameters)
sim_mol = amber_ff.create_prepped_molecule(dehydrated)
In [ ]:
sim_mol.set_energy_model(mdt.models.OpenMMPotential, implicit_solvent='obc', cutoff=8.0*u.angstrom)
sim_mol.set_integrator(mdt.integrators.OpenMMLangevin, timestep=2.0*u.fs)
sim_mol.configure_methods()
In [ ]:
mintraj = sim_mol.minimize()
mintraj.draw()
In [ ]:
traj = sim_mol.run(20*u.ps)
In [ ]:
viewer = traj.draw(display=True)
viewer.autostyle()
In [ ]: