Effective connectivity refers to the causality relationship between different brain regions, or one brain region influences another distant brain region. Generally, there are three approaches for modeling nonlinear brain activity in the study of functional magnetic resonance imaging (fMRI) effective connectivity. These methods can be categorized as white-box, gray-box, and black-box model methods [1]. White-box modeling is an off-building approach where the model is built before quantifying the connectivities. White-box modeling assumes that the structure of the neuron population is known, to address the problem of finding the parameters of the assumed structure. It is difficult to apply this identification method because the assumption of a known, stationary neuron population does not hold. Alternatively, gray-box models are often adopted when limited structural knowledge is available or assumed, and a model structure can be built on physical grounds, with a certain number of parameters to be estimated from the data. This could be a state-space model of a given order and structure. For example, using dynamic causal modeling (DCM) [2–5], an offline state-space model can be built at neuron level, and the effective connectivity can be studied by identifying the known structure of models corresponding to an infinitely precise prior on the model. A major limitation of this method is that it uses one offline model for all the brain regions without considering the variability of the fMRI response from different regions. For the black-box modeling method, the physical insight is not available or assumed, but the chosen model structure belongs to families that are known to have good flexibility and therefore can be adapted more easily to the data. Current black-box model approaches to study effective connectivity include the Granger causality model (GCM) [6–10], the (sparse) multivariate autoregression (MAR) model [11–13], and autoregression and moving average (ARMA) models [14]. These methods, however, do not include the experimental design as input for modeling the fMRI effective connectivity. This is not optimal when knowledge of the experimental design is available [2]. To circumvent this limitation, we have introduced a new nonlinear system identification method (black-box model method) for studying the effective connectivity [15]. We devised a scheme which identifies nonlinear connectivities using a nonlinear autoregressive exogenous model (NARX) method and provides statistics which can be used to test model interactions.

In this chapter, we firstly introduce the basic theory of NARX for fMRI effective connectivity analysis. Second, we show how the model is constructed for the investigation of nonlinear dynamic connectivity. In addition, we present the least squares algorithm to identify the strength of the various connectivities. Thirdly, we give an example of applying this method to the human visual system to study the brain region interactions. Fourthly, we extend this method for resting-state fMRI (rfMRI) data analysis. Finally, we introduce a robust statistical method for studying fMRI effective connectivity analysis at second-level data analysis.

The nonlinear brain black box is depicted in Fig. 4.1a [5, 16, 17]. The input of the black box is the experimental design u(t) (Fig. 4.1b), and the output y(t) (Fig. 4.1c) is the fMRI response with some random noise e(t) added. Figure 4.1b is one example of a brain system input from phase-encoded and standard block designs. The system output or fMRI response y(t) is adaptively changed according to the input u(t) as shown in Fig. 4.1c. Formally, the physiological processes underlying the BOLD response can be modeled as a multiple-input and multiple-output (MIMO) system [16]:

and its discrete form is
where f and g are nonlinear functions and θ represents the set of model parameters. Under some mild assumptions, the discrete-time multivariate system (4.1) with p outputs (e.g., response from different regions) and q inputs (e.g., different types of stimuli) can be described by an autoregressive moving average with exogenous input (NARMAX) as follows [18, 19]:
where , , are the system output, input, and noise, respectively; n _y , n _u , and n _e are the maximum lags in the output, input, and noise; e(t) is a zero mean independent sequence; and f _g is a new nonlinear function which can be obtained from nonlinear functions f and g. A special case of the general NARMAX model (4.3) is the nonlinear autoregressive with exogenous inputs (NARX) model which can be expressed as

