Advances in the sophistication of imaging techniques necessitate the development of techniques to model the neural and cognitive phenomena they represent. Using electrocorticographic (ECoG) data, we propose a data-driven approach using ordinary differential equations and Bayesian differential structural equation modeling (BdSEM) to model effective connectivity related to sensorimotor integration. First, we tested the region-based covariance structure across subjects by each experimental condition to evaluate the tenability of pooling subject data to perform group level versus single-subject analyses. Second, we applied a differential equation approach to model dynamic change originating from regional neuronal states across the data acquisition period. Finally, we employed an information-theoretic search strategy to identify the optimal connectivity model within each experimental condition for a single subject. Results of subject-specific (intra-individual) relationship maps include effective, contemporaneous and delayed effective connections of across different brain regions.

ECoG, electro cortico graphic; BdSEM, bayesian differential structural equation modeling; ODEs, ordinary differential equations; ROI, region of interest; fMRI, functional magnetic resonance imaging; MCMC, markov chain monte carlo; FDA, functional data analysis; STG, superior temporal gyrus; PostSTG, posterior STG; MidSTG, middle STG; BIC, bayesian information criterion; AIC, akaike information criterion; FF, feed forward; FB, feed back; ERP, evoked response potential

To date, imaging neuroscience and electrophysiology have provided a solid foundation for functional specialization as a principle of brain organization in humans. However, accurately modeling functional integration of specialized areas of the brain has proven to be a more difficult task.¹ For example, the human brain can be viewed as a holistic and dynamic system involving functionally specialized areas or regions that are related by effective connections during cognitive processing. Connectivity models (functional or effective) have become useful approximations to holistic and dynamic systems for understanding underlying relationships between regions of neural activation. For example, functional connectivity models provide relational maps based on statistical dependency (correlations) between remote neuro physiological events. Effective connectivity models involve estimating the influence that one neuronal system exerts over another, either at the synaptic or population level.¹ Furthermore, effective connectivity is activity-dependent involving interactions among regions of the brain. Finally, modeling intrinsic connectivity involves capturing dynamic elements of the system reflecting the essential nature of the collective system.

One approach to modeling effective connectivity is based on intracranial electrocorticographic (ECoG) time series data. In ECoG-based studies, data are acquired from multi-contact subdural electrodes implanted during surgical evaluation in epilepsy patients.  ECoG provides a higher signal to noise ratio and temporal and spatial resolution than EEG due to the fact that ECoG signals are not influenced by the low conductivity of the skull and the fact that measurements are acquired in the vicinity of the underlying brain sources.² In clinical practice, ECoG has become the “gold standard” for defining epileptogenic zones.³ The locations of ECoG electrodes implantation are determined using pre-operative clinical data such as ictal semiology, ictal and interictal features on scalp EEG, and structural MRI findings.⁴ The pattern of the subsequent ECoG recordings is then used for further localization of epileptogenic foci.

In this study, we present an innovative modeling approach for the development and estimation of human electrophysiological inferential connectivity maps using ECoG data. The aim of the method is to obtain accurate representations of underlying effective relationships while also considering intrinsic or dynamically changing aspects of regional components of brain activity. Specifically, we use ordinary differential equations within a structural equation modeling framework and a heuristic model search strategy to provide insights into electrode-to-electrode connectivity using ECoG. Here, our interest is based on modeling simultaneous or dynamic change in neural activity so as to provide unique insights into the causal nature of the relationships between brain regions in specific vocalization and auditory processes. Given the high spatial and temporal resolution of ECoG data, inferences about underlying neural functions as mapped on to brain anatomy allow for exquisite representation of both time and space in the human brain.^3,5,6 

The organization of this manuscript is as follows. The first section briefly reviews SEM and the advantages of using Bayesian SEM in studies such as ours. The second section details the application of ordinary differential equations (ODEs) to measure simultaneous change in brain activity propagating from the level of neural activity upward to the region of interest (ROI) level. The third section introduces a model for Bayesian SEM based on differential equations (BdSEM) and describes the heuristic specification search procedure used to arrive at optimal models of effective connectivity. The fourth section outlines a simulation study conducted to evaluate the efficacy of the modeling strategy.

