Comparing hemodynamic models with DCM

1In many models of effective connectivity, like structural equation modeling, psychophysiological interactions, or multivariate autoregression, coupling parameters are estimated directly from the measured time-series. However, the causal system architecture one wishes to identify exists at the level of neuronal dynamics. In modeling dependencies among measured data, one implicitly assumes that neural activity can be observed directly. This is tenable for single neuron recordings, but not for non-invasive techniques like functional magnetic resonance imaging (fMRI) or electroencephalography (EEG). For example, in fMRI, the relationship between neural activity and measured blood oxygen level dependent (BOLD) signals is complicated and not fully understood. In short, the interpretation of models of effective connectivity, which ignore the transformation from the neural activity to the BOLD measurements, can be difficult.

2To enable valid inferences about connectivity at the neural level, effective connectivity models have to combine a model of neural dynamics with a biophysically plausible forward model that describes the transformation from neural activity to measured signals. The inversion of such models furnishes estimates of both neural and forward model parameters (see for review). Dynamic causal modeling (DCM) of fMRI data uses a model that conforms to this principle.1 DCM uses a hemodynamic forward model that is based on work by and its extension by. This model consists of three parts (Fig. 1 ). The first describes the link between neural activity and regional cerebral blood flow (rCBF) and is based on linear differential equations modeling a dampened oscillator: changes in neural activity elicit an exponentially decaying vasodilatory signal that is subject to feedback-regulation by the flow it induces. The second part concerns the dependence of the BOLD signal on rCBF-induced changes of blood volume and deoxyhemoglobin content. This so-called “Balloon model” describes the behavior of the post-capillary venous compartment by analogy to an inflated Balloon, postulating a non-linear dependence of BOLD signal on blood volume ν and deoxyhemoglobin content q. Together, these two components represent the hemodynamic state equations which have been the focus of many empirical and theoretical investigations of the BOLD response (e.g., ) and, as a core component of DCM, find widespread use in analyses of effective connectivity (e.g., ). Finally, the output or BOLD signal change equation λ(q,ν) links ν and q to BOLD signal change (Fig. 1). This paper is only concerned with the last part of the overall hemodynamic forward model; i.e., different variants of the BOLD signal change equation. For brevity, we will refer to the BOLD signal change equation as the “BOLD model” (BM) throughout the paper.

3Over the last years, a variety of extensions (and alternatives) have been proposed for both the hemodynamic state equations and the BM (e.g., ). In particular, and suggested a revision of their original formulation of the BM. Compared to the classical BM, assumed a linear dependence of BOLD signal on blood volume v and deoxyhemoglobin content q. Furthermore, they offered revised definitions of the coefficients, which were based on more recent experimental data and depended explicitly on acquisition-specific parameters like echo time (TE).

4Given the importance of hemodynamic modeling for investigations of BOLD responses and analyses of effective connectivity with DCM, the question arises: which formulation is a more appropriate model of empirically measured BOLD responses? In this article, we address this question by embedding different variants of the BM in a well-established DCM of functional interactions in the ventral stream of the visual cortex, which is fitted to BOLD responses that were measured in a group of twelve subjects. We then compare the DCMs with different BMs by means of Bayesian model selection (BMS; ). Specifically, we addressed the following questions:Is the more recent model by a better model than the classical model by for measured BOLD responses?Can the models be improved further by treating ε, the ratio of intra- and extravascular signal, as a free parameter that is estimated from the data?Can a non-linear variant of the revised model be derived by requiring additional consistency in the approximations and does this variant outperform other models?

5Is the more recent model by a better model than the classical model by for measured BOLD responses?

6Can the models be improved further by treating ε, the ratio of intra- and extravascular signal, as a free parameter that is estimated from the data?

7Can a non-linear variant of the revised model be derived by requiring additional consistency in the approximations and does this variant outperform other models?

8Using fMRI data from a group of 12 subjects, we demonstrate that the best model is a non-linear model with a revised form for the coefficients, in which ε is treated as a free parameter. In what follows, we review the hemodynamic model, dynamic causal modeling and Bayesian model selection. We then consider the BM variants that are evaluated in the final section.

9In the treatment of the hemodynamic state variables below, we follow and in using capital letters for specific dimensional values of a state variable and lower case letters for normalized state variables (normalized to their values at rest).

