I have been testing whether I could use Modeller to minimize a protein-protein complex.
As a first test, I've tried to perform a rigid body minimization of the bound structures of receptor and ligand.
Even though the starting structure has a Lig-RMSD of 3 A from the native complex structure (i.e., after superimposing the receptor subunits and measuring the CA RMSD between the ligand subunits), the final minimized structure does not get closer to the native complex but goes farther away, with Lig-RMSD of 12 A. I've tried with several near-native conformations and different protein complexes with similar results.
The script I am running is the following:
from modeller import * from modeller.scripts import complete_pdb from modeller.optimizers import conjugate_gradients, molecular_dynamics, actions
env = environ() env.io.atom_files_directory = ['../atom_files'] env.edat.dynamic_sphere = True
code = "test.pdb" mdl = complete_pdb(env, code) mdl.write(file=code+'.ini')
# Select rigid bodies # Receptor rigid body atmsel_a = mdl.chains['A'].atoms r_a = rigid_body(atmsel_a)
# Ligand rigid body atmsel_b = mdl.chains['B'].atoms r_b = rigid_body(atmsel_b)
# Select all atoms: atmsel = selection(mdl)
# Generate the restraints: mdl.restraints.rigid_bodies.append(r_a) mdl.restraints.rigid_bodies.append(r_b) mdl.restraints.make(atmsel, restraint_type='STEREO', spline_on_site=False) mdl.restraints.write(file=code+'.rsr')
mpdf = atmsel.energy()
# Create optimizer objects and set defaults for all further optimizations cg = conjugate_gradients(output='REPORT') md = molecular_dynamics(output='REPORT')
# Open a file to get basic stats on each optimization trcfil = open(code+'.D00000001', 'w')
# Run CG on the all-atom selection; write stats every 5 steps cg.optimize(atmsel, max_iterations=20, actions=actions.trace(5, trcfil))
# Run MD; write out a PDB structure (called '1fas.D9999xxxx.pdb') every # 10 steps during the run, and write stats every 10 steps md.optimize(atmsel, temperature=300, max_iterations=20, actions=[actions.write_structure(10, code+'.D9999%04d.pdb'), actions.trace(10, trcfil)])
# Finish off with some more CG, and write stats every 5 steps cg.optimize(atmsel, max_iterations=20, actions=[actions.trace(5, trcfil)])
mpdf = atmsel.energy()
I have also played with several minimization schemes, with and without MD, but the final minimized model does not get closer to the native complex.
Thanks in advance, Miguel
On 2/13/18 7:30 AM, email@example.com wrote: > I have been testing whether I could use Modeller to minimize a > protein-protein complex.
Modeller's strength is satisfying a bunch of distance restraints. It's not a great general purpose dynamics package - for that you'd be better off using something like CHARMM or GROMACS. Your script looks fine but you haven't given Modeller any restraints between the receptor and ligand, so there's nothing to pull them together. In a regular MD package you might expect solvent, dispersion, or electrostatics to do this, but Modeller doesn't model any of those interactions.
Ben Webb, Modeller Caretaker