Structural equation modeling: A frequently used technique used for modeling population-level relationships among regions caused by neural activation in brain regions is structural equation modeling (SEM). SEM was first adapted to imaging data by McIntosh & Gonzales-Lima^7,8 to examine effective⁹ relationships between regions of interests (ROIs). To date, SEM of neural connectivity has predominantly been applied to functional magnetic resonance imaging (fMRI). In part, this trend may reflect difficulties in source localization that are inherent in most electrophysiological data capture methods.

Application of SEM involves developing a set of simultaneous equations to estimate

1. The regression coefficients between measurements on observed variables and associated latent variables (i.e. a measurement model) and
2. Relationships among latent variables comprising a hypothesized model. SEM is based on the general linear and/or generalized nonlinear model and is very useful for complex or longitudinal data structures (e.g., multilevel random coefficients models, growth curve modeling, and differential equation-based manifest or latent change). Perhaps SEM is most often associated with confirmatory (theory confirming) modeling approaches. However, SEM can be used to conduct heuristic search techniques in high dimensional data structures to locate an optimal model from among a set of competing ones. Heuristic search techniques leverage information-theoretic algorithms from the field of artificial intelligence (e.g., mathematically-based approaches based on lagrangian heuristic/incomplete branch-and-bound algorithms). Application of SEM involves measuring the discrepancy between a parameterized causal structure of hypothesized relationships (e.g., between observed or latent variables or a combination of both) and estimated parameters from data as measured through a series of fit indices.^10,11  SEM employing Bayesian probability and statistics^12,13 has recently emerged and provides increased precision and flexibility in modeling scenarios with high dimensional data structures, small sample size and non-normally distributed variables.

Statistical inference and bayesian learning: The history and development of Bayesian statistical methods are substantial and closely related to frequentist statistical methods. In some ways, Bayesian statistical thinking can be viewed as an extension of the traditional (i.e., frequentist) approach, in that it formalizes aspects of the statistical analysis that are left to uninformed judgment by researchers in classical statistical analyses.¹⁴ Bayesian methods include data analytic techniques that are derived from the principles of Bayesian statistical inference. Statistical induction involves learning about the characteristics of a population from a subset of members of a particular population. Numerical values of populations are expressed as parameters (θ) while numerical values of the subset of the population are expressed as (y). Given the numerical values (y) in the subset or sample dataset, uncertainty is reduced about the population parameters. Quantifying this shift in uncertainty is the goal of Bayesian inference. The parameter space Θ is the set of all possible parameter values from which we wish to identify the value(s) that best reflect the true population parameters (e.g., regression weights in connectivity models). Bayesian learning involves a numerical formulation of the joint beliefs about y and θ, expressed as probability distributions over y and θ. In short, Bayesian learning involves the components listed below.

1. Given each value $\theta \in \Theta$ the prior distribution $p(\theta )$  describes the belief that $(\theta )$ represents the true population parameter.
2. Given each value $\theta \in \Theta$ and $y\in Y$  the sampling model $p(y|\theta )$  captures the belief that $y$  will be the outcome of a particular study if we knew $\theta$  to be correct or true.
3. After $y$  is acquired, each numerical value of the posterior distribution $p(y|\theta )$  describes the belief that $\theta$  is the true value having observed dataset $y$ .

The posterior distribution for model parameters is derived using Bayes rule as in Equation 1.

$p(\theta |y)=\frac{p(y|\theta )p(\theta )}{{\displaystyle {\int }_{\Theta }p(y|\tilde{\theta })p(\tilde{\theta })d\tilde{\theta }}}$     (1)

