Model and algorithm
Particle filter [17, 19,20,21,22] is a sequential Monte Carlo methodology based on the Bayes theorem on conditional probability. Particle filters estimate the probability distributions of hidden variables of interest, modeled according to a hypothesized statespace equation. The probability density function (pdf) of the hidden variables is allowed to be timevarying and is therefore sequentially updated when new data become available. Such probability distribution is estimated from the data, modeled according to a hypothesized observation equation. In brain connectivity studies based on fMRI data, the relationship among the timeseries of R different brain Regions of Interest (ROIs) (mathbf{x }_t={x_1(t),dots ,x_R(t)}) can be modeled as a first order linear Vector Autoregressive (VAR) model [12, 21, 23,24,25] as:
$$x_i(t)=sum _j^R a_{ij}(t)x_j(t1)+eta _i(t) quad i=1,ldots ,R $$
(1)
or in matrix notation:
$$begin{aligned} mathbf{x }(t)=mathbf{a }(t)mathbf{x }(t1)+{varvec{eta }}(t) end{aligned}$$
(2)
where
$$begin{aligned} mathbf{a }(t)= begin{bmatrix} a_{11}(t) &{} a_{12}(t) &{} dots &{} a_{1R}(t) \ a_{21}(t) &{} a_{22}(t) &{} dots &{} a_{2R}(t) \ vdots &{} vdots &{} ddots &{} vdots \ a_{R1} (t) &{} a_{R2}(t) &{} dots &{} a_{RR}(t) end{bmatrix} end{aligned}$$
(3)
which is employed as the observation equation describing the relationship between the observations (mathbf{x }(t)) at time t and those at time (t1) (that is, (mathbf{x }(t1))); ({varvec{eta }}(t)) is the vector of observation noise; the matrix of hidden parameters of interest (mathbf{a }(t)) represents the causal influence exerted between different areas, and its elements (a_ij(t)) are the coefficients which represent conditional dependence. In particular, it can be assumed that elements of (mathbf{a }(t)) are allowed to be timevarying:
$$begin{aligned} a_{ij}(t)=a_{ij}(t1)+nu _{ij}(t) end{aligned}$$
(4)
where (a_{ij}(t)) is the ijth element of the matrix (mathbf{a }(t)), describing the influence of the jth region over the ith region, and (nu _{ij}(t)) is the process noise (innovation) term.
The adoption of a linear model was supported by the wellestablished body of literature on fMRI data modelling and brain connectivity analysis at the temporal scales of fMRI data [25, 26]. The adopted autoregressive model was firstorder, which was optimal on the basis of the Schwartz criterion, in accordance with previous findings [12, 25, 27].
The PF algorithm evolves from an initial probability distribution for (a_{ij}(t1)), which we chose to be uniform at (t=1), and through Eq. (4) it generates new possible values for (a_{ij}(t)). The N particles are generated from the probability distributions of the elements of (mathbf{a }_i(t)): the distributions are adapted with each new observation through the mechanism provided by particle filtering, to describe for the set of coefficients (a_{ij}) at every timestep. The algorithm generates N particles, by updating those at the previous time point (initialized to zero at t=0) using a noise innovation term, by following Eq. 4. In our implementation, the innovation values are drawn from a Gaussian distribution. The standard deviation of this distribution follows the absolute difference between the two previous coefficient values estimated by the algorithm and constrained between a minimum of 0.1 and a maximum of 0.4: these values were chosen empirically to prevent divergence and, at the same time, capture nonstationarities.
With Eq. (2), the PF algorithm generates predicted values of the observations at time t. The desired probability density function of the parameters of interest (mathbf{a }(t)) can be estimated via Bayes theorem as follows:
$$begin{aligned} p(mathbf{a }(t)mathbf{x }(1,dots ,t))=frac{p(mathbf{x }(t)mathbf{a }(t)) p(mathbf{a }(t)mathbf{x }(1,dots ,t1))}{p(mathbf{x }(t)mathbf{x }(1,dots ,t1))} end{aligned}$$
(5)
and with the assumption of Gaussian noise we have
$$begin{aligned} p(x_i(t){mathbf {a}}_{i}(t))=frac{1}{(2pi sigma _eta ^2)^{R/2}} text {exp}Big (frac{(x_i(t){hat{x}}_i(t))^2}{2sigma _eta ^2}Big ) end{aligned}$$
(6)
where ({hat{x}}_i(t)) are the data estimated through Eq. (1) at time t for the ith ROI and ({mathbf {a}}_i(t)={a_{i1},dots ,a_{iR}}) is the vector of hidden variables associated with the ith ROI at time t, that is, the ith row of matrix ({mathbf {a}}_t). In Eq. 6 the value of the (mathbf{a }_i(t)) elements determines the value of the estimate ({hat{x}}_i(t)), and the weight of the ith particle.
In most applications, Eq. (5) cannot be solved analytically [28], but it can be computed through the Sequential Monte Carlo sampling scheme, which consists in representing the pdf (p(mathbf{a }(t)mathbf{x }(1,dots ,t))) as a discrete set of N particles:
$$begin{aligned} p(mathbf{a }_i(t)mathbf{x }(1,dots ,t))approx sum _{n=1}^Nw_t^{(n)} delta (mathbf{a }_i(t)mathbf{a }^{(n)}_i(t)) end{aligned}$$
(7)
where (w_t^{(n)}) is the weight associated to the nth particle vector (mathbf{a }^{(n)}_i(t)) for the ith row of matrix (mathbf{a }(t)) at time t. The Sequential Importance Sampling (SIS) [28] methodology provides a strategy to compute the weights. It has been shown [20] that the weights can be sequentially updated as follows:
$$begin{aligned} w_t^{(n)}propto w_{t1}^{(n)} p(mathbf{x }(t)mathbf{a }_i(t)^{(n)}) end{aligned}$$
(8)
where the proportionality takes into account normalization factors. With this approach, at each time instant t we have a sample set ({mathbf{a }_i(t)^{(n)},w_t^{(n)}}) for (n=1,dots ,N) and for (i=1,dots ,R) which can be used to estimate the pdf of the parameters and to infer information about the network. However, after some iterations, most of the particles will have a very low statistical weight, resulting in a lower exploration efficiency of the algorithm. To overcome this typical problem of sequential Monte Carlo methodologies, a step called Resampling is performed. The number of effective particles was defined in [29] as
$$begin{aligned} N_{text{eff}}=frac{1}{sum _{n=1}^N(w_t^{(n)})^2} end{aligned}$$
If (N_{text{eff}}) is below a certain arbitrary threshold the Resampling is performed: particles with weight below a certain threshold are substituted by copies of particles with sufficiently high weights and to each of the new particles set is assigned the same weight 1/N. This results in a more effective exploration of the solution space, because only statistically relevant particles remain after this step. The VAR process estimation problem using SMC was developed in [30], and it was extended and applied to timevarying gene expression networks in [17].
To sum up, the resulting algorithm can be schematically expressed as in Table 1. In our implementation the procedure is repeated (N_r=100) times, all independent from each other, to provide a better exploration of the solution space, and resampling was performed when (N_{text{eff}}<30%) of the total number of particles. The final outputs of the algorithm are the (mathbf{a }_t) computed as the average of the (N_r) repetitions. In this implementation, the running time was proportional to the length T of the time series, the number of particles N, the number of repetitions (N_r), and it increased quadratically with the number of network nodes [31]. The algorithm was implemented in MATLAB (Mathworks, Natick, MA, U.S.A.) R2017b.
Synthetic data
To validate the proposed approach, two different synthetic networks were used.

