Hi, I am trying to model a protein with certain distance restraints and disulphide bonds. I wonder why restraint does not work when I set disulphide bonds. My script looks like.
--------------------------------------------- from modeller import * from modeller.automodel import * # Load the automodel class import sys
log.verbose() env = environ() env.io.hetatm = True env.schedule_scale = physical.values(default=1.0, soft_sphere=0.7)
# directories for input atom files env.io.atom_files_directory = ['.', '../atom_files']
class MyModel(automodel): def special_restraints(self, aln): rsr = self.restraints at = self.atoms
# Residues to be an alpha helix: rsr.add(secondary_structure.alpha(self.residue_range('1:', '22:'))) rsr.add(secondary_structure.alpha(self.residue_range('26:', '32:')))
rsr.add(forms.gaussian(group=physical.xy_distance, feature=features.distance(at['CA:1'], at['CA:269']), mean=4.0, stdev=4.0))
def special_patches(self, aln): # A disulfide between residues X and Y: self.patch(residue_type='DISU', residues=(self.residues['113'], self.residues['117'])) self.patch(residue_type='DISU', residues=(self.residues['125'], self.residues['175'])) self.patch(residue_type='DISU', residues=(self.residues['313'], self.residues['337'])) self.patch(residue_type='DISU', residues=(self.residues['24'],
self.residues['314']))a = MyModel(env, alnfile = 'seq.ali', # alignment filename knowns = ('XYZ','ABC'), # codes of the templates sequence = 'prot', assess_methods=(assess.DOPE, assess.GA341)) # code of the target a.starting_model= 1 # index of the first model a.ending_model = 20 # index of the last model # (determines how many models to calculate)
a.library_schedule = autosched.fast # Very thorough VTFM optimization a.max_var_iterations = 100 a.md_level = refine.fast # Thorough MD optimization a.repeat_optimization = 1 # Repeat the whole cycle 2 times and.. a.max_molpdf = 1e6 # do not stop unless obj.func. > 1E6
#a.use_parallel_job(j) # Use the job for model building a.make() -------------------------------------------------------
If I comment distance restraint seems to work okay. What can be the reason that it does not work with S-S restraints. Moreover, Is there a way to control force constant for harmonic restraints in such cases ?
Best, SM
On 10/3/16 1:47 AM, Sushil Mishra wrote: > I am trying to model a protein with certain distance restraints and > disulphide bonds. I wonder why restraint does not work when I set > disulphide bonds.
What do you mean by "does not work"? Your script looks fine to me.
> What can be the > reason that it does not work with S-S restraints. Moreover, Is there a > way to control force constant for harmonic restraints in such cases ?
Technically, you are not adding an S-S restraint - you are telling Modeller that the two disulfides are bonded, so when it turns on the CHARMM forcefield you'll get a term in the scoring function that corresponds to the S-S bond. So if you want to adjust that force constant (I can't imagine why you'd want to) the way to do it would be to edit the CHARMM parameter file, modlib/par.lib (in the BOND section, which is arranged by CHARMM atom type; the two S atoms are type SM, so you'd be looking around line 291 in Modeller 9.17).
Ben Webb, Modeller Caretaker
On Tue, Oct 4, 2016 at 3:43 AM, Modeller Caretaker modeller-care@salilab.org wrote: > On 10/3/16 1:47 AM, Sushil Mishra wrote: >> >> I am trying to model a protein with certain distance restraints and >> disulphide bonds. I wonder why restraint does not work when I set >> disulphide bonds. > > > What do you mean by "does not work"? Your script looks fine to me.
I mean, in the generated models those restrained atoms (CA of 1 and 169) are far away (>20A) while distance should be somewhere between within 8A. The CA atom of residues 1 and 169 are close enough in the model when function special_patches is commented. However, I am missing a S-S bond without this function. I do not understand why those two CA atoms are >20 in the model when there is restrain set between them:
rsr.add(forms.gaussian(group=physical.xy_distance, feature=features.distance(at['CA:1'], at['CA:269']), mean=4.0, stdev=4.0))
>> What can be the >> reason that it does not work with S-S restraints. Moreover, Is there a >> way to control force constant for harmonic restraints in such cases ? > > > Technically, you are not adding an S-S restraint - you are telling Modeller > that the two disulfides are bonded, so when it turns on the CHARMM > forcefield you'll get a term in the scoring function that corresponds to the > S-S bond. So if you want to adjust that force constant (I can't imagine why > you'd want to) the way to do it would be to edit the CHARMM parameter file, > modlib/par.lib (in the BOND section, which is arranged by CHARMM atom type; > the two S atoms are type SM, so you'd be looking around line 291 in Modeller > 9.17).
I understand it now. I guessed it incorrectly. I thought MODELLER may be using this information to restraint both Cys residues during the modeling so that they are close enough to form S-S bond.
Thanks, SM
> > Ben Webb, Modeller Caretaker > -- > modeller-care@salilab.org https://salilab.org/modeller/ > Modeller mail list: https://salilab.org/mailman/listinfo/modeller_usage
On 10/03/2016 04:38 PM, Sushil Mishra wrote: > On Tue, Oct 4, 2016 at 3:43 AM, Modeller Caretaker > modeller-care@salilab.org wrote: >> What do you mean by "does not work"? Your script looks fine to me. > > I mean, in the generated models those restrained atoms (CA of 1 and > 169) are far away (>20A) while distance should be somewhere between > within 8A.
In that case, most likely either the restraint cannot be satisfied at the same time as satisfying both the S-S bond and your alignment-derived (homology) restraints. You can certainly try building a larger ensemble of models to see if you can sample a conformation that does satisfy all restraints, but more likely you'll have to relax one or more of your restraints (e.g. modify the alignment so that the model can move further from the template).
Ben Webb, Modeller Caretaker
-
Modeller Caretaker -
Sushil Mishra