10In their pioneering work, modeled the BOLD signal at rest, S 0, as a sum of intravascular (S I) and extravascular (S E) signal, weighted by resting venous blood volume fraction V 0:

11Based on theoretical and empirical results, derived the following non-linear equation for BOLD signal change ΔS during activation, relative to resting signal S 0 (see lower part of Fig. 1):

12In this BM, ν and q are the venous blood volume and deoxyhemoglobin content (both normalized to their values at rest), and E 0 is the oxygen extraction fraction at rest. The coefficients k 1...k 3 were estimated assuming a field strength of 1.5 T and a TE of 40 ms. At first glance, the BOLD equation in does not seem to possess anything analogous to ε, the ratio of intra- and extravascular signal. However, the derivation of this equation in the appendix to includes a parameter β which is equivalent to ε, but is substituted in their final equation by an approximation based on results by. Without this approximation, one obtains k 3 = 1 − ε. used fixed values E 0 = 0.4 and ε = 0.4 for which k 3 = 2E 0 − 0.2 = 1 − ε.

13The coefficients k 1...k 3 in Eq. (2), as published in and used subsequently in many studies, can be improved in two ways. The first improvement is to drop any assumptions about TE and use the full expressions that Buxton and colleagues derived initially and then simplified. The second improvement involves correcting an error in the derivation of the coefficient k 2. This error was mentioned in the later article by but its nature was not specified. Additionally, we use a more accurate expression for k 1 that does not rely on the assumption that V 0 is very small. Together, these improvements lead to the following expressions for the coefficients: Here, ϑ 0 = 40.3 s−1 is the frequency offset at the outer surface of the magnetized vessel for fully deoxygenated blood at 1.5 T. In this paper, any results for models based on the classical formulation by use this improved form of the coefficients as shown in Eq. 3.2

14The BM described above requires that one knows the evolution of blood volume ν and deoxyhemoglobin content q. The Balloon model by proposed differential equations for ν and q, based on two main assumptions about the biophysical mechanisms involved: (i) small post-capillary vessels react to an increase in inflowing blood like an inflating balloon, and (ii) oxygen extraction is tightly coupled to blood flow. The first assumption means that changes in (normalized) venous blood volume ν correspond to differences in inflow f in and outflow f out with a time constant τ:Here, τ is the mean transit time of blood, i.e., the average time blood takes to traverse the venous compartment, and corresponds to the ratio of resting blood volume V 0 to resting blood flow F 0: τ = V 0/F 0. Following the results by 3, outflow is modeled as a function of volume with a single parameter α that represents the resistance of the venous balloon, i.e., vessel stiffness:

15The second assumption above (i.e., that oxygen extraction is tightly coupled to blood flow) determines the equation for deoxyhemoglobin content q. Generally, the change in q corresponds to the delivery of deoxyhemoglobin into the venous compartment, minus that expelled. Assuming fully oxygenated hemoglobin in pre-capillary blood, the delivery of deoxyhemoglobin into the venous compartment corresponds to the product of blood inflow and the oxygen extraction fraction E, whereas its clearance is the product of outflow and deoxyhemoglobin concentration q/ν:

16Assuming that any extracted oxygen is metabolized immediately, thereby maintaining a tissue oxygen concentration of zero, the oxygen gradient across the capillary wall and oxygen extraction rate depends entirely on oxygen delivery, and thus on blood flow. showed that a reasonable approximation across a wide range of conditions is

17Together, these considerations led to propose the following state equations for v and q (from now on, blood inflow will simply be denoted as f):

18Eq. (8) shows that the only state variable that ν and q, and thus the BOLD signal, depend on is blood inflow f. The issue of how blood flow depends on neural activity was addressed by. In their completion of the Buxton model, vascular responses to neural activity correspond to a dampened oscillator: changes in neural activity x elicit an exponentially decaying vasodilatory signal s that is also subject to feedback regulation by the flow f it induces (Fig. 1): Here, κ and γ are the rate constants of signal decay and feedback regulation, respectively. Note that f is normalized flow with regard to resting flow F 0, and therefore the feedback regulation term (f − 1) in Eq. (9) becomes zero during rest.