By applying the regression equation, the NARMAX model (4.3) and NARX model (4.4) can be approximated as [20–22]
where P ₀(t) = 1; for M ≥ 1, P _m (t) = y ₁ … y _i u ₁ u ₂ … u _j , i ≥ 1, j ≥ 0; m is the number of nonlinear terms; M is the system order; N is the total number of time point in the time series; i is the number of connected regions; and j is the number of inputs. Equation (4.5) denotes a general case where both input and output terms may be present, but it should be understood that some of the P _m may contain only input or output terms and cross products. For example, for two stationary series of N values, the inputs and y ₂ and output y ₁ of a closed-loop time-invariant nonlinear brain system can be described as [8]
where the coefficients c ₀, {a ₁(i); b ₁(j); c ₁}, and {a ₂(i, j); b ₂(i, j); c ₂(i, j)} denote constant (zeroth order), linear (first order), and nonlinear (second order) contributions to y ₁(t), respectively. represents the experimental input and is the prediction error of y ₁(t). The model orders S ₁ and S ₂ are the maximum lags of the linear and nonlinear autoregressive (AR) influences, respectively, while the maximum lags for linear and nonlinear exogenous effects are determined by the model orders T ₁ and T ₂. The model can be represented in the matrix form as
where the vector Y = [y ₁(1), y ₂(2), … y ₁(N)]^T contains values of output series; is the prediction error series; is the experimental input time series; A ₁, B ₁ and C ₁ are the first-order vector coefficients; A ₂, B ₂, and C ₂ are the second-order vector coefficients; and the matrices and contain the S ₁ linear AR terms and the (T ₁ + 1) linear exogenous terms, respectively:
the matrix contains the S ₂(S ₂ + 1)/2 quadratic AR terms given by the product of the terms of the matrix . In the same way, the matrix contains the (T ₂ + 1)(T ₂ + 2)/2 quadratic exogenous terms and the matrix contains the S ₂(T ₂ + 1) cross terms. Equation (4.7) can be written in the compact matrix format as
where , β=[c ₀, A ₁^T, B ₁^T, A ₂^T, B ₂^T, C ₂^T, c ₁]^T. Coefficient matrix β can be estimated by least squares:
where pinv is the Moore–Penrose pseudoinverse of the matrix. By neglecting the nonlinear terms and experimental input and considering only the first order of AR, i.e., AR(1), these lead to
or

This is the well-known two-connection linear GCM in fMRI data analysis.

Once the coefficients of the model are determined (e.g., Eq. (4.9)), Granger causality tests [6, 23, 24] are derived based on F/T statistics. For simplicity and illustrative purposes, we take the two-connection nonlinear models (4.6) as an example. The same principle can be applied for the linear system (Eqs. (4.11) and (4.12)). The tests for determining Granger causes (GC) are [23]:

T and F statistics are developed to detect significant relations. From Eq. (4.8), and partitioning the coefficients as β = (β ₁ : β ₂) and W = (W ₁ : W ₂) accordingly, we can write this test as
with the maintained hypothesis. Therefore, T statistics given in Chap.​ 2 can be applied for testing the hypothesis. For the F-test, we employ the following equation:
where R _in² is the squared multiple correlation of the model containing all the variables in Eq. (4.8); R _out² is the squared multiple correlation from the reduced model with the term corresponding to β ₂ = 0 removed, i.e., under null hypothesis; m is the number of terms corresponding to β ₂, i.e., number of coefficients being jointly tested; k is the number of predictors in the full regression model, from which R _in² is derived; and n is the number of cases.

Directionality indices are quantified by computation of the absolute or relative predictability improvement obtained by the NARX model compared to the nonlinear autoregressive (NAR) model [8], i.e.,
where represents the residual sum of square (RSS) from its own past (e.g., T ₁ = T ₂ = 0 in Eq. (4.6)) and RSS = . denotes RSS for its own past and the past and present of the input series for a NARX model (e.g., T ₁ ≠ 0, T ₂ ≠ 0 in Eq. (4.6)). The relative causality index for the inputs y ₂ and u ₁ to the output y ₁ is

This index belongs to [0,1]. In the same way, the causality from input y ₁ and u ₂ to output y ₂ can be investigated by reversing the input–output roles of the two series:

Finally, the relative strength of the causal interactions from y ₂ to y ₁ under influence of u ₁ and u ₂ is calculated by the directionality index:

ranges from −1 to 1. A negative value implies that direction of causality is from y ₂ to y ₁, whereas a positive value indicates that the causality is from y ₁ to y ₂, and 0 means balanced bilateral interactions between y ₂ and y ₁.

