Calculating overall protein tumbling time

[Jessica M. Gonzalez-Delgado]

Hi everyone,

I was wondering if Carma can be used to calculate the overall protein
tumbling time from a long MD trajectory. If not, any suggestions on how
to do that? I'd like to compare with some experimental data. Thanks!

--
Jessica M. González-Delgado | she/her
Ph.D. Candidate
Department of Chemistry
Franzen Research Group
North Carolina State University

[Nicholas M. Glykos]

> I was wondering if Carma can be used to calculate the overall protein
> tumbling time from a long MD trajectory.

Probably not. Do you have a published example of exactly what you want to calculate ?




--


Nicholas M. Glykos, Department of Molecular Biology
and Genetics, Democritus University of Thrace, University Campus,
Dragana, 68100 Alexandroupolis, Greece, Tel/Fax (office) +302551030620,
Ext.77620, Tel (lab) +302551030615, https://utopia.duth.gr/glykos/

[Jessica M. Gonzalez-Delgado]

Hi Nicholas,

Thanks for your prompt response. I actually haven't found a paper on how
to calculate the overall tumbling time from MD simulations, through NMR
experiments it is done using the Model Free formalism by Lipari and
Szabo. I did NMR experiments to apply model free analysis and obtain the
S2 order parameter, and other parameter like the time of the motions of
the individual residues (tau_e) and the overall protein tumbling time
(tau_m). For the MD simulations you have to remove the overall protein
tumbling, which is what I use Carma primarily for, in order to get the
S2 and tau_e parameters from the MD simulations. However, I still
haven't seen how to get tau_m from MD. I guess it's just a curiosity at
this point, since I am not sure it is possible to calculate from MD in a
straightforward fashion.

Thanks again,
Jessica

On 2022-03-17 04:43, Nicholas M. Glykos wrote:
>
>> I was wondering if Carma can be used to calculate the overall protein
>> tumbling time from a long MD trajectory.
> Probably not. Do you have a published example of exactly what you want to calculate ?
>
>
>
>

--

[Nicholas M. Glykos]

Is (tau_m) related to the Rotational correlation time (τc)¹ ?
If yes, I could possibly imagine a way of calculating it from carma output.


¹ https://en.wikipedia.org/wiki/Rotational_correlation_time

[Jessica M. Gonzalez-Delgado]

Yes! It actually is tau_c, some people use tau_c and others use tau_m.

-Jessica

[Nicholas M. Glykos]

> Yes! It actually is tau_c, some people use tau_c and others use tau_m.

ok, for what follows assume that it is all nonsense, I don't have the time
right now to validate it :

When you run carma with "carma -v -fit ...", it produces a file named
'carma.fit-rms.dat' which in-between other things contains the rotation
matrix and the translation vector needed to superimpose the frame of
interest to the reference frame (by default the first). To calculate tau_c
you need to calculate the pure rotational offset between the (i) frame and the
frame of reference. This angle is the κ (kappa) angle in the polar angle
convention. So the problem is to convert the rotation matrix to the polar
angles and only keep (κ).

This is possibly easy : The rotation angle (κ) can be derived from the
trace of the rotation matrix :

Trace = R11 + R22 + R33 = 1 + 2 cos( κ )

κ = acos( (R11 + R22 + R33 - 1) / 2 )


To calculate this only once, is a one-liner :

awk 'function acos(x) { return atan2(sqrt(1-x*x), x)} {print acos(($3+$7+$11-1)/2)}' carma.fit-rms.dat

which returns (κ) in rad as a function of frame, using as reference the
first frame.


The thing is that what you want is not one run using as reference the first
frame, but a whole set of runs to calculate averages. This calls for a
small perl script to run carma several times using each time a different
frame of reference.

Verify that all these are correct and make sense, and I think I could find
the time to write the script sometime tomorrow.

[Jessica M. Gonzalez-Delgado]

Thanks, that makes sense to me. So, then the (average number of
frames)*(step_size) that takes kappa to have a value of 1.0 would be an
approximation of tau_c?

-Jessica

[Nicholas M. Glykos]