19This linear model of the relation between neural activity and rCBF concurs with several experimental results. In particular, the results of combined perfusion and BOLD measurements by during various stimulation conditions were “consistent with a model with a non-linear step from stimulus to neural activity, a linear step from neural activity to CBF change, and a non-linear step from CBF change to BOLD signal change.” This is exactly what the hemodynamic model represents. Note that neural activity x, driving the vasodilatory signal in Eq. (9), is the output from the neural state equation in DCM (see Fig. 1 and below). This provides for a flexible model of the link between neural activity and rCBF that can capture a large variety of transients, sustained responses, and adaptation effects—in brief, any kind of neural dynamics that can be modeled with bilinear differential equations as used in DCM (see Fig. 1 in and Fig. 5 in for examples).

20Buxton and colleagues recently proposed a revised form of the classical BM (Eq. (2)). The critical changes they proposed were (i) a linear dependence of BOLD signal on blood volume v and deoxyhemoglobin content q and (ii) new coefficients k 1...k 3 which were computed on the basis of more recent experimental data and which explicitly consider acquisition-specific TE and the ratio of intra- and extravascular signals (ε): According to, for a field strength of 1.5 T, r 0 = 25 s− 1 is the slope of the relation between the intravascular relaxation rate R 2I* and oxygen saturation. Furthermore, they chose fixed values for the resting oxygen extraction fraction (E 0 = 0.4) and for the ratio of intra- and extravascular signal (ε = 1.43).

21Both the classical BM by and the more recent version by can be regarded as special cases of the following generalized BM:

22This model is derived in Appendix B, starting from the same point as, but using more consistent approximations. It retains the non-linear output equation in the classical Buxton model and combines it with the revised coefficients from the Obata model. It is more general than these two models: its coefficients can be reduced to those of the classical model, and its derivation can be simplified to yield the Obata model (see Appendix B). This model is a useful starting point to consider a number of BM variants that differ in terms of (i) classical (Buxton) vs. revised (Obata) forms for the coefficients, (ii) linear vs. non-linear output equations, and (iii) fixed vs. free parameters, ε, for the ratio of intra- and extravascular signals. These variants, all of which we fitted to empirical data and compared using BMS, are described in the next three sections. The nomenclature of these models is summarized in Table 1.

23The classical model obtains by simplifying the revised expressions for the coefficients in Eq. (12) through reduction: it assumes a fixed TE of 40 ms and omits the dependence of k 2 on ε, the ratio of intra- and extravascular signals. In the following, all models based on BMs with classical coefficients will be denoted as “CBM” whereas all models based on BMs with revised coefficients will be referred to as “RBM” (compare Table 1). In short, the distinction between CBM and RBM rests on the substitution:

24Over and above the non-linearities in the state equations for blood volume and deoxyhemoglobin content (Eq. (8)), the Buxton BM and our generalized BM have a non-linear form (Eqs. (2) and (12)) whereas the Obata BM is linear (Eq. (11)). The Obata BM can be regarded as a linearized version of the generalized BM which is due to a simplified approximation to the exact signal change equation (see Appendix B). Concerning the Buxton BM, the non-linearity is introduced by the intravascular term: φ(q,ν) = 1 − q/ν. In order to obtain a linearized version of it, it therefore suffices to complement the reduction of the coefficients described above with a linearization of this term around the resting state (q = ν = 1), using the first-order Taylor approximation; φ(q,ν) ≈ ν − q:

25Throughout this paper, linear and non-linear models are distinguished by the subscripts “L” and “N”, respectively. For example, the Obata model is referred to as “RBML”, whereas the original Buxton model is referred to by “CBMN” (see Table 1).

26A critical parameter in BM is ε, the ratio of intra- and extravascular signals. This quantity is not well characterized because it is difficult to disambiguate intra- and extravascular BOLD signal experimentally. estimated ε = 1.43 by assuming intra- and extravascular relaxation time constants of T 2I* = 90 ms (corresponding to R 2I* = 11.1 s− 1 ) and T 2E* = 50 ms, respectively, and equal spin-densities in intra- and extravascular compartments. The exact basis of these assumptions is not quite clear. For the T 2I* value, referred to the study by who measured T 2* values in large vessels of pigs. However, in their paper, state: “... the R 2* values of arterial (Y ≈ 93%) and venous (Y ≈ 75%) blood measured in our pig studies are approximately 4 and 6 s− 1, respectively, whereas those measured in human volunteers in our previous work are approximately 5 and 10 s− 1, respectively...” For their T 2E* estimate, did not give a reference.

