Abstract
Purpose To develop a generalized signal model for dual-module velocity-selective arterial spin labeling (dm-VSASL) that can integrate arbitrary saturation and inversion profiles.
Theory and Methods A recently developed mathematical framework for single-module VSASL is extended to address the increased complexity of dm-VSASL and to model the use of realistic velocity-selective profiles in the label-control and vascular crushing modules. Expressions for magnetization difference, arterial delivery functions, labeling efficiency, and cerebral blood flow (CBF) estimation error are presented. Sources of error are examined and timing requirements to minimize quantification errors are derived.
Results For ideal velocity-selective profiles, the predicted signals match those of prior work. With realistic profiles, a CBF-dependent estimation error can occur when velocity-selective inversion (VSI) is used for the labeling modules and velocity-selective saturation (VSS) is used for the vascular crushing module. The error reflects a mismatch between the leading and trailing edges of the delivery function for the second bolus and can be minimized by choosing a nominal labeling cutoff velocity that is lower than the nominal saturation cutoff velocity. In the presence of B0 and B1 inhomogeneities, the labeling efficiency of dual-module VSI is more attenuated than that of dual-module VSS.
Conclusion The proposed signal model will enable researchers to more accurately assess and compare the performance of realistic dm-VSASL implementations and improve the quantification of dm-VSASL CBF measures.
1 INTRODUCTION
Dual-module velocity-selective arterial spin labeling (dm-VSASL) has been shown to offer improvements in robustness and SNR efficiency as compared to singlemodule VSASL1,2,3,4,5. The primary gain in SNR efficiency has been shown to reflect the creation of a larger magnetization difference between the control and label conditions through re-labeling of spins by the second module1,2,3. When using modules that invert static tissue, dm-VSASL offers inherent background suppression that can further improve performance and robustness2,6. Although typically negligible, there can also be an SNR efficiency gain due to an increase in bolus duration that reflects labeling by the second module of arterial spins that have entered the coverage of the RF transmit coil after the first module1,2,3.
Prior dm-VSASL signal models have assumed ideal saturation and inversion profiles and are not easily extended to handle realistic profiles. Here we build upon a recently introduced mathematical model for singlemodule VSASL7 to develop a generalized signal model for dm-VSASL that can readily integrate realistic saturation and inversion profiles and thus enable more accurate quantification of cerebral blood flow (CBF).
As the Theory section is rather involved, we start with a guide to the theoretical subsections. In Section 2.1, we review the basic structure and timing of a typical dm-VSASL scan and present the signal model. Then, in Section 2.2, we adopt the concept of passband and saturation reference functions introduced in Ref. [7] to rewrite the signal model in a form that makes explicit the dependence on the amplitudes and shapes of the label-control and saturation functions. In Section 2.3 we demonstrate the application of the framework with ideal rect profiles and compare the resulting expressions with those from prior work2. We present expressions for CBF quantification error and labeling efficiency in Sections 2.4 and 2.5, respectively, and then go on to show the dependence of the CBF quantification error on effective bolus width errors in Section 2.6 and present conditions for minimizing the errors in Section 2.7. Finally, we briefly consider cerebral blood volume effects in Section 2.8. Analyses using realistic Bloch-simulated profiles and taking into account CBF levels, varying timing parameters, and B0 and B1 inhomogeneities are described in the Methods and Results sections, and areas for future work are addressed in the Discussion.
2 THEORY
2.1 Signal Model
We consider the dual-module VSASL diagram shown in Figure 1, with two label-control modules (LCM1 and LCM2), a vascular crushing module (VCM), and timing parameters Tsat, τ1, τ2, and PLD. To simplify the presentation, we assume that LCM1, LCM2, and the VCM are applied instantaneously at t = 0, τ1, and τ1 + τ2, respectively, such that TI = τ1 + τ2 + PLD. The labeling and control profiles are l1(v) and c1(v), respectively, for LCM1, and l2(v) and c2(v) for LCM2, where v denotes velocity. The labeling profiles can represent any velocity-selective implementation, including velocity-selective inversion (VSI), velocity-selective saturation (VSS), or velocity-selective saturation with static spins inverted (VSSI)2. The VCM saturation response is represented by s(v).
For dmVSI and mixed variants (e.g. VSSI+VSS or VSI+VSS), LCM2 uses l2(v) and c2(v) in the control and label conditions, respectively. For dmVSS, LCM2 uses c2(v) and l2(v) in the control and label conditions, respectively.
We focus here on dual-module VSI (dmVSI) and dual-module VSS (dmVSS), as these are currently the most widely used dm-VSASL implementations. In LCM1, c1(v) and l1(v) are applied during the control and label conditions, respectively. For dmVSI, LCM2 uses l2(v) and c2(v) in the control and label conditions, respectively, whereas for dmVSS, c2(v) and l2(v) are applied in the control and label conditions, respectively. The expressions presented below for dmVSI also apply to other combinations2 (which we will refer to as mixed variants), such as VSSI+VSS and VSI+VSS, that swap the label and control responses between LCM1 and LCM2. A glossary of key functions and variables is provided in Table 1.
To simplify the presentation we have omitted subscripts that can be used to indicate different values of the functions and variables for LCM1 and LCM2, e.g. subscripts can be added to l(v) to denote l1(v) and l2(v).
Building on the framework presented in Ref. [7], it is shown in the Appendix that the cumulative magnetization difference (control-label) as a function of time for t > τ1 + τ2 is given by ΔMz(t) = ΔMz,1(t) + ΔMz,2(t) where
and
with longitudinal relaxation weighting functions
and M0b and T1b denoting the equilibrium magnetization and longitudinal relaxation time constant, respectively, of arterial blood. Equivalent arterial delivery functions are defined as
where
and
with Q0 denoting CBF. The function va(t) denotes the initial velocity (at t = 0) of arterial blood that will decelerate to a boundary velocity vb at a later time t, which was previously denoted as va(t, vb) in Ref. [7]. To simplify the presentation, we have made the dependence on vb implicit here, such that va(t) = va(t, vb).
With these definitions, da,1(t) and da,2(t) represent the arterial blood flow that arrives at time t to the boundary velocity vb, where da,1(t) reflects the blood delivered in the bolus created by LCM1, LCM2, and the VCM and da,2(t) reflects the blood delivered in the bolus created by LCM2 and the VCM. We use the term equivalent arterial delivery function to reflect the fact that blood that has already been delivered at some time t can be affected by pulses applied at a later time7, and the effect of the “future” pulses is included in the definition of the delivery function. For example, the value of da,1(t) at some time 0 ≤ t < τ1 + τ2 reflects the application of the VCM at a future time t = τ1 + τ2.
2.2 Passband and Saturation Reference Functions
Following the approach described in Ref. [7], we define the passband function for the ith LCM module as 
where Zi is 1 for modules that invert static spins and 0 otherwise. It is useful to express the passband function as the product pi(v) = Pa,i · pref,i(v) of an amplitude term Pa,i that denotes the steady-state value of the passband function and a reference passband function pref,i(v) that starts at zero and approaches 1.0 for velocities above a labeling cut-off velocity vl,i. Similarly, it is useful to express the VCM saturation function as the product s(v) = S0 · sref(v) of an amplitude term S0 = s(0) and a reference saturation function sref(v) that starts at sref(0) = 1.0 and approaches zero for velocities above a VCM cutoff velocity vc.
To facilitate comparison with prior work2, it is helpful to express the passband amplitude as the product Pa,i = Pa⋆,i · κi · βi where Pa⋆,i denotes the ideal steadystate value (1.0 and 2.0 for VSS and VSI, respectively), κi = |li(0)| = |ci(0)| ≤ 1.0 represents the absolute value of the control and label functions at v = 0 and reflects factors (such as T1, T2, and B1 effects) that equally affect both the label and control, and βi ≥ 0 reflects remaining factors that are not accounted for by κi, such as the steady-state value of the VSI labeling profile7.
With these definitions, we can immediately write
where ζ = 0 for dmVSS and ζ = 1 for dmVSI and mixed variants. With the additional assumption2 that ci(v) = ci(0), we can write for dmVSI (and mixed variants)
For a typical dmVSI implementation where the LCM modules use the same train of pulses this can be written as
where pref(v) is an index-independent reference passband function and κ and β are index-independent scaling terms.
For dmVSS, we have
which for a typical implementation with identical LCM modules reduces to
In practice, pref(va(t−τ1))pref(va(t)) ≈ pref(va(t−τ1)) for τ1 > Δtl, where Δtl is the LCM transit delay7, resulting in
2.3 Ideal Profiles
We first consider delivery functions obtained with ideal passband pref(v) = 1 − rect(v/(2vl)) and saturation sref(v) = rect(v/(2vc)) functions, where vl and vc are the labeling and saturation cutoff velocities7, respectively, and assumed to be equal unless otherwise noted. Examples of da,1(t) and da,2(t) are shown in Figs. 2A and 3A for dmVSI and dmVSS, respectively, with Q0 = 1, τ1 = 1.5 s, τ2 = 0.5 s, vb = 0.1 cm/s, vl = vc = 2 cm/s, and κ = β = S0 = 1.0. Additional examples with realistic Bloch-simulated profiles are discussed in the Results section.
(A, C, E) and cumulative magnetization differences (B, D, F) for ideal rect profile (top row) and Bloch-simulated sinc FTVSI nominal profiles (second and third rows) with Q0 = 1, τ1 = 1.5 s, τ2 = 0.5 s, vb = 0.1 cm/s, vc = 2 cm/s, vl = 2 cm/s (top 2 rows) or vl = 1.2 cm/s (bottom row), and Δtl and Δtc denoting the LCM and VCM transit delays7. For the top two rows, these transit delays are equal and indicate the time needed for arterial blood to decelerate from 2 cm/s to 0.1 cm/s. For the bottom row, Δtl < Δtc since vl < vc, with Δtl indicating the time for arterial blood to decelerate from 1.2 cm/s to 0.1 cm/s. Note that area under the curve of da,2(t) for vl = 1.2 cm/s (panel E) is greater than that for vl = 2.0 cm/s (panel C), reflecting the better matching of the leading and trailing edges in the delivery function for the second bolus. This results in greater amplitudes in the ΔMz,2(t) and ΔMz(t) curves for vl = 1.2 cm/s. The delivery time t′ and measurement time t variables are described in the Appendix. Curves for ΔMz,1(t), ΔMz,2(t), and ΔMz (t) are depicted for t ≥ 0, but reflect what is actually measured only for the period after application of the VCM, i.e. t ≥ τ1 + τ2. For t < τ1 + τ2, they reflect the “virtual” magnetization difference that would be measured if the effects of future events (e.g. LCM2 and VCM) could affect prior periods.
(A, C) and cumulative magnetization differences (B,D) for ideal rect profile (top row) and Bloch-simulated BIR8 nominal profile (bottom row) with Q0 = 1, τ1 = 1.5 s, τ2 = 0.5 s, vb = 0.1 cm/s, vl = vc = 2 cm/s, and Δtl and Δtc denoting the labeling and VCM transit delays7, which are equal for this example and indicate the time needed for arterial blood to decelerate from 2 cm/s to 0.1 cm/s. The delivery time t′ and measurement time t variables are described in the Appendix. Curves for ΔMz,1(t), ΔMz,2(t), and ΔMz (t) are depicted for t ≥ 0, but reflect what is actually measured only for the period after application of the VCM, i.e. t ≥ τ1 + τ2. For t < τ1 + τ2, they reflect the “virtual” magnetization difference that would be measured if the effects of future events (e.g. LCM2 and VCM) could affect prior periods.
For dmVSI, da,1(t) and da,2(t) in Fig. 2A correspond to ideal boluses with amplitudes of 2Q0 spanning the intervals [Δtl, τ1 + Δtl] and [τ1 + Δtl, τ1 + τ2 + Δtc], respectively, with corresponding bolus areas of 2Q0τ1 and 2Q0τ2, where Δtc is the VCM transit delay7 and Δtl and Δtc are equal for this example and represent the time for arterial blood to decelerate from vl = vc = 2 cm/s to 0.1 cm/s. The corresponding cumulative magnetization differences, ΔMz,1(t) and ΔMz,2(t), and the sum ΔMz(t) are shown in Fig. 2B, where all spins are delivered by t = τ1 + τ2 +PLD as the PLD is set equal to Δtc. The selection of PLD is discussed in detail in Section 2.7. For general values of S0, κ, and β, the bolus areas of da,1(t) and da,2(t) are 2Q0S0κ2βτ1 and 2Q0S0κβτ2. If 
(where eTEVS and T2b denote the effective echo time and T2 of arterial blood), then the associated expressions for ΔMz,1(t) and ΔMz,2(t) match those for dmVSI groups 1 and 2, respectively, in Ref. [2]. For dmVSS, da,1(t) and da,2(t) in Fig. 3A correspond to ideal boluses with amplitudes of Q0 spanning the intervals [Δtl, τ1 + τ2 + Δtl] and [τ1 + Δtl, τ1 + τ2 + Δtc], respectively, with corresponding bolus areas of Q0(τ1 + τ2) and Q0τ2. The corresponding cumulative magnetization differences, ΔMz,1(t) and ΔMz,2(t), and the sum ΔMz(t) are shown in Fig. 3B. For general values of S0, κ, and β, the bolus areas of da,1(t) and da,2(t) are Q0S0κ2β(τ1 + τ2 + (1 − β)τ2) and Q0S0κβτ2, respectively. Comparing to Ref. [2], it can be shown that the associated expression for ΔMz,2(t) matches the first two terms of the dmVSS group 2 expression in the prior work, while the expression for ΔMz,1(t) contains a term that matches that of group 1 in the interval [Δtl, τ1 + Δtl] and another term that nearly matches the prior work’s remaining dmVSS group 2 terms in the interval [τ1 + Δtl, τ1 + τ2 + Δtc]. The discrepancy stems from a minor error in the prior derivation2, which is corrected in the current work.
2.4 CBF Quantification
An unbiased estimate of Q0 based on a measurement of the magnetization difference at t = TI can be written as
where
denote effective areas. If the estimate uses values 
and 
that differ from the actual values A1,eff and A2,eff, then the resulting estimate 
is biased. The fractional error of the biased estimate can be written as
where
are the error components for the first and second boluses and
2.5 Labeling Efficiency
In the VSASL white paper3 (n.b. erratum), the difference signal for single-module VSASL has the form
where α and τ were defined as the single-module VSASL labeling efficiency and bolus width, respectively. Comparing this to Eqn. 15, setting τ = τ1 + τ2, and equating like terms leads to
as the expression for dm-VSASL labeling efficiency.
For dmVSI with ideal rect profiles,
and for dmVSS with ideal rect profiles,
where we have assumed that the timing requirements specified below in Section 2.7 are satisfied. There do not appear to be analytical expressions for finding the maxima of αdmVSI,rect and αdmVSS,rect, but they are straightforward to evaluate and optimize numerically.
Examples of ideal rect labeling efficiency versus τ1 are shown by the dashed green lines in Figs. 4A and 5A assuming S0 = κ = 1.0, τ1 + τ2 = 2.0 s, Tsat = 2.5 s, T1b = 1.66 s, and β = 0.88 for dmVSI and β = 1.0 for dmVSS, where the dmVSI β = Pa/2 value corresponds to the amplitude Pa = 1.76 of the sinc Fourier Transform VSI (FTVSI) passband function first examined in Ref. [7] and used later in this paper. Note that since α is referenced to the performance of single-module VSASL (with a maximum value of 1.0 for ideal VSI), the additional labeling of spins provided by dm-VSASL1,2can result in α > 1.0 under ideal conditions, as shown in Fig. 4A.
(A,B) dmVSI labeling efficiency and CBF estimation error as a function of τ1 for sinc FTVSI LCM (vl = 1.2 or 2 cm/s as indicated in legend) and BIR8 VCM with vc = 2, with nominal profiles The dashed green line in (A) indicates the labeling efficiency obtained with ideal rect profiles with S0 = κ = 1.0 and β = 0.88. (C,D) Error contributions from bolus groups 1 and 2. All errors are shown as percentages.
(A,B) dmVSS labeling efficiency and CBF estimation error as a function of τ1 for BIR8 LCM and VCM with vl = vc = 2, with nominal BIR8 profile. The dashed green line in (A) indicates the labeling efficiency obtained with ideal rect profiles with S0 = κ = β = 1.0. (C,D) Error contributions from bolus groups 1 and 2. All errors are shown as percentages.
It is also helpful to note that the fractional CBF quantification error can be written as
where α and α′ denote the actual and assumed labeling efficiencies, consistent with what was previously stated for single-module VSASL7.
2.6 Effective Bolus Width and Error
While the effective areas and related metrics, such as quantification error and labeling efficiency, can be readily calculated for arbitrary profiles, it is useful to write each effective area as the product of an amplitude term and an effective bolus width term7. For dmVSI, we have A1,eff = 2S0κ2βτ1,eff and A2,eff = 2S0κβτ2,eff, where
and
are the effective bolus widths. The fractional CBF estimation error for dmVSI can be written as
where the primed variables represent the assumed values associated with 
and 
. To isolate the estimation error due solely to bolus width errors, we may consider the case where the assumed amplitudes are equal to their actual values (
, κ′ = κ, β′ = β), but the assumed bolus widths are equal to their nominal ideal values (i.e. 
and 
, resulting in
with bolus width errors defined as Δτ1,eff = τ1,eff − τ1 and Δτ2,eff = τ2,eff − τ2.
For dmVSS, A2,eff = S0κβτ2,eff where τ2,eff was defined in Eqn. 28. Using the approximation from Eqn. 14, we can write A1,eff ≈ S0κ2β(τ[1,2],eff + (1 − β)τ2,eff) (see Section S.1) where
The fractional CBF estimation error for dmVSS can be written as
Isolating effects due solely to bolus width errors leads to
where Δτ[1,2],eff = τ[1,2],eff − (τ1 + τ2).
2.7 Minimizing Bolus Width Error
While the effective bolus width errors can be calculated for arbitrary functions using Eqns. 27, 28 and 31, it is also useful to gain some theoretical insight into the conditions needed to obtain (i) Δτ1,eff = 0, (ii) Δτ2,eff = 0, and (iii) Δτ[1,2],eff = 0, which are sufficient for minimizing the CBF quantification errors described in Eqns. 30 and 33.
For case (i), if we assume that
holds, as well as additional conditions min(τ1, τ2) > Δtl and sref(vb) ≈ 1.0, then Eqn. 27 can be approximated as:
where sref,E(v) = 1 − pref(v) is an effective saturation function7. Because Eqn. 35 has the form of Eqn. 25 in Ref. [7] and Eqn. 34 is an extended version of condition W1 from that prior work (which we will denote as W1(i)), we may adopt the arguments from the prior work to state that the sufficient conditions to minimize the bolus width error for case (i) are: condition W1(i), an extended version W2(i): PLD ≥ Δtl − τ2 of condition W2 from the prior work, and the additional conditions.
Similarly, for cases (ii) and (iii), we note that the forms of Eqns. 28 and 31 follow that of Eqn. 25 in Ref. [7]. Sufficient conditions to minimize the bolus width errors for these cases are extended versions of condition W1: sref(v) = 0 for |v| ≥ vc where vc ≤ va(τ2) for condition W1(ii) and vc ≤ va(τ1 + τ2) for condition W1(iii); and an extended version W2(ii, iii): PLD ≥ Δtc of condition W2 that applies for both cases (ii) and (iii).
The expressions for the extended versions of Condition W1 are equivalent to: (i) τ1 ≥ Δtl, (ii) τ2 ≥ Δtc, and (iii) τ1 + τ2 ≥ Δtc, where Δtl and Δtc are the LCM and VCM transit delays, respectively, described above in Sections 2.2 and 2.3. Since most dm-VSASL implementations use τ1 > τ2, the active design constraints are typically:
where we have taken into account the additional conditions for case (i). A summary of the conditions is provided in Table 2. With decreases in CBF, the transit delays are expected to increase7 and therefore the values of τ2 and PLD needed to meet the design constraints will also tend to increase.
In addition, we require min(τ1, τ2) > Δtl and sref(vb) ≈ 1.0 for the approximation in Eqn 35 to hold. Since most implementations use τ1 > τ2 and typically Δtl ∼ Δtc, the active design constraints are: τ2 ≥ max(Δtl, Δtc) and PLD ≥ Δtc.
Note that the conditions derived above for cases (ii) and (iii) assume the matched condition7 sref,E(v) = sref(v), which is applicable for dmVSS where the LCM1, LCM2, and VCM typically all use the same saturation module. For dmVSI, the matched condition holds for case (i) since LCM1 and LCM2 use the same inversion module, but typically does not hold for case (ii) due to differences in the LCM and VCM implementations, leading to the mismatch condition7 sref,E(v) ≠ sref(v) for the delivery of the second bolus da,2(t). As discussed in Ref. [7] for single-module VSASL, this mismatch can lead to a CBF-dependent bolus width error that occurs when vl = vc and can be mitigated by setting vl < vc. In Results, we show that the mismatch-related CBF estimation error for dmVSI (Δτ2,eff ≠ 0) can be similarly reduced with vl < vc. Note that Fig. 6 in Ref. [7] provides a useful graphical interpretation of the mismatch error.
(A,B) dmVSI and dmVSS arterial delivery functions (yellow, purple; vl = vc = 2 cm/s, vb = 0.1 cm/s) obtained using average profiles, with delivery functions for nominal profiles shown in blue and red for comparison (originally shown in Figures 2C and 3C). (C,D) dmVSI and dmVSS labeling efficiencies as a function of τ1 for (C) dmVSI with sinc FTVSI LCM (vl = 1.2 or 2 cm/s as indicated in legend) and (D) dmVSS with BIR8 LCM with vl = 2, with BIR8 VCM with vc = 2 cm/s for both dmVSI and dmVSS. Labeling efficiencies with ideal rect profiles and scaling parameters from the average and nominal responses are shown by the green dashed lines and the black dash-dot lines, respectively.
2.8 Blood Volume Components
As noted in Ref. [7], the arterial delivery functions do not account for either (a) potential labeling of arterial blood in the velocity range [vcap, vb], where vcap denotes capillary velocity, or (b) labeling of capillary and venous blood. Under the reasonable assumption that pi(vcap) ≈ 0, the contribution of the capillary blood volume term can be ignored7. Based on the similarity discussed above of Eqns. 28, 31, and 35 to Eqn. 25 in Ref. [7], the arguments from Appendix A.2 of the prior work can be adopted to show that the arterial blood volume terms are compensated for by negative adjustment terms in the bolus width computations, such that the arterial blood volume contributions can also be ignored. As a result, the prior treatment7 of local versus global model timing requirements for single-module VSASL also applies to dm-VSASL.
Adapting the arguments from Appendix A.6 of the prior work7 and assuming τ1 > τ2, minimization of the venous blood volume contribution requires that either (1) τ2 > ΔtV,c or (2) vv(ΔtV,c − τ2, vcap) is not much greater than vcap where ΔtV,c denotes the time for venous blood to accelerate from vcap to vc and vv(t, vcap) is the velocity at time t of a venous spin that starts with an initial velocity of vcap at t = 0. When these conditions are not satisfied (i.e. for small values of τ2), Eqn. A.13 from Ref. [7] can be adapted to estimate the venous contribution.
3 METHODS
Bloch simulations and laminar flow integration were performed using code from the ISMRM 2022 VSASL Bloch Simulation Tutorial8. The BIR8 option with vl = vc = 2 cm/s was used to simulate the VSS module and the sinc FTVSI option was used to simulate the VSI module with vl values of 1.2 and 2 cm/s (with velocity-dependent profiles shown in Figs. 2 and S1 of Ref. [7]). Key parameters were maximum B1 = 20 µT, maximum gradient 50 mT/m, maximum gradient slew rate 150 T/m/s, T1 and T2 set to ∞, and velocity span of 0 to 60 cm/s with a velocity increment of 0.001 cm/s. The profiles computed under these conditions are referred to as nominal profiles. The arterial acceleration model from Ref. [7] was used to calculate va(t) assuming vb = 0.1 cm/s for both normal (Q0 = 12.5 ml/s) and low flow conditions (Q0 = 7.7 ml/s).
Example equivalent arterial delivery functions and cumulative magnetization difference signals were then calculated assuming Tsat = 2.5 s, τ1 = 1.5 s, τ2 = 0.5 s, and T1b = 1.66 s. For comparison, ideal rect VSI and VSS labeling functions with vc = 2 cm/s and κ = β = S0 = 1.0 were also considered. Profiles were the same across LCM1 and LCM2, i.e. l2(v) = l1(v) and c2(v) = c1(v). We then varied τ1 from 1.0 s to 1.9 s (increment of 0.05 s) with τ1+τ2 = 2.0 s, and computed the delivery functions, magnetization difference signals, labeling efficiency, and CBF estimation error for each value of τ1.
Following the approach of Ref. [9], additional Bloch simulations were performed with T1 = 1660 ms, T2 = 150 ms, B1 scaled by a factor of 70% to 130% of its nominal value with an interval of 10%, and ΔB0 (off-resonance) varied from –200 Hz to 200 Hz with an increment of 50 Hz. Each of the resulting profiles was multiplied by a Gaussian weighting factor, computed from the product of two Gaussian distributions centered around mean values of B1 scale = 100% and ΔB0 = 0 Hz with standard deviations of 
and 
Hz, respectively. Average FTVSI and BIR8 profiles were then calculated as the sum of the weighted profiles. Labeling efficiency and CBF estimation error as a function of τ1 were then computed using the average profiles. Using the definitions from Section 2.2, the values of κ and S0 were set equal to the absolute values of the average labeling profiles of the LCM and VCM, respectively, at v = 0. The value of β was then determined by estimating the scaling factor needed to account for residual reductions in the steady-state value of the average passband function not accounted for by κ.
4 RESULTS
4.1 Nominal profiles
Example arterial delivery functions and cumulative magnetization difference signals computed using nominal profiles are shown in Fig. 2(C-F) for dmVSI (LCM: FTVSI with vl = 2 or 1.2 cm/s; VCM: BIR8) and Fig. 3(C,D) for dmVSS (LCM and VCM: BIR8). As compared to the comparable waveforms obtained with ideal rect profiles, the delivery functions using the FTVSI and BIR8 waveforms exhibit smoothly varying leading and trailing edges. In addition, the value of the FTVSI cumulative magnetization difference at TI = τ1 + τ2 + Δtc is lower, reflecting the fact that Pa < 2.0 for the FTVSI passband function, as noted earlier in Section 2.5.
As shown in Fig. 4A, the dmVSI labeling efficiencies for both normal (blue) and low (yellow) flow with vl = 1.2 cm/s are close to the ideal rect profile labeling efficiency (dashed green line) previously described in Section 2.5, and the corresponding CBF errors (blue and yellow lines Fig. 4B) are close to zero in the region of optimal τ1 ∼ [1.5, 1.6]. On the other hand, the labeling efficiencies with vl = 2.0 cm/s are lower than that of the ideal rect profile, resulting in flow-dependent negative CBF errors (red and purple lines). These errors are largely due to the bolus group 2 error (red and purple lines in Fig. 4D), which reflects the mismatch described in Section 2.7. This error is consistent with the smaller area under the curve of da,2(t) for vl = 2 cm/s (Fig. 2C) as compared to that for vl = 1.2 cm/s (Fig. 2E).
For the low flow condition with vl = 2.0 cm/s, the bolus group 1 error (purple line in Fig. 4C) reflects the interaction between the edges of the pref(va(t)) and pref(va(t − τ1) terms in Eqn. 27, which vary more slowly for lower flow. For larger values of τ1, the corresponding values of τ2 = 2.0 − τ1 are too small to satisfy the condition stated in Eqn. 36, giving rise to the changes in the bolus group 2 errors. It is interesting to note that in this regime of small τ2 values, the errors associated with bolus groups 1 and 2 trend in different directions, such that the errors partially cancel out when computing the overall error. The negative trend for bolus group 1 reflects the interaction between the trailing edges of the LCM1 and VCM profiles, while the upward trend for bolus group 2 reflects the fact that for small τ2 there is less time for spins to decelerate such that the relevant portions of the leading edge of the LCM2 passband function occur at lower velocities and exhibit steeper slopes.
For dmVSS, the labeling efficiencies shown in Fig. 5A, for both normal and low flow with vl = 2.0 cm/s, closely follow the ideal rect profile labeling efficiency (dashed green line) described in Section 2.5. As a result, the corresponding CBF errors (Fig. 5B) are relatively small (< 2%) in the region of optimal τ1 ∼ [1.0, 1.2], reflecting the inherent matching of the LCM and VCM functions described in Section 2.7. For larger values of τ1 where the condition stated in Eqn. 36 is violated, the increase in the bolus group 2 error (Fig. 5D), drives the overall increase in error.
4.2 Average profiles
Delivery functions using average profiles (i.e sum of weighted profiles over range of B1 and ΔB0 values as described in Methods) are shown by the yellow and purple curves in panels A and B of Fig. 6 for dmVSI and dmVSS, respectively, with the corresponding delivery functions using nominal profiles (originally depicted in Figs. 2C and 3C) shown by the blue and red curves for comparison. With the average profiles, the amplitudes of the dmVSI delivery functions are considerably smaller than those for the nominal profiles. These are consistent with the predicted amplitudes of 2S0κ2β and 2S0κβ for da,1(t) and da,2(t), respectively, where S0 = 0.92, κ = 0.68, and β = 1.0 for the average profiles as compared to S0 = κ = 1.0, and β = 0.88 for the nominal profiles. Note that for the average profiles the amplitude of da,1(t) is more attenuated than that of da,2(t), reflecting the difference in their respective dependence on κ (i.e. ∼ κ2 vs. ∼ κ).
For the dmVSS average profiles, the reduction in the delivery function amplitudes are consistent with the predicted amplitudes of ∼ S0κ2β (approximation for β close to 1.0) and S0κβ for da,1(t) and da,2(t), respectively, where S0 = 0.92, κ = 0.92, and β = 0.97 for the average profiles as compared to S0 = κ = β = 1.0 for the nominal profiles. Similar to what was observed for dmVSI, the amplitude of da,1(t) shows greater attenuation than that of da,2(t) due to the κ2 dependence.
Labeling efficiencies (both normal and low flow conditions) for the average profiles are shown in panels C and D of Fig. 6 for dmVSI and dmVSS, respectively, with the corresponding labeling efficiencies using (i) the ideal rect profiles (with the S0, κ, and β values of the average profiles) shown for reference and (ii) the ideal rect profiles with the scaling parameters from the nominal profiles (originally shown in Figs. 4A and 5A) shown for comparison. For dmVSI, the vl = 1.2 cm/s curves show good agreement with the ideal rect reference curve, whereas the vl = 2 cm/s show a flow-dependent negative bias, similar to what was observed for the nominal profiles shown in Figs. 4A. Consistent with the markedly lower dmVSI delivery function amplitudes shown in Fig. 6A, the dmVSI labeling efficiencies are greatly attenuated for the average profiles as compared to the nominal profiles.
The dmVSS labeling efficiencies exhibit a moderate degree of attenuation with the average profiles, reflecting the greater robustness of the BIR8 profile to B0 and B1 inhomogeneities. The maximum dmVSI and dmVSS labeling efficiencies assuming average profiles are both approximately 0.6, in marked contrast to the difference in dmVSI and dmVSS maximum labeling efficiencies obtained for the nominal profiles.
5 DISCUSSION AND CONCLUSIONS
We have presented a generalized signal model for dm-VSASL that can integrate arbitrary saturation and inversion profiles. In addition, we have shown how the framework of passband and saturation reference functions originally developed for the analysis of single-module VSASL7 can facilitate the analysis and interpretation of the signals, including the derivation of concise expressions for labeling efficiency and CBF quantification error.
For the case of scaled ideal rect profiles, the signal expressions match those of the prior work2, with the addition of a minor correction term for dmVSS that addresses a small error in the prior work. To facilitate comparison with the prior work, we defined scaling terms κ and β, where κ denotes the attenuation that is common to l(v) and c(v) at v = 0 and β accounts for additional scaling of the passband function that is not accounted for by κ. This definition supports the derivation of labeling efficiency expressions (Eqns. 24 and 25) for scaled ideal rect profiles that can be used to represent the expected labeling efficiency in the absence of bolus width errors (e.g. green dashed lines in Figs. 4A, 5A, and 6(C,D)). It is important to note, however, that our definition of κ and β differs from the prior work2, in which it was assumed that κ accounted only for the effects of T2 decay and β accounted for all other effects. While the product κβ is the same across definitions, the definition presented here is preferred for future work as it is grounded in the generalized signal model and thereby facilitates a rigorous analysis and interpretation of dm-VSASL labeling efficiency and CBF quantification error.
When using realistic Bloch-simulated profiles, we demonstrated that flow-dependent CBF quantification errors can occur for dmVSI due to the mismatch in the profiles of the LCM and VCM functions, similar to those observed for single-module VSI7. This error is primarily due to the mismatch of leading and trailing edges in the second bolus and can be mitigated by selecting a lower labeling cutoff velocity vl < vc, consistent with the findings for single-module VSI7. For small values of τ2, we observed opposing negative and positive trends in the bolus group 1 and 2 errors that partly canceled out in the overall error. Further work to identify approaches for maximizing the cancelation could be of interest.
In this work, we assumed that LCM1 and LCM2 used the same labeling profiles, with matched VCM for dmVSS and mismatched VCM for dmVSI. In the prior modeling of single-module VSASL7, it was noted that there may be cases where where the mismatch between the LCM and VCM profiles could be of interest, for example to characterize the volume of arterial blood in various velocity ranges. For dm-VSASL, there is an even wider space of potential mismatches to consider due to the presence of three modules (LCM1, LCM2, and VCM), and further work to explore this space would be of interest.
Conditions for minimizing bolus width errors were derived under the assumption of matched profiles. For most implementations, where τ1 > τ2, the active requirements are that τ2 is greater than both the LCM and VCM transit delays and that the PLD is greater than the VCM transit delay. As discussed in Ref. [7], these transit delays are likely to exhibit a dependence on voxel size, with preliminary model-based estimates on the order of 100 to 200 ms when vc = vl = 2 cm/s. While further experimental work is needed to refine the estimates, researchers wishing to minimize bolus width errors may want to consider a minimum value of at least 100 to 200 ms for both PLD and τ2 as an interim design rule.
Consistent with prior work we found that the labeling efficiency of dmVSI was higher than that of dmVSS under the assumption of nominal profiles (with T1 = T2 = ∞ and nominal B0 and B1 values). However, when taking into account B0 and B1 inhomogeneities, the labeling efficiencies were found to be comparable, reflecting the the greater robustness of BIR8 dmVSS to the inhomogeneities9. While this is roughly consistent with prior in vivo findings9 of comparably high ASL signal levels with BIR8 dmVSS and sinc FTVSI VSI (which were found to be about 5% lower than sinc FTVSI dmVSI ASL signal levels2), further work is needed to more accurately assess the empirical effect of the inhomogeneities in vivo (see also Section S.2).
As noted in this paper, dmVSS has the advantage of inherent matching of LCM and VCM profiles, resulting in inherently small CBF quantification errors when using vl = vc. On the other hand, the inversion of static tissue in dmVSI (or a dm-VSASL variant that uses VSSI pulses) provides inherent background suppression that can lead to improved temporal SNR even without additional background suppression pulses6. Due to the label/control condition switching in the second module, dmVSI also offers a more balanced distribution of motion-sensitive gradients, which can reduce artifacts from sources such as eddy currents and diffusion attenuation2. The effect of these sources on the dm-VSASL signal can be examined with the generalized signal model through the inclusion of saturation and inversion profiles that reflect the presence of the confounding factors, and further work in this area would be of interest.
Consistent with the prior work we have assumed that the τ1 + τ2 < BDmax, where BDmax denotes the maximum available physical bolus width and depends on the spatial coverage of the RF transmit coil and the mean velocity of arterial spins1,2. This is equivalent to neglecting the contribution of arterial spins that enter the coverage of the RF transmit coil after the first LCM. Extensions of the current model to model this contribution would be of interest and would most likely entail introducing a joint dependence on velocity and space into the passband function that would allow for an estimation of the physical bolus width.
In this paper we have focused on dm-VSASL, which is currently the primary implementation of multi-module VSASL, as prior work indicates that there is unlikely to be an advantage to using more than two modules1. The challenges of using more than two modules include further decreases in labeling efficiency with the addition of each module and potential concerns with specific absorption rate, especially at higher magnetic field strengths1,2. Nevertheless, if VSASL with more than two modules becomes of interest for future work, the approach described in the Appendix can be readily extended to handle additional modules.
In conclusion, by facilitating the use of realistic saturation and inversion profiles, the proposed signal model can provide a more accurate assessment of the performance (e.g. labeling efficiency) of dm-VSASL implementations and improve the quantification of dm-VSASL CBF measures.
6 DATA AVAILABILITY STATEMENT
Analysis code and files to generate the figures and results presented in this paper will be made available upon publication through the Open Science Framework DOI:10.17605/OSF.IO/WKM54.
SUPPORTING INFORMATION
The following supporting information is available as part of the online article:
Figure S1. Example dmVSI and pCASL normalized ASL images ΔMz/M0b in percentage units.
How to cite this article
Liu TT, Chen C, Guo J, Wong EC, and Bolar DS (2024), A Generalized Signal Model for Dual-Module Velocity-Selective Arterial Spin Labeling, Magn. Reson. Med.,.
APPENDIX
Define 
and the longitudinal recovery function for the blood signal as
Note that Λ (Mb, t) is a general-purpose function that returns the longitudinal magnetization given an initial condition Mb and relaxation time t, whereas the longitudinal relaxation weighting functions L1(t) and L2(t) defined in Section 2.1 can be obtained from the expressions provided in the latter part of this Appendix. We use the notation 
and 
to denote the label and control, respectively, magnetization at experimental time t for blood that is delivered at t′ to the boundary velocity vb, where the subscript (i) indicates the ith step in the processing beginning with LCM1 as the first step (i = 1). With the application of the LCM1 at t = 0, we have
Over the interval t ∈ [0, τ1) the longitudinal magnetization is
At the point of the application of LCM2 at t = τ1, we have
for dmVSI and
for dmVSS. The magnetization difference at this point is
for dmVSI and
for dmVSS. After expanding and reorganizing terms, these expressions may be rewritten as
for dmVSI and
for dmVSS, where term 1 reflects the effects of both LCM1 and LCM2 and term 2 reflects the effect of only LCM2.
After the application of LCM2, both the control and label paths go through the same longitudinal relaxation and VCM, so that the magnetization difference for time t > τ1 + τ2 may be written as
To obtain the cumulative magnetization difference at time t, we integrate over all spins that were delivered over the interval t′ ∈ [0, t], resulting in
which can be expanded and rewritten as the equations presented in Section 2.
Footnotes
minor revisions and added supplementary information
