From 4df67f3141cbb47c86ac54f4fe3892ee13ed3410 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Thu, 14 Dec 2023 15:58:46 -0500 Subject: [PATCH 01/73] Draft model fitting code for GKM. --- aslprep/utils/cbf.py | 118 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 118 insertions(+) diff --git a/aslprep/utils/cbf.py b/aslprep/utils/cbf.py index ec45f7ff8..21735fcc4 100644 --- a/aslprep/utils/cbf.py +++ b/aslprep/utils/cbf.py @@ -763,6 +763,124 @@ def estimate_cbf_pcasl_multipld( return att_arr, cbf +def calculate_deltam(cbf, att, labeleff, alpha_bs, t1blood, m0_a, m0_b, ld, pld, abat, a_bv): + """Specify a model for use with scipy curve fitting. + + The model is fit separately for each voxel. + + Parameters + ---------- + labeleff + Labeling efficiency. + Used for both deltam_tiss and deltam_art calculation. + alpha_bs + Total background suppression inversion efficiency. + Used for both deltam_tiss and deltam_art calculation. + t1blood + Longitudinal relaxation time of arterial blood in seconds. + Used for both deltam_tiss and deltam_art calculation. + m0_a + Equilibrium magnetization of arterial blood, calculated as M0a = SIpd / lambda, + where SIpd is a proton density weighted image and lambda is the tissue/blood partition + coefficient. + Used for deltam_tiss calculation. + m0_b + ??? + Used for deltam_art calculation. + cbf + Cerebral blood flow. + Used for deltam_tiss calculation. + att + Arterial transit time. + Used for deltam_tiss calculation. + ld + Labeling duration. + Used for both deltam_tiss and deltam_art calculation. + pld + Post-labeling delay. + Used for both deltam_tiss and deltam_art calculation. + abat + The arterial bolus arrival time, describing the first arrival of the labeled blood water + in the voxel within the arterial vessel. + Used for deltam_art calculation. + a_bv + Arterial blood volume. The fraction of the voxel occupied by the arterial vessel. + Used for deltam_art calculation. + + Returns + ------- + deltam + """ + if (ld + pld) < att: + # 0 < LD + PLD < ATT + deltam_tiss = 0 + elif att < (ld + pld) < (att + ld): + # ATT < LD + PLD < ATT + LD + deltam_tiss = ( + (2 * labeleff * alpha_bs * t1blood * m0_a * cbf * (np.exp ** (-att / t1blood))) + * (1 - (np.exp ** (-(ld + pld - att) / t1blood))) + / 6000 + ) + else: + # ATT < PLD + deltam_tiss = ( + (2 * labeleff * alpha_bs * t1blood * m0_a * cbf * (np.exp ** (-pld / t1blood))) + * (1 - (np.exp ** (-ld / t1blood))) + / 6000 + ) + + # Intravascular component + if (ld + pld) < abat: + # 0 < LD + PLD < aBAT + deltam_art = 0 + elif abat < (ld + pld) < (abat + ld): + # aBAT < LD + PLD < aBAT + LD + deltam_art = 2 * labeleff * alpha_bs * m0_b * a_bv * (np.exp ** (-abat / t1blood)) + else: + # aBAT < PLD + deltam_art = 0 + + deltam = deltam_tiss + deltam_art + return deltam + + +def fit_deltam(deltam, plds, lds, t1blood, labeleff, scaled_m0data, partition_coefficient): + import scipy + + assert deltam.shape == plds.shape + assert scaled_m0data.shape == deltam.shape[1] + + m0_a = scaled_m0data / partition_coefficient + + # alpha_bs, m0_b, cbf, att, abat, a_bv + + n_voxels = deltam.shape[1] + csf = np.zeros(n_voxels) + att = np.zeros(n_voxels) + for i_voxel in range(n_voxels): + deltam_voxel = deltam[:, i_voxel] + plds_voxel = plds[:, i_voxel] + + popt, cov = scipy.optimize.curve_fit( + calculate_deltam, + deltam_voxel, + bounds=((0, 0), (np.inf, np.inf)), + ld=lds, + pld=plds_voxel, + t1blood=t1blood, + labeleff=labeleff, + m0_a=m0_a, + alpha_bs=1, # TODO: Figure out what this should be. + m0_b=m0_a, # TODO: Figure out what this means and what it should be. + abat=1, # TODO: Figure out what this should be. + a_bv=1, # TODO: Figure out what this should be. + ) + csf[i_voxel] = popt[0] + att[i_voxel] = popt[1] + + return csf, att + + def estimate_t1(metadata): """Estimate the relaxation rates of blood and gray matter based on magnetic field strength. From 1b5357d87cd6444704adbd6eb403071393cab9bb Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Thu, 21 Dec 2023 12:35:02 -0500 Subject: [PATCH 02/73] Keep drafting the fitting code. --- aslprep/utils/cbf.py | 125 ++++++++++++++++++++++++++++--------------- 1 file changed, 81 insertions(+), 44 deletions(-) diff --git a/aslprep/utils/cbf.py b/aslprep/utils/cbf.py index 21735fcc4..e07cf22da 100644 --- a/aslprep/utils/cbf.py +++ b/aslprep/utils/cbf.py @@ -763,54 +763,69 @@ def estimate_cbf_pcasl_multipld( return att_arr, cbf -def calculate_deltam(cbf, att, labeleff, alpha_bs, t1blood, m0_a, m0_b, ld, pld, abat, a_bv): +def calculate_deltam(X, cbf, att, abat, abv): """Specify a model for use with scipy curve fitting. The model is fit separately for each voxel. + This model is used to estimate ``cbf``, ``att``, ``abat``, and ``abv``. + Parameters ---------- - labeleff + cbf : :obj:`float` + Cerebral blood flow. + Used for deltam_tiss calculation. + Estimated by the model. + att : :obj:`float` + Arterial transit time. + Used for deltam_tiss calculation. + Estimated by the model. + abat : :obj:`float` + The arterial bolus arrival time, describing the first arrival of the labeled blood + water in the voxel within the arterial vessel. Used for deltam_art calculation. + Estimated by the model. + abv : :obj:`float` + Arterial blood volume. + The fraction of the voxel occupied by the arterial vessel. + Used for deltam_art calculation. + Estimated by the model. + X : :obj:`tuple` of length 7 + Dependent variables. + + Returns + ------- + deltam + + Notes + ----- + X boils down into seven variables: + + labeleff : :obj:`float` Labeling efficiency. Used for both deltam_tiss and deltam_art calculation. - alpha_bs + alpha_bs : :obj:`float` Total background suppression inversion efficiency. Used for both deltam_tiss and deltam_art calculation. - t1blood + t1blood : :obj:`float` Longitudinal relaxation time of arterial blood in seconds. Used for both deltam_tiss and deltam_art calculation. - m0_a - Equilibrium magnetization of arterial blood, calculated as M0a = SIpd / lambda, + m0_a : :obj:`float` + Equilibrium magnetization of tissue, calculated as M0a = SIpd / lambda, where SIpd is a proton density weighted image and lambda is the tissue/blood partition coefficient. Used for deltam_tiss calculation. - m0_b - ??? + m0_b : :obj:`float` + Equilibrium magnetization of arterial blood. Used for deltam_art calculation. - cbf - Cerebral blood flow. - Used for deltam_tiss calculation. - att - Arterial transit time. - Used for deltam_tiss calculation. - ld - Labeling duration. + ld : :obj:`numpy.ndarray` + Labeling durations. Used for both deltam_tiss and deltam_art calculation. - pld - Post-labeling delay. + pld : :obj:`numpy.ndarray` + Post-labeling delays. Used for both deltam_tiss and deltam_art calculation. - abat - The arterial bolus arrival time, describing the first arrival of the labeled blood water - in the voxel within the arterial vessel. - Used for deltam_art calculation. - a_bv - Arterial blood volume. The fraction of the voxel occupied by the arterial vessel. - Used for deltam_art calculation. - - Returns - ------- - deltam """ + labeleff, alpha_bs, t1blood, m0_a, m0_b, ld, pld = X + if (ld + pld) < att: # 0 < LD + PLD < ATT deltam_tiss = 0 @@ -835,7 +850,7 @@ def calculate_deltam(cbf, att, labeleff, alpha_bs, t1blood, m0_a, m0_b, ld, pld, deltam_art = 0 elif abat < (ld + pld) < (abat + ld): # aBAT < LD + PLD < aBAT + LD - deltam_art = 2 * labeleff * alpha_bs * m0_b * a_bv * (np.exp ** (-abat / t1blood)) + deltam_art = 2 * labeleff * alpha_bs * m0_b * abv * (np.exp ** (-abat / t1blood)) else: # aBAT < PLD deltam_art = 0 @@ -845,6 +860,25 @@ def calculate_deltam(cbf, att, labeleff, alpha_bs, t1blood, m0_a, m0_b, ld, pld, def fit_deltam(deltam, plds, lds, t1blood, labeleff, scaled_m0data, partition_coefficient): + """Estimate CBF, ATT, aBAT, and aBV from multi-PLD data. + + Parameters + ---------- + deltam + plds + lds + t1blood + labeleff + scaled_m0data + partition_coefficient + + Returns + ------- + csf + att + abat + abv + """ import scipy assert deltam.shape == plds.shape @@ -852,33 +886,36 @@ def fit_deltam(deltam, plds, lds, t1blood, labeleff, scaled_m0data, partition_co m0_a = scaled_m0data / partition_coefficient - # alpha_bs, m0_b, cbf, att, abat, a_bv + # alpha_bs, m0_b, cbf, att, abat, abv n_voxels = deltam.shape[1] csf = np.zeros(n_voxels) att = np.zeros(n_voxels) + abat = np.zeros(n_voxels) + abv = np.zeros(n_voxels) for i_voxel in range(n_voxels): deltam_voxel = deltam[:, i_voxel] plds_voxel = plds[:, i_voxel] - popt, cov = scipy.optimize.curve_fit( + # TODO: Figure out what alpha_bs and m0_b should be. + popt = scipy.optimize.curve_fit( calculate_deltam, deltam_voxel, - bounds=((0, 0), (np.inf, np.inf)), - ld=lds, - pld=plds_voxel, - t1blood=t1blood, - labeleff=labeleff, - m0_a=m0_a, - alpha_bs=1, # TODO: Figure out what this should be. - m0_b=m0_a, # TODO: Figure out what this means and what it should be. - abat=1, # TODO: Figure out what this should be. - a_bv=1, # TODO: Figure out what this should be. - ) + X=(labeleff, labeleff, t1blood, m0_a, m0_a, lds, plds_voxel), + # initial estimates for DVs + csf=1, + att=1, + abat=1, + abv=1, + # lower and upper bounds for DVs + bounds=((0, 0, 0, 0), (np.inf, np.inf, np.inf, np.inf)), + )[0] csf[i_voxel] = popt[0] att[i_voxel] = popt[1] + abat[i_voxel] = popt[2] + abv[i_voxel] = popt[3] - return csf, att + return csf, att, abat, abv def estimate_t1(metadata): From 8fdb90d6a6765e9a47268166d936ad582c0ad4bf Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Thu, 21 Dec 2023 12:40:03 -0500 Subject: [PATCH 03/73] It's CBF, not CSF. --- aslprep/utils/cbf.py | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/aslprep/utils/cbf.py b/aslprep/utils/cbf.py index e07cf22da..e576227ff 100644 --- a/aslprep/utils/cbf.py +++ b/aslprep/utils/cbf.py @@ -794,11 +794,12 @@ def calculate_deltam(X, cbf, att, abat, abv): Returns ------- - deltam + deltam : :obj:`float` + Delta-M value for the voxel. Notes ----- - X boils down into seven variables: + X breaks down into seven variables: labeleff : :obj:`float` Labeling efficiency. @@ -874,7 +875,7 @@ def fit_deltam(deltam, plds, lds, t1blood, labeleff, scaled_m0data, partition_co Returns ------- - csf + cbf att abat abv @@ -889,7 +890,7 @@ def fit_deltam(deltam, plds, lds, t1blood, labeleff, scaled_m0data, partition_co # alpha_bs, m0_b, cbf, att, abat, abv n_voxels = deltam.shape[1] - csf = np.zeros(n_voxels) + cbf = np.zeros(n_voxels) att = np.zeros(n_voxels) abat = np.zeros(n_voxels) abv = np.zeros(n_voxels) @@ -903,19 +904,19 @@ def fit_deltam(deltam, plds, lds, t1blood, labeleff, scaled_m0data, partition_co deltam_voxel, X=(labeleff, labeleff, t1blood, m0_a, m0_a, lds, plds_voxel), # initial estimates for DVs - csf=1, + cbf=1, att=1, abat=1, abv=1, # lower and upper bounds for DVs bounds=((0, 0, 0, 0), (np.inf, np.inf, np.inf, np.inf)), )[0] - csf[i_voxel] = popt[0] + cbf[i_voxel] = popt[0] att[i_voxel] = popt[1] abat[i_voxel] = popt[2] abv[i_voxel] = popt[3] - return csf, att, abat, abv + return cbf, att, abat, abv def estimate_t1(metadata): From 462d76e690e937843cb9989de04a9ac3a363e251 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Thu, 21 Dec 2023 12:40:48 -0500 Subject: [PATCH 04/73] Move the functions. --- aslprep/utils/cbf.py | 116 +++++++++++++++++++++---------------------- 1 file changed, 58 insertions(+), 58 deletions(-) diff --git a/aslprep/utils/cbf.py b/aslprep/utils/cbf.py index e576227ff..7e29f1429 100644 --- a/aslprep/utils/cbf.py +++ b/aslprep/utils/cbf.py @@ -763,6 +763,64 @@ def estimate_cbf_pcasl_multipld( return att_arr, cbf +def estimate_t1(metadata): + """Estimate the relaxation rates of blood and gray matter based on magnetic field strength. + + t1blood is set based on the scanner's field strength, + according to :footcite:t:`zhang2013vivo,alsop_recommended_2015`. + If recommended values from these publications cannot be used + (i.e., if the field strength isn't 1.5T, 3T, 7T), + then the formula from :footcite:t:`zhang2013vivo` will be applied. + + t1tissue is set based on the scanner's field strength as well, + according to :footcite:t:`wright2008water`. + At the moment, field strengths other than 1.5T, 3T, and 7T are not supported and will + raise an exception. + + Parameters + ---------- + metadata : :obj:`dict` + Dictionary of metadata from the ASL file. + + Returns + ------- + t1blood : :obj:`float` + Estimated relaxation rate of blood based on magnetic field strength. + t1tissue : :obj:`float` + Estimated relaxation rate of gray matter based on magnetic field strength. + + References + ---------- + .. footbibliography:: + """ + T1BLOOD_DICT = { + 1.5: 1.35, + 3: 1.65, + 7: 2.087, + } + t1blood = T1BLOOD_DICT.get(metadata["MagneticFieldStrength"]) + if not t1blood: + config.loggers.interface.warning( + f"T1blood cannot be inferred for {metadata['MagneticFieldStrength']}T data. " + "Defaulting to formula from Zhang et al. (2013)." + ) + t1blood = (110 * metadata["MagneticFieldStrength"] + 1316) / 1000 + + # TODO: Supplement with formula for other field strengths + T1TISSUE_DICT = { + 1.5: 1.197, + 3: 1.607, + 7: 1.939, + } + t1tissue = T1TISSUE_DICT.get(metadata["MagneticFieldStrength"]) + if not t1tissue: + raise ValueError( + f"T1tissue cannot be inferred for {metadata['MagneticFieldStrength']}T data." + ) + + return t1blood, t1tissue + + def calculate_deltam(X, cbf, att, abat, abv): """Specify a model for use with scipy curve fitting. @@ -917,61 +975,3 @@ def fit_deltam(deltam, plds, lds, t1blood, labeleff, scaled_m0data, partition_co abv[i_voxel] = popt[3] return cbf, att, abat, abv - - -def estimate_t1(metadata): - """Estimate the relaxation rates of blood and gray matter based on magnetic field strength. - - t1blood is set based on the scanner's field strength, - according to :footcite:t:`zhang2013vivo,alsop_recommended_2015`. - If recommended values from these publications cannot be used - (i.e., if the field strength isn't 1.5T, 3T, 7T), - then the formula from :footcite:t:`zhang2013vivo` will be applied. - - t1tissue is set based on the scanner's field strength as well, - according to :footcite:t:`wright2008water`. - At the moment, field strengths other than 1.5T, 3T, and 7T are not supported and will - raise an exception. - - Parameters - ---------- - metadata : :obj:`dict` - Dictionary of metadata from the ASL file. - - Returns - ------- - t1blood : :obj:`float` - Estimated relaxation rate of blood based on magnetic field strength. - t1tissue : :obj:`float` - Estimated relaxation rate of gray matter based on magnetic field strength. - - References - ---------- - .. footbibliography:: - """ - T1BLOOD_DICT = { - 1.5: 1.35, - 3: 1.65, - 7: 2.087, - } - t1blood = T1BLOOD_DICT.get(metadata["MagneticFieldStrength"]) - if not t1blood: - config.loggers.interface.warning( - f"T1blood cannot be inferred for {metadata['MagneticFieldStrength']}T data. " - "Defaulting to formula from Zhang et al. (2013)." - ) - t1blood = (110 * metadata["MagneticFieldStrength"] + 1316) / 1000 - - # TODO: Supplement with formula for other field strengths - T1TISSUE_DICT = { - 1.5: 1.197, - 3: 1.607, - 7: 1.939, - } - t1tissue = T1TISSUE_DICT.get(metadata["MagneticFieldStrength"]) - if not t1tissue: - raise ValueError( - f"T1tissue cannot be inferred for {metadata['MagneticFieldStrength']}T data." - ) - - return t1blood, t1tissue From 77d4f1bdc191c78f9873f3db1debae5bf7131f04 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Thu, 21 Dec 2023 13:09:59 -0500 Subject: [PATCH 05/73] Try using the new method. --- aslprep/data/boilerplate.bib | 7 ++++ aslprep/interfaces/cbf.py | 32 ++++++++++++++++--- aslprep/utils/cbf.py | 62 ++++++++++++++++++++++++------------ aslprep/workflows/asl/cbf.py | 4 +++ 4 files changed, 79 insertions(+), 26 deletions(-) diff --git a/aslprep/data/boilerplate.bib b/aslprep/data/boilerplate.bib index 1acdf6503..3967b123a 100644 --- a/aslprep/data/boilerplate.bib +++ b/aslprep/data/boilerplate.bib @@ -852,3 +852,10 @@ @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 others}, + year={2023}, + publisher={OSF Preprints} +} diff --git a/aslprep/interfaces/cbf.py b/aslprep/interfaces/cbf.py index bd998f4fe..1e8527c70 100644 --- a/aslprep/interfaces/cbf.py +++ b/aslprep/interfaces/cbf.py @@ -27,8 +27,8 @@ from aslprep.utils.cbf import ( _getcbfscore, _scrubcbf, - estimate_cbf_pcasl_multipld, estimate_t1, + fit_deltam_pcasl, ) @@ -342,6 +342,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, @@ -508,15 +518,13 @@ def _run_interface(self, runtime): if is_multi_pld: if is_casl: - att, mean_cbf = estimate_cbf_pcasl_multipld( + cbf, att, abat, abv = fit_deltam_pcasl( deltam_arr, scaled_m0data, plds, tau, labeleff, t1blood=t1blood, - t1tissue=t1tissue, - unit_conversion=UNIT_CONV, partition_coefficient=PARTITION_COEF, ) @@ -526,8 +534,10 @@ def _run_interface(self, runtime): "Multi-delay data are not supported for PASL sequences at the moment." ) - mean_cbf_img = masker.inverse_transform(mean_cbf) + 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,6 +547,18 @@ 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: diff --git a/aslprep/utils/cbf.py b/aslprep/utils/cbf.py index 7e29f1429..312f6512f 100644 --- a/aslprep/utils/cbf.py +++ b/aslprep/utils/cbf.py @@ -821,12 +821,13 @@ def estimate_t1(metadata): return t1blood, t1tissue -def calculate_deltam(X, cbf, att, abat, abv): +def calculate_deltam_pcasl(X, cbf, att, abat, abv): """Specify a model for use with scipy curve fitting. The model is fit separately for each voxel. - This model is used to estimate ``cbf``, ``att``, ``abat``, and ``abv``. + This model is used to estimate ``cbf``, ``att``, ``abat``, and ``abv`` for multi-PLD PCASL + data. Parameters ---------- @@ -876,38 +877,38 @@ def calculate_deltam(X, cbf, att, abat, abv): m0_b : :obj:`float` Equilibrium magnetization of arterial blood. Used for deltam_art calculation. - ld : :obj:`numpy.ndarray` + tau : :obj:`numpy.ndarray` Labeling durations. Used for both deltam_tiss and deltam_art calculation. pld : :obj:`numpy.ndarray` Post-labeling delays. Used for both deltam_tiss and deltam_art calculation. """ - labeleff, alpha_bs, t1blood, m0_a, m0_b, ld, pld = X + labeleff, alpha_bs, t1blood, m0_a, m0_b, tau, pld = X - if (ld + pld) < att: + if (tau + pld) < att: # 0 < LD + PLD < ATT deltam_tiss = 0 - elif att < (ld + pld) < (att + ld): + elif att < (tau + pld) < (att + tau): # ATT < LD + PLD < ATT + LD deltam_tiss = ( (2 * labeleff * alpha_bs * t1blood * m0_a * cbf * (np.exp ** (-att / t1blood))) - * (1 - (np.exp ** (-(ld + pld - att) / t1blood))) + * (1 - (np.exp ** (-(tau + pld - att) / t1blood))) / 6000 ) else: # ATT < PLD deltam_tiss = ( (2 * labeleff * alpha_bs * t1blood * m0_a * cbf * (np.exp ** (-pld / t1blood))) - * (1 - (np.exp ** (-ld / t1blood))) + * (1 - (np.exp ** (-tau / t1blood))) / 6000 ) # Intravascular component - if (ld + pld) < abat: + if (tau + pld) < abat: # 0 < LD + PLD < aBAT deltam_art = 0 - elif abat < (ld + pld) < (abat + ld): + elif abat < (tau + pld) < (abat + tau): # aBAT < LD + PLD < aBAT + LD deltam_art = 2 * labeleff * alpha_bs * m0_b * abv * (np.exp ** (-abat / t1blood)) else: @@ -918,14 +919,26 @@ def calculate_deltam(X, cbf, att, abat, abv): return deltam -def fit_deltam(deltam, plds, lds, t1blood, labeleff, scaled_m0data, partition_coefficient): - """Estimate CBF, ATT, aBAT, and aBV from multi-PLD data. +def fit_deltam_pcasl( + deltam_arr, + scaled_m0data, + plds, + tau, + labeleff, + t1blood, + partition_coefficient, +): + """Estimate CBF, ATT, aBAT, and aBV from multi-PLD PCASL data. + + This function uses the general kinetic model specified by + :footcite:t:`woods2023recommendations`. Parameters ---------- deltam plds - lds + tau + Label duration. May be a single value or may vary across volumes/PLDs. t1blood labeleff scaled_m0data @@ -937,30 +950,37 @@ def fit_deltam(deltam, plds, lds, t1blood, labeleff, scaled_m0data, partition_co att abat abv + + Notes + ----- + This model does not include a separate term for T1,tissue, but we could expand it to do so + in the future. + + References + ---------- + .. footbibliography:: """ import scipy - assert deltam.shape == plds.shape - assert scaled_m0data.shape == deltam.shape[1] + assert deltam_arr.shape == plds.shape + assert scaled_m0data.shape == deltam_arr.shape[1] m0_a = scaled_m0data / partition_coefficient - # alpha_bs, m0_b, cbf, att, abat, abv - - n_voxels = deltam.shape[1] + n_voxels = deltam_arr.shape[1] cbf = np.zeros(n_voxels) att = np.zeros(n_voxels) abat = np.zeros(n_voxels) abv = np.zeros(n_voxels) for i_voxel in range(n_voxels): - deltam_voxel = deltam[:, i_voxel] + deltam_voxel = deltam_arr[:, i_voxel] plds_voxel = plds[:, i_voxel] # TODO: Figure out what alpha_bs and m0_b should be. popt = scipy.optimize.curve_fit( - calculate_deltam, + calculate_deltam_pcasl, deltam_voxel, - X=(labeleff, labeleff, t1blood, m0_a, m0_a, lds, plds_voxel), + X=(labeleff, labeleff, t1blood, m0_a, m0_a, tau, plds_voxel), # initial estimates for DVs cbf=1, att=1, diff --git a/aslprep/workflows/asl/cbf.py b/aslprep/workflows/asl/cbf.py index 648ddcaaf..2c54ffb17 100644 --- a/aslprep/workflows/asl/cbf.py +++ b/aslprep/workflows/asl/cbf.py @@ -240,6 +240,8 @@ def init_cbf_wf( "mean_cbf", "cbf_ts", # Only calculated for single-delay data "att", # Only calculated for multi-delay data + "abat", # Only calculated for multi-delay data + "abv", # Only calculated for multi-delay data "plds", # SCORE/SCRUB outputs "cbf_ts_score", @@ -425,6 +427,8 @@ def _getfiledir(file): ("mean_cbf", "mean_cbf"), ("att", "att"), ("plds", "plds"), + ("abat", "abat"), + ("abv", "abv"), ]), ]) # fmt:skip From 29889858b1ac1301eba2477561614321b38b045f Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Thu, 21 Dec 2023 13:29:33 -0500 Subject: [PATCH 06/73] Update boilerplate.bib --- aslprep/data/boilerplate.bib | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/aslprep/data/boilerplate.bib b/aslprep/data/boilerplate.bib index 3967b123a..7d95da53a 100644 --- a/aslprep/data/boilerplate.bib +++ b/aslprep/data/boilerplate.bib @@ -855,7 +855,10 @@ @article{topup @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 others}, + 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}, - publisher={OSF Preprints} + journal={OSF Preprints}, + publisher={OSF Preprints}, + doi={10.31219/osf.io/4tskr}, + url={https://doi.org/10.31219/osf.io/4tskr}, } From 2f038ef3867d64cdd04ece766c3d56e2a8ace9f0 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Thu, 21 Dec 2023 13:31:52 -0500 Subject: [PATCH 07/73] Update cbf.py --- aslprep/interfaces/cbf.py | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/aslprep/interfaces/cbf.py b/aslprep/interfaces/cbf.py index 1e8527c70..8d57d28d1 100644 --- a/aslprep/interfaces/cbf.py +++ b/aslprep/interfaces/cbf.py @@ -24,12 +24,7 @@ estimate_labeling_efficiency, pcasl_or_pasl, ) -from aslprep.utils.cbf import ( - _getcbfscore, - _scrubcbf, - estimate_t1, - fit_deltam_pcasl, -) +from aslprep.utils.cbf import _getcbfscore, _scrubcbf, estimate_t1, fit_deltam_pcasl class _RefineMaskInputSpec(BaseInterfaceInputSpec): From c4046f13593eed173a23c955c06a15d34864a653 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Thu, 21 Dec 2023 13:56:12 -0500 Subject: [PATCH 08/73] Improve shape checks. --- aslprep/utils/cbf.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/aslprep/utils/cbf.py b/aslprep/utils/cbf.py index 312f6512f..22bb4f7f3 100644 --- a/aslprep/utils/cbf.py +++ b/aslprep/utils/cbf.py @@ -962,8 +962,13 @@ def fit_deltam_pcasl( """ import scipy - assert deltam_arr.shape == plds.shape - assert scaled_m0data.shape == deltam_arr.shape[1] + if deltam_arr.shape != plds.shape: + raise ValueError(f"deltam ({deltam_arr.shape}) != plds ({plds.shape})") + + if scaled_m0data.shape != deltam_arr.shape[1]: + raise ValueError( + f"scaled_m0data ({scaled_m0data.shape}) != deltam_arr ({deltam_arr.shape[1]})" + ) m0_a = scaled_m0data / partition_coefficient From 42125b9d20d80fe7af49c9830aee52af5a78ae0b Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Thu, 21 Dec 2023 14:22:12 -0500 Subject: [PATCH 09/73] Update cbf.py --- aslprep/utils/cbf.py | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/aslprep/utils/cbf.py b/aslprep/utils/cbf.py index 22bb4f7f3..1cc63c3c3 100644 --- a/aslprep/utils/cbf.py +++ b/aslprep/utils/cbf.py @@ -962,24 +962,27 @@ def fit_deltam_pcasl( """ import scipy + n_voxels, n_volumes = deltam_arr.shape if deltam_arr.shape != plds.shape: raise ValueError(f"deltam ({deltam_arr.shape}) != plds ({plds.shape})") - if scaled_m0data.shape != deltam_arr.shape[1]: + if scaled_m0data.size != n_voxels: raise ValueError( - f"scaled_m0data ({scaled_m0data.shape}) != deltam_arr ({deltam_arr.shape[1]})" + f"scaled_m0data ({scaled_m0data.size}) != deltam_arr ({n_voxels})" ) + if tau.size != n_volumes: + raise ValueError(f"tau ({tau.shape}) != n_volumes {n_volumes}") + m0_a = scaled_m0data / partition_coefficient - n_voxels = deltam_arr.shape[1] cbf = np.zeros(n_voxels) att = np.zeros(n_voxels) abat = np.zeros(n_voxels) abv = np.zeros(n_voxels) for i_voxel in range(n_voxels): - deltam_voxel = deltam_arr[:, i_voxel] - plds_voxel = plds[:, i_voxel] + deltam_voxel = deltam_arr[i_voxel, :] + plds_voxel = plds[i_voxel, :] # TODO: Figure out what alpha_bs and m0_b should be. popt = scipy.optimize.curve_fit( From 45a621332a71e07c7e5745c30c29981da3160ad8 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Thu, 21 Dec 2023 14:23:14 -0500 Subject: [PATCH 10/73] Update cbf.py --- aslprep/utils/cbf.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/aslprep/utils/cbf.py b/aslprep/utils/cbf.py index 1cc63c3c3..834733a43 100644 --- a/aslprep/utils/cbf.py +++ b/aslprep/utils/cbf.py @@ -967,9 +967,7 @@ def fit_deltam_pcasl( raise ValueError(f"deltam ({deltam_arr.shape}) != plds ({plds.shape})") if scaled_m0data.size != n_voxels: - raise ValueError( - f"scaled_m0data ({scaled_m0data.size}) != deltam_arr ({n_voxels})" - ) + raise ValueError(f"scaled_m0data ({scaled_m0data.size}) != deltam_arr ({n_voxels})") if tau.size != n_volumes: raise ValueError(f"tau ({tau.shape}) != n_volumes {n_volumes}") From 23b48ca8e5bebfea42940fdf92fe424b11c72d7a Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Thu, 21 Dec 2023 14:42:21 -0500 Subject: [PATCH 11/73] Try adding equations to docs. --- docs/workflows.rst | 31 +++++++++++++++++++++++++++++++ 1 file changed, 31 insertions(+) diff --git a/docs/workflows.rst b/docs/workflows.rst index aab72b1d8..1f6428c76 100644 --- a/docs/workflows.rst +++ b/docs/workflows.rst @@ -438,6 +438,37 @@ Also, multi-delay ASL allows for the estimation of arterial transit time (ATT). Pseudo-Continuous ASL --------------------- +For multi-delay PCASL data, ASLPrep uses :func:`scipy.optimize.curvefit` to estimate +CBF, ATT, aBAT, and aBV values for each voxel using Equations 2 and 4 in +:footcite:t:`woods2023recommendations`. + +Equation 2: + +.. math:: + + \Delta{M} = 0 + + \Delta{M} = \frac{ 2 \cdot \alpha \cdot \alpha_{BS} \cdot T_{1b} \cdot M_{0a} \cdot CBF \cdot + e ^ -\frac{ ATT } { T_{1,blood} } \cdot \left[ + 1 - e ^ {-\frac{ LD + PLD - ATT } { T_{1,tissue} }} + \right] }{ 6000 } + + \Delta{M} = \frac{ 2 \cdot \alpha \cdot \alpha_{BS} \cdot T_{1b} \cdot M_{0a} \cdot CBF \cdot + e ^ -\frac{ PLD } { T_{1,blood} } \cdot \left[ + 1 - e ^ {-\frac{ LD } { T_{1,tissue} }} + \right] }{ 6000 } + +Equation 4: + +.. math:: + + \Delta{M}_{art} = 0 + + \Delta{M}_{art} = 2 \cdot \alpha \cdot \alpha_{BS} \cdot M_{0b} \cdot aBV \cdot + e ^ -\frac{ aBAT } { T_{1,blood} } + + \Delta{M}_{art} = 0 + For multi-delay PCASL data, the following steps are taken: 1. :math:`\Delta{M}` values are first averaged over time for each unique post-labeling delay value. From 1cb0706876bbf57c807fe41756bf0753a70accba Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Thu, 21 Dec 2023 16:45:56 -0500 Subject: [PATCH 12/73] Make tau an array. --- aslprep/interfaces/cbf.py | 2 +- aslprep/utils/cbf.py | 11 ++++++++--- docs/workflows.rst | 6 +++--- 3 files changed, 12 insertions(+), 7 deletions(-) diff --git a/aslprep/interfaces/cbf.py b/aslprep/interfaces/cbf.py index 8d57d28d1..ffd89e00f 100644 --- a/aslprep/interfaces/cbf.py +++ b/aslprep/interfaces/cbf.py @@ -509,7 +509,7 @@ def _run_interface(self, runtime): plds = np.dot(plds[:, None], np.ones((1, deltam_arr.shape[0]))).T if is_casl: - tau = np.array(metadata["LabelingDuration"]) + tau = metadata["LabelingDuration"] if is_multi_pld: if is_casl: diff --git a/aslprep/utils/cbf.py b/aslprep/utils/cbf.py index 834733a43..99ce51136 100644 --- a/aslprep/utils/cbf.py +++ b/aslprep/utils/cbf.py @@ -937,7 +937,7 @@ def fit_deltam_pcasl( ---------- deltam plds - tau + tau : :obj:`float` or :obj:`numpy.ndarray` Label duration. May be a single value or may vary across volumes/PLDs. t1blood labeleff @@ -969,8 +969,13 @@ def fit_deltam_pcasl( if scaled_m0data.size != n_voxels: raise ValueError(f"scaled_m0data ({scaled_m0data.size}) != deltam_arr ({n_voxels})") - if tau.size != n_volumes: - raise ValueError(f"tau ({tau.shape}) != n_volumes {n_volumes}") + if isinstance(tau, float): + tau = np.full((n_volumes,), tau) + else: + tau = np.ndarray(tau) + + if tau.size != n_volumes: + raise ValueError(f"tau ({tau.shape}) != n_volumes {n_volumes}") m0_a = scaled_m0data / partition_coefficient diff --git a/docs/workflows.rst b/docs/workflows.rst index 1f6428c76..d0a6e58cb 100644 --- a/docs/workflows.rst +++ b/docs/workflows.rst @@ -446,14 +446,14 @@ Equation 2: .. math:: - \Delta{M} = 0 + &= 0 - \Delta{M} = \frac{ 2 \cdot \alpha \cdot \alpha_{BS} \cdot T_{1b} \cdot M_{0a} \cdot CBF \cdot + \Delta{M} &= \frac{ 2 \cdot \alpha \cdot \alpha_{BS} \cdot T_{1b} \cdot M_{0a} \cdot CBF \cdot e ^ -\frac{ ATT } { T_{1,blood} } \cdot \left[ 1 - e ^ {-\frac{ LD + PLD - ATT } { T_{1,tissue} }} \right] }{ 6000 } - \Delta{M} = \frac{ 2 \cdot \alpha \cdot \alpha_{BS} \cdot T_{1b} \cdot M_{0a} \cdot CBF \cdot + &= \frac{ 2 \cdot \alpha \cdot \alpha_{BS} \cdot T_{1b} \cdot M_{0a} \cdot CBF \cdot e ^ -\frac{ PLD } { T_{1,blood} } \cdot \left[ 1 - e ^ {-\frac{ LD } { T_{1,tissue} }} \right] }{ 6000 } From 7c9905a609bcb7ed898b3793a75c0a9e96208b9e Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Thu, 21 Dec 2023 16:57:06 -0500 Subject: [PATCH 13/73] Update workflows.rst --- docs/workflows.rst | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/docs/workflows.rst b/docs/workflows.rst index d0a6e58cb..cf8d86947 100644 --- a/docs/workflows.rst +++ b/docs/workflows.rst @@ -446,28 +446,28 @@ Equation 2: .. math:: - &= 0 + &= 0 \hfill \text{(if } 0 < LD + PLD < ATT) - \Delta{M} &= \frac{ 2 \cdot \alpha \cdot \alpha_{BS} \cdot T_{1b} \cdot M_{0a} \cdot CBF \cdot + \Delta{M}_tiss &= \frac{ 2 \cdot \alpha \cdot \alpha_{BS} \cdot T_{1b} \cdot M_{0a} \cdot CBF \cdot e ^ -\frac{ ATT } { T_{1,blood} } \cdot \left[ 1 - e ^ {-\frac{ LD + PLD - ATT } { T_{1,tissue} }} - \right] }{ 6000 } + \right] }{ 6000 } \hfill \text{(if } ATT < LD + PLD < ATT + LD) &= \frac{ 2 \cdot \alpha \cdot \alpha_{BS} \cdot T_{1b} \cdot M_{0a} \cdot CBF \cdot e ^ -\frac{ PLD } { T_{1,blood} } \cdot \left[ 1 - e ^ {-\frac{ LD } { T_{1,tissue} }} - \right] }{ 6000 } + \right] }{ 6000 } \hfill \text{(if } ATT < PLD) Equation 4: .. math:: - \Delta{M}_{art} = 0 + &= 0 - \Delta{M}_{art} = 2 \cdot \alpha \cdot \alpha_{BS} \cdot M_{0b} \cdot aBV \cdot + \Delta{M}_{art} &= 2 \cdot \alpha \cdot \alpha_{BS} \cdot M_{0b} \cdot aBV \cdot e ^ -\frac{ aBAT } { T_{1,blood} } - \Delta{M}_{art} = 0 + &= 0 For multi-delay PCASL data, the following steps are taken: From 59a9e785eab0c0b9004b4761a8e1cdfbb8bad63c Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Thu, 21 Dec 2023 17:03:09 -0500 Subject: [PATCH 14/73] Update workflows.rst --- docs/workflows.rst | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/docs/workflows.rst b/docs/workflows.rst index cf8d86947..6d2c27a3c 100644 --- a/docs/workflows.rst +++ b/docs/workflows.rst @@ -446,28 +446,28 @@ Equation 2: .. math:: - &= 0 \hfill \text{(if } 0 < LD + PLD < ATT) + &= 0 \hfill 0 < LD + PLD < ATT \Delta{M}_tiss &= \frac{ 2 \cdot \alpha \cdot \alpha_{BS} \cdot T_{1b} \cdot M_{0a} \cdot CBF \cdot - e ^ -\frac{ ATT } { T_{1,blood} } \cdot \left[ + e ^ { -\frac{ ATT } { T_{1,blood} } } \cdot \left[ 1 - e ^ {-\frac{ LD + PLD - ATT } { T_{1,tissue} }} - \right] }{ 6000 } \hfill \text{(if } ATT < LD + PLD < ATT + LD) + \right] }{ 6000 } \hfill ATT < LD + PLD < ATT + LD &= \frac{ 2 \cdot \alpha \cdot \alpha_{BS} \cdot T_{1b} \cdot M_{0a} \cdot CBF \cdot - e ^ -\frac{ PLD } { T_{1,blood} } \cdot \left[ + e ^ { -\frac{ PLD } { T_{1,blood} } } \cdot \left[ 1 - e ^ {-\frac{ LD } { T_{1,tissue} }} - \right] }{ 6000 } \hfill \text{(if } ATT < PLD) + \right] }{ 6000 } \hfill ATT < PLD Equation 4: .. math:: - &= 0 + &= 0 \hfill 0 < LD + PLD < aBAT \Delta{M}_{art} &= 2 \cdot \alpha \cdot \alpha_{BS} \cdot M_{0b} \cdot aBV \cdot - e ^ -\frac{ aBAT } { T_{1,blood} } + e ^ { -\frac{ aBAT } { T_{1,blood} } } \hfill aBAT < LD + PLD < aBAT + LD - &= 0 + &= 0 \hfill aBAT < PLD For multi-delay PCASL data, the following steps are taken: From a551dd45182cf63712bfa160e9a119b035db4153 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Thu, 21 Dec 2023 17:07:31 -0500 Subject: [PATCH 15/73] Update workflows.rst --- docs/workflows.rst | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/docs/workflows.rst b/docs/workflows.rst index 6d2c27a3c..f9b30587a 100644 --- a/docs/workflows.rst +++ b/docs/workflows.rst @@ -446,28 +446,28 @@ Equation 2: .. math:: - &= 0 \hfill 0 < LD + PLD < ATT + &= 0 \pushright{ 0 < LD + PLD < ATT } \Delta{M}_tiss &= \frac{ 2 \cdot \alpha \cdot \alpha_{BS} \cdot T_{1b} \cdot M_{0a} \cdot CBF \cdot e ^ { -\frac{ ATT } { T_{1,blood} } } \cdot \left[ 1 - e ^ {-\frac{ LD + PLD - ATT } { T_{1,tissue} }} - \right] }{ 6000 } \hfill ATT < LD + PLD < ATT + LD + \right] }{ 6000 } \pushright{ ATT < LD + PLD < ATT + LD } &= \frac{ 2 \cdot \alpha \cdot \alpha_{BS} \cdot T_{1b} \cdot M_{0a} \cdot CBF \cdot e ^ { -\frac{ PLD } { T_{1,blood} } } \cdot \left[ 1 - e ^ {-\frac{ LD } { T_{1,tissue} }} - \right] }{ 6000 } \hfill ATT < PLD + \right] }{ 6000 } \pushright{ ATT < PLD } Equation 4: .. math:: - &= 0 \hfill 0 < LD + PLD < aBAT + &= 0 \pushright{ 0 < LD + PLD < aBAT } \Delta{M}_{art} &= 2 \cdot \alpha \cdot \alpha_{BS} \cdot M_{0b} \cdot aBV \cdot - e ^ { -\frac{ aBAT } { T_{1,blood} } } \hfill aBAT < LD + PLD < aBAT + LD + e ^ { -\frac{ aBAT } { T_{1,blood} } } \pushright{ aBAT < LD + PLD < aBAT + LD } - &= 0 \hfill aBAT < PLD + &= 0 \pushright{ aBAT < PLD } For multi-delay PCASL data, the following steps are taken: From d2566ed140ad44cc05f66b8321e8f34f1cb80e3f Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Thu, 21 Dec 2023 17:12:03 -0500 Subject: [PATCH 16/73] Update workflows.rst --- docs/workflows.rst | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/docs/workflows.rst b/docs/workflows.rst index f9b30587a..3c7e54079 100644 --- a/docs/workflows.rst +++ b/docs/workflows.rst @@ -446,17 +446,17 @@ Equation 2: .. math:: - &= 0 \pushright{ 0 < LD + PLD < ATT } + &= 0 \hfilll { 0 < LD + PLD < ATT } - \Delta{M}_tiss &= \frac{ 2 \cdot \alpha \cdot \alpha_{BS} \cdot T_{1b} \cdot M_{0a} \cdot CBF \cdot - e ^ { -\frac{ ATT } { T_{1,blood} } } \cdot \left[ + \Delta{M}_tiss &= \frac{ 2 \cdot \alpha \cdot \alpha_{BS} \cdot T_{1b} \cdot M_{0a} \cdot CBF + \cdot \left[ 1 - e ^ {-\frac{ LD + PLD - ATT } { T_{1,tissue} }} - \right] }{ 6000 } \pushright{ ATT < LD + PLD < ATT + LD } + \right] }{ 6000 \cdot e ^ { \frac{ ATT } { T_{1,blood} } } } \hfilll { ATT < LD + PLD < ATT + LD } - &= \frac{ 2 \cdot \alpha \cdot \alpha_{BS} \cdot T_{1b} \cdot M_{0a} \cdot CBF \cdot - e ^ { -\frac{ PLD } { T_{1,blood} } } \cdot \left[ + &= \frac{ 2 \cdot \alpha \cdot \alpha_{BS} \cdot T_{1b} \cdot M_{0a} \cdot CBF + \cdot \left[ 1 - e ^ {-\frac{ LD } { T_{1,tissue} }} - \right] }{ 6000 } \pushright{ ATT < PLD } + \right] }{ 6000 e ^ { \frac{ PLD } { T_{1,blood} } } } \hfilll { ATT < PLD } Equation 4: From 15b8c45575e55a6ecf67147a9cadadbb55e57046 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Thu, 21 Dec 2023 18:00:28 -0500 Subject: [PATCH 17/73] Update. --- aslprep/utils/cbf.py | 4 ++-- docs/workflows.rst | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/aslprep/utils/cbf.py b/aslprep/utils/cbf.py index 99ce51136..7f6638b43 100644 --- a/aslprep/utils/cbf.py +++ b/aslprep/utils/cbf.py @@ -990,8 +990,8 @@ def fit_deltam_pcasl( # TODO: Figure out what alpha_bs and m0_b should be. popt = scipy.optimize.curve_fit( calculate_deltam_pcasl, - deltam_voxel, - X=(labeleff, labeleff, t1blood, m0_a, m0_a, tau, plds_voxel), + xdata=(labeleff, labeleff, t1blood, m0_a, m0_a, tau, plds_voxel), + ydata=deltam_voxel, # initial estimates for DVs cbf=1, att=1, diff --git a/docs/workflows.rst b/docs/workflows.rst index 3c7e54079..3ff13cd78 100644 --- a/docs/workflows.rst +++ b/docs/workflows.rst @@ -448,7 +448,7 @@ Equation 2: &= 0 \hfilll { 0 < LD + PLD < ATT } - \Delta{M}_tiss &= \frac{ 2 \cdot \alpha \cdot \alpha_{BS} \cdot T_{1b} \cdot M_{0a} \cdot CBF + \Delta{M}_{tiss} &= \frac{ 2 \cdot \alpha \cdot \alpha_{BS} \cdot T_{1b} \cdot M_{0a} \cdot CBF \cdot \left[ 1 - e ^ {-\frac{ LD + PLD - ATT } { T_{1,tissue} }} \right] }{ 6000 \cdot e ^ { \frac{ ATT } { T_{1,blood} } } } \hfilll { ATT < LD + PLD < ATT + LD } From 96095aee7cce3965d1e3c7621414f12add1b920b Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Fri, 22 Dec 2023 09:21:46 -0500 Subject: [PATCH 18/73] Try that. --- aslprep/utils/cbf.py | 23 +++++++++++++++++++++-- 1 file changed, 21 insertions(+), 2 deletions(-) diff --git a/aslprep/utils/cbf.py b/aslprep/utils/cbf.py index 7f6638b43..58714b866 100644 --- a/aslprep/utils/cbf.py +++ b/aslprep/utils/cbf.py @@ -884,7 +884,15 @@ def calculate_deltam_pcasl(X, cbf, att, abat, abv): Post-labeling delays. Used for both deltam_tiss and deltam_art calculation. """ - labeleff, alpha_bs, t1blood, m0_a, m0_b, tau, pld = X + # float parameters + labeleff = X[0, 0] + alpha_bs = X[0, 1] + t1blood = X[0, 2] + m0_a = X[0, 3] + m0_b = X[0, 4] + # array parameters + tau = X[:, 5] + pld = X[:, 6] if (tau + pld) < att: # 0 < LD + PLD < ATT @@ -987,10 +995,21 @@ def fit_deltam_pcasl( deltam_voxel = deltam_arr[i_voxel, :] plds_voxel = plds[i_voxel, :] + # The independent variables used to estimate cbf, etc. are either floats or arrays, + # but curve_fit needs them all to be the same size/shape. + xdata = np.zeros((n_volumes, 7)) + xdata[0, 0] = labeleff + xdata[0, 1] = labeleff + xdata[0, 2] = t1blood + xdata[0, 3] = m0_a + xdata[0, 4] = m0_a + xdata[:, 5] = tau + xdata[:, 6] = plds_voxel + # TODO: Figure out what alpha_bs and m0_b should be. popt = scipy.optimize.curve_fit( calculate_deltam_pcasl, - xdata=(labeleff, labeleff, t1blood, m0_a, m0_a, tau, plds_voxel), + xdata=xdata, ydata=deltam_voxel, # initial estimates for DVs cbf=1, From 9f772940ddac4eb69c11a65d8f217cadf06e8294 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Fri, 22 Dec 2023 09:45:59 -0500 Subject: [PATCH 19/73] Update cbf.py --- aslprep/utils/cbf.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/aslprep/utils/cbf.py b/aslprep/utils/cbf.py index 58714b866..e67807c9e 100644 --- a/aslprep/utils/cbf.py +++ b/aslprep/utils/cbf.py @@ -994,6 +994,7 @@ def fit_deltam_pcasl( for i_voxel in range(n_voxels): deltam_voxel = deltam_arr[i_voxel, :] plds_voxel = plds[i_voxel, :] + m0_a_voxel = m0_a[i_voxel] # The independent variables used to estimate cbf, etc. are either floats or arrays, # but curve_fit needs them all to be the same size/shape. @@ -1001,8 +1002,8 @@ def fit_deltam_pcasl( xdata[0, 0] = labeleff xdata[0, 1] = labeleff xdata[0, 2] = t1blood - xdata[0, 3] = m0_a - xdata[0, 4] = m0_a + xdata[0, 3] = m0_a_voxel + xdata[0, 4] = m0_a_voxel xdata[:, 5] = tau xdata[:, 6] = plds_voxel From 2d49598c56fe94043ef43fde52e123e583d88781 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Fri, 22 Dec 2023 10:25:07 -0500 Subject: [PATCH 20/73] More work. --- aslprep/utils/cbf.py | 113 +++++++++++++++++++++++++++++++++++++++---- docs/workflows.rst | 6 +-- 2 files changed, 107 insertions(+), 12 deletions(-) diff --git a/aslprep/utils/cbf.py b/aslprep/utils/cbf.py index e67807c9e..9a0f0e41d 100644 --- a/aslprep/utils/cbf.py +++ b/aslprep/utils/cbf.py @@ -900,15 +900,15 @@ def calculate_deltam_pcasl(X, cbf, att, abat, abv): elif att < (tau + pld) < (att + tau): # ATT < LD + PLD < ATT + LD deltam_tiss = ( - (2 * labeleff * alpha_bs * t1blood * m0_a * cbf * (np.exp ** (-att / t1blood))) - * (1 - (np.exp ** (-(tau + pld - att) / t1blood))) + (2 * labeleff * alpha_bs * t1blood * m0_a * cbf * (np.e ** (-att / t1blood))) + * (1 - (np.e ** (-(tau + pld - att) / t1blood))) / 6000 ) else: # ATT < PLD deltam_tiss = ( - (2 * labeleff * alpha_bs * t1blood * m0_a * cbf * (np.exp ** (-pld / t1blood))) - * (1 - (np.exp ** (-tau / t1blood))) + (2 * labeleff * alpha_bs * t1blood * m0_a * cbf * (np.e ** (-pld / t1blood))) + * (1 - (np.e ** (-tau / t1blood))) / 6000 ) @@ -918,7 +918,7 @@ def calculate_deltam_pcasl(X, cbf, att, abat, abv): deltam_art = 0 elif abat < (tau + pld) < (abat + tau): # aBAT < LD + PLD < aBAT + LD - deltam_art = 2 * labeleff * alpha_bs * m0_b * abv * (np.exp ** (-abat / t1blood)) + deltam_art = 2 * labeleff * alpha_bs * m0_b * abv * (np.e ** (-abat / t1blood)) else: # aBAT < PLD deltam_art = 0 @@ -1013,10 +1013,7 @@ def fit_deltam_pcasl( xdata=xdata, ydata=deltam_voxel, # initial estimates for DVs - cbf=1, - att=1, - abat=1, - abv=1, + p0=(1, 1, 1, 1), # lower and upper bounds for DVs bounds=((0, 0, 0, 0), (np.inf, np.inf, np.inf, np.inf)), )[0] @@ -1026,3 +1023,101 @@ def fit_deltam_pcasl( abv[i_voxel] = popt[3] return cbf, att, abat, abv + + +def calculate_deltam_pasl(X, cbf, att, abat, abv): + """Specify a model for use with scipy curve fitting. + + The model is fit separately for each voxel. + + This model is used to estimate ``cbf``, ``att``, ``abat``, and ``abv`` for multi-PLD PCASL + data. + + Parameters + ---------- + cbf : :obj:`float` + Cerebral blood flow. + Used for deltam_tiss calculation. + Estimated by the model. + att : :obj:`float` + Arterial transit time. + Used for deltam_tiss calculation. + Estimated by the model. + abat : :obj:`float` + The arterial bolus arrival time, describing the first arrival of the labeled blood + water in the voxel within the arterial vessel. Used for deltam_art calculation. + Estimated by the model. + abv : :obj:`float` + Arterial blood volume. + The fraction of the voxel occupied by the arterial vessel. + Used for deltam_art calculation. + Estimated by the model. + X : :obj:`tuple` of length 7 + Dependent variables. + + Returns + ------- + deltam : :obj:`float` + Delta-M value for the voxel. + + Notes + ----- + X breaks down into seven variables: + + labeleff : :obj:`float` + Labeling efficiency. + Used for both deltam_tiss and deltam_art calculation. + alpha_bs : :obj:`float` + Total background suppression inversion efficiency. + Used for both deltam_tiss and deltam_art calculation. + t1blood : :obj:`float` + Longitudinal relaxation time of arterial blood in seconds. + Used for both deltam_tiss and deltam_art calculation. + m0_a : :obj:`float` + Equilibrium magnetization of tissue, calculated as M0a = SIpd / lambda, + where SIpd is a proton density weighted image and lambda is the tissue/blood partition + coefficient. + Used for deltam_tiss calculation. + m0_b : :obj:`float` + Equilibrium magnetization of arterial blood. + Used for deltam_art calculation. + ti : :obj:`numpy.ndarray` + ti1 : ??? + """ + # float parameters + labeleff = X[0, 0] + alpha_bs = X[0, 1] + t1blood = X[0, 2] + m0_a = X[0, 3] + m0_b = X[0, 4] + # array parameters + ti = X[:, 5] + ti1 = X[:, 6] + + if ti < att: + # 0 < TI < ATT + deltam_tiss = 0 + elif att < ti < (att + ti1): + # ATT < TI < ATT + TI1 + deltam_tiss = ( + (2 * labeleff * alpha_bs * m0_a * cbf * (np.e ** (-ti / t1blood))) * (ti - att) + ) / 6000 + else: + # ATT + TI1 < TI + deltam_tiss = ( + (2 * labeleff * alpha_bs * m0_a * cbf * (np.e ** (-ti / t1blood))) * ti1 + ) / 6000 + + # Intravascular component + if ti < abat: + # 0 < TI < aBAT + deltam_art = 0 + elif abat < ti < (abat + ti1): + # aBAT < TI < aBAT + TI1 + deltam_art = 2 * labeleff * alpha_bs * m0_b * abv * (np.e ** (-ti / t1blood)) + else: + # aBAT + TI1 < TI + deltam_art = 0 + + deltam = deltam_tiss + deltam_art + return deltam diff --git a/docs/workflows.rst b/docs/workflows.rst index 3ff13cd78..838049aad 100644 --- a/docs/workflows.rst +++ b/docs/workflows.rst @@ -462,12 +462,12 @@ Equation 4: .. math:: - &= 0 \pushright{ 0 < LD + PLD < aBAT } + &= 0 \hfilll{ 0 < LD + PLD < aBAT } \Delta{M}_{art} &= 2 \cdot \alpha \cdot \alpha_{BS} \cdot M_{0b} \cdot aBV \cdot - e ^ { -\frac{ aBAT } { T_{1,blood} } } \pushright{ aBAT < LD + PLD < aBAT + LD } + e ^ { -\frac{ aBAT } { T_{1,blood} } } \hfilll{ aBAT < LD + PLD < aBAT + LD } - &= 0 \pushright{ aBAT < PLD } + &= 0 \hfilll{ aBAT < PLD } For multi-delay PCASL data, the following steps are taken: From 4cd4d89f2bc6c93569c77c6a64ab8f3a0893543b Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Fri, 22 Dec 2023 11:50:10 -0500 Subject: [PATCH 21/73] Update cbf.py --- aslprep/utils/cbf.py | 70 ++++++++++++++++++++++++-------------------- 1 file changed, 38 insertions(+), 32 deletions(-) diff --git a/aslprep/utils/cbf.py b/aslprep/utils/cbf.py index 9a0f0e41d..c27db77d5 100644 --- a/aslprep/utils/cbf.py +++ b/aslprep/utils/cbf.py @@ -853,8 +853,8 @@ def calculate_deltam_pcasl(X, cbf, att, abat, abv): Returns ------- - deltam : :obj:`float` - Delta-M value for the voxel. + deltam : :obj:`numpy.ndarray` + Delta-M values for the voxel. Notes ----- @@ -891,39 +891,45 @@ def calculate_deltam_pcasl(X, cbf, att, abat, abv): m0_a = X[0, 3] m0_b = X[0, 4] # array parameters - tau = X[:, 5] - pld = X[:, 6] + taus = X[:, 5] + plds = X[:, 6] + + deltam = np.zeros(plds.size) + + for i_pld, pld in enumerate(plds): + tau = taus[i_pld] + + if (tau + pld) < att: + # 0 < LD + PLD < ATT + deltam_tiss = 0 + elif att <= (tau + pld) < (att + tau): + # ATT < LD + PLD < ATT + LD + deltam_tiss = ( + (2 * labeleff * alpha_bs * t1blood * m0_a * cbf * (np.e ** (-att / t1blood))) + * (1 - (np.e ** (-(tau + pld - att) / t1blood))) + / 6000 + ) + else: + # 0 < ATT < PLD + deltam_tiss = ( + (2 * labeleff * alpha_bs * t1blood * m0_a * cbf * (np.e ** (-pld / t1blood))) + * (1 - (np.e ** (-tau / t1blood))) + / 6000 + ) - if (tau + pld) < att: - # 0 < LD + PLD < ATT - deltam_tiss = 0 - elif att < (tau + pld) < (att + tau): - # ATT < LD + PLD < ATT + LD - deltam_tiss = ( - (2 * labeleff * alpha_bs * t1blood * m0_a * cbf * (np.e ** (-att / t1blood))) - * (1 - (np.e ** (-(tau + pld - att) / t1blood))) - / 6000 - ) - else: - # ATT < PLD - deltam_tiss = ( - (2 * labeleff * alpha_bs * t1blood * m0_a * cbf * (np.e ** (-pld / t1blood))) - * (1 - (np.e ** (-tau / t1blood))) - / 6000 - ) + # Intravascular component + if (tau + pld) < abat: + # 0 < LD + PLD < aBAT + deltam_art = 0 + elif abat < (tau + pld) < (abat + tau): + # aBAT < LD + PLD < aBAT + LD + deltam_art = 2 * labeleff * alpha_bs * m0_b * abv * (np.e ** (-abat / t1blood)) + else: + # aBAT < PLD + deltam_art = 0 - # Intravascular component - if (tau + pld) < abat: - # 0 < LD + PLD < aBAT - deltam_art = 0 - elif abat < (tau + pld) < (abat + tau): - # aBAT < LD + PLD < aBAT + LD - deltam_art = 2 * labeleff * alpha_bs * m0_b * abv * (np.e ** (-abat / t1blood)) - else: - # aBAT < PLD - deltam_art = 0 + deltam[i_pld] = deltam_tiss + deltam_art - deltam = deltam_tiss + deltam_art return deltam From b5b2d9f699210031d31d272349c975a15128ce00 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Fri, 22 Dec 2023 13:13:44 -0500 Subject: [PATCH 22/73] Give the new outputs a try. --- aslprep/data/aslprep_bids_config.json | 6 +- aslprep/interfaces/cbf.py | 37 ++++++++- aslprep/utils/cbf.py | 112 +++++++++++++++++++++++++- aslprep/workflows/asl/base.py | 6 +- aslprep/workflows/asl/outputs.py | 6 ++ 5 files changed, 155 insertions(+), 12 deletions(-) diff --git a/aslprep/data/aslprep_bids_config.json b/aslprep/data/aslprep_bids_config.json index 40b1da702..7110fb5b3 100755 --- a/aslprep/data/aslprep_bids_config.json +++ b/aslprep/data/aslprep_bids_config.json @@ -37,9 +37,9 @@ "sub-{subject}[/ses-{session}]/[{datatype|func}/]sub-{subject}[_ses-{session}]_task-{task}[_acq-{acquisition}][_rec-{reconstruction}][_run-{run}][_echo-{echo}][_recording-{recording}]_{suffix}{extension<.tsv.gz|.json>|.tsv.gz}", "sub-{subject}[/ses-{session}]/{datatype|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|timeseries}{extension<.json|.tsv|.csv|>|.tsv}", "sub-{subject}[/ses-{session}]/{datatype|perf}/sub-{subject}[_ses-{session}][_task-{task}][_acq-{acquisition}][_ce-{ceagent}][_dir-{direction}][_rec-{reconstruction}][_run-{run}][_echo-{echo}][_desc-{desc}]_{suffix}{extension<.tsv|.json>|.tsv}", - "sub-{subject}[/ses-{session}]/{datatype|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}{extension<.nii|.nii.gz|.json|.csv|.tsv>|.nii.gz}", - "sub-{subject}[/ses-{session}]/{datatype|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}{extension<.dtseries.nii|.dscalar.nii|.json|.dscalar.json|.dtseries.json|.func.gii|.func.json>|.dtseries.nii}", - "sub-{subject}[/ses-{session}]/{datatype|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}{extension<.svg|.png|.html>|.svg}", + "sub-{subject}[/ses-{session}]/{datatype|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}{extension<.nii|.nii.gz|.json|.csv|.tsv>|.nii.gz}", + "sub-{subject}[/ses-{session}]/{datatype|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}{extension<.dtseries.nii|.dscalar.nii|.json|.dscalar.json|.dtseries.json|.func.gii|.func.json>|.dtseries.nii}", + "sub-{subject}[/ses-{session}]/{datatype|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}{extension<.svg|.png|.html>|.svg}", "sub-{subject}/{datatype}/sub-{subject}[_ses-{session}][_acq-{acquisition}][_ce-{ceagent}][_rec-{reconstruction}][_run-{run}][_space-{space}][_atlas-{atlas}][_cohort-{cohort}][_desc-{desc}]_{suffix}{extension<.html|.svg|.png>}", "sub-{subject}/{datatype}/sub-{subject}[_ses-{session}][_acq-{acquisition}][_ce-{ceagent}][_rec-{reconstruction}][_run-{run}][_space-{space}][_atlas-{atlas}][_cohort-{cohort}][_desc-{desc}]_{suffix}{extension<.html|.svg|.png>}", "sub-{subject}[/ses-{session}]/{datatype|anat}/sub-{subject}[_ses-{session}][_acq-{acquisition}][_ce-{ceagent}][_rec-{reconstruction}][_run-{run}]_from-{from}_to-{to}_mode-{mode|image}_{suffix|xfm}{extension<.txt|.h5|.json>}", diff --git a/aslprep/interfaces/cbf.py b/aslprep/interfaces/cbf.py index ffd89e00f..15ff111d2 100644 --- a/aslprep/interfaces/cbf.py +++ b/aslprep/interfaces/cbf.py @@ -24,7 +24,13 @@ estimate_labeling_efficiency, pcasl_or_pasl, ) -from aslprep.utils.cbf import _getcbfscore, _scrubcbf, estimate_t1, fit_deltam_pcasl +from aslprep.utils.cbf import ( + _getcbfscore, + _scrubcbf, + estimate_t1, + fit_deltam_pasl, + fit_deltam_pcasl, +) class _RefineMaskInputSpec(BaseInterfaceInputSpec): @@ -524,9 +530,32 @@ def _run_interface(self, runtime): ) 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] + + else: + raise ValueError( + f"Unsupported BolusCutOffTechnique ({metadata['BolusCutOffTechnique']}) " + "for multi-PLD data." + ) + + cbf, att, abat, abv = fit_deltam_pasl( + deltam_arr, + scaled_m0data, + plds, + ti1, + labeleff, + t1blood=t1blood, + partition_coefficient=PARTITION_COEF, ) mean_cbf_img = masker.inverse_transform(cbf) diff --git a/aslprep/utils/cbf.py b/aslprep/utils/cbf.py index c27db77d5..43bf3c109 100644 --- a/aslprep/utils/cbf.py +++ b/aslprep/utils/cbf.py @@ -978,10 +978,16 @@ def fit_deltam_pcasl( n_voxels, n_volumes = deltam_arr.shape if deltam_arr.shape != plds.shape: - raise ValueError(f"deltam ({deltam_arr.shape}) != plds ({plds.shape})") + raise ValueError( + f"Number of PostLabelingDelays ({plds.shape}) does not match " + f"number of delta-M volumes ({deltam_arr.shape})." + ) if scaled_m0data.size != n_voxels: - raise ValueError(f"scaled_m0data ({scaled_m0data.size}) != deltam_arr ({n_voxels})") + raise ValueError( + f"Number of voxels in M0 data ({scaled_m0data.size}) does not match " + f"number of delta-M voxels ({n_voxels})." + ) if isinstance(tau, float): tau = np.full((n_volumes,), tau) @@ -989,7 +995,10 @@ def fit_deltam_pcasl( tau = np.ndarray(tau) if tau.size != n_volumes: - raise ValueError(f"tau ({tau.shape}) != n_volumes {n_volumes}") + raise ValueError( + f"Number of values in tau ({tau.shape}) does not match number of " + f"delta-M volumes ({n_volumes})." + ) m0_a = scaled_m0data / partition_coefficient @@ -1018,7 +1027,7 @@ def fit_deltam_pcasl( calculate_deltam_pcasl, xdata=xdata, ydata=deltam_voxel, - # initial estimates for DVs + # initial estimates for DVs (cbf, att, abat, abv) p0=(1, 1, 1, 1), # lower and upper bounds for DVs bounds=((0, 0, 0, 0), (np.inf, np.inf, np.inf, np.inf)), @@ -1127,3 +1136,98 @@ def calculate_deltam_pasl(X, cbf, att, abat, abv): deltam = deltam_tiss + deltam_art return deltam + + +def fit_deltam_pasl( + deltam_arr, + scaled_m0data, + plds, + ti1, + labeleff, + t1blood, + partition_coefficient, +): + """Estimate CBF, ATT, aBAT, and aBV from multi-PLD PCASL data. + + This function uses the general kinetic model specified by + :footcite:t:`woods2023recommendations`. + + Parameters + ---------- + deltam + plds + ti1 + t1blood + labeleff + scaled_m0data + partition_coefficient + + Returns + ------- + cbf + att + abat + abv + + Notes + ----- + This model does not include a separate term for T1,tissue, but we could expand it to do so + in the future. + + References + ---------- + .. footbibliography:: + """ + import scipy + + n_voxels, n_volumes = deltam_arr.shape + if deltam_arr.shape != plds.shape: + raise ValueError( + f"Number of PostLabelingDelays ({plds.shape}) does not match " + f"number of delta-M volumes ({deltam_arr.shape})." + ) + + if scaled_m0data.size != n_voxels: + raise ValueError( + f"Number of voxels in M0 data ({scaled_m0data.size}) does not match " + f"number of delta-M voxels ({n_voxels})." + ) + + m0_a = scaled_m0data / partition_coefficient + + cbf = np.zeros(n_voxels) + att = np.zeros(n_voxels) + abat = np.zeros(n_voxels) + abv = np.zeros(n_voxels) + for i_voxel in range(n_voxels): + deltam_voxel = deltam_arr[i_voxel, :] + plds_voxel = plds[i_voxel, :] + m0_a_voxel = m0_a[i_voxel] + + # The independent variables used to estimate cbf, etc. are either floats or arrays, + # but curve_fit needs them all to be the same size/shape. + xdata = np.zeros((n_volumes, 7)) + xdata[0, 0] = labeleff + xdata[0, 1] = labeleff + xdata[0, 2] = t1blood + xdata[0, 3] = m0_a_voxel + xdata[0, 4] = m0_a_voxel + xdata[0, 5] = ti1 + xdata[:, 6] = plds_voxel + + # TODO: Figure out what alpha_bs and m0_b should be. + popt = scipy.optimize.curve_fit( + calculate_deltam_pasl, + xdata=xdata, + ydata=deltam_voxel, + # initial estimates for DVs (cbf, att, abat, abv) + p0=(1, 1, 1, 1), + # lower and upper bounds for DVs + bounds=((0, 0, 0, 0), (np.inf, np.inf, np.inf, np.inf)), + )[0] + cbf[i_voxel] = popt[0] + att[i_voxel] = popt[1] + abat[i_voxel] = popt[2] + abv[i_voxel] = popt[3] + + return cbf, att, abat, abv diff --git a/aslprep/workflows/asl/base.py b/aslprep/workflows/asl/base.py index 88e46b0df..41cccc131 100644 --- a/aslprep/workflows/asl/base.py +++ b/aslprep/workflows/asl/base.py @@ -232,7 +232,11 @@ def init_asl_wf( cbf_4d_derivs = [] if is_multi_pld: - att_derivs += ["att"] + att_derivs += [ + "att", + "abat", + "abv", + ] else: cbf_4d_derivs += ["cbf_ts"] diff --git a/aslprep/workflows/asl/outputs.py b/aslprep/workflows/asl/outputs.py index e3557d4eb..32527d8ea 100644 --- a/aslprep/workflows/asl/outputs.py +++ b/aslprep/workflows/asl/outputs.py @@ -35,6 +35,12 @@ "att": { "suffix": "att", }, + "abat": { + "suffix": "abat", + }, + "abv": { + "suffix": "abv", + }, # SCORE/SCRUB outputs "cbf_ts_score": { "desc": "scoreTimeseries", From 6b00c3243e75c00f9dfe12d30c8607877c2c9928 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Fri, 22 Dec 2023 13:23:55 -0500 Subject: [PATCH 23/73] Update plotting.py --- aslprep/workflows/asl/plotting.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/aslprep/workflows/asl/plotting.py b/aslprep/workflows/asl/plotting.py index c622edf59..cb4f4a543 100644 --- a/aslprep/workflows/asl/plotting.py +++ b/aslprep/workflows/asl/plotting.py @@ -63,6 +63,8 @@ def init_cbf_reporting_wf( "cifti_cbf_ts", # Multi-delay outputs "att", + "abat", + "abv", # SCORE/SCRUB outputs "cbf_ts_score", # unused "mean_cbf_score", From 2e07bfc471aba21bd6d1132e5a5c3464a977c6cd Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Fri, 22 Dec 2023 13:57:55 -0500 Subject: [PATCH 24/73] reorg --- aslprep/interfaces/cbf.py | 44 +++--- aslprep/utils/cbf.py | 272 ++++++++++++++------------------------ 2 files changed, 112 insertions(+), 204 deletions(-) diff --git a/aslprep/interfaces/cbf.py b/aslprep/interfaces/cbf.py index 15ff111d2..ec23ac707 100644 --- a/aslprep/interfaces/cbf.py +++ b/aslprep/interfaces/cbf.py @@ -24,13 +24,7 @@ estimate_labeling_efficiency, pcasl_or_pasl, ) -from aslprep.utils.cbf import ( - _getcbfscore, - _scrubcbf, - estimate_t1, - fit_deltam_pasl, - fit_deltam_pcasl, -) +from aslprep.utils.cbf import _getcbfscore, _scrubcbf, estimate_t1, fit_deltam_multipld class _RefineMaskInputSpec(BaseInterfaceInputSpec): @@ -514,20 +508,11 @@ 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 = metadata["LabelingDuration"] - if is_multi_pld: + ti1 = None + tau = None if is_casl: - cbf, att, abat, abv = fit_deltam_pcasl( - deltam_arr, - scaled_m0data, - plds, - tau, - labeleff, - t1blood=t1blood, - partition_coefficient=PARTITION_COEF, - ) + tau = metadata["LabelingDuration"] else: if metadata["BolusCutOffTechnique"] == "QUIPSSII": @@ -548,15 +533,17 @@ def _run_interface(self, runtime): "for multi-PLD data." ) - cbf, att, abat, abv = fit_deltam_pasl( - deltam_arr, - scaled_m0data, - plds, - ti1, - labeleff, - t1blood=t1blood, - partition_coefficient=PARTITION_COEF, - ) + 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) @@ -586,6 +573,7 @@ def _run_interface(self, runtime): else: # Single-delay if is_casl: + tau = metadata["LabelingDuration"] denom_factor = t1blood * (1 - np.exp(-(tau / t1blood))) elif not metadata["BolusCutOffFlag"]: diff --git a/aslprep/utils/cbf.py b/aslprep/utils/cbf.py index 43bf3c109..ccf042184 100644 --- a/aslprep/utils/cbf.py +++ b/aslprep/utils/cbf.py @@ -848,7 +848,7 @@ def calculate_deltam_pcasl(X, cbf, att, abat, abv): The fraction of the voxel occupied by the arterial vessel. Used for deltam_art calculation. Estimated by the model. - X : :obj:`tuple` of length 7 + X : :obj:`numpy.ndarray` of shape (n_plds, 6) Dependent variables. Returns @@ -858,13 +858,12 @@ def calculate_deltam_pcasl(X, cbf, att, abat, abv): Notes ----- - X breaks down into seven variables: + X breaks down into six variables: labeleff : :obj:`float` Labeling efficiency. - Used for both deltam_tiss and deltam_art calculation. - alpha_bs : :obj:`float` - Total background suppression inversion efficiency. + If background suppression was performed, then this is already adjusted and combines + alpha_bs as well. Used for both deltam_tiss and deltam_art calculation. t1blood : :obj:`float` Longitudinal relaxation time of arterial blood in seconds. @@ -877,22 +876,21 @@ def calculate_deltam_pcasl(X, cbf, att, abat, abv): m0_b : :obj:`float` Equilibrium magnetization of arterial blood. Used for deltam_art calculation. - tau : :obj:`numpy.ndarray` - Labeling durations. - Used for both deltam_tiss and deltam_art calculation. - pld : :obj:`numpy.ndarray` + plds : :obj:`numpy.ndarray` Post-labeling delays. Used for both deltam_tiss and deltam_art calculation. + taus : :obj:`numpy.ndarray` + Labeling durations. + Used for both deltam_tiss and deltam_art calculation. """ # float parameters labeleff = X[0, 0] - alpha_bs = X[0, 1] - t1blood = X[0, 2] - m0_a = X[0, 3] - m0_b = X[0, 4] + t1blood = X[0, 1] + m0_a = X[0, 2] + m0_b = X[0, 3] # array parameters + plds = X[:, 4] taus = X[:, 5] - plds = X[:, 6] deltam = np.zeros(plds.size) @@ -905,14 +903,14 @@ def calculate_deltam_pcasl(X, cbf, att, abat, abv): elif att <= (tau + pld) < (att + tau): # ATT < LD + PLD < ATT + LD deltam_tiss = ( - (2 * labeleff * alpha_bs * t1blood * m0_a * cbf * (np.e ** (-att / t1blood))) + (2 * labeleff * t1blood * m0_a * cbf * (np.e ** (-att / t1blood))) * (1 - (np.e ** (-(tau + pld - att) / t1blood))) / 6000 ) else: # 0 < ATT < PLD deltam_tiss = ( - (2 * labeleff * alpha_bs * t1blood * m0_a * cbf * (np.e ** (-pld / t1blood))) + (2 * labeleff * t1blood * m0_a * cbf * (np.e ** (-pld / t1blood))) * (1 - (np.e ** (-tau / t1blood))) / 6000 ) @@ -923,7 +921,7 @@ def calculate_deltam_pcasl(X, cbf, att, abat, abv): deltam_art = 0 elif abat < (tau + pld) < (abat + tau): # aBAT < LD + PLD < aBAT + LD - deltam_art = 2 * labeleff * alpha_bs * m0_b * abv * (np.e ** (-abat / t1blood)) + deltam_art = 2 * labeleff * m0_b * abv * (np.e ** (-abat / t1blood)) else: # aBAT < PLD deltam_art = 0 @@ -933,113 +931,6 @@ def calculate_deltam_pcasl(X, cbf, att, abat, abv): return deltam -def fit_deltam_pcasl( - deltam_arr, - scaled_m0data, - plds, - tau, - labeleff, - t1blood, - partition_coefficient, -): - """Estimate CBF, ATT, aBAT, and aBV from multi-PLD PCASL data. - - This function uses the general kinetic model specified by - :footcite:t:`woods2023recommendations`. - - Parameters - ---------- - deltam - plds - tau : :obj:`float` or :obj:`numpy.ndarray` - Label duration. May be a single value or may vary across volumes/PLDs. - t1blood - labeleff - scaled_m0data - partition_coefficient - - Returns - ------- - cbf - att - abat - abv - - Notes - ----- - This model does not include a separate term for T1,tissue, but we could expand it to do so - in the future. - - References - ---------- - .. footbibliography:: - """ - import scipy - - n_voxels, n_volumes = deltam_arr.shape - if deltam_arr.shape != plds.shape: - raise ValueError( - f"Number of PostLabelingDelays ({plds.shape}) does not match " - f"number of delta-M volumes ({deltam_arr.shape})." - ) - - if scaled_m0data.size != n_voxels: - raise ValueError( - f"Number of voxels in M0 data ({scaled_m0data.size}) does not match " - f"number of delta-M voxels ({n_voxels})." - ) - - if isinstance(tau, float): - tau = np.full((n_volumes,), tau) - else: - tau = np.ndarray(tau) - - if tau.size != n_volumes: - raise ValueError( - f"Number of values in tau ({tau.shape}) does not match number of " - f"delta-M volumes ({n_volumes})." - ) - - m0_a = scaled_m0data / partition_coefficient - - cbf = np.zeros(n_voxels) - att = np.zeros(n_voxels) - abat = np.zeros(n_voxels) - abv = np.zeros(n_voxels) - for i_voxel in range(n_voxels): - deltam_voxel = deltam_arr[i_voxel, :] - plds_voxel = plds[i_voxel, :] - m0_a_voxel = m0_a[i_voxel] - - # The independent variables used to estimate cbf, etc. are either floats or arrays, - # but curve_fit needs them all to be the same size/shape. - xdata = np.zeros((n_volumes, 7)) - xdata[0, 0] = labeleff - xdata[0, 1] = labeleff - xdata[0, 2] = t1blood - xdata[0, 3] = m0_a_voxel - xdata[0, 4] = m0_a_voxel - xdata[:, 5] = tau - xdata[:, 6] = plds_voxel - - # TODO: Figure out what alpha_bs and m0_b should be. - popt = scipy.optimize.curve_fit( - calculate_deltam_pcasl, - xdata=xdata, - ydata=deltam_voxel, - # initial estimates for DVs (cbf, att, abat, abv) - p0=(1, 1, 1, 1), - # lower and upper bounds for DVs - bounds=((0, 0, 0, 0), (np.inf, np.inf, np.inf, np.inf)), - )[0] - cbf[i_voxel] = popt[0] - att[i_voxel] = popt[1] - abat[i_voxel] = popt[2] - abv[i_voxel] = popt[3] - - return cbf, att, abat, abv - - def calculate_deltam_pasl(X, cbf, att, abat, abv): """Specify a model for use with scipy curve fitting. @@ -1067,7 +958,7 @@ def calculate_deltam_pasl(X, cbf, att, abat, abv): The fraction of the voxel occupied by the arterial vessel. Used for deltam_art calculation. Estimated by the model. - X : :obj:`tuple` of length 7 + X : :obj:`numpy.ndarray` of shape (n_plds/tis, 6) Dependent variables. Returns @@ -1077,13 +968,12 @@ def calculate_deltam_pasl(X, cbf, att, abat, abv): Notes ----- - X breaks down into seven variables: + X breaks down into six variables: labeleff : :obj:`float` Labeling efficiency. - Used for both deltam_tiss and deltam_art calculation. - alpha_bs : :obj:`float` - Total background suppression inversion efficiency. + In the case of background suppression, this combines alpha and alpha_bs, + so we don't need a separate term for alpha_bs. Used for both deltam_tiss and deltam_art calculation. t1blood : :obj:`float` Longitudinal relaxation time of arterial blood in seconds. @@ -1096,56 +986,59 @@ def calculate_deltam_pasl(X, cbf, att, abat, abv): m0_b : :obj:`float` Equilibrium magnetization of arterial blood. Used for deltam_art calculation. - ti : :obj:`numpy.ndarray` - ti1 : ??? + tis : :obj:`numpy.ndarray` + ti1 : :obj:`float` + Bolus cutoff delay time. """ # float parameters labeleff = X[0, 0] - alpha_bs = X[0, 1] - t1blood = X[0, 2] - m0_a = X[0, 3] - m0_b = X[0, 4] - # array parameters - ti = X[:, 5] - ti1 = X[:, 6] - - if ti < att: - # 0 < TI < ATT - deltam_tiss = 0 - elif att < ti < (att + ti1): - # ATT < TI < ATT + TI1 - deltam_tiss = ( - (2 * labeleff * alpha_bs * m0_a * cbf * (np.e ** (-ti / t1blood))) * (ti - att) - ) / 6000 - else: - # ATT + TI1 < TI - deltam_tiss = ( - (2 * labeleff * alpha_bs * m0_a * cbf * (np.e ** (-ti / t1blood))) * ti1 - ) / 6000 - - # Intravascular component - if ti < abat: - # 0 < TI < aBAT - deltam_art = 0 - elif abat < ti < (abat + ti1): - # aBAT < TI < aBAT + TI1 - deltam_art = 2 * labeleff * alpha_bs * m0_b * abv * (np.e ** (-ti / t1blood)) - else: - # aBAT + TI1 < TI - deltam_art = 0 + t1blood = X[0, 1] + m0_a = X[0, 2] + m0_b = X[0, 3] + tis = X[:, 4] # array parameter + ti1 = X[0, 5] + + deltam = np.zeros(tis.size) + + for i_ti, ti in enumerate(tis): + if ti < att: + # 0 < TI < ATT + deltam_tiss = 0 + elif att < ti < (att + ti1): + # ATT < TI < ATT + TI1 + deltam_tiss = ( + (2 * labeleff * m0_a * cbf * (np.e ** (-ti / t1blood))) * (ti - att) + ) / 6000 + else: + # ATT + TI1 < TI + deltam_tiss = ((2 * labeleff * m0_a * cbf * (np.e ** (-ti / t1blood))) * ti1) / 6000 + + # Intravascular component + if ti < abat: + # 0 < TI < aBAT + deltam_art = 0 + elif abat < ti < (abat + ti1): + # aBAT < TI < aBAT + TI1 + deltam_art = 2 * labeleff * m0_b * abv * (np.e ** (-ti / t1blood)) + else: + # aBAT + TI1 < TI + deltam_art = 0 + + deltam[i_ti] = deltam_tiss + deltam_art - deltam = deltam_tiss + deltam_art return deltam -def fit_deltam_pasl( +def fit_deltam_multipld( deltam_arr, scaled_m0data, plds, - ti1, labeleff, t1blood, partition_coefficient, + is_casl, + tau=None, + ti1=None, ): """Estimate CBF, ATT, aBAT, and aBV from multi-PLD PCASL data. @@ -1156,9 +1049,13 @@ def fit_deltam_pasl( ---------- deltam plds - ti1 + tau : :obj:`float` or :obj:`numpy.ndarray` + Label duration. May be a single value or may vary across volumes/PLDs. t1blood - labeleff + labeleff : :obj:`float` + Labeling efficiency, also referred to as alpha. + In the case of background suppression, this combines alpha and alpha_bs, + so we don't need a separate term for alpha_bs. scaled_m0data partition_coefficient @@ -1193,6 +1090,23 @@ def fit_deltam_pasl( f"number of delta-M voxels ({n_voxels})." ) + if is_casl: + assert ti1 is None + assert tau is not None + if isinstance(tau, float): + tau = np.full((n_volumes,), tau) + else: + tau = np.ndarray(tau) + + if tau.size != n_volumes: + raise ValueError( + f"Number of values in tau ({tau.shape}) does not match number of " + f"delta-M volumes ({n_volumes})." + ) + else: + assert ti1 is not None + assert tau is None + m0_a = scaled_m0data / partition_coefficient cbf = np.zeros(n_voxels) @@ -1204,20 +1118,25 @@ def fit_deltam_pasl( plds_voxel = plds[i_voxel, :] m0_a_voxel = m0_a[i_voxel] + # TODO: Figure out what alpha_bs and m0_b should be. # The independent variables used to estimate cbf, etc. are either floats or arrays, # but curve_fit needs them all to be the same size/shape. - xdata = np.zeros((n_volumes, 7)) + xdata = np.zeros((n_volumes, 6)) xdata[0, 0] = labeleff - xdata[0, 1] = labeleff - xdata[0, 2] = t1blood + xdata[0, 1] = t1blood + xdata[0, 2] = m0_a_voxel xdata[0, 3] = m0_a_voxel - xdata[0, 4] = m0_a_voxel - xdata[0, 5] = ti1 - xdata[:, 6] = plds_voxel + xdata[:, 4] = plds_voxel + + if is_casl: + xdata[:, 5] = tau + func = calculate_deltam_pcasl + else: + xdata[0, 5] = ti1 + func = calculate_deltam_pasl - # TODO: Figure out what alpha_bs and m0_b should be. popt = scipy.optimize.curve_fit( - calculate_deltam_pasl, + func, xdata=xdata, ydata=deltam_voxel, # initial estimates for DVs (cbf, att, abat, abv) @@ -1225,6 +1144,7 @@ def fit_deltam_pasl( # lower and upper bounds for DVs bounds=((0, 0, 0, 0), (np.inf, np.inf, np.inf, np.inf)), )[0] + cbf[i_voxel] = popt[0] att[i_voxel] = popt[1] abat[i_voxel] = popt[2] From 2dcb39453d02532f149ecd6e4906e08552ccca16 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Fri, 22 Dec 2023 14:07:40 -0500 Subject: [PATCH 25/73] Update test_cli.py --- aslprep/tests/test_cli.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/aslprep/tests/test_cli.py b/aslprep/tests/test_cli.py index c0d77302c..2df90e227 100644 --- a/aslprep/tests/test_cli.py +++ b/aslprep/tests/test_cli.py @@ -54,7 +54,7 @@ def test_examples_pasl_multipld(data_dir, output_dir, working_dir): f"{os.path.join(data_dir, 'anatomical/smriprep')}", ] - _run_and_fail(parameters) + _run_and_generate(TEST_NAME, PARTICIPANT_LABEL, parameters, out_dir) @pytest.mark.examples_pcasl_multipld From 0b60bed996fea20a7e5629304866e8793d63971d Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Fri, 22 Dec 2023 15:17:23 -0500 Subject: [PATCH 26/73] Update expected_outputs_examples_pcasl_multipld.txt --- .../tests/data/expected_outputs_examples_pcasl_multipld.txt | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/aslprep/tests/data/expected_outputs_examples_pcasl_multipld.txt b/aslprep/tests/data/expected_outputs_examples_pcasl_multipld.txt index 462ddee77..98512072e 100644 --- a/aslprep/tests/data/expected_outputs_examples_pcasl_multipld.txt +++ b/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 From 5d76dc2d6348cc21161a48408840ede962ec5e59 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Fri, 22 Dec 2023 15:17:48 -0500 Subject: [PATCH 27/73] Update test_cli.py --- aslprep/tests/test_cli.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/aslprep/tests/test_cli.py b/aslprep/tests/test_cli.py index 2df90e227..63413a4cd 100644 --- a/aslprep/tests/test_cli.py +++ b/aslprep/tests/test_cli.py @@ -54,7 +54,8 @@ 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_generate(TEST_NAME, PARTICIPANT_LABEL, parameters, out_dir) + _run_and_fail(parameters) @pytest.mark.examples_pcasl_multipld From 02b999d47a3a8dfb2dddbccd0d169981496f3b0a Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Fri, 22 Dec 2023 15:50:30 -0500 Subject: [PATCH 28/73] Update tests. --- aslprep/tests/test_cli.py | 3 +-- aslprep/tests/test_interfaces_cbf.py | 10 ++++------ 2 files changed, 5 insertions(+), 8 deletions(-) diff --git a/aslprep/tests/test_cli.py b/aslprep/tests/test_cli.py index 63413a4cd..2df90e227 100644 --- a/aslprep/tests/test_cli.py +++ b/aslprep/tests/test_cli.py @@ -54,8 +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) + _run_and_generate(TEST_NAME, PARTICIPANT_LABEL, parameters, out_dir) @pytest.mark.examples_pcasl_multipld diff --git a/aslprep/tests/test_interfaces_cbf.py b/aslprep/tests/test_interfaces_cbf.py index 791819cbe..062687194 100644 --- a/aslprep/tests/test_interfaces_cbf.py +++ b/aslprep/tests/test_interfaces_cbf.py @@ -222,7 +222,7 @@ def test_computecbf_pasl(datasets, tmp_path_factory): m0_file=m0_file, mask=mask_file, ) - with pytest.raises(ValueError, match="Multi-delay data are not supported"): + with pytest.raises(ValueError, match="Unsupported BolusCutOffTechnique"): results = interface.run(cwd=tmpdir) # Scenario 4: QUIPSS PASL with one PostLabelingDelay for each deltam volume (good) @@ -242,7 +242,7 @@ def test_computecbf_pasl(datasets, tmp_path_factory): m0_file=m0_file, mask=mask_file, ) - with pytest.raises(ValueError, match="Multi-delay data are not supported"): + with pytest.raises(ValueError, match="Unsupported BolusCutOffTechnique"): results = interface.run(cwd=tmpdir) # Scenario 5: QUIPSSII PASL with one PostLabelingDelay @@ -289,8 +289,7 @@ def test_computecbf_pasl(datasets, tmp_path_factory): m0_file=m0_file, mask=mask_file, ) - with pytest.raises(ValueError, match="Multi-delay data are not supported"): - results = interface.run(cwd=tmpdir) + results = interface.run(cwd=tmpdir) # Scenario 7: Q2TIPS PASL with one PostLabelingDelay # This should produce CBF time series and mean CBF, but no ATT @@ -336,8 +335,7 @@ def test_computecbf_pasl(datasets, tmp_path_factory): m0_file=m0_file, mask=mask_file, ) - with pytest.raises(ValueError, match="Multi-delay data are not supported"): - results = interface.run(cwd=tmpdir) + results = interface.run(cwd=tmpdir) def test_compare_slicetiming(datasets, tmp_path_factory): From 117aab6f65d66f627121e320f4e2ac44396c4f4e Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Fri, 22 Dec 2023 16:00:14 -0500 Subject: [PATCH 29/73] Update cbf.py --- aslprep/workflows/asl/cbf.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/aslprep/workflows/asl/cbf.py b/aslprep/workflows/asl/cbf.py index 2c54ffb17..ce14c4056 100644 --- a/aslprep/workflows/asl/cbf.py +++ b/aslprep/workflows/asl/cbf.py @@ -182,9 +182,7 @@ def init_cbf_wf( bcut = metadata.get("BolusCutOffTechnique") if is_multi_pld: - raise ValueError( - "Multi-delay data are not supported for PASL sequences at the moment." - ) + workflow.__desc__ += """HAHAHAHA!""" # Single-delay PASL data, with different bolus cut-off techniques if bcut == "QUIPSS": From 7dfe3d96da75b5c85727ac35a1fda33679dbf801 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Tue, 2 Jan 2024 14:25:07 -0500 Subject: [PATCH 30/73] Plot the other metrics. --- aslprep/interfaces/plotting.py | 26 +++++++++++++---- aslprep/utils/cbf.py | 6 ++-- aslprep/workflows/asl/base.py | 1 + aslprep/workflows/asl/plotting.py | 46 ++++++++++++++++++++++++------- 4 files changed, 62 insertions(+), 17 deletions(-) diff --git a/aslprep/interfaces/plotting.py b/aslprep/interfaces/plotting.py index cd3cda097..c48ed9158 100644 --- a/aslprep/interfaces/plotting.py +++ b/aslprep/interfaces/plotting.py @@ -165,8 +165,16 @@ def _run_interface(self, runtime): class _CBFByTissueTypePlotInputSpec(BaseInterfaceInputSpec): - cbf = File(exists=True, mandatory=True, desc="") + in_file = File(exists=True, mandatory=True, desc="CBF, ATT, aBAT, or ABV file.") seg_file = File(exists=True, mandatory=True, desc="Segmentation file") + img_type = traits.Enum( + "cbf", + "att", + "abat", + "abv", + mandatory=True, + desc="Image type", + ) class _CBFByTissueTypePlotOutputSpec(TraitedSpec): @@ -185,12 +193,20 @@ def _run_interface(self, runtime): from nilearn import image, masking self._results["out_file"] = fname_presuffix( - self.inputs.cbf, - suffix="_cbfplot.svg", + self.inputs.in_file, + suffix=f"_{self.inputs.img_type}plot.svg", use_ext=False, newpath=runtime.cwd, ) + column_names = { + "cbf": "CBF\n(mL/100 g/min)", + "att": "ATT\n(s)", + "abat": "aBAT\n(s)", + "abv": "ABV\n(fraction)", + } + unit_str = column_names[self.inputs.img_type] + dfs = [] for i_tissue_type, tissue_type in enumerate(["GM", "WM", "CSF"]): tissue_type_val = i_tissue_type + 1 @@ -200,7 +216,7 @@ def _run_interface(self, runtime): ) tissue_type_vals = masking.apply_mask(self.inputs.cbf, mask_img) df = pd.DataFrame( - columns=["CBF\n(mL/100 g/min)", "Tissue Type"], + columns=[unit_str, "Tissue Type"], data=list( map(list, zip(*[tissue_type_vals, [tissue_type] * tissue_type_vals.size])) ), @@ -214,7 +230,7 @@ def _run_interface(self, runtime): fig, ax = plt.subplots(figsize=(16, 8)) sns.despine(ax=ax, bottom=True, left=True) sns.boxenplot( - y="CBF\n(mL/100 g/min)", + y=unit_str, data=df, width=0.6, showfliers=True, diff --git a/aslprep/utils/cbf.py b/aslprep/utils/cbf.py index ccf042184..dc75b24c9 100644 --- a/aslprep/utils/cbf.py +++ b/aslprep/utils/cbf.py @@ -1140,9 +1140,11 @@ def fit_deltam_multipld( xdata=xdata, ydata=deltam_voxel, # initial estimates for DVs (cbf, att, abat, abv) - p0=(1, 1, 1, 1), + # Values provided by Manuel Taso + p0=(60, 1.2, 1, 0.02), # lower and upper bounds for DVs - bounds=((0, 0, 0, 0), (np.inf, np.inf, np.inf, np.inf)), + # Upper bounds provided by Manuel Taso + bounds=((0, 0, 0, 0), (300, 5, 5, 0.1)), )[0] cbf[i_voxel] = popt[0] diff --git a/aslprep/workflows/asl/base.py b/aslprep/workflows/asl/base.py index 41cccc131..dcb260693 100644 --- a/aslprep/workflows/asl/base.py +++ b/aslprep/workflows/asl/base.py @@ -492,6 +492,7 @@ def init_asl_wf( plot_timeseries=not (is_multi_pld or use_ge or (config.workflow.level == "resampling")), scorescrub=scorescrub, basil=basil, + is_multi_pld=is_multi_pld, name="cbf_reporting_wf", ) workflow.connect([ diff --git a/aslprep/workflows/asl/plotting.py b/aslprep/workflows/asl/plotting.py index cb4f4a543..bfdcbf2fc 100644 --- a/aslprep/workflows/asl/plotting.py +++ b/aslprep/workflows/asl/plotting.py @@ -20,6 +20,7 @@ def init_cbf_reporting_wf( plot_timeseries=True, scorescrub=False, basil=False, + is_multi_pld=False, name="cbf_reporting_wf", ): """Generate CBF reports. @@ -36,6 +37,7 @@ def init_cbf_reporting_wf( "RepetitionTime": 4, "RepetitionTimePreparation": 4, }, + is_multi_pld=True, ) """ from niworkflows.interfaces.images import SignalExtraction @@ -195,11 +197,11 @@ def init_cbf_reporting_wf( workflow.connect([(cbf_summary, ds_report_cbf, [("out_file", "in_file")])]) cbf_by_tt_plot = pe.Node( - CBFByTissueTypePlot(), + CBFByTissueTypePlot(img_type="cbf"), name="cbf_by_tt_plot", ) workflow.connect([ - (inputnode, cbf_by_tt_plot, [("mean_cbf", "cbf")]), + (inputnode, cbf_by_tt_plot, [("mean_cbf", "in_file")]), (warp_t1w_dseg_to_aslref, cbf_by_tt_plot, [("output_image", "seg_file")]), ]) # fmt:skip @@ -216,6 +218,30 @@ def init_cbf_reporting_wf( ) workflow.connect([(cbf_by_tt_plot, ds_report_cbf_by_tt, [("out_file", "in_file")])]) + if is_multi_pld: + for img_type in ["att", "abat", "abv"]: + img_by_tt_plot = pe.Node( + CBFByTissueTypePlot(img_type=img_type), + name=f"{img_type}_by_tt_plot", + ) + workflow.connect([ + (inputnode, img_by_tt_plot, [(img_type, "in_file")]), + (warp_t1w_dseg_to_aslref, img_by_tt_plot, [("output_image", "seg_file")]), + ]) # fmt:skip + + ds_report_img_by_tt = pe.Node( + DerivativesDataSink( + datatype="figures", + desc=f"{img_type}ByTissueType", + suffix=img_type, + keep_dtype=True, + ), + name=f"ds_report_{img_type}_by_tt", + run_without_submitting=True, + mem_gb=config.DEFAULT_MEMORY_MIN_GB, + ) + workflow.connect([(cbf_by_tt_plot, ds_report_img_by_tt, [("out_file", "in_file")])]) + if scorescrub: score_summary = pe.Node( CBFSummaryPlot(label="score", vmax=100), @@ -238,11 +264,11 @@ def init_cbf_reporting_wf( workflow.connect([(score_summary, ds_report_score, [("out_file", "in_file")])]) score_by_tt_plot = pe.Node( - CBFByTissueTypePlot(), + CBFByTissueTypePlot(img_type="cbf"), name="score_by_tt_plot", ) workflow.connect([ - (inputnode, score_by_tt_plot, [("mean_cbf_score", "cbf")]), + (inputnode, score_by_tt_plot, [("mean_cbf_score", "in_file")]), (warp_t1w_dseg_to_aslref, score_by_tt_plot, [("output_image", "seg_file")]), ]) # fmt:skip @@ -280,11 +306,11 @@ def init_cbf_reporting_wf( workflow.connect([(scrub_summary, ds_report_scrub, [("out_file", "in_file")])]) scrub_by_tt_plot = pe.Node( - CBFByTissueTypePlot(), + CBFByTissueTypePlot(img_type="cbf"), name="scrub_by_tt_plot", ) workflow.connect([ - (inputnode, scrub_by_tt_plot, [("mean_cbf_scrub", "cbf")]), + (inputnode, scrub_by_tt_plot, [("mean_cbf_scrub", "in_file")]), (warp_t1w_dseg_to_aslref, scrub_by_tt_plot, [("output_image", "seg_file")]), ]) # fmt:skip @@ -323,11 +349,11 @@ def init_cbf_reporting_wf( workflow.connect([(basil_summary, ds_report_basil, [("out_file", "in_file")])]) basil_by_tt_plot = pe.Node( - CBFByTissueTypePlot(), + CBFByTissueTypePlot(img_type="cbf"), name="basil_by_tt_plot", ) workflow.connect([ - (inputnode, basil_by_tt_plot, [("mean_cbf_basil", "cbf")]), + (inputnode, basil_by_tt_plot, [("mean_cbf_basil", "in_file")]), (warp_t1w_dseg_to_aslref, basil_by_tt_plot, [("output_image", "seg_file")]), ]) # fmt:skip @@ -365,11 +391,11 @@ def init_cbf_reporting_wf( workflow.connect([(pvc_summary, ds_report_pvc, [("out_file", "in_file")])]) pvc_by_tt_plot = pe.Node( - CBFByTissueTypePlot(), + CBFByTissueTypePlot(img_type="cbf"), name="pvc_by_tt_plot", ) workflow.connect([ - (inputnode, pvc_by_tt_plot, [("mean_cbf_gm_basil", "cbf")]), + (inputnode, pvc_by_tt_plot, [("mean_cbf_gm_basil", "in_file")]), (warp_t1w_dseg_to_aslref, pvc_by_tt_plot, [("output_image", "seg_file")]), ]) # fmt:skip From 4e76b3f5eb78b4d8fd486f6befb0aada8b2af357 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Tue, 2 Jan 2024 14:46:09 -0500 Subject: [PATCH 31/73] Update plotting.py --- aslprep/interfaces/plotting.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/aslprep/interfaces/plotting.py b/aslprep/interfaces/plotting.py index c48ed9158..918ce7d67 100644 --- a/aslprep/interfaces/plotting.py +++ b/aslprep/interfaces/plotting.py @@ -214,7 +214,7 @@ def _run_interface(self, runtime): f"(img == {tissue_type_val}).astype(int)", img=self.inputs.seg_file, ) - tissue_type_vals = masking.apply_mask(self.inputs.cbf, mask_img) + tissue_type_vals = masking.apply_mask(self.inputs.in_file, mask_img) df = pd.DataFrame( columns=[unit_str, "Tissue Type"], data=list( From 1cdf6ba19aed312a8aa140afe0a9c28e7fc3a2dc Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Tue, 2 Jan 2024 16:03:44 -0500 Subject: [PATCH 32/73] Update figures. --- aslprep/data/reports-spec.yml | 18 ++++++++++++++++++ aslprep/interfaces/plotting.py | 8 ++++---- 2 files changed, 22 insertions(+), 4 deletions(-) diff --git a/aslprep/data/reports-spec.yml b/aslprep/data/reports-spec.yml index a2491e873..29a75dc1f 100644 --- a/aslprep/data/reports-spec.yml +++ b/aslprep/data/reports-spec.yml @@ -201,6 +201,24 @@ sections: separated by tissue type. The unit is mL/100 g/min. subtitle: Mean CBF by Tissue Type + - bids: {datatype: figures, desc: attByTissueType, suffix: att} + caption: | + The maps plot the distribution of arterial transit time values for the basic output, + separated by tissue type. + The unit is seconds. + subtitle: ATT by Tissue Type + - bids: {datatype: figures, desc: abatByTissueType, suffix: abat} + caption: | + The maps plot the distribution of arterial bolus arrival time (aBAT) values for the basic output, + separated by tissue type. + The unit is seconds. + subtitle: aBAT by Tissue Type + - bids: {datatype: figures, desc: abvByTissueType, suffix: abv} + caption: | + The maps plot the distribution of arterial blood volume (aBV) values for the basic output, + separated by tissue type. + The unit is fraction. + subtitle: aBV by Tissue Type - bids: {datatype: figures, desc: score, suffix: cbf} caption: | The maps plot cerebral blood flow (CBF) for SCORE-corrected CBF. diff --git a/aslprep/interfaces/plotting.py b/aslprep/interfaces/plotting.py index 918ce7d67..5b0f9ddb5 100644 --- a/aslprep/interfaces/plotting.py +++ b/aslprep/interfaces/plotting.py @@ -200,10 +200,10 @@ def _run_interface(self, runtime): ) column_names = { - "cbf": "CBF\n(mL/100 g/min)", - "att": "ATT\n(s)", - "abat": "aBAT\n(s)", - "abv": "ABV\n(fraction)", + "cbf": "Cerebral Blood Flow\n(mL/100 g/min)", + "att": "Arterial Transit Time\n(s)", + "abat": "Arterial Bolus Arrival Time\n(s)", + "abv": "Arterial Blood Volume\n(fraction)", } unit_str = column_names[self.inputs.img_type] From ada7e476f0f9de47d2cf2ddd7f29c4028bab75dc Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Tue, 2 Jan 2024 16:41:19 -0500 Subject: [PATCH 33/73] Fix connection. --- aslprep/workflows/asl/plotting.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/aslprep/workflows/asl/plotting.py b/aslprep/workflows/asl/plotting.py index bfdcbf2fc..9c272553b 100644 --- a/aslprep/workflows/asl/plotting.py +++ b/aslprep/workflows/asl/plotting.py @@ -240,7 +240,7 @@ def init_cbf_reporting_wf( run_without_submitting=True, mem_gb=config.DEFAULT_MEMORY_MIN_GB, ) - workflow.connect([(cbf_by_tt_plot, ds_report_img_by_tt, [("out_file", "in_file")])]) + workflow.connect([(img_by_tt_plot, ds_report_img_by_tt, [("out_file", "in_file")])]) if scorescrub: score_summary = pe.Node( From 93c2a75b68626dcb5e3a7e546b3d8ee92b72d3e8 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Wed, 3 Jan 2024 09:29:36 -0500 Subject: [PATCH 34/73] Drop upper limits for check. --- aslprep/utils/cbf.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/aslprep/utils/cbf.py b/aslprep/utils/cbf.py index dc75b24c9..51e3c927c 100644 --- a/aslprep/utils/cbf.py +++ b/aslprep/utils/cbf.py @@ -1144,7 +1144,11 @@ def fit_deltam_multipld( p0=(60, 1.2, 1, 0.02), # lower and upper bounds for DVs # Upper bounds provided by Manuel Taso - bounds=((0, 0, 0, 0), (300, 5, 5, 0.1)), + bounds=( + (0, 0, 0, 0), + # (300, 5, 5, 0.1), + (np.inf, np.inf, np.inf, np.inf), + ), )[0] cbf[i_voxel] = popt[0] From db989976a7d032c6d988d519518709d7cfcce125 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Wed, 3 Jan 2024 11:13:42 -0500 Subject: [PATCH 35/73] Update formulae and add bounds back in. --- aslprep/utils/cbf.py | 9 +-- aslprep/workflows/asl/cbf.py | 33 +++++------ docs/workflows.rst | 105 +++++++++-------------------------- 3 files changed, 43 insertions(+), 104 deletions(-) diff --git a/aslprep/utils/cbf.py b/aslprep/utils/cbf.py index 51e3c927c..741116a5e 100644 --- a/aslprep/utils/cbf.py +++ b/aslprep/utils/cbf.py @@ -1139,15 +1139,16 @@ def fit_deltam_multipld( func, xdata=xdata, ydata=deltam_voxel, - # initial estimates for DVs (cbf, att, abat, abv) + # Initial estimates for DVs (CBF, ATT, aBAT, aBV) # Values provided by Manuel Taso p0=(60, 1.2, 1, 0.02), - # lower and upper bounds for DVs + # Lower and upper bounds for DVs # Upper bounds provided by Manuel Taso bounds=( (0, 0, 0, 0), - # (300, 5, 5, 0.1), - (np.inf, np.inf, np.inf, np.inf), + # Manuel Taso recommended 5 and 5 for ATT and aBAT, + # but our test data maxed out around 12 when left unbounded. + (300, 15, 15, 0.1), ), )[0] diff --git a/aslprep/workflows/asl/cbf.py b/aslprep/workflows/asl/cbf.py index ce14c4056..993e7cf4c 100644 --- a/aslprep/workflows/asl/cbf.py +++ b/aslprep/workflows/asl/cbf.py @@ -150,28 +150,24 @@ def init_cbf_wf( *ASLPrep* loaded pre-calculated cerebral blood flow (CBF) data from the ASL file. """ - elif is_casl: - if is_multi_pld: - workflow.__desc__ += f"""\ + elif is_multi_pld: + workflow.__desc__ += f"""\ *ASLPrep* calculated cerebral blood flow (CBF) from the multi-delay -{metadata['ArterialSpinLabelingType']} data using the following method. +{metadata['ArterialSpinLabelingType']} data using the general kinetic model (GKM) +:footcite:p:`buxton1998general`, as recommended and extended in +:footcite:t:`woods2023recommendations`. -First, delta-M values were averaged over time for each post-labeling delay (PLD). {m0_str} +The voxel-wise M0 values were used as both M0a and M0b in the GKM. -Next, arterial transit time (ATT) was estimated on a voxel-wise basis according to -@dai2012reduced. - -CBF was then calculated for each delay using the mean delta-M values and the estimated ATT, -according to the formula from @fan2017long. - -CBF was then averaged over delays according to @juttukonda2021characterizing, -in which an unweighted average is calculated for each voxel across all delays in which -PLD + labeling duration > ATT. +CBF, arterial transit time (ATT), arterial bolus arrival time (aBAT), +and arterial blood volume (aBV) were estimated using a nonlinear model fit with SciPy's +``curve_fit``. """ - else: - # Single-delay (P)CASL data - workflow.__desc__ += f"""\ + + elif is_casl: + # Single-delay (P)CASL data + workflow.__desc__ += f"""\ *ASLPrep* calculated cerebral blood flow (CBF) from the single-delay {metadata['ArterialSpinLabelingType']} using a single-compartment general kinetic model [@buxton1998general]. @@ -181,9 +177,6 @@ def init_cbf_wf( else: bcut = metadata.get("BolusCutOffTechnique") - if is_multi_pld: - workflow.__desc__ += """HAHAHAHA!""" - # Single-delay PASL data, with different bolus cut-off techniques if bcut == "QUIPSS": workflow.__desc__ += f"""\ diff --git a/docs/workflows.rst b/docs/workflows.rst index 838049aad..f89224a54 100644 --- a/docs/workflows.rst +++ b/docs/workflows.rst @@ -438,7 +438,7 @@ Also, multi-delay ASL allows for the estimation of arterial transit time (ATT). Pseudo-Continuous ASL --------------------- -For multi-delay PCASL data, ASLPrep uses :func:`scipy.optimize.curvefit` to estimate +For multi-delay PCASL data, *ASLPrep* uses :func:`scipy.optimize.curvefit` to estimate CBF, ATT, aBAT, and aBV values for each voxel using Equations 2 and 4 in :footcite:t:`woods2023recommendations`. @@ -450,13 +450,13 @@ Equation 2: \Delta{M}_{tiss} &= \frac{ 2 \cdot \alpha \cdot \alpha_{BS} \cdot T_{1b} \cdot M_{0a} \cdot CBF \cdot \left[ - 1 - e ^ {-\frac{ LD + PLD - ATT } { T_{1,tissue} }} - \right] }{ 6000 \cdot e ^ { \frac{ ATT } { T_{1,blood} } } } \hfilll { ATT < LD + PLD < ATT + LD } + 1 - e ^ {-\frac{ LD + PLD - ATT } { T_{1,b} }} + \right] }{ 6000 \cdot e ^ { \frac{ ATT } { T_{1,b} } } } \hfilll { ATT < LD + PLD < ATT + LD } &= \frac{ 2 \cdot \alpha \cdot \alpha_{BS} \cdot T_{1b} \cdot M_{0a} \cdot CBF \cdot \left[ - 1 - e ^ {-\frac{ LD } { T_{1,tissue} }} - \right] }{ 6000 e ^ { \frac{ PLD } { T_{1,blood} } } } \hfilll { ATT < PLD } + 1 - e ^ {-\frac{ LD } { T_{1,b} }} + \right] }{ 6000 \cdot e ^ { \frac{ PLD } { T_{1,b} } } } \hfilll { ATT < PLD } Equation 4: @@ -465,96 +465,41 @@ Equation 4: &= 0 \hfilll{ 0 < LD + PLD < aBAT } \Delta{M}_{art} &= 2 \cdot \alpha \cdot \alpha_{BS} \cdot M_{0b} \cdot aBV \cdot - e ^ { -\frac{ aBAT } { T_{1,blood} } } \hfilll{ aBAT < LD + PLD < aBAT + LD } + e ^ { -\frac{ aBAT } { T_{1,b} } } \hfilll{ aBAT < LD + PLD < aBAT + LD } &= 0 \hfilll{ aBAT < PLD } -For multi-delay PCASL data, the following steps are taken: -1. :math:`\Delta{M}` values are first averaged over time for each unique post-labeling delay value. - We shall call these :math:`\Delta{M}` in the following equations for the sake of readability. - -2. Arterial transit time is estimated on a voxel-wise basis according to - :footcite:t:`dai2012reduced`. - - 1. Define a set of possible transit times to evaluate. - The range is defined as the minimum PLD to the maximum PLD, at increments of 0.001. - - 2. Calculate the expected weighted delay (:math:`WD_{E}`) for each possible transit time - (:math:`\delta`), across PLDs (:math:`w`). - - .. math:: - - WD_{E}(\delta_{t}, w_{i}) = e ^ \frac{ -\delta_{t} } { T_{1,blood} } \cdot - \left[ - e ^ {-\frac{ max( 0, w_{i} - \delta_{t} ) } { T_{1,tissue} }} - - e ^ {-\frac{ max( 0, \tau + w_{i} - \delta_{t} ) } { T_{1,tissue} }} - \right] - - WD_{E}(\delta_{t}) = \frac{ \sum_{i=1}^{|w|} w_{i} \cdot - WD_{E}(\delta_{t},w_{i}) } { \sum_{i=1}^{|w|} WD_{E}(\delta_{t},w_{i}) } - - 3. Calculate the observed weighted delay (:math:`WD_{O}`) for the actual data, at each voxel :math:`v`. - - .. math:: - - WD_{O}(v) = \frac{ - \sum_{i=1}^{|w|} w_{i} \cdot \Delta{M}( w_{i},v ) - } - { - \sum_{i=1}^{|w|} \Delta{M}( w_{i},v ) - } - - 4. Truncate the observed weighted delays to valid delay values, - determined based on the expected weighted delays. - - .. math:: +Pulsed ASL +---------- - WD_{O}(v) = max[min(WD_{O}(v), max[WD_{E}]), min(WD_{E})] +For multi-delay Q2TIPS or QUIPSSII PCASL data, +*ASLPrep* uses :func:`scipy.optimize.curvefit` to estimate +CBF, ATT, aBAT, and aBV values for each voxel using Equations 3 and 5 in +:footcite:t:`woods2023recommendations`. - 5. Interpolate the expected weighted delay values to infer the appropriate transit time for each voxel, - based on its observed weighted delay. +Equation 3: -3. CBF is then calculated for each unique PLD value (:math:`w_{i}`) using the 2-compartment model - described in :footcite:t:`fan2017long`. +.. math:: - .. math:: + &= 0 \hfilll { 0 < TI < ATT } - CBF_{i} = 6000 \cdot \lambda \cdot \frac{ \Delta{M}_{i} }{ M_{0} } \cdot - \frac{ - e ^ \frac{ \delta }{ T_{1,blood} } - } - { - 2 \cdot \alpha \cdot T_{1,blood} \cdot - \left[ - e ^ { -\frac{ max(w_{i} - \delta, 0) }{ T_{1,tissue} } } - - - e ^ { -\frac{ max(\tau + w_{i} - \delta, 0) }{ T_{1,tissue} } } - \right] - } + \Delta{M}_{tiss} &= \frac{ 2 \cdot \alpha \cdot \alpha_{BS} \cdot M_{0a} \cdot CBF \cdot (TI - ATT) } + }{ 6000 \cdot e ^ {-\frac{ TI } { T_{1,b} } } \hfilll { ATT < TI < ATT + TI_{1} } - .. note:: + &= \frac{ 2 \cdot \alpha \cdot \alpha_{BS} \cdot M_{0a} \cdot CBF \cdot TI_{1} } + { 6000 \cdot e ^ { \frac{ TI } { T_{1,b} } } } \hfilll { ATT + TI_{1} < TI } - Note that Equation 2 in :footcite:t:`fan2017long` uses different notation. - :math:`T_{1,blood}` is referred to as :math:`T_{1a}`, - :math:`T_{1,tissue}` is referred to as :math:`T_{1t}`, - :math:`\Delta{M}` is referred to as :math:`M`, - :math:`w` is referred to as :math:`PLD`, - :math:`\delta` is referred to as :math:`ATT`, - :math:`\tau` is referred to as :math:`LD`, - and :math:`\alpha` is referred to as :math:`\epsilon`. +Equation 5: -4. CBF is then averaged over PLDs according to :footcite:t:`juttukonda2021characterizing`, - in which an unweighted average is calculated for each voxel across all PLDs (:math:`w`) in which - :math:`w + \tau \gt \delta`. +.. math:: + &= 0 \hfilll{ 0 < TI < aBAT } -Pulsed ASL ----------- + \Delta{M}_{art} &= 2 \cdot \alpha \cdot \alpha_{BS} \cdot M_{0b} \cdot aBV \cdot + e ^ { -\frac{ TI } { T_{1,b} } } \hfilll{ aBAT < TI < aBAT + TI_{1} } -.. warning:: - As of 0.3.0, ASLPrep has disabled multi-delay support for PASL data. - We plan to properly support multi-delay PASL data in the near future. + &= 0 \hfilll{ aBAT + TI_{1} < TI } Additional Denoising Options From d29b4eb951d632e220c30fc8ac3ff402c4810539 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Wed, 3 Jan 2024 11:21:39 -0500 Subject: [PATCH 36/73] Remove unused functions and citations. --- aslprep/data/boilerplate.bib | 46 ------ aslprep/interfaces/cbf.py | 11 +- aslprep/tests/test_utils_cbf.py | 23 --- aslprep/utils/cbf.py | 279 -------------------------------- 4 files changed, 4 insertions(+), 355 deletions(-) diff --git a/aslprep/data/boilerplate.bib b/aslprep/data/boilerplate.bib index 7d95da53a..a162a1d00 100644 --- a/aslprep/data/boilerplate.bib +++ b/aslprep/data/boilerplate.bib @@ -557,27 +557,6 @@ @article{groves2009combined publisher={Elsevier} } -@article{dai2012reduced, - title={Reduced resolution transit delay prescan for quantitative continuous arterial spin labeling perfusion imaging}, - author={Dai, Weiying and Robson, Philip M and Shankaranarayanan, Ajit and Alsop, David C}, - journal={Magnetic resonance in medicine}, - volume={67}, - number={5}, - pages={1252--1265}, - year={2012}, - publisher={Wiley Online Library} -} - -@article{wang2013multi, - title={Multi-delay multi-parametric arterial spin-labeled perfusion MRI in acute ischemic stroke—comparison with dynamic susceptibility contrast enhanced perfusion imaging}, - author={Wang, Danny JJ and Alger, Jeffry R and Qiao, Joe X and Gunther, Matthias and Pope, Whitney B and Saver, Jeffrey L and Salamon, Noriko and Liebeskind, David S and UCLA Stroke Investigators and others}, - journal={NeuroImage: Clinical}, - volume={3}, - pages={1--7}, - year={2013}, - publisher={Elsevier} -} - @article{clement2022asl, title={ASL-BIDS, the brain imaging data structure extension for arterial spin labeling}, author={Clement, Patricia and Castellaro, Marco and Okell, Thomas W and Thomas, David L and Vandemaele, Pieter and Elgayar, Sara and Oliver-Taylor, Aaron and Kirk, Thomas and Woods, Joseph G and Vos, Sjoerd B and others}, @@ -615,31 +594,6 @@ @article{wang2005amplitude doi={10.1148/radiol.2351031663} } -@article{fan2017long, - title={Long-delay arterial spin labeling provides more accurate cerebral blood flow measurements in moyamoya patients: a simultaneous positron emission tomography/MRI study}, - author={Fan, Audrey P and Guo, Jia and Khalighi, Mohammad M and Gulaka, Praveen K and Shen, Bin and Park, Jun Hyung and Gandhi, Harsh and Holley, Dawn and Rutledge, Omar and Singh, Prachi and others}, - journal={Stroke}, - volume={48}, - number={9}, - pages={2441--2449}, - year={2017}, - publisher={Am Heart Assoc}, - url={https://doi.org/10.1161/STROKEAHA.117.017773}, - doi={10.1161/STROKEAHA.117.017773} -} - -@article{juttukonda2021characterizing, - title={Characterizing cerebral hemodynamics across the adult lifespan with arterial spin labeling MRI data from the Human Connectome Project-Aging}, - author={Juttukonda, Meher R and Li, Binyin and Almaktoum, Randa and Stephens, Kimberly A and Yochim, Kathryn M and Yacoub, Essa and Buckner, Randy L and Salat, David H}, - journal={Neuroimage}, - volume={230}, - pages={117807}, - year={2021}, - publisher={Elsevier}, - url={https://doi.org/10.1016/j.neuroimage.2021.117807}, - doi={10.1016/j.neuroimage.2021.117807} -} - @article{noguchi2015technical, title={A technical perspective for understanding quantitative arterial spin-labeling MR imaging using Q2TIPS}, author={Noguchi, Tomoyuki and Nishihara, Masashi and Hara, Yukiko and Hirai, Tetsuyoshi and Egashira, Yoshiaki and Azama, Shinya and Irie, Hiroyuki}, diff --git a/aslprep/interfaces/cbf.py b/aslprep/interfaces/cbf.py index ec23ac707..414dedff2 100644 --- a/aslprep/interfaces/cbf.py +++ b/aslprep/interfaces/cbf.py @@ -365,13 +365,8 @@ class ComputeCBF(SimpleInterface): Single-delay CBF, for both (P)CASL and QUIPSSII PASL is calculated according to :footcite:t:`alsop_recommended_2015`. - Multi-delay CBF is handled using a weighted average, - based on :footcite:t:`dai2012reduced,wang2013multi`. - Multi-delay CBF is calculated according to :footcite:t:`fan2017long`, - although CBF is averaged across PLDs according to the method in - :footcite:t:`juttukonda2021characterizing`. - Arterial transit time is estimated according to :footcite:t:`dai2012reduced`. + Multi-delay CBF is calculated according to :footcite:t:`woods2023recommendations`. If slice timing information is detected, then PLDs will be shifted by the slice times. @@ -381,7 +376,9 @@ class ComputeCBF(SimpleInterface): :func:`~aslprep.utils.asl.determine_multi_pld` :func:`~aslprep.utils.cbf.estimate_t1` :func:`~aslprep.utils.asl.estimate_labeling_efficiency` - :func:`~aslprep.utils.cbf.estimate_cbf_pcasl_multipld` + :func:`~aslprep.utils.cbf.calculate_deltam_pasl` + :func:`~aslprep.utils.cbf.calculate_deltam_pcasl` + :func:`~aslprep.utils.cbf.fit_deltam_multipld` References ---------- diff --git a/aslprep/tests/test_utils_cbf.py b/aslprep/tests/test_utils_cbf.py index 8010d5442..132c5361b 100644 --- a/aslprep/tests/test_utils_cbf.py +++ b/aslprep/tests/test_utils_cbf.py @@ -38,26 +38,3 @@ def test_estimate_att_pcasl_2d(): t1tissue=1.3, ) assert att.shape == (n_voxels,) - - -def test_estimate_cbf_pcasl_multipld(): - """Smoke test aslprep.utils.cbf.estimate_cbf_pcasl_multipld with slice-shifted PLDs.""" - n_voxels, n_unique_plds, n_slice_times = 1000, 6, 50 - base_plds = np.linspace(0.25, 1.5, n_unique_plds) - slice_times = np.linspace(0, 3, n_slice_times) - slice_shifted_plds = slice_times[:, None] + base_plds[None, :] - plds = np.repeat(slice_shifted_plds, n_voxels // n_slice_times, axis=0) - deltam_arr = np.random.random((n_voxels, n_unique_plds)) - scaled_m0data = np.random.random(n_voxels) - - cbf.estimate_cbf_pcasl_multipld( - deltam_arr=deltam_arr, - scaled_m0data=scaled_m0data, - plds=plds, - tau=np.array(1.4), - labeleff=0.7, - t1blood=1.6, - t1tissue=1.3, - unit_conversion=6000, - partition_coefficient=0.9, - ) diff --git a/aslprep/utils/cbf.py b/aslprep/utils/cbf.py index 741116a5e..60364964f 100644 --- a/aslprep/utils/cbf.py +++ b/aslprep/utils/cbf.py @@ -2,7 +2,6 @@ from __future__ import annotations import numpy as np -from scipy.interpolate import interp1d from scipy.stats import median_abs_deviation from aslprep import config @@ -485,284 +484,6 @@ def _scrubcbf(cbf_ts, gm, wm, csf, mask, wfun="huber", thresh=0.7): return newcbf -def estimate_att_pcasl(deltam_arr, plds, lds, t1blood, t1tissue): - """Estimate arterial transit time using the weighted average method. - - The weighted average method comes from :footcite:t:`dai2012reduced`. - - Parameters - ---------- - deltam_arr : :obj:`numpy.ndarray` of shape (S, D) - Delta-M array, averaged by PLD. - plds : :obj:`numpy.ndarray` of shape (S, D) - Post-labeling delays. w in Dai 2012. - In case of a 2D acquisition, PLDs may vary by slice, and thus the plds array will vary - in the spatial dimension. - For 3D acquisitions, or 2D acquisitions without slice timing info, plds will only vary - along the second dimension. - lds : :obj:`numpy.ndarray` - Labeling durations. tau in Dai 2012. - t1blood : :obj:`float` - T1 relaxation rate for blood. - t1tissue : :obj:`float` - T1 relaxation rate for tissue. - - Returns - ------- - att_arr : :obj:`numpy.ndarray` - Arterial transit time array. - - Notes - ----- - This function was originally written in MATLAB by Jianxun Qu and William Tackett. - It was translated to Python by Taylor Salo. - Taylor Salo modified the code to loop over voxels, in order to account for slice timing-shifted - post-labeling delays. - - Please see https://shorturl.at/wCO56 and https://shorturl.at/aKQU3 for the original MATLAB - code. - - The code could probably be improved by operating on arrays, rather than looping over voxels. - It is also overkill for 3D acquisitions, where PLD doesn't vary by voxel. - - License - ------- - MIT License - - Copyright (c) 2023 willtack - - Permission is hereby granted, free of charge, to any person obtaining a copy - of this software and associated documentation files (the "Software"), to deal - in the Software without restriction, including without limitation the rights - to use, copy, modify, merge, publish, distribute, sublicense, and/or sell - copies of the Software, and to permit persons to whom the Software is - furnished to do so, subject to the following conditions: - - The above copyright notice and this permission notice shall be included in all - copies or substantial portions of the Software. - - THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR - IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, - FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE - AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER - LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, - OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE - SOFTWARE. - - References - ---------- - .. footbibliography:: - """ - n_voxels, n_plds = plds.shape - assert deltam_arr.shape == plds.shape, f"{deltam_arr.shape} != {plds.shape}" - - # Beginning of auxil_asl_gen_wsum - assert lds.size == n_plds, f"{lds.size} != {n_plds}" - - att_arr = np.empty(n_voxels) - for i_voxel in range(n_voxels): - deltam_by_voxel = deltam_arr[i_voxel, :] - plds_by_voxel = plds[i_voxel, :] - - # Define the possible transit times to evaluate for this voxel - transit_times = np.arange( - np.round(np.min(plds_by_voxel), 3), - np.round(np.max(plds_by_voxel), 3) + 0.001, - 0.001, - ) - - sig_sum = np.zeros((transit_times.size)) - sig_pld_sum = np.zeros((transit_times.size)) - - for j_pld in range(n_plds): - pld = plds_by_voxel[j_pld] - ld = lds[j_pld] - - # e ^ (-delta / T1a) - exp_tt_T1b = np.exp(-transit_times / t1blood) - # e ^ (-max(w - delta, 0) / T1t) - exp_pld_tt_T1t = np.exp(-np.maximum(0, pld - transit_times) / t1tissue) - # e ^ (-max(tau + w - delta, 0) / T1t) - exp_pld_ld_tt_T1t = np.exp(-np.maximum(0, pld + ld - transit_times) / t1tissue) - - # The combination of exponent terms in Equation 1 - sig = exp_tt_T1b * (exp_pld_tt_T1t - exp_pld_ld_tt_T1t) - - # The elements in Equation 1 that aren't touched here are: - # 2M0t (2 * equilibrium magnetization of brain tissue), - # Beta (a term to compensate for static tissue signal loss caused by vessel supp. - # pulses), - # alpha (labeling efficiency), - # T1t (t1tissue), - # f (CBF; perfusion rate) - # lambda (partition coefficient) - # It seems like they are cancelled out though, so it's probably fine. - - # Numerator in Equation 4, for a range of transit times - sig_pld_sum += sig * pld - # Denominator in Equation 4, for a range of transit times - sig_sum += sig - - # Predicted weighted delay values for a range of transit times - weighted_delay_predicted = sig_pld_sum / sig_sum # TT - # End of auxil_asl_gen_wsum - - # Calculate the observed weighted delay for each voxel - weighted_delay_denom = np.sum(deltam_by_voxel) - weighted_delay_num = np.sum(deltam_by_voxel * plds_by_voxel) - weighted_delay_observed = weighted_delay_num / ( - np.abs(weighted_delay_denom) + np.finfo(float).eps - ) - - # Truncate extreme transit time value to the PLD limits - weighted_delay_min = min(weighted_delay_predicted) - weighted_delay_max = max(weighted_delay_predicted) - weighted_delay_observed = np.maximum(weighted_delay_observed, weighted_delay_min) - weighted_delay_observed = np.minimum(weighted_delay_observed, weighted_delay_max) - - # Use linear interpolation to get the ATT for each weighted delay value (i.e., each voxel), - # using the predicted weighted delay and associated transit time arrays. - interp_func = interp1d(weighted_delay_predicted, transit_times) - - att_arr[i_voxel] = interp_func(weighted_delay_observed) - - return att_arr - - -def estimate_cbf_pcasl_multipld( - deltam_arr, - scaled_m0data, - plds, - tau, - labeleff, - t1blood, - t1tissue, - unit_conversion, - partition_coefficient, -): - """Estimate CBF and ATT for multi-delay PCASL data. - - Parameters - ---------- - deltam_arr : :obj:`numpy.ndarray` of shape (S, P) - Control-label values for each voxel and PLDs. - S = sample (i.e., voxel). - P = Post-labeling delay (i.e., volume). - scaled_m0data : :obj:`numpy.ndarray` of shape (S,) - The M0 volume, after scaling based on the M0-scale value. - plds : :obj:`numpy.ndarray` of shape (S, P) - Post-labeling delays. One value for each volume in ``deltam_arr``. - tau : :obj:`numpy.ndarray` of shape (P,) or (0,) - Label duration. May be a single value or may vary across volumes/PLDs. - labeleff : :obj:`float` - Estimated labeling efficiency. - t1blood : :obj:`float` - The longitudinal relaxation time of blood in seconds. - t1tissue : :obj:`float` - The longitudinal relaxation time of tissue in seconds. - unit_conversion : :obj:`float` - The factor to convert CBF units from mL/g/s to mL/ (100 g)/min. 6000. - partition_coefficient : :obj:`float` - The brain/blood partition coefficient in mL/g. Called lambda in the literature. - - Returns - ------- - att_arr : :obj:`numpy.ndarray` of shape (S,) - Arterial transit time map. - cbf : :obj:`numpy.ndarray` of shape (S,) - Estimated cerebrospinal fluid map, after estimating for each PLD and averaging across - PLDs. - - Notes - ----- - Delta-M values are first averaged over time for each unique post-labeling delay value. - - Arterial transit time is estimated on a voxel-wise basis according to - :footcite:t:`dai2012reduced`. - - CBF is then calculated for each unique PLD value using the mean delta-M values and the - estimated ATT. - - CBF is then averaged over PLDs according to :footcite:t:`juttukonda2021characterizing`, - in which an unweighted average is calculated for each voxel across all PLDs in which - PLD + tau > ATT. - - References - ---------- - .. footbibliography:: - """ - n_voxels, n_volumes = deltam_arr.shape - n_voxels_pld, n_plds = plds.shape - assert n_voxels_pld == n_voxels - first_voxel_plds = plds[0, :] - - if n_plds != n_volumes: - raise ValueError( - f"Number of PostLabelingDelays ({n_plds}) does not match " - f"number of delta-M volumes ({n_volumes})." - ) - - # Formula from Fan 2017 (equation 2) - # Determine unique original post-labeling delays (ignoring slice timing shifts) - unique_first_voxel_plds, unique_pld_idx = np.unique(first_voxel_plds, return_index=True) - unique_plds = plds[:, unique_pld_idx] # S x unique PLDs - n_unique_plds = unique_pld_idx.size - - # tau should be a 1D array, with one volume per unique PLD - if tau.size > 1: - if tau.size != n_plds: - raise ValueError( - f"Number of LabelingDurations ({tau.size}) != " - f"number of PostLabelingDelays ({n_plds})" - ) - - tau = tau[unique_pld_idx] - else: - tau = np.full(n_unique_plds, tau) - - mean_deltam_by_pld = np.zeros((n_voxels, n_unique_plds)) - for i_pld, first_voxel_pld in enumerate(unique_first_voxel_plds): - pld_idx = first_voxel_plds == first_voxel_pld - mean_deltam_by_pld[:, i_pld] = np.mean(deltam_arr[:, pld_idx], axis=1) - - # Estimate ATT for each voxel - att_arr = estimate_att_pcasl( - deltam_arr=mean_deltam_by_pld, - plds=unique_plds, - lds=tau, - t1blood=t1blood, - t1tissue=t1tissue, - ) - - # Start calculating CBF - num_factor = unit_conversion * partition_coefficient - denom_factor = 2 * labeleff * scaled_m0data * t1blood - - # Loop over PLDs and calculate CBF for each, accounting for ATT. - cbf_by_pld = np.zeros((n_voxels, n_unique_plds)) - for i_pld in range(n_unique_plds): - pld_by_voxel = unique_plds[:, i_pld] - tau_for_pld = tau[i_pld] - - pld_num_factor = num_factor * mean_deltam_by_pld[:, i_pld] * np.exp(att_arr / t1blood) - pld_denom_factor = denom_factor * ( - np.exp(-np.maximum(pld_by_voxel - att_arr, 0)) - - np.exp(-np.maximum(tau_for_pld + pld_by_voxel - att_arr, 0)) - ) - cbf_by_pld[:, i_pld] = pld_num_factor / pld_denom_factor - - # Average CBF across PLDs, but only include PLDs where PLD + tau > ATT for that voxel, - # per Juttukonda 2021 (section 2.6). - cbf = np.zeros(n_voxels) # mean CBF - for i_voxel in range(n_voxels): - plds_voxel = unique_plds[i_voxel, :] - cbf_by_pld_voxel = cbf_by_pld[i_voxel, :] - att_voxel = att_arr[i_voxel] - cbf[i_voxel] = np.mean(cbf_by_pld_voxel[(plds_voxel + tau) > att_voxel]) - - return att_arr, cbf - - def estimate_t1(metadata): """Estimate the relaxation rates of blood and gray matter based on magnetic field strength. From ffe80e198af0395ebbeb99a21ba01a0452327350 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Wed, 3 Jan 2024 11:55:18 -0500 Subject: [PATCH 37/73] Update workflows.rst --- docs/workflows.rst | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/docs/workflows.rst b/docs/workflows.rst index f89224a54..1a4ede7e5 100644 --- a/docs/workflows.rst +++ b/docs/workflows.rst @@ -482,24 +482,24 @@ Equation 3: .. math:: - &= 0 \hfilll { 0 < TI < ATT } + &= 0, \hfilll { 0 < TI < ATT } - \Delta{M}_{tiss} &= \frac{ 2 \cdot \alpha \cdot \alpha_{BS} \cdot M_{0a} \cdot CBF \cdot (TI - ATT) } - }{ 6000 \cdot e ^ {-\frac{ TI } { T_{1,b} } } \hfilll { ATT < TI < ATT + TI_{1} } + \Delta{M}_{tiss} &= \frac{ 2 \cdot \alpha \cdot \alpha_{BS} \cdot M_{0a} \cdot CBF \cdot (TI - ATT)} + { 6000 \cdot e ^ { -\frac{ TI } { T_{1,b} } }, \hfilll { ATT < TI < ATT + TI_{1} } &= \frac{ 2 \cdot \alpha \cdot \alpha_{BS} \cdot M_{0a} \cdot CBF \cdot TI_{1} } - { 6000 \cdot e ^ { \frac{ TI } { T_{1,b} } } } \hfilll { ATT + TI_{1} < TI } + { 6000 \cdot e ^ { \frac{ TI } { T_{1,b} } } }, \hfilll { ATT + TI_{1} < TI } Equation 5: .. math:: - &= 0 \hfilll{ 0 < TI < aBAT } + &= 0, \hfilll{ 0 < TI < aBAT } \Delta{M}_{art} &= 2 \cdot \alpha \cdot \alpha_{BS} \cdot M_{0b} \cdot aBV \cdot - e ^ { -\frac{ TI } { T_{1,b} } } \hfilll{ aBAT < TI < aBAT + TI_{1} } + e ^ { -\frac{ TI } { T_{1,b} } }, \hfilll{ aBAT < TI < aBAT + TI_{1} } - &= 0 \hfilll{ aBAT + TI_{1} < TI } + &= 0, \hfilll{ aBAT + TI_{1} < TI } Additional Denoising Options From 2c6ac423ee4e44a0d269a72ef10a4ed1aadf21cd Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Wed, 3 Jan 2024 12:14:59 -0500 Subject: [PATCH 38/73] Update workflows.rst --- docs/workflows.rst | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/docs/workflows.rst b/docs/workflows.rst index 1a4ede7e5..61f704f00 100644 --- a/docs/workflows.rst +++ b/docs/workflows.rst @@ -482,13 +482,15 @@ Equation 3: .. math:: - &= 0, \hfilll { 0 < TI < ATT } + + &= 0, && { 0 < TI < ATT } \\ \Delta{M}_{tiss} &= \frac{ 2 \cdot \alpha \cdot \alpha_{BS} \cdot M_{0a} \cdot CBF \cdot (TI - ATT)} - { 6000 \cdot e ^ { -\frac{ TI } { T_{1,b} } }, \hfilll { ATT < TI < ATT + TI_{1} } + { 6000 \cdot e ^ { -\frac{ TI } { T_{1,b} } }, && { ATT < TI < ATT + TI_{1} } \\ &= \frac{ 2 \cdot \alpha \cdot \alpha_{BS} \cdot M_{0a} \cdot CBF \cdot TI_{1} } - { 6000 \cdot e ^ { \frac{ TI } { T_{1,b} } } }, \hfilll { ATT + TI_{1} < TI } + { 6000 \cdot e ^ { \frac{ TI } { T_{1,b} } } }, && { ATT + TI_{1} < TI } + Equation 5: From f3e9b58501dbbaf915d7b2edc65e48a69a0b8ed1 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Wed, 3 Jan 2024 12:24:46 -0500 Subject: [PATCH 39/73] Update workflows.rst --- docs/workflows.rst | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/docs/workflows.rst b/docs/workflows.rst index 61f704f00..31b78301d 100644 --- a/docs/workflows.rst +++ b/docs/workflows.rst @@ -446,17 +446,17 @@ Equation 2: .. math:: - &= 0 \hfilll { 0 < LD + PLD < ATT } + &= 0 && { 0 < LD + PLD < ATT } \\ \Delta{M}_{tiss} &= \frac{ 2 \cdot \alpha \cdot \alpha_{BS} \cdot T_{1b} \cdot M_{0a} \cdot CBF \cdot \left[ 1 - e ^ {-\frac{ LD + PLD - ATT } { T_{1,b} }} - \right] }{ 6000 \cdot e ^ { \frac{ ATT } { T_{1,b} } } } \hfilll { ATT < LD + PLD < ATT + LD } + \right] }{ 6000 \cdot e ^ { \frac{ ATT } { T_{1,b} } } } && { ATT < LD + PLD < ATT + LD } \\ &= \frac{ 2 \cdot \alpha \cdot \alpha_{BS} \cdot T_{1b} \cdot M_{0a} \cdot CBF \cdot \left[ 1 - e ^ {-\frac{ LD } { T_{1,b} }} - \right] }{ 6000 \cdot e ^ { \frac{ PLD } { T_{1,b} } } } \hfilll { ATT < PLD } + \right] }{ 6000 \cdot e ^ { \frac{ PLD } { T_{1,b} } } } && { ATT < PLD } Equation 4: @@ -482,7 +482,6 @@ Equation 3: .. math:: - &= 0, && { 0 < TI < ATT } \\ \Delta{M}_{tiss} &= \frac{ 2 \cdot \alpha \cdot \alpha_{BS} \cdot M_{0a} \cdot CBF \cdot (TI - ATT)} From b16730c9c12f1573e8240d5edec0ad2d44760c4a Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Wed, 3 Jan 2024 12:25:13 -0500 Subject: [PATCH 40/73] Update workflows.rst --- docs/workflows.rst | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/workflows.rst b/docs/workflows.rst index 31b78301d..d553e291e 100644 --- a/docs/workflows.rst +++ b/docs/workflows.rst @@ -462,12 +462,12 @@ Equation 4: .. math:: - &= 0 \hfilll{ 0 < LD + PLD < aBAT } + &= 0 && { 0 < LD + PLD < aBAT } \Delta{M}_{art} &= 2 \cdot \alpha \cdot \alpha_{BS} \cdot M_{0b} \cdot aBV \cdot - e ^ { -\frac{ aBAT } { T_{1,b} } } \hfilll{ aBAT < LD + PLD < aBAT + LD } + e ^ { -\frac{ aBAT } { T_{1,b} } } && { aBAT < LD + PLD < aBAT + LD } - &= 0 \hfilll{ aBAT < PLD } + &= 0 && { aBAT < PLD } Pulsed ASL From 21098ab639630f8a088ba5d94c163c189a5abf53 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Wed, 3 Jan 2024 12:32:57 -0500 Subject: [PATCH 41/73] Update workflows.rst --- docs/workflows.rst | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/docs/workflows.rst b/docs/workflows.rst index d553e291e..a5abe35c0 100644 --- a/docs/workflows.rst +++ b/docs/workflows.rst @@ -446,12 +446,12 @@ Equation 2: .. math:: - &= 0 && { 0 < LD + PLD < ATT } \\ + &= 0 && { 0 < LD + PLD < ATT } \Delta{M}_{tiss} &= \frac{ 2 \cdot \alpha \cdot \alpha_{BS} \cdot T_{1b} \cdot M_{0a} \cdot CBF \cdot \left[ 1 - e ^ {-\frac{ LD + PLD - ATT } { T_{1,b} }} - \right] }{ 6000 \cdot e ^ { \frac{ ATT } { T_{1,b} } } } && { ATT < LD + PLD < ATT + LD } \\ + \right] }{ 6000 \cdot e ^ { \frac{ ATT } { T_{1,b} } } } && { ATT < LD + PLD < ATT + LD } &= \frac{ 2 \cdot \alpha \cdot \alpha_{BS} \cdot T_{1b} \cdot M_{0a} \cdot CBF \cdot \left[ @@ -482,25 +482,25 @@ Equation 3: .. math:: - &= 0, && { 0 < TI < ATT } \\ + &= 0 && { 0 < TI < ATT } - \Delta{M}_{tiss} &= \frac{ 2 \cdot \alpha \cdot \alpha_{BS} \cdot M_{0a} \cdot CBF \cdot (TI - ATT)} - { 6000 \cdot e ^ { -\frac{ TI } { T_{1,b} } }, && { ATT < TI < ATT + TI_{1} } \\ + \Delta{M}_{tiss} &= \frac{ 2 \cdot \alpha \cdot \alpha_{BS} \cdot M_{0a} \cdot CBF \cdot (TI - ATT) } + { 6000 \cdot e ^ { -\frac{ TI } { T_{1,b} } } && { ATT < TI < ATT + TI_{1} } &= \frac{ 2 \cdot \alpha \cdot \alpha_{BS} \cdot M_{0a} \cdot CBF \cdot TI_{1} } - { 6000 \cdot e ^ { \frac{ TI } { T_{1,b} } } }, && { ATT + TI_{1} < TI } + { 6000 \cdot e ^ { \frac{ TI } { T_{1,b} } } } && { ATT + TI_{1} < TI } Equation 5: .. math:: - &= 0, \hfilll{ 0 < TI < aBAT } + &= 0 && { 0 < TI < aBAT } \Delta{M}_{art} &= 2 \cdot \alpha \cdot \alpha_{BS} \cdot M_{0b} \cdot aBV \cdot - e ^ { -\frac{ TI } { T_{1,b} } }, \hfilll{ aBAT < TI < aBAT + TI_{1} } + e ^ { -\frac{ TI } { T_{1,b} } } && { aBAT < TI < aBAT + TI_{1} } - &= 0, \hfilll{ aBAT + TI_{1} < TI } + &= 0 && { aBAT + TI_{1} < TI } Additional Denoising Options From fa2f7261103d9a900309223553f2a87bea32304b Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Wed, 3 Jan 2024 12:42:52 -0500 Subject: [PATCH 42/73] Update workflows.rst --- docs/workflows.rst | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/docs/workflows.rst b/docs/workflows.rst index a5abe35c0..4e175c3f9 100644 --- a/docs/workflows.rst +++ b/docs/workflows.rst @@ -484,11 +484,11 @@ Equation 3: &= 0 && { 0 < TI < ATT } - \Delta{M}_{tiss} &= \frac{ 2 \cdot \alpha \cdot \alpha_{BS} \cdot M_{0a} \cdot CBF \cdot (TI - ATT) } - { 6000 \cdot e ^ { -\frac{ TI } { T_{1,b} } } && { ATT < TI < ATT + TI_{1} } + \Delta{M}_{tiss} &= \frac{ 2 \cdot \alpha \cdot \alpha_{BS} \cdot M_{0a} \cdot CBF \cdot (TI - ATT) + }{ 6000 \cdot e ^ { -\frac{ TI }{ T_{1,b} } } } && { ATT < TI < ATT + TI_{1} } - &= \frac{ 2 \cdot \alpha \cdot \alpha_{BS} \cdot M_{0a} \cdot CBF \cdot TI_{1} } - { 6000 \cdot e ^ { \frac{ TI } { T_{1,b} } } } && { ATT + TI_{1} < TI } + &= \frac{ 2 \cdot \alpha \cdot \alpha_{BS} \cdot M_{0a} \cdot CBF \cdot TI_{1} + }{ 6000 \cdot e ^ { \frac{ TI }{ T_{1,b} } } } && { ATT + TI_{1} < TI } Equation 5: From b344acd476b237a53ea96fee9e02f08c6ca3e1ef Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Wed, 3 Jan 2024 13:15:16 -0500 Subject: [PATCH 43/73] Update workflows.rst --- docs/workflows.rst | 39 +++++++++++++++++++++++++++++---------- 1 file changed, 29 insertions(+), 10 deletions(-) diff --git a/docs/workflows.rst b/docs/workflows.rst index 4e175c3f9..0ef66a86c 100644 --- a/docs/workflows.rst +++ b/docs/workflows.rst @@ -430,21 +430,34 @@ Multi-Delay ASL =============== In multi-delay ASL, control-label pairs are acquired for multiple post-labeling delay values. -This type of acquisition requires more complicated models, but it also results in more accurate -CBF estimates. -Also, multi-delay ASL allows for the estimation of arterial transit time (ATT). +This type of acquisition requires more complicated models, +but it also results in more accurate CBF estimates. +Also, multi-delay ASL allows for the estimation of arterial transit time (ATT), +arterial blood arrival time (aBAT), and arterial blood volume (aBV). + +For multi-delay PCASL and Q2TIPS/QUIPSSII PASL data, +*ASLPrep* uses :func:`scipy.optimize.curvefit` to estimate +CBF, ATT, aBAT, and aBV values for each voxel according to the two-compartment general kinetic +model described in :footcite:t:`woods2023recommendations`. + +In the two-compartment model, the difference signal (:math:`\Delta{M}`) is decomposed into +a tissue component (:math:`\Delta{M}_{tiss}`) and +a macrovascular component (:math:`\Delta{M}_{art}`). +The combined macrovascular and tissue signal is simply +:math:`\Delta{M} = \Delta{M}_{tiss} + \Delta{M}_{art}`. Pseudo-Continuous ASL --------------------- -For multi-delay PCASL data, *ASLPrep* uses :func:`scipy.optimize.curvefit` to estimate -CBF, ATT, aBAT, and aBV values for each voxel using Equations 2 and 4 in -:footcite:t:`woods2023recommendations`. +The tissue component of the :math:`\Delta{M}` signal is estimated according to +:eq:`woods_equation_02`, while the macrovascular (i.e., arterial) component is estimated using +:eq:`woods_equation_04`. Equation 2: .. math:: + :label: woods_equation_02 &= 0 && { 0 < LD + PLD < ATT } @@ -461,6 +474,7 @@ Equation 2: Equation 4: .. math:: + :label: woods_equation_04 &= 0 && { 0 < LD + PLD < aBAT } @@ -473,14 +487,18 @@ Equation 4: Pulsed ASL ---------- -For multi-delay Q2TIPS or QUIPSSII PCASL data, -*ASLPrep* uses :func:`scipy.optimize.curvefit` to estimate -CBF, ATT, aBAT, and aBV values for each voxel using Equations 3 and 5 in -:footcite:t:`woods2023recommendations`. +The tissue component of the :math:`\Delta{M}` signal is estimated according to +:eq:`woods_equation_03`, while the macrovascular (i.e., arterial) component is estimated using +:eq:`woods_equation_05`. + +.. warning:: + + Multi-delay QUIPSS PASL is not supported. Equation 3: .. math:: + :label: woods_equation_03 &= 0 && { 0 < TI < ATT } @@ -494,6 +512,7 @@ Equation 3: Equation 5: .. math:: + :label: woods_equation_05 &= 0 && { 0 < TI < aBAT } From bd03dffd8f112a0f23fd23e217fefa74d5a36071 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Wed, 3 Jan 2024 13:16:05 -0500 Subject: [PATCH 44/73] Fix tests. --- LICENSE.md | 27 ----------------------- aslprep/tests/test_utils_cbf.py | 39 --------------------------------- 2 files changed, 66 deletions(-) diff --git a/LICENSE.md b/LICENSE.md index c6a850dc5..f68c25e2d 100644 --- a/LICENSE.md +++ b/LICENSE.md @@ -25,30 +25,3 @@ SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - - -## External code - -### aslprep.utils.misc.estimate_att_pcasl - -MIT License - -Copyright (c) 2023 willtack - -Permission is hereby granted, free of charge, to any person obtaining a copy -of this software and associated documentation files (the "Software"), to deal -in the Software without restriction, including without limitation the rights -to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -copies of the Software, and to permit persons to whom the Software is -furnished to do so, subject to the following conditions: - -The above copyright notice and this permission notice shall be included in all -copies or substantial portions of the Software. - -THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -SOFTWARE. diff --git a/aslprep/tests/test_utils_cbf.py b/aslprep/tests/test_utils_cbf.py index 132c5361b..a660edae2 100644 --- a/aslprep/tests/test_utils_cbf.py +++ b/aslprep/tests/test_utils_cbf.py @@ -1,40 +1 @@ """Tests for the aslprep.utils.cbf module.""" -import numpy as np - -from aslprep.utils import cbf - - -def test_estimate_att_pcasl_3d(): - """Smoke test aslprep.utils.cbf.estimate_att_pcasl with 3D data (no slice timing).""" - n_voxels, n_unique_plds = 1000, 6 - base_plds = np.linspace(0.25, 1.5, n_unique_plds)[None, :] - plds = np.repeat(base_plds, n_voxels, axis=0) - deltam_arr = np.random.random((n_voxels, n_unique_plds)) - tau = np.full(n_unique_plds, 1.4) - att = cbf.estimate_att_pcasl( - deltam_arr=deltam_arr, - plds=plds, - lds=tau, - t1blood=1.6, - t1tissue=1.3, - ) - assert att.shape == (n_voxels,) - - -def test_estimate_att_pcasl_2d(): - """Smoke test aslprep.utils.cbf.estimate_att_pcasl with slice-shifted PLDs.""" - n_voxels, n_unique_plds, n_slice_times = 1000, 6, 50 - base_plds = np.linspace(0.25, 1.5, n_unique_plds) - slice_times = np.linspace(0, 3, n_slice_times) - slice_shifted_plds = slice_times[:, None] + base_plds[None, :] - plds = np.repeat(slice_shifted_plds, n_voxels // n_slice_times, axis=0) - deltam_arr = np.random.random((n_voxels, n_unique_plds)) - tau = np.full(n_unique_plds, 1.4) - att = cbf.estimate_att_pcasl( - deltam_arr=deltam_arr, - plds=plds, - lds=tau, - t1blood=1.6, - t1tissue=1.3, - ) - assert att.shape == (n_voxels,) From e069608248090be889ec1ebb7dbe23df9e889f69 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Wed, 3 Jan 2024 14:35:40 -0500 Subject: [PATCH 45/73] Update again. --- aslprep/interfaces/cbf.py | 14 +++++++------- aslprep/utils/cbf.py | 33 ++++++++++++++++++++++---------- aslprep/workflows/asl/outputs.py | 9 ++++----- docs/outputs.rst | 4 +++- 4 files changed, 37 insertions(+), 23 deletions(-) diff --git a/aslprep/interfaces/cbf.py b/aslprep/interfaces/cbf.py index 414dedff2..dc7d5aa50 100644 --- a/aslprep/interfaces/cbf.py +++ b/aslprep/interfaces/cbf.py @@ -502,16 +502,16 @@ def _run_interface(self, runtime): self._results["plds"] = pld_file elif is_multi_pld: - # Broadcast PLDs to voxels by PLDs + # Broadcast PLDs to voxels by PLDs, even though there's no slice timing to account for. plds = np.dot(plds[:, None], np.ones((1, deltam_arr.shape[0]))).T if is_multi_pld: ti1 = None tau = None - if is_casl: + if is_casl: # (P)CASL needs tau, but not ti1 tau = metadata["LabelingDuration"] - else: + else: # PASL needs ti1, but not tau if metadata["BolusCutOffTechnique"] == "QUIPSSII": # PASL + QUIPSSII # Only one BolusCutOffDelayTime allowed. @@ -531,10 +531,10 @@ def _run_interface(self, runtime): ) cbf, att, abat, abv = fit_deltam_multipld( - deltam_arr, - scaled_m0data, - plds, - labeleff, + deltam_arr=deltam_arr, + scaled_m0data=scaled_m0data, + plds=plds, + labeleff=labeleff, t1blood=t1blood, partition_coefficient=PARTITION_COEF, is_casl=is_casl, diff --git a/aslprep/utils/cbf.py b/aslprep/utils/cbf.py index 60364964f..b26ab3aac 100644 --- a/aslprep/utils/cbf.py +++ b/aslprep/utils/cbf.py @@ -708,8 +708,10 @@ def calculate_deltam_pasl(X, cbf, att, abat, abv): Equilibrium magnetization of arterial blood. Used for deltam_art calculation. tis : :obj:`numpy.ndarray` + Inversion times in seconds. ti1 : :obj:`float` - Bolus cutoff delay time. + Bolus cutoff delay time. For Q2TIPS, this is the first delay time. + For QUIPSSII, it's the only one. """ # float parameters labeleff = X[0, 0] @@ -768,17 +770,28 @@ def fit_deltam_multipld( Parameters ---------- - deltam - plds - tau : :obj:`float` or :obj:`numpy.ndarray` - Label duration. May be a single value or may vary across volumes/PLDs. - t1blood + deltam_arr : :obj:`numpy.ndarray` of shape (n_voxels, n_volumes) + Subtracted (control - label) signal. + scaled_m0data : :obj:`numpy.ndarray` of shape (n_voxels,) + M0 data, after any scaling requested by the user. + plds : :obj:`numpy.ndarray` of shape (n_voxels, n_volumes) + Post-labeling delay values. This uses the TI values for PASL data. labeleff : :obj:`float` Labeling efficiency, also referred to as alpha. In the case of background suppression, this combines alpha and alpha_bs, so we don't need a separate term for alpha_bs. - scaled_m0data - partition_coefficient + t1blood : :obj:`float` + Longitudinal relaxation time of arterial blood in seconds. + Used for both deltam_tiss and deltam_art calculation. + partition_coefficient : :obj:`float` + Brain partition coefficient (lambda in Alsop 2015). Generally 0.9. + is_casl : :obj:`bool` + True if is (P)CASL sequence. + tau : :obj:`float` or :obj:`numpy.ndarray` or None, optional + Label duration. May be a single value or may vary across volumes/PLDs. + Only defined for (P)CASL data. + ti1 : :obj:`float` or None, optional + TI1. Only defined for PASL data. Returns ------- @@ -867,9 +880,9 @@ def fit_deltam_multipld( # Upper bounds provided by Manuel Taso bounds=( (0, 0, 0, 0), - # Manuel Taso recommended 5 and 5 for ATT and aBAT, + # Manuel Taso recommended 5, 5, 0.2 for ATT, aBAT, and aBV, # but our test data maxed out around 12 when left unbounded. - (300, 15, 15, 0.1), + (300, 15, 15, 1), ), )[0] diff --git a/aslprep/workflows/asl/outputs.py b/aslprep/workflows/asl/outputs.py index 32527d8ea..ebd9299f1 100644 --- a/aslprep/workflows/asl/outputs.py +++ b/aslprep/workflows/asl/outputs.py @@ -34,12 +34,15 @@ }, "att": { "suffix": "att", + "Units": "s", }, "abat": { "suffix": "abat", + "Units": "s", }, "abv": { "suffix": "abv", + "Units": "fraction", }, # SCORE/SCRUB outputs "cbf_ts_score": { @@ -70,6 +73,7 @@ "att_basil": { "desc": "basil", "suffix": "att", + "Units": "s", }, } @@ -572,10 +576,6 @@ def init_ds_asl_native_wf( workflow.connect([(inputnode, ds_cbf, [(cbf_name, "in_file")])]) for att_name in att: - # TODO: Add EstimationReference and EstimationAlgorithm - att_meta = { - "Units": "s", - } fields = BASE_INPUT_FIELDS[att_name] ds_att = pe.Node( @@ -584,7 +584,6 @@ def init_ds_asl_native_wf( compress=True, dismiss_entities=("echo",), **fields, - **att_meta, ), name=f"ds_{att_name}", run_without_submitting=True, diff --git a/docs/outputs.rst b/docs/outputs.rst index 749fd598a..52b2d591b 100644 --- a/docs/outputs.rst +++ b/docs/outputs.rst @@ -322,6 +322,8 @@ CBF Outputs:: [_space-