firasm / cest Goto Github PK

The simulation now runs through the saturation pulse, allows the system to relax for some time defined by adjustable sequence parameters with a hard flip here and there, and uses the resulting magnetization state as the initial state for the next saturation frequency.

Here I simulated alternating frequencies to try and replicate the double-peak phenomenon we've observed in phantoms. There is only one simulated species to the left of the water peak; the peak on the right is an artifact of residual saturation from the previous frequency. It dissappears if the inter-frequency time is long enough, about 5000 - 70000 ms (in excellent agreement with phantom studies... almost too excellent).

Also note the shoulders on either side of the spectrum as the system approaches a steady state.

(Forgot to say that the legend gives the interfreq delay for each spectrum)

Some parameters:
B0 = 7.0 T (easy upgrade!)
B1 = 1.0e-6
saturation pulse duration = 500 ms
there is one 20 degree flip in the sequence

M0a = 1000.0
M0b = M0ax1e-3x10 (so 10mM, used a high concentration for easily cest effect)
T1a = 2.0 # seconds
T2a = 0.6
T1b = 1.0
T2b = 0.2
kb = 40.
ka = M0b/M0a*kb #Transfer rate of pool a to pool b, s^-1. We want a-->b = b-->a
accFactor = 1 (more on this later)

While you are on the subject of advancing the simulation capabilities, it could be useful to look at saturation schemes where the CW sat pulse is replace by a series (n=100-1000?) of short (ms?), gaussian (or start with block if easier?), separated by short gaps (ms? -> duty cycle 10-100%).

The CERT (chemical exchange rotation transfer or some such) is based around the idea of comparing two spectra with differences in only their sat pulse parameters thereby replacing the questionable Z_asym analysis with something slightly better founded.

Did a quick experiment today with various acceleration factors and RF Spoiling on and off.

I think I have some SNR issues at high acceleration

Some conclusions:

- Cannot reproduce the second peak, with or without RF spoiling
- Noise increases with increased acceleration
- RF spoling reduces noise in the spectrum
- Peak height decreases with increased acceleration
- Line width increases with increased acceleration
- Where the hell is that second peak when you scan without alternating frequencies !?

Hi all,

Here are the CEST Asymmetry maps from today's mouse brain experiment.

They're calculated as follows:

And here's roughly the SNR to expect from a tumour (Source: KD's thesis):

I think we're doing okay!

@DrSAR mentioned to Michael McMahon the oscillation issue that we also see in his simulation.

He knew about it and suggested that it's because the saturation pulse duration isn't long enough and the "phase gets screwy".

@srveale, when you have a chance, can you switch the integrator back to something fast, make sure you can still see the oscillations, and then run another iteration JUST changing the pulse length to 4seconds.

This fixed this in his simulation and it should for us as well.

Firas

This is a report of the two tumour bearing mice I imaged on April 16, 2015.

Tumour type: BT474
Tumour location: mammary fat pad
T1 map:

• The first CEST scan I tried was the one that we didn't get to do with the previous tumour mouse because it died.

Relevant parameters:

Matrix Size: 64x64
Spatial Resolution: 0.469 mm x 0.469 mm x 3 mm
Acc. factor: 8
InterSatFreqDelay: 0ms
Pulse Length: 1000 ms

Lesson learned: Wait time between saturation frequencies matters. First I thought that the only difference would be the mirror peak because that's what we saw in phantoms. But the mirror peak is just one manifestation of what happens if you don't wait. In Vivo, the manifestation of this effect is different.

• I repeated the scan above with InterSatFreqDelay time to be 7000ms

We clearly see one big peak at ~ -3 ppm (aliphatic) but nothing on the other side, where we expect two peaks from the amide and amine protons.

• In an effort to try and eliminate SNR from the picture, the spatial resolution was reduced by cutting in half the matrix size.

Things seemed to have gotten a bit better, but still no peaks on the other side (data not shown).

At this point the mouse had been in the scanner for over 5 hours so it was euthanized and the next mouse put in the scanner.

• The first scan for tumour 3 was the longest possible scan we could acquire - no acceleration, long InterSatFreqDelay (7000 ms), maximum SNR (1mm x 1mm x 5mm).

I also went one step further and first water shifted all the voxels, and then averaged the spectrum to make absolutely sure this wasn't an SNR issue:

Now we see something on the other side, but it’s not clear if there is an “approach to steady state” effect here.