27Overall, there is considerable uncertainty about the value of ε in the literature. In contrast to, assumed that ε = 0.4. Recent data from, obtained with a “vascular space occupancy” technique that nulls intravascular contributions to the BOLD signal, suggest the contribution of extravascular signal to the total change in R 2* during activation is 47% at 1.5 T. This implies that in their measurements ε ≈ 1. Furthermore, the results of imply that ε depends on complex interactions between TE and field strength.

28This suggests that instead of fixing the value of ε, it may be more appropriate to acknowledge the uncertainty about it and make it a free parameter of the BM. Additionally, because different brain regions can differ in their vascularization and therefore in the intravascular contribution to the measured BOLD signal, different values of ε for each brain area in the DCM may be appropriate. The question remains however, how one should choose the prior density of ε. Since ε is the ratio of intra- and extravascular BOLD signal at rest, it is positive. This is assured with a log-normal prior density (see Fig. 2 ). Furthermore, given our uncertainty, this prior should be fairly flat. In our models, we used ε = μ ε exp(ϑ ε), where the scale-parameter ϑ ε had a prior mean of zero and a variance of 0.5. This allows for approximately an order of magnitude variation about μ ε; a range that is substantially larger than the estimates of ε described above. Furthermore, we chose μ ε = 1 (i.e., equivalent intra- and extravascular signal contributions at rest). This is a value roughly in the middle of the ε estimates at 1.5 T reported so far (see above).

29The distinction between treating ε as fixed or free parameter induces the final variant of the BMs we considered in this study. In this paper, a model in which ε is treated as a free parameter for each area is indicated by the suffix “(ε)”, whereas a model in which ε is a fixed parameter does not have a suffix (see Table 1).

30DCM for fMRI uses a simple model of neural dynamics in a system of n interacting brain regions where neural population activity of each region is represented by a single state variable. It models the change of this neural state vector x in time as a bilinear differential equation:Here, the A matrix represents the fixed (context-independent or endogenous) strength of connections between the modeled regions, and the matrices B (j) represent the modulation of these connections (e.g., due to learning, attention, etc.) induced by the jth input u j as an additive change. Finally, the C matrix represents the influence of direct (exogenous) inputs to the system (e.g., sensory stimuli). Note that all parameters are rate constants and are thus in units of s− 1.

31To explain regional BOLD responses, DCM for fMRI combines this model of neural dynamics with the hemodynamic model described above. Together, the neural and hemodynamic state equations yield a deterministic forward model with hidden states. For any given combination of parameters θ and inputs u, the measured BOLD response y is modeled as the predicted BOLD signal h(u,θ) plus a linear mixture of confounds Xβ (e.g., signal drift) and observation error e:

32DCM uses a fully Bayesian approach to parameter estimation, with empirical priors for the hemodynamic parameters and conservative shrinkage priors for the coupling parameters (see and for details). Briefly, the posterior moments are updated iteratively using variational Bayes under a fixed-form Laplace (i.e., Gaussian) approximation, q(θ), to the conditional density p(θ|y). This includes a gradient ascent on a variational free-energy bound on the marginal likelihood to optimize the maximum a posteriori (MAP) estimate of the parameters in the E-step of an EM algorithm, whereas the M-step is concerned with optimizing hyperparameters λ that control the covariance components of observation error.

33Given some observed data, which of several alternative models is optimal? The decision cannot be made solely by comparing the relative fit of competing models. One also needs to account for differences in complexity; i.e., the number of free parameters and the functional form of the generative model. This is important because as model complexity increases, fit increases monotonically, but at some point the model will start fitting noise that is specific to the particular data (i.e., “over-fitting”) and thus becomes less generalizable across multiple realizations of the same underlying generative process. Therefore, the question “What is the optimal model?” can be reformulated as “What is the model that represents the best balance between fit and complexity?” This is the model that maximizes the model evidence:4

