lib.dispersion.ns_mmq

1 ############################################################################### 2 # # 3 # Copyright (C) 2013 Mathilde Lescanne # 4 # Copyright (C) 2013 Dominique Marion # 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 numeric solution for the 2-site Bloch-McConnell equations for MMQ CPMG data, the U{NS MMQ 2-site<http://wiki.nmr-relax.com/NS_MMQ_2-site>} model. 27 28 Description 29 =========== 30 31 This handles proton-heteronuclear SQ, ZQ, DQ and MQ CPMG data. 32 33 34 References 35 ========== 36 37 It uses the equations of: 38 39 - Dmitry M. Korzhnev, Philipp Neudecker, Anthony Mittermaier, Vladislav Yu. Orekhov, and Lewis E. Kay (2005). Multiple-site exchange in proteins studied with a suite of six NMR relaxation dispersion experiments: An application to the folding of a Fyn SH3 domain mutant. I{J. Am. Chem. Soc.}, B{127}, 15602-15611. (U{DOI: 10.1021/ja054550e<http://dx.doi.org/10.1021/ja054550e>}). 40 41 42 Links 43 ===== 44 45 More information on the NS MMQ 2-site model can be found in the: 46 47 - U{relax wiki<http://wiki.nmr-relax.com/NS_MMQ_2-site>}, 48 - U{relax manual<http://www.nmr-relax.com/manual/The_NS_MMQ_2_site_model.html>}, 49 - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#NS_MMQ_2-site>}. 50 """ 51 52 # Python module imports. 53 from math import floor 54 from numpy import array, conj, complex128, dot, einsum, float64, log, multiply 55 from numpy.linalg import matrix_power 56 57 # relax module imports. 58 from lib.float import isNaN 59 from lib.dispersion.matrix_exponential import matrix_exponential 60 61 # Repetitive calculations (to speed up calculations). 62 m_r20a = array([ 63 [-1, 0], 64 [ 0, 0]], float64) 65 66 m_r20b = array([ 67 [ 0, 0], 68 [ 0, -1]], float64) 69 70 m_k_AB = array([ 71 [-1, 0], 72 [ 1, 0]], float64) 73 74 m_k_BA = array([ 75 [ 0, 1], 76 [ 0, -1]], float64) 77 78 m_dw = array([ 79 [ 0, 0], 80 [ 0, 1]], float64) 81 82

83 -def rmmq_2site_rankN(R20A=None, R20B=None, dw=None, k_AB=None, k_BA=None, tcp=None):

84 """The Bloch-McConnell matrix for 2-site exchange, for rank [NE][NS][NM][NO][ND][2][2]. 85 86 @keyword R20A: The transverse, spin-spin relaxation rate for state A. 87 @type R20A: numpy float array of rank [NE][NS][NM][NO][ND] 88 @keyword R20B: The transverse, spin-spin relaxation rate for state B. 89 @type R20B: numpy float array of rank [NE][NS][NM][NO][ND] 90 @keyword dw: The combined chemical exchange difference parameters between states A and B in rad/s. This can be any combination of dw and dwH. 91 @type dw: numpy float array of rank [NE][NS][NM][NO][ND] 92 @keyword k_AB: The rate of exchange from site A to B (rad/s). 93 @type k_AB: float 94 @keyword k_BA: The rate of exchange from site B to A (rad/s). 95 @type k_BA: float 96 @keyword tcp: The tau_CPMG times (1 / 4.nu1). 97 @type tcp: numpy float array of rank [NE][NS][NM][NO][ND] 98 @return: The relaxation matrix. 99 @rtype: numpy float array of rank [NE][NS][NM][NO][ND][2][2] 100 """ 101 102 # Pre-multiply with tcp. 103 r20a_tcp = R20A * tcp 104 r20b_tcp = R20B * tcp 105 k_AB_tcp = k_AB * tcp 106 k_BA_tcp = k_BA * tcp 107 # Complex dw. 108 dw_tcp_C = dw * tcp * 1j 109 110 # Fill in the elements. 111 #matrix[0, 0] = -k_AB - R20A 112 #matrix[0, 1] = k_BA 113 #matrix[1, 0] = k_AB 114 #matrix[1, 1] = -k_BA + 1.j*dw - R20B 115 116 # Multiply and expand. 117 m_r20a_tcp = multiply.outer( r20a_tcp, m_r20a ) 118 m_r20b_tcp = multiply.outer( r20b_tcp, m_r20b ) 119 120 # Multiply and expand. 121 m_k_AB_tcp = multiply.outer( k_AB_tcp, m_k_AB ) 122 m_k_BA_tcp = multiply.outer( k_BA_tcp, m_k_BA ) 123 124 # Multiply and expand. 125 m_dw_tcp_C = multiply.outer( dw_tcp_C, m_dw ) 126 127 # Collect matrix. 128 matrix = (m_r20a_tcp + m_r20b_tcp + m_k_AB_tcp + m_k_BA_tcp + m_dw_tcp_C) 129 130 return matrix

