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
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:
Awas not recognized
Bwas 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) 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 [ ]: