Can it be acceptable to control the b0/dwi ratio?

Hi all,

First of all, best wishes for the new year.

I understand that dwiintensitynorm scales all b0s and dwi’s of a certain dataset using the same scaling factor. Would it be acceptable to use in some cases a separate factor for b0s and dwi’s? I give an example:

For a certain study I have data from two (similar) protocols. They both acquire 10 x b0, 25 x b700, 40 x b1000 and 75 x b2800. The difference is the way the data is bundled in scan sequences. For the first half of the study each b-shell was acquired in a separate sequence, each time with one b0. A separate sequence with 7 b0 was added. After that, the scanner underwent an update (hate it when that happens during a study). During the second half of the study, the data was acquired in two sequences. The first had 5 b0s, 25 x b700 and 40 x b1000. The other one had 5 b0s and 75 x b2800.

Before I start any analysis I want to make sure the (slight) difference in acquisition has no influence on outcome parameters. Now I noticed that median white matter FA values from the ‘new’ protocol were significantly lower compared to results from the ‘old’ protocol (0.35 compared to 0.40). As the effects I’m looking for are assumed small, such bias is too large, so something needs to be done.
I also noticed that the ratio of b0 intensity over dwi intensity of each b-value differs between protocols. To give an example: for the b2800 data the b0/dwi ratio in the ‘old’ protocol was about 40, for the ‘new’ protocol it was about 55.

Initial attempts to work away the bias without altering the b0/dwi ratio (see also this topic), failed, so I think this ratio is the cause of the bias.
Therefore, in short my approach now is to calibrate b0 intensities (similar to dwiintensitynorm) to a common value, but instead of preserving the b0/dwi ratio, I make sure each dataset obtains the same b-value specific ratio.

I derived the ‘ground truth’ b0 intensity and b0/dwi ratio from the new protocol data alone: I computed median b0 intensity in white matter and the average across subjects. This is the target b0 value. Still from the new protocol data alone I computed the b0/dwi ratio of each b-value. These are the target ratios.
Then, all datasets (new and old protocol data) underwent scaling: each b0 was scaled (individually) towards the target b0 value. The dwi’s were scaled (combined per dataset, not individually) to meet the target b0/dwi ratios.
Using this approach, the median white matter FA values of new and old protocol data did not differ significantly anymore, and so the effect of the protocol was worked away. However, I’ld like to have a second opinion, on whether this approach is acceptable.

Looking forward to any opinions on this.
Cheers,

Thibo

Dear Thibo,

This is not really what I would call a ‘slight’ difference in acquisition.

If you acquire each shell as a separate sequence, each sequence will max out the available gradient strength, resulting in different diffusion encoding delta’s to achieve a specific b-value (and potentially also different echo times). On the other hand, if you acquire all shells in one sequence, the delta’s and echo times will be fixed and the different b-values are obtained by changing only the gradient strength.

The first approach invalidates the assumption in the standard DTI and DKI models… Even though these models are formulated as a function of b, they are actually a function of the gradient strength (and its orientation), with all other things assumed constant. So it is not really surprising that you find significantly different FA values.

Luckily, multi-tissue CSD does not come with this assumption. However, you should still be careful when selecting the appropriate response function for each subject, as there will now be two sequence-specific average responses instead of a single study specific average response. In addition, even with the correct response, it is highly likely that the final fODFs will still be biased as a result of the difference in the sequences.

My suggestion is as follows:

  1. Normalise all the scans to a specific b=0 intensity, as you would normally do
  2. Obtain an average response function for the scans before and one for the scans after the update (sequence-specific response functions).
  3. Perform the deconvolution with the sequence-specific response functions.
  4. Include a covariate in your statistical analysis to distinguish the scans before and after the update.

Step 4 will only work if your pre- and post-update data sets contain a decent mix of both controls and patients (or a mix of time points, if you have a longitudinal study). If this is not the case (e.g. most controls are scanned with one sequence and most patients with the other one), I don’t think there is a good way of extracting credible results from these data.

Cheers,
Ben

Hi Ben,

Thanks for your reply.
In the past we have always acquired shells as separate sequences on our Philips scanner. The way we dealt with the different intensity scaling (diffusion encoding deltas) between the sequences was to set the fixed scale factor type to ‘DOTS’, with a fixed value for all sequences. This allowed us to apply the DTI and DKI model on this data.
We still used this same approach when starting to acquired data in 2 instead of 4 sequences, i.e. the ‘new’ protocol. The intensities in each sequence of a dataset should therefore still share a same scale (same deltas) and allow DTI and DKI (obtaining quantitative results really is our goal).

I think that within-protocol the data is fine. It’s only the between-protocol effect that I’m trying to compensate. Similar to a multi-center study, but then on the same scanner. And therefore I was thinking to control the b0/dwi intensity ratios…
(also, luckily ‘patient’ and ‘controls’ are more or less evenly divided over new and old protocol, so I will for sure include a covariate in the statistics later on)

If the ratios are not constant between acquisitions, I would suspect that there are other changes to the sequences than just gradient strength (did you compare all the parameters?) or that the DICOMs somehow did not convert correctly (have you compared the output of mrconvert to that of other tools?).

I agree with @bjeurissen, there’s a potentially quite fundamental problem with the data here. You stated earlier that:

