lib.dispersion.ns_r1rho

1 ############################################################################### 2 # # 3 # Copyright (C) 2000-2001 Nikolai Skrynnikov # 4 # Copyright (C) 2000-2001 Martin Tollinger # 5 # Copyright (C) 2013-2014 Edward d'Auvergne # 6 # Copyright (C) 2014 Troels E. Linnet # 7 # # 8 # This file is part of the program relax (http://www.nmr-relax.com). # 9 # # 10 # This program is free software: you can redistribute it and/or modify # 11 # it under the terms of the GNU General Public License as published by # 12 # the Free Software Foundation, either version 3 of the License, or # 13 # (at your option) any later version. # 14 # # 15 # This program is distributed in the hope that it will be useful, # 16 # but WITHOUT ANY WARRANTY; without even the implied warranty of # 17 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # 18 # GNU General Public License for more details. # 19 # # 20 # You should have received a copy of the GNU General Public License # 21 # along with this program. If not, see <http://www.gnu.org/licenses/>. # 22 # # 23 ############################################################################### 24 25 # Module docstring. 26 """The numerical solution for the 2-site Bloch-McConnell equations for R1rho-type data, the U{NS R1rho 2-site<http://wiki.nmr-relax.com/NS_R1rho_2-site>} model. 27 28 Description 29 =========== 30 31 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 U{https://web.archive.org/web/https://gna.org/support/download.php?file_id=18404} attached to U{https://web.archive.org/web/https://gna.org/task/?7712#comment5}). That code is:: 32 33 function residual = funNumrho(optpar) 34 35 36 global nu_0 x y Rcalc rms nfields offset w1 Tc w1 R1_51 R1_81 37 %keyboard 38 Rcalc=zeros(nfields,size(w1,2)); 39 tau_ex=optpar(1); 40 pb=optpar(2); 41 pa=1-pb; 42 kex=1/tau_ex; k_u=pb*kex; k_f=(1-pb)*kex; 43 44 for k=1:nfields 45 % we assume that A resonates at 0 [s^-1], without loss of generality 46 dw=nu_0(k)*optpar(3)*2*pi; % [s^-1] 47 Wa=0*2*pi; % Larmor frequ. [s^-1] 48 Wb=dw; % Larmor frequ. [s^-1] 49 Wsl=offset*2*pi; % Larmor frequ. of spin lock [s^-1] 50 da=Wa-Wsl; % offset of sl from A 51 db=Wb-Wsl; % offset of sl from B 52 53 54 Rinf=optpar(3+k); 55 56 for kk=1:length(w1) 57 w1x=w1(kk); 58 if k==1; R1=R1_51; end; if k==2; R1=R1_81; end 59 60 R=-[Rinf+k_u -k_f da 0 0 0 ; 61 -k_u Rinf+k_f 0 db 0 0 ; 62 -da 0 Rinf+k_u -k_f w1x 0 ; 63 0 -db -k_u Rinf+k_f 0 w1x ; 64 0 0 -w1x 0 R1+k_u -k_f ; 65 0 0 0 -w1x -k_u R1+k_f ]; 66 % keyboard 67 MAx0= pa; MBx0= pb; MAy0= 0; MBy0= 0; MAz0= 0; MBz0= 0; 68 Mof0=[MAx0 MBx0 MAy0 MBy0 MAz0 MBz0]'; 69 70 % the following lines: rotate the magnetization previous to spin lock into the weff frame 71 % a new Mof0 is otained: Mof0=[sin(thetaa)*pa sin(thetab)*pb 0 0 cos(thetaa)*pa cos(thetab)*pb]; 72 thetaa=atan(w1x/da);thetaa_degrees=thetaa*360/(2*pi); 73 thetab=atan(w1x/db);thetab_degrees=thetab*360/(2*pi); 74 MAxnew=sin(thetaa)*pa; 75 MBxnew=sin(thetab)*pb; 76 MAynew=MAy0; 77 MBynew=MBy0; 78 MAznew=cos(thetaa)*pa; 79 MBznew=cos(thetab)*pb; 80 Mof0=[MAxnew MBxnew MAynew MBynew MAznew MBznew]'; 81 82 Moft(1:6)=(expm3(R*Tc)*Mof0)'; 83 MAx=real(Moft(1))/pa; 84 MAy=real(Moft(3))/pa; 85 MAz=real(Moft(5))/pa; 86 87 MA=sqrt(MAx^2+MAy^2+MAz^2); % for spin A, is equal to 1 if nothing happens (no relaxation) 88 intrat(k,kk)=MA; 89 Rcalc(k,kk)=(-1.0/Tc)*log(intrat(k,kk)); 90 end 91 92 end 93 94 if (size(Rcalc)==size(y)) 95 residual=sum(sum((y-Rcalc).^2)); 96 rms=sqrt(residual/(size(y,1)*size(y,2))); 97 end 98 99 100 References 101 ========== 102 103 The solution has been modified to use the from presented in: 104 105 - 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. I{J. Am. Chem. Soc.}, B{127}, 713-721. (U{DOI: 10.1021/ja0446855<http://dx.doi.org/10.1021/ja0446855>}). 106 107 108 Links 109 ===== 110 111 More information on the NS R1rho 2-site model can be found in the: 112 113 - U{relax wiki<http://wiki.nmr-relax.com/NS_R1rho_2-site>}, 114 - U{relax manual<http://www.nmr-relax.com/manual/The_NS_2_site_R1_rho_model.html>}, 115 - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#NS_R1rho_2-site>}. 116 """ 117 118 # Python module imports. 119 from numpy import array, einsum, float64, isfinite, log, min, multiply, sum 120 from numpy.ma import fix_invalid, masked_less 121 122 # relax module imports. 123 from lib.dispersion.matrix_exponential import matrix_exponential 124 125 # Repetitive calculations (to speed up calculations). 126 m_r1rho_prime = array([ 127 [-1, 0, 0, 0, 0, 0], 128 [ 0, -1, 0, 0, 0, 0], 129 [ 0, 0, 0, 0, 0, 0], 130 [ 0, 0, 0, -1, 0, 0], 131 [ 0, 0, 0, 0, -1, 0], 132 [ 0, 0, 0, 0, 0, 0]], float64) 133 134 m_wA = array([ 135 [ 0, -1, 0, 0, 0, 0], 136 [ 1, 0, 0, 0, 0, 0], 137 [ 0, 0, 0, 0, 0, 0], 138 [ 0, 0, 0, 0, 0, 0], 139 [ 0, 0, 0, 0, 0, 0], 140 [ 0, 0, 0, 0, 0, 0]], float64) 141 142 m_wB = array([ 143 [ 0, 0, 0, 0, 0, 0], 144 [ 0, 0, 0, 0, 0, 0], 145 [ 0, 0, 0, 0, 0, 0], 146 [ 0, 0, 0, 0, -1, 0], 147 [ 0, 0, 0, 1, 0, 0], 148 [ 0, 0, 0, 0, 0, 0]], float64) 149 150 m_w1 = array([ 151 [ 0, 0, 0, 0, 0, 0], 152 [ 0, 0, -1, 0, 0, 0], 153 [ 0, 1, 0, 0, 0, 0], 154 [ 0, 0, 0, 0, 0, 0], 155 [ 0, 0, 0, 0, 0, -1], 156 [ 0, 0, 0, 0, 1, 0]], float64) 157 158 m_k_AB = array([ 159 [-1, 0, 0, 0, 0, 0], 160 [ 0, -1, 0, 0, 0, 0], 161 [ 0, 0, -1, 0, 0, 0], 162 [ 1, 0, 0, 0, 0, 0], 163 [ 0, 1, 0, 0, 0, 0], 164 [ 0, 0, 1, 0, 0, 0]], float64) 165 166 m_k_BA = array([ 167 [ 0, 0, 0, 1, 0, 0], 168 [ 0, 0, 0, 0, 1, 0], 169 [ 0, 0, 0, 0, 0, 1], 170 [ 0, 0, 0, -1, 0, 0], 171 [ 0, 0, 0, 0, -1, 0], 172 [ 0, 0, 0, 0, 0, -1]], float64) 173 174 m_R1 = array([ 175 [ 0, 0, 0, 0, 0, 0], 176 [ 0, 0, 0, 0, 0, 0], 177 [ 0, 0, -1, 0, 0, 0], 178 [ 0, 0, 0, 0, 0, 0], 179 [ 0, 0, 0, 0, 0, 0], 180 [ 0, 0, 0, 0, 0, -1]], float64) 181 182

