[gmx-users] Obtaining PMF for change in domain position

[Abhi Acharya]

Hello Gmx Users,

As suggested previously in this thread, I tried to use enforced rotations
on my system. First I calculated the amount of domain rotation and
translation required to obtain conformation 2 from conformation 1 using
'measure rotation' tool in chimera. I received an output as below:

  Matrix rotation and translation
    -0.96887254  -0.01658882  -0.24700365  59.55727674
     0.17648032  -0.74599360  -0.64214347  85.34692780
    -0.17361074  -0.66574646   0.72570032  11.00502990
  Axis  -0.11353457  -0.35303290   0.92869676     (Axis vector for rotation
)
  Axis point  28.32087779  46.35125870   0.00000000     (point on the axis)
  Rotation angle (degrees) 174.03353949
  Shift along axis -26.67174721  (in Angstroms)

For generating the windows for umbrella sampling, I used the pull code for
generating the translation and enforced rotations for the required rotation
simultaneously. For pull code I used the axis vector as the pull vector and
a pulling rate of -0.0002667 nm/ps for a 10 ns simulation. For enforced
rotations, I used the same Axis vector as the rotation vector and a
rotation rate of 0.0174 deg/ps.  I have provided the relevant parameter
sections from the .mdp file below.

Although I can see the expected movement of the domain of interest, but in
the final conformation it fails to reach the intended final position; falls
short by ~21 Angstroms. Visual analysis revealed that although the domain
rotated correctly, the translation of the domain did not work effectively,
so I suspect the pull code didn't work as I expected. What could be the
possible reasons for such a behavior ? Maybe I am doing some thing wrong.

The link to image containing the superposition of the initial, final
obtained and final expected conformations is provided in the link for
reference.
https://www.dropbox.com/s/7fnhpfq7zmzw1d5/image.jpg?dl=0


The .mdp portion is as follows:

; COM PULLING
; Pull type: no, umbrella, constraint or constant-force
pull                     = umbrella
; Pull geometry: distance, direction, cylinder or position
pull_geometry            = direction
; Select components for the pull vector. default: Y Y Y
pull_dim                 = Y Y Y
; Cylinder radius for dynamic reaction force groups (nm)
pull_r1                  = 1
; Switch from r1 to r0 in case of dynamic reaction force
pull_r0                  = 1.5
pull_constr_tol          = 1e-06
pull_start               = no
pull_nstxout             = 10
pull_nstfout             = 1
; Number of pull groups
pull_ngroups             = 1
; Group name, weight (default all 1), vector, init, rate (nm/ps),
kJ/(mol*nm^2)
pull_group0              = KH-domain
pull_weights0            =
pull_pbcatom0            = 0
pull_group1              = GD-domain
pull_weights1            =
pull_pbcatom1            = 0
pull_vec1                = -0.11353457 -0.35303290 0.92869676
pull_init1               = 0.0
pull_rate1               = -0.0002667
pull_k1                  = 1000
pull_kB1                 = 1000

; ENFORCED ROTATION
; Enforced rotation: No or Yes
rotation                 = Yes
; Output frequency for angle, torque and rotation potential energy for the
whole group
rot_nstrout              = 100
; Output frequency for per-slab data (angles, torques and slab centers)
rot_nstsout              = 1000
; Number of rotation groups
rot_ngroups              = 1
; Rotation group name
rot_group0               = GD-C-alpha
; Rotation potential. Can be iso, iso-pf, pm, pm-pf, rm, rm-pf, rm2,
rm2-pf, flex, flex-t, flex2, flex2-t
rot_type0                = flex
; Use mass-weighting of the rotation group positions
rot_massw0               = no
; Rotation vector, will get normalized
rot_vec0                 = -0.11353457 -0.35303290 0.92869676
; Pivot point for the potentials iso, pm, rm, and rm2 (nm)
rot_pivot0               = 2.832087779  4.635125870   0.00000000
; Rotation rate (degree/ps) and force constant (kJ/(mol*nm^2))
rot_rate0                = -0.0174
rot_k0                   = 400
; Slab distance for flexible axis rotation (nm)
rot_slab_dist0           = 1.5
; Minimum value of Gaussian function for the force to be evaluated (for
flex* potentials)
rot_min_gauss0           = 0.001
; Value of additive constant epsilon' (nm^2) for rm2* and flex2* potentials
rot_eps0                 = 0
; Fitting method to determine angle of rotation group (rmsd, norm, or
potential)
rot_fit_method0          = norm
; For fit type 'potential', nr. of angles around the reference for which
the pot. is evaluated
rot_potfit_nsteps0       = 21
; For fit type 'potential', distance in degrees between two consecutive
angles
rot_potfit_step0         = 0.25