34Here, the numbers of free parameters (as well as the functional form) are subsumed by the integration. Unfortunately, this integral cannot usually be solved analytically; therefore an approximation to the model evidence is used. This is usually the free-energy bound on the log-evidence. In DCM, the negative free energy F is the objective function for inversion:Here, KL denotes the Kullback–Leibler divergence between the approximating posterior density q(θ) and the true posterior p(θ|y,m). After convergence of the estimation scheme, the KL term is minimized and F ≈ ln p(y|m). Rewriting Eq. (17) asadds a useful perspective, The two terms in Eq. (18) can be understood as encoding the two opposing requirements of a good model: that it explains the data and is as simple as possible (i.e., uses minimal number of parameters that deviate minimally from their priors). Eqs. (17) and (18) are derived in Appendix A.

35To quantify the relative goodness of two models m i and m j, the differences in their log-evidences can be transformed into a Bayes factor (BF):

36The group Bayes factor (GBF) for any given model m i, relative to m j, iswhere k indexes subjects. As detailed in, the GBF is equivalent to the product of the subject-specific Bayes factors for a given model comparison. It rests on the assumption that model evidences are independent across subjects; this is tenable if the subjects are statistically independent samples from the population. Furthermore, given a flat prior on the model, the product of model evidences is equivalent to multiplying the posterior probabilities of all models.

37Just as conventions have developed for using p-values in frequentist statistics, there are conventions for Bayes factors. For example, suggests an interpretation of Bayes factors as providing weak (BF < 3), positive (3 ≤ BF < 20), strong (20 ≤ BF < 150), or very strong (BF ≥ 150) evidence of one model over another. A complementary index to the GBF is the positive evidence ratio (PER; ), i.e., the number of subjects where there is positive (or stronger) evidence for model m i divided by the number of subjects with positive (or stronger) evidence for model m j where k = 1,…,N and ∣·∣ denotes set size. The GBF is sensitive to outliers, whereas the PER is not. In contrast, the PER is insensitive to the magnitude of the differences across subjects while the GBF is not. Therefore, the GBF and the PER can be used as complementary measures: the GBF gives a quantitative account of the difference, pooling over all data (i.e., subjects) but ignoring inter-subject variability, whereas the PER describes the qualitative reproducibility of model comparison over subjects.

38We compared the variants of hemodynamic models described above by placing them in a DCM and comparing the resulting models using BMS. We used a well-characterized DCM of interacting visual areas as a vehicle for this comparison. This DCM is the four-area ventral stream model described by, which had emerged as the optimal model from a systematic comparison of sixteen alternatives. Fig. 3A shows the structure of this model. The fMRI data are from the study by and were obtained at 1.5 T and a TE of 66 ms. The advantage of using a full DCM (instead of modeling a single area as in ) is twofold: first, the resulting inference does not depend on the choice of a particular region, and second, the neural state equation in DCM can model a wide range of neuronal transients, sustained responses, and adaptation effects and thus affords much higher realism than other hemodynamic models (see Discussion). Furthermore, our comparisons were performed separately for each of 12 subjects, to avoid conclusions that were specific to a particular individual.

39Our model comparison adopted a factorial approach: we compared (i) classical vs. revised forms for the coefficients, (ii) linear vs. non-linear BMs, and (iii) fixed vs. free ε. All combinations were tested so that we could establish the relative importance of the three model attributes. This resulted in eight DCMs per subject. As an additional reference model for displaying the log-evidences, we chose the model of in its original implementation, where both ε and resting oxygen extraction fraction, E 0 have a fixed value. In our experience, treating E 0 as a free parameter typically improves the evidence of hemodynamic models. Therefore, all of our eight models treated E 0 as a free parameter. Indeed, model comparison showed that all eight models were superior to the reference model (see Fig. 4 ).

40The same DCM, equipped with different BOLD output equations, was fitted to the empirical fMRI data from each subject. Fig. 4 summarizes the results of all model comparisons. For each model, it shows the log of the GBF; i.e., the sum of its log-evidences across the group minus the summed log-evidence for the reference model (cf. Eq. (20)). It can be seen that the best of all models is the non-linear model with revised coefficients and free ε, i.e., RBMN(ε). In the following, we describe the individual comparisons in more detail. We first report the direct comparison between the classical Buxton model (CBMN) and the Obata model (RBML). These two models differ in terms of both reduction (of the coefficients) and linearization; therefore we then investigate whether any differences depend on treating ε as a free or a fixed parameter. Finally, we examine the role of non-linear vs. linear BOLD models.

