Author: bugman
Date: Mon Nov 18 13:01:01 2013
New Revision: 21498
URL: http://svn.gna.org/viewcvs/relax?rev=21498&view=rev
Log:
Implemented the 'MMQ 2-site' CPMG model equations from the Korzhnev et al.,
2005 reference.
The paper reference is:
- 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).
The original code from Mathilde Lescanne and Dominique Marion has only
slightly been modified for
this change as the MQ data treatment in the Korzhnev et al., 2004 reference
is the same as in the
2005 reference, but using a different notation. This has been renamed to
r2eff_mmq_2site_mq().
The new r2eff_mmq_2site_sq_dq_zq() function has been added to the
lib.dispersion.mmq_2site module to
allows the SQ, DQ, and ZQ R2eff data to be calculated. This function follows
the notation of the
2005 paper.
The populate_matrix() function has been modified to only accept one combined
chemical shift
difference value. It can now also accept different values for R20A and R20B,
though the mmq_2site
module defaults to R20A=R20B.
Modified:
branches/relax_disp/lib/dispersion/mmq_2site.py
branches/relax_disp/target_functions/relax_disp.py
Modified: branches/relax_disp/lib/dispersion/mmq_2site.py
URL:
http://svn.gna.org/viewcvs/relax/branches/relax_disp/lib/dispersion/mmq_2site.py?rev=21498&r1=21497&r2=21498&view=diff
==============================================================================
--- branches/relax_disp/lib/dispersion/mmq_2site.py (original)
+++ branches/relax_disp/lib/dispersion/mmq_2site.py Mon Nov 18 13:01:01 2013
@@ -29,35 +29,26 @@
- 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).
"""
-# Dependency check module.
-import dep_check
-
# Python module imports.
-from math import fabs, log
-from numpy import array, conj, dot, float64
+from numpy import array, dot, float64, log
# relax module imports.
-from lib.dispersion.ns_matrices import rcpmg_3d
from lib.float import isNaN
from lib.linear_algebra.matrix_exponential import matrix_exponential
from lib.linear_algebra.matrix_power import square_matrix_power
-def populate_matrix(matrix=None, r20=None, dw=None, dwH=None, k_AB=None,
k_BA=None):
- """Matrix for HMQC experiments.
-
- This corresponds to matrix m1 and m2 matrices of equation 2.2 from:
-
- - 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. J. Am.
Chem. Soc., 126, 3964-3973. (U{DOI:
10.1021/ja039587i<http://dx.doi.org/10.1021/ja039587i>}).
+def populate_matrix(matrix=None, R20A=None, R20B=None, dw=None, k_AB=None,
k_BA=None):
+ """The Bloch-McConnell matrix for 2-site exchange.
@keyword matrix: The matrix to populate.
@type matrix: numpy rank-2, 2D complex64 array
- @keyword r20: The R2 value in the absence of exchange.
- @type r20: float
- @keyword dw: The chemical exchange difference between states
A and B in rad/s.
+ @keyword R20A: The transverse, spin-spin relaxation rate for
state A.
+ @type R20A: float
+ @keyword R20B: The transverse, spin-spin relaxation rate for
state B.
+ @type R20B: float
+ @keyword dw: The combined chemical exchange difference
parameters between states A and B in rad/s. This can be any combination of
dwX and dwH.
@type dw: float
- @keyword dwH: The proton chemical exchange difference between
states A and B in rad/s.
- @type dwH: float
@keyword k_AB: The rate of exchange from site A to B (rad/s).
@type k_AB: float
@keyword k_BA: The rate of exchange from site B to A (rad/s).
@@ -65,14 +56,22 @@
"""
# Fill in the elements.
- matrix[0, 0] = -k_AB - r20
+ matrix[0, 0] = -k_AB - R20A
matrix[0, 1] = k_BA
matrix[1, 0] = k_AB
- matrix[1, 1] = -k_BA - 1.j*(dwH + dw) - r20
-
-
-def r2eff_mmq_2site(M0=None, F_vector=array([1, 0], float64), m1=None,
m2=None, r20=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, n=None):
- """The 2-site numerical solution to the Bloch-McConnell equation for SQ,
ZQ, DQ and MQ data.
+ matrix[1, 1] = -k_BA + 1.j*dw - R20B
+
+
+def r2eff_mmq_2site_mq(M0=None, m1=None, m2=None, r20=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, n=None):
+ """The 2-site numerical solution to the Bloch-McConnell equation for MQ
data.
+
+ The notation used here comes from:
+
+ - 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).
+
+ and:
+
+ - 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).
This function calculates and stores the R2eff values.
@@ -112,18 +111,18 @@
"""
# Populate the m1 and m2 matrices (only once per function call for
speed).
- populate_matrix(matrix=m1, r20=r20, dw=dw, dwH=dwH, k_AB=k_AB, k_BA=k_BA)
- populate_matrix(matrix=m2, r20=r20, dw=-dw, dwH=dwH, k_AB=k_AB,
k_BA=k_BA)
+ populate_matrix(matrix=m1, R20A=r20, R20B=r20, dw=dw+dwH, k_AB=k_AB,
k_BA=k_BA) # D+ matrix component.
+ populate_matrix(matrix=m2, R20A=r20, R20B=r20, dw=-dw+dwH, k_AB=k_AB,
k_BA=k_BA) # Z- matrix component.
# Loop over the time points, back calculating the R2eff values.
for i in range(num_points):
# The M1 and M2 matrices.
- M1 = matrix_exponential(m1*tcp[i])
- M2 = matrix_exponential(m2*tcp[i])
+ M1 = matrix_exponential(m1*tcp[i]) # Equivalent to D+.
+ M2 = matrix_exponential(m2*tcp[i]) # Equivalent to Z-.
# The complex conjugates M1* and M2*
- M1_star = conj(M1)
- M2_star = conj(M2)
+ M1_star = conj(M1) # Equivalent to D+*.
+ M2_star = conj(M2) # Equivalent to Z-*.
# Repetitive dot products (minimised for speed).
M1_M2 = dot(M1, M2)
@@ -176,9 +175,76 @@
# The next lines calculate the R2eff using a two-point
approximation, i.e. assuming that the decay is mono-exponential.
A_B = dot(A, B)
C_D = dot(C, D)
- Mx = dot(dot(F_vector, (A_B + C_D)), M0)
- Mx = Mx.real / (2.0 * pA)
+ Mx = (A_B[0, 0] + C_D[0, 0])*M0[0] + (A_B[0, 1] + C_D[0, 1])*M0[1]
+ Mx = Mx.real
if Mx <= 0.0 or isNaN(Mx):
back_calc[i] = 1e99
else:
- back_calc[i]= -inv_tcpmg * log(Mx)
+ back_calc[i]= -inv_tcpmg * log(Mx / pA)
+
+
+def r2eff_mmq_2site_sq_dq_zq(M0=None, F_vector=array([1, 0], float64),
m1=None, m2=None, r20=None, pA=None, pB=None, dwX=None, k_AB=None, k_BA=None,
inv_tcpmg=None, tcp=None, back_calc=None, num_points=None, power=None):
+ """The 2-site numerical solution to the Bloch-McConnell equation for SQ,
ZQ, and DQ data.
+
+ The notation used here comes from:
+
+ - 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).
+
+ This function calculates and stores the R2eff values.
+
+
+ @keyword M0: This is a vector that contains the initial
magnetizations corresponding to the A and B state transverse magnetizations.
+ @type M0: numpy float64, rank-1, 7D array
+ @keyword F_vector: The observable magnitisation vector. This
defaults to [1, 0] for X observable magnitisation.
+ @type F_vector: numpy rank-1, 2D float64 array
+ @keyword m1: A complex numpy matrix to be populated.
+ @type m1: numpy rank-2, 2D complex64 array
+ @keyword m2: A complex numpy matrix to be populated.
+ @type m2: numpy rank-2, 2D complex64 array
+ @keyword r20: The R2 value in the absence of exchange.
+ @type r20: float
+ @keyword pA: The population of state A.
+ @type pA: float
+ @keyword pB: The population of state B.
+ @type pB: float
+ @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.
+ @type dw: float
+ @keyword k_AB: The rate of exchange from site A to B (rad/s).
+ @type k_AB: float
+ @keyword k_BA: The rate of exchange from site B to A (rad/s).
+ @type k_BA: float
+ @keyword inv_tcpmg: The inverse of the total duration of the CPMG
element (in inverse seconds).
+ @type inv_tcpmg: float
+ @keyword tcp: The tau_CPMG times (1 / 4.nu1).
+ @type tcp: numpy rank-1 float array
+ @keyword back_calc: The array for holding the back calculated R2eff
values. Each element corresponds to one of the CPMG nu1 frequencies.
+ @type back_calc: numpy rank-1 float array
+ @keyword num_points: The number of points on the dispersion curve,
equal to the length of the tcp and back_calc arguments.
+ @type num_points: int
+ @keyword power: The matrix exponential power array.
+ @type power: numpy int16, rank-1 array
+ """
+
+ # Populate the m1 and m2 matrices (only once per function call for
speed).
+ populate_matrix(matrix=m1, R20A=r20, R20B=r20, dw=dwX, k_AB=k_AB,
k_BA=k_BA)
+ populate_matrix(matrix=m2, R20A=r20, R20B=r20, dw=-dwX, k_AB=k_AB,
k_BA=k_BA)
+
+ # Loop over the time points, back calculating the R2eff values.
+ for i in range(num_points):
+ # The A+/- matrices.
+ A_pos = matrix_exponential(m1*tcp[i])
+ A_neg = matrix_exponential(m2*tcp[i])
+
+ # The evolution for one n.
+ evol_block = dot(A_pos, dot(A_neg, dot(A_neg, A_pos)))
+
+ # The full evolution.
+ evol = square_matrix_power(evol_block, power[i])
+
+ # The next lines calculate the R2eff using a two-point
approximation, i.e. assuming that the decay is mono-exponential.
+ Mx = dot(F_vector, dot(evol, M0))
+ Mx = Mx.real
+ if Mx <= 0.0 or isNaN(Mx):
+ back_calc[i] = 1e99
+ else:
+ back_calc[i] = -inv_tcpmg * log(Mx / pA)
Modified: branches/relax_disp/target_functions/relax_disp.py
URL:
http://svn.gna.org/viewcvs/relax/branches/relax_disp/target_functions/relax_disp.py?rev=21498&r1=21497&r2=21498&view=diff
==============================================================================
--- branches/relax_disp/target_functions/relax_disp.py (original)
+++ branches/relax_disp/target_functions/relax_disp.py Mon Nov 18 13:01:01
2013
@@ -38,7 +38,7 @@
from lib.dispersion.m61b import r1rho_M61b
from lib.dispersion.mp05 import r1rho_MP05
from lib.dispersion.mq_cr72 import r2eff_mq_cr72
-from lib.dispersion.mmq_2site import r2eff_mmq_2site
+from lib.dispersion.mmq_2site import r2eff_mmq_2site_mq,
r2eff_mmq_2site_sq_dq_zq
from lib.dispersion.ns_cpmg_2site_3d import r2eff_ns_cpmg_2site_3D
from lib.dispersion.ns_cpmg_2site_expanded import
r2eff_ns_cpmg_2site_expanded
from lib.dispersion.ns_cpmg_2site_star import r2eff_ns_cpmg_2site_star
@@ -264,10 +264,12 @@
if model in [MODEL_MMQ_2SITE, MODEL_NS_CPMG_2SITE_3D,
MODEL_NS_CPMG_2SITE_3D_FULL, MODEL_NS_CPMG_2SITE_EXPANDED,
MODEL_NS_CPMG_2SITE_STAR, MODEL_NS_CPMG_2SITE_STAR_FULL,
MODEL_NS_R1RHO_2SITE]:
self.inv_relax_time = 1.0 / relax_time
- # Special matrices for the multi-quantum CPMG 2-site numerical model.
+ # Special storage matrices for the multi-quantum CPMG 2-site
numerical model.
if model == MODEL_MMQ_2SITE:
self.m1 = zeros((2, 2), complex64)
self.m2 = zeros((2, 2), complex64)
+ self.m3 = zeros((2, 2), complex64)
+ self.m4 = zeros((2, 2), complex64)
# Set up the model.
if model == MODEL_NOREX:
@@ -938,8 +940,19 @@
dw_frq = dw[spin_index] *
self.frqs[exp_index][spin_index][frq_index]
dwH_frq = dwH[spin_index] *
self.frqs[exp_index][spin_index][frq_index]
- # Back calculate the R2eff values.
- r2eff_mmq_2site(M0=self.M0, m1=self.m1, m2=self.m2,
r20=R20[r20_index], pA=pA, pB=pB, dw=dw_frq, dwH=dwH_frq, k_AB=k_AB,
k_BA=k_BA, inv_tcpmg=self.inv_relax_time,
tcp=self.tau_cpmg[exp_index][frq_index],
back_calc=self.back_calc[exp_index][spin_index][frq_index],
num_points=self.num_disp_points[exp_index][frq_index],
n=self.n[exp_index][frq_index])
+ # Back calculate the R2eff values for each experiment
type.
+ if self.exp_types[exp_index] == EXP_TYPE_CPMG:
+ r2eff_mmq_2site_sq_dq_zq(M0=self.M0, m1=self.m1,
m2=self.m2, r20=R20[r20_index], pA=pA, pB=pB, dwX=dw_frq, k_AB=k_AB,
k_BA=k_BA, inv_tcpmg=self.inv_relax_time,
tcp=self.tau_cpmg[exp_index][frq_index],
back_calc=self.back_calc[exp_index][spin_index][frq_index],
num_points=self.num_disp_points[exp_index][frq_index],
power=self.power[exp_index][frq_index])
+ elif self.exp_types[exp_index] ==
EXP_TYPE_PROTON_SQ_CPMG:
+ r2eff_mmq_2site_sq_dq_zq(M0=self.M0, m1=self.m1,
m2=self.m2, r20=R20[r20_index], pA=pA, pB=pB, dwX=dwH_frq, k_AB=k_AB,
k_BA=k_BA, inv_tcpmg=self.inv_relax_time,
tcp=self.tau_cpmg[exp_index][frq_index],
back_calc=self.back_calc[exp_index][spin_index][frq_index],
num_points=self.num_disp_points[exp_index][frq_index],
power=self.power[exp_index][frq_index])
+ elif self.exp_types[exp_index] == EXP_TYPE_DQ_CPMG:
+ r2eff_mmq_2site_sq_dq_zq(M0=self.M0, m1=self.m1,
m2=self.m2, r20=R20[r20_index], pA=pA, pB=pB, dwX=dwH_frq+dw_frq, k_AB=k_AB,
k_BA=k_BA, inv_tcpmg=self.inv_relax_time,
tcp=self.tau_cpmg[exp_index][frq_index],
back_calc=self.back_calc[exp_index][spin_index][frq_index],
num_points=self.num_disp_points[exp_index][frq_index],
power=self.power[exp_index][frq_index])
+ elif self.exp_types[exp_index] == EXP_TYPE_ZQ_CPMG:
+ r2eff_mmq_2site_sq_dq_zq(M0=self.M0, m1=self.m1,
m2=self.m2, r20=R20[r20_index], pA=pA, pB=pB, dwX=dwH_frq-dw_frq, k_AB=k_AB,
k_BA=k_BA, inv_tcpmg=self.inv_relax_time,
tcp=self.tau_cpmg[exp_index][frq_index],
back_calc=self.back_calc[exp_index][spin_index][frq_index],
num_points=self.num_disp_points[exp_index][frq_index],
power=self.power[exp_index][frq_index])
+ elif self.exp_types[exp_index] == EXP_TYPE_MQ_CPMG:
+ r2eff_mmq_2site_mq(M0=self.M0, m1=self.m1,
m2=self.m2, m3=self.m3, m4=self.m4, r20=R20[r20_index], pA=pA, pB=pB,
dw=dw_frq, dwH=dwH_frq, k_AB=k_AB, k_BA=k_BA, inv_tcpmg=self.inv_relax_time,
tcp=self.tau_cpmg[exp_index][frq_index],
back_calc=self.back_calc[exp_index][spin_index][frq_index],
num_points=self.num_disp_points[exp_index][frq_index],
power=self.power[exp_index][frq_index])
+ elif self.exp_types[exp_index] ==
EXP_TYPE_PROTON_MQ_CPMG:
+ r2eff_mmq_2site_mq(M0=self.M0, m1=self.m1,
m2=self.m2, m3=self.m3, m4=self.m4, r20=R20[r20_index], pA=pA, pB=pB,
dw=dwH_frq, dwH=dw_frq, k_AB=k_AB, k_BA=k_BA, inv_tcpmg=self.inv_relax_time,
tcp=self.tau_cpmg[exp_index][frq_index],
back_calc=self.back_calc[exp_index][spin_index][frq_index],
num_points=self.num_disp_points[exp_index][frq_index],
power=self.power[exp_index][frq_index])
# For all missing data points, set the back-calculated
value to the measured values so that it has no effect on the chi-squared
value.
for point_index in
range(self.num_disp_points[exp_index][frq_index]):
r21498 - in /branches/relax_disp: lib/dispersion/mmq_2site.py target_functions/relax_disp.py