On Thu, Jan 8, 2009 at 5:40 AM, Sébastien Morin <sebastien.morin.1@xxxxxxxxx> wrote:
Hi, I thought about the task (Should I add the implementation of relaxation dispersion as a task on the Gna web site ?) and here is how I see the flow-through of an analysis...
If you would like to add it as a task, feel free and then assign it to yourself.
1. A user will want to perform relaxation dispersion analysis. The user will first have to chose whether CPMG or R1rho experiments were recorded, the choice of which will determine if intensities are associated to different CPMG pulse train frequencies or to different R1rho spin lock field strength. Let's say the user types the command: relax_disp.exp_type('cpmg') and chooses the CPMG experiment. The intensities will then be associated to different CPMG pulse train frequencies, or to a null frequency (for the reference spectrum).
Having a value of None or something like a string value of 'ref' for the CPMG frequencies would be enough to identify the reference spectrum.
We can now assume that the user recorded CPMG experiment from which R2eff values are extracted as follows: R2eff = ( 1 / T ) * Ln( Icpmg / Iref ) where T is the constant time relaxation delay, Icpmg is the peak height with CPMG, and Iref is the peak height without CPMG (reference). There should be a function to input the delay T used... Should this function be called something like: relax_disp.cpmg_delayT() and accept a float ?
This should accept a float, and it should be a value in seconds. Other suggestions for the name of the user function: relax_disp.tau_cpmg(), relax_disp.CT_period(), relax_disp.CT_length(), relax_disp.cmpg_delay(), etc. I'm not sure of the equation yet, but If we have to convert to the unitless radial units this can be done later within relax.
From this step, peak intensities will be used to calculate R2eff for the different frequencies used. How should duplicated reference spectra (if present) be treated (since they have an effect in every R2eff calculated) ?
I think they should be treated exactly the same as in the exponential relaxation curve-fitting - i.e. treat each duplicate spectrum as a separate time point. The exception would be the reference spectrum where the average value must, obviously because of the equation you give above, be used. I think a better error analysis for relaxation dispersion would not be replicated spectra but rather measuring the RMSD of the base plane noise and using that as the errors for the minimisation and Monte Carlo simulations. Does anyone know how error analysis is performed in relaxation dispersion? I'm not familiar with it enough.
In the future, we could add a function for selecting this approach or another (if others exist and are of interest).
All this infrastructure is already within relax, and is independent from the analysis type. I developed it all in the 'spectral_errors' branch which was merged back into the mainline at r8093. The only code missing would be for handling duplicated reference spectra, but this should be code specific to relaxation dispersion. And it could be developed later on as I don't know if people collect the reference spectrum in duplicate.
2. The user will then need to chose which mathematical model to use, fast or slow exchange. Let's say the user types a command like: relax_disp.select_model('fast') and chooses fast exchange. (Should this function be renamed to something like: relax_disp.cpmg_select_model() since it is associated to CPMG experiments ?)
For simplicity, I think that this function should handle both the CPMG and the R1rho and simply raise a RelaxError if 'r1rho' is set and a CMPG equation is selected (and probably raise a RelaxError if no experiment type is set). Are there any equations shared between the CMPG and R1rho analyses? If each model identifies the experiment type unambiguously, then the relax_disp.exp_type() user function is not necessary. I'm not sure what is best here until I see all the equations - even those for the accordion type experiments. I saw the code you added here, and I think that putting the equations, references, and descriptions into the user function docstrings as you did in your Mathematica notebook would be very useful.
The R2eff values calculated in part 1 will then be input in the chosen equation and minimised. This is different as in standard curve fitting for R1 and R2, since the R2eff values are first calculated and then are passed to the dispersion equations to be minimised. Am I right ?
This is correct. There should be a method added to the specific_fns.relax_disp.Relax_disp class which calculates R2eff using the equation you gave above. This method can then be called by methods such as calc(), back_calc(), grid_search(), minimise(), etc. to generate the R2eff values.
In the case of CPMG data in the fast-exchange limit, should the relaxation dispersion specific user commands be as follows: relax_disp.exp_type('cpmg') relax_disp.cpmg_delayT(0.040) relax_disp.cpmg_select_model('fast') and should the data go as follows: intensities --> calculated [R2eff] --> minimised [R2, Rex, kex]
Here in the script you would use: spectrum.read_intensities() spectrum.replicated() or spectrum.baseplane_rmsd() spectrum.error_analysis() frq.set() (set the spectrometer frequency) relax_disp.nu_cmpg() (or maybe tcp_cmpg in seconds, and then relax will convert to radial units) grid_search(inc=11) minimise('simplex', scaling=False, constraints=False) The grid_search() and minimise() functions should automatically at the start calculate R2eff.
====================================== ====================================== Moreover, what about the intermediate-exchange limit ? Is that of any importance and should it be implemented also in relax ? I propose yes, although I'll have to check for the equations, and everything...
Again this can be treated as another model - which can be added later. If the infrastructure for the other models is in place, then adding this will be insanely trivial.
Also, as for choosing whether slow-, intermediate- or fast-exchange equations are best suited for the data, there is, of course, the use of model selection implemented into relax, but also the use of the alpha parameter (Millet et al., JACS, 2000, 122: 2867-2877 (equation 14)) which can help determining the exchange limit... How could we combine both ?
If alpha is found within an equation and is optimised - then this is again another model. This model tells you where in the exchange regime you are. This can go directly into AIC and BIC model selection for comparison to the other models such as 'fast', 'slow', etc. When ever there is a different equation, even slightly, then this is a new mathematical model.
Finally, what about group fitting of dispersion data to extract better parameters based on a common exchange rate (kex or kA, for fast- or slow-exchange, respectively) ? Should this be thought about in the beginning of the development or will it be easy to implement this when the rest works ?
Do you mean have one parameter that spans multiple spins? This again would be another model. It would be global though as it affects more than one spin. It will require a different set of equations - ones that loop over the relevant spins. It will be a bit like the global model-free models where the diffusion tensor is optimised affecting the model-free parameters of individual spins. This, like all new models, can be added later with zero redesign of relax. The infrastructure will be in place and you just have to code the equations into a new function in maths_fns (and add some documentation to the relax_disp.select_model() user function). Cheers, Edward