A core component of Bayesian modeling is accurately modeling the generating process (e.g., unknown causes modeled using a prior distribution) and unknown population parameters $\theta$  to observed sensory data values ( $y$ ). Bayesian statistical methods are particularly well-suited for developing generative or recognition models of complex systems^15-17 because the goal in Bayesian learning is to model the generating process that produced the observed data values. For example, in generative and recognition models, functions are applied thereby allowing a mechanism for mapping causes to sensory input. Specifically, the goal of generative modeling is to: “learn representations that are economical to describe but allow the input to be reconstructed accurately”.¹⁸ The goal is to make inferences about the causes and learn the parameters. Bayesian probability provides a natural framework for linking unknown parameters and causes to observed data. Finally, in the classical school of probability, the sample data values are selected randomly with the statistics estimated being fixed point estimates of population parameters. Conversely, in the Bayesian school of probability and inference, the sample is fixed (i.e. not considered random) with the parameters estimated being random (e.g., obtained using Markov chain Monte Carlo [MCMC] resampling methods). Because our estimated parameters are random variables, we can make probabilistic statements about their certainty with a high level of precision.

In the Bayesian modeling approach, we view any unknown quantity (e.g., parameter) as random and these quantities are assigned a probability distribution (e.g., normal, Poisson, multinomial, geometric, etc.) that provides the impetus for generating a particular set of data. In this study, our unknown population parameters were modeled as being random and then assigned to a joint probability distribution. In this way, we were able to summarize our current state of knowledge regarding the model parameters. The sampling-based approach to Bayesian estimation provides a solution for the random parameter vector θ by estimating the posterior density or distribution of a parameter. This posterior distribution is defined as the product of the likelihood function (accumulated over all possible values of θ) and the prior density of θ.¹⁹ In our case, a generative model is specified in terms of a prior distribution relative to the neuronal activity (i.e. the causative mechanism of the observed data).

Modeling long time series relative to small sample size: Developing statistical models for studying populations relevant to neuroscience often poses considerable challenges due to small sample sizes, the issue of low statistical power, and the length of the time series (i.e. number of repeated measurements) being quite long (e.g., > 100 time points). One analytic approach appropriate for the challenging scenarios previously noted are multivariate autoregressive models.^20-22,13 The use of Bayesian statistical modeling is more sensitive for hypothesis testing and interval estimation than frequentist approaches when there are small sample sizes with multivariate structure and a long time series.^23-25,19 In the present study, the number of time series measurements is large (i.e. > 1000) with the length of the time series greater than the sample size (N=1). Although our data acquisition involved multiple trials for each electrode within each subject, the average time series (waveform) for each site within each subject was used for BdSEM modeling.

Modeling fMRI data involves capturing the degree of deoxygenated vascularization relative to baseline for each voxel or congregation of voxels corresponding to ROIs across multiple trials (i.e. time series) nested within subjects. Each case within the time-series data structure for each ROI corresponds to the peak, or averaged peak, of hemodynamic activity for that trial. Furthermore, after slice timing correction, the resulting data structures for each ROI are obtained contemporaneously for any given trial. When extended to multivariate regression or multivariate autoregressive models (or SEMs), mapping hypothesized effective connectivity relationships (including those with theoretical support) becomes fallacious given that inferences of causal influence of one region on another requires temporal precedence. For example, the empirical conditions for inferring a causal relationship between two variables include (a) X is related to Y, (b) X temporally precedes Y, and (c) the relationship between X and Y is not mediated by a third variable - Z.^15,26  Although more complicated, the same rules apply when modeling the relationships in a system of simultaneous equations (e.g., more than two variables).

Kim et al.²⁷ developed a Unified SEM technique of fMRI time-series using data on visual attention using multivariate auto regression to model both contemporaneous and temporal (longitudinal) relationships simultaneously. However, the effective relationships determined from delayed representations and real-time representations of data were separated by a single trial with a time resolution of 3 seconds. Moreover, inferences of interregional causal relationships between multiple regions’ underlying electrophysiology were made from vascular data that was obtained 3 seconds prior. Also, Unified SEM requires relationships between parameters to be specified a priori. Gates et al.²⁸ improved upon this technique, by developing an Extended Unified SEM technique integrating an automatic search procedure based on Lagrangian multiplier tests, or modification indices. In this paper, we improve on these procedures by the estimation of electrophysiological connectivity maps to obtain more accurate representations of effective relationships by leveraging the high temporal and spatial resolution of (ECoG), a form of electrophysiological data collected from electrodes placed directly on the brain cortical surface during neurosurgery.^2,3,5,6  ECoG remedies issues of localization frequently encountered in EEG while retaining the advantage of high temporal resolution (0.5 – 1.0ms) and spatial resolution (diameter 2.0-3.0mm). The improved temporal resolution relative to fMRI makes ECoG ideal for reducing bias to very small or optimal levels specific to inferences based on underlying neurophysiology. Leveraging the high quality of ECoG data, we present a novel method for modeling dynamic data structures of ROIs to be input in a data-driven method of modeling the relationships between regions of the cortical voice network. 