To apply NSIM for fMRI effective connectivity study, the first step is to build the structure of brain network, although this step is not a requirement to apply NSIM method, because we can identify the system parameters and network structure using model selection method. The reason for using predefined brain network structure is because if we directly apply the method to the whole brain, we have to deal with the computational demanding problem, e.g., selecting the network structure from huge amount of voxels in the brain. If prior brain network structure can be obtained, it will greatly simplify the calculation. As an example for studying fMRI effective connectivity, we build the visual dorsal pathway network as displayed in Fig. 4.2.

Once the network has been constructed, we need to define the region of interests (ROI) from the brain image for the effective connectivity study. ROI can be defined by using standard template method or defined individually by stimulus-driven method. The later method is more accurate and considers individual variability, although it is time consuming to carry out this additional experiment. For instance, we can define the early visual cortex regions, i.e., V1–V4, by means of retinotopic mapping experiments. Because each brain’s ROI contains more than one fMRI time series, we need to choose one representative fMRI response from each ROI for connectivity study. There are two methods to select typical fMRI response for each region, e.g., random sample method and directly average method or using SVD method to average the time series within each region. For the random sample method, we can randomly select fMRI time series from the ROI and use it as representative response for the study. For the average method, we apply the singular value decomposition (SVD) method to get the representative time course for each ROI. In this method, assuming we have m time series (m voxels) in an ROI with n observations (fMRI image frames), we can form a matrix m × n matrix X and apply SVD method as
where SVD is the MATLAB function, and then we keep the largest singular value S(1,1) and set all other elements in S to be 0. Then we do a reverse calculation to get the average representative signal Y as

To demonstrate how to use NSIM for effective connectivity study, we adopt the prior-defined dorsal visual network structure as shown in Fig. 4.2. We employ a random sample method to obtain the representative fMRI time course for each region from the retinotopic mapping experiment to overcome the limitation of averaging which can lead to blur the response delay. To use Eq. (4.5) for effective study, we need to select the order of nonlinearity and autoregression (AR) to model the nonlinear brain system. For reason of simplicity, although automatic AR order selection method is available to achieve this goal, we choose the nonlinearity of 2 and AR order of 1 for the study. Because the nonlinear system structure is defined (Fig. 4.2), we can apply Eq. (4.5) to calculate the system parameters. We employ the method to the retinotopic mapping data for one random selected run. The fMRI response (dotted curve) from V1, V2, and MT is displayed in Fig. 4.3b, d, f, respectively. Using Fourier transform method, we obtain the corresponding fundamental frequency (thin curve) as shown in Fig. 4.3a, c, e, respectively. Since system input (boxcar function denoted by the thick curves in Fig. 4.3a, c, e), output (dotted curves in Fig. 4.3b, d, f), and structure are determined, applying Eq. (4.9), we obtain system parameters as shown in the following equation:

where u ₁(t) (Fig. 4.3a), u ₂(t) (Fig. 4.3c), and u ₃(t) (Fig. 4.3e) are the system inputs. System output is displayed in Fig. 4.3b (y ₁(t)), Fig. 4.3d (y ₂(t)), and Fig. 4.3f (y ₃(t)), respectively. Using the GC tests from Sect. 4.1.3, we calculate the F-test as shown in Table 4.1 for the fMRI time series in Fig. 4.3. In Table 4.1, 4.7106 is the V1 regional self-influence, 0.1604 denotes the regional influence from MT to V2, etc.

Based on Eqs. (4.15, 4.16, 4.17, 4.18, and 4.19), the directionality indices we obtained from V1 to V2 is 0.2562, from V2 to MT is 0.0584, and from MT to V1 is −0.3099. To investigate the input difference, we used the sinusoidal as shown in Fig. 4.3 (thin curve in Fig. 4.3a, c, e); the GC test results are given in Table 4.2. Similarly, the directionality indices we obtained from V1 to V2 is 0.2135; from V2 to MT is 0.2094, and from MT to V1 is −0.4048. Compare Table 4.1 with Table 4.2, we found that F values change using different input function. As we can see, the maximum influence is V1 self-influence in both tables, and the minimum influence is MT region self-influence using both boxcar function and sinusoidal function as inputs.