Andy, you need to be a bit more fair here! If Troels can come up with a proof of the Carver and Richards 1972 equations reducing to the Luz and Meiboom 1963 equations, I would expect at least a few pints ;)
Regards, Edward On 6 May 2014 09:25, Andrew Baldwin <[email protected]> wrote: > Hiya Troels, > > Very commendable tenacity with the fluorescence result :) Good stuff indeed! > > 1) Definition of alpha- : > > That's a fair point. I agree. Alpha- is later squared so either: > > zeta(B14) = 2.0*dw (Delta R2 + k_eg - k_ge) > zeta(B14) = 2.0*dw (-Delta R2 - k_eg + k_ge) > > Will work equally well. Which you use depends on how you choose to > parametrise the Eigenvalues. h1 can be plus or minus depending on your mood. > So my mood says keeping the positive deltaR2 makes the equation prettier. > > 2) You might need to give a bit more thought to the numerical > implementation. Notice that I add the factor of R2g back in, during the > function Numdisp. Turns out you can factor anything out from the diagonal as > long as you put it back in at the end, hence the form of the matrix in the > code I sent. > > You can prove this by looking at the Eigenvalues, and how they propagate > when you take things like matrix exponentials. In the same vein, notice that > In the paper I subtract R2G from both f00 and f11. I add it back in again > later in the constants out the front. This is following from the lead of > Nikolai in his reconstructing invisible states paper. Implicitly looking at > deltaR2s shows clearly where differences in relaxation have effects, so this > brings out the physics. > > So look at the 'error' reported in the outputs. When you made the Troels > version of the numerical code, the error for the Baldwin formula was 10s-1. > This is actually the exact value you set R2G to. So while you added in R2G > in the evolution matrix (which is not wrong to do), you didn't then remove > it in numdisp, so your solution is out by one unit of R2G. > > Looking at the outputs, the code I sent for CR had the carver richard's > formula set in the wrong 'minus' configuration, added during testing while > writing the previous email. This was certainly my foobar. However, note that > the maximum 'error' reported by the code was 4989 s-1 when you ran the code > for the first time. In my previous email: > > >> with the alpha- sign made consistent with the Carver Richard's paper gives >> an error of 4989 s-1 in the same test > > > That should have set off alarm bells for you. Setting either: > > zeta=2*dw*(deltaR2+keg-kge) #zeta=g1 > zeta=2*dw*(-deltaR2-keg+kge) #zeta=g1 > > In the CR function as you note (and as you do in relax) is good. If you > subtract the 'Baldwin error' from the 'CR error' you get 8. 8 s-1 is the > biggest error incurred using CR in this parameter range (as I stated in the > last email). So it's fair to call the CR equation an approximation. Fig 1 > speaketh the truth. Calculate it yourself first using your code before > wielding the accusing stick! > > So I'm deducting one half pint from your score for not checking the results > from your modified numerical code more carefully. > > 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. > > > Pint for the typo in keg/kge. That's a genuine foobar in the manu - the > equations should match code i've written, and I know the code is right > because i've matched it to simulation. I'll keep my Eigenvalue definitions > though. > > > 2) zeta typo in eq 70. z=zeta > > Agreed! +1. > > > 3) Watch alpha_ in eq 70. Get Delta R straight > > See 1. No extra pint here. > > 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. > > See 1. No extra pint here either. > > 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??? > > > -1/2 pint for an over-reached conclusion from your modified code. > > So current score = 1.5 > > > Meiboom challenge: I'd be interested in a proof for how Carver Richard's > reduces to Meiboom. I'd pay a pint for that. To get you started, note that > in the fast exchange limit, E0 and E2 can be simplified using the fast > exchange limits in supp section 1. > > Best, > > Andy. > > > > On 05/05/2014 15:27, Troels Emtekær Linnet 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
