Setting up a kinase:inhibitor constant pH simulation

Setting up a kinase simulation

Please post any errors you might have on the constph slack channel so I can look into them when I get back.

Note: You may want to try this without a ligand first. While the code should support it in theory, the ligand inclusion from a structure perspective is in no way streamlined yet. I'll add instructions at the end on how to give that a try.

Simulating a protein with constant-pH using protons

The most complete current example of setting up a protein system using the OpenMM ForceField class can be found in tests/test_ff.py in the test_create_peptide_simulation_using_protons_xml test function, but the required steps are detailed below.

Retrieve a structure

For a kinase, you can grab a kinase structure from the pdbfinder repo.

Rename residues to their titratable forms

This currently must be done manually.

Rename the following residues:

• all histidines to HIP
• all ASP to AS4
• all GLU to GL4.

Model missing atoms

Use app.modeler to add any missing atoms, especially the hydrogens present in the titratable residue templates (HIP, AS4, GL4).

To help you, paths to useful OpenMM format xml files are provided as variables in the protons.ff module:

• protonsff
• hydrogens
• bonds
• ions
• gaff

Be sure to load the provided bond definitions into topology before reading a pdb/pdbx file to be able to recognize residues correctly:

from protons import ff
app.Topology.loadBondDefinitions(ff.bonds)

You will also need to load the hydrogen definitions into app.Modeller before adding hydrogens; use can feed the protons.ff.hydrogens file path to Modeller.loadHydrogenDefinitions() for this:

from protons import ff
app.Modeller.loadHydrogenDefinitions(ff.hydrogens)

Solvate

Use app.Modeller to add solvent, and include some salt. You can use the standard water models included with OpenMM (tip4pew, tip3p, spc).

The ion models are provided through the ff module (e.g. protons.ff.ions_tip3p).

Minimize the system

Simulate using the ForceFieldProtonDrive driver

Rough example:

    from openmmtools.integrators import BAOABIntegrator

    integrator = BAOABIntegrator(temperature, friction, dt)
    integrator.setConstraintTolerance(constraintTolerance)

    driver = ForceFieldProtonDrive(
        system, temperature, 7.0, [ff.protonsff],
        topology, integrator, debug=False, pressure=pressure,
        ncmc_steps_per_trial=1, implicit=False, cationName="NA",
        anionName="CL")

    simulation = app.Simulation(
        topology,
        system,
        driver.compound_integrator,
        platform)
    simulation.context.setPositions(positions)
    simulation.context.setVelocitiesToTemperature(temperature)

    # Run simulation by alternating between integration and protonation state update
    for iteration in range(niterations):
        simulation.step(nsteps_per_iteration)
        driver.update(simulation.context, nattempts=nattempts_per_iteration)

Ligands

The protons.ligands module provides tools to parametrize a ligand.

This requires schrodinger tools (for Epik) and openeye toolkit (for small molecule manipulation, parameterization, and charging). We are working to reduce the number of non-open-source external requirements for ligand simulation, but this is hindered by the lack of methods for protonation state prediction.

Use the protons.ligands.generate_protons_ffxml() function to generate an enhanced ffxml file containing ligand parameters for all protonation states, as well as protonation state penalties for the pH of interest.

It is recommended that you use a Maestro mae file as input, but a mol2 file should also work. Make sure that all ligand atoms have unique names, and match the atom names in the pdb file you will use to simulate (or the corresponding Topology object). You will also need to make sure that all possible protons are included in the input file, so this file should represent the fully-protonated form of the small molecule. You can run epik once to find out which ones need to be there.

(Question: Does the PDB file need to also contain the fully-protonated form of the ligand, or will missing protons be added automatically by some tool?)

(Question: How do you currently generate the fully-protonated ligand input form?)

Note that it might be easiest to add these atoms to the pdb file first, and subsequently delete the ones that are no longer necessary after you've seen the output from Epik/protons.ligands.

If you run the prepared ligand structure, it should give you the location of the ffxml file that was generated. Double-check it for errors and please report when you see them.

(Question: What is meant by "run" above?)

If you prepared the same atom names in the pdb file, it should be relatively straightforward from here. Provide the forcefield template to the app.ForceField class, and it should recognize your ligand. Proceed as you would with a protein simulation.

I recommend you try imatinib-abl first. It's been parametrized successfully before.

(Question: Why is this not in protons/examples/?)