[modeller_usage] Help with rigid body minimization of protein-protein complex

Hello,

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

env.libs.topology.read(file='$(LIB)/top_heav.lib')
env.libs.parameters.read(file='$(LIB)/par.lib')

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()

mdl.write(file=code+'.B')

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

http://bsc.es/disclaimer