183 -def 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):

184 """Definition of the multidimensional 3D exchange matrix, of rank [NE][NS][NM][NO][ND][6][6]. 185 186 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). 187 188 189 @keyword R1: The longitudinal, spin-lattice relaxation rate. 190 @type R1: numpy float array of rank [NE][NS][NM][NO][ND] 191 @keyword r1rho_prime: The R1rho transverse, spin-spin relaxation rate in the absence of exchange. 192 @type r1rho_prime: numpy float array of rank [NE][NS][NM][NO][ND] 193 @keyword dw: The chemical exchange difference between states A and B in rad/s. 194 @type dw: numpy float array of rank [NE][NS][NM][NO][ND] 195 @keyword omega: The chemical shift for the spin in rad/s. 196 @type omega: numpy float array of rank [NE][NS][NM][NO][ND] 197 @keyword offset: The spin-lock offsets for the data. 198 @type offset: numpy float array of rank [NE][NS][NM][NO][ND] 199 @keyword w1: The spin-lock field strength in rad/s. 200 @type w1: numpy float array of rank [NE][NS][NM][NO][ND] 201 @keyword k_AB: The forward exchange rate from state A to state B. 202 @type k_AB: float 203 @keyword k_BA: The reverse exchange rate from state B to state A. 204 @type k_BA: float 205 @keyword relax_time: The total relaxation time period for each spin-lock field strength (in seconds). 206 @type relax_time: numpy float array of rank [NE][NS][NM][NO][ND] 207 @return: The relaxation matrix. 208 @rtype: numpy float array of rank [NE][NS][NM][NO][ND][6][6] 209 """ 210 211 # Wa: The chemical shift offset of state A from the spin-lock. Larmor frequency [s^-1]. 212 Wa = omega 213 # Wb: The chemical shift offset of state A from the spin-lock. Larmor frequency [s^-1]. 214 Wb = omega + dw 215 216 # Population-averaged Larmor frequency [s^-1]. 217 #W = pA*Wa + pB*Wb 218 219 # Offset of spin-lock from A. 220 dA = Wa - offset 221 222 # Offset of spin-lock from B. 223 dB = Wb - offset 224 225 # Offset of spin-lock from population-average. 226 #d = W - offset 227 228 # Alias to original parameter name. 229 wA = dA 230 wB = dB 231 232 # Multiply and expand. 233 mat_r1rho_prime = multiply.outer( r1rho_prime * relax_time, m_r1rho_prime ) 234 235 mat_wA = multiply.outer( wA * relax_time, m_wA ) 236 mat_wB = multiply.outer( wB * relax_time, m_wB ) 237 238 mat_w1 = multiply.outer( w1 * relax_time, m_w1 ) 239 240 mat_k_AB = multiply.outer( k_AB * relax_time, m_k_AB ) 241 mat_k_BA = multiply.outer( k_BA * relax_time, m_k_BA ) 242 243 mat_R1 = multiply.outer( R1 * relax_time, m_R1 ) 244 245 # Collect matrix. 246 matrix = (mat_r1rho_prime + mat_wA + mat_wB 247 + mat_w1 + mat_k_AB + mat_k_BA 248 + mat_R1) 249 250 # Return the matrix. 251 return matrix

