Welcome to MRI analysis!

Nonlinear system identification theory for fMRI effective connectivity study is mainly based on my previous paper. In brief, a dynamic brain system can be expressed using a black box as shown in Figure 1.

Figure 1. Structure of a nonlinear dynamic brain system. A, brain system block diagram, the grey box represents the nonlinear brain system; B, brain system input for experimental design; C, brain system output (fMRI response) curve from MRI scanner in response to the system input .

In this nonlinear system, the nonlinear brain black-box is depicted in Figure 1A and the input of the black-box is the experimental design (Figure 1B), and the output (Figure 1C) is the fMRI response with some random noise added. Figure 1B is one example of a brain system input from standard block designs. The system output or fMRI response is adaptively changed according to the input as shown in Figure 1C. Formally, the physiological processes underlying the BOLD response can be modelled as a multiple-input and multiple-output (MIMO) system:

, (1)

and its discrete form is

(2)

where and are nonlinear functions, and represents the set of model parameters. Under some mild assumptions the discrete-time multivariate system (1) with outputs (e.g., response from different regions) and inputs (e.g., different types of stimuli) can be described by an autoregressive moving average with exogenous input (NARMAX) as follows:

(3)

where ,,, are the system output, input and noise, respectively; ,,and are the maximum lags in the output, input, and noise; is a zero mean independent sequence; is a new nonlinear function which can be obtained from nonlinear functions and . A special case of the general NARMAX model (3) is the nonlinear autoregressive with exogenous inputs (NARX) model which can be expressed as:

(4)

By applying the regression equation, the NARMAX model (3) and NARX model (4) can be approximated as:

, (5)

where ; for ,,; is the number of nonlinear terms;is the system order; is the total number of time point in the time series; is the number of connected regions; is the number of inputs. Equation (5) denotes a general case where both input and output terms may be present, but it should be understood that some of the may contain only input or output terms and cross-products (if using polynomials base to approximate the output). For example, for two stationary series of values, the inputs and, output of a closed-loop time–invariant nonlinear brain system can be described as:

+ (6)

where the coefficients, , and denote constant (zero-th order), linear (first order), and nonlinear (second order) contributions to , respectively. represents the experimental input, and is the prediction error of . The model orders and 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 and . The model can be represented in the matrix form as :

(7)

where the vector contains values of output series, is the prediction error series. is the experimental input time series; and are the first order vector coefficients; , and are the second order vector coefficients. The matrices and contain the linear AR terms and the (+1) linear exogenous terms respectively:

the matrix contains of the quadratic AR terms given by the product of the terms of the matrix . Similarly, the matrix contains the ()()/2 quadratic exogenous terms, and the matrix contains the cross-terms. Equation (7) can be written in the compact matrix format as :

(8)

where , . Coefficient matrix can be estimated by least squares:

(9)

where is the Moore-Penrose pseudoinverse of the matrix. By neglecting the nonlinear terms , experimental input , and considering only the first order of AR, i.e. AR(1), this leads to:

(10)

(11)

(12)

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

Once the coefficients of the model are determined (e.g., equation (9)), Granger causality tests are derived based on F/T statistics. For simplicity and illustrative purposes, we take the two-connection nonlinear models (6) as an example. The same principle can be applied for the linear system (equations (11) and (12)). The test for determining Granger-cause (GC) is:

(i) is GC of if in equation (6) is not true. Given the data, we reach this conclusion if is rejected.

(ii) Similarly, Granger causes of can be investigated by reversing the input-output roles of the two series.

T and statistics are developed to detect significant relations. From equation (8), and partitioning the coefficients as and accordingly, we can write this test as:

versus , (13)

with the maintained hypothesis. Therefore, T statistics can be applied for testing the hypothesis. For the F test, the following equation:

(14)

can be used, where is the squared multiple correlation of the model containing all the variables in equation (8); is the squared multiple correlation from the reduced model with the term corresponding to removed, i.e., under null hypothesis; is the number of terms corresponding to , i.e., number of coefficients being jointly tested;is the number of predictors in the full regression model, from which is derived; is the number of cases.

Figure 2. A 3-connection brain network system

Nonlinear model for fMRI effective connectivity study: an example for a 3-connection system.

Before we apply model selection algorithm, we analyze the nonlinear model from equation (5). Equation (5) can also represent a standard Volterra series expansion to approximate a nonlinear biological system. For example, for a 3-connection network (3-connection regions , , and as shown in Figure 2) with 2 order of nonlinearity and lags (AR order = 2) on the input of 1, the output can be expressed as:

(15)

where is the coefficient for the constant drift (intercept term), and

(16)

is the nonlinear function for the mapping between output and AR inputand. Similarly,

(17)

(18)

(19)

(20)

(21)

and

(22)

in this study. Apart from the linear covariates as shown in the first two terms of the equations (16), this model also includes nonlinear covariates. For a second order nonlinearity with an AR order of 2, equation (15) includes 31 covariates, i.e., 6 covariates of linear AR, 21 covariates of the nonlinearity which is from the combination of these 6 linear AR covariates, 3 inputs, and one constant drift. Therefore it is necessary to develop model selection methods to determine the linear/nonlinear function in equation (3), efficiently and accurately.

Model selection for NSIM in effective connectivity study

Consider the term selection problem for the linear-in-parameters model (5). Let be a vector of measured fMRI response (system output) at time instants, and be a vector formed by the th candidate model term, where = 1, 2,…, . Let be a dictionary of the candidate bases from equation (5). From the viewpoint of practical modelling and identification, the finite dimensional set is often redundant as shown in equation (4). The model term selection problem is equivalent to finding a full dimensional subset of () bases, from the library, where , and , so that system output can be satisfactorily approximated using a linear combination of ,,,as below :

(23)

or in a compact matrix form

(24)

where the matrix is assumed to be of full column rank, is a parameter vector and is the approximation error vector; where is the total variables/covariates in equation (24), i.e., for a 3-connection model with second order of AR (for ,, and) and nonlinearity (equation (15)). For a 3-connection network, we set, , and equals to the corresponding AR covariates in the equations (15-22), i.e., ,,;,,;;, which is an matrix, where is the number of fMRI time frames after excluding the order of AR.