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 <[email protected]> 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 <[email protected]>: >>> 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 <[email protected]>: >>>> >>>> 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 <[email protected]>: >>>> >>>> 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. >>>> >>>> >>>> >>>> >> _______________________________________________ relax (http://www.nmr-relax.com) This is the relax-devel mailing list [email protected] To unsubscribe from this list, get a password reminder, or change your subscription options, visit the list information page at https://mail.gna.org/listinfo/relax-devel
