GROMACS version: 2024.1
GROMACS modification: Yes/No
Here post your question
I want to apply flat-bottomed position restraints. Unfortunately, I could not find any document explaining how to apply flat-bottomed position restraints in Gromacs. I’d be grateful if someone provides me a link to that or some detailed explanation.
I want to keep the protein channel limited to a cylindrical space in the lipid bilayer on the XY plane.
Dear @Dulaj
Flat-bottomed position restraints are defined in the topology file of what you want to restraint, and here you can find the details. The formalism is like the following
#ifdef flat_bottomed_restraint_name
[ position_restraints ]
; id funct g r k
atom_number 2 direction radius force_constant
atom_number 2 direction radius force_constant
atom_number 2 direction radius force_constant
...
...
#endif
where atom_number is the number of the atom to which the restraint is applied (the one in the topology), 2 is the flat bottom restraints, direction is the axis and the symmetry of the potential (the types are explained in the link I put before), radius is the length of the region of activation for the potential and force_constant is the force constant of the spring potential you are applying.
Your symmetry is having something that has an axis along z and that you want to constrain along x and y, so you should go for a cylindrical potential with the axis along z, that is, direction is 8. Then you want a flat bottomed potential that is inverse, as you want everything to stay WITHIN the cylinder, which means that you will need a negative radius. So, if you want everything to stay within a 2.5 radius cylinder then the radius will be -2.5. The force constant and the atoms to which you are applying the potential are system dependent and you have to chose them. You will also have to prepare a special restraint file (to be given to GROMACS with the -r flag during the pre-processing step with gmx grompp) where you prepare a configuration where all atoms of interest of your protein are projected on the axis of the cylinder, that is, the center of the pore in the bilayer where you want to keep your protein.
However, I would say that a protein embedded in a membrane shouldn’t unfold/get scattered in the bilayer. If this is what is happening, I guess you have to really take a look at the system because you need a good justification to force it in place!
Dear @obZehn ,
Thank you very much for the valuable information you provided!
Could you please tell me how I set the flat-bottomed position restrains using one of the gromacs commands such as gmx genrestr ? Or should I use a script (eg; Python script) ?
Also, as you have mentioned in this part, my justification is that I’m applying flat-bottomed position restraints as a centre of mass position restraints (I mean that I want to keep the protein channel without moving everywhere in the lipid membrane by applying a force to the centre of mass.). Even for a part of the equilibration to avoid moving the protein due to unequilibrated lipids.
I did some MD simulations with the protein embedded in a POPE bilayer. I observed sometimes protein had moved through the bilayer. Since I’m going to use smaller X and Y lengths of the lipid bilayer to decrease the computational cost, I’m afraid that the protein would be moved through the bilayer to reach the box boundaries. Even with PBC, it was harder for me to centre the protein in that bilayer system than in the process of centring the protein when it is simulated in a water box (for example when calculating RMSD.
So, a few points.
My understanding of how topology-set restraints work is the following. The reference is provided with the -r flag, which usually is a structure file (e.g. .gro) that is basically a copy of your system that it is used by GROMACS not as a starting point to update but as a reference. For example, in the gmx grompp processing you provide a structure with -c, which is the starting point, and you can also add a reference with -r. My understanding is that GROMACS will read the positions from the -cfile and, if restraints are set, will read the reference for the restraint from the structure file provided in the -r flag. As an example, if you are restraining the position of a certain atom with a typical harmonic restraint, then GROMACS will read the starting position from -c and will apply the bias, which is purely distance dependent (k parameter is fixed by you), by calculating how much the given restrained atom moves away from the reference, which is the position of itself in the -r file. That’s why the -rflag must contain, structurally speaking, a copy of your system. The positions can be different, but the number of atoms and their ordering must be the same as there will be a one-to-one correspondence between the structures in -c and -r to calculate the restraints.
Going back to the flat bottomed with a negative radius. In this case it is a bit tricky, as you are defining a region for your atoms, not a specific point for each one of them. So you want ALL the atoms of the protein to stay within an axis of a cylinder. If you just provide the same file in -c and -r, then every atom that you list in the #ifdef ... #endif section of the topology will stay within a certain radius with respect to its own starting position. This is nearly useless basically, it is basically equivalent to request to each atom to stay within a cylinder which z axis it’s going through the starting position of the atom. It is a bit hard maybe to visualize, but this is completely different from requiring that the WHOLE protein stays within a cylindrical volume. To achieve this, you have to modify the .gro file and project each atom of the protein that you want to restraint to the (x,y) value that defines the z axis along which you want to define the cylinder. As far as I know, there are no tools in GROMACS to do so, but you can do it easily with whatever scripting you like, python/bash/etc.
If you want just to restrain the center of mass, then a much easier way to go is to use the pull module of GROMACS. You can define the two groups in the .mdp file, e.g. protein and membrane, and put an harmonic restraint that keeps the center of mass of the two at a constant distance.
A general word of warning. Restraints are tricky things. You are applying forces to the system, so very roughly speaking you are injecting/extracting energy to modify the dynamics of the system to follow your will. Restraints can give rise to artifacts, like torque on the components, that may be nearly invisible but can do a lot of harm. You should always justify very well the WHY you are applying a restraint!
Again, I would really consider why you are applying the restraints. If the system pops out of the membrane, then the problem is probably the embedding/the forcefield/the equilibration/etc. Similarly if the protein is unfolding or more in general showing structural defects. If the system is diffusing laterally (which seems a bit weird as I would expect the time scale for something like that to be very large, but this is not really my expertise) and is breaking across boundaries, then you can always reconstruct it with pbc tools and recenter it, as PBC are only a representation problem. I would rather keep my system virgin from biases like a restraint as long as possible, especially if it is just for a recentering problem. For analysis, you have tools like gmx trjconv that give you the opportunity to fit and remove the translation along the xy plane, for example.
Similarly, saving computational time is important, but be careful to not make the box too small, otherwise you will end up with fast systems that have systematic artifacts. As a rule of thumb, I think you should have at least a distance that is at least your cut-off for electrostatics from the sides of the box, but taking a look at previous literature from big groups that study things similar to your will for sure guide you better than me. My point here is that I do not think it may be safe for you to justify a restraint just by saying ‘I do not want my protein to touch the borders of the box’, since as you said PBC are there and virtually there is no boundary. Mine are just words of caution though, as you are the expert of your system and you should decide how to proceed.
Hope this helps!
Dear @obZehn ,
Thank you so much for your time spent on providing these valuable things.
Could you please explain how to apply the pull module? Even for future needs, I think it will be important for me to know about it. :)
I can understand the points you provided on how these restraints will be suitable for my system. I am also thinking about how to avoid these types of things that may cause artifacts.
The problem is that once I tried to recenter the system I couldn’t do it easily as in water box simulations. I remember I could do it for one frame using gmx editconf or gmx trajconv. However, I couldn’t center the trajectory in that way using gmx trajconv tools like -center / -pbc mol … etc.
Thank you!
Dear @Dulaj,
the options for the pull code in the manual. It is well documented. Ideally in your case you can define as the two pull groups the group of the protein and that of the membrane and apply an harmonic restraint to the distance of the two centres of mass. If you set a null rate then the distance between the COMs will be kept restrained by an umbrella potential. This is different though from a cylindrical restraint, though.
Yes, recentering can notoriously be a bit of a pain, but I guess the gmx trjconv -fit should really come in handy if you want to remove the xy diffusion!