41First, we investigated the effect of freeing E0 in the Obata model. Our implementation, which treated E 0 as a free parameter (i.e., RBML), was slightly superior to the original formulation by in which E 0 had a fixed value (GBF = 3.49 × 102; PER = 1:0; second column of Table 2 ). Note that the latter is used as a reference model for plotting the results in Fig. 4. In a next step, we compared the linear RBM with the classical Buxton model (CBMN). This comparison indicated that CBMN was superior in each and every subject (PER = 12:0), giving a GBF of 8.41 × 1025 (Table 2, third column). We next investigated the two possible reasons why RBML performed worse than CBMN: (i) a suboptimal value of ε, and/or (ii) the failure to model non-linearities in the BOLD equation.

42For the linear RBML, freeing ε dramatically improved the model evidence: RBML was clearly inferior to RBML(ε) (GBF = 2.12 × 10− 30; PER = 0:12; Table 2, fourth column). RBML(ε) also performed better, albeit only slightly, than CBMN (GBF = 5.60 × 103, PER = 2:1).

43In contrast, for CBMN, freeing ε did not improve it and actually slightly decreased its log-evidence: CBMN(ε) performed worse than the original CBMN with fixed ε (GBF = 5.24 × 10− 2, PER = 1:3). These comparisons illustrate that å plays a different role within the linear RBML(ε) and the non-linear CBMN(ε), respectively: in the latter, including ε as an additional free parameter per region led to an increase in model complexity (by increasing the overall number of model parameters and changing how other parameters had to deviate from their priors) that outweighed the improvement in data fit, whereas in the former it did not. Would the same results also hold for the linear CBML(ε) and the non-linear RBMN(ε)?

44Fitting the linear CBML to the data and evaluating its log-evidence showed that it performed worse than the other CBM variants described above (see Fig. 4). With fixed ε, although it was still clearly superior to the linear RBML (GBF = 1.57 × 1021, PER = 12:0), it was not as good as the original CBMN (GBF = 1.87 × 10− 5, PER = 0:3), and the linear RBML(ε) (GBF = 3.34 × 10− 9, PER = 0:3). Freeing ε further decreased its log-evidence: the linear CBML(ε) turned out to be the worst of all CBM variants tested here (see Fig. 4).

45When the non-linear RBMN was evaluated using BMS, using a fixed ε, it was superior to the linear RBML with fixed ε (GBF = 3.63 × 103; PER = 3:0), but inferior to the CBMN (GBF = 4.32 × 10−23; PER = 0:11) (see Table 3 ). However, when ε was treated as a free parameter in RBMN, it became superior to all other models (see Fig. 4). It even performed slightly better than the previously best model, the linear RBML(ε) (GBF = 1.92 × 102). However, examination of the subject-specific BFs in Table 3 (rightmost column) shows that this improvement is fairly subtle: the individual BFs are only in the range 0.97 to 2.32. This indicates that the differences are very small at the level of individual subjects. However, given the consistent direction of the individual BFs over the group, the GBF reflects clear evidence in favor of the non-linear RBMN(ε). In summary, freeing ε decreased the goodness of both linear and non-linear CBMs, but improved both linear and non-linear RBMs, with the non-linear RBMN(ε) emerging as the best model from our comparisons.

46An interesting issue is whether estimates of the new parameter ε allow for straightforward interpretation, or whether this parameter exhibits conditional dependencies with other parameters. The region-specific parameter estimates for ε, averaged across subjects, were 0.63 ± 0.43 (left lingual gyrus, LG), 0.58 ± 0.36 (right LG), 1.98 ± 0.45 (left fusiform gyrus, FG), and 1.38 ± 0.46 (right FG). It is striking that for the regions receiving direct inputs (left and right LG; compare Fig. 3A) the estimates are lower than the prior mean of ε. In contrast, for the regions not receiving direct inputs (left and right FG), they are higher.

