1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24 """The Baldwin (2014) 2-site exact solution model for all time scales U{B14<http://wiki.nmr-relax.com/B14>}.
25
26 Description
27 ===========
28
29 This module is for the function, gradient and Hessian of the U{B14<http://wiki.nmr-relax.com/B14>} model.
30
31
32 References
33 ==========
34
35 The model is named after the reference:
36
37 - 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>}).
38
39
40 Equations
41 =========
42
43 The equation used is::
44
45 R2A0 + R2B0 + kex Ncyc 1 ( 1+y 1-y )
46 R2eff = ------------------ - ------ * cosh^-1 * v1c - ------ ln( --- + ------------------ * (v2 + 2*kAB*pD ) ) ,
47 2 Trel Trel ( 2 2*sqrt(v1c^2 -1 ) )
48
49 1 ( 1+y 1-y )
50 = R2eff(CR72) - ------ ln( --- + ------------------ * (v2 + 2*kAB*pD ) ) ,
51 Trel ( 2 2*sqrt(v1c^2 -1 ) )
52
53 Which have these following definitions::
54
55 v1c = F0 * cosh(tauCP * E0)- F2 * cosh(tauCP * E2) ,
56 v1s = F0 * sinh(tauCP * E0)- F2 * sinh(tauCP * E2) ,
57 v2*N = v1s * (OB-OA) + 4OB * F1^a * sinh(tauCP * E1) ,
58 pD N = v1s + (F1^a + F1^b) * sinh(tauCP * E1) ,
59 v3 = ( v2^2 + 4 * kBA * kAB * pD^2 )^1/2 ,
60 y = ( (v1c-v3)/(v1c+v3) )^NCYC ,
61
62 Note, E2 is complex. If |x| denotes the complex modulus::
63
64 cosh(tauCP * E2) = cos(tauCP * |E2|) ,
65 sinh(tauCP * E2) = i sin(tauCP * |E2|) ,
66
67 The term pD is based on product of the off diagonal elements in the CPMG propagator (Supplementary Section 3).
68
69 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::
70
71 sqrt(v1c^2-1) ~ v2 + 2*kAB * pD ,
72
73 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.
74
75 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.
76
77 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.
78
79
80 Links
81 =====
82
83 More information on the B14 model can be found in the:
84
85 - U{relax wiki<http://wiki.nmr-relax.com/B14>},
86 - U{relax manual<http://www.nmr-relax.com/manual/reduced_B14_2_site_CPMG_model.html>},
87 - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#B14>}.
88
89 More information on the B14 full model can be found in the:
90
91 - U{relax wiki<http://wiki.nmr-relax.com/B14_full>},
92 - U{relax manual<http://www.nmr-relax.com/manual/full_B14_2_site_CPMG_model.html>},
93 - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#B14_full>}.
94
95
96 Comparison to CR72 model
97 ========================
98
99 Comparison to CR72 model can be found in the:
100
101 - U{relax wiki<http://wiki.nmr-relax.com/CR72>},
102 - U{relax manual<http://www.nmr-relax.com/manual/reduced_CR72_2_site_CPMG_model.html>},
103 - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#CR72>}.
104
105 Comparison to CR72 full model can be found in the:
106
107 - U{relax wiki<http://wiki.nmr-relax.com/CR72_full>},
108 - U{relax manual<http://www.nmr-relax.com/manual/full_CR72_2_site_CPMG_model.html>},
109 - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#CR72_full>}.
110 """
111
112
113 from numpy import arccosh, arctan2, array, cos, cosh, isfinite, log, max, power, sin, sinh, sqrt, sum
114
115
116 g_fact = 1/sqrt(2)
117
118 -def r2eff_B14(r20a=None, r20b=None, pA=None, pB=None, dw=None, kex=None, k_AB=None, k_BA=None, ncyc=None, inv_tcpmg=None, tcp=None, back_calc=None, num_points=None):
119 """Calculate the R2eff values for the CR72 model.
120
121 See the module docstring for details.
122
123
124 @keyword r20a: The R20 parameter value of state A (R2 with no exchange).
125 @type r20a: float
126 @keyword r20b: The R20 parameter value of state B (R2 with no exchange).
127 @type r20b: float
128 @keyword pA: The population of state A.
129 @type pA: float
130 @keyword pB: The population of state B.
131 @type pB: float
132 @keyword dw: The chemical exchange difference between states A and B in rad/s.
133 @type dw: float
134 @keyword kex: The kex parameter value (the exchange rate in rad/s).
135 @type kex: float
136 @keyword k_AB: The rate of exchange from site A to B (rad/s).
137 @type k_AB: float
138 @keyword k_BA: The rate of exchange from site B to A (rad/s).
139 @type k_BA: float
140 @keyword ncyc: The matrix exponential power array. The number of CPMG blocks.
141 @type ncyc: numpy int16, rank-1 array
142 @keyword inv_tcpmg: The inverse of the total duration of the CPMG element (in inverse seconds).
143 @type inv_tcpmg: float
144 @keyword tcp: The tau_CPMG times (1 / 4.nu1).
145 @type tcp: numpy rank-1 float array
146 @keyword back_calc: The array for holding the back calculated R2eff values. Each element corresponds to one of the CPMG nu1 frequencies.
147 @type back_calc: numpy rank-1 float array
148 @keyword num_points: The number of points on the dispersion curve, equal to the length of the cpmg_frqs and back_calc arguments.
149 @type num_points: int
150 """
151
152
153 if dw == 0.0 or pA == 1.0 or k_AB == 0.0:
154 back_calc[:] = array([r20a]*num_points)
155 return
156
157
158 deltaR2 = r20a - r20b
159
160
161 alpha_m = deltaR2 + k_AB - k_BA
162 zeta = 2.0 * dw * alpha_m
163 Psi = alpha_m**2 + 4.0 * k_BA * k_AB - dw**2
164
165
166 dw2 = dw**2
167 two_tcp = 2.0 * tcp
168
169
170
171 quad_zeta2_Psi2 = (zeta**2 + Psi**2)**0.25
172 g3 = cos(0.5 * arctan2(-zeta, Psi)) * quad_zeta2_Psi2
173 g4 = sin(0.5 * arctan2(-zeta, Psi)) * quad_zeta2_Psi2
174
175
176 g32 = g3**2
177 g42 = g4**2
178
179
180
181 N = g3 + g4*1j
182
183 NNc = g32 + g42
184
185
186 F0 = (dw2 + g32) / NNc
187
188
189 F2 = (dw2 - g42) / NNc
190
191
192
193
194 F1b = (dw + g4) * (dw - g3*1j) / NNc
195
196
197 F1a_plus_b = (2. * dw2 + zeta*1j) / NNc
198
199
200
201 E0 = two_tcp * g3
202
203
204
205 if max(E0) > 700:
206 back_calc[:] = array([r20a]*num_points)
207 return
208
209
210 E2 = two_tcp * g4
211
212
213 E1 = (g3 - g4*1j) * tcp
214
215
216 v1s = F0 * sinh(E0) - F2 * sin(E2)*1j
217
218
219 v4 = F1b * (-alpha_m - g3 ) + F1b * (dw - g4)*1j
220
221
222 ex1c = sinh(E1)
223
224
225 v5 = (-deltaR2 + kex + dw*1j) * v1s - 2. * (v4 + k_AB * F1a_plus_b) * ex1c
226
227
228 v1c = F0 * cosh(E0) - F2 * cos(E2)
229
230
231 v3 = sqrt(v1c**2 - 1.)
232
233 y = power( (v1c - v3) / (v1c + v3), ncyc)
234
235 Tog = 0.5 * (1. + y) + (1. - y) * v5 / (2. * v3 * N )
236
237
238
239
240
241
242
243
244
245
246
247
248 R2eff = (r20a + r20b + kex) / 2.0 - inv_tcpmg * ( ncyc * arccosh(v1c.real) + log(Tog.real) )
249
250
251
252 if not isfinite(sum(R2eff)):
253 R2eff = array([1e100]*num_points)
254
255 back_calc[:] = R2eff
256