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,2019 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 3-site Bloch-McConnell equations for R1rho-type data, the U{NS R1rho 3-site linear<http://wiki.nmr-relax.com/NS_R1rho_3-site_linear>} and U{NS R1rho 3-site<http://wiki.nmr-relax.com/NS_R1rho_3-site>} model. 27 28 Description 29 =========== 30 31 This is the model of the numerical solution for the 3-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}). 32 33 34 References 35 ========== 36 37 The solution has been modified to use the from presented in: 38 39 - 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>}). 40 41 42 Links 43 ===== 44 45 More information on the NS R1rho 3-site linear model can be found in the: 46 47 - U{relax wiki<http://wiki.nmr-relax.com/NS_R1rho_3-site_linear>}, 48 - U{relax manual<http://www.nmr-relax.com/manual/The_NS_3_site_linear_R1_rho_model.html>}, 49 - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#NS_R1rho_3-site_linear>}. 50 51 More information on the NS R1rho 3-site model can be found in the: 52 53 - U{relax wiki<http://wiki.nmr-relax.com/NS_R1rho_3-site>}, 54 - U{relax manual<http://www.nmr-relax.com/manual/The_NS_3_site_R1_rho_model.html>}, 55 - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#NS_R1rho_3-site>}. 56 """ 57 58 # Python module imports. 59 from numpy import array, einsum, float64, isfinite, log, min, multiply, sum 60 from numpy.ma import fix_invalid, masked_less 61 62 # relax module imports. 63 from lib.dispersion.matrix_exponential import matrix_exponential 64 65 # Repetitive calculations (to speed up calculations). 66 m_R1 = array([ 67 [ 0, 0, 0, 0, 0, 0, 0, 0, 0], 68 [ 0, 0, 0, 0, 0, 0, 0, 0, 0], 69 [ 0, 0, -1, 0, 0, 0, 0, 0, 0], 70 [ 0, 0, 0, 0, 0, 0, 0, 0, 0], 71 [ 0, 0, 0, 0, 0, 0, 0, 0, 0], 72 [ 0, 0, 0, 0, 0, -1, 0, 0, 0], 73 [ 0, 0, 0, 0, 0, 0, 0, 0, 0], 74 [ 0, 0, 0, 0, 0, 0, 0, 0, 0], 75 [ 0, 0, 0, 0, 0, 0, 0, 0, -1]], float64) 76 77 m_r1rho_prime = array([ 78 [-1, 0, 0, 0, 0, 0, 0, 0, 0], 79 [ 0, -1, 0, 0, 0, 0, 0, 0, 0], 80 [ 0, 0, 0, 0, 0, 0, 0, 0, 0], 81 [ 0, 0, 0, -1, 0, 0, 0, 0, 0], 82 [ 0, 0, 0, 0, -1, 0, 0, 0, 0], 83 [ 0, 0, 0, 0, 0, 0, 0, 0, 0], 84 [ 0, 0, 0, 0, 0, 0, -1, 0, 0], 85 [ 0, 0, 0, 0, 0, 0, 0, -1, 0], 86 [ 0, 0, 0, 0, 0, 0, 0, 0, 0]], float64) 87 88 m_wA = array([ 89 [ 0, -1, 0, 0, 0, 0, 0, 0, 0], 90 [ 1, 0, 0, 0, 0, 0, 0, 0, 0], 91 [ 0, 0, 0, 0, 0, 0, 0, 0, 0], 92 [ 0, 0, 0, 0, 0, 0, 0, 0, 0], 93 [ 0, 0, 0, 0, 0, 0, 0, 0, 0], 94 [ 0, 0, 0, 0, 0, 0, 0, 0, 0], 95 [ 0, 0, 0, 0, 0, 0, 0, 0, 0], 96 [ 0, 0, 0, 0, 0, 0, 0, 0, 0], 97 [ 0, 0, 0, 0, 0, 0, 0, 0, 0]], float64) 98 99 m_wB = array([ 100 [ 0, 0, 0, 0, 0, 0, 0, 0, 0], 101 [ 0, 0, 0, 0, 0, 0, 0, 0, 0], 102 [ 0, 0, 0, 0, 0, 0, 0, 0, 0], 103 [ 0, 0, 0, 0, -1, 0, 0, 0, 0], 104 [ 0, 0, 0, 1, 0, 0, 0, 0, 0], 105 [ 0, 0, 0, 0, 0, 0, 0, 0, 0], 106 [ 0, 0, 0, 0, 0, 0, 0, 0, 0], 107 [ 0, 0, 0, 0, 0, 0, 0, 0, 0], 108 [ 0, 0, 0, 0, 0, 0, 0, 0, 0]], float64) 109 110 m_wC = array([ 111 [ 0, 0, 0, 0, 0, 0, 0, 0, 0], 112 [ 0, 0, 0, 0, 0, 0, 0, 0, 0], 113 [ 0, 0, 0, 0, 0, 0, 0, 0, 0], 114 [ 0, 0, 0, 0, 0, 0, 0, 0, 0], 115 [ 0, 0, 0, 0, 0, 0, 0, 0, 0], 116 [ 0, 0, 0, 0, 0, 0, 0, 0, 0], 117 [ 0, 0, 0, 0, 0, 0, 0, -1, 0], 118 [ 0, 0, 0, 0, 0, 0, 1, 0, 0], 119 [ 0, 0, 0, 0, 0, 0, 0, 0, 0]], float64) 120 121 m_w1 = array([ 122 [ 0, 0, 0, 0, 0, 0, 0, 0, 0], 123 [ 0, 0, -1, 0, 0, 0, 0, 0, 0], 124 [ 0, 1, 0, 0, 0, 0, 0, 0, 0], 125 [ 0, 0, 0, 0, 0, 0, 0, 0, 0], 126 [ 0, 0, 0, 0, 0, -1, 0, 0, 0], 127 [ 0, 0, 0, 0, 1, 0, 0, 0, 0], 128 [ 0, 0, 0, 0, 0, 0, 0, 0, 0], 129 [ 0, 0, 0, 0, 0, 0, 0, 0, -1], 130 [ 0, 0, 0, 0, 0, 0, 0, 1, 0]], float64) 131 132 m_k_AB = array([ 133 [-1, 0, 0, 0, 0, 0, 0, 0, 0], 134 [ 0, -1, 0, 0, 0, 0, 0, 0, 0], 135 [ 0, 0, -1, 0, 0, 0, 0, 0, 0], 136 [ 1, 0, 0, 0, 0, 0, 0, 0, 0], 137 [ 0, 1, 0, 0, 0, 0, 0, 0, 0], 138 [ 0, 0, 1, 0, 0, 0, 0, 0, 0], 139 [ 0, 0, 0, 0, 0, 0, 0, 0, 0], 140 [ 0, 0, 0, 0, 0, 0, 0, 0, 0], 141 [ 0, 0, 0, 0, 0, 0, 0, 0, 0]], float64) 142 143 m_k_BA = array([ 144 [ 0, 0, 0, 1, 0, 0, 0, 0, 0], 145 [ 0, 0, 0, 0, 1, 0, 0, 0, 0], 146 [ 0, 0, 0, 0, 0, 1, 0, 0, 0], 147 [ 0, 0, 0, -1, 0, 0, 0, 0, 0], 148 [ 0, 0, 0, 0, -1, 0, 0, 0, 0], 149 [ 0, 0, 0, 0, 0, -1, 0, 0, 0], 150 [ 0, 0, 0, 0, 0, 0, 0, 0, 0], 151 [ 0, 0, 0, 0, 0, 0, 0, 0, 0], 152 [ 0, 0, 0, 0, 0, 0, 0, 0, 0]], float64) 153 154 m_k_BC = array([ 155 [ 0, 0, 0, 0, 0, 0, 0, 0, 0], 156 [ 0, 0, 0, 0, 0, 0, 0, 0, 0], 157 [ 0, 0, 0, 0, 0, 0, 0, 0, 0], 158 [ 0, 0, 0, -1, 0, 0, 0, 0, 0], 159 [ 0, 0, 0, 0, -1, 0, 0, 0, 0], 160 [ 0, 0, 0, 0, 0, -1, 0, 0, 0], 161 [ 0, 0, 0, 1, 0, 0, 0, 0, 0], 162 [ 0, 0, 0, 0, 1, 0, 0, 0, 0], 163 [ 0, 0, 0, 0, 0, 1, 0, 0, 0]], float64) 164 165 m_k_CB = array([ 166 [ 0, 0, 0, 0, 0, 0, 0, 0, 0], 167 [ 0, 0, 0, 0, 0, 0, 0, 0, 0], 168 [ 0, 0, 0, 0, 0, 0, 0, 0, 0], 169 [ 0, 0, 0, 0, 0, 0, 1, 0, 0], 170 [ 0, 0, 0, 0, 0, 0, 0, 1, 0], 171 [ 0, 0, 0, 0, 0, 0, 0, 0, 1], 172 [ 0, 0, 0, 0, 0, 0, -1, 0, 0], 173 [ 0, 0, 0, 0, 0, 0, 0, -1, 0], 174 [ 0, 0, 0, 0, 0, 0, 0, 0, -1]], float64) 175 176 m_k_AC = array([ 177 [-1, 0, 0, 0, 0, 0, 0, 0, 0], 178 [ 0, -1, 0, 0, 0, 0, 0, 0, 0], 179 [ 0, 0, -1, 0, 0, 0, 0, 0, 0], 180 [ 0, 0, 0, 0, 0, 0, 0, 0, 0], 181 [ 0, 0, 0, 0, 0, 0, 0, 0, 0], 182 [ 0, 0, 0, 0, 0, 0, 0, 0, 0], 183 [ 1, 0, 0, 0, 0, 0, 0, 0, 0], 184 [ 0, 1, 0, 0, 0, 0, 0, 0, 0], 185 [ 0, 0, 1, 0, 0, 0, 0, 0, 0]], float64) 186 187 m_k_CA = array([ 188 [ 0, 0, 0, 0, 0, 0, 1, 0, 0], 189 [ 0, 0, 0, 0, 0, 0, 0, 1, 0], 190 [ 0, 0, 0, 0, 0, 0, 0, 0, 1], 191 [ 0, 0, 0, 0, 0, 0, 0, 0, 0], 192 [ 0, 0, 0, 0, 0, 0, 0, 0, 0], 193 [ 0, 0, 0, 0, 0, 0, 0, 0, 0], 194 [ 0, 0, 0, 0, 0, 0, -1, 0, 0], 195 [ 0, 0, 0, 0, 0, 0, 0, -1, 0], 196 [ 0, 0, 0, 0, 0, 0, 0, 0, -1]], float64) 197 198

199 -def rr1rho_3d_3site_rankN(R1=None, r1rho_prime=None, dw_AB=None, dw_AC=None, omega=None, offset=None, w1=None, k_AB=None, k_BA=None, k_BC=None, k_CB=None, k_AC=None, k_CA=None, relax_time=None):

200 """Definition of the 3D exchange matrix. 201 202 @keyword R1: The longitudinal, spin-lattice relaxation rate. 203 @type R1: numpy float array of rank [NE][NS][NM][NO][ND] 204 @keyword r1rho_prime: The R1rho transverse, spin-spin relaxation rate in the absence of exchange. 205 @type r1rho_prime: numpy float array of rank [NE][NS][NM][NO][ND] 206 @keyword omega: The chemical shift for the spin in rad/s. 207 @type omega: numpy float array of rank [NS][NM][NO][ND] 208 @keyword offset: The spin-lock offsets for the data. 209 @type offset: numpy float array of rank [NE][NS][NM][NO][ND] 210 @keyword dw_AB: The chemical exchange difference between states A and B in rad/s. 211 @type dw_AB: numpy float array of rank [NE][NS][NM][NO][ND] 212 @keyword dw_AC: The chemical exchange difference between states A and C in rad/s. 213 @type dw_AC: numpy float array of rank [NE][NS][NM][NO][ND] 214 @keyword w1: The spin-lock field strength in rad/s. 215 @type w1: numpy float array of rank [NE][NS][NM][NO][ND] 216 @keyword k_AB: The forward exchange rate from state A to state B. 217 @type k_AB: float 218 @keyword k_BA: The reverse exchange rate from state B to state A. 219 @type k_BA: float 220 @keyword k_BC: The forward exchange rate from state B to state C. 221 @type k_BC: float 222 @keyword k_CB: The reverse exchange rate from state C to state B. 223 @type k_CB: float 224 @keyword k_AC: The forward exchange rate from state A to state C. 225 @type k_AC: float 226 @keyword k_CA: The reverse exchange rate from state C to state A. 227 @type k_CA: float 228 @keyword relax_time: The total relaxation time period for each spin-lock field strength (in seconds). 229 @type relax_time: numpy float array of rank [NE][NS][NM][NO][ND] 230 """ 231 232 # Repetitive calculations (to speed up calculations). 233 # The chemical shift offset of state A from the spin-lock. Larmor frequency for state A [s^-1]. 234 Wa = omega 235 236 # The chemical shift offset of state B from the spin-lock. Larmor frequency for state B [s^-1]. 237 Wb = omega + dw_AB 238 239 # The chemical shift offset of state C from the spin-lock. Larmor frequency for state C [s^-1]. 240 Wc = omega + dw_AC 241 242 # Population-averaged Larmor frequency [s^-1]. 243 #W = pA*Wa + pB*Wb + pC*Wc 244 245 # Offset of spin-lock from A. 246 dA = Wa - offset 247 248 # Offset of spin-lock from B. 249 dB = Wb - offset 250 251 # Offset of spin-lock from C. 252 dC = Wc - offset 253 254 # Offset of spin-lock from population-average. 255 #d = W - offset_i 256 257 # Parameter alias. 258 wA = dA 259 wB = dB 260 wC = dC 261 262 # Multiply and expand. 263 mat_R1 = multiply.outer( R1 * relax_time, m_R1 ) 264 mat_r1rho_prime = multiply.outer( r1rho_prime * relax_time, m_r1rho_prime ) 265 266 mat_wA = multiply.outer( wA * relax_time, m_wA ) 267 mat_wB = multiply.outer( wB * relax_time, m_wB ) 268 mat_wC = multiply.outer( wC * relax_time, m_wC ) 269 mat_w1 = multiply.outer( w1 * relax_time, m_w1 ) 270 271 mat_k_AB = multiply.outer( k_AB * relax_time, m_k_AB ) 272 mat_k_BA = multiply.outer( k_BA * relax_time, m_k_BA ) 273 mat_k_BC = multiply.outer( k_BC * relax_time, m_k_BC ) 274 275 mat_k_CB = multiply.outer( k_CB * relax_time, m_k_CB ) 276 mat_k_AC = multiply.outer( k_AC * relax_time, m_k_AC ) 277 mat_k_CA = multiply.outer( k_CA * relax_time, m_k_CA ) 278 279 # Collect matrix. 280 matrix = (mat_R1 + mat_r1rho_prime 281 + mat_wA + mat_wB + mat_wC + mat_w1 282 + mat_k_AB + mat_k_BA + mat_k_BC 283 + mat_k_CB + mat_k_AC + mat_k_CA ) 284 285 # Return the matrix. 286 return matrix

287 288

289 -def ns_r1rho_3site(M0=None, M0_T=None, r1rho_prime=None, omega=None, offset=None, r1=0.0, pA=None, pB=None, dw_AB=None, dw_BC=None, kex_AB=None, kex_BC=None, kex_AC=None, spin_lock_fields=None, relax_time=None, inv_relax_time=None, back_calc=None, num_points=None):

290 """The 3-site numerical solution to the Bloch-McConnell equation for R1rho data. 291 292 This function calculates and stores the R1rho values. 293 294 295 @keyword M0: This is a vector that contains the initial magnetizations corresponding to the A and B state transverse magnetizations. 296 @type M0: numpy float array of rank [NE][NS][NM][NO][ND][9][1] 297 @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. 298 @keyword r1rho_prime: The R1rho_prime parameter value (R1rho with no exchange). 299 @type r1rho_prime: numpy float array of rank [NE][NS][NM][NO][ND][1][9] 300 @keyword omega: The chemical shift for the spin in rad/s. 301 @type omega: numpy float array of rank [NS][NM][NO][ND] 302 @keyword offset: The spin-lock offsets for the data. 303 @type offset: numpy float array of rank [NS][NM][NO][ND] 304 @keyword r1: The R1 relaxation rate. 305 @type r1: numpy float array of rank [NS][NM][NO][ND] 306 @keyword pA: The population of state A. 307 @type pA: float 308 @keyword pB: The population of state B. 309 @type pB: float 310 @keyword dw_AB: The chemical exchange difference between states A and B in rad/s. 311 @type dw_AB: numpy float array of rank [NS][NM][NO][ND] 312 @keyword dw_BC: The chemical exchange difference between states B and C in rad/s. 313 @type dw_BC: numpy float array of rank [NS][NM][NO][ND] 314 @keyword kex_AB: The exchange rate between sites A and B for 3-site exchange with kex_AB = k_AB + k_BA (rad.s^-1) 315 @type kex_AB: float 316 @keyword kex_BC: The exchange rate between sites A and C for 3-site exchange with kex_AC = k_AC + k_CA (rad.s^-1) 317 @type kex_BC: float 318 @keyword kex_AC: The exchange rate between sites B and C for 3-site exchange with kex_BC = k_BC + k_CB (rad.s^-1) 319 @type kex_AC: float 320 @keyword spin_lock_fields: The R1rho spin-lock field strengths (in rad.s^-1). 321 @type spin_lock_fields: numpy float array of rank [NS][NM][NO][ND] 322 @keyword relax_time: The total relaxation time period for each spin-lock field strength (in seconds). 323 @type relax_time: numpy float array of rank [NS][NM][NO][ND] 324 @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. 325 @type inv_relax_time: numpy float array of rank [NS][NM][NO][ND] 326 @keyword back_calc: The array for holding the back calculated R2eff values. Each element corresponds to one of the CPMG nu1 frequencies. 327 @type back_calc: numpy float array of rank [NS][NM][NO][ND] 328 @keyword num_points: The number of points on the dispersion curve, equal to the length of the tcp and back_calc arguments. 329 @type num_points: numpy int array of rank [NS][NM][NO] 330 """ 331 332 # Once off parameter conversions. 333 dw_AC = dw_AB + dw_BC 334 pC = 1.0 - pA - pB 335 pA_pB = pA + pB 336 pA_pC = pA + pC 337 pB_pC = pB + pC 338 k_BA = pA * kex_AB / pA_pB 339 k_AB = pB * kex_AB / pA_pB 340 if pB_pC != 0.0: 341 k_CB = pB * kex_BC / pB_pC 342 k_BC = pC * kex_BC / pB_pC 343 elif pB == 0.0 and pC == 0.0: 344 k_CB = 0.0 345 k_BC = 0.0 346 elif pB == 0.0: 347 k_CB = 0.0 348 k_BC = 1e100 349 elif pC == 0.0: 350 k_CB = 1e100 351 k_BC = 0.0 352 else: 353 k_CB = 1e100 354 k_BC = 1e100 355 k_CA = pA * kex_AC / pA_pC 356 k_AC = pC * kex_AC / pA_pC 357 358 # Extract shape of experiment. 359 NE, NS, NM, NO = num_points.shape 360 361 # The matrix that contains all the contributions to the evolution, i.e. relaxation, exchange and chemical shift evolution. 362 R_mat = rr1rho_3d_3site_rankN(R1=r1, r1rho_prime=r1rho_prime, omega=omega, offset=offset, dw_AB=dw_AB, dw_AC=dw_AC, w1=spin_lock_fields, k_AB=k_AB, k_BA=k_BA, k_BC=k_BC, k_CB=k_CB, k_AC=k_AC, k_CA=k_CA, relax_time=relax_time) 363 364 # This matrix is a propagator that will evolve the magnetization with the matrix R. 365 Rexpo_mat = matrix_exponential(R_mat) 366 367 # Magnetization evolution. 368 Rexpo_M0_mat = einsum('...ij, ...jk', Rexpo_mat, M0) 369 370 # Magnetization evolution, which include all dimensions. 371 MA_mat = einsum('...ij, ...jk', M0_T, Rexpo_M0_mat)[:, :, :, :, :, 0, 0] 372 373 # Insert safe checks. 374 if min(MA_mat) < 0.0: 375 mask_min_MA_mat = masked_less(MA_mat, 0.0) 376 # Fill with high values. 377 MA_mat[mask_min_MA_mat.mask] = 1e100 378 379 # Do back calculation. 380 back_calc[:] = -inv_relax_time * log(MA_mat) 381 382 # Catch errors, taking a sum over array is the fastest way to check for 383 # +/- inf (infinity) and nan (not a number). 384 if not isfinite(sum(back_calc)): 385 # Replaces nan, inf, etc. with fill value. 386 fix_invalid(back_calc, copy=False, fill_value=1e100)

387

Source Code for Module lib.dispersion.ns_r1rho_3site