Module ns_r1rho_2site

The numerical solution for the 2-site Bloch-McConnell equations for R1rho-type data, the NS R1rho 2-site model.

Description

This is the model of the numerical solution for the 2-site Bloch-McConnell equations. It originates from the funNumrho.m file from the Skrynikov & Tollinger code (the sim_all.tar file https://web.archive.org/web/https://gna.org/support/download.php?file_id=18404 attached to https://web.archive.org/web/https://gna.org/task/?7712#comment5). That code is:

   function residual = funNumrho(optpar)


   global nu_0 x y Rcalc rms nfields offset w1 Tc w1 R1_51 R1_81
   %keyboard
   Rcalc=zeros(nfields,size(w1,2));
   tau_ex=optpar(1);
   pb=optpar(2);
   pa=1-pb;
   kex=1/tau_ex; k_u=pb*kex; k_f=(1-pb)*kex;

   for k=1:nfields
       % we assume that A resonates at 0 [s^-1], without loss of generality
       dw=nu_0(k)*optpar(3)*2*pi;      % [s^-1]
       Wa=0*2*pi;                      % Larmor frequ. [s^-1]
       Wb=dw;                          % Larmor frequ. [s^-1]
       Wsl=offset*2*pi;                        % Larmor frequ. of spin lock [s^-1]
       da=Wa-Wsl;                              % offset of sl from A
       db=Wb-Wsl;                              % offset of sl from B


       Rinf=optpar(3+k);

       for kk=1:length(w1)
         w1x=w1(kk);
         if k==1; R1=R1_51; end; if k==2; R1=R1_81; end

         R=-[Rinf+k_u   -k_f       da       0       0       0     ;
              -k_u    Rinf+k_f      0      db       0       0     ;
              -da         0     Rinf+k_u  -k_f     w1x      0     ;
              0          -db      -k_u  Rinf+k_f    0      w1x    ;
              0           0       -w1x      0    R1+k_u   -k_f    ;
              0           0         0     -w1x    -k_u   R1+k_f  ];
       % keyboard
        MAx0= pa; MBx0= pb; MAy0= 0; MBy0= 0; MAz0= 0; MBz0= 0;
        Mof0=[MAx0 MBx0 MAy0 MBy0 MAz0 MBz0]';

   % the following lines: rotate the magnetization previous to spin lock into the weff frame
   % a new Mof0 is otained: Mof0=[sin(thetaa)*pa sin(thetab)*pb 0 0 cos(thetaa)*pa cos(thetab)*pb];
   thetaa=atan(w1x/da);thetaa_degrees=thetaa*360/(2*pi);
   thetab=atan(w1x/db);thetab_degrees=thetab*360/(2*pi);
   MAxnew=sin(thetaa)*pa;
   MBxnew=sin(thetab)*pb;
   MAynew=MAy0;
   MBynew=MBy0;
   MAznew=cos(thetaa)*pa;
   MBznew=cos(thetab)*pb;
   Mof0=[MAxnew MBxnew MAynew MBynew MAznew MBznew]';

        Moft(1:6)=(expm3(R*Tc)*Mof0)';
        MAx=real(Moft(1))/pa;
        MAy=real(Moft(3))/pa;
        MAz=real(Moft(5))/pa;

        MA=sqrt(MAx^2+MAy^2+MAz^2); % for spin A, is equal to 1 if nothing happens (no relaxation)
        intrat(k,kk)=MA;
        Rcalc(k,kk)=(-1.0/Tc)*log(intrat(k,kk));
       end

   end

   if (size(Rcalc)==size(y))
       residual=sum(sum((y-Rcalc).^2));
       rms=sqrt(residual/(size(y,1)*size(y,2)));
   end

References

The solution has been modified to use the from presented in:

• Korzhnev, D. M., Orekhov, V. Y., and Kay, L. E. (2005). Off-resonance R(1rho) NMR studies of exchange dynamics in proteins with low spin-lock fields: an application to a Fyn SH3 domain. J. Am. Chem. Soc., 127, 713-721. (DOI: 10.1021/ja0446855).

Links

More information on the NS R1rho 2-site model can be found in the:

• relax wiki,
• relax manual,
• relaxation dispersion page of the relax website.

Functions

numpy float array of rank [NE][NS][NM][NO][ND][6][6] rr1rho_3d_2site_rankN(R1=None, r1rho_prime=None, dw=None, omega=None, offset=None, w1=None, k_AB=None, k_BA=None, relax_time=None)
Definition of the multidimensional 3D exchange matrix, of rank [NE][NS][NM][NO][ND][6][6].

ns_r1rho_2site(M0=None, M0_T=None, r1rho_prime=None, omega=None, offset=None, r1=0.0, pA=None, dw=None, kex=None, spin_lock_fields=None, relax_time=None, inv_relax_time=None, back_calc=None)
The 2-site numerical solution to the Bloch-McConnell equation for R1rho data.

Variables

m_r1rho_prime = array([[-1., 0., 0., 0., 0., 0...

m_wA = array([[ 0., -1., 0., 0., 0., 0...

m_wB = array([[ 0., 0., 0., 0., 0., 0...

m_w1 = array([[ 0., 0., 0., 0., 0., 0...

m_k_AB = array([[-1., 0., 0., 0., 0., 0...

m_k_BA = array([[ 0., 0., 0., 1., 0., 0...

m_R1 = array([[ 0., 0., 0., 0., 0., 0...

__package__ = 'lib.dispersion'

Imports: array, einsum, float64, isfinite, log, min, multiply, sum, fix_invalid, masked_less, matrix_exponential

Function Details

rr1rho_3d_2site_rankN(R1=None, r1rho_prime=None, dw=None, omega=None, offset=None, w1=None, k_AB=None, k_BA=None, relax_time=None)

Definition of the multidimensional 3D exchange matrix, of rank [NE][NS][NM][NO][ND][6][6].

This code originates from the funNumrho.m file from the Skrynikov & Tollinger code (the sim_all.tar file https://web.archive.org/web/https://gna.org/support/download.php?file_id=18404 attached to https://web.archive.org/web/https://gna.org/task/?7712#comment5).

Parameters:

• R1 (numpy float array of rank [NE][NS][NM][NO][ND]) - The longitudinal, spin-lattice relaxation rate.
• r1rho_prime (numpy float array of rank [NE][NS][NM][NO][ND]) - The R1rho transverse, spin-spin relaxation rate in the absence of exchange.
• dw (numpy float array of rank [NE][NS][NM][NO][ND]) - The chemical exchange difference between states A and B in rad/s.
• omega (numpy float array of rank [NE][NS][NM][NO][ND]) - The chemical shift for the spin in rad/s.
• offset (numpy float array of rank [NE][NS][NM][NO][ND]) - The spin-lock offsets for the data.
• w1 (numpy float array of rank [NE][NS][NM][NO][ND]) - The spin-lock field strength in rad/s.
• k_AB (float) - The forward exchange rate from state A to state B.
• k_BA (float) - The reverse exchange rate from state B to state A.
• relax_time (numpy float array of rank [NE][NS][NM][NO][ND]) - The total relaxation time period for each spin-lock field strength (in seconds).

Returns: numpy float array of rank [NE][NS][NM][NO][ND][6][6]

The relaxation matrix.

ns_r1rho_2site(M0=None, M0_T=None, r1rho_prime=None, omega=None, offset=None, r1=0.0, pA=None, dw=None, kex=None, spin_lock_fields=None, relax_time=None, inv_relax_time=None, back_calc=None)

The 2-site numerical solution to the Bloch-McConnell equation for R1rho data.

This function calculates and stores the R1rho values.

Parameters:

• M0 (numpy float array of rank [NE][NS][NM][NO][ND][6][1]) - This is a vector that contains the initial magnetizations corresponding to the A and B state transverse magnetizations.
• M0_T (numpy float array of rank [NE][NS][NM][NO][ND][1][6]) - This is a vector that contains the initial magnetizations corresponding to the A and B state transverse magnetizations, where the outer two axis has been swapped for efficient dot operations.
• r1rho_prime (numpy float array of rank [NE][NS][NM][NO][ND]) - The R1rho_prime parameter value (R1rho with no exchange).
• omega (numpy float array of rank [NE][NS][NM][NO][ND]) - The chemical shift for the spin in rad/s.
• offset (numpy float array of rank [NE][NS][NM][NO][ND]) - The spin-lock offsets for the data.
• r1 (numpy float array of rank [NE][NS][NM][NO][ND]) - The R1 relaxation rate.
• pA (float) - The population of state A.
• dw (numpy float array of rank [NE][NS][NM][NO][ND]) - The chemical exchange difference between states A and B in rad/s.
• kex (float) - The kex parameter value (the exchange rate in rad/s).
• spin_lock_fields (numpy float array of rank [NE][NS][NM][NO][ND]) - The R1rho spin-lock field strengths (in rad.s^-1).
• relax_time (numpy float array of rank [NE][NS][NM][NO][ND]) - The total relaxation time period for each spin-lock field strength (in seconds).
• inv_relax_time (numpy float array of rank [NE][NS][NM][NO][ND]) - The inverse of the relaxation time period for each spin-lock field strength (in inverse seconds). This is used for faster calculations.
• back_calc (numpy float array of rank [NE][NS][NM][NO][ND]) - The array for holding the back calculated R2eff values. Each element corresponds to one of the CPMG nu1 frequencies.

Variables Details

m_r1rho_prime

Value:

array([[-1.,  0.,  0.,  0.,  0.,  0.],
       [ 0., -1.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0., -1.,  0.,  0.],
       [ 0.,  0.,  0.,  0., -1.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.]])

m_wA

Value:

array([[ 0., -1.,  0.,  0.,  0.,  0.],
       [ 1.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.]])

m_wB

Value:

array([[ 0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0., -1.,  0.],
       [ 0.,  0.,  0.,  1.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.]])

m_w1

Value:

array([[ 0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0., -1.,  0.,  0.,  0.],
       [ 0.,  1.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0., -1.],
       [ 0.,  0.,  0.,  0.,  1.,  0.]])

m_k_AB

Value:

array([[-1.,  0.,  0.,  0.,  0.,  0.],
       [ 0., -1.,  0.,  0.,  0.,  0.],
       [ 0.,  0., -1.,  0.,  0.,  0.],
       [ 1.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  1.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  1.,  0.,  0.,  0.]])

m_k_BA

Value:

array([[ 0.,  0.,  0.,  1.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  1.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  1.],
       [ 0.,  0.,  0., -1.,  0.,  0.],
       [ 0.,  0.,  0.,  0., -1.,  0.],
       [ 0.,  0.,  0.,  0.,  0., -1.]])

m_R1

Value:

array([[ 0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0., -1.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0., -1.]])