Package lib :: Package dispersion :: Module ns_mmq_3site
[hide private]
[frames] | no frames]

Source Code for Module lib.dispersion.ns_mmq_3site

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2013 Mathilde Lescanne                                        # 
  4  # Copyright (C) 2013 Dominique Marion                                         # 
  5  # Copyright (C) 2013-2014 Edward d'Auvergne                                   # 
  6  #                                                                             # 
  7  # This file is part of the program relax (http://www.nmr-relax.com).          # 
  8  #                                                                             # 
  9  # This program is free software: you can redistribute it and/or modify        # 
 10  # it under the terms of the GNU General Public License as published by        # 
 11  # the Free Software Foundation, either version 3 of the License, or           # 
 12  # (at your option) any later version.                                         # 
 13  #                                                                             # 
 14  # This program is distributed in the hope that it will be useful,             # 
 15  # but WITHOUT ANY WARRANTY; without even the implied warranty of              # 
 16  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               # 
 17  # GNU General Public License for more details.                                # 
 18  #                                                                             # 
 19  # You should have received a copy of the GNU General Public License           # 
 20  # along with this program.  If not, see <http://www.gnu.org/licenses/>.       # 
 21  #                                                                             # 
 22  ############################################################################### 
 23   
 24  # Module docstring. 
 25  """The numeric solution for the 3-site Bloch-McConnell equations for MMQ CPMG data, the U{NS MMQ 3-site linear<http://wiki.nmr-relax.com/NS_MMQ_3-site_linear>} and U{NS MMQ 3-site<http://wiki.nmr-relax.com/NS_MMQ_3-site>} models. 
 26   
 27  Description 
 28  =========== 
 29   
 30  This handles proton-heteronuclear SQ, ZQ, DQ and MQ CPMG data. 
 31   
 32   
 33  References 
 34  ========== 
 35   
 36  It uses the equations of: 
 37   
 38      - 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>}). 
 39   
 40   
 41  Links 
 42  ===== 
 43   
 44  More information on the NS MMQ 3-site linear model can be found in the: 
 45   
 46      - U{relax wiki<http://wiki.nmr-relax.com/NS_MMQ_3-site_linear>}, 
 47      - U{relax manual<http://www.nmr-relax.com/manual/NS_MMQ_3_site_linear_model.html>}, 
 48      - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#NS_MMQ_3-site_linear>}. 
 49   
 50  More information on the NS MMQ 3-site model can be found in the: 
 51   
 52      - U{relax wiki<http://wiki.nmr-relax.com/NS_MMQ_3-site>}, 
 53      - U{relax manual<http://www.nmr-relax.com/manual/NS_MMQ_3_site_model.html>}, 
 54      - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#NS_MMQ_3-site>}. 
 55  """ 
 56   
 57  # Python module imports. 
 58  from math import floor 
 59  from numpy import array, conj, dot, float64, log 
 60   
 61  # relax module imports. 
 62  from lib.float import isNaN 
 63  from lib.linear_algebra.matrix_exponential import matrix_exponential 
 64  from lib.linear_algebra.matrix_power import square_matrix_power 
 65   
 66   
