Distance restraints applied to multiple proteins

User discussions
1

GROMACS version: 2020, 2018
GROMACS modification: Yes/No
Here post your question

Dear all,

I am trying to simulate a system with 4 identical proteins and restrain the C-Alpha of each of them (there are also some small molecules in the box too), so that the proteins wouldn’t unfold during the thermal annealing stage. I found few discussions on the subject (i.e. https://mailman-1.sys.kth.se/pipermail/gromacs.org_gmx-users/2015-July/099219.html), and tested it on small molecules, and it worked just fine. Nevertheless, for proteins the very same approach doesn’t seem to work.

Say, I have the following system:

#include “Protein1.itp”
#include “Protein2.itp”; exactly the same as Protein1, except the molecule type name
#include “Protein3.itp”; exactly the same as Protein1, except the molecule type name
#include “Protein4.itp”; exactly the same as Protein1, except the molecule type name

#include “API1.itp”
#include “API2.itp”; exactly the same as API1, except the molecule type name

[ system ]
; Name
Protein with API

[ molecules ]
; Compound #mols
Protein1 1
Protein2 1
Protein3 1
Protein4 1
API1 5
API2 5

If I make some 5 APIs restrained (and specify it in the API.itp [ position_restraints ] section) and keep the original topology of the other 5, everything works well. Also, if I restrain only the Protein1 it works also fine. However, if I try to restrain Protein1 and Protein2 (either as two different types of files Protein1.itp + restr_Protein1.itp and Protein2.itp + restr_Protein2.itp, or as “Protein1 2”) it either explodes with segfault, or says:
Assertion failed:
Condition: type_max - type_min + 1 == dd->nres
All distance restraint parameter entries in the topology should be consecutive.

As it works fine for smaller molecules and even 1 restrained Protein plus several restrained APIs, I assume it has something to do with the numbering of the residues for the protein (which has 153 amino acids).

I would be extremely grateful if someone could help with a hint on how to set up intra-protein distance restraints properly for all 4 molecules.

Thank you!
Aleksei

2

Hi Alex, I am having a similar problem, did you manage to resolve it in the meantime?

3

Hi srdjan,

I sort of found an alternative solution. Instead of the actual restraints (in a separate section of the .itp) I added inter-protein bonds with type 10. You can look it up here:
https://manual.gromacs.org/documentation/current/reference-manual/topologies/topology-file-formats.html
To make a list of the pairs I used the gmx genrestr applied to C-Alpha of a single protein. Then just copied it (and probably edited manually…) and inserted into the [bonds] section of the protein topology file.

Hope it helps!

Best,
Aleksei

4

Hi Alex,

Thanks a lot for your help, I am trying it now.

May I just ask one more questions, so you have created file by “gmx genrestr -disre” option which gives you something like this:

1 4 1 0 1 0.197768 0.397768 1.39777 1
1 5 1 1 1 0.281579 0.481579 1.48158 1
1 10 1 2 1 0.618802 0.818802 1.8188 1
1 11 1 3 1 0.499951 0.699951 1.69995 1

Before copying it into Bonds section of topology file, I need to change to bond type 10:

1 4 1 0 10 0.197768 0.397768 1.39777 1
1 5 1 1 10 0.281579 0.481579 1.48158 1
1 10 1 2 10 0.618802 0.818802 1.8188 1
1 11 1 3 10 0.499951 0.699951 1.69995 1

After that do you copy all these lines into bonds section or you need to omit some section, I see that bond section should look like this:
bonds ]
; ai aj funct c0 c1 c2 c3
1 2 1
2 3 1
2 4 1
4 5 1

5

The formats of the two bond types are different, but that is not an issue.
For instance, the “type 10” binds will look like:
5 24 10 0.276033 0.476033 1.47603 1000
, whereas the normal harmonic bonds are
1258 1259 1
However you need to throw away the labels (columns 3 and 4 in your example). You can do that with awk or something similar. Editing spaces to treat only numbers as columns might take some time, but the rest should go pretty smoothly.

6

Hi, thank you for help.

Yes, in the meantime I have supposed that columns 3 and 4 should disappear (I did it in Excel but it is smarter to write some python code for that).

But when I tried to run simulation I got an error: invalid pair type 10.

Which gromacs are you using? Or maybe I have copied them in wrong part of the topology file? I have copied it in “bonds” section, at the end after the section with “normal” bonds of type 1 have finished.

Where exactly in topology file have you copied it?

7

Are gromacs topology files space sensitive (like PDB files) and needs to be exact spaces between each columns etc?

8

I have it in the beginning of the [bonds] section. I think type 10 was introduced long ago, but I’m using 2018.

9

Sorry, it was my mistake, I have tried it now and it works :)))))). Before I have copied parameters to [pairs] section, even I intended to copy them into [bonds] section, hence I have got an error “no such pair type 10”.

Thank you very much for your help!

At the end is 1000, which is force constant or whatever is called that 4tht parameter for potential to restraint bonds of type 10.
Did you choose 1000 after testing what gives the best results? (Or maybe somewhere is recommended to use that value).

Greetings.

10

No worries!
1000 is a default value, I think, but I would rather pick it based on your specific molecules and how rigid you want it to be.
Good luck with your simulations!

Best,
Aleksei

11

Thanx :). For me even worked when I have put “force constant 10” although I was not checking thoroughly by RMSD how stable is conformation, but visually in VMD it looked pretty fine.