252 253

254 -def 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):

255 """The 2-site numerical solution to the Bloch-McConnell equation for R1rho data. 256 257 This function calculates and stores the R1rho values. 258 259 260 @keyword M0: This is a vector that contains the initial magnetizations corresponding to the A and B state transverse magnetizations. 261 @type M0: numpy float array of rank [NE][NS][NM][NO][ND][6][1] 262 @keyword M0_T: 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. 263 @type M0_T: numpy float array of rank [NE][NS][NM][NO][ND][1][6] 264 @keyword r1rho_prime: The R1rho_prime parameter value (R1rho with no exchange). 265 @type r1rho_prime: numpy float array of rank [NE][NS][NM][NO][ND] 266 @keyword omega: The chemical shift for the spin in rad/s. 267 @type omega: numpy float array of rank [NE][NS][NM][NO][ND] 268 @keyword offset: The spin-lock offsets for the data. 269 @type offset: numpy float array of rank [NE][NS][NM][NO][ND] 270 @keyword r1: The R1 relaxation rate. 271 @type r1: numpy float array of rank [NE][NS][NM][NO][ND] 272 @keyword pA: The population of state A. 273 @type pA: float 274 @keyword dw: The chemical exchange difference between states A and B in rad/s. 275 @type dw: numpy float array of rank [NE][NS][NM][NO][ND] 276 @keyword kex: The kex parameter value (the exchange rate in rad/s). 277 @type kex: float 278 @keyword spin_lock_fields: The R1rho spin-lock field strengths (in rad.s^-1). 279 @type spin_lock_fields: numpy float array of rank [NE][NS][NM][NO][ND] 280 @keyword relax_time: The total relaxation time period for each spin-lock field strength (in seconds). 281 @type relax_time: numpy float array of rank [NE][NS][NM][NO][ND] 282 @keyword inv_relax_time: The inverse of the relaxation time period for each spin-lock field strength (in inverse seconds). This is used for faster calculations. 283 @type inv_relax_time: numpy float array of rank [NE][NS][NM][NO][ND] 284 @keyword back_calc: The array for holding the back calculated R2eff values. Each element corresponds to one of the CPMG nu1 frequencies. 285 @type back_calc: numpy float array of rank [NE][NS][NM][NO][ND] 286 """ 287 288 # Once off parameter conversions. 289 pB = 1.0 - pA 290 k_BA = pA * kex 291 k_AB = pB * kex 292 293 # The matrix that contains all the contributions to the evolution, i.e. relaxation, exchange and chemical shift evolution. 294 R_mat = rr1rho_3d_2site_rankN(R1=r1, r1rho_prime=r1rho_prime, dw=dw, omega=omega, offset=offset, w1=spin_lock_fields, k_AB=k_AB, k_BA=k_BA, relax_time=relax_time) 295 296 # This matrix is a propagator that will evolve the magnetization with the matrix R. 297 Rexpo_mat = matrix_exponential(R_mat) 298 299 # Magnetization evolution. 300 Rexpo_M0_mat = einsum('...ij, ...jk', Rexpo_mat, M0) 301 302 # Magnetization evolution, which include all dimensions. 303 MA_mat = einsum('...ij, ...jk', M0_T, Rexpo_M0_mat)[:, :, :, :, :, 0, 0] 304 305 # Insert safe checks. 306 if min(MA_mat) < 0.0: 307 mask_min_MA_mat = masked_less(MA_mat, 0.0) 308 # Fill with high values. 309 MA_mat[mask_min_MA_mat.mask] = 1e100 310 311 # Do back calculation. 312 back_calc[:] = -inv_relax_time * log(MA_mat) 313 314 # Catch errors, taking a sum over array is the fastest way to check for 315 # +/- inf (infinity) and nan (not a number). 316 if not isfinite(sum(back_calc)): 317 # Replaces nan, inf, etc. with fill value. 318 fix_invalid(back_calc, copy=False, fill_value=1e100)

319

Source Code for Module lib.dispersion.ns_r1rho_2site