131 132

133 -def r2eff_ns_mmq_2site_mq(M0=None, F_vector=array([1, 0], float64), R20A=None, R20B=None, pA=None, dw=None, dwH=None, kex=None, inv_tcpmg=None, tcp=None, back_calc=None, num_points=None, power=None):

134 """The 2-site numerical solution to the Bloch-McConnell equation for MQ data. 135 136 The notation used here comes from: 137 138 - Dmitry M. Korzhnev, Philipp Neudecker, Anthony Mittermaier, Vladislav Yu. Orekhov, and Lewis E. Kay (2005). Multiple-site exchange in proteins studied with a suite of six NMR relaxation dispersion experiments: An application to the folding of a Fyn SH3 domain mutant. J. Am. Chem. Soc., 127, 15602-15611. (doi: http://dx.doi.org/10.1021/ja054550e). 139 140 and: 141 142 - Dmitry M. Korzhnev, Philipp Neudecker, Anthony Mittermaier, Vladislav Yu. Orekhov, and Lewis E. Kay (2005). Multiple-site exchange in proteins studied with a suite of six NMR relaxation dispersion experiments: An application to the folding of a Fyn SH3 domain mutant. J. Am. Chem. Soc., 127, 15602-15611. (doi: http://dx.doi.org/10.1021/ja054550e). 143 144 This function calculates and stores the R2eff values. 145 146 147 @keyword M0: This is a vector that contains the initial magnetizations corresponding to the A and B state transverse magnetizations. 148 @type M0: numpy float64, rank-1, 7D array 149 @keyword F_vector: The observable magnitisation vector. This defaults to [1, 0] for X observable magnitisation. 150 @type F_vector: numpy rank-1, 2D float64 array 151 @keyword R20A: The transverse, spin-spin relaxation rate for state A. 152 @type R20A: numpy float array of rank [NS][NM][NO][ND] 153 @keyword R20B: The transverse, spin-spin relaxation rate for state B. 154 @type R20B: numpy float array of rank [NS][NM][NO][ND] 155 @keyword pA: The population of state A. 156 @type pA: float 157 @keyword dw: The chemical exchange difference between states A and B in rad/s. 158 @type dw: numpy float array of rank [NS][NM][NO][ND] 159 @keyword dwH: The proton chemical exchange difference between states A and B in rad/s. 160 @type dwH: numpy float array of rank [NS][NM][NO][ND] 161 @keyword kex: The kex parameter value (the exchange rate in rad/s). 162 @type kex: float 163 @keyword inv_tcpmg: The inverse of the total duration of the CPMG element (in inverse seconds). 164 @type inv_tcpmg: numpy float array of rank [NS][NM][NO][ND] 165 @keyword tcp: The tau_CPMG times (1 / 4.nu1). 166 @type tcp: numpy float array of rank [NS][NM][NO][ND] 167 @keyword back_calc: The array for holding the back calculated R2eff values. Each element corresponds to one of the CPMG nu1 frequencies. 168 @type back_calc: numpy float array of rank [NS][NM][NO][ND] 169 @keyword num_points: The number of points on the dispersion curve, equal to the length of the tcp and back_calc arguments. 170 @type num_points: numpy int array of rank [NS][NM][NO] 171 @keyword power: The matrix exponential power array. 172 @type power: numpy int array of rank [NS][NM][NO][ND] 173 """ 174 175 # Once off parameter conversions. 176 pB = 1.0 - pA 177 k_BA = pA * kex 178 k_AB = pB * kex 179 180 # This is a vector that contains the initial magnetizations corresponding to the A and B state transverse magnetizations. 181 M0[0] = pA 182 M0[1] = pB 183 184 # Extract shape of experiment. 185 NS, NM, NO = num_points.shape 186 187 # Populate the m1 and m2 matrices (only once per function call for speed). 188 # D+ matrix component. 189 m1_mat = rmmq_2site_rankN(R20A=R20A, R20B=R20B, dw=-dw - dwH, k_AB=k_AB, k_BA=k_BA, tcp=tcp) 190 # Z- matrix component. 191 m2_mat = rmmq_2site_rankN(R20A=R20A, R20B=R20B, dw=dw - dwH, k_AB=k_AB, k_BA=k_BA, tcp=tcp) 192 193 # The M1 and M2 matrices. 194 # Equivalent to D+. 195 M1_mat = matrix_exponential(m1_mat, dtype=complex128) 196 # Equivalent to Z-. 197 M2_mat = matrix_exponential(m2_mat, dtype=complex128) 198 199 # The complex conjugates M1* and M2* 200 # Equivalent to D+*. 201 M1_star_mat = conj(M1_mat) 202 # Equivalent to Z-*. 203 M2_star_mat = conj(M2_mat) 204 205 # Repetitive dot products (minimised for speed). 206 M1_M2_mat = einsum('...ij, ...jk', M1_mat, M2_mat) 207 M2_M1_mat = einsum('...ij, ...jk', M2_mat, M1_mat) 208 M1_M2_M2_M1_mat = einsum('...ij, ...jk', M1_M2_mat, M2_M1_mat) 209 M2_M1_M1_M2_mat = einsum('...ij, ...jk', M2_M1_mat, M1_M2_mat) 210 M1_M2_star_mat = einsum('...ij, ...jk', M1_star_mat, M2_star_mat) 211 M2_M1_star_mat = einsum('...ij, ...jk', M2_star_mat, M1_star_mat) 212 M1_M2_M2_M1_star_mat = einsum('...ij, ...jk', M1_M2_star_mat, M2_M1_star_mat) 213 M2_M1_M1_M2_star_mat = einsum('...ij, ...jk', M2_M1_star_mat, M1_M2_star_mat) 214 215 # Loop over spins. 216 for si in range(NS): 217 # Loop over the spectrometer frequencies. 218 for mi in range(NM): 219 # Loop over offsets: 220 for oi in range(NO): 221 num_points_i = num_points[si, mi, oi] 222 223 # Loop over the time points, back calculating the R2eff values. 224 for i in range(num_points_i): 225 # Extract data from array. 226 power_i = int(power[si, mi, oi, i]) 227 M1_M2_i = M1_M2_mat[si, mi, oi, i] 228 M1_M2_star_i = M1_M2_star_mat[si, mi, oi, i] 229 M2_M1_i = M2_M1_mat[si, mi, oi, i] 230 M2_M1_star_i = M2_M1_star_mat[si, mi, oi, i] 231 M1_M2_M2_M1_i = M1_M2_M2_M1_mat[si, mi, oi, i] 232 M2_M1_M1_M2_star_i = M2_M1_M1_M2_star_mat[si, mi, oi, i] 233 M2_M1_M1_M2_i = M2_M1_M1_M2_mat[si, mi, oi, i] 234 M1_M2_M2_M1_star_i = M1_M2_M2_M1_star_mat[si, mi, oi, i] 235 236 # Special case of 1 CPMG block - the power is zero. 237 if power_i == 1: 238 # M1.M2. 239 A = M1_M2_i 240 241 # M1*.M2*. 242 B = M1_M2_star_i 243 244 # M2.M1. 245 C = M2_M1_i 246 247 # M2*.M1*. 248 D = M2_M1_star_i 249 250 # Matrices for even number of CPMG blocks. 251 elif power_i % 2 == 0: 252 # The power factor (only calculate once). 253 fact = int(floor(power_i / 2)) 254 255 # (M1.M2.M2.M1)^(n/2). 256 A = matrix_power(M1_M2_M2_M1_i, fact) 257 258 # (M2*.M1*.M1*.M2*)^(n/2). 259 B = matrix_power(M2_M1_M1_M2_star_i, fact) 260 261 # (M2.M1.M1.M2)^(n/2). 262 C = matrix_power(M2_M1_M1_M2_i, fact) 263 264 # (M1*.M2*.M2*.M1*)^(n/2). 265 D = matrix_power(M1_M2_M2_M1_star_i, fact) 266 267 # Matrices for odd number of CPMG blocks. 268 else: 269 # The power factor (only calculate once). 270 fact = int(floor((power_i - 1) / 2)) 271 272 # (M1.M2.M2.M1)^((n-1)/2).M1.M2. 273 A = matrix_power(M1_M2_M2_M1_i, fact) 274 A = dot(A, M1_M2_i) 275 276 # (M1*.M2*.M2*.M1*)^((n-1)/2).M1*.M2*. 277 B = matrix_power(M1_M2_M2_M1_star_i, fact) 278 B = dot(B, M1_M2_star_i) 279 280 # (M2.M1.M1.M2)^((n-1)/2).M2.M1. 281 C = matrix_power(M2_M1_M1_M2_i, fact) 282 C = dot(C, M2_M1_i) 283 284 # (M2*.M1*.M1*.M2*)^((n-1)/2).M2*.M1*. 285 D = matrix_power(M2_M1_M1_M2_star_i, fact) 286 D = dot(D, M2_M1_star_i) 287 288 # The next lines calculate the R2eff using a two-point approximation, i.e. assuming that the decay is mono-exponential. 289 A_B = dot(A, B) 290 C_D = dot(C, D) 291 Mx = dot(dot(F_vector, (A_B + C_D)), M0) 292 Mx = Mx.real / 2.0 293 if Mx <= 0.0 or isNaN(Mx): 294 back_calc[si, mi, oi, i] = 1e99 295 else: 296 back_calc[si, mi, oi, i]= -inv_tcpmg[si, mi, oi, i] * log(Mx / pA)

