Implementing an iterative delay correction strategy to improve CVR estimation in clinical populations

[beccaclements99]

This code implements the iterative delay correction approach described in this preprint. This approach is useful for clinical cohorts when the appropriate maximum shift is not known. We found that using a maximum shift that is too small can result in many delay estimates at the boundary (meaning that they are probably not optimized), but a conservatively large maximum shift can result in spurious negative CVR values. To address this, the delay search range is selectively expanded only for voxels with estimated delays at a boundary (i.e., within 1 lag step of the minimum or maximum shift) until the estimated delay is no longer at a boundary or a predefined value is reached.

Proposed Changes

• Two new parameters were introduced: a starting maximum lag and an increment for increasing the maximum shift
• If these parameters are defined, and a voxel's lag value is within 1 lag step of the minimum or starting maximum lag, the maximum lag for that voxel will be iteratively increased by the amount specified with --lag-increment until the optimal lag is not at within 1 lag step of the boundary or --lag-max (the final lag maximum) is reached, whichever comes first.
• If a starting maximum lag is not defined, the iterative approach is not applied and the standard (fixed maximum) approach is used
• Note that only the maximum lag is expanded, not the minimum lag. If the bulk shift is estimated from normal-appearing tissue or, in cases of pathology such as Moyamoya disease or carotid stenosis, from unaffected vascular territories, hemodynamic responses more than 9 s earlier than the bulk shift are not expected to be physiologically plausible (and thus an iterative expansion of the minimum lag would not be necessary). Rare exceptions like arteriovenous shunting would require additional modifications to the code.
• An error is returned if:
  - A starting_lag_max is provided and no lag_increment is specified
  - A starting_lag_max is provided and it is greater than or equal to lag_max
  - A lag map is provided AND starting_lag_max is specified. (The lag map option does not work with the iterative lag search)

Change Type

• bugfix (+0.0.1)
• minor (+0.1.0)
• major (+1.0.0)
• refactoring (no version update)
• test (no version update)
• infrastructure (no version update)
• documentation (no version update)
• other

Checklist before review

• I added everything I wanted to add to this PR.
• [Code or tests only] I wrote/updated the necessary docstrings.
• [Code or tests only] I ran and passed tests locally.
• [Documentation only] I built the docs locally.
• My contribution is harmonious with the rest of the code: I'm not introducing repetitions.
• My code respects the adopted style, especially linting conventions.
• The title of this PR is explanatory on its own, enough to be understood as part of a changelog.
• I added or indicated the right labels.

• I added information regarding the timeline of completion for this PR.

• Please, comment on my PR while it's a draft and give me feedback on the development!

[smoia]

@beccaclements99 I saw that pre-commit didn't run properly - I would suggest you install it locally for it to run when you commit (see here).

[smoia]

I also have a couple of initial questions for clarification:

1. Could you explain the rationale behind the fact that the maximum lag expands when the minimum lag is close to boundaries?
   2. Is it necessary to have an increment for the maximum lag that is different from the delta lag already specified by the user? Scratch that, I understood why. Maybe changing the name of the flag and the variable would be better though.

[smoia]

Suggested change

The empty line will be flagged by Ruff

[smoia]

Suggested change

	lag = (lag_idx * step) / freq + (mask * lag_min)
	else:
	lag_idx = np.argmax(r_square_all, axis=-1)
	lag = (lag_idx * step) / freq + (mask * lag_min)
	else:
	lag_idx = np.argmax(r_square_all, axis=-1)
	lag = (lag_idx * step) / freq + (mask * lag_min)

Line 795 and 798 are the same, so we can just bring them out of the if statement ;)

[smoia]

Suggested change

	return final_lag_idx
	return final_lag_idx

Ruff will flag not having 2 empty lines

[smoia]

It doesn't look like you are using freq in the function - can you have another look at it please?

[smoia]

@beccaclements99 besides the couple of feedback points I already left, there are a couple of things that I'd like to discuss regarding the implementation - bearing in mind that we might have to modularise the workflow a bit more than previously planned.

1. I find a little confusing that you are specifying the possible extra lag analysis as you would the "normal" one, and you are adding a first soft stop to emulate what is the normal behaviour of phys2cvr. I understand why you are doing it at the moment (it's a neat trick), but from a user perspective (and a maintenance perspective), since this approach is an addition to the standard workflow, it would make more sense to treat it as such extra. That means, keep lag max as classically implemented (also because in case the lag min is not specified, its negative would still be the minimum lag - I have a suspicion that you relied on the asymmetric lag the whole time, right?), and have the extra parameter allowing more exploration if and when required, as a true expansion. Granted, that would mean having to run both the regressor creation and the lag regression a second time, but it would only be for the new exploratory space alone.

2. Going down this path of "extra exploration" would allow you to only work within a very confined set of voxels, those that hit the lag boundary, thus making this second execution far faster and computationally lighter than keep solving OLS for the whole volume. It would also make the use of the recursive incremental exploration (if that is the idea behind lag_increment) far easier and better integrated in the existing codebase, without having to change masks the whole time.

3. It would, however, mean to break the current workflow into different blocks (it is not a bad idea in general, this workflow is massive after all), one of which deals with the recursive lag estimation. I'm happy to help you do this though, and i think it would become very beneficial on the long run.

4. I would change the name of the parameter used to increase the exploration space (lag_increment) to make it very distinct from lag_step, but I am also wondering if we really need it, or if we couldn't use trial-length (or one of its fractions), an already existing parameter that is used to indicate the length of one trial repetition, to express the maximum expansion allowed/used.

Let me know what do you think and how I can be of help!