Troubles in the Transverse Plane

[srveale]

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.

[DrSAR]

What is t=0 in these plots? Broadly speaking the signal in the xy-plane is an oscillation, modulated by a decay times the amplitude (maybe an offset). So the 1st point (t=0) should give you the required amplitude.

[DrSAR]

problem: t=0 is the start of saturation - I think I now understand better. disregard previous comment

[andrewcyung]

what is the reference frame of your simulation rotating at, relative to the laboratory frame: Is it at zero offset or at the frequency of the second pool? Have you tried simulating right at zero offset, with no exchange to the small pool? If things don't make sense, I try to simulate the simplest case possible to see if that makes sense.

[DrSAR]

I think the thing is that the oscillation do make sense. There is a low-power saturation pulse applied off-resonance for water (somewhere in the spectrum) and as a result you see the water magnetization oscillating around (rotating around an effective B1,eff). The problem is that for the long saturation pulse the water whizzes around many times and hence the value of any of the components at the end of that pulse becomes very sensitive to the power, off-resonance and duration of that pulse. In reality we don't see that.

My suspicion is that this sensitivity is 'damped' since a spread of off-resonance conditions is encountered smearing out the oscillation resulting in an average value. I think Scott is trying to see whether averaging over an ensemble of magnetizations with slightly varied frequency offsets (as could be expected within a pixel) gives a reasonable average value. [This is how simulation times blow up quickly ...]

[srveale]

Just an update into how things are going. For each saturation offset, delta, I told the simulation to run for a range of delta ± 20 rad/s, and used the average of the signals as the data point. I've tried using a variety of ranges and the graph shows similar features. I've also tried using a range over omega1 (the resonant frequency of water nuclei under the influence of B1), and the graph is still the same.

I've also tried using different flip angles with results that made no sense.

In the next figure, I simply told the simulation that the M components after the saturation pulse should be the average of their values over the last cycle of rotation (see the first figure in this issue to see what I mean by rotation). There are still troubles but it looks much better.

I will have to determine: -why the shoulder is so sharp -why there is a mirror peak even though the frequencies are linear (there might be ghosts in the code) -what is happening around delta = 0

I think these are manageable, and I think progress has been made on this issue.

[firasm]

For the last simulation, can you repeat with a shorter T1? Try 1/3 as high

On Sat, May 23, 2015 at 8:50 PM, srveale notifications@github.com wrote:

Just an update into how things are going. For each saturation offset, delta, I told the simulation to run for a range of delta ± 20 rad/s, and used the average of the signals as the data point. I've tried using a variety of ranges and the graph shows similar features. I've also tried using a range over omega1 (the resonant frequency of water nuclei under the influence of B1), and the graph is still the same. I've also tried using different flip angles with results that made no sense. In the next figure, I simply told the simulation that the M components after the saturation pulse should be the average of their values over the last cycle of rotation (see the first figure in this issue to see what I mean by rotation). There are still troubles but it looks much better. I will have to determine: -why the shoulder is so sharp -why there is a mirror peak even though the frequencies are linear (there might be ghosts in the code) -what is happening around delta = 0

I think these are manageable, and I think progress has been made on this issue.

Reply to this email directly or view it on GitHub: https://github.com/firasm/CEST/issues/37#issuecomment-104975893

[srveale]

Divided T1 by 3 and ran it again:

[firasm]

Brilliant !

That's as @DrSAR predicted. What's the T1 value after dividing by 3?

[srveale]

Water T1: 666 ms Species T1: 333 ms I should also say that this was a 90 degree flip angle

[DrSAR]

This occurred to me a little late: When you average over a period (integrating the vector of magnetization over time) you are doing the same as averaging over a range of frequencies for a fixed time (.e.g. the end of the sat pulse). That frequency range, naturally, is the inverse of the period from above. The advantage is that for a period-average you don't need to simulate a range of frequencies.

[DrSAR]

