1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
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
58 from math import floor
59 from numpy import array, conj, dot, float64, log
60
61
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
97 matrix[0, 0] = -k_AB - k_AC - R20A
98 matrix[0, 1] = k_BA
99 matrix[0, 2] = k_CA
100
101
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
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
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)
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)
181
182
183 for i in range(num_points):
184
185 M1 = matrix_exponential(m1*tcp[i])
186 M2 = matrix_exponential(m2*tcp[i])
187
188
189 M1_star = conj(M1)
190 M2_star = conj(M2)
191
192
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
203 if power[i] == 1:
204
205 A = M1_M2
206
207
208 B = M1_M2_star
209
210
211 C = M2_M1
212
213
214 D = M2_M1_star
215
216
217 elif power[i] % 2 == 0:
218
219 fact = int(floor(power[i] / 2))
220
221
222 A = square_matrix_power(M1_M2_M2_M1, fact)
223
224
225 B = square_matrix_power(M2_M1_M1_M2_star, fact)
226
227
228 C = square_matrix_power(M2_M1_M1_M2, fact)
229
230
231 D = square_matrix_power(M1_M2_M2_M1_star, fact)
232
233
234 else:
235
236 fact = int(floor((power[i] - 1) / 2))
237
238
239 A = square_matrix_power(M1_M2_M2_M1, fact)
240 A = dot(A, M1_M2)
241
242
243 B = square_matrix_power(M1_M2_M2_M1_star, fact)
244 B = dot(B, M1_M2_star)
245
246
247 C = square_matrix_power(M2_M1_M1_M2, fact)
248 C = dot(C, M2_M1)
249
250
251 D = square_matrix_power(M2_M1_M1_M2_star, fact)
252 D = dot(D, M2_M1_star)
253
254
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
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
324 for i in range(num_points):
325
326 A_pos = matrix_exponential(m1*tcp[i])
327 A_neg = matrix_exponential(m2*tcp[i])
328
329
330 evol_block = dot(A_pos, dot(A_neg, dot(A_neg, A_pos)))
331
332
333 evol = square_matrix_power(evol_block, power[i])
334
335
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