47As demonstrated in Fig. 5, increasing ε increases the amplitude of the evoked BOLD response (although note that the effect of changing å depends on the value of E 0). This effect is similar to that achieved by increasing the driving inputs (C in Eq. (14)). In order to fit a given BOLD response, any change in ε or C, respectively, would therefore have to be compensated, to some degree, by changing the other parameter in the opposite direction. This means that the two parameters are likely to show inverse conditional correlations. This was confirmed by an analysis of the estimated posterior covariances of the parameters in the RBMN(ε) model. Since covariance matrices are difficult to interpret visually, we normalized the posterior covariance matrix from each subject to provide a posterior correlation matrix. Fig. 6 shows the average posterior correlation matrix over subjects. The rectangles in this figure highlight three results: first, ε shows a fairly strong negative correlation with the input parameters to the system, C (see r ε,C in Fig. 6). As expected, this posterior correlation was particularly pronounced (up to − 0.61) for the ε’s in the regions that received driving inputs (i.e., left and right lingual gyrus, LG). Second, ε also shows a negative correlation (up to − 0.53) with the fixed strengths (A in Eq. (14)) of those connections that convey activity elicited in input areas (left and right LG) to those areas not receiving direct inputs (left and right FG). These connections play a similar role in determining the amplitude of the BOLD response in left and right FG as do the driving inputs for left and right LG. Most importantly, however, ε is not strongly correlated with the main parameters of interest within DCM, that is, with the parameters encoding the context-dependent modulation of connection strengths (B in Eq. (14)). These correlations did not go beyond − 0.29 and were mostly close to zero (see r ε,B in Fig. 6).

48Tables 4 and 5 show how a change in the hemodynamic model affects the average neuronal parameter estimates and their standard error across subjects. It can be seen that in models with free ε the conditional dependencies between ε and the neuronal parameters described above have only modest effects on the parameter estimates. As expected, the changes are most profound for the C parameters (Table 5), which are usually not of interest. Concerning the A/B parameters (Table 4), in no case were differences induced by freeing ε large enough to affect second-level statistical inference about the parameter estimates5. Finally, as can be seen in the tables, other changes in the hemodynamic model affect the parameter estimates even less than freeing ε.

