Implement general kinetic model for multi-PLD data

aslprep/data/aslprep_bids_config.json

@@ -37,9 +37,9 @@
     "sub-{subject}[/ses-{session}]/[{datatype<func|beh>|func}/]sub-{subject}[_ses-{session}]_task-{task}[_acq-{acquisition}][_rec-{reconstruction}][_run-{run}][_echo-{echo}][_recording-{recording}]_{suffix<physio|stim>}{extension<.tsv.gz|.json>|.tsv.gz}",
     "sub-{subject}[/ses-{session}]/{datatype<func|perf>|func}/sub-{subject}[_ses-{session}][_task-{task}][_acq-{acquisition}][_ce-{ceagent}][_dir-{direction}][_rec-{reconstruction}][_run-{run}][_echo-{echo}][_space-{space}][_atlas-{atlas}][_cohort-{cohort}][_desc-{desc}]_{suffix<regressors|timeseries>|timeseries}{extension<.json|.tsv|.csv|>|.tsv}",
     "sub-{subject}[/ses-{session}]/{datatype<perf>|perf}/sub-{subject}[_ses-{session}][_task-{task}][_acq-{acquisition}][_ce-{ceagent}][_dir-{direction}][_rec-{reconstruction}][_run-{run}][_echo-{echo}][_desc-{desc}]_{suffix<confounds>}{extension<.tsv|.json>|.tsv}",
-    "sub-{subject}[/ses-{session}]/{datatype<perf>|perf}/sub-{subject}[_ses-{session}][_task-{task}][_acq-{acquisition}][_ce-{ceagent}][_dir-{direction}][_rec-{reconstruction}][_run-{run}][_echo-{echo}][_space-{space}][_res-{res}][_atlas-{atlas}][_cohort-{cohort}][_desc-{desc}]_{suffix<asl|aslref|att|cbf|coverage|mask>}{extension<.nii|.nii.gz|.json|.csv|.tsv>|.nii.gz}",
-    "sub-{subject}[/ses-{session}]/{datatype<perf>|perf}/sub-{subject}[_ses-{session}][_task-{task}][_acq-{acquisition}][_ce-{ceagent}][_dir-{direction}][_rec-{reconstruction}][_run-{run}][_echo-{echo}][_space-{space}][_res-{res}][_den-{den}][_cohort-{cohort}][_desc-{desc}]_{suffix<asl|att|cbf>}{extension<.dtseries.nii|.dscalar.nii|.json|.dscalar.json|.dtseries.json|.func.gii|.func.json>|.dtseries.nii}",
-    "sub-{subject}[/ses-{session}]/{datatype<figures>|figures}/sub-{subject}[_ses-{session}][_task-{task}][_acq-{acquisition}][_ce-{ceagent}][_dir-{direction}][_rec-{reconstruction}][_run-{run}][_echo-{echo}][_space-{space}][_res-{res}][_den-{den}][_cohort-{cohort}][_desc-{desc}]_{suffix<asl|aslref|att|cbf|basil|score|scrub|pvc|mask>}{extension<.svg|.png|.html>|.svg}",
+    "sub-{subject}[/ses-{session}]/{datatype<perf>|perf}/sub-{subject}[_ses-{session}][_task-{task}][_acq-{acquisition}][_ce-{ceagent}][_dir-{direction}][_rec-{reconstruction}][_run-{run}][_echo-{echo}][_space-{space}][_res-{res}][_atlas-{atlas}][_cohort-{cohort}][_desc-{desc}]_{suffix<asl|aslref|abat|abv|att|cbf|coverage|mask>}{extension<.nii|.nii.gz|.json|.csv|.tsv>|.nii.gz}",
+    "sub-{subject}[/ses-{session}]/{datatype<perf>|perf}/sub-{subject}[_ses-{session}][_task-{task}][_acq-{acquisition}][_ce-{ceagent}][_dir-{direction}][_rec-{reconstruction}][_run-{run}][_echo-{echo}][_space-{space}][_res-{res}][_den-{den}][_cohort-{cohort}][_desc-{desc}]_{suffix<asl|abat|abv|att|cbf>}{extension<.dtseries.nii|.dscalar.nii|.json|.dscalar.json|.dtseries.json|.func.gii|.func.json>|.dtseries.nii}",
+    "sub-{subject}[/ses-{session}]/{datatype<figures>|figures}/sub-{subject}[_ses-{session}][_task-{task}][_acq-{acquisition}][_ce-{ceagent}][_dir-{direction}][_rec-{reconstruction}][_run-{run}][_echo-{echo}][_space-{space}][_res-{res}][_den-{den}][_cohort-{cohort}][_desc-{desc}]_{suffix<asl|aslref|abat|abv|att|cbf|basil|score|scrub|pvc|mask>}{extension<.svg|.png|.html>|.svg}",
     "sub-{subject}/{datatype<figures>}/sub-{subject}[_ses-{session}][_acq-{acquisition}][_ce-{ceagent}][_rec-{reconstruction}][_run-{run}][_space-{space}][_atlas-{atlas}][_cohort-{cohort}][_desc-{desc}]_{suffix<T1w|T2w|T1rho|T1map|T2map|T2star|FLAIR|FLASH|PDmap|PD|PDT2|inplaneT[12]|angio|dseg|mask|dwi|epiref|fieldmap>}{extension<.html|.svg|.png>}",
     "sub-{subject}/{datatype<figures>}/sub-{subject}[_ses-{session}][_acq-{acquisition}][_ce-{ceagent}][_rec-{reconstruction}][_run-{run}][_space-{space}][_atlas-{atlas}][_cohort-{cohort}][_desc-{desc}]_{suffix<dseg|mask|dwi|epiref|fieldmap>}{extension<.html|.svg|.png>}",
     "sub-{subject}[/ses-{session}]/{datatype<anat>|anat}/sub-{subject}[_ses-{session}][_acq-{acquisition}][_ce-{ceagent}][_rec-{reconstruction}][_run-{run}]_from-{from}_to-{to}_mode-{mode<image|points>|image}_{suffix<xfm>|xfm}{extension<.txt|.h5|.json>}",

