The underlying vascular and metabolic changes are provided by a non-linear inverse solution to our vascular forward model and estimated by curve-fitting the observed hemodynamic responses. In order to illustrate the performance of this inverse model, we demonstrate the characteristics of our model using three examples that examine the local precision, global accuracy, and physiological consistency of our model when used to analyze fMRI, NIRS, or multimodal data.
3.1 Example 1. Parameter sensitivity and local stability of the inverse problem
A requirement of a well-poised inverse model is that the number of independent measurements is equal to or greater than the number of free model parameters. This is informally referred to as the “counting rule” and can provide a qualitative common-sense check for such a model (e.g. (Bamber and van Santen)). To begin our investigation of our windkessel model, we will first discuss this qualitative rule. Our model used for analysis of the fMRI signals has a total of 15 free parameters. Although the ASL and BOLD responses (recovered at 2Hz) contain a greater number of sample points, each sample point does not represent an independent measurement due to the inherent temporal autocorrelation of the evoked response and partial correlation between measurements. In order to apply the counting rule, we must estimate the number of independent degrees-of-freedom to these measurements. The BOLD response is believed to have up to three distinct phases; an initial dip period, a main response, and a post-stimulus undershoot. In order to correctly model the overall BOLD response, for example in a general linear model, three independent canonical functions should be used with one for each of these phases. Since each canonical function requires three degrees-of-freedom for the amplitude, first order timing (e.g. onset or time-to-peak) and second order timing (e.g. temporal width), it follows that a total of nine parameters should be included to model all three phases of the BOLD response. Similar arguments can be applied to determine that oxy- and deoxy-hemoglobin are expected to contain nine degrees-of-freedom each and flow to contain six (with an assumption of no “initial dip” period). By this argument, we can expect that the BOLD / ASL fMRI data should have the minimum number of independent degrees-of-freedom needed for our model. This argument, however, is at best only a qualitative deduction. In order to more rigorously investigate the ability to determine all model parameters from fMRI measurements, we must turn to more formal mathematical definitions of model sensitivity and parameter identifiability.
The stability of our inverse model indicates the feasibility of our model to estimate the model parameters from experimental data and requires that measurements be sensitive to changes in each of the estimated model parameters (reviewed in (Bamber and van Santen)). Model sensitivity can be examined using a variety of numerical methods and has been recently discussed in the context of a similar vascular model by (Deneux and Faugeras 2006
). However, model sensitivity does not guarantee that the estimated values are accurate or unique. Instead, it allows us to test whether a solution to the inverse problem exists. Because of the non-linear nature of our model, we focused on local data perturbation methods to examine the theoretical sensitivity of measurements to changes in the model parameters.
Model sensitivity quantifies the influence of each parameter on the output of the forward model. Two parameters; namely, the mean vascular transit time (τ) and the Windkessel vascular compliance factor (β), are primarily involved in modeling the dynamic relationship between flow and volume. We begin our investigation of model sensitivity by qualitatively examining whether these two parameters affect the temporal dynamics of the hemodynamic response and if it is plausible to estimate these values based on relative temporal dynamics. Based on our model, we predicted that these two parameters would have distinguishable effects on the relative temporal dynamics of the BOLD and ASL measurements. For example, in we see that if vascular transit time is increased (i.e. with constant baseline flow and thus increased baseline volume), the BOLD response amplitude is reduced and the initial dip component is removed. In particular, a longer transit time increases the temporal delay between arterial and venous responses. Likewise, increases in the value of β lead to larger and quicker washout changes that produce larger magnitude changes in BOLD and deoxy-hemoglobin and reduce the relative time-to-peak value of these responses. The compliance parameter also has pronounced effects on the relative magnitude (i.e. ratio) of the oxy-, deoxy-, and total-hemoglobin changes that are measured by NIRS reflecting differing degrees of washout. We note that the relationships between multimodal measurements are uniquely related to these parameters when compared to the individual evoked responses considered independently.
Parameter sensitivity of the evoked hemodynamic response
he sensitivity of our model each parameter in the model can be numerically quantified by examining the Jacobian of the model with respect to the model parameters. We examined the rank and condition of the Jacobian throughout the bounds of the parameter space to determine if the set of measurements was sensitive to each model parameter. We examined these properties for the full multimodal, BOLD and ASL, and NIRS observation models. We found that the Jacobian was full rank for the parameters estimated in each of these models, which implies that an inverse of the model around a localized point is theoretically possible. Singular value decomposition of the Jacobian provided similar conclusions and showed that the number of independent components was equal to the number of model unknowns (for further discussion refer to (Huppert 2007
)). In comparison, when we examined the model with BOLD only measurements, we found that it was not possible to uniquely estimate many of the parameters of the model, including the separation of flow and consumption changes (i.e. the Jacobian was no longer full rank). This finding was consistent with the previous results from Deneux and Faugeras (Deneux and Faugeras 2006
), which attempted similar analysis of BOLD signals alone and found that many but not all of their model parameters could be independently determined. Here, we find that if we use both ASL and BOLD measurements, the model is better defined because information about the model parameters can be additionally determined from the interrelationships of these signals. To further explore the limits of the model, we examined the Jacobian of each model as a function of the simulated sample frequency and determined that the minimum acquisition speed needed was approximately half of the mean transit time. This implies that the estimation of this model will require a minimum sampling rate of above 1 Hz for both flow and BOLD signals, which is unfortunately faster than is possible in most current fMRI studies and was only possible in this study using jittered stimulation presentation.
A further metric of model sensitivity is the variance inflation factor (Γ), which quantifies the minimum uncertainty that can be theoretically expected in the estimate of each model parameter based on an analysis of the slope (Jacobian) of the minimization function around a given local model solution. This analysis allows us to examine how distinguishable the effects of one parameter are from the other parameters. Using the notation discussed in (Deneux and Faugeras 2006
), for a small change in one parameter (θi
), the perturbation of the observation varies by Jidθi
, where Ji
denotes the Jacobian of the forward model with respect to the ith
parameter. If the Jacobian associated with the parameter under consideration is not orthogonal to the Jacobian matrix with respect to the remaining parameters (denoted J
), then the cross talk between parameters results in non-uniqueness of the solution because an error in one parameter can be compensated by errors in other parameters. The maximum precision of each parameter can be estimated for a given level of noise in the measurements from the variance inflation factors as described in (Deneux and Faugeras 2006
) using the following equation.
represents the transpose of the Jacobian matrix. The smaller the corresponding diagonal element of the variance inflation matrix (Γi
), the more precisely the θi
parameter can be determined for a given variance in the observation (Y) and the more sensitive the data set is to changes in the parameter under consideration.
In , we examine the sensitivity of the measurements to changes in model parameters used in each fitting procedure. Note that the full model (ASL, BOLD, and NIRS) has more degrees-of-freedom than fMRI alone, NIRS alone, or NIRS and BOLD models since the value of the baseline blood oxygen saturation is fixed in these models (refer ). The analysis was conducted for a linearization about the parameter set determined from the model fits to the empirical data (). The variance inflation matrix was inspected for other linearization points (not shown) and was found to be in quantitative agreement to these results.
Estimation of Model Parameters
3.2 Example 2. Simulation results and global identifiability
The analysis of the Jacobian matrix described in our previous example indicated that the observation sets were theoretically sensitive to the parameters being estimated and thus the model was expected to be locally invertible. However, the Jacobian can only be defined about a local linearization point in non-linear systems and such sensitivity analysis can be used to probe only the local precision of the model but not its global identifiability. Using this analysis we know that a local solution can be determined but we do not know if this solution is unique.
To investigate the uniqueness of an estimated solution, it is necessary to examine the ability of the model to correctly determine the parameter set from a global search of the parameter space. We ran forward simulations of the model using a sampling of the physiological states and parameters from within the expected physiological ranges (refer to ). To emulate the fMRI and NIRS experimental data, the measurements were simulated at a 2Hz sample frequency with a random additive instrumental noise term (10:1 contrast-to-noise ratio to approximately match the data). We reconstructed system parameters for each simulation using the fMRI data alone, the NIRS data alone, and the full multimodality data set to explore the accuracy of the model to infer the physiological states. For all reconstructions, the Levenberg-Marquardt minimization algorithm was initialized at a random position within the physiological range of the parameters. The same initialization was used for each modality set. Reconstructions were iterated until a set convergence criterion was met (approximately 106 iterations that took around six hours per model estimation). A total of 350 simulations were run for each of the three data sets.
In , we show parametric plots of the simulated (truth) values for the parameters used in the forward simulation against the reconstructed values obtained from each of the three data sets (NIRS only, BOLD and ASL, or NIRS, BOLD and ASL data). The ability of a model to predict one subset of a data set (in this case the NIRS measurements) from analysis of a different subset (the fMRI data) is often used as a strategy to asses if a model has been over-parameterized. We found that relative CMRO2
and arterial resistance changes (ΔRA
) can be estimated accurately with both single and multi- modality measurement sets, consistent with the tests of local sensitivity. We examined the accuracy of the estimation of key static and response-related parameters in the model. We found that the estimates of the Windkessel vascular reserve parameter (β
) had the most variance but could still be estimated fairly accurately in spite of a slight systematic model under-estimation bias at high values of β
with the optical or fMRI only data sets. In , the covariance in the error of the estimates for the parameters is plotted. The lack of high correlation indicated a low interdependence between variables and minimal cross-talk between these parameters in the model. We noted the largest cross-talk (R2
=0.07) between the vascular transit time through the pial-venous compartment and the relative change in CMRO2
. We also observed negative correlation between CMRO2
and the blood oxygen saturation parameters in the BOLD/ASL/NIRS model as also noted in (Huppert et al. 2007
) (not shown). Recall that baseline oxygen saturation was not varied in the NIRS or fMRI alone models since we found these parameters were not identifiable in the reduced data sets.
Examination of global parameter identifiability
3.3. Example 3. Estimating relative CMRO2 changes from empirical measurements
As a final example of the inverse procedure for our vascular model, the model was used to analyze empirical measurements of hemodynamic changes in the human motor cortex following a 2-second finger-tapping task. In order to examine the estimate of the CMRO2 changes and other model parameters given in , the model was first fit using the combined NIRS, ASL, and BOLD data set. The values of baseline blood oxygen saturation from the model fit to the full data set was then fixed and the model was applied to the ASL and BOLD data and finally to the NIRS only data.
The fits of the fMRI only and NIRS only data sets closely predicted the NIRS and fMRI data respectively as shown in . The 75% confidence bounds for the model fit to each measurement and the model predictions of the opposing measurement set are shown as dotted lines in . To verify that the results of the model fit were independent of the initial position of the Levenberg-Marquardt algorithm, we repeated the minimization for different initial seed positions (see ). The model fit using the group-averaged response curves was consistent with the mean of the fits to the five individual subjects (). Inter-subject variations in the parameter estimates were observed, but remained consistent whether the reduced or full data sets were used in the fitting. The values of the estimated parameters and uncertainties are given in . In , we graphically present the estimates for the key model parameters and uncertainties for both the group and individual subject analysis. The estimated dynamic changes in arterial diameter and CMRO2 are shown in .
Dynamic CMRO2 and arterial diameter changes
3.4. Calculated baseline physiological properties
Based on the values of the estimated biomechanical parameters presented in , several additional physiological quantities can be determined and are presented in . Importantly, the vascular transit time and the vascular volume fraction, needed to scale the BOLD signal were both precisely determined. Vascular volume fraction, in units of milliliters of blood per milliliter of tissue, is easily converted to blood volume, in units of milliliters of blood per 100 grams of tissue, given the approximate density of brain tissue (1.04 g/ml (DiResta et al. 1991
)). This baseline blood volume divided by the vascular transit time provides an estimate of baseline blood flow. Likewise, baseline CMRO2
can be calculated from the baseline blood flow and oxygen saturation within the venous compartment. We were able to precisely estimate the venous oxygen saturation only in the model using the combined fMRI and NIRS data and we found that this estimation was not as accurate using only the fMRI only data. Thus, we fixed the value of baseline vascular oxygen saturations in the model estimates using the fMRI only or NIRS only data. Prior studies have noted that the oxygen extraction fraction varies only a few percent between different regions of the brain and this justifies constraining the range of expected cerebral venous oxygen saturation to the range of 60±9% based on previous measurements in the healthy adult brain (He and Yablonskiy 2007
; Raichle et al
2001). The estimated baseline blood volume, blood flow, and oxygen consumption using different combinations of fMRI and NIRS data are shown in . Although we assumed a value for the resting oxygen extraction in the fMRI-only and NIRS-only models, we found that the estimate of baseline blood flow and oxygen consumption based on the estimated scaling parameters and vascular transit time from the fMRI data alone was consistent with the estimate from the NIRS alone fit and also with the combined fMRI and NIRS fit. The values of the estimated baseline physiology for the group data and the mean of the results from the five individual subjects are presented in .
Estimates of baseline physiological parameters