67 -def populate_matrix(matrix=None, R20A=None, R20B=None, R20C=None, dw_AB=None, dw_AC=None, k_AB=None, k_BA=None, k_BC=None, k_CB=None, k_AC=None, k_CA=None):
68 """The Bloch-McConnell matrix for 3-site exchange. 69 70 @keyword matrix: The matrix to populate. 71 @type matrix: numpy rank-2, 3D complex64 array 72 @keyword R20A: The transverse, spin-spin relaxation rate for state A. 73 @type R20A: float 74 @keyword R20B: The transverse, spin-spin relaxation rate for state B. 75 @type R20B: float 76 @keyword R20C: The transverse, spin-spin relaxation rate for state C. 77 @type R20C: float 78 @keyword dw_AB: The combined chemical exchange difference parameters between states A and B in rad/s. This can be any combination of dw and dwH. 79 @type dw_AB: float 80 @keyword dw_AC: The combined chemical exchange difference parameters between states A and C in rad/s. This can be any combination of dw and dwH. 81 @type dw_AC: float 82 @keyword k_AB: The rate of exchange from site A to B (rad/s). 83 @type k_AB: float 84 @keyword k_BA: The rate of exchange from site B to A (rad/s). 85 @type k_BA: float 86 @keyword k_BC: The rate of exchange from site B to C (rad/s). 87 @type k_BC: float 88 @keyword k_CB: The rate of exchange from site C to B (rad/s). 89 @type k_CB: float 90 @keyword k_AC: The rate of exchange from site A to C (rad/s). 91 @type k_AC: float 92 @keyword k_CA: The rate of exchange from site C to A (rad/s). 93 @type k_CA: float 94 """ 95 96 # The first row. 97 matrix[0, 0] = -k_AB - k_AC - R20A 98 matrix[0, 1] = k_BA 99 matrix[0, 2] = k_CA 100 101 # The second row. 102 matrix[1, 0] = k_AB 103 matrix[1, 1] = -k_BA - k_BC + 1.j*dw_AB - R20B 104 matrix[1, 2] = k_CB 105 106 # The third row. 107 matrix[2, 0] = k_AC 108 matrix[2, 1] = k_BC 109 matrix[2, 2] = -k_CB - k_CA + 1.j*dw_AC - R20C
110 111
112 -def r2eff_ns_mmq_3site_mq(M0=None, F_vector=array([1, 0, 0], float64), m1=None, m2=None, R20A=None, R20B=None, R20C=None, pA=None, pB=None, pC=None, dw_AB=None, dw_AC=None, dwH_AB=None, dwH_AC=None, k_AB=None, k_BA=None, k_BC=None, k_CB=None, k_AC=None, k_CA=None, inv_tcpmg=None, tcp=None, back_calc=None, num_points=None, power=None):
113 """The 3-site numerical solution to the Bloch-McConnell equation for MQ data. 114 115 The notation used here comes from: 116 117 - 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). 118 119 and: 120 121 - 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). 122 123 This function calculates and stores the R2eff values. 124 125 126 @keyword M0: This is a vector that contains the initial magnetizations corresponding to the A and B state transverse magnetizations. 127 @type M0: numpy float64, rank-1, 7D array 128 @keyword F_vector: The observable magnitisation vector. This defaults to [1, 0] for X observable magnitisation. 129 @type F_vector: numpy rank-1, 3D float64 array 130 @keyword m1: A complex numpy matrix to be populated. 131 @type m1: numpy rank-2, 3D complex64 array 132 @keyword m2: A complex numpy matrix to be populated. 133 @type m2: numpy rank-2, 3D complex64 array 134 @keyword R20A: The transverse, spin-spin relaxation rate for state A. 135 @type R20A: float 136 @keyword R20B: The transverse, spin-spin relaxation rate for state B. 137 @type R20B: float 138 @keyword R20C: The transverse, spin-spin relaxation rate for state C. 139 @type R20C: float 140 @keyword pA: The population of state A. 141 @type pA: float 142 @keyword pB: The population of state B. 143 @type pB: float 144 @keyword pC: The population of state C. 145 @type pC: float 146 @keyword dw_AB: The chemical exchange difference between states A and B in rad/s. 147 @type dw_AB: float 148 @keyword dw_AC: The chemical exchange difference between states A and C in rad/s. 149 @type dw_AC: float 150 @keyword dwH_AB: The proton chemical exchange difference between states A and B in rad/s. 151 @type dwH_AB: float 152 @keyword dwH_AC: The proton chemical exchange difference between states A and C in rad/s. 153 @type dwH_AC: float 154 @keyword k_AB: The rate of exchange from site A to B (rad/s). 155 @type k_AB: float 156 @keyword k_BA: The rate of exchange from site B to A (rad/s). 157 @type k_BA: float 158 @keyword k_BC: The rate of exchange from site B to C (rad/s). 159 @type k_BC: float 160 @keyword k_CB: The rate of exchange from site C to B (rad/s). 161 @type k_CB: float 162 @keyword k_AC: The rate of exchange from site A to C (rad/s). 163 @type k_AC: float 164 @keyword k_CA: The rate of exchange from site C to A (rad/s). 165 @type k_CA: float 166 @keyword inv_tcpmg: The inverse of the total duration of the CPMG element (in inverse seconds). 167 @type inv_tcpmg: float 168 @keyword tcp: The tau_CPMG times (1 / 4.nu1). 169 @type tcp: numpy rank-1 float array 170 @keyword back_calc: The array for holding the back calculated R2eff values. Each element corresponds to one of the CPMG nu1 frequencies. 171 @type back_calc: numpy rank-1 float array 172 @keyword num_points: The number of points on the dispersion curve, equal to the length of the tcp and back_calc arguments. 173 @type num_points: int 174 @keyword power: The matrix exponential power array. 175 @type power: numpy int16, rank-1 array 176 """ 177 178 # Populate the m1 and m2 matrices (only once per function call for speed). 179 populate_matrix(matrix=m1, R20A=R20A, R20B=R20B, R20C=R20C, dw_AB=-dw_AB-dwH_AB, dw_AC=-dw_AC-dwH_AC, k_AB=k_AB, k_BA=k_BA, k_BC=k_BC, k_CB=k_CB, k_AC=k_AC, k_CA=k_CA) # D+ matrix component. 180 populate_matrix(matrix=m2, R20A=R20A, R20B=R20B, R20C=R20C, dw_AB=dw_AB-dwH_AB, dw_AC=dw_AC-dwH_AC, k_AB=k_AB, k_BA=k_BA, k_BC=k_BC, k_CB=k_CB, k_AC=k_AC, k_CA=k_CA) # Z- matrix component. 181 182 # Loop over the time points, back calculating the R2eff values. 183 for i in range(num_points): 184 # The M1 and M2 matrices. 185 M1 = matrix_exponential(m1*tcp[i]) # Equivalent to D+. 186 M2 = matrix_exponential(m2*tcp[i]) # Equivalent to Z-. 187 188 # The complex conjugates M1* and M2* 189 M1_star = conj(M1) # Equivalent to D+*. 190 M2_star = conj(M2) # Equivalent to Z-*. 191 192 # Repetitive dot products (minimised for speed). 193 M1_M2 = dot(M1, M2) 194 M2_M1 = dot(M2, M1) 195 M1_M2_M2_M1 = dot(M1_M2, M2_M1) 196 M2_M1_M1_M2 = dot(M2_M1, M1_M2) 197 M1_M2_star = dot(M1_star, M2_star) 198 M2_M1_star = dot(M2_star, M1_star) 199 M1_M2_M2_M1_star = dot(M1_M2_star, M2_M1_star) 200 M2_M1_M1_M2_star = dot(M2_M1_star, M1_M2_star) 201 202 # Special case of 1 CPMG block - the power is zero. 203 if power[i] == 1: 204 # M1.M2. 205 A = M1_M2 206 207 # M1*.M2*. 208 B = M1_M2_star 209 210 # M2.M1. 211 C = M2_M1 212 213 # M2*.M1*. 214 D = M2_M1_star 215 216 # Matrices for even number of CPMG blocks. 217 elif power[i] % 2 == 0: 218 # The power factor (only calculate once). 219 fact = int(floor(power[i] / 2)) 220 221 # (M1.M2.M2.M1)^(n/2). 222 A = square_matrix_power(M1_M2_M2_M1, fact) 223 224 # (M2*.M1*.M1*.M2*)^(n/2). 225 B = square_matrix_power(M2_M1_M1_M2_star, fact) 226 227 # (M2.M1.M1.M2)^(n/2). 228 C = square_matrix_power(M2_M1_M1_M2, fact) 229 230 # (M1*.M2*.M2*.M1*)^(n/2). 231 D = square_matrix_power(M1_M2_M2_M1_star, fact) 232 233 # Matrices for odd number of CPMG blocks. 234 else: 235 # The power factor (only calculate once). 236 fact = int(floor((power[i] - 1) / 2)) 237 238 # (M1.M2.M2.M1)^((n-1)/2).M1.M2. 239 A = square_matrix_power(M1_M2_M2_M1, fact) 240 A = dot(A, M1_M2) 241 242 # (M1*.M2*.M2*.M1*)^((n-1)/2).M1*.M2*. 243 B = square_matrix_power(M1_M2_M2_M1_star, fact) 244 B = dot(B, M1_M2_star) 245 246 # (M2.M1.M1.M2)^((n-1)/2).M2.M1. 247 C = square_matrix_power(M2_M1_M1_M2, fact) 248 C = dot(C, M2_M1) 249 250 # (M2*.M1*.M1*.M2*)^((n-1)/2).M2*.M1*. 251 D = square_matrix_power(M2_M1_M1_M2_star, fact) 252 D = dot(D, M2_M1_star) 253 254 # The next lines calculate the R2eff using a two-point approximation, i.e. assuming that the decay is mono-exponential. 255 A_B = dot(A, B) 256 C_D = dot(C, D) 257 Mx = dot(dot(F_vector, (A_B + C_D)), M0) 258 Mx = Mx.real / 2.0 259 if Mx <= 0.0 or isNaN(Mx): 260 back_calc[i] = 1e99 261 else: 262 back_calc[i]= -inv_tcpmg * log(Mx / pA)
263 264
265 -def r2eff_ns_mmq_3site_sq_dq_zq(M0=None, F_vector=array([1, 0, 0], float64), m1=None, m2=None, R20A=None, R20B=None, R20C=None, pA=None, pB=None, pC=None, dw_AB=None, dw_AC=None, dwH_AB=None, dwH_AC=None, k_AB=None, k_BA=None, k_BC=None, k_CB=None, k_AC=None, k_CA=None, inv_tcpmg=None, tcp=None, back_calc=None, num_points=None, power=None):
266 """The 3-site numerical solution to the Bloch-McConnell equation for SQ, ZQ, and DQ data. 267 268 The notation used here comes from: 269 270 - 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). 271 272 This function calculates and stores the R2eff values. 273 274 275 @keyword M0: This is a vector that contains the initial magnetizations corresponding to the A and B state transverse magnetizations. 276 @type M0: numpy float64, rank-1, 7D array 277 @keyword F_vector: The observable magnitisation vector. This defaults to [1, 0] for X observable magnitisation. 278 @type F_vector: numpy rank-1, 3D float64 array 279 @keyword m1: A complex numpy matrix to be populated. 280 @type m1: numpy rank-2, 3D complex64 array 281 @keyword m2: A complex numpy matrix to be populated. 282 @type m2: numpy rank-2, 3D complex64 array 283 @keyword R20A: The transverse, spin-spin relaxation rate for state A. 284 @type R20A: float 285 @keyword R20B: The transverse, spin-spin relaxation rate for state B. 286 @type R20B: float 287 @keyword R20C: The transverse, spin-spin relaxation rate for state C. 288 @type R20C: float 289 @keyword pA: The population of state A. 290 @type pA: float 291 @keyword pB: The population of state B. 292 @type pB: float 293 @keyword pC: The population of state C. 294 @type pC: float 295 @keyword dw_AB: 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. 296 @type dw_AB: float 297 @keyword dw_AC: The combined chemical exchange difference between states A and C 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. 298 @type dw_AC: float 299 @keyword dwH_AB: Unused - this is simply to match the r2eff_mmq_3site_mq() function arguments. 300 @type dwH_AB: float 301 @keyword dwH_AC: Unused - this is simply to match the r2eff_mmq_3site_mq() function arguments. 302 @type dwH_AC: float 303 @keyword k_AB: The rate of exchange from site A to B (rad/s). 304 @type k_AB: float 305 @keyword k_BA: The rate of exchange from site B to A (rad/s). 306 @type k_BA: float 307 @keyword inv_tcpmg: The inverse of the total duration of the CPMG element (in inverse seconds). 308 @type inv_tcpmg: float 309 @keyword tcp: The tau_CPMG times (1 / 4.nu1). 310 @type tcp: numpy rank-1 float array 311 @keyword back_calc: The array for holding the back calculated R2eff values. Each element corresponds to one of the CPMG nu1 frequencies. 312 @type back_calc: numpy rank-1 float array 313 @keyword num_points: The number of points on the dispersion curve, equal to the length of the tcp and back_calc arguments. 314 @type num_points: int 315 @keyword power: The matrix exponential power array. 316 @type power: numpy int16, rank-1 array 317 """ 318 319 # Populate the m1 and m2 matrices (only once per function call for speed). 320 populate_matrix(matrix=m1, R20A=R20A, R20B=R20B, R20C=R20C, dw_AB=dw_AB, dw_AC=dw_AC, k_AB=k_AB, k_BA=k_BA, k_BC=k_BC, k_CB=k_CB, k_AC=k_AC, k_CA=k_CA) 321 populate_matrix(matrix=m2, R20A=R20A, R20B=R20B, R20C=R20C, dw_AB=-dw_AB, dw_AC=-dw_AC, k_AB=k_AB, k_BA=k_BA, k_BC=k_BC, k_CB=k_CB, k_AC=k_AC, k_CA=k_CA) 322 323 # Loop over the time points, back calculating the R2eff values. 324 for i in range(num_points): 325 # The A+/- matrices. 326 A_pos = matrix_exponential(m1*tcp[i]) 327 A_neg = matrix_exponential(m2*tcp[i]) 328 329 # The evolution for one n. 330 evol_block = dot(A_pos, dot(A_neg, dot(A_neg, A_pos))) 331 332 # The full evolution. 333 evol = square_matrix_power(evol_block, power[i]) 334 335 # The next lines calculate the R2eff using a two-point approximation, i.e. assuming that the decay is mono-exponential. 336 Mx = dot(F_vector, dot(evol, M0)) 337 Mx = Mx.real 338 if Mx <= 0.0 or isNaN(Mx): 339 back_calc[i] = 1e99 340 else: 341 back_calc[i] = -inv_tcpmg * log(Mx / pA)
342