These splittings worry me. There are two peaks where there should be one and there is one peak where there should be none. Are there three new little peaks? Can you change the chemical shift of the species, i.e. move it closer or farther from water and leave the rest the same?

[srveale]

I've been trying to unravel the mystery of these ghost peaks. Failing that, let me see if I can define the problem properly. There's a lot going on so this might be better in person. Narrowing it down:

-The peaks remain unaffected if exchange is turned off -They also remain unaffected if the resonance frequency of the species is varied -They disappear with increasing inter-frequency delay

This is the z-spectrum around the ghost peak. -3.35 ppm corresponds to ~-6300rad/s offset frequency. Remember that exchange is zero for this:

Normally I would take the average of a full cycle of oscillation, by finding the previous two peaks. The plateau at ~6300 is where that averaging function failed, and I just plotted some placeholder values. Why? There weren't two peaks to find.

This is the time evolution of magnetization of each component, plotted for a few deltas about where the ghost peak is located:

Look closely at Mya. For an interval of deltas, the frequency of oscillations is too slow to go through two full cycles, so using the average over the last two peaks fails.

It can't be that the frequency of these oscillations is sqrt(delta^2 + omega1^2). Maybe it should be, but that's not what the sim is giving. The oscillation frequency slows down around the -6300 region. Here's what the oscillation periods look like as a function of delta:

The gap is where there weren't enough peaks to get a period measurement.

So something is special about delta = ±6250 and I don't know what it is. My gut tells me that the small scale oscillations (that don't go away unless I average over the previous cycle and a range of deltas for each point) are related to the ghost peaks. I think they must be an artefact of these weirdly behaving oscillations.

[firasm]

Whoa - weird. The fact that this isn't exchange related is relieving and also strange.

You're right this may be better in person. @DrSAR are you available this afternoon?

I'm just gonna drop Hira off at the airport and make my way in. Should be there by 1230 today. If that complicates things, feel free to meet without me.

[srveale]

@DrSAR and I got to the bottom of this yesterday. The integrator we were using had sampling frequency issues. Building and using the simplest integrator possible enabled the production of this graph:

This is what we are looking for. Unfortunately, this graph took an hour to make. I will start producing graphs to show the dependencies on inter-frequency delay and B1, but I think this issue has been resolved.

[firasm]

So we can't use any other integrator period? I'm assuming it took an hour because of your integrator - is that correct? Ultimately I'm wondering if every full simulation experiment take an hour to run?

Also can you comment on How much of our old simulation work will we need to revisit?

[firasm]

P.S. Nice work solving the issue!

[srveale]

Yes, it took a long time because of the rudimentary integrator. I'm going to upgrade it to an RK4 method and that should bring the time down, but I'm now wary of using pre-built integrators with adaptive step size because I think that is the source of this problem.

[DrSAR]

You will have to chose a step-size for RK4 (or any other integrator) and somehow avoid the trap that the stock algorithm fell into. Good work and I hope this can be sped up in the fullness of time. S

srveale mailto:notifications@github.com 2015 May 30, at 13:01

Yes, it took a long time because of the rudimentary integrator. I'm going to upgrade it to an RK4 method and that should bring the time down, but I'm now wary of using pre-built integrators with adaptive step size because I think that is the source of this problem.

— Reply to this email directly or view it on GitHub https://github.com/firasm/CEST/issues/37#issuecomment-107079867.

srveale mailto:notifications@github.com 2015 May 21, at 9:37

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.

each component vs time https://cloud.githubusercontent.com/assets/7737922/7753490/9ffd19dc-ff9a-11e4-8b5f-da95837686a7.png

The oscillating xy-signal looks like this.

xy signal vs time https://cloud.githubusercontent.com/assets/7737922/7753675/d5557308-ff9b-11e4-8fdc-711423a4774c.png

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:

xy z-spectrum https://cloud.githubusercontent.com/assets/7737922/7753743/5c831ff6-ff9c-11e4-9283-a230a2349cdc.png

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.

— Reply to this email directly or view it on GitHub https://github.com/firasm/CEST/issues/37.