Thanks in advance.


Abhishek Acharya
Shasara Research Foundation

On Tue, Dec 30, 2014 at 5:15 PM, Abhi Acharya <abhi117acharya at gmail.com>
wrote:

>
>
> Abhishek Acharya
> Shasara Research Foundation
>
> On Tue, Dec 30, 2014 at 10:42 AM, Justin Lemkul <jalemkul at vt.edu> wrote:
>
>>
>>
>> On 12/29/14 11:34 AM, Abhi Acharya wrote:
>>
>>> On Mon, Dec 29, 2014 at 8:59 PM, Justin Lemkul <jalemkul at vt.edu> wrote:
>>>
>>>
>>>>
>>>> On 12/29/14 6:57 AM, Abhi Acharya wrote:
>>>>
>>>>  Hello GROMACS Users,
>>>>>
>>>>> This is a problem I am facing for the first time. Kindly guide be to
>>>>> the
>>>>> best options.
>>>>>
>>>>> I have a protein which has two large domains connected by a flexible
>>>>> linker
>>>>> peptide (~10 aa). The two domains seem to interact with each other and
>>>>> have
>>>>> been crystallized in three different conformations. I want to calculate
>>>>> the
>>>>> change in binding energy of the two domains wrt change in their
>>>>> relative
>>>>> position, i.e. keeping position of one domain constant what is change
>>>>> in
>>>>> binding energy as the other domain moves from conformation 1, through
>>>>> conformation 2 to finally, conformation 3. What is the best way to do
>>>>> so?
>>>>>
>>>>>
>>>>>  You need to describe what these three conformations are.  Do they
>>>> involve
>>>> rotations of the domains with respect to one another, or is it a simple
>>>> linear distance that varies?
>>>>
>>>
>>>
>>> Yes, there are significant rotations+translations along different axes
>>> involved which makes it challenging. Is there a way to model such
>>> movements
>>> using the pull code?
>>>
>>>
>> Perhaps the enforced rotation options can be useful here.
>
>
> I had a brief look at these options; seem to be appropriate for my
> problem. I will try them first
>
>>
>>
>>
>>>> At first, I considered using umbrella sampling to address the problem.
>>>> Here
>>>>
>>>>> too, I was not sure how to write a pull code for such a complex
>>>>> reaction
>>>>> coordinate. Also, even if I can generate the conformations, what would
>>>>> be
>>>>> the ideal number of windows that would be needed to correctly generate
>>>>> a
>>>>> PMF ? My feeling is it would be a large number, but again it just a
>>>>> guess.
>>>>> Is there is a better way than this?
>>>>>
>>>>>
>>>>>  Hard to know a priori.  You may need some trial and error here, but
>>>> with
>>>> umbrella sampling it's trivial to just go back and add windows, since
>>>> the
>>>> simulations are still all independent of one another.
>>>>
>>>>
>>> My concern is that since the domain movement from initial to final
>>> conformation is of about 10 nm, the number of windows required maybe too
>>> large. I read that the windows need to be spaced close enough so that
>>> each
>>> one samples some part of the next window. My system is also very large
>>> (200K atoms), so it seems to be computationally very expensive.
>>>
>>>
>> The number of windows needed is largely a function of the force constant,
>> e.g. how strong the biasing potential is.  But that's also a function of
>> how strong the interactions are that are intrinsic to the system.
>>
>>  I was just now thinking to reduce the problem to  conducting umbrella
>>> sampling for each conformation, simply using COM distance of the two
>>> domains as the reaction coordinate to obtain the PMF in each case. From
>>> the
>>> PMF, the binding energy of the domains in each conformation can be
>>> obtain.
>>> This will allow me to circumvent the complications introduced by domain
>>> rotations and maybe reduce the number of simulations required a bit .
>>> Will
>>> this be a correct approach?
>>>
>>>
>> Possibly, but COM distance can be a degenerate measure in the case of
>> rotation.  The PLUMED plugin may be useful here rather than standard
>> Gromacs options, otherwise try to find an example in the literature and
>> follow a known procedure.
>
>
> I was suggested the same by someone else too, just thought to try options
> in gromacs before I try something unknown to me. I would try PLUMED in case
> the enforced rotations don't work.
>
> Thank you for the suggestions, always appreciated.
>
> Regards,
> Abhishek
>