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

Source Code for Module lib.dispersion.ns_mmq_2site

  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 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. 
 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 2-site model can be found in the: 
 45   
 46      - U{relax wiki<http://wiki.nmr-relax.com/NS_MMQ_2-site>}, 
 47      - U{relax manual<http://www.nmr-relax.com/manual/NS_MMQ_2_site_model.html>}, 
 48      - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#NS_MMQ_2-site>}. 
 49  """ 
 50   
 51  # Python module imports. 
 52  from math import floor 
 53  from numpy import array, conj, dot, float64, log 
 54   
 55  # relax module imports. 
 56  from lib.float import isNaN 
 57  from lib.linear_algebra.matrix_exponential import matrix_exponential 
 58  from lib.linear_algebra.matrix_power import square_matrix_power 
 59   
 60   
61 -def populate_matrix(matrix=None, R20A=None, R20B=None, dw=None, k_AB=None, k_BA=None):
62 """The Bloch-McConnell matrix for 2-site exchange. 63 64 @keyword matrix: The matrix to populate. 65 @type matrix: numpy rank-2, 2D complex64 array 66 @keyword R20A: The transverse, spin-spin relaxation rate for state A. 67 @type R20A: float 68 @keyword R20B: The transverse, spin-spin relaxation rate for state B. 69 @type R20B: float 70 @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. 71 @type dw: float 72 @keyword k_AB: The rate of exchange from site A to B (rad/s). 73 @type k_AB: float 74 @keyword k_BA: The rate of exchange from site B to A (rad/s). 75 @type k_BA: float 76 """ 77 78 # Fill in the elements. 79 matrix[0, 0] = -k_AB - R20A 80 matrix[0, 1] = k_BA 81 matrix[1, 0] = k_AB 82 matrix[1, 1] = -k_BA + 1.j*dw - R20B
83 84
85 -def r2eff_ns_mmq_2site_mq(M0=None, F_vector=array([1, 0], float64), m1=None, m2=None, R20A=None, R20B=None, pA=None, pB=None, dw=None, dwH=None, k_AB=None, k_BA=None, inv_tcpmg=None, tcp=None, back_calc=None, num_points=None, power=None):
86 """The 2-site numerical solution to the Bloch-McConnell equation for MQ data. 87 88 The notation used here comes from: 89 90 - 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). 91 92 and: 93 94 - 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). 95 96 This function calculates and stores the R2eff values. 97 98 99 @keyword M0: This is a vector that contains the initial magnetizations corresponding to the A and B state transverse magnetizations. 100 @type M0: numpy float64, rank-1, 7D array 101 @keyword F_vector: The observable magnitisation vector. This defaults to [1, 0] for X observable magnitisation. 102 @type F_vector: numpy rank-1, 2D float64 array 103 @keyword m1: A complex numpy matrix to be populated. 104 @type m1: numpy rank-2, 2D complex64 array 105 @keyword m2: A complex numpy matrix to be populated. 106 @type m2: numpy rank-2, 2D complex64 array 107 @keyword R20A: The transverse, spin-spin relaxation rate for state A. 108 @type R20A: float 109 @keyword R20B: The transverse, spin-spin relaxation rate for state B. 110 @type R20B: float 111 @keyword pA: The population of state A. 112 @type pA: float 113 @keyword pB: The population of state B. 114 @type pB: float 115 @keyword dw: The chemical exchange difference between states A and B in rad/s. 116 @type dw: float 117 @keyword dwH: The proton chemical exchange difference between states A and B in rad/s. 118 @type dwH: float 119 @keyword k_AB: The rate of exchange from site A to B (rad/s). 120 @type k_AB: float 121 @keyword k_BA: The rate of exchange from site B to A (rad/s). 122 @type k_BA: float 123 @keyword inv_tcpmg: The inverse of the total duration of the CPMG element (in inverse seconds). 124 @type inv_tcpmg: float 125 @keyword tcp: The tau_CPMG times (1 / 4.nu1). 126 @type tcp: numpy rank-1 float array 127 @keyword back_calc: The array for holding the back calculated R2eff values. Each element corresponds to one of the CPMG nu1 frequencies. 128 @type back_calc: numpy rank-1 float array 129 @keyword num_points: The number of points on the dispersion curve, equal to the length of the tcp and back_calc arguments. 130 @type num_points: int 131 @keyword power: The matrix exponential power array. 132 @type power: numpy int16, rank-1 array 133 """ 134 135 # Populate the m1 and m2 matrices (only once per function call for speed). 136 populate_matrix(matrix=m1, R20A=R20A, R20B=R20B, dw=-dw-dwH, k_AB=k_AB, k_BA=k_BA) # D+ matrix component. 137 populate_matrix(matrix=m2, R20A=R20A, R20B=R20B, dw=dw-dwH, k_AB=k_AB, k_BA=k_BA) # Z- matrix component. 138 139 # Loop over the time points, back calculating the R2eff values. 140 for i in range(num_points): 141 # The M1 and M2 matrices. 142 M1 = matrix_exponential(m1*tcp[i]) # Equivalent to D+. 143 M2 = matrix_exponential(m2*tcp[i]) # Equivalent to Z-. 144 145 # The complex conjugates M1* and M2* 146 M1_star = conj(M1) # Equivalent to D+*. 147 M2_star = conj(M2) # Equivalent to Z-*. 148 149 # Repetitive dot products (minimised for speed). 150 M1_M2 = dot(M1, M2) 151 M2_M1 = dot(M2, M1) 152 M1_M2_M2_M1 = dot(M1_M2, M2_M1) 153 M2_M1_M1_M2 = dot(M2_M1, M1_M2) 154 M1_M2_star = dot(M1_star, M2_star) 155 M2_M1_star = dot(M2_star, M1_star) 156 M1_M2_M2_M1_star = dot(M1_M2_star, M2_M1_star) 157 M2_M1_M1_M2_star = dot(M2_M1_star, M1_M2_star) 158 159 # Special case of 1 CPMG block - the power is zero. 160 if power[i] == 1: 161 # M1.M2. 162 A = M1_M2 163 164 # M1*.M2*. 165 B = M1_M2_star 166 167 # M2.M1. 168 C = M2_M1 169 170 # M2*.M1*. 171 D = M2_M1_star 172 173 # Matrices for even number of CPMG blocks. 174 elif power[i] % 2 == 0: 175 # The power factor (only calculate once). 176 fact = int(floor(power[i] / 2)) 177 178 # (M1.M2.M2.M1)^(n/2). 179 A = square_matrix_power(M1_M2_M2_M1, fact) 180 181 # (M2*.M1*.M1*.M2*)^(n/2). 182 B = square_matrix_power(M2_M1_M1_M2_star, fact) 183 184 # (M2.M1.M1.M2)^(n/2). 185 C = square_matrix_power(M2_M1_M1_M2, fact) 186 187 # (M1*.M2*.M2*.M1*)^(n/2). 188 D = square_matrix_power(M1_M2_M2_M1_star, fact) 189 190 # Matrices for odd number of CPMG blocks. 191 else: 192 # The power factor (only calculate once). 193 fact = int(floor((power[i] - 1) / 2)) 194 195 # (M1.M2.M2.M1)^((n-1)/2).M1.M2. 196 A = square_matrix_power(M1_M2_M2_M1, fact) 197 A = dot(A, M1_M2) 198 199 # (M1*.M2*.M2*.M1*)^((n-1)/2).M1*.M2*. 200 B = square_matrix_power(M1_M2_M2_M1_star, fact) 201 B = dot(B, M1_M2_star) 202 203 # (M2.M1.M1.M2)^((n-1)/2).M2.M1. 204 C = square_matrix_power(M2_M1_M1_M2, fact) 205 C = dot(C, M2_M1) 206 207 # (M2*.M1*.M1*.M2*)^((n-1)/2).M2*.M1*. 208 D = square_matrix_power(M2_M1_M1_M2_star, fact) 209 D = dot(D, M2_M1_star) 210 211 # The next lines calculate the R2eff using a two-point approximation, i.e. assuming that the decay is mono-exponential. 212 A_B = dot(A, B) 213 C_D = dot(C, D) 214 Mx = dot(dot(F_vector, (A_B + C_D)), M0) 215 Mx = Mx.real / 2.0 216 if Mx <= 0.0 or isNaN(Mx): 217 back_calc[i] = 1e99 218 else: 219 back_calc[i]= -inv_tcpmg * log(Mx / pA)
220 221
222 -def r2eff_ns_mmq_2site_sq_dq_zq(M0=None, F_vector=array([1, 0], float64), m1=None, m2=None, R20A=None, R20B=None, pA=None, pB=None, dw=None, dwH=None, k_AB=None, k_BA=None, inv_tcpmg=None, tcp=None, back_calc=None, num_points=None, power=None):
223 """The 2-site numerical solution to the Bloch-McConnell equation for SQ, ZQ, and DQ data. 224 225 The notation used here comes from: 226 227 - 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). 228 229 This function calculates and stores the R2eff values. 230 231 232 @keyword M0: This is a vector that contains the initial magnetizations corresponding to the A and B state transverse magnetizations. 233 @type M0: numpy float64, rank-1, 7D array 234 @keyword F_vector: The observable magnitisation vector. This defaults to [1, 0] for X observable magnitisation. 235 @type F_vector: numpy rank-1, 2D float64 array 236 @keyword m1: A complex numpy matrix to be populated. 237 @type m1: numpy rank-2, 2D complex64 array 238 @keyword m2: A complex numpy matrix to be populated. 239 @type m2: numpy rank-2, 2D complex64 array 240 @keyword R20A: The transverse, spin-spin relaxation rate for state A. 241 @type R20A: float 242 @keyword R20B: The transverse, spin-spin relaxation rate for state B. 243 @type R20B: float 244 @keyword pA: The population of state A. 245 @type pA: float 246 @keyword pB: The population of state B. 247 @type pB: float 248 @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. 249 @type dw: float 250 @keyword dwH: Unused - this is simply to match the r2eff_ns_mmq_2site_mq() function arguments. 251 @type dwH: float 252 @keyword k_AB: The rate of exchange from site A to B (rad/s). 253 @type k_AB: float 254 @keyword k_BA: The rate of exchange from site B to A (rad/s). 255 @type k_BA: float 256 @keyword inv_tcpmg: The inverse of the total duration of the CPMG element (in inverse seconds). 257 @type inv_tcpmg: float 258 @keyword tcp: The tau_CPMG times (1 / 4.nu1). 259 @type tcp: numpy rank-1 float array 260 @keyword back_calc: The array for holding the back calculated R2eff values. Each element corresponds to one of the CPMG nu1 frequencies. 261 @type back_calc: numpy rank-1 float array 262 @keyword num_points: The number of points on the dispersion curve, equal to the length of the tcp and back_calc arguments. 263 @type num_points: int 264 @keyword power: The matrix exponential power array. 265 @type power: numpy int16, rank-1 array 266 """ 267 268 # Populate the m1 and m2 matrices (only once per function call for speed). 269 populate_matrix(matrix=m1, R20A=R20A, R20B=R20B, dw=dw, k_AB=k_AB, k_BA=k_BA) 270 populate_matrix(matrix=m2, R20A=R20A, R20B=R20B, dw=-dw, k_AB=k_AB, k_BA=k_BA) 271 272 # Loop over the time points, back calculating the R2eff values. 273 for i in range(num_points): 274 # The A+/- matrices. 275 A_pos = matrix_exponential(m1*tcp[i]) 276 A_neg = matrix_exponential(m2*tcp[i]) 277 278 # The evolution for one n. 279 evol_block = dot(A_pos, dot(A_neg, dot(A_neg, A_pos))) 280 281 # The full evolution. 282 evol = square_matrix_power(evol_block, power[i]) 283 284 # The next lines calculate the R2eff using a two-point approximation, i.e. assuming that the decay is mono-exponential. 285 Mx = dot(F_vector, dot(evol, M0)) 286 Mx = Mx.real 287 if Mx <= 0.0 or isNaN(Mx): 288 back_calc[i] = 1e99 289 else: 290 back_calc[i] = -inv_tcpmg * log(Mx / pA)
291