Differential structural equation modeling: Modeling dynamic intraindividual change using differential equations based on long time series is a rapidly evolving technique. Functional data analysis (FDA)²⁹ provides an approach for fitting differential equations directly from acquired measurements in studies of human growth and EEG. Using differential equations within SEM provides an approach for applying the FDA approach to model the effective connectivity expressed as dynamic coupling between the neuronal states of various brain systems when exposed to experimental conditions. Typically, in neuroimaging studies, experimental conditions are modeled as inputs via boxcar or stick functions. Here we approximate the (reciprocal) relationships among brain region activity using a bilinear approximation.

Ordinary differential equations (ODEs): Using ODEs to model neuronal activity provides a way to express the rate of change of the states as a parameterized function of the states and inputs. In Equation 2, $\text{u}(\text{t})$  represents a particular experimental stimulus as an input.  Here we use ODEs to model instantaneous changes in neuronal activity in two steps.  First, we employ neuronal state equations thereby linking the derivatives of neuronal states $\text{x}(t)=(x_1(t),\dots ,x_d(t))\text{'}$  of d brain regions to themselves under the influence of an experimental condition.² The ODE-based dynamic model assumes a Markov property whereby instantaneous changes of the system depend only on system states and experimental inputs at the same moment in time. Importantly, the Markov property is tenable for the brain system performing simple auditory and vocalization experimental tasks as in our study. The ODE-based state equation detailing the dynamic changes of neural states is provided in Equation 2.

$\frac{\text{d}\text{x}\text{(t)}}{\text{d}t}={\text{F}}_1(\text{x}\text{(t)},\text{u}(\text{t}),{\text{θ}}_1)$     (2)

Where F1 is a set of nonlinear basis functions capturing neuronal influences that specific brain regions $\text{x}\text{(t)}$  and experimental stimuli $\text{u}(\text{t})$ exert on other specific regions. The vector ${\text{θ}}_1$  contains the unknown parameters in the system. In Equation 2, $\text{x}\text{(t)}$  is a continuous function reflecting an average of the neuronal activity in a specific region. Equation 3 conveys the observation-level (output) equation that enables a description of how the underlying neuronal activity causes changes in the observed data vector $y$  in each ROI.

  $y\text{(t)=}{\text{F}}_2(\text{x}(\text{t}),{\theta }_2,\text{e}(\text{t}))$    (3)

Where ${\text{F}}_2$  is an unknown function, ${\theta }_2$ are parameters to be estimated, and $\text{e}(\text{t})$ are error terms. Equations 2 and 3 are requisite to modeling any dynamic system (i.e. they link underlying or hidden states to observable outputs).³⁰ Particularly relevant to our study is that the causal relationships between the outputs and the inputs conform to a Volterra series, which expresses the outputs as a generalized convolution of the input – without reference to the hidden states $\text{x(t)}$  .¹ In short, the Volterra series is a functional Taylor expansion of the outputs with respect to the inputs.

Next, in Equation 4 we are able to include a bilinear approximation to F1 by combining elements in Equations 2 and 3. In Equation 4, causal influences among regions are possible to elucidate since bidirectional estimates of parameters are included in the ${\text{A}}_{i_1{i}_2}$ and ${\text{A}}_{i_1{i}_2}$ portion of the equation. These bidirectional influences include the dynamic effect of the time series. The parameters in the $B$  matrix enable the estimation of stimulus-dependent effective connectivity between component regions. Parameters are allowed to vary over the time course of the data capture (ECoG signal) in order to ensure accurate approximation to F1.

  $\frac{\text{d}\mathrm{x}(\text{t})}{\text{d}(\text{t})}=Ax(\text{t})+{\displaystyle \sum _{j=1}^J\text{u}}(\text{t})\cdot B_{j}x(\text{t})+Cu(\text{t})+D$     (4)

Here $A$  = ${({\text{A}}_{i_1{i}_2})}_{\text{dxd}}$  with ${\text{A}}_{i_1{i}_2}$ denotes the effect of component ${\text{i}}_2$  on component ${\text{i}}_1$  exerted at the current state; $B_j$ = ${({\text{B}}_{\text{j,}{\text{i}}_{\text{1}}{\text{i}}_{\text{2}}})}_{\text{dxd}}\text{,j=1,…,J,}$  couples the j^th stimulus with the neuronal states and nonzero ${\text{B}}_{\text{j,}{\text{i}}_{\text{1}}{\text{i}}_{\text{2}}}$  implies that the effect exerted by component  ${\text{i}}_2$ on component ${\text{i}}_1$ depends on stimulus j; $C={({\text{C}}_{\text{ij}})}_{\text{dxJ}}$  with ${\text{C}}_{\text{ij}}$  indicated the effect of stimulus j on component i; and $D=({\text{D}}_{\text{1,…,}}{\text{D}}_{\text{d}})\text{'}$  with ${\text{D}}_{\text{i}}$  representing the intercept for component i.²

To estimate parameters at the ROI (observation) level, we use Equation 5.

$y\text{(t)=}x\text{(t)+}e\text{(t)}$    (5)

Where $e\text{(t)}$  is a d-dimensional vector of measurement errors with mean error of zeros.  Further, we assume that the distribution of errors is normally distributed with mean zero and variance ${\sigma }_{\text{i}}^{\text{2}}$ . Finally, to estimate the differential equations, we use basis function expansion as described in Ramsay & Silverman.³¹ The basis function expansion approach is appropriate due to the temporally dense observations of brain regions and the average of a large number of neuron activities in the region of interest. The basis function approach allows for closed form solutions based on a large number of neural activity and the associated derivatives of the regional activity. Estimating parameters in the bilinear model is straightforward and reduces to solving d linear regression equations.

Analytic strategy: Here we use a three-step approach to achieve an optimal model of effective connectivity.  First, we fit differential equations to ECoG signal data to capture instantaneous change of neuronal activity within each region of interest based on data acquired via target subdural electrode sites. Second, we developed a BdSEM to model effective connectivity of the brain system using derivatives and observational data of each ROI. Third, we use a Bayesian information-theoretic algorithm to identify the optimal network model (i.e. the model with the highest posterior probability of being optimal). Our approach is flexible in that one can (a) use the entire ECoG electrode array to model signals from a large amount of regional activity (e.g., 48 to 96 electrodes within a specific area of the brain) propagating from dense neuronal activity, or (b) use one ECoG electrode (single channel) within the array to model a targeted brain region. The first example above begins by modeling the neuronal dynamics and is similar to Dynamic Causal Modeling (DCM, Friston et al.³²) because the brain is viewed as a continuous-time dynamic system of neuronal activity where signals propagate upward to an observational or population level. DCM incorporates two equations: a bilinear equation capturing neuronal activity in a bidirectional, nonlinear fashion, and an observation-level equation enabling effective connectivity modeling. The approach presented here is viewed as a special case of a DCM.

Modeling complex relationships of neural systems based on functional and effective connectivity poses unique challenges to neuroscientists. In this paper we presented an approach to modeling effective and intrinsic relationships between ROIs simultaneously within a functional network. A primary goal of our work included systematically addressing issues specific to effective and intrinsic connectivity in a way that yielded an accurate map of vocal sensorimotor integration within and between ROIs. This is particularly important in dealing with electrophysiological data from neurosurgical patients because of the exquisite anatomical localization and temporal resolution.

The modeling approach presented here is aimed toward increasing an understanding of the causal relationships observed among electrodes during different conditions that implicate various stages of sensorimotor vocal control. Our approach is geared to the development of network models that can explain the neural control of human vocalization.  Specifically, our work is focused on understanding feed forward (FF) and feedback (FB) control mechanisms using a paradigm by which we alter auditory feedback to subjects in real time while they are vocalizing. Such data inform specific processes such as self-voice identification and FB based error detection.

Our work proceeded by first demonstrating an approach for estimating delayed data structures relative to each ROI using differential equations. Second, real-time and delayed representations were input for each region into a baseline or null model that assumed no relationships. Next, an heuristic iterative search algorithm was used to iteratively specify directional paths between regions. Our method was demonstrated on 5 regions of the voice network using ECoG data. The approach presented here offers a systematic and comprehensive way to model effective and intrinsic relationships and appears to work well for electrocorticography (ECoG), a form of electrophysiological data collected from electrodes placed on the surface of patents' cortex during neurosurgery.

While this is not an applications paper, it is useful to determine the nature of the relations among the ROIs modeled here. In this work we used data for which subjects vocalized an “ah” at a comfortable loudness and pitch feel or listened to their own vocalizations (passive listening). Based on our previous work with fMRI, ERP and ECoG we selected regions that are implicated in sensory control of the human voice. What is novelabout the use of ECoG data is that we were able to divide STG into multiple regions in order to determine functional connectivity associated with the percolation of this region, a task that is highly challenging with fMRI or scalp based ERP data. Since we know that the STG plays a central role in the sensory motor control of the voice, this region likely serves as the hub for modulation of other brain regions.^51,52 Given this, by understanding how different regions within the area relate to one another and to the rest of the brain will shed light on this complex process. That the middle and posterior STG show differential connectivity patterns allows for the generation of neural control hypotheses that can drive our future work in a large group of subjects and help provide mechanistic information about vocalization in healthy and disease.  For example, we know that patients with Parkinson’s disease have hypophonia in which they vocalize with very low loudness levels. Importantly, we know that they under scale motor output such that they report that they are speaking at normal or greater than normal loudness levels. These patients have abnormal responses to altered auditory feedback and this may be associated with abnormal STG function within the broader network of brain regions controlling vocalization.

Figure 1 Listen and Vocalization Conditions – Right Hemisphere (Subject 186).

In our fMRI and evoked response potential (ERP) studies we have shown key connectivity differences associated with alterations in auditory feedback during vocalization.  Our dynamic causal modeling of ERP signals has shown important coupling properties associated with bilateral STG, inferior frontal and pre-motor cortices. We have extended this observations to study musicians with absolute pitch, those with relative pitch and non-musicians showing that connectivity variations in these regions differentiated the three subject groups. While these findings are important, the ability to use ECoG data and capitalize on its temporal and spatial resolution will allow for far greater insight into the role of different regions of STG in vocal control.

Limitations: The method introduced here includes a novel and powerful approach to modeling ECoG data. However, we acknowledge several limitations that warrant investigation in future studies. First, in the present study the selection of electrodes was anatomically driven. In a future study we plan to investigate the utility of our modeling approach using a physiological response-based electrode location strategy. Second, a larger sample of subjects will allow us to conduct a more comprehensive evaluation of our modeling framework in measuring neuronal dynamics focusing on intra-individual (single-subject) versus inter individual (group-level). Third, in a related manner, using a larger sample of subjects exhibiting language dominant versus non-dominant network differences, we will be able to evaluate the utility of our modeling approach in accurately capturing brain dynamics. 

None.

Author declares that there are no conflicts of interest.

1. Friston KJ.  Functional Integration. In: K Friston, Ashburner SJ, Editors. Statistical Parametric Mapping: The Analysis of Functional Brain Images. London, UK. Academic Press; 2007.
2. Zhang T, Wu J, Li F, et al. A dynamic directional model for effective brain connectivity using electrocorticographic (ECoG) time series. Journal of the American Statistical Association. 2015;110(509):93‒106.
3. Freestone DR, Karoly PJ, Dragan N, et al. Estimation of effective connectivity via data-driven neural modeling. Frontiers in Neuroscience. 2014;8:383.
4. Bénar CG, Grova C, Kobayashi E, et al. EEG-fMRI of epileptic spikes: Concordance with EEG source localization and intracranial EEG. NeuroImage. 2006;30(4):1161‒1170.
5. Greenlee JDW, Jackson AW, Chen F, et al.  Human Auditory Cortical Activation during Self-Vocalization. PLoS One. 2011;6(3):e14744.
6. Greenlee JDW, Behroozmand R, Larson CR, et al. Sensory-Motor Interactions for Vocal Pitch Monitoring in Non-Primary Human Auditory Cortex. PLoS One. 2013;8(4):e60783.
7. McIntosh AR, Gonzales Lima F. Structural Equation Modeling and its Application to Network Analysis in Functional Brain Imaging. Human Brain Mapping. 1994;2(1-2):2‒22.
8. Price LR. Analysis of Imaging Data. In: Little TD, Editor. Oxford Handbook of Quantitative Methods. UK. Oxford University Press; 1994(2):175‒197.
9. Friston KJ. Functional and Effective Connectivity in Neuroimaging: A Synthesis. Human Brain Mapping. 1994;2:56‒78.
10. Jöreskog K, Sörbom D. Advances in Factor Analysis and Structural Equation Models. Massachusetts, USA. Abt Books;1979.
11. Bentler IM, Weeks DG. Linear Structural Equations with Latent Variables. Psychometrika. 1980;45(3):289‒308.
12. Lee SY. Structural Equation Modeling: A Bayesian Approach. New York, USA. Wiley Series in Probability & Statistics; 2007.
13. Price LR. Small sample properties of Bayesian multivariate autoregressive time series models. Structural Equation Modeling. 2012;19(1):51‒64.
14. Price LR, Peter TF, Roger JI, et al. Modeling dynamic functional neuroimaging data using structural equation modeling. Structural Equation Modeling. 2009;16(1):147‒162.
15. Pearl J. Causality: Models, reasoning and inference. Cambridge, New York, USA. University Press; 2000.
16. Dayan P, Abbott LF. Theoretical neuroscience. Computational and mathematical modeling of neural systems. MIT Press. 2001;36:3.
17. Dayan P, Hinton GE, Neal NM, et al. The Helmholtz machine. Neural Computation. 1995;7(5):889‒904.
18. Lee SY, Song XY. Evaluation of Bayesian and maximum likelihood approaches in analyzing structural equation models with small sample sizes. Multivariate Behavioral Research. 2004;39(4):653‒686.
19. Box GEP, Jenkins GM, Reinsel GC. Time series analysis: Forecasting and control. 4th edn. 2008:784.
20. Lütkepohl H. New introduction to multiple time series analysis. New York, USA. Springer; 2006.
21. West M, Harrison J. Bayesian forecasting and dynamic control.2nd edn. NY: New York, USA. Springer-Verlag; 1997.
22. Soares TM, Gonclaves FB, Gamerman D.  An integrated Bayesian model for DIF analysis. Journal of Educational and Behavioral Statistics. 2009;34(3):348‒377.
23. Dunson DB. Bayesian latent variable models for clustered mixed outcomes. Journal of the Royal Statistical Society Series B. 2000;62(2):355‒366.
24. Scheines R, Hoijtink H, Boomsma A. Bayesian estimation and testing of structural equation models. Psychometrika. 1999;64(1):37‒52.
25. Cohen J, Cohen P, West SG, et al. Applied multiple regression/correlation analysis for the behavioral sciences. Journal of Educational and Behavioral Statistic. 2003;30(2):227‒229.
26. Kim J, Zhu W, Chang L, et al. Unified Structural Equation Modeling Approach for the Analysis of Multisubject, Multivariate Functional MRI data. Human Brain Mapping. 2007;28(2):85‒93.
27. Gates KM, Molenaar PC, Hillary FG, et al. Extended unified SEM approach for modeling event-related fMRI data. Neuroimage. 2011;54(2):1151‒1158.
28. Ramsey JO, Hooker G, Campbell D, et al. Parameter estimation for differential equations: a generalized smoothing approach. Journal of the Royal Statistical Society Series B. 2007;69(5):741‒796.
29. Fleiss M, Lamnabhi M, Lamnabhi-Lagarrigue F. An algebraic approach to nonlinear functional expansions. IEEE Trans Circuits Syst. 1983;30(8):554‒579.
30. Ramsay JO, Silverman BW. Applied functional data analysis. New York, USA. Springer Series in Statistics; 2005.
31. Friston KJ, Harrison L, Penny W. Dynamic causal modelling. NeuroImage. 2003:19(4):1273‒1302.
32. Molenaar PCM, Sinclair KO, Rovine MJ, et al. Analyzing developmental processes on an individual level using nonstationary time series modeling. Developmental Psychology. 2009;45(1): 260‒271.
33. Nesselroade JR, Molenaar PCM. Pooling Lagged Covariance Structures Based on Short, Multivariate Time Series for Dynamic Factor Analysis. In: Hoyle, Editor. Statistical Strategies for Small Sample Research. Thousand Oaks, CA: Newbury Park, California. Sage Publications; 1999:223‒251.
34. An HZ, Chen SG. A note on the ergodicity of non-linear autoregressive models. Statistics & Probability Letters. 1997;34(4):365‒372.
35. Salhi S. Heuristic Search Methods. In: Marcoulides GA, Editor. Modern Methods for Business Research, Quantitative Methods Series. New Jersey, USA. Lawrence Erlbaum Associates; 2009:147‒176.
36. Hastie T, Tibshirani R, Freidman. The elements of statistical learning: Data mining, inference and prediction. 2nd edn. Springer Series in Statistics. New York, USA. Springer Nature; 2009:1‒745.
37. Burnham KP, Anderson DR. Model selection and multimodal inference: A practical information-theoretic approach (2nd edn) Statistical Theory and Methods. New York, USA. Springer Nature; 2002:1‒485.
38. Kullback S, Leibler RA. On information and sufficiency. Annals of Mathematical Statistics. 1951;22(1):79-86.
39. Raftery AE. Bayesian model selection in structural equation models. In: Bollen KA,  Long JS, Ediotors. Newbury Park, California. Sage Publications;1993:163‒180.
40. Browne MW, Cudeck R. Single-sample cross-validation indices for covariance structures. Multivariate Behavioral Research. 1989;24(4):445‒455.
41. Gelman A, Carlin JB, Stern HS, et al. Bayesian data analysis.2nd edn. Boca Raton FL, USA. Chapman & Hall; 1990..
42. Steiger JH. Structural model evaluation and modification: An interval estimation approach. Multivariate Behav Res. 1990;25(2):173‒180.
43. Hoff, PD. A first course in Bayesian statistical methods. New York, NY, USA. Springer; 2009.
44. Sevinç V, Ergün G. Usage of different prior distributions in Bayesian vector autoregressive models. Hacettepe Journal of Mathematics and Statistics. 2009;38(1):85‒93.
45. Gelman A, Jakulin A, Pittau MG, et al. A weakly informative default prior distribution for logistic and other regression models. The Annals of Applied Statistics. 2008;2(4):1360‒1383.
46. Anderson TW. An introduction to multivariate statistical analysis. 3rd edn. New York, USA. Wiley; 2003.
47. http://www.statmodel.com.
48. Muthén BO,  Muthén L. Mplus computer program version 7.3. Mplus, Los Angeles, California. 2012:1‒711.
49. Brooks S, Gelman A, Jones G, et al. Handbook of Markov Chain Monte Carlo. Chapman & Hall/CRC. Taylor & Francis, Boca Raton, Florida. 2011:1‒575.
50. Cudeck R, Kleve KJ, Henly SJ.  A Simple Gauss-Newton Procedure for Covariance Structure Analysis with High-Level Computer Languages. Psychometrika. 1993;58(2): 211‒232.
51. Laird AR, Robbins JM, Li K, et al. Modeling Motor Connectivity using TMS/PET and Structural Equation Modeling. Neuroimage. 2008;41(2):424‒436.
52. Parkinson AL,  Korzyukov O. Modulation of effective connectivity during vocalization with perturbed auditory feedback. Neuropsychologia. 2013;51(8):1471‒1480.