Switching Models of Oscillatory Networks Greatly Improve Inference of Dynamic
Functional Connectivity

Wan-Chi Hsin1, Uri T. Eden1 and Emily P. Stephen1
1Department of Mathematics and Statistics, Boston University, Boston, MA
02215, U.S.A.

Keywords: State-space models, functional networks, functional connectivity,
switching Kalman filters

Abstract

Functional brain networks can change rapidly as a function of stimuli or
cognitive shifts. Tracking dynamic functional connectivity is particularly
challenging as it requires estimating the structure of the network at each
moment as well as how it is shifting through time. In this paper, we describe
a general modeling framework and a set of specific models that provides
substantially increased statistical power for estimating rhythmic dynamic
networks, based on the assumption that for a particular experiment or task,
the network state at any moment is chosen from a discrete set of possible
network modes. Each model is comprised of three components: (1) a set of
latent switching states that represent transitions between the expression of
each network mode; (2) a set of latent oscillators, each characterized by an
estimated mean oscillation frequency and an instantaneous phase and amplitude
at each time point; and (3) an observation model that relates the observed
activity at each electrode to a linear combination of the latent oscillators.
We develop an expectation-maximization procedure to estimate the network
structure for each switching state and the probability of each state being
expressed at each moment. We conduct a set of simulation studies to illustrate
the application of these models and quantify their statistical power, even in
the face of model misspecification.

##  1 Introduction

Neural oscillations are critical for coordinating information across neurons
and brain areas (Singer, 1999; Fries et al., 2001; Fries, 2005; Buzsáki et
al., 2012). When groups of neurons engage in rhythmic activity, they can
become functionally connected across large brain distances. This phenomenon
can be viewed through the lens of coupled oscillators, where dynamics of
neural synchronization reflect changes in information processing and cognitive
function (Horwitz, 2003; Rogers et al., 2007). Changes in brain-wide
oscillatory coupling networks can be used to track changes in healthy and
diseased brain states and enhance our understanding of the role played by
neuronal oscillations (Schnitzler and Gross, 2005; Başar, 2013; Cocchi et al.,
2017; Wang et al., 2018). To analyze complex rhythmic interactions in the
brain, researchers often rely on descriptive statistics, in particular,
spectral estimators, such as coherence and cross-spectral estimates (Buzsáki
and Draguhn, 2004; Fries, 2005). More advanced statistical modeling methods,
such as global coherence analysis, have been developed to understand complex
brain activity and functional networks (Mitra, 2007; Cimenser et al., 2011;
Wong et al., 2011). However, these descriptive methodologies are difficult to
incorporate into statistical models of brain dynamics and are typically
restricted to single narrow-band rhythms. More importantly, these approaches
are limited in statistical efficiency and temporal resolution in estimating
dynamic functional connectivity.

Recently, state-space models have been successful in characterizing rhythmic
neural activity in the time-domain (Matsuda and Komaki, 2017; Beck et al.,
2018; Soulat et al., 2022). State-space models provide a powerful analytic
framework across diverse disciplines including control engineering, signal
processing, and neuroscience. They are typically characterized by two key
equations: the state equation, which describes how the latent state of the
system evolves over time, and the observation equation, which relates the
observed data to this latent state (Chen and Brown, 2013). State-space models
are particularly effective in analyzing neural signal dynamics because of
their generality and ease of interpretation (Kass et al., 2014), and have been
widely used in neuroscience to analyze neural dynamics (Smith and Brown, 2003;
Eden et al., 2004; Chen et al., 2010; Paninski et al., 2010; Akram et al.,
2016; Eden et al., 2018; Wodeyar et al., 2021). Many of these analyses assume
that the statistical properties of the models are time-invariant, but neural
activity can shift rapidly between different states such as sleep and
wakefulness (Steriade et al., 1993). In this article, we introduce a switching
state-space modeling framework that consists of multiple linear state-space
models and transitions between these models through time (Murphy, 1998).

The switching Kalman filter, with its capability to capture dynamic
transitions by having a bank of different linear models and switching between
them, has been applied in various fields. In (Veeraraghavan et al., 2005), it
was applied to track vehicles and detect events at traffic intersections, with
the goal of improving traffic management systems. In meteorology, it was used
to monitor and forecast storm locations in meteorological networks by
dynamically adjusting to changing weather conditions (Manfredi et al., 2005).
It has also been used to predict potential failures for mechanical systems and
proven to be a promising tool to support maintenance decision-making in
aerospace, manufacturing, and transportation (Lim and Mba, 2015).

Our work is motivated by the challenges associated with conventional windowed
methods of estimating dynamic functional connectivity, which have limited
statistical power. Time-frequency analyses using sliding windows have been
widely applied by the neuroimaging and electrophysiology community to
understand how brain dynamics relate to our cognitive behavior (Kucyi et al.,
2013; Kucyi and Davis, 2014; Elton and Gao, 2015; O’Neill et al., 2015; Chang
and Glover, 2010; Preti et al., 2017). Consider a scenario where multiple
neural sources are described by coupled oscillators. Researchers often compute
coherograms, which measure the coherence between pairs of signals as a
function of time and frequency, to assess the degree of linear correlation
between these signal pairs in the frequency domain. However, the statistical
power of coherograms is limited because it relies on relatively small temporal
windows for analysis. This constraint impedes our ability to identify
connections with confidence over time and, consequently, prevents us from
accurately inferring the dynamics of the network structure.

Rhythmic networks in brain activity may reflect several different mechanisms
of coupling: for example, networks may be driven by a common source
influencing rhythms in multiple areas, or rhythms may propagate through local
networks either by driving each other or through a shared source of noise. In
any case, we posit that specific cognitive states and behavioral tasks may be
mediated by a small number of network modes whose expression changes
dynamically. For instance, during general anesthesia, as a subject transitions
from a conscious to unconscious state and back, this change in state is
reflected in stereotypical patterns of oscillatory coherence across brain
areas (Cimenser et al., 2011; Purdon et al., 2013). In our previous work (Hsin
et al., 2022), we modeled the first mechanism using a “Common Oscillator Model
(COM),” in which the activity across multiple electrodes are linked through a
small number of latent oscillators, with each oscillator describing the
associated activity of a single network mode. We used an EM algorithm to
estimate the composition of each network mode, with a switching Kalman filter
and smoother to estimate the instantaneous state of the oscillations driving
each network mode. In contrast to windowed methods, such as the coherogram,
which makes no assumptions about which links are present at any given time,
our approach assumes a discrete set of networks and by switching between them,
greatly improves the inference of dynamic functional connectivity.

In this work, we first introduce a broad framework for analyzing dynamically
coupled oscillators in the brain. In particular, we describe three specific
modeling structures that reflect common types of coupling. The fundamental
problem we aim to address is the changing functional connectivity among
observations, characterized by a small number of distinct functional networks.
This modeling framework enables us to estimate both the structures of these
functional networks and the times at which switches occur. In addition to the
Common Oscillator Model, we develop two other models to capture propagating
rhythms. The first we call the Correlated Noise Model (CNM), which defines
networks linked by correlation structure in the noise driving the oscillations
on each electrode. The second we call the Directed Influence Model (DIM),
which allows the instantaneous state of one oscillating source to directly
influence another. All these methods use the full data to estimate a discrete
set of functional network modes and use local data over short time scales to
estimate which network modes are being expressed at each moment. Each model is
comprised of three components: (1) a set of latent switching states that
represent transitions between the expression of each network mode; (2) a set
of latent oscillators, each characterized by an estimated mean oscillation
frequency and an instantaneous phase and amplitude at each time point; and (3)
an observation model that relates the observed activity at each electrode to a
linear combination of the latent oscillators. We use switching Kalman filters
and smoothers to estimate the instantaneous phase and amplitude of each
oscillator and the probability of each switching state at every time step. The
estimated switching state characterizes which of the network modes is being
expressed at each time point.

##  2 Methods

We begin by introducing a broad analytic framework that uses latent processes
to represent the oscillators that drive brain dynamics, and subsequently
derive three specific model structures to show different ways that network
structure can be captured. In our approach, we extend the oscillation
decomposition framework described in (Matsuda and Komaki, 2017) and use state
models representing sets of oscillators. Each oscillator comprises a pair of
latent variables with coupled dynamics. For example, we model a single
oscillator as a random vector process
$x_{t}=[x_{t,\text{Re}},x_{t,\text{Im}}]^{\prime}$, with one-step update
equation $x_{t}=Ax_{t-1}+\epsilon_{t}$, where $A$ is a $2\times 2$ rotation
matrix of the form
$A=a\begin{pmatrix}\text{cos}(\omega)&-\text{sin}(\omega)\\\
\text{sin}(\omega)&\text{cos}(\omega)\end{pmatrix}$, $a$ is an autoregressive
parameter, and $\omega$ is the rotation frequency. Note that we write
$x_{t,\text{Re}}$ and $x_{t,\text{Im}}$ to evoke the real and imaginary parts
of the analytic representation of a narrow-band signal, $z(t)=x(t)+ix_{H}(t)$,
where $H$ represents the Hilbert transform (See (Matsuda and Komaki, 2017) for
more details), even though all random variables are real valued.

More generally, let $K$ denote the number of oscillators. The $2K\times 1$
vector
$x_{t}=[x^{1}_{t,\text{Re}},x^{1}_{t,\text{Im}}\ldots,x^{K}_{t,\text{Re}},x^{K}_{t,\text{Im}}]^{\prime}$
is defined as the hidden state variable with $K$ oscillators at time $t$,
where $[x^{k}_{t,\text{Re}},x^{k}_{t,\text{Im}}]^{\prime}$ is the state of
oscillator $k$ at time $t$.

Next, we define a discrete switching state $S_{t}$ which takes on one of a
small number of discrete states at each time step. The switching state
determines the dynamics of the oscillatory state $x_{t}$ at each time step and
the influence of the state on the observed data. Importantly, $S_{t}$ will
control which functional network is expressed at time $t$. We assume that
$S_{t}$ is a Markov chain, and denote its transition matrix $Z$, whose
$(i,j)^{th}$ element defines the probability of a switch between time steps
$t-1$ and $t$.

$P(S_{t}=j|S_{t-1}=i)=Z_{ij}.$ (2.1)

The oscillatory latent process model is defined so as to depend on the current
value of the switching state:

$(x_{t}|x_{t-1},S_{t}=j)\sim N(A_{j}x_{t-1},\Sigma_{j}),$ (2.2)

where $A_{j}$ is a $2K\times 2K$ state transition matrix defining the dynamics
of all $K$ oscillators in the current switching state, and $\Sigma_{j}$ is a
$2K\times 2K$ process noise covariance matrix.

The data recorded at each of $N$ neural sources is modeled as an
$N$-dimensional multivariate time series,
$y_{t}=[y_{t}^{1},\ldots,y_{t}^{N}]^{\prime}$. The observation equation is
also defined so as to depend on the current value of the switching state:

$(y_{t}|x_{t},S_{t}=j)\sim N(B_{j}x_{t},R),$ (2.3)

where $B_{j}$ is an $N\times 2K$ observation matrix defining the influence of
each of the two components of each of the $K$ oscillators on each of the $N$
neural sources, and $R$ is an $N\times N$ observation noise covariance matrix.

There are multiple ways that the $N$ neural sources can be coupled in this
general model framework. In particular, options include: (1) using the
observation matrices, $B_{j}$, to express how the activity of a small number
of latent oscillators are reflected across many neural sources, (2) using the
state covariance matrices, $\Sigma_{j}$, to express shared noise structure
propagating through the network, or (3) using the state update matrices,
$A_{j}$, to express the direct influence of one source on each of the others.
In the following three subsections, we derive specific models to capture
functional networks based on these mechanisms.

###  2.1 Common Oscillator Model

For this model, we posit that neural sources are linked by being driven by a
common set of oscillators. The network structure is reflected in the
observation matrix, which determines which of the oscillators is driving each
node and the relation between the instantaneous phase and amplitudes of the
oscillators and each observed signal. A schematic for this model is shown in
Figure 1, where the red arrows going from switching state variable $S_{t}$ to
the observations $y_{t}$ indicate that the observation matrix $B_{j}$ changes
over time depending on the switching state. Note that in this model, the state
transition matrix $A$ and the state noise covariance $\Sigma$ do not depend on
the switching state $S_{t}$, so are assumed to be constant in time. When $B$
switches, the functional network driven by each oscillator changes as well.
This model structure is particularly efficient when there are a small number
of latent oscillators $K$ relative to the number of sources $N$, so we show
the arrow from $x_{t}$ to $y_{t}$ splitting. This schematic displays an
example of one oscillator driving a pair of observed signals.

Figure 1: Schematic representation of the Common Oscillator Model. The
observations $y_{t}$ are a linear function of the continuous latent
oscillators $x_{t}$, and the observation matrix $B$ depends on the discrete
latent state $S_{t}$. This schematic shows an example with one oscillator
driving two neural sources.

The latent process model is defined as

$(x_{t}|x_{t-1})\sim N(Ax_{t-1},\Sigma),$ (2.4)

where

$A=\begin{pmatrix}A_{1}&0&\cdots&0\\\ 0&A_{2}&\cdots&0\\\
\vdots&\vdots&\ddots&\vdots\\\ 0&0&\cdots&A_{K}\end{pmatrix},$

and

$\Sigma=\begin{pmatrix}\sigma_{1}I_{2\times 2}&0&\cdots&0\\\
0&\sigma_{2}I_{2\times 2}&\cdots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&0&\cdots&\sigma_{K}I_{2\times 2}\end{pmatrix}.$

$A$ is a $2K\times 2K$ block diagonal matrix composed of $A_{1},\ldots,A_{K}$.
For $k=1,\ldots,K$,

