Protein-ligand Binding free Energy calculation ????

[rijuaxt]

I was wondering if there any code or package available to calculate binding free energy of a ligand employing alchemical free energy calculation method or MM/PBSA method? There are many MD builds available for generating good trajectory but not analysing the data sufficiently. I haven't found any reliable opensource software to calculate MM/PBSA, MM/GBSA from single PDB file. If i find any help from here, I will be anticipating to employ the code on a series of PDB file using some loop functions.

[rijuaxt]

Hello, this question is yet to be answered?
Could anyone please respond?

[jchodera]

Are you looking for endpoint methods only, where you run simulations in explicit solvent and then strip off the solvent and reprocess with PB or GB? Or are you looking for alchemical free energy methods (like yank for absolute and perses for relative)?

[peastman]

Is this still an open question, or have you resolved it?

[peastman]

I'm closing this, since there's been no response since October.

[rafalbachorz]

Hello - I would like to reopen that and kindly ask for any hints related to the possibilities of getting MMPBSA/MMGBSA binding energies with openMM, preferably with AMBER trajectories. Thank in advance for any hint.

[swails]

Amber has a utility for doing that in MMPBSA.py.

It's reasonably straightforward to use mdtraj with OpenMM to loop through a trajectory and compute energies (less straightforward to do the decomposition that Amber can do). But there is no PB force or any plug-in I'm aware of so you'd be limited to GB.

[pguillem]

Hello @swails
Thanks for sharing.

Could you point to some examples to compute GB energies, looping the traj with mdtraj/openmm?
Thanks in advance
Pedro

[peastman]

I assume you've already created a System and Simulation in OpenMM, and that you've loaded the trajectory with MDTraj. You can compute the energy of each frame like this.

for i in range(trajectory.n_frames):
    simulation.context.setPositions(trajectory.xyz[i])
    energy = simulation.context.getState(getEnergy=True).getPotentialEnergy()
    # Do something with it...

[pguillem]

Thank you @peastman

So, if I use "atom_slice()" to narrow the trajectory.xyz() output to atoms of a protein-protein interface (a mask), and plug those coords to context.setPositions() ... then State.getPotentialEnergy() should return the isolated potential energy of such interface, right?.

For clarification:
I'm only studying energy fluctuations within a protein-protein interface, frame by frame.

I'm also trying to do the same on a per-residue basis... one by one.
Amber MMPBSA/MMGBSA already has this functionality+E.decomposition... but its annoyingly complicated to use.

[swails]

You can't "slice" the system in OpenMM -- you have to create a new one with just the subset of atoms you want. ParmEd has this kind of functionality, but OpenMM doesn't out-of-the-box.

[pguillem]

Thanks @swails!
I'll run straight for that approach.

I think per-residue energies (using atom masks) would be a super handy feature for structural bioinformaticians in OpenMM. I'll try to make it into a plugin and share the code (...still an openmm newbie here... maybe in a few months).

[mainguyenanhvu]

Do you have any python code for MMGBSA by using openMM?

Please help me. Thank you very much.

Besides, I found a sample code:

from simtk.openmm.app import *
from simtk.openmm import *
from simtk.unit import *
import os

# Load the PDB file of your protein-ligand complex.
pdb = PDBFile('complex.pdb')

# Create an OpenMM System object from the PDB file.
forcefield = ForceField('amber14-all.xml', 'amber14/tip3pfb.xml')
system = forcefield.createSystem(pdb.topology, nonbondedMethod=PME, nonbondedCutoff=1*nanometer)

# Add the MMPBSA force to the OpenMM System object.
mmpbsa = MMPBSAForce()
system.addForce(mmpbsa)

# Create an OpenMM Simulation object from the OpenMM System object.
integrator = LangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
simulation = Simulation(pdb.topology, system, integrator)

# Run the simulation and calculate the binding free energy.
simulation.context.setPositions(pdb.positions)
simulation.minimizeEnergy()
simulation.reporters.append(StateDataReporter('output.txt', 1000, step=True, potentialEnergy=True))
simulation.step(10000)

Is it true?

How to use GBSAOBCForce within python script?

[peastman]

See the user guide for instructions on how to use implicit solvent.