Conclusion: Need to do some more reading of Kim's thesis to understand what's going on here.

Here's the plot generated by the simulation. Some parameters:
B0 = 3.0 Tesla
Species chemical shift = 3000 Hz (~3.7 ppm)

Pool a is water, b is the species
T1a = 0.6 s
T2a = 0.1 s
T2b = 0.33 s
T1b = 0.770 s
ka = 1 (Transfer rate of pool a to pool b, s^-1)
kb = 200.

I'll try to walk through this here, but the notebook CEST/simulation/Acc Factor might be good to look at.

First I plotted each component of magnetization, starting just after saturation finished. It was an off-resonance pulse so the magnetization was still close to equilibrium. Then the flip of 20 degrees around the y direction is applied, putting some of the z -magnetization into the x-direction. Free precession is allowed for 10 ms, then another flip, repeating like this for accfac = 72 times, producing a jagged exponential decay type curve. When Mz is small enough, T1 relaxation is fast enough to go back to the same level after each flip, and a steady state is reached.

The steady state should depend on T1 and the length of time between flips. It should not depend on the state of the magnetization just after saturation.

Next is a plot showing several z-spectra. The topmost curve corresponds to the first acquisition after the saturation, and the order continues downwards. CEST contrast is the highest for the first curve, as expected.

Different weighting schemes could be used, depending on T1 and the time in between flips. Or, we could use the steady-state saturation to define T1.... Or we could use a known T1 to determine the best times between sequence pulses.

housed here:
http://godzilla.kennedykrieger.org/

The data processing tools might also be useful but are written in MatLab. Might work in octave (the open-source pendant to this)

For DCE scans we have a VTC (voxel-based time curve) that would be handy when looking at Z Spectra.

Here's the first crack at it using the second mouse tumour from this week.

Any suggestions for modifications and enhancements ?

Yesterday @andrewcyung emailed an image of the first successful scan that exhibited CEST spectra. We've been trying to reproduce this effect size in all of our subsequent imaging sessions.

I decided to have a closer look at this scan, and it turns out that it isn't as great as we thought.

Firstly, this effect isn't actually representative. For instance, this example just two pixels away doesn't show 3 peaks nearly as well as that one.

Second, here's a fit of the data ; as you can see, peak size is very comparable to that in issue #46 with the EPI scan.

Third, that peak at -3.5 ppm is very strange ...

Fourth, our shim on this day wasn't that good and the "good" pixel is right in the middle of the bad shim region. Here is a map of the water shift (in ppm)

This map was obtained by fitting just the water peak to the CEST spectrum, and plotting the ppm shift (from 0) spatially.

The fits aren't great, but the water peak fits reasonably well.

Stefan asked me to prepare this for the group meeting.

Animal: scn = sarpy.Scan('HPCoil1.tQ2/14')

This might actually be due to a bug in atleast_4d (the way I have programmed it...) Why would you not get this error?

In [1]: run NecS3_LL_roi_flattened.py
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
/usr/local/lib/python2.7/dist-packages/IPython/utils/py3compat.pyc in execfile(fname, *where)
    176             else:
    177                 filename = fname
--> 178             __builtin__.execfile(filename, *where)

/home/stefan/firas-analysis/NecS3/NecS3_LL_roi_flattened.py in <module>()
     49         roi2 = sarpy.Scan(master_sheet[k][key_list[5]][0]).adata[key_list[7]]
     50 
---> 51         roi_resample1 = sarpy.ImageProcessing.resample_onto.resample_onto_pdata(roi1,LL_1)
     52         roi_resample2 = sarpy.ImageProcessing.resample_onto.resample_onto_pdata(roi2,LL_2)
     53 

/home/stefan/sarpy/sarpy/ImageProcessing/resample_onto.pyc in resample_onto_pdata(source_pdata, target_pdata)
    157     # loop over all non-3D dimensions:
    158     # We have to flatten all higher dimensions and then go through them.
--> 159     flat_input_img = atleast_4d(atleast_4d(source_pdata.data)[:,:,:,...])
    160     dims = source_pdata.data.shape
    161     extra_dims = int(numpy.prod(dims[3:]))

