Using SMIRNOFF with Amber on BRD4:inhibitor complexes: Exporting parameterized complexes to Amber, Gromacs, and CHARMM

This example applies SMIRNOFF-format parameters to BRD4 inhibitors from the living review on binding free energy benchmark systems by Mobley and Gilson. The BRD4 system comes from the accompanying GitHub repository.

This example uses ParmEd to combine a protein parameterized with Amber with a ligand parameterized with SMIRNOFF. This example is meant to illustrate how to apply parameters to a single ligand, but it's also easy to process many ligands.



In [1]:

    
# Retrieve protein and ligand files for BRD4 and a docked inhibitor from the benchmark systems GitHub repository
# https://github.com/MobleyLab/benchmarksets
import requests
repo_url = 'https://raw.githubusercontent.com/MobleyLab/benchmarksets/master/input_files/'
sources = {
    'receptor.pdb' : repo_url + 'BRD4/pdb/BRD4.pdb',
    'ligand.pdb'   : repo_url + 'BRD4/pdb/ligand-1.pdb',
    'ligand.sdf'   : repo_url + 'BRD4/sdf/ligand-1.sdf',
}
for (filename, url) in sources.items():
    r = requests.get(url)
    open(filename, 'w').write(r.text)

Parametrize a molecule with the SMIRNOFF-format "Parsley" force field

First, we parametrize the ligand with the Parsley force field.

We use both a PDB file and an SDF file for the ligand.



In [2]:

    
# Create an OpenFF Molecule object from the ligand SDF fiel
from openforcefield.topology import Molecule
ligand_off_molecule = Molecule('ligand.sdf')

# Load the SMIRNOFF-format Parsley force field
from openforcefield.typing.engines.smirnoff import ForceField
force_field = ForceField('openff_unconstrained-1.0.0.offxml')

# Parametrize the ligand molecule by creating a Topology object from it
ligand_system = force_field.create_openmm_system(ligand_off_molecule.to_topology())

and we convert the OpenMM System to a ParmEd Structure that we'll be able to mix with the protein

Warning: ParmEd's Structure model is inspired by AMBER. Some information in an OpenMM System are not directly translatable into a Structure. In particular, long-range interaction treatment method (e.g., PME, CutoffPeriodic) and parameters (e.g., cutoff and cutoff switching distance, PME error tolerance) are known to be lost during the conversion.



In [3]:

    
# Read in the coordinates of the ligand from the PDB file
from simtk.openmm.app import PDBFile
ligand_pdbfile = PDBFile('ligand.pdb')

# Convert OpenMM System object containing ligand parameters into a ParmEd Structure.
import parmed
ligand_structure = parmed.openmm.load_topology(ligand_pdbfile.topology,
                                                ligand_system,
                                                xyz=ligand_pdbfile.positions)

Create a ParmEd `Structure` of an AMBER-parametrized receptor

We have to create a ParmEd Structure containing positions and parameters of the receptor (BRD4) that we will then combine with the positions and parameters in the ligand Structure we created above.

First, we use OpenMM to assign Amber14 parameters using OpenMM, but you can also create a Structure from Amber prmtop and inpcrd files.

Note: If you already have AMBER (prmtop/inpcrd), GROMACS (top/gro), or any other file specifying the protein parameters supported by ParmEd, you can simply load the files directly into a Structure using ParmEd's functionalities. See https://parmed.github.io/ParmEd/html/readwrite.html .



In [4]:

    
receptor_pdbfile = PDBFile('receptor.pdb')

# Load the AMBER protein force field through OpenMM.
from simtk.openmm import app
omm_forcefield = app.ForceField('amber14-all.xml')

# Parameterize the protein.
receptor_system = omm_forcefield.createSystem(receptor_pdbfile.topology)

# Convert the protein System into a ParmEd Structure.
receptor_structure = parmed.openmm.load_topology(receptor_pdbfile.topology,
                                                 receptor_system,
                                                 xyz=receptor_pdbfile.positions)

Combine receptor and ligand structures

We can then merge the receptor and ligand Structure objects into a single Structure containing both positions and parameters.



In [5]:

    
complex_structure = receptor_structure + ligand_structure

Export to OpenMM

Once we have the Structure for the complex, containing both positions and parameters, we can chose to create an OpenMM System object that we can simulate using OpenMM:



In [6]:

    
from simtk.openmm.app import NoCutoff, HBonds
from simtk import unit

# Convert the Structure to an OpenMM System in vacuum.
complex_system = complex_structure.createSystem(nonbondedMethod=NoCutoff,
                                                nonbondedCutoff=9.0*unit.angstrom,
                                                constraints=HBonds,
                                                removeCMMotion=False)



In [7]:

    
# Export the System to an OpenMM System XML and PDB file
complex_structure.save('complex.pdb', overwrite=True)
from simtk.openmm import XmlSerializer
with open('complex.xml', 'w') as f:
    f.write(XmlSerializer.serialize(complex_system))

Export to Amber

We can also export the Structure to Amber prmtop and inpcrd files using ParmEd:



In [8]:

    
# Export the Structure to AMBER files
complex_structure.save('complex.prmtop', overwrite=True)
complex_structure.save('complex.inpcrd', overwrite=True)

Export to gromacs

We can export the Structure to gromacs gro and top files using ParmEd:



In [9]:

    
# Export the Structure to Gromacs files
complex_structure.save('complex.gro', overwrite=True)
complex_structure.save('complex.top', overwrite=True)