PD calculation without OLS broken if nr_echoes_forA does not equal nr_echoes4avg

[lukeje]

Background

In versions prior to v0.2.5, the PD (A) map was computed using R1 (rather, T1) and the T1w data (i.e. the PDw data was only used in the calculation of the R1 map). The map calculation was implemented in such a way as to allow a different number of echoes to be used for the A calculation using the hmri_def.PDproc.nr_echoes_forA variable if extrapolation to TE=0 was turned off. As only the T1w data and not the PDw data was needed, a special T1w_forA volume was created which was averaged over nr_echoes_forA echoes and input into the A calculation.

Problem

From v0.2.5 onwards, the PD calculation uses both T1w and PDw data. Therefore in order to correctly account for R2* effects, the PDw data would also need to use only nr_echoes_forA echoes if extrapolation to TE=0 was turned off.

Possible fixes

• Remove nr_echoes_forA entirely and just use the averaged PDw, T1w and MTw used for registration (which averages the same number of echoes for all contrasts)
• Introduce PDw_forA (and MTw_forA) which are also averaged over nr_echoes_forA and used for A calculation
• Keep things as they are, but apply R2* correction on averaged PDw, MTw, and T1w_forA data separately rather than correcting A, using the echoes actually used for each average
• Go back to computing A just using R1 map and T1w data like in v0.2.5

[lukeje]

@mfcallaghan @ChristophePhillips @antoinelutti Do you have a strong opinion which solution would be best?

[ChristophePhillips]

@lukeje Honestly, no idea...
My answer consists in two questions:

• why was the PD calculation changed, i.e. what is the gain in using both T1w and PDw?
• if that was useful (which I assume it is...), then are the implications of using your proposed solutions 1, 2 or 3?

[lukeje]

@ChristophePhillips There were a couple of overlapping reasons for the code to change:

1. to match what is in the hMRI toolbox paper
2. to increase modularisation and make unit testing easier (the R1 and PD calculations were moved to external files)
3. to implement a version of the equations allowing for larger flip angles

In principle 2. and 3. could also have been done within the old approach (rewriting the PD calculation function as necessary), but 1. is something that we can't easily get around.

Pros and cons of the different solutions to follow in the next comment...

[lukeje]

As promised, here's a list of what I see as the pros and cons of the different options.

I'm ignoring here that only option C (or option B modified in line with C) would allow for different echo times between the contrasts (rather than just different numbers of echoes acquired with the same echo spacing), as I think that is a separate issue.

A: Remove nr_echoes_forA entirely and just use the averaged PDw, T1w and MTw used for registration (which averages the same number of echoes for all contrasts)

Pros:

• simple to implement
• simplifies options for user

Cons:

• user may want to boost SNR for registration (use max echoes) but minimise R2* bias in PD map (use fewer echoes)
• user may want to use more echoes to compute PD than are found in the MTw acquisition

B: Introduce PDw_forA (and MTw_forA) which are also averaged over nr_echoes_forA and used for A calculation

Pros:

• keeps all contrasts symmetric

Cons:

• extra complexity in the code
• number of echoes to use restricted by minimum number between T1w, PDw and MTw

C: Keep things as they are, but apply R2* correction on averaged PDw, MTw, and T1w_forA data separately rather than correcting A, using the echoes actually used for each average

Pros:

• relatively simple to rearrange the code
• different numbers of echoes for each contrast accounted for

Cons:

• would be going beyond what is in the papers
• would add additional complexity
• would contaminate R1 and MTsat with any R2* mapping errors (rather than just A)

D: Go back to computing A just using R1 map and T1w data like in v0.2.5

Pros:

• only need to change PD mapping function to optionally accept R1 map instead of PDw data
• T1w_forA stays as it was

Cons:

• PD mapping would not represent what was written in the hMRI toolbox paper
• could be unknown numerical issues

E: Reduce echo selection to a single parameter nr_echoes4avg specified by the user

(Suggested by @antoinelutti)
Pros:

• simplifies code
• simple to implement

Cons (similar to D):

• user may want to boost SNR for registration (use max echoes) but minimise R2* bias in PD map (use fewer echoes)
• user may want to use more echoes to compute PD than are found in the MTw acquisition
• need to rearrange the code to use this user-defined parameter

[antoinelutti]

Hi @lukeje , @ChristophePhillips and @mfcallaghan
I wasn't aware that the PD estimation now involved both the PDw and T1w data.
The original motivation for the specific number of echo images in the PD calculation is that the R2* bias due to the averaging is expected to be stronger for this map. In the case of e.g. R1 estimation, the R2* bias is approximately consistent between the PDw and T1w contrasts, which reduces the bias on the R1 estimates. That's not the case for the PD estimates - or at least wasn't when only the PDw data was used.
Despite this, my personal preference would go to keeping only one variable to determine the number of echoes used for map estimation when interpolation to TE=0 is turned off (nr_echoes4avg) to facilitate the management of the toolbox, code readability and interpretability of results.

[Barisevrenugur]

unfortunately, I think I (still) did not understand this one. But I came across the following lines in create_MTprot which seems related:

hMRI-toolbox/hmri_create_MTProt.m

Line 1399 in 894a38b

% consistency check for number of echoes averaged for A calculation:

and

hMRI-toolbox/hmri_create_MTProt.m

Line 1220 in 894a38b

mpm_params.nr_echoes4avg = min(length(find(mpm_params.input(1).TE<maxTEval4avg))+1,ncommonTEvals);

now, in the code following L1399 nr_echoes_forA is set to the (smaller number) of T1w echoes if there is such a problem.

and in the code following L1220 it simply

hMRI-toolbox/hmri_create_MTProt.m

Line 1218 in 894a38b

% find maximum number of echoes that are common to all available contrasts

why do not the code following the above hinted lines fix the problem of non-equal number of echoes for averaging? but they seem to? with respect to processing, they happen pretty early in the code, namely already in get_mpm_params, which reads off the params at the beginning of mt_prot

[lukeje]

why do not the code following the above hinted lines fix the problem of non-equal number of echoes for averaging?

Because we are talking about nr_echoes_forA, which is not necessarily equal to nr_echoes4avg. My suggested solution A would behave as you suggested, but the current code does not.

[lukeje]

@antoinelutti

Despite this, my personal preference would go to keeping only one variable to determine the number of echoes used for map estimation when interpolation to TE=0 is turned off (nr_echoes4avg) to facilitate the management of the toolbox, code readability and interpretability of results.

I guess this is an alternative "E" closely related to "A": rather than defining nr_echoes4avg automatically, let the user set it, and remove nr_echoes_forA. This would also work...