{ "cells": [ { "cell_type": "markdown", "id": "9afbaa7d", "metadata": {}, "source": [ "## HSP90 with ADP:Mg2+ simulation\n", "\n", "*Simulating a protein:ligand complex with multisite solvated ions.*\n", "\n", "### Files\n", " - All input files can be found in [hsp90_adp_mg.zip](https://openmm.org/tutorials_/hsp90_adp_mg/files/hsp90_adp_mg.zip)\n", "\n", "\n", "### Introduction\n", "\n", "OpenMM allows users to model their systems using Amber and provide `prmtop`, and `inpcrd` files as input. This allows users more familiar with the [Amber](http://ambermd.org/) modeling environment to continue using their setup tools, while harnessing the speed and versatility of OpenMM. This also allows the use of non-standard force fields that have been published for use with Amber, such as metal ion models where dummy atoms are applied to mimic particular coordination geometries. Furthermore, this is facilitated by OpenMM's support for Extra Particles---particles that are not ordinary atoms, such as these metal dummy atoms, virtual sites in many water models, etc.\n", "\n", "\n", "This example illustrates the use of Amber's tleap to set up a simulation of the Heat Shock Protein 90 (HSP90, Uniprot: [P07900](http://www.uniprot.org/uniprot/P07900), including the ADP-Mg2+ complex, where we will model the Mg2+ cation using the Sept Lab octahedral multisite model ([DOI:10.1021/ct400177g](https://doi.org/10.1021/ct400177g), [Sept Lab](http://septlab.engin.umich.edu/multisite-ions.html)), and use the Carlson parameters for ADP ([DOI:10.1002/jcc.10262](https://doi.org/10.1002/jcc.10262), downloaded from [Bryce group Amber parameter database](http://research.bmh.manchester.ac.uk/bryce/amber)).\n", "\n", "We begin from the [1BYQ](https://www.rcsb.org/pdb/explore.do?structureId=1BYQ) PDB file, add missing residues (only those in the middle of the chain) and missing heavy atoms using PDBFixer. Using [MDTraj](http://mdtraj.org/), residue and atom naming of the Mg2+ ion and the ADP are fixed to match those in the parameter files, and dummy atoms of the multisite Mg2+ are added. [CONECT](http://www.wwpdb.org/documentation/file-format-content/format33/sect10.html#CONECT) records for bonds between Mg2+ and ADP are deleted. Crystallographic waters are preserved for coordination to the Mg2+.\n", "\n", "\n", "\n", "\n" ] }, { "cell_type": "markdown", "id": "60d7cfce", "metadata": {}, "source": [ "Finally, `tleap` is run to add hydrogens and parametrize. We use the `ff99SBildn` force field, mag.lib and frcmod_mg.mmg files from the multisite Mg2+ model (downloaded from [Sept Lab](http://septlab.engin.umich.edu/multisite-ions.html)), ADP.prep and frcmod.phos files for the ADP parameters (downloaded from [Bryce group Amber parameter database](http://research.bmh.manchester.ac.uk/bryce/amber)). The prmtop and inpcrd files are saved for simulation in OpenMM. The multisite Mg2+ model requires a correction to the Lennard-Jones B-coefficients in the prmtop file, this is done by a short Python script (obtained by email from Prof. David Sept, [Sept Lab](http://septlab.engin.umich.edu/))." ] }, { "cell_type": "code", "execution_count": null, "id": "68bd18e5", "metadata": {}, "outputs": [], "source": [ "!wget https://openmm.org/tutorials_/hsp90_adp_mg/files/hsp90_adp_mg.zip\n", "!unzip -o hsp90_adp_mg.zip" ] }, { "cell_type": "code", "execution_count": null, "id": "e830a81a", "metadata": {}, "outputs": [], "source": [ "from pdbfixer import PDBFixer\n", "from openmm.app import PDBFile\n", "import mdtraj as md\n", "import os\n", "import numpy as np\n", "\n", "# clean up the original PDB file and add missing residues and heavy atoms\n", "fixer = PDBFixer('pdb1byq.ent')\n", "fixer.findMissingResidues()\n", "\n", "# only add missing residues in the middle of the chain, do not add terminal ones\n", "chains = list(fixer.topology.chains())\n", "keys = fixer.missingResidues.keys()\n", "missingResidues = dict()\n", "for key in keys:\n", " chain = chains[key[0]]\n", " if not (key[1] == 0 or key[1] == len(list(chain.residues()))):\n", " missingResidues[key] = fixer.missingResidues[key]\n", "fixer.missingResidues = missingResidues\n", "\n", "fixer.findMissingAtoms()\n", "fixer.addMissingAtoms()\n", "\n", "PDBFile.writeFile(fixer.topology, fixer.positions, open('1byq_fixed.pdb', 'w'))\n", "\n", "traj = md.load('1byq_fixed.pdb')\n", "\n", "# fix residue names\n", "# in adp - change ` to * in atom names\n", "# add Mg dummy atoms\n", "for residue in traj.top.residues:\n", " if residue.name == 'MG':\n", " residue.name = 'MMG'\n", " residue.atom(0).name = 'XZ'\n", " for i in range(6):\n", " traj.top.add_atom('D%s' % (i+1), md.element.Element.getBySymbol('Mg'), residue)\n", " # add dummy atoms\n", " central_index = residue.atom(0).index\n", " central_positions = traj.xyz[0,central_index]\n", " # positions from magnesium_ms.pdb, Sept Lab: http://septlab.engin.umich.edu/multisite-ions.html\n", " dummy_positions = central_positions + [\n", " [0.0, 0.0, 0.09],\n", " [0.09, 0.0, 0.0],\n", " [0.0, 0.0, -0.09],\n", " [-0.09, 0.0, 0.0],\n", " [0.0, -0.09, 0.0],\n", " [0.0, 0.09, 0.0]\n", " ] \n", " traj.xyz = np.array([np.insert(traj.xyz[0], central_index+1, dummy_positions, axis=0)]) \n", " elif residue.name == 'ADP':\n", " residue.name = 'adp'\n", " for atom in residue.atoms:\n", " if atom.name[-1] == \"'\":\n", " atom.name = atom.name[:-1] + \"*\"\n", "\n", "# remove adp - mg CONECT bonds\n", "bonds = []\n", "\n", "for bond in traj.top._bonds:\n", " if not (bond[0].residue.name == 'MMG' or bond[1].residue.name == 'MMG'):\n", " bonds.append(bond)\n", "\n", "traj.top._bonds = bonds\n", "\n", "traj.save('1byq_fixed2.pdb')\n", "\n", "# save the tleap script to file\n", "with open('leaprc.hsp90', 'w') as f:\n", " f.write('''\n", "source oldff/leaprc.ff99SBildn\n", "loadOff mag.lib\n", "loadamberparams frcmod_mg.mmg\n", "loadAmberPrep ADP.prep\n", "loadamberparams frcmod.phos\n", "x = loadPdb 1byq_fixed2.pdb\n", "addIons x Na+ 0\n", "solvateBox x TIP3PBOX 10.0\n", "savePdb x topology.pdb\n", "saveAmberParm x input.prmtop input.inpcrd\n", "quit\n", "''')\n", "\n", "# run tleap\n", "# if you have amber tools installed on your system you can uncomment and run this line\n", "#os.system('tleap -f leaprc.hsp90')\n", "\n", "# we have included the output files that get produced: \n", "# input.inpcrd and input.prmtop\n", "\n", "#!echo input.prmtop | python zeroBvalues.py\n" ] }, { "cell_type": "markdown", "id": "64da9756", "metadata": {}, "source": [ "We load the `prmtop` and `inpcrd` files in, by creating `AmberPrmtopFile` and `AmberInpcrdFile` objects. Next, the `System` is created by calling the `createSystem` method on the `AmberPrmtopFile` object. Next, the `LangevinIntegrator` and the `Simulation` are set up, using the topology from the `AmberPrmtopFile` and positions from the `AmberInpcrdFile`. The simulation is energy minimized and equilibrated for 100 steps. Reporters are attached and the production simulation propagated for 50 ns." ] }, { "cell_type": "code", "execution_count": null, "id": "e3edf284", "metadata": {}, "outputs": [], "source": [ "from openmm import app\n", "import openmm as mm\n", "from openmm import unit\n", "from sys import stdout\n", "\n", "# load in Amber input files\n", "prmtop = app.AmberPrmtopFile('input.prmtop.mod')\n", "inpcrd = app.AmberInpcrdFile('input.inpcrd')\n", "\n", "# prepare system and integrator\n", "system = prmtop.createSystem(nonbondedMethod=app.PME,\n", " nonbondedCutoff=1.0*unit.nanometers, constraints=app.HBonds, rigidWater=True,\n", " ewaldErrorTolerance=0.0005)\n", "integrator = mm.LangevinIntegrator(300*unit.kelvin, 1.0/unit.picoseconds,\n", " 2.0*unit.femtoseconds)\n", "integrator.setConstraintTolerance(0.00001)\n", "\n", "# prepare simulation\n", "simulation = app.Simulation(prmtop.topology, system, integrator)\n", "simulation.context.setPositions(inpcrd.positions)\n", "\n", "# minimize\n", "print('Minimizing...')\n", "simulation.minimizeEnergy()\n", "\n", "# equilibrate for 100 steps\n", "simulation.context.setVelocitiesToTemperature(300*unit.kelvin)\n", "print('Equilibrating...')\n", "simulation.step(100)\n", "\n", "# append reporters\n", "simulation.reporters.append(app.DCDReporter('trajectory.dcd', 1000))\n", "simulation.reporters.append(app.StateDataReporter(stdout, 1000, step=True,\n", " potentialEnergy=True, temperature=True, progress=True, remainingTime=True,\n", " speed=True, totalSteps=25000000, separator='\\t'))\n", "\n", "# run 50 ns of production simulation\n", "print('Running Production...')\n", "simulation.step(25000000)\n", "print('Done!')" ] } ], "metadata": { "tags": [ "tutorial" ], "kernelspec": { "display_name": "Python 3.9.13 ('openmm')", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.13" }, "nbsphinx": { "execute": "never" }, "vscode": { "interpreter": { "hash": "16b2d2c1789d035bceb6d775bd7ffc39b805c8f0529038638d98b11c7a85ade5" } } }, "nbformat": 4, "nbformat_minor": 5 }