One network with (R=6) nodes, each with (T=100) time points, stationary coefficients generated with the MATLAB function varm() with a SignaltoNoise Ratio (SNR) set to either (infty ) ((sigma _eta =0), ideal case) or 6dB.

Another network with (R=2) and (T=250) was used to assess the PF capability to capture timevarying hidden parameters. In this case, (a_{ij}) coefficients were zero except for coefficient (a_{21}), whose value switched from 1 to (1) with a period of 125 time points. The SNR was 10dB.
The first synthetic dataset allowed us to verify the reliability of the results of the PF and to decide the optimal values of the parameters; the second one allowed us to address the capability of the PF to track variations through time of the AR coefficients.
Real fMRI data
The proposed approach was also retrospectively applied to real fMRI data acquired on healthy volunteers in two different experimental setups, with two and four participants respectively, acquired with twodimensional singleshot echoplanar imaging (EPI) on a 7T MRI system (MR950, GE Healthcare, Chicago, IL, U.S.A.).

Motor task Timeseries consisting of 240 time points with a temporal resolution of 2s were acquired on two subjects with the following acquisition parameters: Time of Echo ((hbox {TE}) = 23,hbox {ms}), Flip Angle ((hbox {FA}) = 60^{circ }), Field of View ((hbox {FoV}) = (192,hbox {mm})^2), acquisition matrix size = (128times 128), 32 slices of thickness = 1.5mm, resulting in isotropic spatial resolution of ((1.5 ,hbox {mm})^3) and Time of Repetition ((hbox {TR}) = 2,hbox {s}). During acquisition, the subjects’ thumb and indexfingertips were stimulated via a pneumatic device (Linari Engineering, Pisa, Italy). The subjects’ task was to move the finger whenever it was stimulated. The fMRI data were motioncorrected using MCFLIRT [32]. Spatial smoothing was applied by using a Gaussian kernel of FWHM 3.0mm; each 4D dataset was demeaned and normalized by a single multiplicative factor; highpass temporal filtering was applied to remove slow temporal drifts of the fMRI signal. Four ROIs were studied covering primary somatosensory (S1), primary motor (M1), supplementary motor (SM) and parietal (P) cortices. All ROIs consisted in four voxels and were manually drawn on each subject on one slice only, to avoid potential slice timing confounds (Fig. 1). The resultant timeseries of each ROI were obtained by averaging the four timeseries of individual voxels.

