That CarverRichards number now looks what I would personally expect
from the theory!
I'd agree with that. This is noteworthy - we're agreeing on something! Had to
happen eventually :)
I don't like the Nikolai results, but I also don't
like the fact that I have two completely different implementations of
this code (from the sim_all software and the compareTroells.py
script).
The different versions of code is a good point. I copied and pasted what I
think is your function here:
http://www.nmr-relax.com/api/3.1/lib.dispersion.ns_cpmg_2site_expanded-pysrc.html
And have added this as the fun mode (see function nikolaiFun). My Nikolai
version came from 'reiko_ishima_v1.tar' and 'Exact_Rex_mod_v1.m' as i'm sure
he'll explain shortly. I ran my python version against the matlab original to
confirm that the transcription was good.
Note in the grid search in this version of testdisp, the grid is no longer
going over R2E, as the funmode is R2g=R2e. Note also that I've jigged up the
Pb in the attached code to cover higher pbs (we now go to 40%) to increase
the apparent errors in the CR equation. This is well outside anywhere
experimentally accessible, but fun to observe.
So your function is already faster by a factor of 2 which is good. Please can
you check the function to make sure i've wired your definitions up to my
function correctly? When I use a narrow grid search, the function seems to
give the exact result so I think it's probably right.
Here's the o/p:
####
Relative times of the various approaches:
Baldwin 0.194165 1.0
Numerical 21.562522 111.052568692
Meiboom 0.070515 0.363170499318
Carver 0.105395 0.542811526279
compareTroells.py:396: RuntimeWarning: invalid value encountered in log
R2eff=(1/Tc)*numpy.log(intensity0/intensity); # we did not factor
out (Ka+Kb)/2 here
Nikolai 0.615217 3.16852676847
compareTroells.py:274: RuntimeWarning: invalid value encountered in log
R2eff=-inv_relax_time * numpy.log(Mx)
NikolaiFun 0.376738 1.94029819998
NikolaiLong 81.447843 419.477470193
NikolaiLong 82.098139 422.82666289
Maximum error over grid (s-1):
Meiboom: 237214256.582
Baldwin: 6.57168541807e-10
CarverRichards: 29.9287668504
Nikolai fun: 2964.41300775
Nikolai dougle (9): 142.46676955
Nikolai long double (18): 0.063449549192819563320223
Nikolai long double (23): 5.7775487944482434059452e-7
####
I was hoping that the modification would make it more stable. If anything the
reverse is true. I would have included a note to this effect in the paper
(and still will if you can show it's better). Note the error is very
temperamental. You can miss it even by setting up the grid search
differently. Give it a narrow grid, and you won't see the error. Only
relatively specific, but seemingly hard to predict sets of parameters (all
with large deltaOs) lead it to explode.
Best,
Andy.
On 06/05/2014 08:40, Edward d'Auvergne wrote:
Hi Troels and Andy,
Troels, for your modification to compareTroells.py
(compareTroells_w_Troels_mod.py,
https://gna.org/support/download.php?file_id=20642), the change to the
ExchangeMat2x2() function requires a further change to the code:
$ diff -u compareTroells_w_Troels_mod.py compareTroells_w_Troels_mod_ed.py
--- compareTroells_w_Troels_mod.py 2014-05-05 16:06:06.522388615 +0200
+++ compareTroells_w_Troels_mod_ed.py 2014-05-06 09:17:00.328524062 +0200
@@ -148,10 +148,9 @@
if(verb=='n'):
return
else:
- #add on R2g to R2eff
arrayNew=[]
for i in range(len(array)):
- arrayNew.append((array[i][0],array[i][1]+R2g))
+ arrayNew.append((array[i][0],array[i][1]))
array=arrayNew
if(outfile!='Null'):
outy=open(outfile,'w')
The results I see are then, for compareTroells.py:
Relative times of the various approaches:
Baldwin 1.135749 1.0
Numerical 150.26362 132.303545942
Meiboom 0.40303 0.354858335777
Carver 0.610853 0.537841547736
Nikolai 3.33659 2.93778819088
Maximum error over grid (s-1):
Meiboom: 93015865.9296
Baldwin: 7.20705273238e-10
CarverRichards: 4989.8296812
Nikolai dougle (9): 204.761583779
Nikolai long double (18):
For compareTroells_w_Troels_mod.py:
Relative times of the various approaches:
Baldwin 1.159897 1.0
Numerical 150.864266 130.066950772
Meiboom 0.399363 0.344309020542
Carver 0.609061 0.525099211395
Nikolai 3.302521 2.84725367856
NikolaiLong 382.673914 329.920599846
Maximum error over grid (s-1):
Meiboom: 22664942.3388
Baldwin: 10.0000000007
CarverRichards: 18.2238852463
Nikolai dougle (9): 214.761583779
Nikolai long double (23): 10.000000206721384181558
And for the above change (compareTroells_w_Troels_mod_ed.py,
https://gna.org/support/download.php?file_id=20646):
Relative times of the various approaches:
Baldwin 1.157454 1.0
Numerical 151.098834 130.544137391
Meiboom 0.391843 0.338538723785
Carver 0.603573 0.521466079861
Nikolai 3.328743 2.87591817904
NikolaiLong 385.616089 333.158889252
Maximum error over grid (s-1):
Meiboom: 22664952.3388
Baldwin: 7.41007255556e-10
CarverRichards: 8.22388524626
Nikolai dougle (9): 204.761583779
Nikolai long double (23): 2.0672138595791449231557e-7
That CarverRichards number now looks what I would personally expect
from the theory! I don't like the Nikolai results, but I also don't
like the fact that I have two completely different implementations of
this code (from the sim_all software and the compareTroells.py
script).
Regards,
Edward
On 5 May 2014 16:27, Troels Emtekær Linnet <tlinnet@xxxxxxxxxxxxx> wrote:
Hi Andy.
Heh, the funniest part with the fluorescence was to battle with the
moderators of wikipedia to change 9000 to 9 in R06.
http://en.wikipedia.org/wiki/F%C3%B6rster_resonance_energy_transfer#Theoretical_basis
That was a brawl, to convince them that the bible was wrong. But
referring to the old dusty library book of Forster (Ref 14) made the
deal.
Back to the beer.
You are swapping two things in your code, according to CR72.
If I write up CR72: (Directly from cr72.py)
#######
fact = r20a - r20b - k_BA + k_AB
zeta = 2.0*dw * fact
So in total:
%%%%%
zeta(CR72) = 2.0*dw (R2A - R2B + k_AB - k_BA)
%%%%
Then we take your code:
#######
deltaR2=R2e-R2g
g1=2*dw*(deltaR2+keg-kge) #same as carver richards zeta
Change notation, to make it easier to see:
%%%%%
zeta(B14) = 2.0*dw (R2B - R2A + k_BA - k_AB)
%%%%
But huh, is that not just:
zeta(B14) = - zeta(CR72)
= 2.0*dw ( - R2A + R2B - k_AB + k_BA)
= 2.0*dw ( R2B - R2A + k_BA - k_AB)
Well, that is interesting.
So your g1 is the negative of carver richards zeta.
Fortunately we are lucky, that zeta is to the power **2.
So that leaves it out.
How about psi?
It boils down to that
fact(B14) = - fact(CR72) <=>
alpha_ (B14) = - alpha_ (CR72)
But alpha_ is also to the power **2. Phew...
So what is up.
If we look at: (compareTroells.py)
https://gna.org/support/download.php?file_id=20636
##############
def CarverRichards(kex,pb,dw,ncyc,Trelax,R2g,R2e,outfile,verb='n'):
deltaR2=R2e-R2g
keg=(1-pb)*kex
kge=pb*kex
zeta=2*dw*(-deltaR2+keg-kge) #zeta=g1
psi=(deltaR2+keg-kge)**2.0+4*keg*kge-dw**2 #psi=g2
#############
Well here we have something !
zeta = 2.0 * dw ( R2A - R2B - k_AB + k_BA )
So this is not about how you make the difference of ( R2A - R2B )
it is the swapping of + k_AB - k_BA << - k_AB + k_BA
kex(1 - 2pa) << kex( -1 + 2pa )
You have implemented CR72 wrong ?
Hmmmm. What is up?
You get a difference between CR72, and numerical.
You cannot get CR72 to work, unless you modify alpha_.
You invert the signs of forward / backward exchange rate.
You say that numerical is true, so CR72 is wrong.
Lets inspect numerical.
If we look at: (compareTroells.py)
https://gna.org/support/download.php?file_id=20636
Now we look up the numerical solution.
#####
NumDisp(kex,pb,dw,ncyc,Trelax,R2g,R2e,outfile,verb='n'):
array=CPMG2x2(dw,(R2e-R2g),kex,pb,Trelax,ncyc)
CPMG2x2(dOmega,DeltaR2,kex,pb,Trel,ncyc):
L=ExchangeMat2x2(dOmega,DeltaR2,kex,pb)
ExchangeMat2x2(dOmega,DeltaR2,kex,pb):
L = mat(zeros((2, 2),dtype=complex))
L[1,1] -= DeltaR2
Hmmm.
According to your paper at equation 3.
R+ =
[[ -kAB - R2A] [kBA] ]
[ kAB] [ -kBA - R2B - i dw] ]
But I see that it gets populated like
R+ =
[[ -kAB ] [kBA] ]
[ kAB] [ -kBA - R2B + R2A - i dw] ]
Has this something to do with it:
????
So I took the liberty to change your code, hoping the Theory is right.
############# Normal
Baldwin 1.183188 1.0
Numerical 58.556647 49.4905687008
Meiboom 0.391873 0.331200958766
Carver 0.693676 0.586277075156
compareTroells_w_Troels_mod.py:260: RuntimeWarning: invalid value
encountered in log
R2eff=(1/Tc)*numpy.log(intensity0/intensity); # we did not
factor out (Ka+Kb)/2 here
Nikolai 3.434938 2.90312105938
NikolaiLong 585.485493 494.837247335
Maximum error over grid (s-1):
Meiboom: 93015865.9296
Baldwin: 1.49539403083e-09
CarverRichards: 4989.8296812
Nikolai dougle (9): 204.761587913
Nikolai long double (23): 2.0692254415912185547869e-7
############ Modifying Meiboom to 4 NCYC
Maximum error over grid (s-1):
Meiboom: 22664952.3388
Baldwin: 1.49539403083e-09
CarverRichards: 4989.8296812
########### Modifying ExchangeMat2x2 to your paper version.
Meiboom: 22664942.3388
Baldwin: 10.0000000012
CarverRichards: 4979.8296812
########## Modify CR72 zeta=2*dw*(-deltaR2-keg+kge)
Baldwin 1.216869 1.0
Numerical 60.302727 49.555644034
Meiboom 0.397395 0.326571718073
Carver 0.700189 0.575402118059
compareTroells_w_Troels_mod.py:261: RuntimeWarning: invalid value
encountered in log
R2eff=(1/Tc)*numpy.log(intensity0/intensity); # we did not
factor out (Ka+Kb)/2 here
Nikolai 3.569577 2.93341107383
NikolaiLong 591.459824 486.050531323
Maximum error over grid (s-1):
Meiboom: 22664942.3388
Baldwin: 10.0000000012
CarverRichards: 18.2238852463
Nikolai dougle (9): 214.761587913
Nikolai long double (23): 10.000000207062814176945
WOOOOOWWWW
The change I made was:
##############################################
37c37
< def ExchangeMat2x2(dOmega,DeltaR2,kex,pb):
---
def ExchangeMat2x2(dOmega,R2e,R2g,kex,pb):
47c47
< L[1,1] -= DeltaR2
---
L[1,1] -= R2e
49a50
L[0, 0] -= R2g
64,65c65,66
< def CPMG2x2(dOmega,DeltaR2,kex,pb,Trel,ncyc):
< L=ExchangeMat2x2(dOmega,DeltaR2,kex,pb)
---
def CPMG2x2(dOmega,R2e,R2g,kex,pb,Trel,ncyc):
L=ExchangeMat2x2(dOmega,R2e,R2g,kex,pb)
86,87c87,88
< tcp=1/(2*nu_cpmg)
<
R2eff=(1-pb)*R2g+pb*R2e+R2Fast*(1.-(2./(kex*tcp))*numpy.tanh(kex*tcp/2.0))
---
tcp=1/(4*nu_cpmg)
R2eff=(1-pb)*R2g+pb*R2e+R2Fast*(1.-(4./(kex*tcp))*numpy.tanh(kex*tcp/4.0))
112c113
< zeta=2*dw*(-deltaR2+keg-kge) #zeta=g1
---
zeta=2*dw*(-deltaR2-keg+kge) #zeta=g1
145c146
< array=CPMG2x2(dw,(R2e-R2g),kex,pb,Trelax,ncyc)
---
array=CPMG2x2(dw,R2e,R2g,kex,pb,Trelax,ncyc)
546c547
< #tom,jarr=TestDisp(ncyc,Trelax,R2g,'NikolaiLong',tim)
---
tom,jarr=TestDisp(ncyc,Trelax,R2g,'NikolaiLong',tim)
560c561
< print 'Nikolai long double (18):
',numpy.max(numpy.abs(narr[:,:,1]-larr[:,:,1]))
---
#print 'Nikolai long double (18):
',numpy.max(numpy.abs(narr[:,:,1]-larr[:,:,1]))
[tlinnet@tomat ~/Downloads]$ diff compareTroells.py
compareTroells_w_Troels_mod.py
37c37
< def ExchangeMat2x2(dOmega,DeltaR2,kex,pb):
---
def ExchangeMat2x2(dOmega,R2e,R2g,kex,pb):
47c47
< L[1,1] -= DeltaR2
---
L[1,1] -= R2e
49a50
L[0, 0] -= R2g
64,65c65,66
< def CPMG2x2(dOmega,DeltaR2,kex,pb,Trel,ncyc):
< L=ExchangeMat2x2(dOmega,DeltaR2,kex,pb)
---
def CPMG2x2(dOmega,R2e,R2g,kex,pb,Trel,ncyc):
L=ExchangeMat2x2(dOmega,R2e,R2g,kex,pb)
86,87c87,88
< tcp=1/(2*nu_cpmg)
<
R2eff=(1-pb)*R2g+pb*R2e+R2Fast*(1.-(2./(kex*tcp))*numpy.tanh(kex*tcp/2.0))
---
tcp=1/(4*nu_cpmg)
R2eff=(1-pb)*R2g+pb*R2e+R2Fast*(1.-(4./(kex*tcp))*numpy.tanh(kex*tcp/4.0))
112c113
< zeta=2*dw*(-deltaR2+keg-kge) #zeta=g1
---
zeta=2*dw*(-deltaR2-keg+kge) #zeta=g1
145c146
< array=CPMG2x2(dw,(R2e-R2g),kex,pb,Trelax,ncyc)
---
array=CPMG2x2(dw,R2e,R2g,kex,pb,Trelax,ncyc)
546c547
< #tom,jarr=TestDisp(ncyc,Trelax,R2g,'NikolaiLong',tim)
---
tom,jarr=TestDisp(ncyc,Trelax,R2g,'NikolaiLong',tim)
560c561
< print 'Nikolai long double (18):
',numpy.max(numpy.abs(narr[:,:,1]-larr[:,:,1]))
---
#print 'Nikolai long double (18):
',numpy.max(numpy.abs(narr[:,:,1]-larr[:,:,1]))
[tlinnet@tomat ~/Downloads]$ diff compareTroells.py
compareTroells_w_Troels_mod.py
37c37
< def ExchangeMat2x2(dOmega,DeltaR2,kex,pb):
---
def ExchangeMat2x2(dOmega,R2e,R2g,kex,pb):
47c47
< L[1,1] -= DeltaR2
---
L[1,1] -= R2e
49a50
L[0, 0] -= R2g
64,65c65,66
< def CPMG2x2(dOmega,DeltaR2,kex,pb,Trel,ncyc):
< L=ExchangeMat2x2(dOmega,DeltaR2,kex,pb)
---
def CPMG2x2(dOmega,R2e,R2g,kex,pb,Trel,ncyc):
L=ExchangeMat2x2(dOmega,R2e,R2g,kex,pb)
86,87c87,88
< tcp=1/(2*nu_cpmg)
<
R2eff=(1-pb)*R2g+pb*R2e+R2Fast*(1.-(2./(kex*tcp))*numpy.tanh(kex*tcp/2.0))
---
tcp=1/(4*nu_cpmg)
R2eff=(1-pb)*R2g+pb*R2e+R2Fast*(1.-(4./(kex*tcp))*numpy.tanh(kex*tcp/4.0))
112c113
< zeta=2*dw*(-deltaR2+keg-kge) #zeta=g1
---
zeta=2*dw*(-deltaR2-keg+kge) #zeta=g1
145c146
< array=CPMG2x2(dw,(R2e-R2g),kex,pb,Trelax,ncyc)
---
array=CPMG2x2(dw,R2e,R2g,kex,pb,Trelax,ncyc)
546c547
< #tom,jarr=TestDisp(ncyc,Trelax,R2g,'NikolaiLong',tim)
---
tom,jarr=TestDisp(ncyc,Trelax,R2g,'NikolaiLong',tim)
560c561
< print 'Nikolai long double (18):
',numpy.max(numpy.abs(narr[:,:,1]-larr[:,:,1]))
---
#print 'Nikolai long double (18):
',numpy.max(numpy.abs(narr[:,:,1]-larr[:,:,1]))
####################################
File can be downloaded from:
https://gna.org/support/download.php?file_id=20642
What have we (maybe) learned?????
1) Baldwin has created a code, which is as good as the numerical of
Nicolai. Thanks !!! Wuhuuu !
2) CR72 is still valid though. Doing well.
And now for beer counting:
1) Appendix 1 – recipe for exact calculation of R2,eff
is: h1 = 2.0*dw (R2B - R2A + k_BA - k_AB) = - CR72 (But it have no
influence, but please keep history! Change it back to CR72. :-))
same for h2.
2) zeta typo in eq 70. z=zeta
3) Watch alpha_ in eq 70. Get Delta R straight
4) Comments in CODE:
http://svn.gna.org/viewcvs/*checkout*/relax/trunk/lib/dispersion/b14.py?revision=HEAD
g1=2*dw*(deltaR2+keg-kge) #same as carver richards zeta
That is not true. That is minus Zeta. But it does not matter, since order **2.
g2=(deltaR2+keg-kge)**2.0+4*keg*kge-dw**2 #same as carver richards psi
Not totally correct, but since deltaR2 is swapped, it gets minus
alpha, and goes to order 2.
5) Did you actually implement ExchangeMat2x2 wrong??
If yes, theeeeeen, no wonder figure 1 shows what it does. :-)
What is up for the deltaOmegas story???
This is my 2 cents. :-)
Let me hear if my hours have been wasted.
Best
Troels
2014-05-05 11:00 GMT+02:00 Andrew Baldwin <andrew.baldwin@xxxxxxxxxxxxx>:
I had to struggle long to find out the Forster distance calculation in
Lakowicz was wrong. The bible of fluorescence...
Very good stuff :) Most would not take the time to check :) The finest
biophysics really does demand both first rate biology and first rate
physics.
I don't think anyone has ever used the CR equation to analyse data outside
of the limit where R2A=R2B, so it's very unlikely that this has error has
had any real consequence. Note also that in general, it's actually much
easier to code up a numerical solution - ie once you've typed up an
evolution matrix, all you have left to do is write some lines taking matrix
exponentials. Whereas with an equation, not only are you making implicit
approximations that aren't necessarily obvious, there are many, many, more
lines to type, thus increasing the probability of a foobar. So numerical
solutions for spin dynamics = no mistakes, and also easy to code! This has
been the approach in the Kay group, arguably cultivated by Dimitri and
Nikolai, and then pushed hard more recently by Flemming, Pramodh and
Guillaume.
This paper was truly groundbreaking (in my view):
Journal of biomolecular NMR. 2000 18(1):49-63
On this note, there are a few suggestions that Flemming and I are putting
together for another thread sensibly started by Edward a few days ago for
analysing R1rho and CPMG data. I think that with just a few extra tweaks,
the backend to your code can be as powerful and as versatile as your front
end, for millisecond data analysis. I agree that it really is in the
interests of the field to have people able to perform the experiments and
get information out of their data as efficiently (and as rigorously) as
possible.
This means more citations for those that made the experiments, and the
community can keep showing more examples where Xray people have
inadvertently missed functionally relevant conformers :)
Best,
Andy.
On 04/05/2014 23:09, Troels Emtekær Linnet wrote:
Well Andrew.
Your paper popped out of the Mail alerts. :-)
My supervisor Kaare Teilum guided my attention to it, and asked me to
go through it to check.
I discussed it with a lab colleague, and we were thrived.
The math derivations still kills me, but I thought that the essential
results was alarming.
I mean, half of our lab do CPMG or R1rho.
Any changes to those equations gets our attention.
The "worst" problem we have in our lab, is the ability to analyse data.
With relax, we now have a GUI, which make sure that at least the
Bachelor Students still survive.
And what else can I then do than get dirty with your code, to test if
it is right?
If you suspicion about the CR72 is correct, then I really really
wonder why not errata article is out there?
This situation reminds me of my master thesis.
I had to struggle long to find out the Forster distance calculation in
Lakowicz was wrong. The bible of fluorescence...
Jesus....
but thanks for the beer!
Best
Your biggest fan. ;-)
2014-05-04 22:59 GMT+02:00 Andrew Baldwin <andrew.baldwin@xxxxxxxxxxxxx>:
Dear Troels,
The existence of typos is exceedingly probable! The editors have added
more
errors in the transcrpition, so it'll be curious to see how accurate the
final draft will be. Your comments come just in time for alterations.
So I'll pay a pint of beer for every identified error. I didn't use that
particular supplementary section to directly calculate anything, so this
is
the place in the manu where errors are most likely to occur. I do note
the
irony that this is probably one of few bits in the article that people
are
likely to read.
1) Yep - that's a typo. alpha- should definitely be KEG-KGE (score 1
pint).
2) Yep! zs should be zetas (score 1 pint).
3) This one is more interesting. Carver and Richard's didn't define a
deltaR2. To my knowledge, I think Nikolai was the first to do that in his
2002 paper on reconstructing invisible states. So here, I think you'll
find
that my definition is the correct one. There is a typo in the Carver
Richard's equation. I've just done a brief lit review, and as far as I
can
see, this mistake has propagated into every transcription of the equation
since. I should note that in the manu in this section. Likely this has
gone
un-noticed as differences in R2 are difficult to measure by CPMG, instead
R2A=R2B being assumed. Though thanks to Flemming, it is now possible to
break this assumption:
J Am Chem Soc. 2009 Sep 9;131(35):12745-54. doi: 10.1021/ja903897e.
Also, in the Ishima and Torchia paper looking at the effects of
relaxation
differences:
Journal of biomolecular NMR. 2006 34(4):209-19
I understand their derivations started from Nikolai's formula, which is
exact, so the error wouldn't come to light there either, necessarily.
To prove the point try this: numerically evaluate R2,effs over a big grid
of
parameters, and compare say Nikolai's formula to your implementation of
Carver Richard's with a +deltaR2 and a -deltaR2 in the alpha- parameter.
Which agrees better? You'll see that dR2=R2E-R2G being positive wins.
The attached should also help. It's some test code I used for the paper
that
has been briefly sanitised. In it, the numerical result, Carver Richards,
Meiboom, Nikolai's and my formula are directly compared for speed
(roughly)
and precision.
The code added for Nikolai's formula at higher precision was contributed
by
Nikolai after the two of us discussed some very odd results that can come
from his algorithm. Turns out evaluating his algorithm at double
precision
(roughly 9 decimal places, default in python and matlab) is
insuffficient.
About 23 decimal places is necessary for exactness. This is what is
discussed in the supplementary info in the recent paper. In short, on my
laptop:
formula / relative speed / biggest single point error in s-1 over grid
Mine 1 7.2E-10
Meiboom 0.3 9E7
Carver Richards 0.53 8
Nikolai (double precision 9dp) 3.1 204
Nikolai (long double 18dp) 420 3E-2
Nikolai (23 dp) 420 2E-7
Numerical 118 0
The error in Carver Richard's with the alpha- sign made consistent with
the
Carver Richard's paper gives an error of 4989 s-1 in the same test. Ie,
it's
a bad idea.
Note that the precise calculations in python at arbitrary precision are
slow
not because of the algorithm, but because of the mpmath module, which
deals
with C code via strings. Accurate but slow. This code has a small amount
of
additional overhead that shouldn't strictly be included in the time. I've
no
doubt also that algorithms here could be sped up further by giving them a
bit of love. Those prefaces aside, I'm surprised that your benchmarking
last
Friday showed Nikolai's to be so much slower than mine? In my hands
(attached) at double precision, Nikolai's formula is only 3x slower than
mine using numpy. Though this is not adequate precision for mistake free
calculations.
Perhaps more physically compelling, to get my equation to reduce to
Carver
Richard's, alpha- has to be defined as I show it. Note that the origin of
this parameter is its role in the Eigenvalue of the 2x2 matrix. This is
the
one advantage of a derivation - you can see where the parameters come
from
and get a feel for what their meaning is.
So this typo, and the fact that Carver and Richards inadvertently change
their definition of tau_cpmg half way through the paper, are both errors
in
the original manuscript. Though given that proof reading and type-setting
in
1972 would have been very very painful, that there are this few errors is
still actually quite remarkable. I notice in your CR72 code, you are
currently using the incorrect definition, so I would suggest you change
that.
So overall, I think you're still ahead, so I'm calling that 2 pints to
Troels.
I would actually be very grateful if you can bring up any other
inconsistencies in the paper! Also from an editorial perspective, please
feel free to note if there are parts of the paper that were particularly
unclear when reading. If you let me know, there's still time to improve
it.
Thusfar, I think you are the third person to read this paper (the other
two
being the reviewers who were forced to read it), and you're precisely its
target demographic!
Best,
Andy.
On 03/05/2014 13:54, Troels Emtekær Linnet wrote:
Dear Andrew Baldwin.
I am in the process of implementing your code in relax.
I am getting myself into a terrible mess by trying to compare
equations in papers, code and conversions.
But I hope you maybe have a little time to clear something up?
I am looking at the current version of your paper:
dx.doi.org/10.1016/j.jmr.2014.02.023
which still is still "early version of the manuscript. The manuscript
will undergo copyediting, typesetting, and review of the resulting
proof before it is published in its final form."
1) Equations are different ?
In Appendix 1 – recipe for exact calculation of R2,eff
h1 = 2 dw (dR2 + kEG - kGE)
But when I compare to:
Supplementary Section 4. Relation to Carver Richards equation. Equation
70.
h1 = zeta = 2 dw alpha_ = 2 dw (dR2 + kGE - kEG)
Is there a typo here? The GE and EG has been swapped.
2) Missing zeta in equation 70
There is missing \zeta instead of z in equation 70, which is probably
a typesetting error.
3) Definition of " Delta R2" is opposite convention?
Near equation 10, 11 you define: delta R2 = R2e - R2g
if we let e=B and g=A, then delta R2 = R2B - R2A
And in equation 70:
alpha_ = delta R2 + kGE - kEG
but i CR original work, A8,
alpha_ = R2A - R2B + kA - kB = - delta R2 + kGE - kEG
So, that doesn't match?
Sorry if these questions are trivial.
But I hope you can help clear them up for me. :-)
Best
Troels
2014-05-01 10:07 GMT+02:00 Andrew Baldwin <andrew.baldwin@xxxxxxxxxxxxx>:
The Carver and Richards code in relax is fast enough, though Troels
might have an interest in squeezing out a little more speed. Though
it would require different value checking to avoid NaNs, divide by
zeros, trig function overflows (which don't occur for math.sqrt), etc.
Any changes would have to be heavily documented in the manual, wiki,
etc.
The prove that the two definitions are equivalent is relatively
straightforward. The trig-to-exp command in the brilliant (and free)
symbolic program sage might prove helpful in doing that.
For the tau_c, tau_cp, tau_cpmg, etc. definitions, comparing the relax
code to your script will ensure that the correct definition is used.
That's how I've made sure that all other dispersion models are
correctly handled - simple comparison to other software. I'm only
aware of two definitions though:
tau_cpmg = 1 / (4 nu_cpmg)
tau_cpmg = 1 / (2 nu_cpmg)
tau_cpmg = 1/(nu_cpmg).
Carver and Richard's switch definitions half way through the paper
somewhere.
What is the third? Internally in relax, the 1 / (4 nu_cpmg)
definition is used. But the user inputs nu_cpmg. This avoids the
user making this mistake as nu_cpmg only has one definition.
I've seen people use nu_cpmg defined as the pulse frequency. It's just an
error that students make when things aren't clear. I've often seen brave
student from a lab that has never tried CPMG before do this. Without
someone
to tell them that this is wrong, it's not obvious to them that they've
made
a mistake. I agree with you that this is solved with good documentation.
You guys are free to use my code (I don't mind the gnu license) or of
course
implement from scratch as needed.
Cheers! For a valid copyright licensing agreement, you'll need text
something along the lines of:
"I agree to licence my contributions to the code in the file
http://gna.org/support/download.php?file_id=20615 attached to
http://gna.org/support/?3154 under the GNU General Public Licence,
version three or higher."
Feel free to copy and paste.
No problem:
"I agree to licence my contributions to the code in the file
http://gna.org/support/download.php?file_id=20615 attached to
http://gna.org/support/?3154 under the GNU General Public Licence,
version three or higher."
I'd like to note again though that anyone using this formula to fit data,
though exact in the case of 2site exchange/inphase magnetisation,
evaluated
at double floating point precision should not be doing so! Neglecting
scalar
coupling/off resonance/spin flips/the details of the specific pulse
sequence
used will lead to avoidable foobars. I do see value in this, as described
in
the paper, as being a feeder for initial conditions for a more fancy
implemenation. But beyond that, 'tis a bad idea. Ideally this should
appear
in big red letters in your (very nice!) gui when used.
In relax, we allow the user to do anything though, via the
auto-analysis (hence the GUI), we direct the user along the best path.
The default is to use the CR72 and 'NS CPMG 2-site expanded' models
(Carver & Richards, and Martin Tollinger and Nikolai Skrynnikov's
Maple expansion numeric model). We use the CR72 model result as the
starting point for optimisation of the numeric models, allowing a huge
speed up in the analysis. The user can also choose to not use the
CR72 model results for the final model selection - for determining
dispersion curve significance.
Here's the supp from my CPMG formula paper (attached). Look at the last
figure. Maybe relevant. Nikolai's algorithm blows up sometimes when you
evaluate to double float point precision (as you will when you have a
version in python or matlab). The advantage of Nicolai's formula, or mine
is
that they won't fail when Pb starts to creep above a per cent or so.
Using the simple formula as a seed for the more complex on is a good
idea.
The most recent versions of CATIA have something analogous.
As for the scalar coupling and spin flips, I am unaware of any
dispersion software that handles this.
CATIA. Also cpmg_fig I believe. In fact, I think we may have had this
discussion before?
https://plus.google.com/s/cpmg%20glove
If you consider experimental validation a reasonable justification for
inclusion of the effects then you might find this interesting:
Spin flips are also quite relevant to NH/N (and in fact any spin system).
The supplementary section of Flemming and Pramodh go into it here for
NH/N
http://www.pnas.org/content/105/33/11766.full.pdf
And this:
JACS (2010) 132: 10992-5
Figure 2:
r2 0.97, rmsd 8.0 ppb (no spin flips)
r2 0.99, rmsd 5.7 ppb (spin flips).
The improvements are larger than experimental uncertainties.
When we design these experiments and test them, we need to think about
the
details. This is in part because Lewis beats us if we don't. You can
imagine
that it comes as a surprise then when we see people neglecting this. In
short, the parameters you extract from fitting data will suffer if the
details are not there. In the case of spin flips, the bigger the protein,
the bigger the effect. In your code, you have the opportunity to do
things
properly. This leaves the details behind the scenes, so the naive user
doesn't have to think about them.
And only Flemming's CATIA
handles the CPMG off-resonance effects. This is all explained in the
relax manual. I have asked the authors of most dispersion software
about this too, just to be sure. I don't know how much of an effect
these have though. But one day they may be implemented in relax as
well, and then the user can perform the comparison themselves and see
if all these claims hold.
Myint, W. & Ishima, R. Chemical exchange effects during refocusing pulses
in
constant-time CPMG relaxation dispersion experiments. J Biomol Nmr 45,
207-216, (2009).
also:
Bain, A. D., Kumar Anand, C. & Nie, Z. Exact solution of the CPMG pulse
sequence with phase variation down the echo train: application to R(2)
measurements. J Magn Reson 209, 183- 194, (2011).
Or just simulate the off-resonance effects yourself to see what happens.
For
NH you see the effects clearly for glycines and side chains, if the
carrier
is in the centre around 120ppm. The problem generally gets worse the
higher
field you go though this of course depends when you bought your probe.
You
seem to imply that these effects are almost mythical. I assure you that
they
come out happily from the Bloch equations.
Out of curiosity, from a philosophical perspective, I wonder if you'd
agree
with this statement:
"the expected deviations due to approximations in a model should be lower
than the experimental errors on each datapoint."
Anyway, the warnings about analytic versers numeric are described in
the manual. But your new model which adds a new factor to the CR72
model, just as Dimitry Korzhnev's cpmg_fit software does for his
multiple quantum extension of the equation (from his 2004 and 2005
papers), sounds like it removes the major differences between the
analytic and numeric results anyway. In any case, I have never seen
an analytic result which is more than 10% different in parameter
values (kex specifically) compared to the numeric results. I am
constantly on the lookout for a real or synthetic data set to add to
relax to prove this wrong though.
I think there's a misunderstanding in what I mean by numerical modelling.
For the 2x2 matrix (basis: I+, ground and excited) from the Bloch
McConnell
equation, you can manipulate this to get an exact solution. Nikolai's
approach does this, though his algorithm can trip up sometimes for
relatively exotic parameters when you use doubles (see attached). My
formula
also does this. I agree with you entirely that in this instance,
numerically
solving the equations via many matrix exponentials is irrelevant as
you'll
get an identical result to using a formula.
My central thesis here though is that to get an accurate picture of the
experiment you need more physics. This means a bigger basis. To have off
resonance effects, you need a minimum 6x6 (Ix,Iy,Iz, ground and excited).
To
include scalar coupling and off resonance, you need a 12x12
(Ix,Iy,Iz,IxHz,IyHz,IzHz, ground and excited). Including R1 means another
element and so on. The methyl group, for example, means you need 24.
So when we use the term numerical treatment, we generally mean a
calculation
in a larger basis, as is necessary to take into account the appropriate
spin
physics. There aren't formulas to deal with these situations. In the 6x6
case for example, you need 6 Eigenvalues, which is going to make life
very
rough for someone brave enough to attempt a close form solution. The
Palmer
and Trott trick used in 2002 for R1rho is a smart way of ducking the
problem
of having 6 Eigenvalues, but for CPMG unfortunately you need several
Eigenvalues, not just the smallest.
The 2x2 matrix that Nikolai and Martin, Carver Richard's and I analyse
does
not include scalar coupling, as magnetisation is held in-phase (in
addition
to neglecting all the other stuff detailed above). So it is a good
representation for describing the continuous wave in-phase experiments
introduced here (neglecting relaxation effects and off resonance):
Vallurupalli, P.; Hansen, D. F.; Stollar, E. J.; Meirovitch, E.; Kay, L.
E.
Proc. Natl. Acad. Sci. U.S.A. 2007, 104, 18473–18477.
and here:
Baldwin, A. J.; Hansen, D. F.; Vallurupalli, P.; Kay, L. E. J. Am. Chem.
Soc. 2009, 131, 11939–48.
But I think are the only two where this formula is directly applicable.
Only
if you have explicit decoupling during the CPMG period do you satisfy
this
condition. So this is not the case for all other experiments and
certainly
not true for those used most commonly.
Anyhow. Best of luck with the software. I would recommend that you
consider
implementing these effects and have a look at some of the references. The
physics are fairly complex, but the implementations are relatively
straightforward and amount to taking many matrix exponentials. If you do
this, I think you'd end up with a solution that really is world-leading.
As it stands though, in your position, I would worry that on occasion,
users
will end up getting slightly wrong parameters out from your code by
neglecting these effects. If a user trusts this code then, in turn, they
might lead themselves to dodgy biological conclusions. For the time
being,
I'll stick to forcing my students to code things up themselves.
All best wishes,
Andy.
#!/usr/bin/python
######################################################################
# Script to simulate different forms of the CPMG intensity,
# and to compare approaches
#
# Started
# A.Baldwin 20th October 2011
# Simplified 30th Jan 2014
#
# NikolaiLong function coded up using mpmath module. This was kindly
# contributed by Nikolai. The python implementation is going via
# strings and so is hell-of-slow. It does give the user full control
# of precision. Note mp.mp.dps>=23 required for numerical precision
# 18 is long double precision and gives acceptable errors.
# 9 is double precision, and can give foobars.
# March 2014.
#
# Copyright (c) A.Baldwin.
# University of Oxford.
# Please do not circulate without permission.
#
import numpy,sys
from math import cos,sin,atan2,log10
from datetime import datetime,timedelta
from scipy import mat,zeros,log
from scipy.linalg import expm
import mpmath as mp
mp.mp.dps=23
mp.mp.dps=18
#######################################################################
#Numerically take eigenvalues and propagate the 2x2 matrix
#
#setup 2x2 matrix
def ExchangeMat2x2(dOmega,DeltaR2,kex,pb,R2g):
#rates in s-1, domage in rad s-1, dfrq in MHz
k_ab=pb*kex
k_ba=(1-pb)*kex
L = mat(zeros((2, 2),dtype=complex))
#add chemical shift evolution
L[0,0]=complex(0,0.0)
L[1,1]=complex(0,-dOmega)
#add intrinsic relaxation
L[0,0] -= R2g
L[1,1] -= DeltaR2
#add exchange
L[0, 0] -= k_ab
L[1, 0] += k_ab
L[1, 1] -= k_ba
L[0, 1] += k_ba
return L
def make_free_precess_propagator(L, time):
return mat(expm(L * time))
def Make2x2Equil(pb):
I0=mat(zeros((2, 1),dtype=complex))
I0[0,0]=complex((1-pb),0.0)
I0[1,0]=complex(pb,0.0)
return I0
def CPMG2x2(dOmega,DeltaR2,kex,pb,Trel,ncyc,R2g):
L=ExchangeMat2x2(dOmega,DeltaR2,kex,pb,R2g)
R2eff=[]
nu_cpmg=(ncyc)/Trel
for i in range(len(ncyc)):
I0=Make2x2Equil(pb)
if(i==0):
Ifin=I0
tcp=Trel/(4*ncyc[i])
LP=make_free_precess_propagator(L,tcp) #make the tau propagator
LPncyc=(LP*LP.conjugate()*LP.conjugate()*LP)**ncyc[i] #make the
(tau/180/tau/tau/180/tau)*ncyc propagator
Ifin=LPncyc*I0
R2eff.append((nu_cpmg[i],(1/Trel)*log((1-pb)/Ifin[0].real)))
return R2eff
#######################################################################
#Meiboom formula
def Meiboom(kex,pb,dw,ncyc,Trelax,R2g,R2e,outfile,verb='y'):
R2Fast=pb*(1-pb)*dw**2.0/kex
nu_cpmg=ncyc/Trelax
tcp=1/(2*nu_cpmg)
R2eff=(1-pb)*R2g+pb*R2e+R2Fast*(1.-(2./(kex*tcp))*numpy.tanh(kex*tcp/2.0))
if(verb=='n'):
return
else:
array=[]
for i in range(len(ncyc)):
array.append((nu_cpmg[i],R2eff[i]))
if(outfile!='Null'):
outy=open(outfile,'w')
for i in range(len(array)):
outy.write('%f\t%f\n' % (array[i][0],array[i][1]))
outy.close()
return array
#######################################################################
#Carver Richard's equation
def CarverRichards(kex,pb,dw,ncyc,Trelax,R2g,R2e,outfile,verb='n'):
deltaR2=R2e-R2g
keg=(1-pb)*kex
kge=pb*kex
zeta=2*dw*(deltaR2+keg-kge) #zeta=g1
psi=(deltaR2+keg-kge)**2.0+4*keg*kge-dw**2 #psi=g2
NpFac=numpy.sqrt( psi+numpy.sqrt(psi**2.0+zeta**2.0))*2./numpy.sqrt(2.0)
NmFac=numpy.sqrt(-psi+numpy.sqrt(psi**2.0+zeta**2.0))*2./numpy.sqrt(2.0)
Dp=0.5*( 1.+(psi+2.*dw**2.0)/(numpy.sqrt(psi**2.0+zeta**2.0)))
Dm=0.5*(-1.+(psi+2.*dw**2.0)/(numpy.sqrt(psi**2.0+zeta**2.0)))
nu_cpmg=ncyc/Trelax
tcp=1/(4*nu_cpmg)
Np=NpFac*tcp
Nm=NmFac*tcp
R2eff=(R2g+R2e+kex)/2.0-(1/(4*tcp))*(numpy.arccosh(Dp*numpy.cosh(Np)-Dm*numpy.cos(Nm)))
if(verb=='n'):
return
else:
array=[]
for i in range(len(ncyc)):
array.append((nu_cpmg[i],R2eff[i]))
if(outfile!='Null'):
outy=open(outfile,'w')
for i in range(len(array)):
outy.write('%f\t%f\n' % (array[i][0],array[i][1]))
outy.close()
return array
#######################################################################
#numerical solution
def NumDisp(kex,pb,dw,ncyc,Trelax,R2g,R2e,outfile,verb='n'):
array=CPMG2x2(dw,(R2e),kex,pb,Trelax,ncyc,R2g)
if(verb=='n'):
return
else:
#add on R2g to R2eff
arrayNew=[]
for i in range(len(array)):
arrayNew.append((array[i][0],array[i][1]))
array=arrayNew
if(outfile!='Null'):
outy=open(outfile,'w')
for i in range(len(array)):
outy.write('%f\t%f\n' % (array[i][0],array[i][1]))
outy.close()
return array
#New fun mode.
def NikolaiFun(kex,pb,dw,ncyc,Trelax,R2g,R2e,outfile,verb='n'):
#rename variables
relax_time=Trelax;
inv_relax_time=1/Trelax
r20=R2g
pa=(1-pb)
pb=pb
pA=pa
pB=pb
k_AB=kex*pB;
k_BA=kex*pA;
nu_cpmg=ncyc/Trelax; # get cpmg frequencies
tcp=1./(2*nu_cpmg); # tcp is the interval between consequitive 180
pulses
num_cpmg=ncyc
j=complex(0,1)
# Repeditive calculations.
half_tcp = 0.5 * tcp
k_AB_plus_k_BA = k_AB + k_BA
k_BA_minus_k_AB = k_BA - k_AB
# The expansion factors (in numpy array form for all dispersion points).
t3 = 1.j
t4 = t3 * dw
two_t4 = 2.0 * t4
t5 = k_BA**2
t8 = two_t4 * k_BA
t10 = 2.0 * k_BA * k_AB
t11 = dw**2
t14 = two_t4 * k_AB
t15 = k_AB**2
t5_t10_t11_t15 = t5 + t10 - t11 + t15
t8_t14 = t8 - t14
t17 = numpy.sqrt(t5_t10_t11_t15 - t8_t14)
k_AB_plus_k_BA_minus_t4 = k_AB_plus_k_BA - t4
t21 = numpy.exp((t17 - k_AB_plus_k_BA_minus_t4) * half_tcp)
t22 = 1.0/t17
t28 = numpy.exp(-(t17 + k_AB_plus_k_BA_minus_t4) * half_tcp)
t31 = t22*k_AB * (t21 - t28)
t33 = numpy.sqrt(t5_t10_t11_t15 + t8_t14)
k_AB_plus_k_BA_plus_t4 = k_AB_plus_k_BA + t4
k_BA_minus_k_AB_plus_t4 = k_BA_minus_k_AB + t4
t34 = k_BA_minus_k_AB_plus_t4 + t33
t37 = numpy.exp((t33 - k_AB_plus_k_BA_plus_t4) * tcp)
t39 = 1.0/t33
t41 = k_BA_minus_k_AB_plus_t4 - t33
t44 = numpy.exp(-(t33 + k_AB_plus_k_BA_plus_t4) * tcp)
t47 = 0.5*t39 * (t34*t37 - t41*t44)
k_BA_minus_k_AB_minus_t4 = k_BA_minus_k_AB - t4
t49 = k_BA_minus_k_AB_minus_t4 - t17
t51 = t21 * t49 * t22
t52 = k_BA_minus_k_AB_minus_t4 + t17
t54 = t28 * t52 * t22
t55 = -t51 + t54
t60 = 0.5*t39*k_AB * (t37 - t44)
t62 = t31*t47 + t55*t60
t63 = 1.0/k_AB
t68 = 0.5*t63 * (t49*t54 - t52*t51)
t69 = 0.5*t62 * t68
t72 = t37 * t41 * t39
t76 = t44 * t34 * t39
t78 = 0.5*t63 * (t41*t76 - t34*t72)
t80 = 0.5 * (t76 - t72)
t82 = 0.5 * (t31*t78 + t55*t80)
t83 = t82 * t55/2.0
t88 = 0.5 * t22 * (t52*t21 - t49*t28)
t91 = t88 * t47 + t68*t60
t92 = t91 * t88
t95 = 0.5 * (t88*t78 + t68*t80)
t96 = t95 * t31
t97 = t69 + t83
t98 = t97**2
t99 = t92 + t96
t102 = t99**2
t108 = t62 * t88 + t82 * t31
t112 = numpy.sqrt(t98 - 2.0 * t99 * t97 + t102 + 2.0 * (t91 * t68 + t95 *
t55) * t108)
t97_t99 = t97 + t99
t97_nt99 = t97 - t99
t113 = t97_nt99 - t112
t115 = num_cpmg
t116 = numpy.power(0.5*(t97_t99 + t112), t115)
t118 = 1.0/t112
t120 = t97_nt99 + t112
t122 = numpy.power(0.5*(t97_t99 - t112), t115)
t127 = 0.5/t108
t120_t122 = t120*t122
t139 = 0.5/(k_AB + k_BA) * ((t120_t122 - t113*t116)*t118*k_BA +
(t120_t122 - t116*t120)*t127*t113*t118*k_AB)
# Calculate the initial and final peak intensities.
intensity0 = pA
intensity = t139.real * numpy.exp(-relax_time * r20)
# The magnetisation vector.
Mx = intensity / intensity0
# Calculate the R2eff using a two-point approximation, i.e. assuming that
the decay is mono-exponential, and store it for each dispersion point.
R2eff=-inv_relax_time * numpy.log(Mx)
#for i in range(num_points):
# if Mx[i] <= 0.0 or isNaN(Mx[i]):
# back_calc[i] = 1e99
# else:
# back_calc[i]= -inv_relax_time * log(Mx[i])
if(verb=='n'):
return
else:
array=[]
for i in range(len(ncyc)):
array.append((nu_cpmg[i],R2eff[i]))
if(outfile!='Null'):
outy=open(outfile,'w')
for i in range(len(array)):
outy.write('%f\t%f\n' % (array[i][0],array[i][1]))
outy.close()
return array
#######################################################################
#coding up nikolai's formula
def Nikolai(kex,pb,dw,ncyc,Trelax,R2g,R2e,outfile,verb='n'):
#rename variables
Tc=Trelax;
Ra=R2g
Rb=R2e
dR=Ra-Rb;
pa=1-pb;
Ka=kex*pb;
Kb=kex*pa;
nu_cpmg=ncyc/Trelax; # get cpmg frequencies
tcp=1./(2*nu_cpmg); # tcp is the interval between consequitive 180
pulses
N=ncyc*2; # ncyc*2 (no need for round with these variables)
#off we go.
Rinf=(Ra+Rb)/2;
t3 = complex(0,1)/2;
t4 = t3*dw;
t5 = Ka/2;
t6 = Ra/2;
t7 = Kb/2;
t8 = Rb/2;
t9 = -2*(complex(0,1));
t14 = 2*(complex(0,1));
t20 = 2*Kb*Ka;
t22 = 2*Ka*Rb;
t24 = 2*Ra*Kb;
t26 = 2*Ra*Rb;
t28 = 2*Kb*Rb;
t30 = 2*Ka*Ra;
t31 = Rb**2;
t32 = Kb**2;
t33 = Ka**2;
t34 = Ra**2;
t35 = dw**2;
t36 =
t9*dw*Kb+t9*Rb*dw+t14*Ka*dw+t14*Ra*dw+t20-t22-t24-t26+t28+t30+t31+t32+t33+t34-t35;
t37 = numpy.sqrt(t36);
t38 = t37/2;
t42 = numpy.exp((t4-t5-t6-t7-t8+t38)*tcp/2);
t43 = 1/t37;
t49 = numpy.exp((t4-t5-t6-t7-t8-t38)*tcp/2);
t52 = t42*t43*Ka-t49*t43*Ka;
t61 =
t14*Rb*dw+t20-t22-t24-t26+t28+t30+t31+t32+t33+t34+t9*Ka*dw+t9*Ra*dw-t35+t14*dw*Kb;
t62 = numpy.sqrt(t61);
t63 = t62/2;
t64 = t4-t5-t6+t7+t8+t63;
t65 = -complex(0,1)/2;
t66 = t65*dw;
t69 = numpy.exp((t66-t5-t6-t7-t8+t63)*tcp);
t71 = 1/t62;
t73 = t4-t5-t6+t7+t8-t63;
t76 = numpy.exp((t66-t5-t6-t7-t8-t63)*tcp);
t79 = t64*t69*t71-t73*t76*t71;
t81 = complex(0,1);
t82 = t81*dw;
t85 = t42*(t82+Ka+Ra-Kb-Rb+t37)*t43;
t88 = t49*(t82+Ka+Ra-Kb-Rb-t37)*t43;
t89 = t85-t88;
t94 = t69*t71*Ka-t76*t71*Ka;
t96 = t52*t79+t89*t94/2;
t97 = t66-t5-t6+t7+t8+t38;
t98 = 1/Ka;
t101 = t66-t5-t6+t7+t8-t38;
t104 = t97*t98*t85-t101*t98*t88;
t105 = t96*t104/2;
t109 = t69*(t82-Ka-Ra+Kb+Rb-t62)*t71;
t114 = t76*(t82-Ka-Ra+Kb+Rb+t62)*t71;
t116 = -t64*t98*t109+t73*t98*t114;
t118 = -t109+t114;
t120 = t52*t116/2+t89*t118/4;
t121 = t120*t89/2;
t126 = t97*t42*t43-t101*t49*t43;
t129 = t126*t79+t104*t94/2;
t130 = t129*t126;
t133 = t126*t116/2+t104*t118/4;
t134 = t133*t52;
t135 = t105+t121;
t136 = t135**2;
t137 = t130+t134;
t140 = t137**2;
t146 = t96*t126+t120*t52;
t150 = numpy.sqrt(t136-2*t137*t135+t140+4*(t129*t104/2+t133*t89/2)*t146);
t151 = t105+t121-t130-t134-t150;
t153 = N/2;
t154 = (t105/2+t121/2+t130/2+t134/2+t150/2)**t153;
t156 = 1./t150;
t158 = t105+t121-t130-t134+t150;
t160 = (t105/2+t121/2+t130/2+t134/2-t150/2)**t153;
t165 = 1./t146;
t177 =
1/(Ka+Kb)*((-t151*t154*t156/2+t158*t160*t156/2)*Kb+(-t151*t165*t154*t158*t156/2+t158*t165*t160*t151*t156/2)*Ka/2);
intensity0 = pa ; # pA
intensity = numpy.real(t177);
R2eff=(1/Tc)*numpy.log(intensity0/intensity); # we did not factor
out (Ka+Kb)/2 here
if(verb=='n'):
return
else:
array=[]
for i in range(len(ncyc)):
array.append((nu_cpmg[i],R2eff[i]))
if(outfile!='Null'):
outy=open(outfile,'w')
for i in range(len(array)):
outy.write('%f\t%f\n' % (array[i][0],array[i][1]))
outy.close()
return array
def NikolaiLong(kex,pb,dw,ncyc,Trel,R2g,R2e,outfile,verb='n'):
#print 'Running test for exactness'
dfrq=mp.mpf('200.0') #spectrometer frequency of nuclei (MHz)
kex=mp.mpf(str(kex))
pb=mp.mpf(str(pb))
dw=mp.mpf(str(dw))
#rename variables
Tc=mp.mpf(str(Trel));
Ra=R2g
Rb=R2e
dR=Ra-Rb;
pa=1-pb;
Ka=kex*pb;
Kb=kex*pa;
R2eff=[]
nu_cpmg=[]
for i in range(len(ncyc)):
ncycV=mp.mpf(str(ncyc[i])) #number of cpmg cycles
nu_cpmgV=mp.mpf(str(ncycV))/Trel; # get cpmg frequencies
tcp=1/(2*nu_cpmgV); # tcp is the interval between consequitive
180 pulses
N=ncycV*2; # ncyc*2 (no need for round with these
variables)
nu_cpmg.append(nu_cpmgV.real) #off we go.
t3 = mp.mpc(0,1)/2;
t4 = t3*dw;
t5 = Ka/2;
t6 = Ra/2;
t7 = Kb/2;
t8 = Rb/2;
t9 = -2*(mp.mpc(0,1));
t14 = 2*(mp.mpc(0,1));
t20 = 2*Kb*Ka;
t22 = 2*Ka*Rb;
t24 = 2*Ra*Kb;
t26 = 2*Ra*Rb;
t28 = 2*Kb*Rb;
t30 = 2*Ka*Ra;
t31 = Rb**2;
t32 = Kb**2;
t33 = Ka**2;
t34 = Ra**2;
t35 = dw**2;
t36 =
t9*dw*Kb+t9*Rb*dw+t14*Ka*dw+t14*Ra*dw+t20-t22-t24-t26+t28+t30+t31+t32+t33+t34-t35;
t37 = mp.sqrt(t36);
t38 = t37/2;
t42 = mp.exp((t4-t5-t6-t7-t8+t38)*tcp/2);
t43 = 1/t37;
t49 = mp.exp((t4-t5-t6-t7-t8-t38)*tcp/2);
t52 = t42*t43*Ka-t49*t43*Ka;
t61 =
t14*Rb*dw+t20-t22-t24-t26+t28+t30+t31+t32+t33+t34+t9*Ka*dw+t9*Ra*dw-t35+t14*dw*Kb;
t62 = mp.sqrt(t61);
t63 = t62/2;
t64 = t4-t5-t6+t7+t8+t63;
t65 = -1*mp.mpc(0,1)/2;
t66 = t65*dw;
t69 = mp.exp((t66-t5-t6-t7-t8+t63)*tcp);
t71 = 1/t62;
t73 = t4-t5-t6+t7+t8-t63;
t76 = mp.exp((t66-t5-t6-t7-t8-t63)*tcp);
t79 = t64*t69*t71-t73*t76*t71;
t81 = mp.mpc(0,1);
t82 = t81*dw;
t85 = t42*(t82+Ka+Ra-Kb-Rb+t37)*t43;
t88 = t49*(t82+Ka+Ra-Kb-Rb-t37)*t43;
t89 = t85-t88;
t94 = t69*t71*Ka-t76*t71*Ka;
t96 = t52*t79+t89*t94/2;
t97 = t66-t5-t6+t7+t8+t38;
t98 = 1/Ka;
t101 = t66-t5-t6+t7+t8-t38;
t104 = t97*t98*t85-t101*t98*t88;
t105 = t96*t104/2;
t109 = t69*(t82-Ka-Ra+Kb+Rb-t62)*t71;
t114 = t76*(t82-Ka-Ra+Kb+Rb+t62)*t71;
t116 = -t64*t98*t109+t73*t98*t114;
t118 = -t109+t114;
t120 = t52*t116/2+t89*t118/4;
t121 = t120*t89/2;
t126 = t97*t42*t43-t101*t49*t43;
t129 = t126*t79+t104*t94/2;
t130 = t129*t126;
t133 = t126*t116/2+t104*t118/4;
t134 = t133*t52;
t135 = t105+t121;
t136 = t135**2;
t137 = t130+t134;
t140 = t137**2;
t146 = t96*t126+t120*t52;
t150 = mp.sqrt(t136-2*t137*t135+t140+4*(t129*t104/2+t133*t89/2)*t146);
t151 = t105+t121-t130-t134-t150;
t153 = N/2;
t154 = (t105/2+t121/2+t130/2+t134/2+t150/2)**t153;
t156 = 1/t150;
t158 = t105+t121-t130-t134+t150;
t160 = (t105/2+t121/2+t130/2+t134/2-t150/2)**t153;
t165 = 1/t146;
t177 =
1/(Ka+Kb)*((-t151*t154*t156/2+t158*t160*t156/2)*Kb+(-t151*t165*t154*t158*t156/2+t158*t165*t160*t151*t156/2)*Ka/2);
intensity0 = mp.mpf(str(pa)); # pA
intensity = t177.real;
R2eff.append( ((1/Tc)*mp.log(intensity0/intensity))); # we did
not factor out (Ka+Kb)/2 here
#print 'Precision of calculations mp.dps: ', mp.dps
#print 'R2eff: ', R2eff
# R2eff.append(2.)
array=[]
for i in range(len(ncyc)):
array.append((nu_cpmg[i],R2eff[i]))
if(outfile!='Null'):
outy=open(outfile,'w')
for i in range(len(array)):
outy.write('%f\t%f\n' % (array[i][0],array[i][1]))
outy.close()
return array
############################################################################
#numpy approximate formula
def Baldwin(kex,pb,dw,ncyc,Trel,R2g,R2e,outfile,verb='n'):
pa=(1-pb)
keg=kex*(1-pb)
kge=kex*pb
deltaR2=R2e-R2g
nu_cpmg=ncyc/Trel
tcp=Trel/(4.0*ncyc) #time for one free precession element
#########################################################################
#get the real and imaginary components of the exchange induced shift
g1=2*dw*(deltaR2+keg-kge) #same as carver richards zeta
g2=(deltaR2+keg-kge)**2.0+4*keg*kge-dw**2 #same as carver richards psi
g3=cos(0.5*atan2(g1,g2))*(g1**2.0+g2**2.0)**(1/4.0) #trig faster than
square roots
g4=sin(0.5*atan2(g1,g2))*(g1**2.0+g2**2.0)**(1/4.0) #trig faster than
square roots
#########################################################################
#time independent factors
N=complex(kge+g3-kge,g4) #N=oG+oE
NNc=(g3**2.+g4**2.)
f0=(dw**2.+g3**2.)/(NNc) #f0
f2=(dw**2.-g4**2.)/(NNc) #f2
#t1=(-dw+g4)*(complex(-dw,-g3))/(NNc) #t1
t2=(dw+g4)*(complex(dw,-g3))/(NNc) #t2
t1pt2=complex(2*dw**2.,-g1)/(NNc) #t1+t2
oGt2=complex((deltaR2+keg-kge-g3),(dw-g4))*t2 #-2*oG*t2
Rpre=(R2g+R2e+kex)/2.0 #-1/Trel*log(LpreDyn)
#do calc in numpy
E0= 2.0*tcp*g3 #derived from relaxation #E0=-2.0*tcp*(f00R-f11R)
E2= 2.0*tcp*g4 #derived from chemical shifts
#E2=complex(0,-2.0*tcp*(f00I-f11I))
E1=(complex(g3,-g4))*tcp #mixed term (complex) (E0-iE2)/2
ex0b=(f0*numpy.cosh(E0)-f2*numpy.cos(E2)) #real
ex0c=(f0*numpy.sinh(E0)-f2*numpy.sin(E2)*complex(0,1.)) #complex
ex1c=(numpy.sinh(E1)) #complex
v3=numpy.sqrt(ex0b**2.-1) #exact result for v2v3
y=numpy.power((ex0b-v3)/(ex0b+v3),ncyc)
v2pPdN=(( complex(deltaR2+kex,dw) )*ex0c+(-oGt2-kge*t1pt2)*2*ex1c)
#off diagonal common factor. sinh fuctions
Tog=(((1+y)/2+(1-y)/(2*v3)*(v2pPdN)/N))
Minty=Rpre-ncyc/(Trel)*numpy.arccosh((ex0b).real)-1/Trel*numpy.log((Tog.real))
#estimate R2eff
if(verb=='n'):
return
else:
array=[]
for i in range(len(ncyc)):
array.append((nu_cpmg[i],Minty[i]))
if(outfile!='Null'):
outy=open(outfile,'w')
for i in range(len(array)):
outy.write('%f\t%f\n' % (array[i][0],array[i][1]))
outy.close()
return array
#######################################################################
#Loop over a bunch of exchange parameters
def TestDisp(ncyc,Trel,R2g,flg,tim):
kexMax=1.0 #in s-1
kexMin=5000.0
dwMax=0.1 #in ppm
dwMin=1000.0
pbMax=0.004
pbMin=0.4
dRMin=1.
dRMax=1E4
Grid=11
now1=datetime.now()
rar=[]
for k in range(Grid):
for i in range(Grid):
for j in range(Grid):
#for l in range(Grid):
kex=kexMin*10**(+(1.0*i/(Grid-1.0))*log10(kexMax/kexMin))
pb= pbMin*10**(+(1.0*j/(Grid-1.0))*log10(pbMax/pbMin))
dw= dwMin*10**(+(1.0*k/(Grid-1.0))*log10(dwMax/dwMin))
#R2e= dRMin*10**(+(1.0*l/(Grid-1.0))*log10(dRMax/dRMin))
R2e=R2g
if(flg=='Carver'):
arr=CarverRichards(kex,pb,ppm_to_rads(dw,dfrq),ncyc,Trel,R2g,R2e,'Null',verb='y')
if(flg=='Baldwin'):
arr=Baldwin(kex,pb,ppm_to_rads(dw,dfrq),ncyc,Trel,R2g,R2e,'Null',verb='y')
if(flg=='Numerical'):
arr=NumDisp(kex,pb,ppm_to_rads(dw,dfrq),ncyc,Trel,R2g,R2e,'Null',verb='y')
if(flg=='Nikolai'):
arr=Nikolai(kex,pb,ppm_to_rads(dw,dfrq),ncyc,Trel,R2g,R2e,'Null',verb='y')
if(flg=='NikolaiLong'):
arr=NikolaiLong(kex,pb,ppm_to_rads(dw,dfrq),ncyc,Trel,R2g,R2e,'Null',verb='y')
if(flg=='NikolaiFun'):
arr=NikolaiFun(kex,pb,ppm_to_rads(dw,dfrq),ncyc,Trel,R2g,R2e,'Null',verb='y')
if(flg=='Meiboom'):
arr=Meiboom(kex,pb,ppm_to_rads(dw,dfrq),ncyc,Trel,R2g,R2e,'Null',verb='y')
rar.append(arr)
now2=datetime.now()
val=timedelta.total_seconds(now2-now1)
if(flg=='Baldwin'):
tim=timedelta.total_seconds(now2-now1)
print flg,val,val/tim
return val,numpy.array(rar)
def ppm_to_rads(ppm,dfrq):
return ppm*2*numpy.pi*dfrq
if __name__ == "__main__": #run this if this file is run as standalone
print 'Relative times of the various approaches:'
dfrq=200. #spectrometer frequency of nucleci (MHz)
ncyc=numpy.array((2,4,8,10,12,14,16,18,20,22,24,26,28,30,32,34,36,38,40))
#number of cpmg cycles
R2g=10. #relaxation rate of ground
R2e=10. #relaxation rate of excited
Trelax=0.02 #total time of CPMG block
tim,barr=TestDisp(ncyc,Trelax,R2g,'Baldwin',0.)
tom,narr=TestDisp(ncyc,Trelax,R2g,'Numerical',tim)
tom,marr=TestDisp(ncyc,Trelax,R2g,'Meiboom',tim)
tom,carr=TestDisp(ncyc,Trelax,R2g,'Carver',tim)
tom,sarr=TestDisp(ncyc,Trelax,R2g,'Nikolai',tim)
tom,garr=TestDisp(ncyc,Trelax,R2g,'NikolaiFun',tim)
mp.mp.dps=18
tom,larr=TestDisp(ncyc,Trelax,R2g,'NikolaiLong',tim)
mp.mp.dps=23
tom,jarr=TestDisp(ncyc,Trelax,R2g,'NikolaiLong',tim)
print 'Maximum error over grid (s-1):'
print 'Meiboom: ',numpy.max(numpy.abs(narr[:,:,1]-marr[:,:,1]))
print 'Baldwin: ',numpy.max(numpy.abs(narr[:,:,1]-barr[:,:,1]))
print 'CarverRichards: ',numpy.max(numpy.abs(narr[:,:,1]-carr[:,:,1]))
#remove nan's from Nikolai's double output.
test=numpy.abs(narr[:,:,1]-garr[:,:,1])
for i in range(len(test)):
for j in range(len(test[i])):
if(test[i][j]!=test[i][j]):
test[i][j]=0.
print 'Nikolai fun: ',numpy.max(test)
#remove nan's from Nikolai's double output.
test=numpy.abs(narr[:,:,1]-sarr[:,:,1])
for i in range(len(test)):
for j in range(len(test[i])):
if(test[i][j]!=test[i][j]):
test[i][j]=0.
print 'Nikolai dougle (9): ',numpy.max(test)
print 'Nikolai long double (18):
',numpy.max(numpy.abs(narr[:,:,1]-larr[:,:,1]))
print 'Nikolai long double (23):
',numpy.max(numpy.abs(narr[:,:,1]-jarr[:,:,1]))
"""
print
print 'Short test for exactness'
kex=1.
pb=0.1
dw=1000
ncyc=numpy.array((2,4,6,8,10,12,14,16,20,24,28,32,36,40)) #number of cpmg
cycles
Narr=NumDisp(kex,pb,ppm_to_rads(dw,dfrq),ncyc,Trelax,R2g,R2e,'Null',verb='y')
Carr=CarverRichards(kex,pb,ppm_to_rads(dw,dfrq),ncyc,Trelax,R2g,R2e,'Null',verb='y')
Aarr=Baldwin(kex,pb,ppm_to_rads(dw,dfrq),ncyc,Trelax,R2g,R2e,'Null',verb='y')
Warr=Nikolai(kex,pb,ppm_to_rads(dw,dfrq),ncyc,Trelax,R2g,R2e,'Null',verb='y')
Marr=Meiboom(kex,pb,ppm_to_rads(dw,dfrq),ncyc,Trelax,R2g,R2e,'Null')
larr=NikolaiLong(kex,pb,ppm_to_rads(dw,dfrq),ncyc,Trelax,R2g,R2e,'Null')
col=len(ncyc)-1 #comparing ncyc=40
print 'R2eff for ncyc=40:'
print 'Numerical result :',Narr[col][1]
print 'Nikolai result : ',Warr[col][1]
print 'CarverRichards : ',Carr[col][1]
print 'Baldwin result : ',Aarr[col][1]
print
print 'Errors for ncyc=40, R2(calc)-R2(numerical)'
print 'Nikolai Error : ',Warr[col][1]-Narr[col][1]
print 'Baldwin Error : ',Aarr[col][1]-Narr[col][1]
print 'CarverRichards err: ',Carr[col][1]-Narr[col][1]
print 'To compare directly with matlab code:'
print 'Python implementation of Nikolais formula gives (nu_cpmg,R2eff):'
for i in range(len(Warr)):
print Warr[i]
print Marr
for i in range(len(Aarr)):
print larr[i]
"""
Re: [sr #3154] Implementation of Baldwin (2014) B14 model - 2-site exact solution model for all time scales