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
Re: Relaxation dispersion