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 Baldwin (2014) 2-site exact solution model for all time scales U{B14<http://wiki.nmr-relax.com/B14>}.
26
27 Description
28 ===========
29
30 This module is for the function, gradient and Hessian of the U{B14<http://wiki.nmr-relax.com/B14>} model.
31
32
33 References
34 ==========
35
36 The model is named after the reference:
37
38 - Andrew J. Baldwin (2014). An exact solution for R2,eff in CPMG experiments in the case of two site chemical exchange. I{J. Magn. Reson.}. (U{DOI: 10.1016/j.jmr.2014.02.023 <http://dx.doi.org/10.1016/j.jmr.2014.02.023>}).
39
40
41 Equations
42 =========
43
44 The equation used is::
45
46 R2A0 + R2B0 + kex Ncyc 1 ( 1+y 1-y )
47 R2eff = ------------------ - ------ * cosh^-1 * v1c - ------ ln( --- + ------------------ * (v2 + 2*kAB*pD ) ) ,
48 2 Trel Trel ( 2 2*sqrt(v1c^2 -1 ) )
49
50 1 ( 1+y 1-y )
51 = R2eff(CR72) - ------ ln( --- + ------------------ * (v2 + 2*kAB*pD ) ) ,
52 Trel ( 2 2*sqrt(v1c^2 -1 ) )
53
54 Which have these following definitions::
55
56 v1c = F0 * cosh(tauCP * E0)- F2 * cosh(tauCP * E2) ,
57 v1s = F0 * sinh(tauCP * E0)- F2 * sinh(tauCP * E2) ,
58 v2*N = v1s * (OB-OA) + 4OB * F1^a * sinh(tauCP * E1) ,
59 pD N = v1s + (F1^a + F1^b) * sinh(tauCP * E1) ,
60 v3 = ( v2^2 + 4 * kBA * kAB * pD^2 )^1/2 ,
61 y = ( (v1c-v3)/(v1c+v3) )^NCYC ,
62
63 Note, E2 is complex. If |x| denotes the complex modulus::
64
65 cosh(tauCP * E2) = cos(tauCP * |E2|) ,
66 sinh(tauCP * E2) = i sin(tauCP * |E2|) ,
67
68 The term pD is based on product of the off diagonal elements in the CPMG propagator (Supplementary Section 3).
69
70 It is interesting to consider the region of validity of the Carver Richards result. The two results are equal when the correction is zero, which is true when::
71
72 sqrt(v1c^2-1) ~ v2 + 2*kAB * pD ,
73
74 This occurs when 2*kAB * pD tends to zero, and so v2=v3. Setting "kAB * pD" to zero, amounts to neglecting magnetisation that starts on the ground state ensemble and end on the excited state ensemble and vice versa. This will be a good approximation when pA >> p_B.
75
76 In practise, significant deviations from the Carver Richards equation can be incurred if pB > 1 %. Incorporation of the correction term into equation (50), results in an improved description of the CPMG experiment over the Carver Richards equation.
77
78 kex is the chemical exchange rate constant, pA and pB are the populations of states A and B, and delta_omega is the chemical shift difference between the two states in ppm.
79
80
81 Links
82 =====
83
84 More information on the B14 model can be found in the:
85
86 - U{relax wiki<http://wiki.nmr-relax.com/B14>},
87 - U{relax manual<http://www.nmr-relax.com/manual/The_reduced_B14_2_site_CPMG_model.html>},
88 - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#B14>}.
89
90 More information on the B14 full model can be found in the:
91
92 - U{relax wiki<http://wiki.nmr-relax.com/B14_full>},
93 - U{relax manual<http://www.nmr-relax.com/manual/The_full_B14_2_site_CPMG_model.html>},
94 - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#B14_full>}.
95
96
97 Comparison to CR72 model
98 ========================
99
100 Comparison to CR72 model can be found in the:
101
102 - U{relax wiki<http://wiki.nmr-relax.com/CR72>},
103 - U{relax manual<http://www.nmr-relax.com/manual/The_reduced_CR72_2_site_CPMG_model.html>},
104 - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#CR72>}.
105
106 Comparison to CR72 full model can be found in the:
107
108 - U{relax wiki<http://wiki.nmr-relax.com/CR72_full>},
109 - U{relax manual<http://www.nmr-relax.com/manual/The_full_CR72_2_site_CPMG_model.html>},
110 - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#CR72_full>}.
111 """
112
113
114 from numpy import any, arccosh, arctan2, cos, cosh, fabs, isfinite, log, max, min, power, sin, sinh, sqrt, sum
115 from numpy.ma import fix_invalid, masked_greater_equal, masked_where
116
117
118 g_fact = 1/sqrt(2)
119
120
121 -def r2eff_B14(r20a=None, r20b=None, pA=None, dw=None, dw_orig=None, kex=None, ncyc=None, inv_tcpmg=None, tcp=None, back_calc=None):
122 """Calculate the R2eff values for the CR72 model.
123
124 See the module docstring for details.
125
126
127 @keyword r20a: The R20 parameter value of state A (R2 with no exchange).
128 @type r20a: numpy float array of rank [NE][NS][NM][NO][ND]
129 @keyword r20b: The R20 parameter value of state B (R2 with no exchange).
130 @type r20b: numpy float array of rank [NE][NS][NM][NO][ND]
131 @keyword pA: The population of state A.
132 @type pA: float
133 @keyword dw: The chemical exchange difference between states A and B in rad/s.
134 @type dw: numpy float array of rank [NE][NS][NM][NO][ND]
135 @keyword dw_orig: The chemical exchange difference between states A and B in ppm. This is only for faster checking of zero value, which result in no exchange.
136 @type dw_orig: numpy float array of rank-1
137 @keyword kex: The kex parameter value (the exchange rate in rad/s).
138 @type kex: float
139 @keyword ncyc: The matrix exponential power array. The number of CPMG blocks.
140 @type ncyc: numpy int16 array of rank [NE][NS][NM][NO][ND]
141 @keyword inv_tcpmg: The inverse of the total duration of the CPMG element (in inverse seconds).
142 @type inv_tcpmg: numpy float array of rank [NE][NS][NM][NO][ND]
143 @keyword tcp: The tau_CPMG times (1 / 4.nu1).
144 @type tcp: numpy float array of rank [NE][NS][NM][NO][ND]
145 @keyword back_calc: The array for holding the back calculated R2eff values. Each element corresponds to one of the CPMG nu1 frequencies.
146 @type back_calc: numpy float array of rank [NE][NS][NM][NO][ND]
147 """
148
149
150 t_dw_zero = False
151 t_max_e = False
152 t_v3_N_zero = False
153 t_log_tog_neg = False
154 t_v1c_less_one = False
155
156
157
158 if kex == 0.0 or pA == 1.0:
159 back_calc[:] = r20a
160 return
161
162
163 if min(fabs(dw_orig)) == 0.0:
164 t_dw_zero = True
165 mask_dw_zero = masked_where(dw == 0.0, dw)
166
167
168 pB = 1.0 - pA
169 k_BA = pA * kex
170 k_AB = pB * kex
171
172
173 deltaR2 = r20a - r20b
174 dw2 = dw**2
175 two_tcp = 2.0 * tcp
176
177
178 alpha_m = deltaR2 + k_AB - k_BA
179 zeta = 2.0 * dw * alpha_m
180 Psi = alpha_m**2 + 4.0 * k_BA * k_AB - dw2
181
182
183
184 quad_zeta2_Psi2 = (zeta**2 + Psi**2)**0.25
185 fact = 0.5 * arctan2(-zeta, Psi)
186 g3 = cos(fact) * quad_zeta2_Psi2
187 g4 = sin(fact) * quad_zeta2_Psi2
188
189
190 g32 = g3**2
191 g42 = g4**2
192
193
194
195 N = g3 + g4*1j
196
197 NNc = g32 + g42
198
199
200 F0 = (dw2 + g32) / NNc
201
202
203 F2 = (dw2 - g42) / NNc
204
205
206
207
208 F1b = (dw + g4) * (dw - g3*1j) / NNc
209
210
211 F1a_plus_b = (2. * dw2 + zeta*1j) / NNc
212
213
214
215 E0 = two_tcp * g3
216
217
218
219 if max(E0) > 700:
220 t_max_e = True
221 mask_max_e = masked_greater_equal(E0, 700.0)
222
223 E0[mask_max_e.mask] = 1.0
224
225
226 E2 = two_tcp * g4
227
228
229 E1 = (g3 - g4*1j) * tcp
230
231
232 v1s = F0 * sinh(E0) - F2 * sin(E2)*1j
233
234
235 v4 = F1b * (-alpha_m - g3 ) + F1b * (dw - g4)*1j
236
237
238 ex1c = sinh(E1)
239
240
241 v5 = (-deltaR2 + kex + dw*1j) * v1s - 2. * (v4 + k_AB * F1a_plus_b) * ex1c
242
243
244 v1c = F0 * cosh(E0) - F2 * cos(E2)
245
246
247
248 mask_v1c_less_one = v1c < 1.0
249 if any(mask_v1c_less_one):
250 t_v1c_less_one = True
251 v1c[mask_v1c_less_one] = 1.0
252
253
254 v3 = sqrt(v1c**2 - 1.)
255
256 y = power( (v1c - v3) / (v1c + v3), ncyc)
257
258 Tog_div = 2. * v3 * N
259
260
261
262 mask_v3_N_zero = Tog_div == 0.0
263 if any(mask_v3_N_zero):
264 t_v3_N_zero = True
265 Tog_div[mask_v3_N_zero] = 1.0
266
267 Tog = 0.5 * (1. + y) + (1. - y) * v5 / Tog_div
268
269
270
271
272
273
274
275
276
277
278
279
280
281 mask_log_tog_neg = Tog.real < 0.0
282 if any(mask_log_tog_neg):
283 t_log_tog_neg = True
284 Tog.real[mask_log_tog_neg] = 1.0
285
286
287 back_calc[:] = (r20a + r20b + kex) / 2.0 - inv_tcpmg * ( ncyc * arccosh(v1c.real) + log(Tog.real) )
288
289
290
291 if t_dw_zero:
292 back_calc[mask_dw_zero.mask] = r20a[mask_dw_zero.mask]
293
294
295 if t_max_e:
296 back_calc[mask_max_e.mask] = r20a[mask_max_e.mask]
297
298
299 if t_v1c_less_one:
300 back_calc[mask_v1c_less_one] = 1e100
301
302
303 if t_v3_N_zero:
304 back_calc[mask_v3_N_zero] = 1e100
305
306
307 if t_log_tog_neg:
308 back_calc[mask_log_tog_neg] = 1e100
309
310
311
312 if not isfinite(sum(back_calc)):
313
314 fix_invalid(back_calc, copy=False, fill_value=1e100)
315