Phi-angle problem?

[spreus]

Dear Xiang-Jun,

First of all, your work with 3DNA has been very inspiring and useful for me.

The problem:
I am currently trying to implement the CEHS scheme in a (Matlab) program of my own which I would like to be compatible with the base-pair step parameters output of a 3DNA-analysis. Following almost the exact matrix equations provided in the SCHNArP paper (J Mol Biol, 1997) and the original CEHS paper by El Hassan & Calladine (J Mol Biol, 1995) I've made a script that rebuilds a structure based on bp step parameters as the input. The program produces a .mol file so I can check the output structure with the output structure from the 3DNA-rebuilding option using the same input parameters.
However, it now appears that the produced structure is identical to the corresponding 3DNA output only for small roll/tilt values, whereas for large roll/tilt it is not. The deviation seems to be depending on the position in the sequence: E.g. using the input parameters provided in the attached file the 3DNA/my_script structures virtually deviates only at the last three base-pairs. Then, making only the last three base-pairs (bottom in attached file) the relative orientation and position of these three base pairs is not the same as in the 10 bp structure, despite their local bp step parameters are set to be the same.

The question:
Since the problem seems to be related to roll/tilt I'm suspecting it lies in my calculation of the angle 'Phi'. The value I'm refering to is the 'Phi' in equation (18-20) in the SCHNArP paper and also in equation (1) in the 3DNA paper (NAR, 2003). Therefore, I hope the answer to the following question might solve my problem:
- How exactly does SCHNArP and 3DNA calculate the 'Phi'-angle of each bp-step when rebuilding a structure based on shift/slide/rise/tilt/roll/twist inputs? Does it use an identity matrix to represent the MST at every base pair step or just the first one?

I hope you have time to help with this despite it's not directly related to the use of 3DNA.

All the best,
Søren

PS. If it is of any help, this is the relevant piece of my Matlab code:

Code: [Select]

...

 % Make first bp:
    r_i = [0 0 0];
    T_i = eye(3,3);
    MST = eye(3,3);
    XYZ(:,1-3) = xyz(:,1-3); % Here XYZ is coordinates in global ref. frame, xyz is in the local bp ref. frame.

 % Make consecutive bp's:
for j = 2:NoOfBasePairs
  ***(it first calls input parameters and coordinates in the bp ref. frames). Then:***

    RollTilt = sqrt(roll^2+tilt^2);
    rolltiltaxis = normr((tilt/RollTilt)*MST(:,1)'+(roll/RollTilt)*MST(:,2)');

    Phi = acos(dot(rolltiltaxis,MST(:,2)));
     if dot(cross(rolltiltaxis,MST(:,2)),MST(:,3)) > 0
         Phi = sqrt(Phi^2);
     elseif dot(cross(rolltiltaxis,MST(:,2)),MST(:,3)) < 0
         Phi = -sqrt(Phi^2);
     end

     % Make T matrix and translation vector:
     T_i_mst = Rz(twist/2-Phi)*Ry(RollTilt/2)*Rz(Phi);
     T_i_ii = Rz(twist/2-Phi)*Ry(RollTilt)*Rz(twist/2+Phi);
     r_i_ii = [shift slide rise]*T_i_mst';
     
     T_g_ii = T_i_ii*T_g_i;
     r_g_ii = r_g_i+r_i_ii*T_g_ii';
     
    % Rotation and translation of base pair ii:
     run = xyz(:,(j*2+j-2):(j*2+j))*T_g_ii';
     XYZ(:,(j*2+j-2):(j*2+j)) = [run(:,1)+r_g_ii(1) run(:,2)+r_g_ii(2) run(:,3)+r_g_ii(3)];
   
     T_g_i = T_g_ii;
     r_g_i = r_g_ii;
   
end

% Then, visualization of coordinates in XYZ
...


[spreus]

Hi again,

I found the problem. It was not in my calculation of Phi but in the calculation of T_g_ii. I see now it is supposed to be T_g_ii = T_g_i*T_i_ii (and not T_i_ii*T_g_i). Sorry if you have spent time looking into it.

I guess there is a typo in the SCHNArP paper equation 21 that confused me.

Best,
Søren

[xiangjun]

Hi Søren,

Thanks for your interest in SCHNAaP/SCHNArP and 3DNA. I am impressed by your diligence in digging into the bottom of seemingly trivial details; in my experience, they are what count the most.

Regarding the definition of the phi-angle, reading El Hassan and Calladine's original JMB paper, and my SCHNAaP/SCHNArP JMB publications would certainly help. However, they are no substitutes of going through a worked example -- that's exactly what I put into Section #5 of the 3DNA user's manual (http://3dna.rutgers.edu:8080/3DNA_v1.5/x3dna.pdf). In 3DNA v2.0 distribution, I extracted that section into a separate file named "tech_details.pdf". For any serious 3DNA users, I strongly suggest them to work through it carefully in order to understand each step thoroughly. If you want to see my SCHNAaP/SCHNArP implementation in ANSI C, please let me know.

I am glad that you "found the problem". Regarding "a typo in the SCHNArP paper equation 21" that confused you, I agree that there is indeed a typo -- thank for bring this issue to me attention:
[pre:2wbo9dyp]T^g_{i+1} = [red:2wbo9dyp]T^T_i[/red:2wbo9dyp] * T^i_{i+1}   -- (21) from the SCHNArP paper
should be
T^g_{i+1} = [red:2wbo9dyp]T^g_i[/red:2wbo9dyp] * T^i_{i+1}[/pre:2wbo9dyp]
 I.e., the [red:2wbo9dyp]red[/red:2wbo9dyp] term should have a superscript [red:2wbo9dyp]g[/red:2wbo9dyp] instead of [red:2wbo9dyp]T[/red:2wbo9dyp]. Note that two terms following the equal sign is the correct order.

There are many subtle details that must be properly accounted for in order to get the correct result. Note that the equations may be derived alternatively in a mathematically equivalent form; this may explain why you have to swap the two right hand terms in equation 21 to solve your puzzle.

HTH,

Xiang-Jun