$A_{k}=a_{k}\begin{pmatrix}\text{cos}(\omega_{k})&-\text{sin}(\omega_{k})\\\
\text{sin}(\omega_{k})&\text{cos}(\omega_{k})\end{pmatrix},$

and $\omega_{k}=2\pi f_{k}/F_{s},$ where $f_{k}$ is the oscillation frequency
and $F_{s}$ is the sampling frequency. Each $A_{k}$ is a 2-dimensional
rotation matrix multiplied by an autoregressive parameter $a_{k}\in(0,1)$.
$\Sigma$ is a $2K\times 2K$ diagonal noise covariance matrix.

For an $N$-dimensional multivariate time series sampled at frequency $F_{s}$,
the observation at time $t$ is $y_{t}=[y_{t}^{1},\ldots,y_{t}^{N}]^{\prime}$,
where $y_{t}^{n}$ is the data recorded at node $n$ at time $t$. The
observation equation is given by

$(y_{t}|x_{t},S_{t}=j)\sim N(B_{j}x_{t},R).$ (2.5)

The observation process $y_{t}$ is the weighted sum of contributions from one
or more oscillatory latent states. These contributions are determined by the
$N\times 2K$ observation matrix $B_{j}$, where the subscript $j$ indicates the
matrix is associated with switching state $S_{t}=j$. Each switching state has
its own observation matrix, and each of these matrices represents a distinct
network mode. The observation noise is Gaussian with $N\times N$ covariance
matrix $R$. We can write the observation matrix for state $j$ elementwise:

$B_{j}=\begin{pmatrix}b_{j,\text{Re}}^{1,1}&b_{j,\text{Im}}^{1,1}&\cdots&b_{j,\text{Re}}^{1,K}&b_{j,\text{Im}}^{1,K}\\\
\vdots&\vdots&\ddots&\vdots&\vdots\\\
b_{j,\text{Re}}^{N,1}&b_{j,\text{Im}}^{N,1}&\cdots&b_{j,\text{Re}}^{N,K}&b_{j,\text{Im}}^{N,K}\end{pmatrix}$
(2.6)

The elements $b_{j,\text{Re}}^{n,k}$ and $b_{j,\text{Im}}^{n,k}$ represent the
respective influence of the first and second components of the $k^{\text{th}}$
oscillator on the $n^{\text{th}}$ node. The amplitude
$\sqrt{(b_{j,\text{Re}}^{n,k})^{2}+(b_{j,\text{Im}}^{n,k})^{2}}$ describes the
degree to which the oscillator $k$ that is expressed at node $n$ in state $j$.
For any two nodes, $n_{1}$ and $n_{2}$ driven by the same oscillator $k$, the
phase difference between them is equal to
$\theta_{n_{1},n_{2}}=\theta_{n_{1}}-\theta_{n_{2}}$, where $\theta_{n_{i}}$
is the argument of the complex number
$b_{j,\text{Re}}^{n_{i},k}+ib_{j,\text{Im}}^{n_{i},k}$.

By comparing the state equation 2.4 and the observation model 2.5 of the
Common Oscillator Model with those from the global model, equations 2.2 and
2.3, the key difference is that we suppress the the subscript $j$ in the state
transition matrix $A$ and the process noise covariance $\Sigma$. In the COM,
we use the observation model to express the network structure. This is
represented by the observation matrix $B_{j}$, which changes over time
depending on the switching state variable $S_{t}$.

###  2.2 Correlated Noise Model

In contrast to the COM, the Correlated Noise Model (CNM) assumes that each
node is characterized by its own oscillator, any combination of which can be
linked through correlated noise structure. The network structure is then
reflected in the process noise covariance matrix. In Figure 2, the red arrows
from the switching state variable $S_{t}$ to the latent oscillatory state
$x_{t}$ indicate that the process noise covariance $\Sigma_{j}$ changes
through time, depending on the switching state $S_{t}$. In this model, the
state transition matrix $A$ and the observation matrix $B$ do not depend on
the switching state $S_{t}$. Note that the number of layers in $y_{t}$ is same
as the number of layers in $x_{t}$, since the number of observations is equal
to the number of oscillators, with one oscillator per node.

Figure 2: Schematic representation of the Correlated Noise Model. The latent
oscillators $x_{t}$ is directly influenced by the discrete latent state
$S_{t}$ through the process noise covariance matrix $\Sigma$.

The latent process model is defined as

$(x_{t}|x_{t-1},S_{t}=j)\sim N(Ax_{t-1},\Sigma_{j}),$ (2.7)

where the $A$ matrix shares the same structure as the $A$ matrix described in
the Common Oscillator Model in Section 2.1.

The process noise covariance matrix is