/home/stefan/sarpy/sarpy/ImageProcessing/resample_onto.pyc in atleast_4d(arr)
     27         (120, 1, 1, 1)
     28     '''
---> 29     arr.shape += (1,) * (4 - arr.ndim)
     30     return arr
     31 

/usr/lib/python2.7/dist-packages/numpy/core/memmap.pyc in __array_finalize__(self, obj)
    255         if hasattr(obj, '_mmap'):
    256             self._mmap = obj._mmap

--> 257             self.filename = obj.filename
    258             self.offset = obj.offset
    259             self.mode = obj.mode

AttributeError: 'memmap' object has no attribute 'filename'

Seems like it is possible for a smaller peak to be hidden by its neighbours.

From Andrew's email:

This is from a single pixel in the cortex. At first glance, the data looks much better than before, but I am also using a low resolution comparable to what Kim was using. I repeated with a variety of powers. Looks like 1 uT gave us the best results on that day.

Tried to fit this with multiple lorentzians

To start off, here is what the signal history looks like for one frequency offset (0.2ppm), for acceleration factor of 8:

And for acceleration factor of 24:

For brevity I won't post all of them, but I'll go straight to the point. The z-spectra at the centre of k-space:

Yes, the order is out of whack in places. AccFactor= 72 is above accFactor = 1 in signal magnitude, but has a smaller peak. My suspicion is that this is related to what was happening before with shorter saturation durations - notice the small scale oscillations. This was done using saturation duration = 500ms, and I'll try with 4000ms, but these take a long time because each frequency offset has as many as 72 saturation pulses.

Let me know if there is anything I should be looking at differently.

Here's an arginine NMR spectrum. It doesn't say what the reference is so I'm assuming it's TMS.

Here's a 3D viewer for the molecule:

http://bmrb.wisc.edu/metabolomics/mol_summary/jmol_display.php?path=/metabolomics/standards/L_arginine/lit/3362.mol&mol=L-Arginine&bmrbid=bmse000029#

And here's which atoms they assigned to which peak (there's an option in the 3D viewer to show atom ID):

"Assigned Chemical shifts for bmse000029, L-Arginine"

Shift ID Atom ID Chemical shift Ambiguity
1 C2 57.002 1
2 C3 177.238 1
3 C5 30.281 1
4 C6 26.577 1
5 C7 43.201 1
6 C9 159.504 1
7 H15 3.764 1
8 H16 1.909 1
9 H17 1.909 1
10 H18 1.679 1
11 H19 1.679 1
12 H20 3.236 1
13 H21 3.236 1

So my guess is we are interested in the main chain hydrogens, ie H16 and H17, and the second peak that appears at low pH presumably corresponds to H18 and H19. Is this the right direction?

By increasing our InterSatFreq Delay time, it seems that we may as well vastly increase our Spatial Resolution without much of a time cost (there is a ghosting artifact that is introduced)!! Of course, this makes sense in retrospect, but who knew!

P.S. This data was acquired at an Acceleration Factor of 32.

Here are the images:

And the CEST spectra:

Ignore the empty plots, that data was not acquired.

To address the issue of the mirror peak (#18), currently we need to increase the InterSaturationDelay time to be much more than the T1 of the sample.

This is very inefficient because we end up spending AccFactor time as much as we need to for relaxation.

I seem to remember us thinking about this already and there is a delay time present, setting it is just buried deep in the method editor.

Results of an experiment I did last night with a variety of Flip angles (5 degrees - 25 degrees), two powers (0.5 uT and 1uT) and two acceleration factors (2 and 4).

The gold standard is the ubcCEST Flash sequence, in green.

Summary: you can go pretty damn fast, but there's an SNR penalty, and you must pay it and be above the threshold otherwise your data will be crap (but it will be acquired very fast).

For some inexplicable reason, a new peak pops up at -3 ppm if you image too fast.

Overall, good news...

Would be handy to check on status of scans remotely...

Now this is a cool result. At pH 5, we now see two peaks! Another interesting result is how little the spectrum changes due to acceleration factors

Scott and I (mostly Scott) has spent most of the night scanning and here is the result of a particularly interesting experiment.

We kept all parameters the same, and just changed the Acceleration Factor from 1 all the way to 32. The scan speed-up and the Acceleration Factors are shown in the legend.

Updated to have thinner lines so you can see more differences

More results from 'CEST_SV.uQ1/' experiment that Scott conducted

For a T1 of around 500ms, the question was how long does our InterSatFreq Delay (wait time between two successive offset freq. acquisitions).

Increasing InterSatFreq Delay reduces the impact of the mirror peak, but it would be nice to sort this out. It seems related to T1, but it can't just be T1 related because 7500ms is 15 times our T1

Now that we have working phantoms, we finally did the experiment to see how the InterSaturation delay time affects the peaks.

The hypothesis was that the second peak (at -3 ppm) was appearing because of insufficient relaxation between two successive saturation pulses.

Here is an image of the phantom used for this experiment, the left phantom has a high Gd concentration (5x the right one) and the right one has a low Gd concentration. The two stars are the pixels used for the plot below.

The left phantom has a T1 of ~ 527 ms and the right phantom has a T1 of ~ 1441 ms.

At acceleration factor 8, we acquired data at different Intersaturation delay times and the plot is here:

Conclusion: We need to wait a sufficiently long time (5xT1 ?) in order to eliminate the second peak

This is the output I get:

In [1]: run NecS3_LL_roi_flattened.py
Expected error, please ignore NecS3Hs02
Expected error, please ignore NecS3Hs03
Expected error, please ignore NecS3Hs13
Expected error, please ignore NecS3Hs12
Expected error, please ignore NecS3Hs06
Expected error, please ignore NecS3Hs07
Expected error, please ignore NecS3Hs04
Expected error, please ignore NecS3Hs05
Expected error, please ignore NecS3Hs08
Expected error, please ignore NecS3Hs09
Expected error, please ignore NecS3Hs01a
Expected error, please ignore NecS3Hs15
Expected error, please ignore NecS3Hs10
Expected error, please ignore NecS3Hs11
Expected error, please ignore NecS3Hs14

You are catching all exceptions in your promiscuous try/except block and I don't know why it's failing... In the end the relevant array come up empty:

In [2]: T1_vals_day1.shape
Out[2]: (0,)
In [4]: T1_vals_day2.shape
Out[4]: (0,)

which is a little disappointing.

Here's a handy list of all the exceptions, prune this list with the most common ones and float them to the top.

Source: http://docs.python.org/2/library/exceptions.html

The following exceptions are the exceptions that are actually raised.

exception AssertionError
Raised when an assert statement fails.

exception AttributeError
Raised when an attribute reference (see Attribute references) or assignment fails. (When an object does not support attribute references or attribute assignments at all, TypeError is raised.)

exception EOFError
Raised when one of the built-in functions (input() or raw_input()) hits an end-of-file condition (EOF) without reading any data. (N.B.: the file.read() and file.readline() methods return an empty string when they hit EOF.)

exception FloatingPointError
Raised when a floating point operation fails. This exception is always defined, but can only be raised when Python is configured with the --with-fpectl option, or the WANT_SIGFPE_HANDLER symbol is defined in the pyconfig.h file.

exception GeneratorExit
Raise when a generator‘s close() method is called. It directly inherits from BaseException instead of StandardError since it is technically not an error.

New in version 2.5.

Changed in version 2.6: Changed to inherit from BaseException.

exception IOError
Raised when an I/O operation (such as a print statement, the built-in open() function or a method of a file object) fails for an I/O-related reason, e.g., “file not found” or “disk full”.

This class is derived from EnvironmentError. See the discussion above for more information on exception instance attributes.

Changed in version 2.6: Changed socket.error to use this as a base class.

exception ImportError
Raised when an import statement fails to find the module definition or when a from ... import fails to find a name that is to be imported.

exception IndexError
Raised when a sequence subscript is out of range. (Slice indices are silently truncated to fall in the allowed range; if an index is not a plain integer, TypeError is raised.)

exception KeyError
Raised when a mapping (dictionary) key is not found in the set of existing keys.

exception KeyboardInterrupt
Raised when the user hits the interrupt key (normally Control-C or Delete). During execution, a check for interrupts is made regularly. Interrupts typed when a built-in function input() or raw_input() is waiting for input also raise this exception. The exception inherits from BaseException so as to not be accidentally caught by code that catches Exception and thus prevent the interpreter from exiting.

Changed in version 2.5: Changed to inherit from BaseException.

exception MemoryError
Raised when an operation runs out of memory but the situation may still be rescued (by deleting some objects). The associated value is a string indicating what kind of (internal) operation ran out of memory. Note that because of the underlying memory management architecture (C’s malloc() function), the interpreter may not always be able to completely recover from this situation; it nevertheless raises an exception so that a stack traceback can be printed, in case a run-away program was the cause.

exception NameError
Raised when a local or global name is not found. This applies only to unqualified names. The associated value is an error message that includes the name that could not be found.

exception NotImplementedError
This exception is derived from RuntimeError. In user defined base classes, abstract methods should raise this exception when they require derived classes to override the method.

New in version 1.5.2.

exception OSError
This exception is derived from EnvironmentError. It is raised when a function returns a system-related error (not for illegal argument types or other incidental errors). The errno attribute is a numeric error code from errno, and the strerror attribute is the corresponding string, as would be printed by the C function perror(). See the module errno, which contains names for the error codes defined by the underlying operating system.

For exceptions that involve a file system path (such as chdir() or unlink()), the exception instance will contain a third attribute, filename, which is the file name passed to the function.

New in version 1.5.2.

exception OverflowError
Raised when the result of an arithmetic operation is too large to be represented. This cannot occur for long integers (which would rather raise MemoryError than give up) and for most operations with plain integers, which return a long integer instead. Because of the lack of standardization of floating point exception handling in C, most floating point operations also aren’t checked.

exception ReferenceError
This exception is raised when a weak reference proxy, created by the weakref.proxy() function, is used to access an attribute of the referent after it has been garbage collected. For more information on weak references, see the weakref module.

New in version 2.2: Previously known as the weakref.ReferenceError exception.

exception RuntimeError
Raised when an error is detected that doesn’t fall in any of the other categories. The associated value is a string indicating what precisely went wrong. (This exception is mostly a relic from a previous version of the interpreter; it is not used very much any more.)

exception StopIteration
Raised by an iterator‘s next() method to signal that there are no further values. This is derived from Exception rather than StandardError, since this is not considered an error in its normal application.

New in version 2.2.

exception SyntaxError
Raised when the parser encounters a syntax error. This may occur in an import statement, in an exec statement, in a call to the built-in function eval() or input(), or when reading the initial script or standard input (also interactively).

Instances of this class have attributes filename, lineno, offset and text for easier access to the details. str() of the exception instance returns only the message.

exception IndentationError
Base class for syntax errors related to incorrect indentation. This is a subclass of SyntaxError.

exception TabError
Raised when indentation contains an inconsistent use of tabs and spaces. This is a subclass of IndentationError.

exception SystemError
Raised when the interpreter finds an internal error, but the situation does not look so serious to cause it to abandon all hope. The associated value is a string indicating what went wrong (in low-level terms).

You should report this to the author or maintainer of your Python interpreter. Be sure to report the version of the Python interpreter (sys.version; it is also printed at the start of an interactive Python session), the exact error message (the exception’s associated value) and if possible the source of the program that triggered the error.

exception SystemExit
This exception is raised by the sys.exit() function. When it is not handled, the Python interpreter exits; no stack traceback is printed. If the associated value is a plain integer, it specifies the system exit status (passed to C’s exit() function); if it is None, the exit status is zero; if it has another type (such as a string), the object’s value is printed and the exit status is one.

Instances have an attribute code which is set to the proposed exit status or error message (defaulting to None). Also, this exception derives directly from BaseException and not StandardError, since it is not technically an error.

A call to sys.exit() is translated into an exception so that clean-up handlers (finally clauses of try statements) can be executed, and so that a debugger can execute a script without running the risk of losing control. The os._exit() function can be used if it is absolutely positively necessary to exit immediately (for example, in the child process after a call to fork()).

The exception inherits from BaseException instead of StandardError or Exception so that it is not accidentally caught by code that catches Exception. This allows the exception to properly propagate up and cause the interpreter to exit.

Changed in version 2.5: Changed to inherit from BaseException.

exception TypeError
Raised when an operation or function is applied to an object of inappropriate type. The associated value is a string giving details about the type mismatch.

exception UnboundLocalError
Raised when a reference is made to a local variable in a function or method, but no value has been bound to that variable. This is a subclass of NameError.

New in version 2.0.

exception UnicodeError
Raised when a Unicode-related encoding or decoding error occurs. It is a subclass of ValueError.

UnicodeError has attributes that describe the encoding or decoding error. For example, err.object[err.start:err.end] gives the particular invalid input that the codec failed on.

encoding
The name of the encoding that raised the error.

reason
A string describing the specific codec error.

object
The object the codec was attempting to encode or decode.

start
The first index of invalid data in object.

end
The index after the last invalid data in object.

New in version 2.0.

exception UnicodeEncodeError
Raised when a Unicode-related error occurs during encoding. It is a subclass of UnicodeError.

New in version 2.3.

exception UnicodeDecodeError
Raised when a Unicode-related error occurs during decoding. It is a subclass of UnicodeError.

New in version 2.3.

exception UnicodeTranslateError
Raised when a Unicode-related error occurs during translating. It is a subclass of UnicodeError.

New in version 2.3.

exception ValueError
Raised when a built-in operation or function receives an argument that has the right type but an inappropriate value, and the situation is not described by a more precise exception such as IndexError.

exception VMSError
Only available on VMS. Raised when a VMS-specific error occurs.

exception WindowsError
Raised when a Windows-specific error occurs or when the error number does not correspond to an errno value. The winerror and strerror values are created from the return values of the GetLastError() and FormatMessage() functions from the Windows Platform API. The errno value maps the winerror value to corresponding errno.h values. This is a subclass of OSError.

New in version 2.0.

Changed in version 2.5: Previous versions put the GetLastError() codes into errno.

exception ZeroDivisionError
Raised when the second argument of a division or modulo operation is zero. The associated value is a string indicating the type of the operands and the operation.

The following exceptions are used as warning categories; see the warnings module for more information.

exception Warning
Base class for warning categories.

exception UserWarning
Base class for warnings generated by user code.

exception DeprecationWarning
Base class for warnings about deprecated features.

exception PendingDeprecationWarning
Base class for warnings about features which will be deprecated in the future.

exception SyntaxWarning
Base class for warnings about dubious syntax

exception RuntimeWarning
Base class for warnings about dubious runtime behavior.

exception FutureWarning
Base class for warnings about constructs that will change semantically in the future.

exception ImportWarning
Base class for warnings about probable mistakes in module imports.

New in version 2.5.

exception UnicodeWarning
Base class for warnings related to Unicode.

• Explain effect of flip angle (and why MT effect should be avoided)
• Figure out why saturation pulse power of 0.57 µT works best for tumour and 1 µT works best for normal brain (Source: KD)

I would stick with saturation pulse power between 0.5 uT and around 3 uT
if you are using a standard "CW" type sequence. I found that 0.6 uT was
better for resolving small changes in CEST effect size in tumour, and 1 uT
was better in brain tissue.

• What's the harm in doing pulsed CEST?

To walk through this, the first step is looking at the magnetization of each component of each pool in time. The first, third, and fifth plots are the x,y, and z components for the water pool. They oscillate at a frequency proportional to B1; as the saturation pulse rotates M from the z-direction, the y signal increases and so on.

The oscillating xy-signal looks like this.

Because it is oscillating, it is hard to define it's magnitude. If I just say it is the magnitude at some point in time, I could stop the saturation at a minimum or a maximum or in between, leading to inconsistencies when I apply the imaging sequence. Here is what a z-spectrum looks like with this method:

Not good.

An alternative could be to define the xy-signal as the height of the envelope, which should solve the problem but I'm not sure if that's an accurate representation of reality.

I might be getting ahead of myself, but I've upgraded the simulation to accommodate 6 proton pools in addition to water.

Let me know if there are any specific situations, in terms of concentrations, exchange rates etc, that would be good to run through the multi-pool model.

It works!

Here is my arginine spectrum fitted. I had to give it some reasonable initial parameters, but once I did, I got a pretty good fit (I think).

The trouble is that I had to tell it how many peaks there were, so that's maybe not so useful in vivo.

Next step:

1. construct maps of the arginine data.
2. Shift the spectrum to centre it at 0 (this one is pretty good but others aren't)
3. two rounds of fitting?

a) to get the water peak
b) shift the entire spectrum to centre the water peak at 0
c) subtract the water peak and then try fitting just the species of interest

In Kim's thesis, there is a plot that shows a bi-exponential decay of S/So vs. time. This image is acquired while saturating at a frequency much far off-resonance (20 kHz) - ideally this saturation pulse should have no effect on the signal.

Here's our signal drift:

Here's Kim's signal drift plot from the Bruker scanner in Toronto:

As I've mentioned to a couple of you, I made a function that accepts relaxation times, saturation pulse power, and exchange rates, and tells you how long until the system reaches a steady state saturation. Here are some results using that function.

Some interesting shapes of the curves. If pressed I could come up with physical justification for why these shapes make sense, but if they don't seem right to everybody else, then there must be something wrong in the code.

As a side note, this could be useful for determining exchange rates if the time until steady state saturation is known.

Here's a single voxel and the resulting CEST spectrum:

I see a few peaks: 3ppm, 2.2ppm,
I will be trying to average a few more pixels together to see if I can improve the peaks, and then trying to fit multiple lorentzians to it.

Perhaps @piotrkoz, @DrSAR, or @andrewcyung can comment on whether I should select the whole tumour as an roi and average over it (after shifting the water peak to 0) or get someone to outline a sub-region in the tumour corresponding to a brain region.

I will update this page a I do more work on this... @DrSAR: perhaps you could give some direction on what you're looking for in the talk?

@andrewcyung did a CEST brain scan of a mouse before coming so I thought I'd try and visualize it so we can make sense of it

These data are normalized to the 300 ppm freq. offset so they should be comparable to each other.

Here are andrew's notes on the experiment:

repeat successful mouse brain CEST experiment where I got amide/amine peaks
duplicate successful experiment with turboCEST sequence method
investigate the effect of acceleration factor, with respect to total saturation time over the course of an entire image acquisition for one particular offset

** Data Key **
2,3,4: RARE localizers
5,6: fieldmap and shim voxel
9. quick 3-offset CEST scan to characterize B0 homogeneity
10. reproduce successful experiment with ubcCEST_FLASH (1.0 uT, 72x72, 3.6cm FOV, 490 ms sat pulse, TR=500 ms (i.e. ~10 ms imaging time), -1 to 5 ppm with 0.4 ppm stepsize,
11. try to reproduce same thing but with turboCEST acceleration factor 1 (this time, TR has new meaning; does not include sat pulse.  So TR=10ms)
12. same turboCEST sequence but now with 8 sec of inter-experiment delay (between  frequency offset images)
13. now try #12 but with acceleration factor 2
15. acceleration factor of 1, but saturation pulse length is cut in half.   Therefore, time is same as #13 (accel factor 2) and will be about the same saturation duration from start of scan to centre of k-space
16. same as #15 but sat pulse length is 350 ms (narrow range)
17. same as 15 but reverse encoding so that centre of k-space is acquired last
18. repeat of 12 but with narrower frequency range
19. repeat of 18 but accel = 2
20. repeat of 18 but accel = 4
21. repeat of 18 but accel = 8
22. repeat of 12 but narrower frequency range

For 1 pixel (blue star)

Hey Scott,

Just had a thought...

For our meeting today, it will be worth knowing how long it takes for a single run of the simulation. This will help us determine if there's any need to optimize the code to make it run faster, or parallellize it.

Remember the timeit module I showed you on Python?

Can you run it again and have a breakdown of the potential bottlenecks ready for us?

Firas

Hey stefan,

the "Number of inserted freq offset" feature is inexplicably broken. I cannot enter in any number, whatever number I enter in, it gets converted to 0.

help?

(not urgent)

more for convenience than anything else

also on the wishlist:

-- Sat Freq List --
Start_val = -4.5 ppm
End_val = 4.5 ppm
step_size = 0.2 ppm

--Additional Freqs--
0, 200000, some close to 0

On second thought, this should just automatically create the array, and if you want to insert reference scans in the middle it should be able to do that.

Firas

@andrewcyung has suggested a method to identify peaks by fitting a spline/polynomial to the baseline and then subtracting points from the spline and creating a plot of the difference from the peak and the baseline.

@srveale and I have done that, and here are the results for one scan from last Wednesday ("dCEST_brain.wk2/12")

I'll update this with more details soon, and for the other scans

I need to create a function that fits just the water peak in the CEST spectrum pixel by pixel and shift the spectrum accordingly.

The new spectrum and x-values should be saved as adata.

This is from Piotr's suggestion to average neighbouring voxels to try and increase SNR.

Reproducing the Pulse Power experiment...

Again, I think we're doing well with our choice at 1µT

I've coded the shift_water_peak but it needs some testing in different situations. You can find the code here: https://github.com/firasm/CEST/blob/master/cest.py#L60-L119

@srveale: Since you're fitting data, can you have a look at how well this function works? We need to make sure it's stress tested as we'll be relying on it a lot. Once the function is working properly, we should be using this function instead of the raw data.

Here are some things to watch out for:

• What's the output in extremely noisy situations when there isn't a water peak?
• what happens when the peak is really broad or really narrow?
• I've coded it to only fit peaks that in the range -1 ppm to +1 ppm; do we ever have situations when the peaks are beyond that?
• The normalization condition is a bit awkward, it scales to the highest value in the array because we haven't always acquired a 200 ppm frequency.
• Rather than give a value that means nothing, I'd rather have this function return a numpy.nan or Inf so that we can easily identify it and address it later so if there are any other fringe cases we should consider, let me know

Here are the maps for the scans done on June 24th...

scans_of_interest['EPI_4s_water']=8
scans_of_interest['EPI_4s_1uT']=9
scans_of_interest['EPI_4s_3avg']=10
scans_of_interest['EPI_5s_1uT']=11
scans_of_interest['EPI_0.5s_4avg']=12 #(TR=1s down from 4s)
scans_of_interest['EPI_0.5s_8avg']=13
scans_of_interest['EPI_7s_1uT']=14
scans_of_interest['EPI_10s_1us']=15
scans_of_interest['EPI_5s_1uT_sl2']=16
scans_of_interest['EPI_5s_1uT_sl3']=17

2.5 ppm peak

3.4 ppm peak

In the literature, most people fit CEST spectra to a Super Lorentzian, rather than a lorentzian.

This paper (Xavier Golay is in the list of authors) describes some reasons why the Super Lorentzian is preferred and is a better model for the line shape.

where θ is the dipole orientation angle with respect to the external magnetic field

This is something for us to consider before we get too invested in the Lorentzian line shape when fitting a model to the data.

Sinclair, C. D. J., Samson, R. S., Thomas, D. L., Weiskopf, N., Lutti, A., Thornton, J. S., & Golay, X. (2010). Quantitative magnetization transfer in in vivo healthy human skeletal muscle at 3 T. Magnetic Resonance in Medicine, 64(6), 1739–1748. doi:10.1002/mrm.22562.

Hey Andrew,

Just checking in to solidify plans for tomorrow’s experiment.

Unfortunately I have a lecture from 11-12, so I can be there right after 12 PM - would it be possible you to get started without me so we have the maximum amount of time with you around ?

I can also be there from 10:15 - 10:50 to help set up if you need it. Scott will be there and can help if you if you need it.

Here’s what I was thinking we should do:

• Localize and prescribe a slice (If we could find another place in the body that doesn't move, I would be willing to try that as well)
• CEST with -100ppm, 0 ppm, 100ppm
• Shim
• CEST with -100ppm, 0 ppm, 100ppm
• T1 map
• T2 map
• Coarse CEST spectrum acquisition (also doubles as an SNR check)
• full CEST spectrum

Firas

mouse study name: 'cest_FM.w01'

This is the result from today's brain experiment that Barry and @srveale did for us. @DrSAR, @andrewcyung, and I got some ideas on things to try from some people we've talked to and we were impatient so we thought we'd get started.

What we mean by accelerating of 72x, is that for 72 phase encodes, only one saturation pulse was used. The "delay" in the plot below refers to the waiting time between subsequent saturation pulses (in this case, equivalent to subsequent frequency offsets).

This is still the first mouse, the second mouse is currently being imaged.

For reference, here is the CEST spectrum for the non-accelerated scan (one sat pulse per k-space line) with a 4s pulse length. I've actually shown two reference datasets, one with a low FA (5) and one with a high FA (80) to demonstrate the importance of SNR in the CEST Spectrum. You'll notice that the images look comparable but in reality aren't and the spectra show the noise

The main questions this poses/solves:

1. InterSaturationFreq Delay doesn't seem to matter, so we will likely ditch it.
2. Long Pulse length seems to result in CEST peaks around 3.5 ppm
3. Outstanding question is: how many peaks are there!?
4. It does seem to make a difference whether you do linear offsets or alternating offsets, but the question is which one is the correct one to do? The linear offsets results in 3 (!) discernible peaks around 3.5 ppm whereas the alternating offsets results in a single (maybe 2?) broad peak
5. If you have a long sat pulse (4s) before each k-space line, your CEST signal disappear!? Why !?!?

The plot thickens...

Something is wrong in the code for the new sequence somewhere.

I prescribed 100 offsets in the scan, and then I hit run and started doing stuff.

At some point, I realized that the scanner was no longer scanning, but there was still 48 minutes or so left on the scanner! I called over Stefan and Piotr but they were about to leave. At this point, I shut down ParaVision and restated it. When I came back, my scan was "Not Reconstructed".

Rather than play around, I decided instead to use the old sequence and type in the frequencies manually.

Need to try multiple slices and see what happens

Looks like we were okay using 490 ms for our purse length, but we could go for higher and get a bit better for our peaks!