> Thanks, that makes sense to me. So, then the (average number of
> frames)*(step_size) that takes kappa to have a value of 1.0 would be an
> approximation of tau_c?

Correct. You can either average across the number of steps needed to
cross a value of 1.0, or you can calculate mean + standard deviation
of (κ) for each step (in a 'sliding' window across the trajectory).

Both methods should give identical results in the case of a sufficiently
long trajectory.

I have a lot of teaching these days, I'll try to cook-up the perl script
sometime tomorrow or the day after tomorrow.

[Nicholas M. Glykos]

You can try this perl script :

________________________________________________________________________________________________

use Math::Trig;

for ( $offset=0 ; $offset <= 100000 ; $offset += 1000 )
{
$first = $offset+1;
$reference = $offset+1;
$last = 30000 + $offset;
`carma -v -fit -segid A -segid B -first $first -ref $reference -last $last ../equi_out.dcd ../../notails.psf`;

open( RMS, "carma.fit-rms.dat");
while( $line = <RMS>)
{
@column = split( ' ', $line);
$kappa = acos( ($column[2] + $column[6] + $column[10] -1 ) / 2 );
if ( $kappa > 1.0 )
{
printf "%8d\n", $column[0] - $first;
last;
}
}
if ( $kappa <= 1.0 )
{
print "Increase maximal expected tau_c !\n";
}
close( RMS );
}

________________________________________________________________________________________________


You'll have to change :
- Paths for DCD & PSF (in the carma command)
- The SEGID selection in the carma command
- In the line $last = 30000 + $offset; the number 30000 should be a large enough
value so that the program does not produce the error "Increase maximal expected tau_c".
If it is too large, you'll waste time.
- How fine and how long will be the search along the trajectory are defined by the for()

The output is the list of frames at which the transition to κ>1 occurred.
You'll do the averaging yourself.

This is written late in the evening. I'd be surprised if it works for you.

[Jessica M. Gonzalez-Delgado]

Hi Nicholas,

Thanks for the script! I appreciate your help with this tau_c curiosity.
I tried the perl script but I was unsuccessful. However, I understand
what the script is doing, so I will make a bash script later and see how
it goes. Using the first frame as a reference I got tau_c equal to ~50
ns, which is about 5 times higher than the experimental value so I'll
take a look at rotation matrix and see if there is a factor of 5 somewhere.

-Jessica

[Nicholas M. Glykos]

> Thanks for the script! I appreciate your help with this tau_c curiosity. I
> tried the perl script but I was unsuccessful.

What was the error message ? I have tried it with a homodimeric protein of
about 120 residues and it was performing as expected.

[Nicholas M. Glykos]

If you get no output, then the carma command is wrong for your system.


Replace the line

`carma -v -fit -segid A -segid B -first $first -ref $reference -last $last ../equi_out.dcd ../../notails.psf`;

with

system("carma -v -fit -segid A -segid B -first $first -ref $reference -last $last ../equi_out.dcd ../../notails.psf");



so that you can see the carma error message.

[Jessica M. Gonzalez-Delgado]

Yup, that did it! Now carma is running. I'll try it on a longer
trajectory, but it should work as expected. Thanks again!

-Jessica

On 2022-03-18 15:02, Nicholas M. Glykos wrote:
>
> If you get no output, then the carma command is wrong for your system.
>
>
> Replace the line
>
> `carma -v -fit -segid A -segid B -first $first -ref $reference -last $last ../equi_out.dcd ../../notails.psf`;
>
> with
>
> system("carma -v -fit -segid A -segid B -first $first -ref $reference -last $last ../equi_out.dcd ../../notails.psf");
>
>
>
> so that you can see the carma error message.
>
>
>
>
>
>
> On Fri, Mar 18, 2022 at 08:47:20PM +0200, Nicholas M. Glykos wrote:
>>
>>> Thanks for the script! I appreciate your help with this tau_c curiosity. I
>>> tried the perl script but I was unsuccessful.
>> What was the error message ? I have tried it with a homodimeric protein of
>> about 120 residues and it was performing as expected.
>

--