HSP90 with ADP:Mg2+ simulation¶
Simulating a protein:ligand complex with multisite solvated ions.
Files¶
All input files can be found in hsp90_adp_mg.zip
Introduction¶
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 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.
This example illustrates the use of Amber’s tleap to set up a simulation of the Heat Shock Protein 90 (HSP90, 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, Sept Lab), and use the Carlson parameters for ADP (DOI:10.1002/jcc.10262, downloaded from Bryce group Amber parameter database).
We begin from the 1BYQ PDB file, add missing residues (only those in the middle of the chain) and missing heavy atoms using PDBFixer. Using MDTraj, 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 records for bonds between Mg2+ and ADP are deleted. Crystallographic waters are preserved for coordination to the Mg2+.
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), ADP.prep and frcmod.phos files for the ADP parameters (downloaded from Bryce group Amber parameter database). 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).
[ ]:
!wget https://openmm.org/tutorials_/hsp90_adp_mg/files/hsp90_adp_mg.zip
!unzip -o hsp90_adp_mg.zip
[ ]:
from pdbfixer import PDBFixer
from openmm.app import PDBFile
import mdtraj as md
import os
import numpy as np
# clean up the original PDB file and add missing residues and heavy atoms
fixer = PDBFixer('pdb1byq.ent')
fixer.findMissingResidues()
# only add missing residues in the middle of the chain, do not add terminal ones
chains = list(fixer.topology.chains())
keys = fixer.missingResidues.keys()
missingResidues = dict()
for key in keys:
chain = chains[key[0]]
if not (key[1] == 0 or key[1] == len(list(chain.residues()))):
missingResidues[key] = fixer.missingResidues[key]
fixer.missingResidues = missingResidues
fixer.findMissingAtoms()
fixer.addMissingAtoms()
PDBFile.writeFile(fixer.topology, fixer.positions, open('1byq_fixed.pdb', 'w'))
traj = md.load('1byq_fixed.pdb')
# fix residue names
# in adp - change ` to * in atom names
# add Mg dummy atoms
for residue in traj.top.residues:
if residue.name == 'MG':
residue.name = 'MMG'
residue.atom(0).name = 'XZ'
for i in range(6):
traj.top.add_atom('D%s' % (i+1), md.element.Element.getBySymbol('Mg'), residue)
# add dummy atoms
central_index = residue.atom(0).index
central_positions = traj.xyz[0,central_index]
# positions from magnesium_ms.pdb, Sept Lab: http://septlab.engin.umich.edu/multisite-ions.html
dummy_positions = central_positions + [
[0.0, 0.0, 0.09],
[0.09, 0.0, 0.0],
[0.0, 0.0, -0.09],
[-0.09, 0.0, 0.0],
[0.0, -0.09, 0.0],
[0.0, 0.09, 0.0]
]
traj.xyz = np.array([np.insert(traj.xyz[0], central_index+1, dummy_positions, axis=0)])
elif residue.name == 'ADP':
residue.name = 'adp'
for atom in residue.atoms:
if atom.name[-1] == "'":
atom.name = atom.name[:-1] + "*"
# remove adp - mg CONECT bonds
bonds = []
for bond in traj.top._bonds:
if not (bond[0].residue.name == 'MMG' or bond[1].residue.name == 'MMG'):
bonds.append(bond)
traj.top._bonds = bonds
traj.save('1byq_fixed2.pdb')
# save the tleap script to file
with open('leaprc.hsp90', 'w') as f:
f.write('''
source oldff/leaprc.ff99SBildn
loadOff mag.lib
loadamberparams frcmod_mg.mmg
loadAmberPrep ADP.prep
loadamberparams frcmod.phos
x = loadPdb 1byq_fixed2.pdb
addIons x Na+ 0
solvateBox x TIP3PBOX 10.0
savePdb x topology.pdb
saveAmberParm x input.prmtop input.inpcrd
quit
''')
# run tleap
# if you have amber tools installed on your system you can uncomment and run this line
#os.system('tleap -f leaprc.hsp90')
# we have included the output files that get produced:
# input.inpcrd and input.prmtop
#!echo input.prmtop | python zeroBvalues.py
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.
[ ]:
from openmm import app
import openmm as mm
from openmm import unit
from sys import stdout
# load in Amber input files
prmtop = app.AmberPrmtopFile('input.prmtop.mod')
inpcrd = app.AmberInpcrdFile('input.inpcrd')
# prepare system and integrator
system = prmtop.createSystem(nonbondedMethod=app.PME,
nonbondedCutoff=1.0*unit.nanometers, constraints=app.HBonds, rigidWater=True,
ewaldErrorTolerance=0.0005)
integrator = mm.LangevinIntegrator(300*unit.kelvin, 1.0/unit.picoseconds,
2.0*unit.femtoseconds)
integrator.setConstraintTolerance(0.00001)
# prepare simulation
simulation = app.Simulation(prmtop.topology, system, integrator)
simulation.context.setPositions(inpcrd.positions)
# minimize
print('Minimizing...')
simulation.minimizeEnergy()
# equilibrate for 100 steps
simulation.context.setVelocitiesToTemperature(300*unit.kelvin)
print('Equilibrating...')
simulation.step(100)
# append reporters
simulation.reporters.append(app.DCDReporter('trajectory.dcd', 1000))
simulation.reporters.append(app.StateDataReporter(stdout, 1000, step=True,
potentialEnergy=True, temperature=True, progress=True, remainingTime=True,
speed=True, totalSteps=25000000, separator='\t'))
# run 50 ns of production simulation
print('Running Production...')
simulation.step(25000000)
print('Done!')
