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 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
52 from math import floor
53 from numpy import array, conj, dot, float64, log
54
55
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
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
136 populate_matrix(matrix=m1, R20A=R20A, R20B=R20B, dw=-dw-dwH, k_AB=k_AB, k_BA=k_BA)
137 populate_matrix(matrix=m2, R20A=R20A, R20B=R20B, dw=dw-dwH, k_AB=k_AB, k_BA=k_BA)
138
139
140 for i in range(num_points):
141
142 M1 = matrix_exponential(m1*tcp[i])
143 M2 = matrix_exponential(m2*tcp[i])
144
145
146 M1_star = conj(M1)
147 M2_star = conj(M2)
148
149
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
160 if power[i] == 1:
161
162 A = M1_M2
163
164
165 B = M1_M2_star
166
167
168 C = M2_M1
169
170
171 D = M2_M1_star
172
173
174 elif power[i] % 2 == 0:
175
176 fact = int(floor(power[i] / 2))
177
178
179 A = square_matrix_power(M1_M2_M2_M1, fact)
180
181
182 B = square_matrix_power(M2_M1_M1_M2_star, fact)
183
184
185 C = square_matrix_power(M2_M1_M1_M2, fact)
186
187
188 D = square_matrix_power(M1_M2_M2_M1_star, fact)
189
190
191 else:
192
193 fact = int(floor((power[i] - 1) / 2))
194
195
196 A = square_matrix_power(M1_M2_M2_M1, fact)
197 A = dot(A, M1_M2)
198
199
200 B = square_matrix_power(M1_M2_M2_M1_star, fact)
201 B = dot(B, M1_M2_star)
202
203
204 C = square_matrix_power(M2_M1_M1_M2, fact)
205 C = dot(C, M2_M1)
206
207
208 D = square_matrix_power(M2_M1_M1_M2_star, fact)
209 D = dot(D, M2_M1_star)
210
211
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
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
273 for i in range(num_points):
274
275 A_pos = matrix_exponential(m1*tcp[i])
276 A_neg = matrix_exponential(m2*tcp[i])
277
278
279 evol_block = dot(A_pos, dot(A_neg, dot(A_neg, A_pos)))
280
281
282 evol = square_matrix_power(evol_block, power[i])
283
284
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