Visual task Timeseries with a temporal resolution of 0.8s were acquired on four trials of either 300 time points (two subjects) or 600 (two subjects). Scanning parameters were (hbox {TE} = 21,hbox {ms}), (hbox {FA} = 48^{circ }), (hbox {FoV} = (192,hbox {mm})^2), acquisition matrix size (= 64times 64), 22 slices of thickness (= 3,hbox {mm}), resulting in isotropic spatial resolution of (3 mm)(^3) and (hbox {TR} = 800,hbox {ms}). The same preprocessing steps described above for the motor task, including motion correction, spatial smoothing, normalization and signal drift removal, were adopted. Subjects underwent a periodic visual stimulation alternating between black and white dots moving along spiral trajectories over a gray background (stimulation ON) and presentation of the gray background alone (stimulation OFF). Four ROIs were studied covering the Lateral Geniculate Nucleus (LGN), the Middle temporal cortex (MT), the Primary Visual area (V1) and one control ROI in the tempoparietal cortex (CTRL). Four voxels wide ROIs were manually drawn on each subject on one slice only.
The optimal order of the autoregressive model describing the time series was 1, as estimated by the Schwartz criterion [12, 25, 27].
The particle filter was applied not only to the original time series but also to the same fMRI data shuffled in time. This test allowed to address the effective dependency of the results from the temporal order of the data, i.e. to address causal dependency through time in the data. Subsequently, the (a_{ij}) coefficients estimated by particle filtering were compared to the delayed correlation (DC) (c_{ij}) between signals (x_i(t)) and (x_j(t1)) which reflect the timeinvariant conditional correlation between network node (ROI) j over node i. Furthermore, the results of the PF were also compared to coefficients estimated in a stationary framework fitting a multivariate linear regression model with stationary coefficients to the data.
With a t test between coefficient values in the presence and in the absence of stimulation we searched statistically significant (pvalue (<0.05)) changes in connectivity following the underlying stimulation pattern. Also, as a control analysis, a second ttest was performed, simulating a different stimulation pattern from the actual one, which allowed to exclude variations deriving from spurious fluctuations of the results, as explained in Fig. 2.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Disclaimer:
This article is autogenerated using RSS feeds and has not been created or edited by OA JF.
Click here for Source link (https://www.springeropen.com/)