49In this study, we implemented a variety of different BOLD signal models within the hemodynamic forward model of DCM. We chose an established four-area DCM of the visual cortex which was fitted to empirical fMRI data from 12 subjects. The relative goodness of the different BOLD signal models was assessed by means of BMS. This comparison gave three main results. First, we showed that, in its original formulation, the linear RBML (which was recently introduced by ) is not superior to the CBMN (the classical BM by ). Second, we have demonstrated that RBML can be made to outperform the CBMN if ε, the ratio between intra- and extravascular signal components, is allowed to vary as a free parameter. Third, we have derived a generalized non-linear BM with revised coefficients that outperforms any other model tested, when ε is treated as a free parameter (the RBMN(ε) variant in this study). This model was founded on the previous observation that a non-linear output equation improves modeling of the BOLD signal. Because of its superior performance in this study, this hemodynamic model will replace the CBM in the next update of the DCM routines in the open source software package SPM5 (http://www.fil.ion.ucl.ac.uk/spm). The other hemodynamic models tested here will be included as options, allowing users to compare alternatives and find the most appropriate hemodynamic model for their particular experimental set-up.

50The posterior covariances among parameters suggest that most parameters of interest in a DCM, i.e., the connectivity among regions and particularly its context-dependent modulation (A and B matrices in Eq. (14)), are not strongly correlated with the hemodynamic parameters, including the new parameter ε (Fig. 6). This relative independence of the neuronal and hemodynamic parameters is an important feature of DCM, which largely results from normalizing the coupling parameters with regard to the system’s dampening rate-constant (i.e., the diagonal of the A matrix; see for details). This renders them less dependent on the amplitude of the hemodynamic response and the neuronal input parameters (C in Eq. (14)). If the two sets of parameters were strongly coupled, it could be difficult to disambiguate the influence of neuronal connectivity from neurovascular coupling mechanisms on the resulting BOLD response. Tables 4 and 5 lists all mean parameter estimates in our model and their standard error across subjects for each of the hemodynamic models tested in this study. It can be seen that in no case changes in coupling parameters were large enough that the nature of statistical inference (as obtained by a second-level t-test) would have changed.

51From a broader perspective, these results suggest that the particular form of the hemodynamic model is not terribly important for inference about the underlying neuronal dynamics and how they were caused. Although one might have intuited this from the rather simple form of empirically derived hemodynamic impulse response functions, it is reassuring to see that the specific parameterizations employed by the biophysical models in DCM provide inference on neuronal parameters that are robust to changes or misspecification of the hemodynamic model. Although the nature of the hemodynamic model may not affect inferences at the neuronal level, it is clearly important for inference on the hemodynamic parameters that may be of interest when studying regional variations in physiology or pathophysiology per se.

52Our characterization of conditional dependencies among the parameters is related to measures of system identifiability. Identifiability can be addressed from various angles. For example, sensitivity analyses investigate the partial derivatives of model output with regard to the parameters. A system is not identifiable if, for any change of a given parameter, an identical effect on model output can be achieved by changing one or several other parameters. In Bayesian approaches to system identification, conventional sensitivity analyses and analyses of posterior covariances have a close relationship. An example of this was given by who performed a sensitivity analysis for the hemodynamic model of, deriving a sensitivity factor that was proportional to the inverse of the posterior covariance of the parameters (under the Laplace approximation this inverse is also known as Fishers Information matrix). We used the conditional covariances directly (after normalizing them to form conditional correlation matrices). These matrices encode the degree of interdependence between the parameters (see Fig. 6).

53Several considerations should be kept in mind when interpreting the results of this study. First, the results were obtained using a particular data set and BOLD time-series from visual cortex. Even though we used four different visual areas and replicated the results in twelve subjects, we cannot ensure that the optimal model is invariant across brain regions and data sets.

54Second, the various hemodynamic models in this study were investigated as an integral part of DCM, a causal model that links experimentally designed manipulations (e.g., presentation of sensory stimuli) via predicted neural and vascular dynamics to observed BOLD responses of discrete brain areas. Although this model is fairly abstract at the neural level, representing neuronal population activity by a single state variable for each region, it allows for a much more flexible representation of neural responses to external perturbations than other hemodynamic models available. For example, it is not constrained to adaptation effects of a particular form as in, but can represent any kind of neural dynamics that can be modeled with bilinear differential equations. This comprises a wide range of transients, sustained responses and adaptation effects (see Fig. 1 in and Fig. 5 in for examples).

55Third, the mathematical form of the dependence of the BOLD signal on deoxyhemoglobin content q and blood volume v remains an interesting research question. Several previous studies suggested that the CBF–BOLD relation has significant non-linear components, and the work presented here supports this view. These non-linearities could enter the system at various levels of the biophysical process. In fact, the state equations for both q and ν are non-linear (although, given the typical range of values for q and ν, these non-linearities are rather weak). Thus, even a model with a linear BOLD equation as, for example, the model by can show a non-linear CBF–BOLD relationship. However, in agreement with, the present results suggest that the output non-linearity of the Balloon model is important, over and above the non-linearity of the state equations for q and ν. The reason for this is that the estimation scheme in DCM uses a bilinear approximation to the states while retaining the output non-linearity (cf. ).

56Fourth, the hemodynamic models in the present study assumed a field strength of 1.5 T. The finding that the neuronal connectivity parameters are fairly robust to changes in the hemodynamic output equation (Table 4), suggests that the current implementation of DCM can also be used in situations for which the coefficients in the hemodynamic model are not optimized (e.g., field strengths higher than 1.5 T). It should also be possible in the future, however, to adapt the model to higher field strengths since most constants in Eq. (12) are known or can be computed for higher field strengths. The main problem is ε for which, to our knowledge, there is a lack of reliable empirical estimates, particularly at high field strengths. However, the present study demonstrated that one can treat ε as a free parameter, and using a flat prior for this parameter properly accommodates the uncertainty about its value.

57Finally, the importance of non-linearity in the BOLD output equation depends on the exact stimulation conditions; in particular it is likely to increase with shorter stimulus-onset asynchronies (SOA) than the ones used in this study (1.5–2.5 s) and to decrease with longer SOAs. In any case, the procedure we introduced here (and the code that will be made available in SPM5) enables the user to examine this question for her or his specific data, by comparing different hemodynamic models.