297 298

299 -def r2eff_ns_mmq_2site_sq_dq_zq(M0=None, F_vector=array([1, 0], float64), R20A=None, R20B=None, pA=None, dw=None, dwH=None, kex=None, inv_tcpmg=None, tcp=None, back_calc=None, num_points=None, power=None):

300 """The 2-site numerical solution to the Bloch-McConnell equation for SQ, ZQ, and DQ data. 301 302 The notation used here comes from: 303 304 - Dmitry M. Korzhnev, Philipp Neudecker, Anthony Mittermaier, Vladislav Yu. Orekhov, and Lewis E. Kay (2005). Multiple-site exchange in proteins studied with a suite of six NMR relaxation dispersion experiments: An application to the folding of a Fyn SH3 domain mutant. J. Am. Chem. Soc., 127, 15602-15611. (doi: http://dx.doi.org/10.1021/ja054550e). 305 306 This function calculates and stores the R2eff values. 307 308 309 @keyword M0: This is a vector that contains the initial magnetizations corresponding to the A and B state transverse magnetizations. 310 @type M0: numpy float64, rank-1, 7D array 311 @keyword F_vector: The observable magnitisation vector. This defaults to [1, 0] for X observable magnitisation. 312 @type F_vector: numpy rank-1, 2D float64 array 313 @keyword R20A: The transverse, spin-spin relaxation rate for state A. 314 @type R20A: numpy float array of rank [NS][NM][NO][ND] 315 @keyword R20B: The transverse, spin-spin relaxation rate for state B. 316 @type R20B: numpy float array of rank [NS][NM][NO][ND] 317 @keyword pA: The population of state A. 318 @type pA: float 319 @keyword dw: The combined chemical exchange difference between states A and B in rad/s. It should be set to dwH for 1H SQ data, dw for heteronuclear SQ data, dwH-dw for ZQ data, and dwH+dw for DQ data. 320 @type dw: numpy float array of rank [NS][NM][NO][ND] 321 @keyword dwH: Unused - this is simply to match the r2eff_ns_mmq_2site_mq() function arguments. 322 @type dwH: numpy float array of rank [NS][NM][NO][ND] 323 @keyword kex: The kex parameter value (the exchange rate in rad/s). 324 @type kex: float 325 @keyword inv_tcpmg: The inverse of the total duration of the CPMG element (in inverse seconds). 326 @type inv_tcpmg: numpy float array of rank [NS][NM][NO][ND] 327 @keyword tcp: The tau_CPMG times (1 / 4.nu1). 328 @type tcp: numpy float array of rank [NS][NM][NO][ND] 329 @keyword back_calc: The array for holding the back calculated R2eff values. Each element corresponds to one of the CPMG nu1 frequencies. 330 @type back_calc: numpy float array of rank [NS][NM][NO][ND] 331 @keyword num_points: The number of points on the dispersion curve, equal to the length of the tcp and back_calc arguments. 332 @type num_points: numpy int array of rank [NS][NM][NO] 333 @keyword power: The matrix exponential power array. 334 @type power: numpy int array of rank [NS][NM][NO][ND] 335 """ 336 337 # Once off parameter conversions. 338 pB = 1.0 - pA 339 k_BA = pA * kex 340 k_AB = pB * kex 341 342 # This is a vector that contains the initial magnetizations corresponding to the A and B state transverse magnetizations. 343 M0[0] = pA 344 M0[1] = pB 345 346 # Extract shape of experiment. 347 NS, NM, NO = num_points.shape 348 349 # Populate the m1 and m2 matrices (only once per function call for speed). 350 m1_mat = rmmq_2site_rankN(R20A=R20A, R20B=R20B, dw=dw, k_AB=k_AB, k_BA=k_BA, tcp=tcp) 351 m2_mat = rmmq_2site_rankN(R20A=R20A, R20B=R20B, dw=-dw, k_AB=k_AB, k_BA=k_BA, tcp=tcp) 352 353 # The A+/- matrices. 354 A_pos_mat = matrix_exponential(m1_mat, dtype=complex128) 355 A_neg_mat = matrix_exponential(m2_mat, dtype=complex128) 356 357 # The evolution for one n. 358 evol_block_mat = einsum('...ij, ...jk', A_neg_mat, A_pos_mat) 359 evol_block_mat = einsum('...ij, ...jk', A_neg_mat, evol_block_mat) 360 evol_block_mat = einsum('...ij, ...jk', A_pos_mat, evol_block_mat) 361 362 # Loop over spins. 363 for si in range(NS): 364 # Loop over the spectrometer frequencies. 365 for mi in range(NM): 366 # Loop over offsets: 367 for oi in range(NO): 368 # Extract number of points. 369 num_points_i = num_points[si, mi, oi] 370 371 # Loop over the time points, back calculating the R2eff values. 372 for i in range(num_points_i): 373 # Extract data from array. 374 power_i = int(power[si, mi, oi, i]) 375 evol_block_i = evol_block_mat[si, mi, oi, i] 376 377 # The full evolution. 378 evol = matrix_power(evol_block_i, power_i) 379 380 # The next lines calculate the R2eff using a two-point approximation, i.e. assuming that the decay is mono-exponential. 381 Mx = dot(F_vector, dot(evol, M0)) 382 Mx = Mx.real 383 if Mx <= 0.0 or isNaN(Mx): 384 back_calc[si, mi, oi, i] = 1e99 385 else: 386 back_calc[si, mi, oi, i] = -inv_tcpmg[si, mi, oi, i] * log(Mx / pA)

387

Source Code for Module lib.dispersion.ns_mmq_2site