Presumably this is on a Philips system? These scanners have a habit of ‘optimising’ everything behind your back - which makes them simple to operate most of the time, but can also screw things up quite badly. In this case, I’d be very surprised if your echo times matched between b-values, and as @bjeurissen already mentioned, your diffusion gradient timings (δ & Δ) will also most likely be completely different. I’d also expect your repetition times to differ. There’s no amount of fixing the scale factor with DOTS or anything else that is going to get around the fundamental changes in contrast that these variations in δ, Δ, TE & TR will introduce between shells…

Personally, I don’t think there’s any way you could defensibly analyse these data together - certainly not with a model that doesn’t also account for this (and pretty much nothing does, as far as I know). If your different shells inherently have different δ, Δ, TE & TR, their relative signal scaling will be off in ways that are going to be very difficult to correct for. This is particularly the case anywhere you have contributions to the signal from a mixture of components with different T2 values, since their relative contributions to the signal will change with TE (and a similar argument holds for T1 with varying TR). For example, the signal fraction that gets assigned to CSF will be different with a different echo time, simply because the tissue signal has a shorter T2 than the CSF - nothing to do with diffusion as such.

As @bjeurissen also mentioned, theoretically you could nonetheless deal with this to some extent with multi-tissue CSD, since it doesn’t really care about this - all it needs to know is what is the expected signal for each component within each shell. This is where the ability to ‘train’ the algorithm (i.e. get the response functions) specifically for your data really helps. But that would only work if your protocol had remained constant - it looks like your parameters will be quite different again after the upgrade: at the very least, I’d expect the first half of your sequence to now share the same δ, Δ, TE & TR, when these shells wouldn’t have before the upgrade. Assuming this is all correct, it’s not at all surprising you end up with different results.

While you could try to correct for this by scaling relative to the b=0 images, I have a strong suspicion even that is problematic. If your b=0 volumes were all acquired with different TE & TR, the intensity of the tissue will vary relative to the CSF, purely because CSF has very different T1 & T2 values from parenchyma. So even that simple operation becomes ill-defined. At the very least, you need to restrict the voxels being used to measure the appropriate scale factor to pure white matter (no CSF, GM, background, or anything else). Currently, your b=0 / DWI ratios seem far too large, I’d expect values less than 10 in white matter - so I suspect you’re including the CSF in that measurement, and that will definitely skew things differently due to TE & TR differences…

So I’m not sure what to tell you here, other than I don’t envy your position… There is no way to ‘fix’ this that I can think of. Anything you do will require a lot of effort, and in the end you’re unlikely to get the data published - I know what my verdict would be if I got asked to review something like this… I guess the first thing would be to look at what the TE & TR values were in these scans - if they match across shells, and pre & post upgrade, then you might be in with a chance - but frankly, the evidence so far suggests otherwise. Sorry…

Hi Donald,

Thanks for your elaborate reply.
I hope I do not need to recall all previous publications we had with the ‘old’ protocol :slight_smile:

I just checked the dicom files and TE and TR of all sequences of new and old protocol are the same (TE=90ms, TR=7800ms).
Also, I mentioned the wrong b0/dwi ratios in my previous post. Those were results based on voxels in CSF (one of my earlier attempts was based on CSF masks). The ones based on WM voxels are between 1 and 5, but still somewhat higher in the new protocol compared to the old protocol.

Can I retrieve the gradient diffusion timings from the dicom file? Could not find them immediately…

In any case, thanks @bjeurissen for the suggestion of using protocol-specific response function for analysis of AFD. For DTI (I know, old model, but still asked for) I could probably just use the highest b-shell. Although not ideal for DTI, it’s the only one that was consistently acquired as a separate sequence.

Well that’s great news! There should be no problem analysing your data in that case. Sounds like you can safely ignore my previous doomsday speculations… Not sure why you’d be having so much trouble re-scaling all the data in that case though?

I don’t think you can get them - there’s certainly no DICOM standard tags for this. But then Philips do use a lot of custom tags to store all kinds of other information, so it might be in there somewhere… Realistically though, I don’t anticipate the gradient timings would be any different if you used the same TE, at least for the same software version - but I can’t say for sure…

OK, that makes more sense. Is there any obvious difference in SNR? That would certainly influence that ratio for the high b-values.

That makes less sense - unless there is an obvious difference in SNR. The best thing to do would be to take the ratio of the mean b=0 to the mean DW images over each shell - that tends to produce much more stable and homogeneous values over the white matter. You can use a combination of dwiextract and mrmath for this, e.g.:

dwiextract dwi.mif -bzero - | mrmath - mean -axis 3 mean_b0.mif
dwiextract dwi.mif -shell 2800 - | mrmath - mean -axis 3 mean_b2800.mif
mrcalc mean_b0.mif mean_b2800.mif -div - | mrstats - -mask white_matter_mask.mif -output mean

I’d still make a concerted effort at figuring out these differences in scaling. If there is an issue, you’d want to know now rather than later… Given what you’ve shown, I can’t think of any reasons why these intensities would differ particularly - but you’d need to have rescaled the data properly first, not sure what procedure you’ve been using so far. But given that you have b=0 images within each scan, whatever procedure you use must scale each acquisition with a single global scale factor (i.e. you shouldn’t scale the b=0 volumes differently from the DW volumes if they were acquired in the same scan), and should provide you with data where the resulting b=0 intensities come out constant across acquisitions, particularly within the parenchyma - any other outcome would signal some kind of failure somewhere…