$\begin{split}&\Sigma_{j}=\begin{pmatrix}\Sigma_{j}^{1}&\Sigma_{j}^{1,2}&\cdots&\Sigma_{j}^{1,K}\\\
\Sigma_{j}^{2,1}&\Sigma_{j}^{2}&\cdots&\Sigma_{j}^{2,K}\\\
\vdots&\vdots&\ddots&\vdots\\\
\Sigma_{j}^{K,1}&\Sigma_{j}^{K,2}&\cdots&\Sigma_{j}^{K}\end{pmatrix},\text{where}\\\
\\\ &\Sigma_{j}^{k}=\sigma_{j}^{k}I_{2\times 2},\text{and}\\\
&\Sigma_{j}^{n_{1},n_{2}}=\left\\{\begin{array}[]{ll}\rho_{j}^{n_{1},n_{2}}\begin{pmatrix}\mbox{cos}(\theta_{j}^{n_{1},n_{2}})&-\mbox{sin}(\theta_{j}^{n_{1},n_{2}})\\\
\mbox{sin}(\theta_{j}^{n_{1},n_{2}})&\mbox{cos}(\theta_{j}^{n_{1},n_{2}})&\end{pmatrix},&\text{if
$n_{1}$ and $n_{2}$ are linked.}\\\ 0_{2\times
2},&\text{otherwise.}\end{array}\right.\end{split}$ (2.8)

The amplitude $\rho_{j}^{n_{1},n_{2}}$ describes the degree to which the nodes
$n_{1}$ and $n_{2}$ are linked and $\theta_{j}^{n_{1},n_{2}}$ is the phase
difference between these two nodes in the switching state $j$. The rotational
structure in the off-diagonal blocks indicates that when one node receives a
random perturbation, the other node to which it is linked receives a related
random perturbation in a phase-dependent manner. Note that the process noise
covariance $\Sigma_{j}$ is symmetric. For example, if $n_{1}$ and $n_{2}$ are
linked with a 90-degree phase difference in state $j$
($\theta_{j}^{n_{1},n_{2}}=\pi/2$), then $n_{2}$ and $n_{1}$ are linked with a
270-degree (-90-degree) phase difference ($\theta_{j}^{n_{2},n_{1}}=-\pi/2$).

Under this model, the observation equation is given by

$(y_{t}|x_{t})\sim N(Bx_{t},R),$ (2.9)

where $B$ is a $K\times 2K$ block diagonal matrix

$B=\begin{pmatrix}B_{1}&0&\cdots&0\\\ 0&B_{2}&\cdots&0\\\
\vdots&\vdots&\ddots&\vdots\\\ 0&0&\cdots&B_{K}\end{pmatrix}$

and $B_{k}=[b^{k}_{\text{Re}},b^{k}_{\text{Im}}]$ determines the extent to
which the real and imaginary components of the oscillator are expressed in the
data. For our simulations we opted to set $b^{k}_{\text{Re}}=1$ and
$b^{k}_{\text{Im}}=0$ so that the signal at each node is determined by only
the real component of the oscillator.

The main difference between the global model descried in equations 2.2 and 2.3
and the specific CNM structure, is the suppression the subscript $j$ in the
state transition matrix $A$ and the observation matrix $B$. In the CNM, the
state noise covariance $\Sigma_{j}$ determines the network structure, and it
changes over time depending on the switching state variable $S_{t}$.

###  2.3 Directed Influence Model

The last model structure we develop is the Directed Influence Model (DIM), in
which each neural source is characterized by its own oscillator, and the state
of one oscillator directly influences another. The network structure is then
determined by the state transition matrix. In Figure 3, the red arrows going
from the switching state variable $S_{t}$ to the latent oscillatory state
$x_{t}$ illustrates that the switching state variable determines which state
transition matrix $A_{j}$ is used at any given time. In this model, the
observation noise covariance $\Sigma$ and the observation matrix $B$ do not
depend on the switching state $S_{t}$. Note that, similar to the CNM, the
number of oscillators is same as the number of neural sources.

Figure 3: Schematic representation of the Directed Influence Model. The latent
oscillators $x_{t}$ are directly influenced by the discrete latent state
$S_{t}$ through the state transition matrix $A$.

The latent process model is defined as

$(x_{t}|x_{t-1},S_{t}=j)\sim N(A_{j}x_{t-1},\Sigma),$ (2.10)

where the process noise covariance matrix $\Sigma$ has the same structure as
in the Common Oscillator Model in Section 2.1. The state transition matrix
$A_{j}$ is defined as

$\begin{split}&A_{j}=\begin{pmatrix}A_{j}^{1}&A_{j}^{1,2}&\cdots&A_{j}^{1,K}\\\
A_{j}^{2,1}&A_{j}^{2}&\cdots&A_{j}^{2,K}\\\ \vdots&\vdots&\ddots&\vdots\\\
A_{j}^{K,1}&A_{j}^{K,2}&\cdots&A_{j}^{K}\end{pmatrix},\text{where}\\\ \\\
&A_{j}^{k}=\alpha_{j}^{k}\begin{pmatrix}\mbox{cos}(2\pi
f_{n}/F_{s})&-\mbox{sin}(2\pi f_{n}/F_{s})\\\ \mbox{sin}(2\pi
f_{n}/F_{s})&\mbox{cos}(2\pi
f_{n}/F_{s})&\end{pmatrix}-\displaystyle\sum_{n_{1}=k}\alpha_{j}^{n_{1},n_{2}}\begin{pmatrix}1&0\\\
0&1\end{pmatrix},\text{and}\\\
&A_{j}^{n_{1},n_{2}}=\left\\{\begin{array}[]{ll}\alpha_{j}^{n_{1},n_{2}}\begin{pmatrix}\mbox{cos}(\phi_{j}^{n_{1},n_{2}})&-\mbox{sin}(\phi_{j}^{n_{1},n_{2}})\\\
\mbox{sin}(\phi_{j}^{n_{1},n_{2}})&\mbox{cos}(\phi_{j}^{n_{1},n_{2}})\end{pmatrix},&\text{if
$n_{1}$ and $n_{2}$ are linked.}\\\ 0_{2\times
2},&\text{otherwise.}\end{array}\right.\end{split}$ (2.11)

The amplitude $\alpha_{j}^{n_{1},n_{2}}$ represents the coupling strength
between nodes $n_{1}$ and $n_{2}$, and $\phi_{j}^{n_{1},n_{2}}$ is the phase
difference between these two nodes in switching state $j$. The diagonal blocks
include one term that is a rotation matrix, representing the oscillator for
the corresponding node, and a second adjustment term that is a row sum of
coupling strengths present on the off-diagonal blocks, multiplied by an
identity matrix. This adjustment term comes from the influences received from
other nodes. Note that the state transition matrix $A_{j}$ is not necessarily
symmetric.

As with the CNM, under this model, the observation equation is given by

$(y_{t}|x_{t})\sim N(Bx_{t},R),$ (2.12)

where $B$ is a $K\times 2K$ block diagonal matrix

$B=\begin{pmatrix}B_{1}&0&\cdots&0\\\ 0&B_{2}&\cdots&0\\\
\vdots&\vdots&\ddots&\vdots\\\ 0&0&\cdots&B_{K}\end{pmatrix}$

and $B_{k}=[b^{k}_{\text{Re}},b^{k}_{\text{Im}}]$ determines the extent to
which the real and imaginary components of the oscillator are expressed at
each node. For the sake of interpretability, we opted to set
$b^{k}_{\text{Re}}=\frac{1}{\sqrt{2}}$ and
$b^{k}_{\text{Im}}=\frac{1}{\sqrt{2}}$ so that the signal at each node is
determined equally by the real and imaginary components of the oscillator. If
we were to use an observation matrix that takes only the real parts, it would
imply that the imaginary component of each oscillator remains unconstrained
across oscillators, allowing for a mixing of oscillators when we estimate the
A matrix.

In the DIM, we remove the subscript $j$ from the process noise covariance
$\Sigma$ and the observation matrix $B$ in the global model within equations
2.2 and 2.3. Here, the state transition matrix $A_{j}$ describes how different
latent sources influence each other in a directed manner, leading to directed
functional connectivity between these nodes. We note that this is the only one
of these three models that captures directional influences.

###  2.4 Oscillatory and Switching State Estimation

We aim to estimate the latent oscillatory state $x_{t}$ and the switching
state $S_{t}$ given the observations $y_{t}$. However, the posterior
distribution of $x_{t}$ is a mixture of Gaussians where the number of mixture
components grows exponentially in time so that there would be $M^{t}$
components by time step $t$. To deal with this exponential growth, we utilize
a collapse approach, which uses the Generalized Pseudo Bayesian algorithm
(Bar-Shalom et al., 2001; Li and Zhang, 2000) to approximate the mixture of
$M^{t}$ Gaussians using only $M$ components at each time step $t$.

We compute the posterior smoothing distributions $P(x_{t}|y_{1:T})$ and
$P(S_{t}=j|y_{1:T})$, of the latent states up to time $t$, given all of the
observations up to time $T$ in two steps. We first apply the Kalman filter
algorithm (Kalman, 1960; Sarkka, 2013) to compute the filter distributions,
$P(x_{t}|y_{1:t})$ and $P(S_{t}=j|y_{1:t})$, at each time step. The
propagation from time $t-1$ to $t$ is given by

$P(x_{t}|y_{1:t},S_{t-1}=i,S_{t}=j)=N(x_{t}^{ij},V_{t}^{ij}),$

where

$\begin{split}x_{t}^{ij}&=\text{E}(x_{t}|y_{1:t},S_{t-1}=i,S_{t}=j),\\\
V_{t}^{ij}&=\text{Cov}(x_{t}|y_{1:t},S_{t-1}=i,S_{t}=j).\end{split}$

Let $L_{t}^{ij}=P(y_{t}|y_{1:t-1},S_{t-1}=i,S_{t}=j)$ be the likelihood of
observing $y_{t}$, and denote

$\begin{split}w_{t}^{ij}&=P(S_{t-1}=i,S_{t}=j|y_{1:t}),\\\
w_{t}^{j}&=P(S_{t}=j|y_{1:t}),\\\
w_{t}^{i|j}&=P(S_{t-1}=i|S_{t}=j,y_{1:t}).\end{split}$

We then compute the following weights:

$\begin{split}w_{t}^{ij}&=\frac{L_{t}^{ij}Z_{ij}w_{t-1}^{i}}{\sum_{k}\sum_{l}L_{t}^{kl}Z_{kl}w_{t-1}^{k}},\\\
w_{t}^{j}&=\sum_{i}w_{t}^{ij},\\\
w_{t}^{i|j}&=\frac{w_{t}^{ij}}{w_{t}^{j}}.\end{split}$

The posterior filter distribution is

$P(x_{t}|y_{1:t})=\sum_{i=1}^{M}\sum_{j=1}^{M}w_{t}^{ij}P(x_{t}|y_{1:t},S_{t-1}=i,S_{t}=j),$

which is a mixture of $M^{2}$ Gaussians. To collapse the mixture of $M^{2}$
Gaussians into a mixture of $M$ Gaussians, we apply the collapse operator as
described in (Murphy, 1998) to combine the estimates from the filters in the
previous time step with the weight $w_{t}^{i|j}$. That is, for each $j$,

$P(x_{t}|y_{1:t},S_{t}=j)=\sum_{i=1}^{M}w_{t}^{i|j}P(x_{t}|y_{1:t},S_{t-1}=i,S_{t}=j),$

is a mixture of $M$ components and we approximate this as a single Gaussian by
applying the collapse operator detailed in the Appendix, so that we have

$P(x_{t}|y_{1:t},S_{t}=j)\approx N(x_{t}^{j},V_{t}^{j}),$

where

$\begin{split}x_{t}^{j}&=\text{E}(x_{t}|y_{1:t},S_{t}=j),\\\
V_{t}^{j}&=\text{Cov}(x_{t}|y_{1:t},S_{t}=j).\end{split}$

Then, the posterior probability of the hidden state conditioned on the
observations can be computed by

$P(x_{t}|y_{1:t})\approx\sum_{j=1}^{M}w_{t}^{j}P(x_{t}|y_{1:t},S_{t}=j)$

where the posterior distribution $P(x_{t}|y_{1:t})$ is approximated by a
mixture of $M$ Gaussians at each time step $t$. Note that if the goal is to
compute $P(x_{t}|y_{1:T})$, we do not have to apply the collapse operator
again for $\hat{x_{t}}$ and $V_{t}$ where
$\hat{x_{t}}=\text{E}(x_{t}|y_{1:t})$ and $V_{t}=\text{Cov}(x_{t}|y_{1:t})$,
instead, we will pass $x_{t}^{j}$ and $V_{t}^{j}$ on to the smoothing step.

In the smoothing step, we estimate $P(x_{t}|y_{1:T})$ and $P(S_{t}=j|y_{1:T})$
by using all of the data including both past and future observations. Similar
to the filter step, we have

$P(x_{t}|y_{1:T},S_{t}=j,S_{t+1}=k)=N(\tilde{x}_{t}^{jk},\tilde{V}_{t}^{jk}),$

where

$\begin{split}\tilde{x}_{t}^{jk}&=\text{E}(x_{t}|y_{1:T},S_{t}=j,S_{t+1}=k),\\\
\tilde{V}_{t}^{jk}&=\text{Cov}(x_{t}|y_{1:T},S_{t}=j,S_{t+1}=k).\end{split}$

To obtain $\tilde{x}_{t}^{jk}$ and $\tilde{V}_{t}^{jk}$, we perform the Rauch-
Tung-Strieber (RTS) smoothing algorithm (Rauch et al., 1965; Sarkka, 2013)
starting at the last time step $T$ and stepping backwards to the first time
step. We denote

$\begin{split}g_{t}^{jk}&=P(S_{t}=j,S_{t+1}=k|y_{1:T}),\\\
g_{t}^{j}&=P(S_{t}=j|y_{1:T}),\\\
g_{t}^{k|j}&=P(S_{t+1}=k|S_{t}=j,y_{1:T}).\end{split}$

We then compute the following weights:

$\begin{split}g_{t}^{jk}&\approx\frac{w_{t}^{j}Z_{jk}}{\sum_{j^{\prime}}w_{t}^{j^{\prime}}Z_{j^{\prime}k}},\\\
g_{t}^{j}&=\sum_{k}g_{t}^{jk},\\\
g_{t}^{k|j}&=\frac{g_{t}^{jk}}{g_{t}^{j}}.\end{split}$

The approximation of $g_{t}^{k|j}$ follows the standard approach for dynamic
linear models (DLMs) with Markov-switching (Kim, 1994; Murphy, 1998). We use
the collapse operator to combine the estimates from the smoothers and the
weight $g_{t}^{k|j}$ so that we have

$P(x_{t}|y_{1:T},S_{t}=j)\approx N(\tilde{x}_{t}^{j},\tilde{V}_{t}^{j}),$

where

$\begin{split}\tilde{x}_{t}^{j}&=\text{E}(x_{t}|y_{1:T},S_{t}=j),\\\
\tilde{V}_{t}^{j}&=\text{Cov}(x_{t}|y_{1:T},S_{t}=j).\end{split}$

Then, the posterior probability of the hidden oscillatory state can be
computed by

$P(x_{t}|y_{1:T})\approx\sum_{j=1}^{M}g_{t}^{j}P(x_{t}|y_{1:T},S_{t}=j),$

where $P(x_{t}|y_{1:T})$ is approximated by a mixture of $M$ Gaussians at
every time step $t$. Finally, to obtain the posterior smoothing distribution

$P(x_{t}|y_{1:T})=N(\tilde{x}_{t},\tilde{V}_{t}),$

where

$\begin{split}\tilde{x}_{t}&=\text{E}(x_{t}|y_{1:T}),\\\
\tilde{V}_{t}&=\text{Cov}(x_{t}|y_{1:T}),\end{split}$

we perform the collapse operator again to obtain a single estimate
$\tilde{x}_{t}$ and its uncertainty $\tilde{V}_{t}$.

In summary, we apply the Kalman filter and RTS smoother algorithms to estimate
the oscillatory and switching states in sequence as follows. The details of
the Filter, Smoother, and Collapse operators are outlined in the Appendix.

The forward step for filtering from time $t=1$ to $t=T$:

$\displaystyle(x_{t}^{ij},V_{t}^{ij},L_{t}^{ij})$
$\displaystyle=\text{Filter}(x_{t-1}^{i},V_{t-1}^{i},y_{t},A_{j},\Sigma_{j},B_{j},R)$
$\displaystyle w_{t}^{ij}$
$\displaystyle=\frac{L_{t}^{ij}Z_{ij}w_{t-1}^{i}}{\sum_{i}\sum_{j}L_{t}^{ij}Z_{ij}w_{t-1}^{i}},$
$\displaystyle w_{t}^{j}$ $\displaystyle=\sum_{i}w_{t}^{ij},$ $\displaystyle
w_{t}^{i|j}$ $\displaystyle=w_{t}^{ij}/w_{t}^{j},$
$\displaystyle(x_{t}^{j},V_{t}^{j})$
$\displaystyle=\text{Collapse}(x_{t}^{ij},V_{t}^{ij},w_{t}^{i|j})$

The backward step for smoothing from time $t=T$ to $t=1$:

$\displaystyle(\tilde{x}_{t}^{jk},\tilde{V}_{t}^{jk},\tilde{V}_{t+1,t}^{jk})$
$\displaystyle=\text{Smoother}(\tilde{x}_{t+1}^{k},\tilde{V}_{t+1}^{k},x_{t}^{j},V_{t}^{j},A_{j},\Sigma_{j})$
$\displaystyle g_{t}^{jk}$
$\displaystyle\approx\frac{w_{t}^{j}Z_{jk}}{\sum_{j^{\prime}}w_{t}^{j^{\prime}}Z_{j^{\prime}k}},$
$\displaystyle g_{t}^{j}$ $\displaystyle=\sum_{k}g_{t}^{jk},$ $\displaystyle
g_{t}^{k|j}$ $\displaystyle=g_{t}^{jk}/g_{t}^{j},$
$\displaystyle(\tilde{x}_{t}^{j},\tilde{V}_{t}^{j})$
$\displaystyle=\text{Collapse}(\tilde{x}_{t}^{jk},\tilde{V}_{t}^{jk},g_{t}^{k|j})$
$\displaystyle(\tilde{x}_{t},\tilde{V}_{t})$
$\displaystyle=\text{Collapse}(\tilde{x}_{t}^{j},\tilde{V}_{t}^{j},g_{t}^{j})$

To obtain the one-step covariance term,
$\tilde{V}_{t+1,t}=\text{Cov}(x_{t+1},x_{t}|y_{1:T})$, we can perform the
cross-collapse operator (details in the Appendix) twice on
$\tilde{V}_{t+1,t}^{jk}=\text{Cov}(x_{t+1},x_{t}|y_{1:T},S_{t}=j,S_{t+1}=k)$,
which is computed in the smoother step. Note that depending on which one of
these three model structures is applied, $A_{j}$, $\Sigma_{j}$, and $B_{j}$
will be fixed as $A$, $\Sigma$, and $B$, accordingly.

###  2.5 Model Parameter Estimation

We develop an expectation-maximization (EM) algorithm (Dempster et al., 1977)
to estimate simultaneously the model parameters
$\theta=\\{A_{j},\Sigma_{j},B_{j},R,Z\\}$ and the values of the latent
switching state $S_{t}$ and oscillatory state $x_{t}$ at each time $t$. The
likelihood of the complete data is given by

$P(x_{1:T},S_{1:T},y_{1:T})=\prod_{t=1}^{T}P(y_{t}|x_{t},S_{t})\prod_{t=2}^{T}P(x_{t}|x_{t-1},S_{t})\prod_{t=2}^{T}P(S_{t}|S_{t-1})P(x_{1})P(S_{1}),$

where

$\begin{split}&P(x_{t}|x_{t-1},S_{t}=j)=(2\pi)^{-\frac{1}{2}}|\Sigma_{j}|^{-\frac{1}{2}}\text{exp}\left[-\frac{1}{2}(x_{t}-A_{j}x_{t-1})^{{}^{\prime}}\Sigma_{j}^{-1}(x_{t}-A_{j}x_{t-1})\right],\\\
&P(y_{t}|x_{t},S_{t}=j)=(2\pi)^{-\frac{1}{2}}|R|^{-\frac{1}{2}}\text{exp}\left[-\frac{1}{2}(y_{t}-B_{j}x_{t})^{{}^{\prime}}R^{-1}(y_{t}-B_{j}x_{t})\right],\\\
&P(S_{t}=j|S_{t-1}=i)=Z_{ij},\\\
&P(x_{1})=(2\pi)^{-\frac{1}{2}}|C_{1}|^{-\frac{1}{2}}\text{exp}\left[-\frac{1}{2}(x_{1}-\mu_{1})^{{}^{\prime}}C_{1}^{-1}(x_{1}-\mu_{1})\right].\end{split}$
(2.13)

Note that $\mu_{1}$ and $C_{1}$ are the initial mean and covariance of the
initial latent oscillatory state $x_{1}$.

In the EM algorithm, we iteratively maximize the expected value of the
complete data log likelihood with respect to $\theta$. The E step involves
computing the expected value of the complete data log likelihood function of
$\theta$, which is also called the $Q$ function:

$\begin{split}Q(\theta)&=\text{E}[\text{log}(P(x_{1:T},S_{1:T},y_{1:T}))]\\\
&\circeq-\frac{1}{2}\sum_{t=2}^{T}\sum_{S_{t}=j}g_{t}^{j}\left[\text{E}\left[(x_{t}-A_{j}x_{t-1})^{\prime}\Sigma_{j}^{-1}(x_{t}-A_{j}x_{t-1})-\text{log}|\Sigma_{j}|\right]\right]\\\
&\quad-\frac{1}{2}\sum_{t=1}^{T}\sum_{S_{t}=j}g_{t}^{j}\left[\text{E}\left[(y_{t}-B_{j}x_{t})^{\prime}R^{-1}(y_{t}-B_{j}x_{t})\right]\right]-\frac{T}{2}\text{log}|R|,\end{split}$
(2.14)

where $g_{t}^{j}=P(S_{t}=j|y_{1:T})$. The details of the derivation process of
the E-step of the switching Kalman filter are outlined in (Murphy, 1998). The
M-step involves maximizing the $Q$ function by taking its derivative with
respect to $\theta$ and setting it to 0. The closed-form solutions for
$A_{j}$, $\Sigma_{j}$, and $B_{j}$ can be found in (Murphy, 1998), which are
essentially weighted versions of the standard closed-form solutions of the
linear Gaussian state-space models in (Ghahramani and Hinton, 1996). We
derived the analytic solution for $R$ in our previous work (Hsin et al.,
2022). The solution for the transition matrix $Z$ can be solved by using
Lagrange multipliers, which was derived in (Rabiner, 1989). Below, we
summarize the analytic solutions for updating each parameter:

$\displaystyle A_{j}$
$\displaystyle=\left(\sum_{t=2}^{T}g_{t}^{j}\hat{P}_{t,t-1}\right)\left(\sum_{t=2}^{T}g_{t}^{j}\hat{P}_{t-1}\right)^{-1}$
(2.15) $\displaystyle\Sigma_{j}$
$\displaystyle=\frac{1}{\sum_{t=2}^{T}g_{t}^{j}}\left(\sum_{t=2}^{T}g_{t}^{j}\hat{P_{t}}-A\sum_{t=2}^{T}g_{t}^{j}\hat{P}_{t,t-1}^{{}^{\prime}}\right)$
(2.16) $\displaystyle B_{j}$
$\displaystyle=\left(\sum_{t=1}^{T}g_{t}^{j}y_{t}\hat{x_{t}}^{\prime}\right)\left(\sum_{t=1}^{T}g_{t}^{j}\hat{P_{t}}\right)^{-1}$
(2.17) $\displaystyle R$
$\displaystyle=\frac{1}{\sum_{t=1}^{T}\sum_{S_{t}=j}g_{t}^{j}}\left(\sum_{t=1}^{T}\sum_{S_{t}=j}g_{t}^{j}(y_{t}y_{t}^{\prime}-B_{j}\hat{x_{t}}y_{t}^{\prime})\right)$
(2.18) $\displaystyle Z_{ij}$
$\displaystyle=\frac{1}{\sum_{t=2}^{T}g_{t}^{i}}\left(\sum_{t=2}^{T}P(S_{t-1}=i,S_{t}=j|y_{1:T})\right)$
(2.19)

where $i,j=1,\ldots,M$,

$\displaystyle\hat{x_{t}}$
$\displaystyle=\text{E}[x_{t}|y_{1:T}]=\tilde{x}_{t},$
$\displaystyle\hat{P_{t}}$
$\displaystyle=\text{E}[x_{t}x_{t}^{\prime}|y_{1:T}]=\tilde{x}_{t}\tilde{x}_{t}^{\prime}+\tilde{V}_{t},$
$\displaystyle\hat{P}_{t,t-1}$
$\displaystyle=\text{E}[x_{t}x_{t-1}^{\prime}|y_{1:T}]=\tilde{x}_{t}\tilde{x}_{t-1}^{\prime}+\tilde{V}_{t,t-1}.$

Note that $\tilde{x}_{t},\tilde{V}_{t}$, and $\tilde{V}_{t,t-1}$ are obtained
from the Smoother and Collapse operators as described in 2.4. Closed-form
solutions for the mean and covariance, $\mu_{1}$ and $C_{1}$ for the initial
oscillatory state $x_{1}$ follow the standard derivation for a Gaussian
mixture model (Hamilton, 1990; Xu and Jordan, 1996).

In Sections 2.2 and 2.3, as we introduced the structures of $\Sigma_{j}$ and
$A_{j}$; these matrices are composed of scaled rotation matrices in each
$2\times 2$ off-diagonal block. The analytic solutions for these parameters,
as shown in equations 2.15 and 2.16, do not guarantee this specific rotational
structure. We therefore added a step to impose the desired rotational
structure in $\Sigma_{j}$ and $A_{j}$. Recall that a scaled rotational matrix
is a $2\times 2$ matrix of the form $\gamma\cdot R(\psi)$ , where $\gamma$ is
the scaling factor, and $R(\psi)$ is a rotation matrix of the form:

$R(\psi)=\begin{pmatrix}\text{cos}(\psi)&-\text{sin}(\psi)\\\
\text{sin}(\psi)&\text{cos}(\psi).\end{pmatrix}$

In the M-step of the EM algorithm, we start with the analytic solution for
$\Sigma_{j}$ or $A_{j}$, depending on the model being estimated, either the
Correlated Noise Model or the Directed Influence Model. For both models, our
goal is to preserve the rotational structure. We consider the projection space
for each off-diagonal block in $\Sigma_{j}$ or $A_{j}$ to be the space of
scaled rotation matrices. Let the singular value decomposition (SVD) of
$A_{j}^{n_{1},n_{2}}$ be
$U_{j}^{n_{1}n_{2}}S_{j}^{n_{1}n_{2}}V_{j}^{n_{1}n_{2}*}$, where $*$ is the
conjugate transpose (Beaver, 2006). Then the projection to the space of scaled
rotation matrices is the product $U_{j}^{n_{1}n_{2}}V_{j}^{n_{1}n_{2}*}$. To
incorporate scaling, we use the geometric mean of the singular values,
$\sqrt{s_{j,1}^{n_{1}n_{2}}s_{j,2}^{n_{1}n_{2}}}$, where
$s_{j,1}^{n_{1}n_{2}}$ and $s_{j,2}^{n_{1}n_{2}}$ are the singular values on
the diagonal of $S_{j}^{n_{1}n_{2}}$. Thus, the scaled rotation matrix becomes
$\sqrt{s_{j,1}^{n_{1}n_{2}}s_{j,2}^{n_{1}n_{2}}}U_{j}^{n_{1}n_{2}}V_{j}^{n_{1}n_{2}*}$.

The diagonal blocks of the state transition matrix, $A_{j}^{k}$ in equation
2.11 are not purely scaled rotational matrices; instead, each one is a scaled
rotational matrix from which a scaling factor times an identity matrix has
been subtracted. The scaling factors associated with the identity matrix can
be obtained, as they derive from the off-diagonal scaling factors.
Consequently, performing SVD on
$A_{j}^{k}+\sum_{n_{1}=k}\sqrt{s_{j,1}^{n_{1}n_{2}}s_{j,2}^{n_{1}n_{2}}}\cdot
I_{2\times 2}$ yields $U_{j}^{k}S_{j}^{k}V_{j}^{k*}$. The updated diagonal
block then becomes
$\sqrt{s_{j,1}^{k}s_{j,2}^{k}}U_{j}^{k}V_{j}^{k*}-\sum_{n_{1}=k}\sqrt{s_{j,1}^{n_{1}n_{2}}s_{j,2}^{n_{1}n_{2}}}$,
where $s_{j,1}^{k}s_{j,2}^{k}$ are the singular values on the diagonal of
$S_{j}^{k}$. This approach ensures that, with each iteration of the EM
algorithm, the updated A matrices maintain a desired rotational structure.

Similarly, we perform the same process on each off-diagonal block of
$\Sigma_{j}$ as outlined in equation 2.16 during the M-step to ensure that the
rotational structure is maintained. Note that there is no concern about the
diagonal blocks in $\Sigma_{j}$, as they are diagonal matrices that simply
represent the covariance of each oscillator.

###  2.6 Theoretical Cross-spectrum of the State-Space Model

Understanding the relationships in the frequency domain between time series is
important for unraveling the latent dynamics that govern a system’s behavior.
Frequency domain methods decompose the signal into sinusoids of specific
frequencies, and describe the relative contribution of each of these in
explaining the total variability in the signal (Kass et al., 2014). Frequency
domain methods are extremely common in the analysis of rhythmic neural data
(Greenhouse et al., 1987; Logothetis et al., 2001; Chaumon et al., 2009; Al-
Fahoum and Al-Fraihat, 2014). The theoretical cross-spectrum is defined as the
Fourier transform of the theoretical cross-covariance function. We derive a
parametric solution for the cross-spectrum of our general state-space
oscillator model, which allows us to visualize the coordinated behavior of any
pair of neural sources. Recall that the latent process model is

$(x_{t}|x_{t-1},S_{t}=j)\sim N(A_{j}x_{t-1},\Sigma_{j})$

and the observation model is

$(y_{t}|x_{t},S_{t}=j)\sim N(B_{j}x_{t},R).$

Conditioning on $S_{t}=j$, we suppress subscript $j$ in the following
derivation. The latent oscillatory state $x_{t}$ follows an AR(1) process. The
AR($k$) model is given by

$x_{t}+a_{1}x_{t-1}+\ldots+a_{k}x_{t-k}=u_{t}.$ (2.20)

We set $a_{1}=-A$, $a_{2},\ldots,a_{k}=0$. Then, the transfer function matrix
of this model (Priestley, 1981) is given by

$\Gamma_{x}(\omega)=\left[\alpha(e^{-i\omega})\right]^{-1},$ (2.21)

where the matrix polynomial $\alpha(z)$ is defined by
$\alpha(z)=\sum_{p=0}^{k}a_{p}z^{p}$ and $a_{0}=I$. By plugging in $a_{1}=-A$,
$a_{2},\ldots,a_{k}=0$, we have

$\Gamma_{x}(\omega)=(I-Ae^{-i\omega})^{-1}.$ (2.22)

Then, following (Priestley, 1981), we obtain the spectral matrix of $x_{t}$ as

$h_{x}(\omega)=\frac{1}{f_{s}}\Gamma_{x}(\omega)\Sigma\Gamma_{x}^{*}(\omega),$
(2.23)

where $f_{s}$ is the sampling frequency, $\omega$ is the frequency of
interest, and $*$ denotes the conjugate transpose.

To compute the spectral matrix of observations $y_{t}$, we first consider the
spectral matrix of $Bx_{t}$, which can be represented as
$\sum_{p=-\infty}^{\infty}\mathbf{g}(p)x_{t-p}$ by setting $g(0)=B$ and
$g(p)=0$ for $p\neq 0$, where $g(p)=\\{\mathbf{g}_{ij}(p)\\}$ is referred to
as the impulse response matrix. Following (Priestley, 1981), we derive the
transfer function matrix as

$\Gamma_{Bx,ij}(\omega)=\sum_{p=-\infty}^{\infty}\mathbf{g}_{ij}(p)e^{-i\omega
p}=B_{ij}e^{-i\omega\cdot 0}=B_{ij}.$ (2.24)

That is, the transfer function matrix $\Gamma=B$, then the spectral matrix of
$Bx_{t}$ is given by

$h_{Bx}(\omega)=Bh_{x}(\omega)B^{T}.$ (2.25)

Given the observation equation $y_{t}=Bx_{t}+v_{t}$, where $v_{t}$ represents
independent Gaussian noise, the spectral matrix of the observations $y_{t}$
can be expressed as the sum of the spectral matrix of $Bx_{t}$ and the
spectral contribution of the noise $v_{t}$. Hence, the spectral matrix of the
observations $y_{t}$ is

$h_{y}(\omega)=Bh_{x}(\omega)B^{T}+\frac{1}{f_{s}}R\quad\text{for
all}\,\omega.$ (2.26)

With the spectral matrix of the observations, we can then compute the
normalized coherence of the observations as follows:

$C_{ij}(\omega)=\frac{h_{ij}(\omega)}{\sqrt{h_{ii}(\omega)h_{jj}(\omega)}},$
(2.27)

where $h_{ij}(\omega)$ is the $ij$-th element of $h_{y}(\omega)$, for
$i,j=1,\ldots,N$.

##  3 Simulation Study

We conducted a simulation study to demonstrate the application and
interpretation of the three model structures introduced in Sections 2.1, 2.2,
and 2.3. In section 3.1, we generate example data from small networks
generated from each model class and demonstrate how inference from the
corresponding model is able to capture the true structure of each network and
the times when the networks switch. In section 3.2, we simulate data from
multiple larger networks under each model structure and compare inferences
from both the correctly specified and misspecified models as well as standard
time-frequency estimators. This process aims to highlight each model’s
properties and potential applications, rather than to perform a systematic
comparative analysis among these three model structures.

###  3.1 Model Visualization

We simulate data for each model structure described in Sections 2.1, 2.2, and
2.3. For each model, we generate data with four nodes and three different
switching states. For illustration purposes, we force two switches to occur
only at times 80 seconds and 200 seconds, instead of simulating the switching
states directly from the model. Note that the estimation procedures do not
know when these switches occur and still attempt to estimate switch times. We
utilize the EM algorithm to estimate the model parameters as described in
Section 2.5. For each set of simulated data, we fit it to its corresponding
model, assuming all the parameters are known except for the one that reflects
the network structures. Specifically, for the COM, the CNM, and the DIM, we
fit only $B_{j}$, $\Sigma_{j}$, and $A_{j}$, respectively. The initial
parameters, whether $B_{j}$, $\Sigma_{j}$, and $A_{j}$ are set without any
specific linkage structure or phase relationship.

####  3.1.1 Common Oscillator Model

Starting with the Common Oscillator Model introduced in Section 2.1, the
network structure is captured by the observation matrix $B_{j}$, with nodes
driven by the same oscillators linked in a network mode. In this simulation,
there are four electrodes, two latent states oscillating with a center
frequency of 7 Hz, and three switching states, each representing a distinct
network structure. The system switches from state 1 to state 2 at 80 seconds
and from state 2 to state 3 at 200 seconds. In state 1, node 1 is driven by
oscillator 1, and node 2 by oscillator 2, which means these two nodes are not
connected. In state 2, nodes 1 and 2 are linked through oscillator 1 with a
90-degree phase difference, while nodes 3 and 4 are connected through
oscillator 2, also with a 90-degree phase difference. In state 3, nodes 1, 2,
and 3 are driven by the same oscillator, with nodes 1 and 2 being 180 degrees
out of phase, nodes 2 and 3 are also 180 degrees out of phase, and nodes 1 and
3 are in phase. Figure 4A shows the diagram of the three states. If any pair
of nodes are connected, they are linked with an edge. We color only the edge
between nodes 1 and 2 in orange and label it with weights, indicating linkage
strength between nodes 1 and 2, which decreases from 0.3 in state 2 to 0.25 in
state 3. This specific coloring is chosen because Figure 4B focuses on
presenting the coherogram for the pair of nodes 1 and 2.

Figure 4: Common Oscillator Model. (A) Diagram of the network structures for
states 1, 2, and 3, illustrating how nodes are interconnected in each state.
The edge between nodes 1 and 2 is highlighted in orange and annotated with
weights, as this pair is analyzed in the coherogram. (B) Coherograms for nodes
1 and 2. The top panel is the theoretical coherogram based on the true network
structures and the actual occurrences of switches; the second is generated
from a non-parametric multi-taper method estimate of coherence over a sliding
window; and the third is the theoretical coherogram based on the estimated
parameter, $B_{j}$, and the estimated switching states. (C) Visualization of
the true switching states (red) and the estimated switching states (blue). (D)
Visualization of the magnitudes from both the true (left two panels) and
estimated (right two panels) observation matrices, $B_{j}$, highlighting the
true and estimated network structure for each state. (E) Visualization of the
phase relationship from the true (left two panels) and estimated (right two
panels) observation matrices, $B_{j}$.

In Figure 4B, the top panel shows the true theoretical coherogram, accounting
for the actual network structure in each state and the occurrences of
switches. In states 2 and 3, nodes 1 and 2 are coherent at around 7 Hz, and
from state 2 to state 3, the strength of their linkage decreases. The
coherogram in the middle panel was generated using a non-parametric multi-
taper estimator with a moving window of 1 s and a bandwidth of 2 Hz. We use 1
second windows to have a hope of estimating network switches from second to
second, but the resulting estimator is extremely noisy, making it difficult to
tell whether the nodes are linked at any frequency. The bottom panel is the
theoretical coherogram, constructed based on the estimated parameters of the
model and the estimated switching states. The estimates from the COM clearly
indicate that nodes 1 and 2 are linked at around 7 Hz in states 2 and 3 and
are able to detect the change in linkage strength at time 200 seconds. In
Figure 4C, the true switching states are colored in red, and the estimated
switching states are in blue. The estimated states align almost perfectly with
the ground truth.

In Figure 4D, we visualize the magnitude of the parameters that describe the
degree to which each oscillator is expressed at each node from equation 2.6 in
Section 2.1. The true network structure is shown in the left two panels, where
the first and second columns show the extent to which oscillators 1 and 2 are
expressed at each node, respectively. The right two panels are the estimated
network structures, with two columns representing the two oscillators. Each
row represents a single switching state. The estimated model correctly infers
that in switching state 1, the first oscillator drives the first node, the
second oscillator drives the second node, and these nodes are not linked.
Similarly in switching state 2, the model correctly estimates that nodes 1 and
2 are linked through oscillator 1 and that nodes 3 and 4 are linked through
oscillator 2. In switching state 3, the model correctly infers that states 1,
2, and 3 are all linked through a common oscillator. In Figure 4E, we
visualize the relative phases of each pair of nodes in a similar manner. For
nodes that are driven by the same oscillator, we can compute their phase
difference as described in equation 2.6. From state 2 to state 3, nodes 1 and
2 are linked with different magnitudes and phases, which can be effectively
captured by the COM. By comparing the panels for ground truth with the
estimated results, we find that the COM is effective at not only detecting the
switches but also estimating the linkage structures and the phase
relationships.

In a regime with a low signal-to-noise ratio, the conventional windowed
method, the multi-taper coherogram, fails to detect the linkage structure, as
shown in the second panel of Figure 4B. As previously mentioned, the multi-
taper coherogram does not assume any specific network structures and allows
for any type of link at any time. In contrast, the modeling structures we
utilize assume that the system switches between a small number of distinct
networks. With switching components, the network structure for each state is
estimated over all the time points when that switching state value is
correctly estimated. By doing so, we improve the statistical power for
estimating network structures and can identify which network is present at any
given time, as illustrated in the third panel of Figure 4B.

####  3.1.2 Correlated Noise Model

Moving to the Correlated Noise Model introduced in Section 2.2, the network
structure is captured by the process noise covariance matrix $\Sigma_{j}$,
with the connectivity of nodes arising from shared noise structure. In this
simulation there are four electrodes, and therefore, four latent states
oscillating with a central frequency of 7 Hz, where each node is associated
with its own oscillator, and three switching states, each representing a
unique network structure. Similar to the COM, the system transitions from
state 1 to state 2 at 80 seconds, and from state 2 to state 3 at 200 seconds.

Figure 5: Correlated Noise Model. (A) Diagram of the network structures for
states 1, 2, and 3, illustrating how nodes are interconnected in each state.
The edge between nodes 1 and 2 is highlighted in orange and annotated with
weights, as this pair is analyzed in the coherogram. (B) Coherograms for nodes
1 and 2. The top panel is the theoretical coherogram based on the true network
structures and the actual occurrences of switches; the second is generated
from a non-parametric multi-taper method estimate of coherence over a sliding
window; and the third is the theoretical coherogram based on the estimated
parameter, $\Sigma_{j}$, and the estimated switching state. (C) Visualization
of the true switching states (red) and the estimated switching states (blue)
(D) Visualization of linkage magnitudes from the true (top panel) and
estimated (bottom panel) noise covariance matrices, $\Sigma_{j}$, indicating
the network structures for each state. (E) Visualization of the phase
relationships from the true (top panel) and estimated (bottom panel) noise
covariance matrices, $\Sigma_{j}$.

In state 1, there are no connections among nodes. In state 2, nodes 1 and 2
are linked with a 90-degree phase difference, and nodes 3 and 4 are connected
with a 270-degree phase difference. In state 3, nodes 1, 2, and 3 are
interconnected, all in phase. Figure 5A shows the diagram of the three states.
Once again, we focus on the pair of nodes 1 and 2, with the corresponding
coherograms shown in Figure 5B. The linkage strength for this pair decreases
from 0.5 in state 2 to 0.4 in state 3.

In Figure 5B, we display the true theoretical coherogram (top), the multi-
taper coherogram (middle), and the theoretical coherogram from the estimated
model (bottom). The multi-taper coherogram is noisy and it is impossible to
tell that nodes 1 and 2 are linked at 7 Hz. In contrast, the estimate from the
CNM clearly shows that nodes 1 and 2 are linked at around 7 Hz in states 2 and
3, and successfully detects the change in linkage strength at 200 seconds. In
Figure 5C, the estimated switching states (blue) nearly perfectly align with
the true values of the switching states (red).

In Figure 5D, the true and estimated linkage strengths from equation 2.8 are
presented in the top and bottom panels, respectively. From left to right, each
column represents a different switching state. The true links are accurately
detected by the model estimates. In Figure 5E, we display the true and
estimated phase relationships for each pair of nodes, with the true
relationships at the top and the estimated ones at the bottom. From state 2 to
state 3, nodes 1 and 2 are linked, with their linkage strength decreasing from
0.5 to 0.4 and phase relationship changing from being 90-degree out of phase
to in phase. These changes are captured by the CNM.

####  3.1.3 Directed Influence Model

For the Directed Influence Model (DIM) introduced in Section 2.3, we define
the network structure via the process transition matrix $A_{j}$, simulating
direct influences among oscillators. Similar to the CNM, each node corresponds
to an individual oscillator. In this simulation there are four electrodes, and
hence four latent states oscillating with a central frequency of 7 Hz, and
three distinct switching states. the system switches from state 1 to state 2
at 80 seconds, and from state 2 to state 3 at 200 seconds.

Figure 6: Directed Influence Model. (A) Diagram of the network structures for
states 1, 2, and 3, illustrating how nodes are interconnected in each state.
The edge between nodes 1 and 2 is highlighted in orange and annotated with
weights, as this pair is analyzed in the coherogram. (B) Coherograms for nodes
1 and 2. The top panel is the theoretical coherogram based on the true network
structures and the actual occurrences of switches; the second is generated
from a non-parametric multi-taper method estimate of coherence over a sliding
window; and the third is the theoretical coherogram based on the estimated
parameter, $A_{j}$, and the estimates of switching states. (C) The true
switching states are (red) and the estimated switching states (blue). (D)
Visualization of magnitude from the true (top panel) and estimated (bottom
panel) process transition matrix $A_{j}$, indicating the network structure for
each state. (E) Visualization of the phase relationship from the true (top
panel) and estimated (bottom panel) process transition matrix, $A_{j}$.

In state 1, all nodes remain unconnected. Transitioning to state 2, a directed
link is established from node 1 to 2 with a 90-degree phase difference. At the
same time, nodes 3 influences 4 with a 270-degree phase difference. In state
3, node 1 influences 2 and 3, and node 2 influences 4, all in phase. The
diagram in Figure 6A illustrates these networks, with a particular emphasis on
the dynamic between nodes 1 and 2, as detailed in the coherograms of Figure
6B, where we observe a decrease in linkage strength from 0.4 in state 2 to 0.2
in state 3. As previously mentioned, the DIM is the only model structure we
introduce that has directionality. The arrow in the diagram indicates the
direction of influence.

Figure 6B presents a series of coherograms: the true theoretical, the multi-
taper, and the theoretical based on the estimated model. The multi-taper
estimate does not indicate a linkage at 7 Hz between nodes 1 and 2 at any
time. On the contrary, the estimates from the DIM effectively reveal the
connection at 7 Hz between these two nodes in states 2 and 3, capturing the
change in linkage strength. The true switching states and the estimated
switching states are shown in Figure 6C. The estimated states align almost
perfectly with the ground truth.

In Figure 6D, we visualize both the actual and estimated linkage strengths in
the network structures from equation 2.11 at the top and bottom, respectively.
From left to right, each column represents a different switching state. Figure
6E shows the true and estimated phase relationships for each pair of nodes,
again with the true relationships at the top and the estimated ones at the
bottom. From state 2 to state 3, nodes 1 and 2 are linked, with their linkage
magnitude decreasing from 0.4 to 0.2 and phase relationship changing from
90-degree out of phase to in phase. These are captured by the DIM.

The consistent feature in each of these simulations is that the model
estimates are able to capture both the network structure in each state and the
switches between these networks. The sliding window estimator is able to
estimate neither, because short windows do not provide enough data to estimate
linkages, and long windows cannot track rapid changes in the network
structure. Conversely, by positing that the system switches between a small
number of distinct networks, we notably improve our ability to detect the
presence and evolution of links between nodes.

###  3.2 Model Performance Evaluations

In this subsection, we analyze simulations from larger networks, simulating
the switching state $S_{t}$, the latent oscillatory state $x_{t}$, and the
observation $y_{t}$ under each model structure, and comparing the estimates
using both correctly specified and misspecified models. For each simulation,
the number of neural sources is set to ten, the oscillation frequency is
centered at 7 Hz, and the sampling frequency is set to 100 Hz. Similar to the
previous subsection, we assume all the parameters are known except for the
ones that reflect the network structures, and we use the EM algorithm to
estimate these parameters. The initial parameters, whether $B_{j}$,
$\Sigma_{j}$, or $A_{j}$ are set without any specific linkage structure or
phase relationship.

We examine two evaluation metrics. First, we use the estimated parameters from
each model to compute the theoretical cross-spectrum of the fitted model, as
described in Section 2.6. We first compute the element-wise difference between
the estimated and true cross-spectral matrices at the frequency of interest,
and then calculate the root mean square of the magnitudes of all the off-
diagonal elements in the difference matrix:

$\left\|h_{y}^{\text{est}}(w)-h_{y}^{\text{true}}(w)\right\|=\sqrt{\frac{1}{N(N-1)}\sum_{i\neq
j}|h_{y}^{\text{est}}(w)-h_{y}^{\text{true}}(w)|_{ij}^{2}},$ (3.1)

where $h_{y}^{\text{est}}(w)$ is the estimated cross-spectral matrix,
$h_{y}^{\text{true}}(w)$ is the true cross-spectral matrix, $N$ is the number
of nodes, and $w$ is the frequency of interest (specifically, $w$ = 7 Hz in
simulations). The normalization term is $1/N(N-1)$ because we focus on the
off-diagonal terms in the $N\times N$ cross-spectrum matrix. Secondly, we use
the fitted parameters to compute the coherence matrix, as described in Section
2.6. We then conduct a hypothesis test on every pair of nodes to determine if
they are significantly coherent. This is achieved by fitting an empirical
gamma distribution to the off-diagonal elements in the estimated coherence
matrix, using a significance level of 0.05. We report the sensitivity and
false positive rate. Sensitivity is calculated by computing the proportion of
actual coherent links that are correctly identified by our model (true
positives) out of all the true coherent links present, according to the true
coherence matrix. This metric evaluates how effectively our model identifies
coherent pairs of nodes when they are indeed coherent. The corresponding
formula for sensitivity is given by:

$\text{Sensitivity}=\frac{\text{True Positives (TP)}}{\text{True Positives
(TP)}+\text{False Negatives (FN)}}$

On the other hand, the false positive rate (FPR) computes the ratio of non-
coherent links that are wrongly identified as coherent (false positives) out
of all the actual non-coherent pairs, measuring the model’s tendency to
incorrectly label non-coherent pairs as coherent. The formula for FPR is:

$\text{FPR}=\frac{\text{False Positives (FP)}}{\text{\text{False Positives
(FP)}}+\text{True Negatives (TN)}}$

For assessing the significance of the multi-taper coherence, we perform the
F-test with test statistics equal to $(v_{0}/2-1)|C|^{2}/(1-|C|^{2})$, which
is distributed as $F_{2,v_{0}-2}$ under the null hypothesis of zero population
coherence (Hannan, 1970; Jarvis and Mitra, 2001). Here, $v_{0}/2-1$ is the
number of tapers and $|C|$ is the magnitude of the multi-taper coherence at
the frequency of interest. The derivation of the sample coherence distribution
is given in (Hannan, 1970), and details of the test statistics can be found in
(Jarvis and Mitra, 2001). Similarly, we compute the sensitivity and false
positive rate for this test.

For example, if the generative model is the COM, we obtain three estimates of
the cross-spectral matrix: one from the COM itself, one from the CNM, and one
from the DIM. Additionally, we compute the multi-taper cross-spectral matrix
in 1 second windows throught the simulation. To compare the cross-spectral
estimates from the COM, CNM, and DIM with the true cross-spectral matrix, we
compute the norm of the difference matrix, as defined in Equation 3.1, at
every time point. We then report the mean and standard deviation of these norm
values in Table 1. Recall that, by conditioning on the switching state
$S_{t}=j$, we can compute the cross-spectral matrix as described in Section
2.6. In our example, with three switching states, this yields three estimated
cross-spectral matrices. At every time step, we obtain a single estimated
cross-spectral matrix by calculating the weighted sum of cross-spectral
matrices from the three states, with the weights based on the probability of
the estimated switching state. Similarly, to compare the multi-taper cross-
spectral estimates to the true cross-spectral matrix, we compute the norm of
the difference matrix at every window and report the mean and standard
deviation of these norm values, as shown in the last row in Table 1. The
results of the sensitivity and false positive rate for testing coherence
between any pair of nodes in the estimated coherence matrix are presented in
Tables 2 and 3, respectively.

Table 1: Mean and standard deviation (in parentheses) of the element-wise error between the estimated and true cross-spectral matrices (off-diagonal elements). The values are only comparable within each column. Fitted model | Generative model
---|---
DIM  | CNM  | COM
DIM | 0.0069 (0.0007)  | 0.0192 (0.0066)  | 0.0030 (0.0008)
CNM | 0.0093 (0.0013)  | 0.0085 (0.0011)  | 0.0030 (0.0011)
COM | 0.0145 (0.0014)  | 0.0180 (0.0018)  | 0.0015 (0.0003)
Multi-taper | 0.0221 (0.0021)  | 0.0313 (0.0040)  | 0.0072 (0.0014)
Table 2: Sensitivity of tests for coherence between pair of nodes with true links, with fractions shown in parentheses. Fitted model | Generative model
---|---
DIM | CNM | COM
DIM | 1.000 (40/40) | 1.000 (54/54) | 0.917 (22/24)
CNM | 1.000 (40/40) | 1.000 (54/54) | 0.917 (22/24)
COM | 0.900 (36/40) | 0.963 (52/54) | 1.000 (24/24)
Multi-taper | 0.345 (1404/4066) | 0.364 (1966/5396) | 0.446 (1114/2496)
Table 3: False positive rate of tests for coherence between pairs of nodes without a link, with fractions shown in parentheses. Fitted model | Generative model
---|---
DIM | CNM | COM
DIM | 0.052 (12/230) | 0.056 (12/216) | 0.057 (14/246)
CNM | 0.043 (10/230) | 0.046 (10/216) | 0.049 (12/246))
COM | 0.052 (12/230) | 0.046 (10/216) | 0.041 (10/246)
Multi-taper | 0.097 (2230/22934) | 0.129 (2786/21604) | 0.074 (1816/24504)

In Table 1, for data generated by any one of these models—DIM, CNM, or COM —
estimates from and of these three models all have smaller error values than
those of the multi-taper method, even when the estimated model is
misspecified. In Table 2, we observe a massive improvement in the ability to
identify the linkage structures in the data when compared to the multi-taper
method, again, even for misspecified models. Note that for the multi-taper
method, we perform the coherence test on each window of the tapered segments
of the signals. In both tables, it is evident that when the data is fitted
with its corresponding generative model, these models successfully capture the
structures they are designed to detect and have the smallest error values in
each set of the simulated data. Table 3 shows that the false positive rates
for all of the models are well controlled. However, the false positive rates
for the multi-taper method are less well controlled. For example, the FPR for
multitaper estimates from data generated by the CNM is 0.129, which does not
adequately control the FPR to the desired significance level of 0.05. The
potential reason for this is that each test is based on a small window,
resulting in a lack of accuracy in the approximated F distribution.

As expected, we consistently find that estimation of network structure using
the true model that generated the data performs as well or better than
estimation from the misspecified models. This suggests that none of these
models is superior to the others irrespective of the structure in the data.
However, what is more striking is the fact that estimation of the network
structure is nearly perfect in these simulations, even when using misspecified
models. This highlights the flexibility and general applicability of the
models, even in cases where there is minimal prior knowledge about the
mechanisms underlying functional connectivity. Each of these model structures
takes advantage of the increased statistical power gained by combining
information over time. Conversely, if one has specific hypotheses about the
mechanisms of functional connectivity, drawing inferences from the
corresponding model leads to even better statistical power.

For each simulated dataset and fitted model, we also report the accuracy of
the switching state estimate by calculating the percentage of time steps where
the state is correctly estimated, as shown in Table 4. These models are able
to detect the switching states almost perfectly when using the corresponding
generative model and still perform well, with accuracy over 95%, when using a
misspecified model. Note that the denominator is 30,000 because we simulate
300 seconds of data with a sampling frequency of 100 Hz, resulting in 30,000
data points.

Table 4: Fraction of switching states that are correctly identified, with frequencies shown in parentheses. Fitted model | Generative model
---|---
DIM | CNM | COM
DIM | 0.999 (29970/30000) | 0.967 (29001/30000) | 0.999 (29970/30000)
CNM | 0.998 (29940/30000) | 0.999 (29977/30000) | 0.999 (29967/30000)
COM | 0.997 (29906/30000) | 0.995 (29844/30000) | 0.999 (29968/30000)

Figure 7: Model evaluation metrics vs the signal to noise ratio (dB) of the
simulated data.

The results from Tables 1-4 are in a simulation regime with low signal-to-
noise, where tests based on the multi-taper coherogram consisently fail to
properly detect the linkage structure in the data. Next, we examined the
models’ performance by adjusting the signal-to-noise ratio (SNR) and observing
how error in the estimate cross-spectrum, sensitivity for detecting links, and
switching accuracy are affected. In Figure 7, each column corresponds to a
different generative model for the simulated data. The top row shows the error
in the estimated cross spectrum when fitting with each model or using the
multitaper estimator. As expected, the error values decrease as the SNR
increases for every estimation method. We observe that the multi-taper
estimator consistently has higher error values than any of the three models.
When the fitted model is correctly specified, it consistently has the smallest
error values, but the differences between the models tends to be small
compared to the error of the multitaper estimator. The middle row shows the
sensitivity for detecting true links using each model and the multitaper
estimator. While all of the estimators have low sensitivity in our lowest SNR
regime, we observe that the sensitivity of the multi-taper method increases
much more slowly as the SNR increases than do the model-based estimates. The
bottom row shows the accuracy of the estimated switch state for each of the
models. Note that the multitaper method does not provide a simple method for
estimating which network is expressed at each time. Not surprisingly, we find
that switch detection is inaccurate at the lowest SNR regimes and increases
rapidly with SNR. Interestingly, we find that for the correctly specified
model, switching state estimation starts increasing at lower SNR levels than
sensitivity of estimating network structure. this suggests that even when
completely accurate inference of the networks is impossible, it may still be
possible to detect when changes in network structure occur.

##  4 Discussion

In this paper, we introduce a broad modeling framework for analyzing dynamical
functional connectivity in rhythmic neural data, and investigate three
specific model structures, each of which describes different mechanisms of
coupling. Broadly, the model framework assumes that a small number of discrete
network structures are present at any given time, and maximizes statistical
power by estimating network structure over all periods where the same network
is determined to be expressed, and determining these periods using the
improved network estimates. These approaches allow us to more effectively
analyze dynamically coupled oscillators by estimating the underlying
functional network structure and identifying changes in their expression. Our
purpose in proposing multiple model structures is not to suggest that any one
of these is superior to the others, but to demonstrate the inferential value
of each and to highlight the improvement in statistical power derived from the
broader model framework.

We believe that the distinct model structures proposed are suited to
addressing distinct physiological questions and have the potential to provide
different inferential insights. From a physiological perspective, researchers
might be interested in specific questions about the mechanisms by which
functional oscillatory networks are formed in the brain. Possible explanations
include generation from a single source that propagates outward, generation
across multiple locations influenced by shared, unobserved inputs, or
generation through the direct action of anatomically linked sources. These
scenarios might guide them towards the Common Oscillator Model (COM), the
Correlated Noise Model (CNM), and the Directed Influence Model (DIM),
respectively. From an inferential perspective, researchers might be interested
in describing high-dimensional activity within a large connected network using
a low-dimensional manifold (i.e., a small set of oscillators). Conversely,
researchers might be interested in making inferences about directional or
causal links in a network and their role in shaping coordinated activity. In
such cases, the choice of a model structure would reflect the inferential
goals of the analysis rather than physiological prior knowledge about data
generation mechanisms. The choice of which model to apply depends on the
specific questions being asked. However, our simulation results suggest that
the inferences from these models are robust, even when there is a mismatch
between the model used for fitting and the true data generating mechanism.

Another advantage of this modeling framework is its ability to accommodate
oscillators at multiple frequencies. While our simulation study focused on
networks linked at a single central frequency, it is simple to include
oscillators each with its own central frequency to model transitions between
networks that are dominated by distinct rhythms. For example, in the analysis
of dynamic brain networks during propofol anesthesia, both alpha rhythms and
slow waves are though to entrain distinct networks during different stages of
unconsciousness (Stephen et al., 2020; Mukamel et al., 2014; Purdon et al.,
2013). Another example is the analysis of signals in the basal ganglia of
Parkinson’s disease (PD) patients, where distinct networks driven by beta and
gamma oscillations have been observed to predict features of motor control
tasks (Hirschmann et al., 2011; Litvak et al., 2011; Swann et al., 2016;
McGregor and Nelson, 2019). The most commonly reported change is an increase
in beta band (13-30 Hz) activity in the motor cortex and basal ganglia when
the patient is experiencing motor impairments such as rigidity and slowing of
movement (bradykinesia) (Little and Brown, 2014; Vinding et al., 2020).
Conversely, during movement initiation or execution, there is a notable
decrease in beta oscillations and an increase in gamma band activity,
reflecting the kinetic processes underlying motor control (Brown, 2003; van
Wijk et al., 2017). The multiplicity of frequencies is crucial for modeling
and understanding the physiological basis of different brain states. By
allowing for the simultaneous representation of these different frequencies,
this modeling framework enables us to better interpret the complex patterns of
brain activity.

There are a number of ways in which this work may be extended. For the COM,
the number of oscillators can be much smaller than the number of neural
sources being recorded. In each of our simulations, we fixed this number and
used the correct number of oscillators when using the COM to estimate network
structure. In prior work, (Hsin et al., 2022) we analyzed EEG data during
anesthesia induction and return to consciousness using the COM with four
oscillators - one for the alpha rhythm and three for the slow waves - to
describe the data across 64 electrodes. It is noteworthy that tools have been
developed for estimating the number of oscillators from time series data based
on this model framework. This includes the oscillator decomposition method
(Matsuda and Komaki, 2017), which has been employed to extract oscillation
components from functional near-infrared spectroscopy (fNIRS) data (Matsuda et
al., 2022). The authors highlight that this method is applicable to local
field potential (LFP), EEG, magnetoencephalography (MEG), and
electrocorticography (ECoG) data as well. In addition, a probabilistic
parametric generative model (Das and Purdon, 2022) offers a Bayesian inference
procedure for identifying the optimal number of oscillation sources through
empirical Bayes model selection. The effectiveness of this approach has been
demonstrated on both simulated and real EEG data. If there is no prior
knowledge regarding the optimal number of oscillators for use in the COM, the
aforementioned approaches can be referred to determine the appropriate number
of oscillators.

The broad model framework we proposed was based on the assumption that one
network state is expressed at a time, and the system transitions between
different network structures nearly instantaneously when a switch occurs. Both
the latent oscillatory model and the observation model are conditioned on a
single switching state at any given moment. A natural extension would be to
model the data as a superposition of multiple states, where all network modes
are present at any time, but the degree of expression of each network
structure changes over time. In this case, instead of using a latent discrete
random variable to describe a switching state, we could use a multi-
dimensional continuous random vector to represent expression of each network
at each time step. Then, the latent oscillatory state and the observation
model would no longer be conditioned on a single state; instead, they would
depend on the expression level of all the network modes at the current time.
For example, it is reasonable to model different states of motor control in
Parkinson’s disease as a changing balance between beta- and gamma-driven
networks rather than as an instantaneous transition between the two.

Modeling high dimensional rhythmic neural activity using a state space
paradigm, with latent variables that take advantage of known oscillatory
dynamics, provides a general and robust framework for understanding shared
oscillatory dynamics. This broader analytical framework allows us to capture
various functional network structures, each offering insights into the
mechanisms and structures underlying coordinated oscillations. These latent
process models can be effective in accurately estimating the functional
network structure across a wide range of coupling mechanisms and identifying
dynamic changes in network expression related to clinical or cognitive
context.

## Code Availability

The code for all the methods described and simulation studies performed in the
paper is publicly available through a MATLAB code repository on GitHub at
https://github.com/Stephen-Lab-BU/Switching_Oscillator_Networks/

## Acknowledgments

This work was supported by the Simons Collaboration on the Global Brain
(521921 and NC-GB-CULM-00002730-06). The authors would like to thank Proloy
Das for helpful discussions on simulation studies.

## Appendix

### A1: Filter operator

In the standard Kalman filter framework, we assume the latent state $x_{t}$
and the observations $y_{t}$ satisfy the following linear equations.

$\begin{split}x_{t}&=Ax_{t-1}+u_{t}\\\ y_{t}&=Bx_{t}+v_{t},\end{split}$ (4.1)

where $u_{t}\sim N(0,\Sigma)$ and $v_{t}\sim N(0,R)$. Here, we summarize the
filter operator

$(x_{t},V_{t},L_{t})=\text{Filter}(x_{t-1},V_{t-1},y_{t},A,\Sigma,B,R).$

At each time step $t$, it has two steps: (1) prediction of the state and its
error covariance estimates and (2) update these estimates with new measurement
data. The computation is based on a recursive approach for the incorporation
of the new observation $y_{t}$. The steps are as follows. We first compute the
predicted state $x_{t|t-1}^{-}$ and its error covariance matrix
$V_{t|t-1}^{-}$ as follows:

$\begin{split}x_{t|t-1}^{-}&=Ax_{t-1}\\\
V_{t|t-1}^{-}&=AV_{t-1}A^{{}^{\prime}}+\Sigma\end{split}$ (4.2)

We then compute the innovation $e_{t}$ (the state error in the prediction),
the Kalman gain matrix $K_{t}$, the updated state estimate $\hat{x}_{t}$, and
the updated state error covariance matrix $V_{t}$.

$\begin{split}e_{t}&=y_{t}-Bx_{t|t-1}^{-}\\\
K_{t}&=V_{t|t-1}^{-}B^{{}^{\prime}}(BV_{t|t-1}^{-}B^{{}^{\prime}}+R)^{-1}\\\
\hat{x}_{t}&=x_{t|t-1}^{-}+K_{t}e_{t}\\\
V_{t}&=(I-K_{t}B)V_{t|t-1}^{-}\end{split}$

In the prediction step, the equations use the system’s model to forecast the
next state and its error covariance. The update step adjusts these predictions
based on the new measurement $y_{t}$, with the Kalman Gain $K_{t}$ determining
the weight of the adjustment. This process iteratively refines the state
estimated by balancing between the predictions and the measurements. The
likelihood $L_{t}$ of observing the observation $y_{t}$ can be computed by

$L_{t}=N(e_{t};0,BV_{t|t-1}^{-}B^{{}^{\prime}}+R).$

### A2: Smoother operator

We denote the smoother operator:

$(\tilde{x}_{t},\tilde{V}_{t},\tilde{V}_{t,t-1})=\text{Smoother}(\tilde{x}_{t+1},\tilde{V}_{t+1},x_{t},V_{t},A,\Sigma).$

To apply the smoother algorithm, we start with the final estimates from the
Kalman filter for the state $x_{T}$ and covaraince $V_{T}$ at the last time
step $T$. Then, moving backward from time step $t=T$ to $t=1$, for the
prediction step, we can pass the estimates in from the filtering stage. Then,
we compute the smoother gain matrix $G_{t}$, the updated state estimate
$\tilde{x}_{t}$, the updated state covariance estimate $\tilde{V}_{t}$, and
the one-step covariance estimate $\tilde{V}_{t,t-1}$ as follows:

$\begin{split}G_{t}&=V_{t}A^{{}^{\prime}}V_{t+1}^{-1}\\\
\tilde{x}_{t}&=\hat{x}_{t}+G_{t}(\tilde{x}_{t+1}-\hat{x}_{t+1})\\\
\tilde{V}_{t}&=V_{t}+G_{t}(\tilde{V}_{t+1}-V_{t+1})G_{t}^{{}^{\prime}}\\\
\tilde{V}_{t,t-1}&=V_{t}G_{t-1}^{{}^{\prime}}+G_{t}(\tilde{V}_{t+1,t}-Av_{t})G_{t-1}^{{}^{\prime}}\end{split}$

The computation of the one-step covariance estimates shown above is an
alternative way that does not required the corresponding filtered terms
(Shumway and Stoffer, 1991), where the initial condition is

$\tilde{V}_{T,T-1}=(I-K_{T}B)AV_{T-1}.$

The smoother gain matrix $G_{t}$ is used to adjust the state estimate
$\hat{x}_{t}$ with information from the future state $\tilde{x}_{t+1}$. The
estimates $\tilde{x}_{t}$, $\tilde{V}_{t}$ and $\tilde{V}_{t,t-1}$ represent
the smoothed estimates by taking into account all observations up to time $T$.
This backward pass through all the data allows the smoother algorithm to
provide a more accurate estimate of the state at each time step by
incorporating information from the entire dataset, rather than just the
information available up to that point in time as in the filtering stage.

### A3: Collapse and Cross-Collapse operators

The Collapse and Cross-Collapse operators are used to approximate a mixture of
Gaussians as a single Gaussian by moment matching. Given the conditional mean
for each state is $x_{t}^{j}=\text{E}(x_{t}|S_{t}=j)$, its covariance
$V_{t}^{j}=\text{Cov}(x_{t}|S_{t}=j)$, and the mixing coefficient
$g_{t}^{j}=P(S_{t}=j)$ as the inputs of the operator, we can compute the
unconditional moments as follows. For the collapse operator

$(x_{t},V_{t})=\text{Collapse}(x_{t}^{j},V_{t}^{j},g_{t}^{j}),$

we compute

$\displaystyle x_{t}$ $\displaystyle=\sum_{j}g_{t}^{j}x_{t}^{j}$
$\displaystyle V_{t}$
$\displaystyle=\sum_{j}g_{t}^{j}(V_{t}^{j}-(x_{t}^{j}-x_{t})(x_{t}^{j}-x_{t})^{{}^{\prime}}).$

For the cross-collapse operator

$(x_{t+1},x_{t},V_{t+1,t})=\text{Cross-
Collapse}(x_{t+1}^{j},x_{t}^{j},V_{t+1,t}^{j},g_{t}^{j}),$

we compute

$\displaystyle x_{t+1}$ $\displaystyle=\sum_{j}g_{t+1}^{j}x_{t+1}^{j}$
$\displaystyle x_{t}$ $\displaystyle=\sum_{j}g_{t}^{j}x_{t}^{j}$
$\displaystyle V_{t+1,t}$
$\displaystyle=\sum_{j}g_{t}^{j}(V_{t+1,t}^{j}-(x_{t+1}^{j}-x_{t+1})(x_{t}^{j}-x_{t})^{{}^{\prime}}).$

It has been demonstrated that a Gaussian with these moments is the nearest
possible Gaussian approximation to the original mixture distribution
(Lauritzen, 1996; Murphy, 1998).

## References

  * Akram et al. (2016) Akram, S., Presacco, A., Simon, J. Z., Shamma, S. A., & Babadi, B. (2016).  Robust decoding of selective auditory attention from MEG in a competing-speaker environment via state-space modeling.  _NeuroImage_ , _124_ , 906–917.
  * Al-Fahoum and Al-Fraihat (2014) Al-Fahoum, A. S., & Al-Fraihat, A. A. (2014).  Methods of EEG signal features extraction using linear analysis in frequency and time-frequency domains.  _International Scholarly Research Notices_ , _2014_.
  * Bar-Shalom et al. (2001) Bar-Shalom, Y., et al. (2001).  Estimation with Applications to Tracking and Navigation: Theory, Algorithms and Software.  Wiley.
  * Başar (2013) Başar, E. (2013).  Brain oscillations in neuropsychiatric disease.  In _Dialogues in Clinical Neuroscience_ , _15_(3), 291–300.  doi:10.31887/DCNS.2013.15.3/ebasar
  * Beaver (2006) Beaver, S. (2006).  The singular value decomposition in symmetric orthogonalization and data compression.  Western Oregon University.
  * Beck et al. (2018) Beck, A. M., Stephen, E. P., & Purdon, P. L. (2018).  State Space Oscillator Models for Neural Data Analysis.  In _Annual International Conference of the IEEE Engineering in Medicine and Biology Society_. IEEE Engineering in Medicine and Biology Society. Annual International Conference, 2018, 4740–4743.  doi:10.1109/EMBC.2018.8513215
  * Berger (1929) Berger, H. (1929).  Über das Elektroenkephalogramm des Menschen.  _Archiv für Psychiatrie und Nervenkrankheiten_ , _87_ , 527–570.
  * Brown (2003) Brown, P. (2003).  Oscillatory nature of human basal ganglia activity: Relationship to the pathophysiology of Parkinson’s disease.  _Movement Disorders_ , _18_(4), 357–363.  doi:10.1002/mds.10358
  * Buzsáki (2006) Buzsáki, G. (2006).  _Rhythms of the Brain_.  Oxford University Press.
  * Buzsáki et al. (2012) Buzsáki, G., Anastassiou, C. A., & Koch, C. (2012).  The origin of extracellular fields and currents–EEG, ECoG, LFP and spikes.  _Nature Reviews Neuroscience_ , _13_(6), 407–420.  doi:10.1038/nrn3241
  * Buzsáki and Draguhn (2004) Buzsáki, G., & Draguhn, A. (2004).  Neuronal oscillations in cortical networks.  _Science_ (New York, N.Y.), _304_(5679), 1926–1929.  doi:10.1126/science.1099745
  * Chang and Glover (2010) Chang, C., & Glover, G. H. (2010).  Time–frequency dynamics of resting-state brain connectivity measured with fMRI.  _Neuroimage_ , _50_(1), 81–98.
  * Chaumon et al. (2009) Chaumon, M., Schwartz, D., & Tallon-Baudry, C. (2009).  Unconscious learning versus visual perception: dissociable roles for gamma oscillations revealed in MEG.  _Journal of Cognitive Neuroscience_ , _21_(12), 2287–2299.
  * Chen et al. (2010) Chen, Z., Barbieri, R., & Brown, E. N. (2010).  State space modeling of neural spike train and behavioral data.  In _Statistical Signal Processing for Neuroscience and Neurotechnology_ (pp. 175-218). Academic Press.
  * Chen and Brown (2013) Chen, Z., & Brown, E. (2013).  State Space Model.  _Scholarpedia_ , _8_(3), 30868.  doi:10.4249/scholarpedia.30868
  * Cimenser et al. (2011) Cimenser, A., Purdon, P. L., Pierce, E. T., Walsh, J. L., Salazar-Gomez, A. F., Harrell, P. G., Tavares-Stoeckel, C., Habeeb, K., & Brown, E. N. (2011).  Tracking brain states under general anesthesia by using global coherence analysis.  _Proceedings of the National Academy of Sciences of the United States of America_ , _108_(21), 8832 – 8837.  doi:10.1073/pnas.1017041108
  * Cocchi et al. (2017) Cocchi, L., Gollo, L. L., Zalesky, A., & Breakspear, M. (2017).  Criticality in the brain: A synthesis of neurobiology, models and cognition.  _Progress in Neurobiology_ , _158_ , 132–152.  doi:10.1016/j.pneurobio.2017.07.002
  * Das and Purdon (2022) Das, P., & Purdon, P. L. (2022).  Extracting common oscillatory time-courses from multichannel recordings: Oscillation Component Analysis.  In _2022 56th Asilomar Conference on Signals, Systems, and Computers_ (pp. 602-606), Pacific Grove, CA, USA.  doi:10.1109/IEEECONF56349.2022.10052084
  * Dempster et al. (1977) Dempster, A. P., Laird, N. M., & Rubin, D. B. (1977).  Maximum Likelihood from Incomplete Data via the EM Algorithm.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , _39_(1), 1–38.
  * Eden et al. (2004) Eden, U. T., Frank, L. M., Barbieri, R., Solo, V., & Brown, E. N. (2004).  Dynamic Analysis of Neural Encoding by Point Process Adaptive Filtering.  _Neural Computation_ , _16_(5), 971–998.  doi:10.1162/089976604773135069
  * Eden et al. (2018) Eden, U. T., Frank, L. M., & Tao, L. (2018).  Characterizing complex, multi-scale neural phenomena using state-space models.  In _Dynamic Neuroscience: Statistics, Modeling, and Control_ (pp. 29-52).
  * Elton and Gao (2015) Elton, A., & Gao, W. (2015).  Task‐related modulation of functional connectivity variability and its behavioral correlations.  _Human Brain Mapping_ , _36_(8), 3260–3272.
  * Fries et al. (2001) Fries, P., et al. (2001).  Modulation of oscillatory neuronal synchronization by selective visual attention.  _Science_ , _291_ , 1560 – 1563.  doi:10.1126/science.1055465
  * Fries (2005) Fries, P. (2005).  A mechanism for cognitive dynamics: neuronal communication through neuronal coherence.  _Trends in Cognitive Sciences_ , _9_ , 474 – 480.  doi:10.1016/j.tics.2005.08.011
  * Ghahramani and Hinton (1996) Ghahramani, Z., & Hinton, G. E. (1996).  Parameter estimation for linear dynamical systems.
  * Greenhouse et al. (1987) Greenhouse, J. B., Kass, R. E., & Tsay, R. S. (1987).  Fitting nonlinear models with ARMA errors to biological rhythm data.  _Statistics in Medicine_ , _6_(2), 167–183.
  * Hamilton (1990) Hamilton, James D. (1990).  Analysis of time series subject to changes in regime.  _Journal of Econometrics_ , _45_(1–2), 39–70.  doi:10.1016/0304-4076(90)90093-9
  * Hannan (1970) Hannan, E. J. (1970).  _Multiple Time Series_.  New York: Wiley.
  * Hirschmann et al. (2011) Hirschmann, J., Özkurt, T. E., Butz, M., Homburger, M., Elben, S., Hartmann, C. J., Vesper, J., Wojtecki, L., & Schnitzler, A. (2011).  Distinct oscillatory STN-cortical loops revealed by simultaneous MEG and local field potential recordings in patients with Parkinson’s disease.  _NeuroImage_ , _55_(3), 1159–1168.  doi:10.1016/j.neuroimage.2010.11.063
  * Hsin et al. (2022) Hsin, W.-C., Eden, U. T., & Stephen, E. P. (2022).  Switching Functional Network Models of Oscillatory Brain Dynamics.  In _2022 56th Asilomar Conference on Signals, Systems, and Computers_ (pp. 607-612), Pacific Grove, CA, USA.  doi:10.1109/IEEECONF56349.2022.10052077
  * Horwitz (2003) Horwitz, B. (2003).  The elusive concept of brain connectivity.  In _NeuroImage_ , _19_(2 Pt 1), 466–470.  doi:10.1016/s1053-8119(03)00112-5
  * Jarvis and Mitra (2001) Jarvis, M. R., & Mitra, P. P. (2001).  Sampling properties of the spectrum and coherency of sequences of action potentials.  _Neural Computation_ , _13_(4), 717–749.  doi:10.1162/089976601300014312
  * Kalman (1960) Kalman, R. E. (1960).  A New Approach to Linear Filtering and Prediction Problems.  _ASME. Journal of Basic Engineering_ , _82_(1), 35–45.  doi:10.1115/1.3662552
  * Kass et al. (2014) Kass, R. E., Eden, U. T., & Brown, E. N. (2014).  _Analysis of Neural Data_ (Vol. 491).  Springer.
  * Kim (1994) Kim, Chang-Jin (1994).  Dynamic linear models with Markov-switching.  _Journal of Econometrics_ , _60_(1–2), 1–22.  doi:10.1016/0304-4076(94)90036-1
  * Kucyi and Davis (2014) Kucyi, A., & Davis, K. D. (2014).  Dynamic functional connectivity of the default mode network tracks daydreaming.  _Neuroimage_ , _100_ , 471–480.
  * Kucyi et al. (2013) Kucyi, A., Salomons, T. V., & Davis, K. D. (2013).  Mind wandering away from pain dynamically engages antinociceptive and default mode brain networks.  _Proceedings of the National Academy of Sciences_ , _110_(46), 18692–18697.
  * Lauritzen (1996) Lauritzen, Steffen L. (1996).  _Graphical Models_.  Oxford University Press.
  * Little and Brown (2014) Little, S., & Brown, P. (2014).  The functional role of beta oscillations in Parkinson’s disease.  _Parkinsonism & Related Disorders_, _20 Suppl 1_ , S44–S48.  doi:10.1016/S1353-8020(13)70013-0
  * Litvak et al. (2011) Litvak, V., Jha, A., Eusebio, A., Oostenveld, R., Foltynie, T., Limousin, P., Zrinzo, L., Hariz, M. I., Friston, K., & Brown, P. (2011).  Resting oscillatory cortico-subthalamic connectivity in patients with Parkinson’s disease.  _Brain_ , _134_(2), 359–374.  doi:10.1093/brain/awq332
  * Li and Zhang (2000) Li, X. R., & Zhang, Y. (2000).  Multiple-model estimation with variable structure. V. Likely-model set algorithm.  _IEEE Transactions on Aerospace and Electronic Systems_ , _36_(2), 448-466.  doi:10.1109/7.845222
  * Lim and Mba (2015) Lim, C. K. R., & Mba, D. (2015).  Switching Kalman Filter for Failure Prognostic.  _Mechanical Systems and Signal Processing_ , _Volumes 52–53_ , 426 – 435.  doi:10.1016/j.ymssp.2014.08.006
  * Logothetis et al. (2001) Logothetis, N. K., Pauls, J., Augath, M., Trinath, T., & Oeltermann, A. (2001).  Neurophysiological investigation of the basis of the fMRI signal.  _Nature_ , _412_(6843), 150–157.
  * Manfredi et al. (2005) Manfredi, V., Mahadevan, S., & Kurose, J. (2005).  Switching Kalman Filters for Prediction and Tracking in an Adaptive Meteorological Sensing Network.  In _2005 Second Annual IEEE Communications Society Conference on Sensor and Ad Hoc Communications and Networks (IEEE SECON 2005)_ (pp. 197-206), Santa Clara, CA, USA.  doi:10.1109/SAHCN.2005.1557075
  * Matsuda and Komaki (2017) Matsuda, T., & Komaki, F. (2017).  Multivariate Time Series Decomposition into Oscillation Components.  _Neural Computation_ , _29_(8), 2055–2075.  doi:10.1162/NECO.a.00981
  * Matsuda and Komaki (2017) Matsuda, T., & Komaki, F. (2017).  Time Series Decomposition into Oscillation Components and Phase Estimation.  _Neural Computation_ , _29_(2), 332–367.
  * Matsuda et al. (2022) Matsuda, T., Homae, F., Watanabe, H., Taga, G., & Komaki, F. (2022).  Oscillator decomposition of infant fNIRS data.  _PLOS Computational Biology_ , _18_(3), e1009985.  doi:10.1371/journal.pcbi.1009985
  * McGregor and Nelson (2019) McGregor, M. M., & Nelson, A. B. (2019).  Circuit Mechanisms of Parkinson’s Disease.  _Neuron_ , _101_(6), 1042–1056.  doi:10.1016/j.neuron.2019.03.004
  * Mitra (2007) Mitra, P. (2007).  _Observed Brain Dynamics_.  Oxford University Press.
  * Mukamel et al. (2014) Mukamel, E. A., Pirondini, E., Babadi, B., Wong, K. F. K., Pierce, E. T., Harrell, P. G., Walsh, J. L., Salazar-Gomez, A. F., Cash, S. S., Eskandar, E. N., Weiner, V. S., Brown, E. N., & Purdon, P. L. (2014).  A Transition in Brain State during Propofol-Induced Unconsciousness.  _Journal of Neuroscience_ , _34_(3), 839–845.  doi:10.1523/JNEUROSCI.5813-12.2014
  * Murphy (1998) Murphy, K. P. (1998).  Switching Kalman filter.  Technical Report 98-10, Compaq Cambridge Research Laboratory, Cambridge, MA.
  * O’Neill et al. (2015) O’Neill, G. C., Barratt, E. L., Hunt, B. A., Tewarie, P. K., & Brookes, M. J. (2015).  Measuring electrophysiological connectivity by power envelope correlation: a technical review on MEG methods.  _Physics in Medicine and Biology_ , _60_(21), R271–R295.  doi:10.1088/0031-9155/60/21/R271
  * Paninski et al. (2010) Paninski, L., Ahmadian, Y., Ferreira, D. G., Koyama, S., Rahnama Rad, K., Vidne, M., et al. (2010).  A new look at state-space models for neural data.  _Journal of Computational Neuroscience_ , _29_ , 107–126.
  * Priestley (1981) Priestley, M.B. (1981).  _Spectral Analysis and Time Series_.  Academic Press Inc., New York.
  * Preti et al. (2017) Preti, M. G., Bolton, T. A. W., & Van De Ville, D. (2017).  The dynamic functional connectome: State-of-the-art and perspectives.  _NeuroImage_ , _160_ , 41–54.  doi:10.1016/j.neuroimage.2016.12.061
  * Purdon et al. (2013) Purdon, P. L., Pierce, E. T., Mukamel, E. A., Prerau, M. J., Walsh, J. L., Wong, K. F., Salazar-Gomez, A. F., Harrell, P. G., Sampson, A. L., Cimenser, A., Ching, S., Kopell, N. J., Tavares-Stoeckel, C., Habeeb, K., Merhar, R., & Brown, E. N. (2013).  Electroencephalogram signatures of loss and recovery of consciousness from propofol.  _Proceedings of the National Academy of Sciences of the United States of America_ , _110_(12), E1142–E1151.  doi:10.1073/pnas.1221180110
  * Rabiner (1989) Rabiner, L. R. (1989).  A tutorial on hidden Markov models and selected applications in speech recognition.  In _Proceedings of the IEEE_ , _77_(2), 257–286.  doi:10.1109/5.18626
  * Rauch et al. (1965) Rauch, H. E., Tung, F., & Striebel, C. T. (1965).  Maximum likelihood estimates of linear dynamic systems.  _AIAA Journal_ , _3_(8), 1445–1450.
  * Rogers et al. (2007) Rogers, B. P., Morgan, V. L., Newton, A. T., & Gore, J. C. (2007).  Assessing functional connectivity in the human brain by fMRI.  _Magnetic Resonance Imaging_ , _25_(10), 1347–1357.
  * Sarkka (2013) Särkkä, S. (2013).  _Bayesian Filtering and Smoothing_.  Cambridge University Press.
  * Schnitzler and Gross (2005) Schnitzler, A., & Gross, J. (2005).  Normal and pathological oscillatory communication in the brain.  _Nature Reviews Neuroscience_ , _6_ , 285–296.  doi:10.1038/nrn1650
  * Shumway and Stoffer (1991) Shumway, R. H., & Stoffer, D. S. (1991).  Dynamic Linear Models With Switching.  _Journal of the American Statistical Association_ , _86_(415), 763–769.  doi:10.2307/2290410
  * Singer (1999) Singer, W. (1999).  Neuronal synchrony: A versatile code for the definition of relations?  _Neuron_ , _24_ , 49 – 65, 111 – 125.  doi:10.1016/S0896-6273(00)80821-1
  * Smith and Brown (2003) Smith, A. C., & Brown, E. N. (2003).  Estimating a state-space model from point process observations.  _Neural Computation_ , _15_(5), 965–991.  doi:10.1162/089976603765202622
  * Soulat et al. (2022) Soulat, H., Stephen, E.P., Beck, A.M., et al. (2022).  State space methods for phase amplitude coupling analysis.  _Scientific Reports_ , _12_ , 15940.  doi:10.1038/s41598-022-18475-3
  * Stephen et al. (2020) Stephen, E.P., Hotan, G.C., Pierce, E.T., et al. (2020).  Broadband slow-wave modulation in posterior and anterior cortex tracks distinct states of propofol-induced unconsciousness.  _Scientific Reports_ , _10_ , 13701.  doi:10.1038/s41598-020-68756-y
  * Steriade et al. (1993) Steriade, M., McCormick, D. A., & Sejnowski, T. J. (1993).  Thalamocortical oscillations in the sleeping and aroused brain.  _Science_ , _262_(5134), 679–685.  doi:10.1126/science.8235588
  * Swann et al. (2016) Swann, N. C., De Hemptinne, C., Miocinovic, S., Qasim, S., Wang, S. S., Ziman, N., et al. (2016).  Gamma oscillations in the hyperkinetic state detected with chronic human brain recordings in Parkinson’s disease.  _Journal of Neuroscience_ , _36_(24), 6445–6458.
  * van Wijk et al. (2017) van Wijk, B. C. M., Neumann, W.-J., Schneider, G.-H., Sander, T. H., Litvak, V., & Kühn, A. A. (2017).  Low-beta cortico-pallidal coherence decreases during movement and correlates with overall reaction time.  _NeuroImage_ , _159_ , 1–8.  doi:10.1016/j.neuroimage.2017.07.024
  * Veeraraghavan et al. (2005) Veeraraghavan, H., Schrater, P., & Papanikolopoulos, N. (2005).  Switching Kalman Filter-Based Approach for Tracking and Event Detection at Traffic Intersections.  In _Proceedings of the 2005 IEEE International Symposium on Mediterranean Conference on Control and Automation Intelligent Control_ (pp. 1167-1172), Limassol, Cyprus.  doi:10.1109/.2005.1467180
  * Vinding et al. (2020) Vinding, M. C., Tsitsi, P., Waldthaler, J., Oostenveld, R., Ingvar, M., Svenningsson, P., & Lundqvist, D. (2020).  Reduction of spontaneous cortical beta bursts in Parkinson’s disease is linked to symptom severity.  _Brain Communications_ , _2_(1), fcaa052.  doi:10.1093/braincomms/fcaa052
  * Wang et al. (2018) Wang, D. D., de Hemptinne, C., Miocinovic, S., Ostrem, J. L., Galifianakis, N. B., San Luciano, M., & Starr, P. A. (2018).  Pallidal Deep-Brain Stimulation Disrupts Pallidal Beta Oscillations and Coherence with Primary Motor Cortex in Parkinson’s Disease.  _The Journal of Neuroscience : The Official Journal of the Society for Neuroscience_ , _38_(19), 4556–4568.  doi:10.1523/JNEUROSCI.0431-18.2018
  * Wodeyar et al. (2021) Wodeyar, A., Schatza, M., Widge, A. S., Eden, U. T., & Kramer, M. A. (2021).  A state space modeling approach to real-time phase estimation.  _eLife_ , _10_ , e68803.
  * Wong et al. (2011) Wong, K. F., Mukamel, E. A., Salazar, A. F., Pierce, E. T., Harrell, P. G., Walsh, J. L., Sampson, A., Brown, E. N., Purdon, P. L. (2011).  Robust time-varying multivariate coherence estimation: application to electroencephalogram recordings during general anesthesia.  In _Annual International Conference of the IEEE Engineering in Medicine and Biology Society_ (pp. 4725-4728).  doi:10.1109/IEMBS.2011.6091170  PMID: 22255393; PMCID: PMC3293406.
  * Xu and Jordan (1996) Xu, Lei, & Jordan, Michael I. (1996).  On Convergence Properties of the EM Algorithm for Gaussian Mixtures.  _Neural Computation_ , _8_(1), 129–151.  doi:10.1162/neco.1996.8.1.129