1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23 """The CR72 model extended for MMQ CPMG data, called the U{MMQ CR72<http://wiki.nmr-relax.com/MMQ_CR72>} model.
24
25 Description
26 ===========
27
28 This module is for the function, gradient and Hessian of the U{MMQ CR72<http://wiki.nmr-relax.com/MMQ_CR72>} model.
29
30
31 References
32 ==========
33
34 The Carver and Richards (1972) 2-site model for all times scales was extended for multiple-MQ (MMQ) CPMG data by:
35
36 - Korzhnev, D. M., Kloiber, K., Kanelis, V., Tugarinov, V., and Kay, L. E. (2004). Probing slow dynamics in high molecular weight proteins by methyl-TROSY NMR spectroscopy: Application to a 723-residue enzyme. I{J. Am. Chem. Soc.}, B{126}, 3964-3973. (U{DOI: 10.1021/ja039587i<http://dx.doi.org/10.1021/ja039587i>}).
37
38
39 Links
40 =====
41
42 More information on the MMQ CR72 model can be found in the:
43
44 - U{relax wiki<http://wiki.nmr-relax.com/MMQ_CR72>},
45 - U{relax manual<http://www.nmr-relax.com/manual/MMQ_CR72_model.html>},
46 - U{relaxation dispersion page of the relax website<http://www.nmr-relax.com/analyses/relaxation_dispersion.html#MMQ_CR72>}.
47 """
48
49
50 from numpy import arccosh, array, cos, cosh, isfinite, log, max, sin, sqrt, sum
51
52
53 -def r2eff_mmq_cr72(r20=None, pA=None, pB=None, dw=None, dwH=None, kex=None, k_AB=None, k_BA=None, cpmg_frqs=None, inv_tcpmg=None, tcp=None, back_calc=None, num_points=None, power=None):
54 """The CR72 model extended to MMQ CPMG data.
55
56 This function calculates and stores the R2eff values.
57
58
59 @keyword r20: The R2 value in the absence of exchange.
60 @type r20: float
61 @keyword pA: The population of state A.
62 @type pA: float
63 @keyword pB: The population of state B.
64 @type pB: float
65 @keyword dw: The chemical exchange difference between states A and B in rad/s.
66 @type dw: float
67 @keyword dwH: The proton chemical exchange difference between states A and B in rad/s.
68 @type dwH: float
69 @keyword kex: The kex parameter value (the exchange rate in rad/s).
70 @type kex: float
71 @keyword k_AB: The rate of exchange from site A to B (rad/s).
72 @type k_AB: float
73 @keyword k_BA: The rate of exchange from site B to A (rad/s).
74 @type k_BA: float
75 @keyword cpmg_frqs: The CPMG nu1 frequencies.
76 @type cpmg_frqs: numpy rank-1 float array
77 @keyword inv_tcpmg: The inverse of the total duration of the CPMG element (in inverse seconds).
78 @type inv_tcpmg: float
79 @keyword tcp: The tau_CPMG times (1 / 4.nu1).
80 @type tcp: numpy rank-1 float array
81 @keyword back_calc: The array for holding the back calculated R2eff values. Each element corresponds to one of the CPMG nu1 frequencies.
82 @type back_calc: numpy rank-1 float array
83 @keyword num_points: The number of points on the dispersion curve, equal to the length of the tcp and back_calc arguments.
84 @type num_points: int
85 @keyword power: The matrix exponential power array.
86 @type power: numpy int16, rank-1 array
87 """
88
89
90 if (dw == 0.0 and dwH == 0.0) or pA == 1.0 or k_AB == 0.0:
91 back_calc[:] = array([r20]*num_points)
92 return
93
94
95 dw2 = dw**2
96 r20_kex = r20 + kex/2.0
97 pApBkex2 = k_AB * k_BA
98 isqrt_pApBkex2 = 1.j*sqrt(pApBkex2)
99 sqrt_pBpA = sqrt(pB/pA)
100 ikex = 1.j*kex
101
102
103 d = dwH + dw
104 dpos = d + ikex
105 dneg = d - ikex
106
107
108 z = dwH - dw
109 zpos = z + ikex
110 zneg = z - ikex
111
112
113 fact = 1.j*dwH + k_BA - k_AB
114 Psi = fact**2 - dw2 + 4.0*pApBkex2
115 zeta = -2.0*dw * fact
116
117
118 sqrt_psi2_zeta2 = sqrt(Psi**2 + zeta**2)
119
120
121 D_part = (Psi + 2.0*dw2) / sqrt_psi2_zeta2
122 Dpos = 0.5 * (1.0 + D_part)
123 Dneg = 0.5 * (-1.0 + D_part)
124
125
126 eta_scale = 2.0**(-3.0/2.0)
127 etapos_part = eta_scale * sqrt(Psi + sqrt_psi2_zeta2)
128 etaneg_part = eta_scale * sqrt(-Psi + sqrt_psi2_zeta2)
129
130
131 etapos = etapos_part / cpmg_frqs
132
133
134
135 if max(etapos) > 700:
136 back_calc[:] = array([r20]*num_points)
137 return
138
139
140 etaneg = etaneg_part / cpmg_frqs
141
142
143 mD = isqrt_pApBkex2 / (dpos * zpos) * (zpos + 2.0*dw*sin(zpos*tcp)/sin((dpos + zpos)*tcp))
144
145
146 mZ = -isqrt_pApBkex2 / (dneg * zneg) * (dneg - 2.0*dw*sin(dneg*tcp)/sin((dneg + zneg)*tcp))
147
148
149 Q = 1 - mD**2 + mD*mZ - mZ**2 + 0.5*(mD + mZ)*sqrt_pBpA
150 Q = Q.real
151
152
153 lambda1 = r20_kex - cpmg_frqs * arccosh(Dpos * cosh(etapos) - Dneg * cos(etaneg))
154
155
156 R2eff = lambda1.real - inv_tcpmg * log(Q)
157
158
159
160 if not isfinite(sum(R2eff)):
161 R2eff = array([1e100]*num_points)
162
163 back_calc[:] = R2eff
164