aslprep/data/boilerplate.bib

@@ -852,3 +852,13 @@ @article{topup
     doi = {10.1016/S1053-8119(03)00336-7},
     url = {https://www.sciencedirect.com/science/article/pii/S1053811903003367}
 }
+
+@article{woods2023recommendations,
+  title={Recommendations for quantitative cerebral perfusion MRI using multi-timepoint arterial spin labeling: Acquisition, quantification, and clinical applications},
+  author={Woods, Joseph G and Achten, Eric and Asllani, Iris and Bolar, Divya S and Dai, Weiying and Detre, John A and Fan, Audrey P and Fern{\'a}ndez-Seara, Maria and Golay, Xavier and G{\"u}nther, Matthias and Guo, Jia and Hernandez-Garcia, Luis and Ho, Mai-Lan and Juttukonda, Meher R. and Lu, Hanzhang and MacIntosh, Bradley J. and Madhuranthakam, Ananth J. and Mutsaerts, Henk J.M.M. and Okell, Thomas W. and Parkes, Laura M. and Pinter, Nandor and Pinto, Joana and Qin, Qin and Smits, Marion and Suzuki, Yuriko and Thomas, David L. and van Osch, Matthias J.P. and Wang, Danny J.J. and Warnert, Esther A.H. and Zaharchuk, Greg and Zelaya, Fernando and Zhao, Moss and Chappell, Michael A.},
+  year={2023},
+  journal={OSF Preprints},
+  publisher={OSF Preprints},
+  doi={10.31219/osf.io/4tskr},
+  url={https://doi.org/10.31219/osf.io/4tskr},
+}

aslprep/interfaces/cbf.py

@@ -24,12 +24,7 @@
     estimate_labeling_efficiency,
     pcasl_or_pasl,
 )
-from aslprep.utils.cbf import (
-    _getcbfscore,
-    _scrubcbf,
-    estimate_cbf_pcasl_multipld,
-    estimate_t1,
-)
+from aslprep.utils.cbf import _getcbfscore, _scrubcbf, estimate_t1, fit_deltam_multipld
 
 
 class _RefineMaskInputSpec(BaseInterfaceInputSpec):
@@ -342,6 +337,16 @@ class _ComputeCBFOutputSpec(TraitedSpec):
         None,
         desc="Arterial transit time map, in seconds. Only generated for multi-delay data.",
     )
+    abat = traits.Either(
+        File(exists=True),
+        None,
+        desc="Arterial bolus arrival time map, in seconds. Only generated for multi-delay data.",
+    )
+    abv = traits.Either(
+        File(exists=True),
+        None,
+        desc="Arterial blood volume map. Only generated for multi-delay data.",
+    )
     plds = traits.Either(
         File(exists=True),
         None,
@@ -503,31 +508,47 @@ def _run_interface(self, runtime):
             # Broadcast PLDs to voxels by PLDs
             plds = np.dot(plds[:, None], np.ones((1, deltam_arr.shape[0]))).T
 
-        if is_casl:
-            tau = np.array(metadata["LabelingDuration"])
-
         if is_multi_pld:
+            ti1 = None
+            tau = None
             if is_casl:
-                att, mean_cbf = estimate_cbf_pcasl_multipld(
-                    deltam_arr,
-                    scaled_m0data,
-                    plds,
-                    tau,
-                    labeleff,
-                    t1blood=t1blood,
-                    t1tissue=t1tissue,
-                    unit_conversion=UNIT_CONV,
-                    partition_coefficient=PARTITION_COEF,
-                )
+                tau = metadata["LabelingDuration"]
 
             else:
-                # Dai's approach can't be used on PASL data, so we'll need another method.
-                raise ValueError(
-                    "Multi-delay data are not supported for PASL sequences at the moment."
-                )
+                if metadata["BolusCutOffTechnique"] == "QUIPSSII":
+                    # PASL + QUIPSSII
+                    # Only one BolusCutOffDelayTime allowed.
+                    assert isinstance(metadata["BolusCutOffDelayTime"], Number)
+                    ti1 = metadata["BolusCutOffDelayTime"]
+
+                elif metadata["BolusCutOffTechnique"] == "Q2TIPS":
+                    # PASL + Q2TIPS
+                    # Q2TIPS should have two BolusCutOffDelayTimes.
+                    assert len(metadata["BolusCutOffDelayTime"]) == 2
+                    ti1 = metadata["BolusCutOffDelayTime"][0]
 
-            mean_cbf_img = masker.inverse_transform(mean_cbf)
+                else:
+                    raise ValueError(
+                        f"Unsupported BolusCutOffTechnique ({metadata['BolusCutOffTechnique']}) "
+                        "for multi-PLD data."
+                    )
+
+            cbf, att, abat, abv = fit_deltam_multipld(
+                deltam_arr,
+                scaled_m0data,
+                plds,
+                labeleff,
+                t1blood=t1blood,
+                partition_coefficient=PARTITION_COEF,
+                is_casl=is_casl,
+                tau=tau,  # defined for (P)CASL
+                ti1=ti1,  # defined for PASL
+            )
+
+            mean_cbf_img = masker.inverse_transform(cbf)
             att_img = masker.inverse_transform(att)
+            abat_img = masker.inverse_transform(abat)
+            abv_img = masker.inverse_transform(abv)
 
             # Multi-delay data won't produce a CBF time series
             self._results["cbf_ts"] = None
@@ -537,9 +558,22 @@ def _run_interface(self, runtime):
                 newpath=runtime.cwd,
             )
             att_img.to_filename(self._results["att"])
+            self._results["abat"] = fname_presuffix(
+                self.inputs.deltam,
+                suffix="_abat",
+                newpath=runtime.cwd,
+            )
+            abat_img.to_filename(self._results["abat"])
+            self._results["abv"] = fname_presuffix(
+                self.inputs.deltam,
+                suffix="_abv",
+                newpath=runtime.cwd,
+            )
+            abv_img.to_filename(self._results["abv"])
 
         else:  # Single-delay
             if is_casl:
+                tau = metadata["LabelingDuration"]
                 denom_factor = t1blood * (1 - np.exp(-(tau / t1blood)))
 
             elif not metadata["BolusCutOffFlag"]:

aslprep/tests/data/expected_outputs_examples_pcasl_multipld.txt

@@ -149,6 +149,10 @@ sub-01/perf/sub-01_atlas-Tian_desc-basil_cbf.tsv
 sub-01/perf/sub-01_atlas-Tian_desc-basil_coverage.tsv
 sub-01/perf/sub-01_att.json
 sub-01/perf/sub-01_att.nii.gz
+sub-01/perf/sub-01_abat.json
+sub-01/perf/sub-01_abat.nii.gz
+sub-01/perf/sub-01_abv.json
+sub-01/perf/sub-01_abv.nii.gz
 sub-01/perf/sub-01_cbf.json
 sub-01/perf/sub-01_cbf.nii.gz
 sub-01/perf/sub-01_desc-basilGM_cbf.json

aslprep/tests/test_cli.py

@@ -54,6 +54,7 @@ def test_examples_pasl_multipld(data_dir, output_dir, working_dir):
         f"{os.path.join(data_dir, 'anatomical/smriprep')}",
     ]
 
+    # _run_and_generate(TEST_NAME, PARTICIPANT_LABEL, parameters, out_dir)
     _run_and_fail(parameters)