# Network inference in a stochastic multi-population neural mass model via
approximate Bayesian computation

Susanne Ditlevsen11footnotemark: 1, Massimiliano Tamborrino22footnotemark: 2,
Irene Tubikanec33footnotemark: 3
Department of Mathematical Sciences, University of Copenhagen
(susanne@math.ku.dk)Department of Statistics, University of Warwick
(massimiliano.tamborrino@warwick.ac.uk)Department of Statistics, University of
Klagenfurt<EMAIL_ADDRESS>
## Abstract

The aim of this article is to infer the connectivity structures of brain
regions before and during epileptic seizure. Our contributions are fourfold.
First, we propose a $6N$-dimensional stochastic differential equation for
modelling the activity of $N$ coupled populations of neurons in the brain.
This model further develops the (single population) stochastic Jansen and Rit
neural mass model, which describes human electroencephalography (EEG) rhythms,
in particular signals with epileptic activity. Second, we construct a reliable
and efficient numerical scheme for the model simulation, extending a splitting
procedure proposed for one neural population. Third, we propose an adapted
Sequential Monte Carlo Approximate Bayesian Computation algorithm for
simulation-based in-ference of both the relevant real-valued model parameters
as well as the $\\{0,1\\}$-valued network parameters, the latter describing
the coupling directions among the $N$ modelled neural populations. Fourth,
after illustrating and validating the proposed statistical approach on
different types of simulated data, we apply it to a set of multi-channel EEG
data recorded before and during an epileptic seizure. The real data
experiments suggest, for example, a larger activation in each neural
population and a stronger connectivity on the left brain hemisphere during
seizure.

#### Keywords

Jansen and Rit neural mass model, Stochastic Hamiltonian systems, Splitting
methods, Approximate Bayesian computation, Simulation-based inference, EEG
data

#### AMS subject classifications

60H10, 60H35, 65C30

#### Authorship contribution statement

S.D., M.T. and I.T. designed the study. I.T. performed the analysis, wrote the
code, and produced all figures. I.T. wrote the first draft. S.D., M.T. and
I.T. discussed the results and implications, commented on the manuscript at
all stages, edited the text and approved the final version.

#### Acknowledgements

S.D. was supported by the Novo Nordisk Foundation NNF20OC0062958 and the Independent Research Fund Denmark | Natural Sciences 9040-00215B. A part of this paper was written while I.T. was member of the Institute of Stochastics, Johannes Kepler University Linz, 4040 Linz, Austria. During this time, I.T. was supported by the Austrian Science Fund (FWF): W1214-N15, project DK14. The authors are also grateful for the generous mobility funding provided by the University of Klagenfurt.

##  1 Introduction

Estimating the connectivity in a network of units is important in a vast
variety of applications, ranging from biology over climate to finance,
physics, sociology and other fields. This is a statistically challenging task,
as these interacting units typically follow some complex underlying stochastic
dynamics, which may be only partially observed. Moreover, the detection of
directed connections and the distinction from spurious ones is particularly
difficult. This article is motivated by neuroscience, where the study of
neural activity and the underlying connectivity between brain regions is
essential for understanding brain function and its implications in various
neurological conditions. In particular, we are interested in inferring the
connectivity structure of the brain before and during epileptic seizure, a
critical area of study for understanding and managing epilepsy.
Electroencephalography (EEG) is a widely used technique to measure and analyze
brain activity, providing insights into the complex dynamics of the brain.
Here, we aim at estimating the directed connections between $N$ coupled neural
populations whose activity corresponds to the simultaneous measurements from
$N$ electrodes in the EEG recordings, one per population.

In this article, each neural population is modelled with a stochastic version
of the Jansen and Rit neural mass model (JR-NMM) originally proposed in [18].
This model has shown useful for reproducing human EEG rhythms, including those
associated with epileptic activity. The original JR-NMM is a $6$-dimensional
system of ordinary differential equations and describes the average activity
of one neural population, corresponding to the measured activity of one
electrode in EEG recordings. The model includes a term representing noisy
extrinsic input from the neighborhood or more distant regions. As this term
can be interpreted as a stochastic process, the solution of the dynamical
system is a stochastic process as well, inheriting the analytical properties
of this stochastic input function. Therefore, the model was reformulated as a
stochastic differential equation (SDE) in [1], and proved to be geometrically
ergodic, a property guaranteeing that the distribution of the solution
converges to a unique limit distribution, exponentially fast and for any
initial value. This has two important statistical implications: First, the
choice of the (unknown) initial value is negligible, since its impact on the
distribution of the process decreases exponentially fast. Second, quantities
related to the distribution of the process can be estimated from a
(sufficiently long) single path, instead of requiring many repeated path
simulations. While the original JR-NMM [18] has been extended to model $N$
coupled neural populations in [42], no such extension has been proposed for
the SDE version [1] of this model. We fill in this gap by proposing a
$6N$-dimensional SDE model accounting for directed connections between $N$
neural populations (first contribution). In contrast to [42], our model
contains $\\{0,1\\}$-valued coupling direction parameters, describing the
underlying network structure that we aim to infer from EEG data.

As neither the underlying transition density nor exact simulation schemes are
available for the $6N$-dimensional SDE model, a suitable numerical
approximation is required. The commonly applied Euler-Maruyama discretisation
is not suitable for this SDE, as it has been shown to fail in preserving
crucial structural properties for the single population model, such as the
dynamics of the modelled neural oscillations [1, 8]. The interested reader is
also referred to [7, 10, 11, 20, 37], where similar issues of the Euler-
Maruyama method applied to oscillatory SDE models are reported. Consequently,
when embedded into a statistical inference procedure, this standard numerical
scheme may yield wrong estimation results, make the inference algorithm
computationally infeasible or lead to ill-conditioned estimation methods [8,
14, 29]. Therefore, to obtain a reliable and efficient numerical simulation
method, we further develop the structure-preserving splitting procedure
proposed in [1] for $N=1$ population to our $N$-population SDE model (second
contribution). In contrast to the Euler-Maruyama scheme, which is based on
truncating a stochastic Taylor series, the idea behind the splitting approach
is to divide the unsolvable equation into solvable subequations, and to
compose their solutions in a proper way [5, 24]. The constructed splitting
scheme for the stochastic $N$-population neural mass model is based on a
Hamiltonian type re-formulation of the SDE and successfully handles multiple
interacting components, allowing to accurately simulate the complex dynamics
of coupled neural populations.

As a next step, we aim to estimate the $\\{0,1\\}$-valued coupling direction
parameters between the $N$ neural populations, as well as relevant real-valued
model parameters of the $6N$-dimensional SDE from $N$ simultaneously recorded
EEG signals. This is particularly challenging, since this SDE falls into the
class of hypoelliptic models [14, 21, 25, 29] and is only partially observed
via $N$ one-dimensional linear functions of a subset of the $6N$ model
components. These issues, combined with the lack of a tractable underlying
likelihood, make this problem naturally suitable for likelihood-free inference
approaches. We focus on the simulation-based Approximate Bayesian Computation
(ABC) method [34], which has become one of the leading tools for parameter
estimation in complex mathematical models in the last decades. The idea behind
the classical acceptance-rejection ABC algorithm is to: (a) sample a parameter
candidate from a prior distribution; (b) simulate a dataset conditioned on
this value; (c) keep the candidate as a sample from the desired approximate
posterior if the distance between the summary statistics of the simulated and
observed datasets is smaller than some threshold level [4, 23, 34]. This basic
algorithm is computationally expensive, because parameter candidates are
sampled from the prior distribution throughout, whose mass typically
concentrates in regions “far away” from those of the posterior distribution,
leading to low acceptance rates and high computational costs. To tackle this,
here we consider the Sequential Monte Carlo (SMC) ABC approach, which
represents the state-of-the-art sampler within ABC. This sequential algorithm
relies on targeted proposal samplers (constructed across iterations) that
efficiently avoid improbable parameter regions, yielding intermediate
approximate posterior distributions that move towards the desired posterior
during several iterations [3, 13, 35]. In particular, we propose an adapted
SMC-ABC algorithm to infer both the $\\{0,1\\}$-valued network parameters and
the real-valued model parameters (third contribution). The algorithm builds
upon the work of [8] for Hamiltonian-type SDEs observed via univariate time
series, further developing it at the following three levels. First, we rely on
the numerical splitting scheme for synthetic data generation from the
$6N$-dimensional SDE. Second, we construct effective summary statistics by
mapping the $N$-dimensional time series to their estimated $N$ marginal
densities and spectral densities, as well as to their estimated cross-spectral
densities and cross-correlation functions to account for possible dependencies
among the populations. Third and foremost, to handle the high-dimensional
parameter space and the presence of $\\{0,1\\}$-valued parameters, we consider
two independent kernels for the SMC-ABC proposal, one Gaussian for the real-
valued and one Bernoulli for the $\\{0,1\\}$-valued parameters. By estimating
the latter, we uncover the directed connectivity structure of the neural
populations, enabling the identification of the intricate network connections
within brain regions.

After validating the performance of the proposed statistical approach on
different types of simulated data, we apply it on real multi-channel EEG data
measured before and during epileptic seizure (fourth contribution). We
successfully obtain unimodal posteriors of a relatively large number of
continuous model parameters and a clear network estimate during seizure.
Comparing this network with that obtained based on the EEG recordings before
seizure, we observe, e.g., a stronger connectivity on the left brain
hemisphere during seizure.

The proposed algorithm, along with the choice of its key ingredients, holds
promise beyond the specific application, as it can be used for other coupled
ergodic SDE models, provided that a reliable numerical simulation method can
be derived. This flexibility broadens the scope of our work, offering a
potential avenue for parameter and network estimation in various complex
stochastic systems involving coupled units.

The paper is organised as follows. In Section 2, we introduce the model. In
Section 3, we construct the numerical splitting method. In Section 4, we
detail the proposed adapted SMC-ABC algorithm for network inference. In
Section 5, we illustrate its performance on simulated data. In Section 6, we
apply our method to real EEG data with epileptic activity. Conclusions and
discussion are reported in Section 7. Sample code is available at
https://github.com/IreneTubikanec/networkABC.

##  2 Multi-population stochastic Jansen and Rit neural mass model

We first present the stochastic JR-NMM of one neural population [1]. We then
extend it to a system of multiple coupled neural populations, following the
strategies proposed in [19, 42]. Finally, inspired by [1], we derive the
corresponding Hamiltonian type formulation of the multi-population model.

######  Notation.

We denote by $0_{d}$ the $d$-dimensional zero vector, by $\mathbb{O}_{d}$ the
$d\times d$-dimensional zero matrix, by $\mathbb{I}_{d}$ the $d\times
d$-dimensional identity matrix, and by $\textrm{diag}[a_{1},\ldots,a_{d}]$ a
$d\times d$-dimensional diagonal matrix with diagonal entries
$a_{1},\ldots,a_{d}$. The transpose is denoted by ⊤ and the Euclidean norm by
$\left\lVert\cdot\right\rVert$. We sometimes omit the time index of a
stochastic process and use, e.g., $(X(t))_{t\in[0,T]}$ and $X$
interchangeably.

###  2.1 Modelling one neural population

Let $[0,T]$, with $T>0$, be the time interval of interest. Let
$(\Omega,\mathcal{F},\mathbb{P})$ be a complete probability space with a
complete and right-continuous filtration $(\mathcal{F}(t))_{t\in[0,T]}$, and
let $(W_{i}(t))_{t\in[0,T]}$, $i=4,5,6$, be independent Wiener processes on
$(\Omega,\mathcal{F},\mathbb{P})$ and adapted to
$(\mathcal{F}(t))_{t\in[0,T]}$.

Let $X_{i}$, $i=1,2,3$, model the average postsynaptic potentials of three
groups of neurons, i.e., the main neurons, excitatory interneurons and
inhibitory interneurons, respectively. Define the sigmoid function
$\textrm{sig}\colon\mathbb{R}\to[0,\nu_{\max}]$, with $\nu_{\max}>0$, by

$\textrm{sig}(x):=\frac{\nu_{\max}}{1+e^{\gamma(v_{0}-x)}},$

where $x$ is a potential, $\nu_{\max}$ describes the maximum firing rate of
the neural population, $v_{0}\in\mathbb{R}$ is the potential value for which
$50\%$ of the maximum firing rate is attained and $\gamma>0$ is proportional
to the slope of the sigmoid function at $v_{0}$. The sigmoid function is used
to transform the average membrane potential of a neural subpopulation into an
average firing rate. The stochastic version of the JR-NMM [1] is given by

$\displaystyle\begin{split}dX_{1}(t)&=X_{4}(t)dt\\\ dX_{2}(t)&=X_{5}(t)dt\\\
dX_{3}(t)&=X_{6}(t)dt\\\
dX_{4}(t)&=\left[Aa\bigl{(}\textrm{sig}\left(X_{2}(t)-X_{3}(t)\right)\bigr{)}-2aX_{4}(t)-a^{2}X_{1}(t)\right]dt+\bar{\epsilon}dW_{4}(t)\\\
dX_{5}(t)&=\left[Aa\bigl{(}\mu+C_{2}\textrm{sig}\bigl{(}C_{1}X_{1}(t)\bigr{)}\bigr{)}-2aX_{5}(t)-a^{2}X_{2}(t)\right]dt+\sigma
dW_{5}(t)\\\
dX_{6}(t)&=\left[BbC_{4}\textrm{sig}\left(C_{3}X_{1}(t)\right)-2bX_{6}(t)-b^{2}X_{3}(t)\right]dt+\tilde{\epsilon}dW_{6}(t),\end{split}$
(1)

with $\mathcal{F}(0)$-measurable initial value
$X_{0}=(X_{1}(0),\ldots,X_{6}(0))^{\top}$ which is independent of
$W=(W_{4},W_{5},W_{6})^{\top}$ and satisfies $\mathbb{E}\bigl{[}\left\lVert
X_{0}\right\rVert^{2}\bigr{]}<\infty$. The meaning and typical values of the
model parameters $A$, $B$, $a$, $b$, $C_{1}$, $C_{2}$, $C_{3}$, $C_{4}$,
$\nu_{\max}$, $\gamma$ and $v_{0}$ are reported in Table 1 (see also [1, 18,
42]). The parameters $\mu$ and $\sigma$ scale the deterministic and stochastic
input, respectively, coming from neighboring regions in the brain. Together
with the Wiener process $W_{5}$, they replace the stochastic input function of
the original model. While usually $\sigma\gg 1$, weak noise acts on the
components $X_{4}$ and $X_{6}$ with noise intensities
$\bar{\epsilon},\tilde{\epsilon}\ll\sigma$. Throughout, we set
$\epsilon=\bar{\epsilon}=\tilde{\epsilon}=1\ s^{-1}$, since their role is
marginal. See [1] for further details.

The stochastic JR-NMM (1) is an additive noise SDE with globally Lipschitz
drift coefficient, and thus has a pathwise unique solution which is adapted to
$(\mathcal{F}(t))_{t\in[0,T]}$ [2, 22]. Moreover, system (1) is hypoelliptic
and geometrically ergodic [1]. The $6$-dimensional solution
${X}=({X_{1}},{X_{2}},{X_{3}},{X_{4}},{X_{5}},{X_{6}})^{\top}$ is only
partially observed through the difference of two of its coordinates

$Y(t):=X_{2}(t)-X_{3}(t),\quad t\in[0,T].$ (2)

The process $Y=(Y(t))_{t\in[0,T]}$ is the observed output of the model and
describes the average membrane potential of the main neurons as measured with
EEG. An illustration of a simulated trace with parameters chosen to produce
$\alpha$-waves (neural oscillations in the 8–12 Hz frequency band) is provided
in Figure 1, Panel A. The model can also produce more complex behaviour, such
as brain signals occurring before and during epileptic seizures. This is
achieved by increasing the excitation-inhibition-ratio $A/B$ [42]. In Figure 1
Panel B, regular EEG activity is produced by using $\mu=90$, $\sigma=500$,
$A=3.25$ and the standard values reported in Table 1. Increasing the
excitation-inhibition-ratio by increasing $A$ causes the model to generate
sporadic and frequently occurring spikes ($A=3.5$ and $A=3.6$) and rhythmic
discharge of spikes ($A=4.3$), as observed in Panels C-E, respectively.

Figure 1: Traces of one neural population with different types of activity.
Simulated paths of the process $Y$ in eq. (2) of model (1). Panel A: Trace
showing $\alpha$-rhythmic activity using the standard values of Table 1,
$A=3.25,C=134.263,\mu=202.547$, and $\sigma=1859.211$ (values taken from [8]).
Panels B-E: Traces describing activity occurring during epileptic seizures
using the standard values of Table 1, $\mu=90,\sigma=500$ and different values
of $A$. Panel B: Regular EEG, $A=3.25$. Panel C: Sporadic spikes, $A=3.5$.
Panel D: Frequently occurring spikes, $A=3.6$. Panel E: Rhythmic discharge of
spikes, $A=4.3$. The signals resemble experimental stereo EEG recordings (cf.
Figure 3 in [42]).  Table 1: Standard parameter values for the Jansen and Rit
Neural Mass Model [1, 18, 42].

Parameter | Meaning | Standard value
---|---|---
$A$ | Average excitatory synaptic gain |  $3.25$ mV
$B$ | Average inhibitory synaptic gain |  $22$ mV
$a$ | Membrane time constant of excitatory postsynaptic potential |  $100$ s-1
$b$ | Membrane time constant of inhibitory postsynaptic potential |  $50$ s-1
$C$ | Average number of synapses between the subpopulations | $135$
$C_{1}$, $C_{2}$ | Avg. no. of synaptic contacts in the excitatory feedback loop |  $C$, $0.8\ C$
$C_{3}$, $C_{4}$ | Avg. no. of synaptic contacts in the inhibitory feedback loop |  $0.25\ C$, $0.25\ C$
$\nu_{\max}$ | Maximum firing rate (Maximum of the sigmoid function) |  $5$ s-1
$v_{0}$ |  Value for which $50\%$ of the maximum firing rate is attained |  $6$ mV
$\gamma$ |  Determines the slope of the sigmoid function at $v_{0}$ |  $0.56$ mV-1

###  2.2 Modelling multiple coupled neural populations

Here, we define the $6N$-dimensional SDE describing a system of $N>1$ neural
populations, where the $k$-th population $X^{k}=(X_{1}^{k},\ldots,X_{6}^{k})$,
$k=1,\ldots,N$, satisfies (1) with suitable index $k$, except for the fifth
equation, which is instead given by

$dX^{k}_{5}(t)=\Bigl{[}A_{k}a_{k}\Bigl{(}\mu_{k}+C_{2,k}\textrm{sig}\left(C_{1,k}X^{k}_{1}(t)\right)+\sum\limits_{j=1,j\neq
k}^{N}\rho_{jk}K_{jk}X_{1}^{j}(t)\Bigr{)}-2a_{k}X^{k}_{5}(t)-a_{k}^{2}X^{k}_{2}(t)\Bigr{]}dt+\sigma_{k}dW^{k}_{5}(t).$

In particular, the $k$-th population $X^{k}$ receives inputs from populations
$j\neq k$, $j=1,\ldots,N$, through its coordinate $X^{k}_{5}$, where
$K_{jk}>0$ model the coupling strengths and $\rho_{jk}\in\\{0,1\\}$ determine
whether there is ($\rho_{jk}=1$) or not ($\rho_{jk}=0$) a directed coupling
from the $j$-th to the $k$-th population, defining thus the functional network
of relations among the neural populations. Since the main pyramidal cells are
excitatory and have axons reaching other brain areas, the average action
potential $X_{1}^{j}(t)$ of population $j\neq k$ provides excitatory input to
population $k$, and thus $K_{jk}>0$. Note that the coupling direction
parameters $\rho_{jk}$ where not considered in the modelling approach
discussed in [42].

The $N$-dimensional observed output is

$Y(t):=(Y^{1}(t),\ldots,Y^{N}(t))^{\top}=(X^{1}_{2}(t)-X^{1}_{3}(t),\ldots,X^{N}_{2}(t)-X^{N}_{3}(t))^{\top},\quad
t\in[0,T].$ (3)

Not only the excitation-inhibition-ratio $A/B$ is relevant for epileptic
behaviour, but also the coupling strengths and directions between neural
groups play a crucial role [42]. This is illustrated in Figure 2, where
simulated activity of $N=4$ neural populations under different coupling
regimes is shown. In the left panels, no coupling occurs, i.e., all
$\rho_{jk}$ are set to zero. In the middle and right panels, there is a
unidirectional (cascade) coupling structure (illustrated in Figure 3a), i.e,
$\rho_{12}=\rho_{23}=\rho_{34}=1$, for different coupling strengths $K$. The
activity of a passive site (no epileptic spikes occur without input from other
populations) strongly depends on that of an active site (epileptic spikes
occur without input from other populations). In Population $1$, the
excitation-inhibition-ratio is increased by setting $A_{1}=3.6$ to obtain
spiking activity. In the remaining populations, the standard values of Table
1, $\mu=90$ and $\sigma=500$ are used for each component. When there is no
coupling among the populations, no activation of Populations 2–4 occurs (left
panels). Introducing coupling and setting the coupling strength parameters to
$K_{12}=K_{23}=K_{34}=300$ (central panels) leads to a dependence of the
activity of Populations 2–4 on that of Population $1$. When the coupling is
strong enough ($K_{12}=K_{23}=K_{34}=500$, right panels), rhythmic
synchronization occurs. A similar behaviour for two populations has also been
observed for the original JR-NMM (cf. Figure 6 in [42]).

Figure 2: Cascade network of four populations with one active population.
Simulated paths of $Y$ in eq. (3) of $N=4$ neural populations. The standard
values of Table 1, $\mu=90$ and $\sigma=500$ are used, except for Population
$1$ where $A_{1}=3.6$ to make it active. Panel A: $\rho_{jk}=0$,
$j,k=1,\ldots,4$, $j\neq k$. Panels B-C: $\rho_{12}$, $\rho_{23},\rho_{34}=1$.
The coupling strength parameters $K_{12}$, $K_{23}$ and $K_{34}$ equal $300$
(Panel B) and $500$ (Panel C).

###  2.3 Formulation as a stochastic Hamiltonian type system

The multi-population JR-NMM can be formulated as a damped stochastic
Hamiltonian system with non-linear displacement, similarly to [1] for $N=1$.
This allows for a more compact formulation and constitutes the basis for the
derivation of the numerical method in Section 3.

First, let us decompose the $k$-th process as $X^{k}:=(Q^{k},P^{k})^{\top}$
with $3$-dimensional components $Q^{k}=(X_{1}^{k},X_{2}^{k},X_{3}^{k})^{\top}$
and $P^{k}=(X_{4}^{k},X_{5}^{k},X_{6}^{k})^{\top}$. Then, let us denote the
$3$-dimensional Wiener process in the $k$-th population by
$W^{k}=(W_{4}^{k},W_{5}^{k},W_{6}^{k})^{\top}$ and the corresponding $3\times
3$-dimensional diffusion matrix by
$\Sigma_{k}=\text{diag}[\epsilon_{k},\sigma_{k},\epsilon_{k}]$. The
Hamiltonian type formulation of the $k$-th population is

$d\begin{pmatrix}Q^{k}(t)\\\
P^{k}(t)\end{pmatrix}=\begin{pmatrix}\nabla_{P}H{{}_{k}}\bigl{(}Q^{k}(t),P^{k}(t)\bigr{)}\\\
-\nabla_{Q}H{{}_{k}}\bigl{(}Q^{k}(t),P^{k}(t)\bigr{)}-2\Gamma_{k}P^{k}(t)+G_{k}({Q(t)})\end{pmatrix}dt\
+\ \begin{pmatrix}\mathbb{O}_{3}\\\ \Sigma_{k}\end{pmatrix}dW^{k}(t).$ (4)

System (4) consists of a Hamiltonian part defined by the Hamiltonian function
$H{{}_{k}}:\mathbb{R}^{6}\to\mathbb{R}_{0}^{+}$ given by

$H{{}_{k}}(Q^{k},P^{k}):=\frac{1}{2}\left(\left\lVert
P^{k}\right\rVert^{2}+\left\lVert\Gamma_{k}Q^{k}\right\rVert^{2}\right),$

with gradients $\nabla_{P}H{{}_{k}}\bigl{(}Q^{k}(t),P^{k}(t)\bigr{)}=P^{k}(t)$
and
$\nabla_{Q}H{{}_{k}}\bigl{(}Q^{k}(t),P^{k}(t)\bigr{)}=\Gamma_{k}^{2}Q^{k}(t)$,
a damping part determined by the $3\times 3$-dimensional diagonal matrix
$\Gamma_{k}=\text{diag}[a_{k},a_{k},b_{k}]$, and a non-linear displacement and
coupling term $G_{k}:\mathbb{R}^{3{N}}\to\mathbb{R}^{3}$ given by

$G_{k}({Q(t)})=\begin{pmatrix}A_{k}a_{k}\textrm{sig}\left(X^{k}_{2}(t)-X^{k}_{3}(t)\right)\\\
A_{k}a_{k}\Bigl{(}\mu_{k}+C_{2,k}\textrm{sig}\left(C_{1,k}X_{1}^{k}(t)\right)+\sum\limits_{j=1,j\neq
k}^{N}\rho_{jk}K_{jk}X_{1}^{j}(t)\Bigr{)}\\\
B_{k}b_{k}C_{4,k}\textrm{sig}(C_{3,k}X_{1}^{k}(t))\end{pmatrix},$

where
$Q=(Q^{1},\ldots,Q^{N})^{\top}=(X_{1}^{1},X_{2}^{1},X_{3}^{1},\ldots,X_{1}^{N},X_{2}^{N},X_{3}^{N})^{\top}$.
A compact Hamiltonian type formulation for the $N$ coupled populations is then
obtained by the process $X:=(Q,P)^{\top}$ with $Q$ as above, $3N$-dimensional
component
$P=(P^{1},\ldots,P^{N})^{\top}=(X_{4}^{1},X_{5}^{1},X_{6}^{1},\ldots,X_{4}^{N},X_{5}^{N},X_{6}^{N})^{\top}$
and $3N$-dimensional Wiener process
$W=(W^{1},\ldots,W^{N})^{\top}=(W_{4}^{1},W_{5}^{1},W_{6}^{1},\ldots,W_{4}^{N},W_{5}^{N},W_{6}^{N})^{\top}$.
In particular, the whole $6N$-dimensional stochastic Hamiltonian system is
given by

$d\begin{pmatrix}Q(t)\\\ P(t)\end{pmatrix}=\begin{pmatrix}P(t)\\\
-\Gamma^{2}Q(t)-2\Gamma P(t)+G(Q(t))\end{pmatrix}dt\ +\
\begin{pmatrix}\mathbb{O}_{3N}\\\ \Sigma\end{pmatrix}dW(t),$ (5)

where $\Gamma=\text{diag}[a_{1},a_{1},b_{1},\ldots,a_{N},a_{N},b_{N}]$ and
$\Sigma=\text{diag}[\epsilon_{1},\sigma_{1},\epsilon_{1},\ldots,\epsilon_{N},\sigma_{N},\epsilon_{N}]$
are $3N\times 3N$-dimensional diagonal matrices, and
$G:\mathbb{R}^{3N}\to\mathbb{R}^{3N}$ is given by
$G(Q)=(G_{1}({Q}),\ldots,G_{N}({Q}))^{\top}$.

######  Remark 2.1.

Including the coupling terms into the displacement function
$G:\mathbb{R}^{3N}\to\mathbb{R}^{3N}$ enables closed-form expressions of all
required components of the splitting framework discussed in Section 3 for
arbitrary $N$.

##  3 Simulation of the multi-population JR-NMM: Splitting integrators

Let $0=t_{0}<\ldots<t_{m}=T$ be a partition of the time interval $[0,T]$ with
equidistant time steps $\Delta=t_{i+1}-t_{i}$ for $i=0,\ldots,m-1$,
$m\in\mathbb{N}$. As the solution $X(t_{i})$ of system (5) at the discrete
time points $t_{i}$ cannot be simulated exactly, we derive reliable
approximations $\widetilde{X}(t_{i})$ of $X(t_{i})$, extending the numerical
splitting method proposed in [1] for single neural populations to the multi-
population case.

The splitting approach consists of the following three steps [6, 24]:

  * (i)

Split the equation for $X(t)$ into $d\in\mathbb{N}$ exactly solvable
subequations for $X^{[l]}(t)$, $l=1,\ldots,d$;

  * (ii)

Derive the exact solutions (flows) of the $d$ subequations over an increment
of length $\Delta$, i.e., for initial value $X^{[l]}(t_{i})$, the
$\Delta$-flow is $\varphi_{\Delta}^{[l]}(X^{[l]}(t_{i}))=X^{[l]}(t_{i+1})$;

  * (iii)

Compose the $d$ exact solutions in a suitable way. Prominent methods are the
Lie-Trotter [40] and Strang [36] compositions, given by

$\displaystyle\widetilde{X}^{\textrm{LT}}(t_{i+1})$
$\displaystyle=\left(\varphi_{\Delta}^{[1]}\circ\ldots\circ\varphi_{\Delta}^{[d]}\right)\left(\widetilde{X}^{\textrm{LT}}(t_{i})\right),$
$\displaystyle\widetilde{X}^{\textrm{S}}(t_{i+1})$
$\displaystyle=\left(\varphi_{\Delta/2}^{[1]}\circ\ldots\circ\varphi_{\Delta/2}^{[d-1]}\circ\varphi_{\Delta}^{[d]}\circ\varphi_{\Delta/2}^{[d-1]}\circ\ldots\circ\varphi_{\Delta/2}^{[1]}\right)\left(\widetilde{X}^{\textrm{S}}(t_{i})\right),$

respectively. Note that the order of the compositions can be changed.

Despite having the same strong order of convergence, Strang compositions have
been shown to outperform Lie-Trotter methods, yielding a better approximation
of the true solution and preserving the model dynamics for larger time steps
[1, 7, 9, 10, 17, 28, 41]. This is possibly due to the symmetry of the Strang
splitting [1], its smaller mean biases [41] and its higher-order one-step
predictions [28]. We therefore apply the Strang splitting for system (5),
generalizing the Ornstein-Uhlenbeck integrator presented in [1]. Their Wiener
integrator can be extended analogously.

##### Step (i): Choice of subequations

We separate the non-linear term $G(Q(t))$ of system (5), and consider the
subequations

$\displaystyle d\begin{pmatrix}Q^{[1]}(t)\\\ P^{[1]}(t)\end{pmatrix}$
$\displaystyle=\begin{pmatrix}P^{[1]}(t)\\\ -\Gamma^{2}Q^{[1]}(t)-2\Gamma
P^{[1]}(t)\end{pmatrix}dt\ +\ \begin{pmatrix}\mathbb{O}_{3N}\\\
\Sigma\end{pmatrix}dW(t),$ (6) $\displaystyle d\begin{pmatrix}Q^{[2]}(t)\\\
P^{[2]}(t)\end{pmatrix}$ $\displaystyle=\begin{pmatrix}0_{3N}\\\
G(Q^{[2]}(t))\end{pmatrix}dt.$ (7)

##### Step (ii): Derivation of exact solutions

Write subequation (6) as

$dX^{[1]}(t)=FX^{[1]}(t)dt+\Sigma_{0}dW(t),$

where

$F=\begin{pmatrix}\mathbb{O}_{3N}&\mathbb{I}_{3N}\\\
-\Gamma^{2}&-2\Gamma\end{pmatrix},\quad\Sigma_{0}=\begin{pmatrix}\mathbb{O}_{3N}\\\
\Sigma\end{pmatrix}.$

Let $X^{[1]}(t_{i})=(Q^{[1]}(t_{i}),P^{[1]}(t_{i}))^{\top}$ denote the
solution of system (6) at time $t_{i}$. The exact solution at time $t_{i+1}$
is then given by

$X^{[1]}(t_{i+1})=\varphi_{\Delta}^{[1]}\left(X^{[1]}(t_{i})\right)=e^{F\Delta}X^{[1]}(t_{i})+\xi_{i}(\Delta),$

where $\xi_{i}(\Delta)$, $i=1,\ldots,m$, are independent $6N$-dimensional
Gaussian random vectors with zero mean $\mathbb{E}[\xi_{i}(\Delta)]=0_{6N}$
and covariance matrix given by

$\textrm{Cov}(\Delta)=\int\limits_{0}^{\Delta}e^{F(\Delta-s)}\Sigma_{0}\Sigma_{0}^{\top}\left(e^{F(\Delta-s)}\right)^{\top}ds.$

For general matrices $F$, the exponential matrix $e^{F\Delta}$ and the
covariance matrix $\textrm{Cov}(\Delta)$ may be costly to compute. However,
due to the sparseness of $F$ coming from the damped Hamiltonian term, both
$e^{F\Delta}$ and $\textrm{Cov}(\Delta)$ can be expressed in closed-form for
arbitrary $N$ as

$e^{F\Delta}=\begin{pmatrix}e^{-\Gamma\Delta}\left(\mathbb{I}_{3N}+\Gamma\Delta\right)&e^{-\Gamma\Delta}\Delta\\\
-\Gamma^{2}e^{-\Gamma\Delta}\Delta&e^{-\Gamma\Delta}\left(\mathbb{I}_{3N}-\Gamma\Delta\right)\end{pmatrix}=:\begin{pmatrix}\vartheta(\Delta)&\kappa(\Delta)\\\
\vartheta^{\prime}(\Delta)&\kappa^{\prime}(\Delta)\end{pmatrix}$

and

$\textrm{Cov}(\Delta)=\begin{pmatrix}\frac{1}{4}\Gamma^{-3}\Sigma^{2}\left(\mathbb{I}_{3N}+\kappa(\Delta)\vartheta^{\prime}(\Delta)-\vartheta^{2}(\Delta)\right)&\frac{1}{2}\Sigma^{2}\kappa^{2}(\Delta)\\\
\frac{1}{2}\Sigma^{2}\kappa^{2}(\Delta)&\frac{1}{4}\Gamma^{-1}\Sigma^{2}\left(\mathbb{I}_{3N}+\kappa(\Delta)\vartheta^{\prime}(\Delta)-\kappa^{\prime
2}(\Delta)\right)\end{pmatrix}.$

Denote the solution of the second subequation (7) at time $t_{i}$ by
$X^{[2]}(t_{i})=(Q^{[2]}(t_{i}),P^{[2]}(t_{i}))^{\top}$. Since the first
component of the right side of system (7) is zero and the second component
only depends on $Q$, the exact solution at time $t_{i+1}$ is given by

$X^{[2]}(t_{i+1})=\varphi_{\Delta}^{[2]}\left(X^{[2]}(t_{i})\right)=X^{[2]}(t_{i})+\Delta\begin{pmatrix}0_{3N}\\\
G(Q^{[2]}(t_{i}))\end{pmatrix}.$

##### Step (iii): Composition of exact solutions

The Strang splitting integrator for (5) is then given by

$\widetilde{X}^{\textrm{S}}(t_{i+1})=\left(\varphi_{\Delta/2}^{[2]}\circ\varphi_{\Delta}^{[1]}\circ\varphi_{\Delta/2}^{[2]}\right)\left(\widetilde{X}^{\textrm{S}}(t_{i})\right).$
(8)

Using (8), a path of system (5) can be simulated by applying Algorithm 1. An
analogous analysis as in [1] yields that the derived splitting scheme is mean-
square convergent of order $1$ and preserves the qualitative behaviour of the
solution of (5), e.g., amplitudes of oscillations, marginal invariant
densities and spectral densities.

Algorithm 1 Strang splitting scheme for the $N$-population stochastic JR-NMM
Input: Initial value $X_{0}$, step size $\Delta$, number of time steps $m$ in
$[0,T]$ and model parameters
Output: Approximated path of $(X(t))_{t\in[0,T]}$ at discrete times
$t_{i}=i\Delta$, $i=0,\ldots,m$, $t_{m}=T$.

1:Set $\widetilde{X}^{\textrm{S}}(t_{0})=X_{0}$

2:for $i=0:({m-1})$ do

3: Set
$X^{\textrm{[2]}}=\widetilde{X}^{\textrm{S}}(t_{i})+\frac{\Delta}{2}\begin{pmatrix}0_{3N}\\\
G(\widetilde{Q}^{\textrm{S}}(t_{i}))\end{pmatrix}$

4: Set $X^{\textrm{[1]}}=e^{F\Delta}X^{[2]}+\xi_{i}(\Delta)$

5: Set
$\widetilde{X}^{S}(t_{i+1})=X^{[1]}+\frac{\Delta}{2}\begin{pmatrix}0_{3N}\\\
G(Q^{[1]})\end{pmatrix}$

6:end for

7:Return $\widetilde{X}^{\textrm{S}}(t_{i})$, $i=0,\ldots,m$.

######  Remark 3.1.

The Strang splitting integrator

$\widetilde{X}^{\textrm{S2}}(t_{i+1})=\left(\varphi_{\Delta/2}^{[1]}\circ\varphi_{\Delta}^{[2]}\circ\varphi_{\Delta/2}^{[1]}\right)\left(\widetilde{X}^{\textrm{S2}}(t_{i})\right)$

is another possible choice. However, this composition requires to evaluate the
more costly stochastic subsystem twice at every iteration step, generating
twice as many pseudo-random numbers.

##  4 Adjusted SMC-ABC for network inference

In this section, we first revise the acceptance-rejection, reference-table
acceptance-rejection and SMC-ABC schemes. Then, as the considered parameters
of interest are both real and binary, we propose an adjusted SMC-ABC scheme,
where the Gaussian proposal sampler within the “canonical” SMC-ABC algorithm
is modified to account for both continuous (real-valued) model parameters and
discrete ($\\{0,1\\}$-valued) network parameters. Finally, we describe all
underlying key ingredients (numerical scheme for data generation, summary
statistics, proposal samplers, etc.). The resulting algorithm further develops
the spectral density-based and measure-preserving ABC method proposed in [8]
for ergodic SDEs with a univariate observed output process to multi-
dimensional SDEs with coupled ergodic components, partially observed via $N$
simultaneously recorded univariate time series. On the one hand, we replace
the initially proposed reference table acceptance-rejection ABC with the more
computationally efficient adjusted SMC-ABC (cf. Section 4.2). On the other
hand, we generalise the summaries and distances to cope with the $N>1$
simultaneously observed coupled time series (cf. Section 4.3).

###  4.1 ABC schemes

Let $\theta=(\theta_{1},\ldots,\theta_{d})$ denote the parameter vector to be
inferred from the observed data $y$. Denoting by $\pi(\theta)$ the prior and
by $\mathcal{L}(y|\theta)$ the likelihood function, the posterior distribution
$\pi(\theta|y)$ satisfies

$\pi(\theta|y)\propto\pi(\theta)\mathcal{L}(y|\theta).$

In general, the likelihood function, and thus the true posterior, are not
available for SDEs like (5). The idea of ABC is to replace the likelihood via
a large amount of “synthetic” datasets simulated from the model, obtaining an
approximate posterior $\pi_{\textrm{ABC}}(\theta|y)$ targeting the true
posterior $\pi(\theta|y)$.

Several ABC algorithms have been proposed throughout the years, see [34] for
an overview. Among all, the simplest is acceptance-rejection ABC [4, 23, 34],
consisting of three steps: (a) Sample $\theta^{{}^{\prime}}$ from the prior
$\pi(\theta)$; (b) Conditioned on $\theta^{{}^{\prime}}$, simulate a synthetic
dataset $\tilde{y}_{\theta^{{}^{\prime}}}$ (here a path) from the observed
output $Y$; (c) Keep the sampled value $\theta^{{}^{\prime}}$ if the distance
$D(\cdot,\cdot)$ between a vector of summary statistics $s(\cdot)$ of the
observed and simulated data is smaller than a threshold $\delta>0$, i.e.,
$D(s(y),s(\tilde{y}_{\theta^{\prime}}))<\delta$. Steps (a)-(c) are then
repeated until $M$ draws are accepted, which typically happens after $n\gg M$
drawings. This leads to the approximate acceptance-rejection ABC posterior

$\pi(\theta|y)\approx\pi_{\textrm{ABC}}^{\delta}(\theta|s(y))\propto\int\mathbbm{1}_{\\{D(s(y),s(\tilde{y}_{\theta^{\prime}}))<\delta\\}}\pi(\theta)\mathcal{L}(s(\tilde{y}_{\theta^{\prime}})|\theta)ds.$

Instead of keeping only the samples whose distance is smaller than some
apriori fixed threshold $\delta$, the reference table acceptance-rejection ABC
scheme [8, 12] first produces a reference table $\\{\theta_{j},D_{j}\\}$,
$j=1\ldots,n$, and then selects the threshold level $\delta$ as the $q$-th
percentile of the calculated distances $D_{j}$. This procedure has the
computational advantage of fixing the number of drawings $n$ in advance.

Acceptance-rejection ABC and its variant are computationally inefficient by
construction, as the proposals $\theta^{\prime}$ are sampled from the prior
distribution throughout. To tackle this, here we consider SMC-ABC [3, 13, 35],
a sequential algorithm using proposal samplers constructed based on the kept
sampled values (called particles) at the previous iteration
$r=1,\ldots,r_{\textrm{last}}$, yielding samplers defined on more likely
parameter regions. This increases the acceptance probability of each particle,
yielding intermediate approximate posterior distributions during several
iterations which move closer and closer to the desired posterior. In
particular, SMC-ABC works as follows: At iteration one, acceptance-rejection
ABC is run, sampling from the prior distribution until $M$ particles
$\Theta_{1}=(\theta_{1}^{(1)},\ldots,\theta_{1}^{(M)})$ have been accepted,
i.e., have yielded a distance smaller than an initial threshold $\delta_{1}$.
Then, the initial weights are set to
$w_{1}=(w_{1}^{(1)},\ldots,w_{1}^{(M)})=(1/M,\ldots,1/M)$. At iteration $r>1$,
a particle $\theta$ is initially sampled from the set of kept candidates
$\Theta_{r-1}$ of the previous iteration with the corresponding weights
$w_{r-1}$ and then perturbed to a value
$\theta^{*}\sim\mathcal{K}(\cdot|\theta)$, where $\mathcal{K}$ is a suitable
perturbation kernel. Such kernel is also known as proposal sampler or
important sampler, and is commonly assumed to be Gaussian [13, 16], even if
other possibilities have been proposed [27]. Here, we consider the optimized
multivariate normal sampler as proposed in [16]. In particular, for given
$j,l\in\\{1,\ldots,M\\}$ at iteration $r$, the perturbation kernel is given by

$\mathcal{K}_{r}\left(\theta_{r}^{(j)}\Big{|}\theta_{r-1}^{(l)}\right)=(2\pi)^{-d/2}\left(\det\hat{\Sigma}_{r}\right)^{-1/2}\exp\left(-\frac{1}{2}\left[\theta_{r}^{(j)}-\theta_{r-1}^{(l)}\right]^{\top}\hat{\Sigma}_{r}^{-1}\left[\theta_{r}^{(j)}-\theta_{r-1}^{(l)}\right]\right),$
(9)

where $\hat{\Sigma}_{r}$ is twice the weighted empirical covariance matrix
obtained from the previous population $\Theta_{r-1}$ and
$\hat{\Sigma}^{-1}_{r}$ denotes its inverse. Synthetic data
$\tilde{y}_{\theta^{*}}$ are then simulated conditioned on the perturbed
$\theta^{*}$, which is accepted if
$d(s(y),s(\tilde{y}_{\theta^{*}}))<\delta_{r}$, with
$\delta_{r}<\delta_{r-1}$. This is repeated until $M$ particles
$\Theta_{r}=(\theta_{r}^{(1)},\ldots,\theta_{r}^{(M)})$ have been accepted,
after which the corresponding important weights
$(\tilde{w}_{r}^{(1)},\ldots,\tilde{w}_{r}^{(M)})$ are computed as

$\tilde{w}_{r}^{(j)}=\pi\left(\theta_{r}^{(j)}\right)/\sum\limits_{l=1}^{{M}}w_{r-1}^{(l)}\mathcal{K}_{r}\left(\theta_{r}^{(j)}\Big{|}\theta_{r-1}^{(l)}\right)$

and then normalised via

${w}_{r}^{(j)}=\tilde{w}_{r}^{(j)}/\sum\limits_{l=1}^{M}\tilde{w}_{r}^{(l)}.$

This procedure is then repeated over several iterations until a suitable
stopping criterion is reached. The final SMC-ABC posterior is obtained by
sampling from the particles kept at the last iteration $r_{\textrm{last}}$
with probabilities given by the corresponding normalized weights.

###  4.2 Adjusted SMC-ABC algorithm for real and binary parameters

Here, we are interested in estimating both continuous (real-valued) model
parameters and discrete ($\\{0,1\\}$-valued) network parameters. To tackle
this, we adapt the classical SMC-ABC algorithm with Gaussian proposal samplers
to the specific case of additional $\\{0,1\\}$-valued parameters, noting that
the scheme can be adjusted to other distributions with finite support as well.
The method is reported in Algorithm 2 and abbreviated as nSMC-ABC, where the
“n” indicates that it includes network estimation.

We split the parameter vector $\theta$ in continuous and discrete parts as

$\theta:=(\theta_{c},\theta_{d}),$

where $\theta_{c}$ contains continuous (real-valued) model parameters, while
$\theta_{d}$ consists of the discrete ($\\{0,1\\}$-valued) network parameters.
We denote by $c_{n}$ and $d_{n}$ the dimensions of these two vectors,
respectively, and their entries by $\theta_{c}^{k}$, $k=1,\ldots,c_{n}$, and
$\theta_{d}^{k}$, $k=1,\ldots,d_{n}$. We write $\Theta_{c,r}$ and
$\Theta_{d,r}$ to denote the $M$ kept continuous (resp. discrete) particles at
iteration $r$, where each particle $\theta_{r}^{(j)}$ is represented as

$\theta_{r}^{(j)}=\bigl{(}\theta_{c,r}^{(j)},\theta_{d,r}^{(j)}\bigr{)},\quad
j=1,\ldots,M.$

We denote by $\pi^{c}$ and $\pi^{d}$ the prior distributions of $\theta_{c}$
and $\theta_{d}$, and by $\mathcal{K}_{r}^{c}$ and $\mathcal{K}_{r}^{d}$ the
perturbation kernels applied at the $r$-th iteration, respectively.

At iteration one, we run the canonical acceptance-rejection ABC, sampling from
the prior (lines $2-10$ of Algorithm 2). At iteration $r>1$, we draw a
particle $\theta_{c,r}^{(j)}$ for the continuous parameters, which we sample
from the weighted set $\\{\Theta_{c,r-1},{w}_{r-1}\\}$ (line 17 of Algorithm
2). However, differently from the classical SMC-ABC sampler, each entry of a
particle $\theta_{d,r}^{(j)}$ for the binary parameters is drawn from a
Bernoulli distribution with probability given by the sample mean derived from
the previous population (line 19 of Algorithm 2). Then, instead of jointly
perturbing $\theta_{r}^{(j)}=(\theta_{c,r}^{(j)},\theta_{d,r}^{(j)})$ with a
joint Gaussian proposal kernel $\mathcal{K}_{r}$, as done in the classical
SMC-ABC, here we perturb them independently with respect to some continuous
and discrete perturbation kernels $\mathcal{K}_{r}^{c}$ and
$\mathcal{K}_{r}^{d}$ (lines $18$ and $20$ of Algorithm 2). This independence
assumption comes from the fact that there is no natural relation between the
network and continuous parameters, and it turns out to be beneficial, as it
reduces the computational cost required to perturb a particle.

Algorithm 2 Adjusted SMC-ABC for network inference (nSMC-ABC)
Input: Summaries $s(y)$ of the observed data $y$, prior distributions
$\pi^{c}$ and $\pi^{d}$, perturbation kernels $\mathcal{K}_{r}^{c}$ and
$\mathcal{K}_{r}^{d}$, number of kept samples per iteration $M$, initial
threshold $\delta_{1}$
Output: Samples from the nSMC-ABC posterior

1:Set $r=1$

2:for $j=1:M$ do

3: repeat

4: Sample $\theta_{d}$ from $\pi^{d}$ and $\theta_{c}$ from $\pi^{c}$, and set
$\theta=(\theta_{c},\theta_{d})$

5: Conditioned on $\theta$, simulate a synthetic dataset $\tilde{y}_{\theta}$
from the observed output $Y$

6: Compute the summaries $s(\tilde{y}_{\theta})$

7: Calculate the distance $D=d\bigl{(}s(y),s(\tilde{y}_{\theta})\bigr{)}$

8: until $D<\delta_{1}$

9: Set $\theta_{d,1}^{(j)}=\theta_{d}$ and $\theta_{c,1}^{(j)}=\theta_{c}$

10:end for

11:Initialize the weights by setting each entry of
$w_{1}=(w_{1}^{(1)},\ldots,w_{1}^{({M})})$ to $1/{M}$

12:repeat

13: Set $r=r+1$

14: Determine $\delta_{r}<\delta_{r-1}$

15: for $j=1:M$ do

16: repeat

17: Sample $\theta_{c}$ from the weighted set $\\{\Theta_{c,r-1},{w}_{r-1}\\}$

18: Perturb $\theta_{c}$ to obtain $\theta_{c}^{*}$ from
$\mathcal{K}^{c}_{r}(\cdot|\theta_{c})$

19: Sample $\theta_{d}^{k}$, $k=1,\ldots,d_{n}$, from
$\textrm{Bernoulli}(\hat{p}_{r}^{k})$, where
$\hat{p}_{r}^{k}=\frac{1}{M}\sum\limits_{l=1}^{M}\theta^{k,(l)}_{d,r-1}$

20: Perturb $\theta_{d}=(\theta_{d}^{1},\ldots,\theta_{d}^{d_{n}})$ to obtain
$\theta_{d}^{*}$ from $\mathcal{K}^{d}_{r}(\cdot|\theta_{d})$

21: Conditioned on $\theta^{*}=(\theta_{c}^{*},\theta_{d}^{*})$, simulate a
dataset $\tilde{y}_{\theta^{*}}$ from the observed output $Y$

22: Compute the summaries $s(\tilde{y}_{\theta^{*}})$

23: Calculate the distance $D=d\bigl{(}s(y),s(\tilde{y}_{\theta^{*}})\bigr{)}$

24: until $D<\delta_{r}$

25: Set $\theta_{d,r}^{(j)}=\theta_{d}^{*}$ and
$\theta_{c,r}^{(j)}=\theta_{c}^{*}$

26: Set
$\tilde{w}_{r}^{(j)}=\pi^{c}\left(\theta_{c,r}^{(j)}\right)/\sum\limits_{l=1}^{{M}}w_{r-1}^{(l)}\mathcal{K}_{r}^{c}\left(\theta_{c,r}^{(j)}\Big{|}\theta_{c,r-1}^{(l)}\right)$

27: end for

28: Normalise the weights
$w_{r}^{(j)}=\tilde{w}_{r}^{(j)}/\sum\limits_{l=1}^{{M}}\tilde{w}_{r}^{(l)}$,
for $j=1,\ldots,{M}$

29:until stopping criterion is reached

30:Return the final $\Theta_{d,r_{\textrm{last}}}$ and
$\\{\Theta_{c,r_{\textrm{last}}},w_{r_{\textrm{last}}}\\}$.

######  Remark 4.1.

The proposed nSMC-ABC Algorithm 2 is related to ABC for model selection as
discussed in [39], in the sense that each possible network (obtained for a
given combination of the binary parameters) may be interpreted as a model.
Their algorithm samples real-valued parameter candidates conditioned on a
given model and obtains the posterior distribution of the model based on the
number of particles kept under each model. Thus, it is suitable for a low
number of models, while in our case we would obtain a prohibitive number of
possible models, namely $2^{N(N-1)}$ (which equals, e.g., $4096$ for $N=4$).

###  4.3 Ingredients of the nSMC-ABC algorithm

The accuracy and performance of any of the discussed ABC algorithms depend on
various aspects, such as the numerical method used to simulate the data, the
data summaries, the distances, the threshold values, the proposal kernels,
etc. In the following, we describe the choice of these key ingredients.

##### Choice of simulation method

We apply the Strang splitting integrator (8) described in Algorithm 1 to
simulate datasets from the observed output process $Y$ (3) of SDE (5).

##### Choice of summary statistics

For the marginal behaviour of each population $Y^{k}$, $k=1,\ldots,N$, we use
the summary statistics proposed in [8]. In particular, taking advantage of the
underlying geometric ergodicity, we map the $N$-dimensional time series $y$ to
the $N$ estimated marginal invariant densities and spectral densities, the
latter describing the dynamics in the frequency domain. Moreover, we include
summaries capturing the interactions between populations, specifically the
magnitude squared coherence (MSC), a statistic used to examine the similarity
between two signals at various frequencies [38], and the cross-correlation as
a function of lag.

Let $R_{k}$ denote the auto-correlation function of $Y^{k}$ and $R_{jk}$ the
cross-correlation function of $Y_{j}$ and $Y_{k}$, that is

$\displaystyle R_{k}(\tau)$
$\displaystyle=\mathbb{E}[Y^{k}(t)Y^{k}(t+\tau)],\quad k\in\\{1,\ldots,N\\},$
$\displaystyle R_{jk}(\tau)$
$\displaystyle=\mathbb{E}[Y^{j}(t)Y^{k}(t+\tau)],\quad
j,k\in\\{1,\ldots,N\\},\ j\neq k,$

where $\tau$ denotes the time-lag. It holds that $R_{jk}(\tau)=R_{kj}(-\tau)$.
The spectral density $S_{k}$ of $Y^{k}$ is given by the Fourier transform of
$R_{k}(\tau)$, that is

$S_{k}(\nu)=\mathcal{F}\\{R_{k}\\}(\nu)=\int\limits_{-\infty}^{\infty}R_{k}(\tau)e^{-i2\pi\nu\tau}\
d\tau,\quad k\in\\{1,\ldots,N\\},$

where $\nu$ denotes the frequency. Similarly, the cross-spectral density
$S_{jk}$ of $Y^{j}$ and $Y^{k}$ is given by

$S_{jk}(\nu)=\mathcal{F}\\{R_{jk}\\}(\nu)=\int\limits_{-\infty}^{\infty}R_{jk}(\tau)e^{-i2\pi\nu\tau}\
d\tau.$

The MSC is defined as

$Z_{jk}(\nu):=\frac{|S_{jk}(\nu)|^{2}}{S_{j}(\nu)S_{k}(\nu)},\quad
j,k\in\\{1,\ldots,N\\},\ j\neq k,$

where $|\cdot|$ denotes the magnitude.

The summary functions are thus the marginal densities, denoted by $f_{k}$, the
marginal spectral densities $S_{k}$, $k=1,\ldots,N$, the symmetric MSCs
$Z_{jk}$, $j,k=1,\ldots,N$, $k>j$, and the non-symmetric cross-correlations
$R_{jk}$, $j,k=1,\ldots,N$, $j\neq k$. These functions are estimated from the
available observations, see Section 4.4 for further details. Denoting the
estimates by $\hat{f}_{k},\hat{S}_{k}$, $\hat{Z}_{jk}$ and $\hat{R}_{jk}$,
respectively, the proposed summaries $s(\cdot)$ of a dataset $y$ are defined
as

$s(y):=\left\\{\hat{f}_{k},\hat{S}_{k},\hat{Z}_{jk},\hat{R}_{jk}\right\\}.$
(10)

##### Choice of distance measure

Following [8], as distances between the summary functions (10) of the observed
and simulated datasets, we use the integrated absolute error (IAE) given by

$\textrm{IAE}(g_{1},g_{2}):=\int\limits_{\mathbb{R}}|g_{1}(x)-g_{2}(x)|\ dx\
\in\ \mathbb{R}^{+},$

approximated by rectangular integration. The distance between the summaries of
the observed dataset $y$ and a simulated dataset $\tilde{y}_{\theta}$ is then
defined as

$\displaystyle D(s(y),s(\tilde{y}_{\theta}))$ $\displaystyle:=$ $\displaystyle
v_{1}\frac{1}{N}\sum\limits_{k=1}^{N}\textrm{IAE}(\hat{S}_{k},\tilde{S}_{k})+v_{2}\frac{1}{N(N-1)/2}\sum\limits_{j=1,k>j}^{N}\textrm{IAE}(\hat{Z}_{jk},\tilde{Z}_{jk})$
(11) $\displaystyle+v_{3}\frac{1}{N(N-1)}\sum\limits_{j,k=1,j\neq
k}^{N}\textrm{IAE}(\hat{R}_{jk},\tilde{R}_{jk})+v_{4}\frac{1}{N}\sum\limits_{k=1}^{N}\textrm{IAE}(\hat{f}_{k},\tilde{f}_{k}),$

where $\\{\hat{f}_{k},\hat{S}_{k},\hat{Z}_{jk},\hat{R}_{jk}\\}$ and
$\\{\tilde{f}_{k},\tilde{S}_{k},\tilde{Z}_{jk},\tilde{R}_{jk}\\}$ denote the
summaries of the observed and simulated datasets, respectively. The values
$v_{l}\geq 0$, $l=1,2,3,4$, are weights guaranteeing that the different
summary functions have a comparable impact on the distance measure (cf.
Section 4.4).

##### Choice of prior distributions

We use independent uniform priors

$\theta_{c}^{k}\sim\textrm{U}\left(\theta_{c,\textrm{min}}^{k},\theta_{c,\textrm{max}}^{k}\right),\quad
k=1,\ldots,c_{n},$

for the continuous parameters and independent Bernoulli priors

$\theta_{d}^{k}\sim\textrm{Bernoulli}\left(p\right),\quad p=\frac{1}{2},\quad
k=1,\ldots,d_{n},$

for the discrete parameters. In particular,

$\displaystyle\pi^{c}(\theta_{c})$
$\displaystyle=\prod_{k=1}^{c_{n}}\pi(\theta_{c}^{k})=\prod_{k=1}^{c_{n}}\frac{1}{\theta_{c,\textrm{max}}^{k}-\theta_{c,\textrm{min}}^{k}},$
$\displaystyle\pi^{d}(\theta_{d})$
$\displaystyle=\prod_{k=1}^{d_{n}}\pi(\theta_{d}^{k})=\prod_{k=1}^{d_{n}}\left(\frac{1}{2}\right)^{\theta_{d}^{k}}\left(\frac{1}{2}\right)^{1-\theta_{d}^{k}}=\frac{1}{2^{d_{n}}}.$

##### Choice of proposal samplers

For the continuous parameters, we choose $\mathcal{K}_{r}^{c}$ as in (9),
i.e., the optimized Gaussian sampler [16]. A particle $\theta_{c}$ sampled at
iteration $r$ according to the normalized weights of $w_{r-1}$ (line 17 of
Algorithm 2) is thus perturbed to $\theta_{c}^{*}$, a realization of the
multivariate normal distribution $\mathcal{N}(\theta_{c},\hat{\Sigma}_{c,r})$,
where $\hat{\Sigma}_{c,r}$ is twice the weighted covariance matrix obtained
from the previous population $\Theta_{c,r-1}$ (line 18 of Algorithm 2).

For the discrete parameters, a value $\theta_{d}^{k}$, $k=1,\ldots,d_{n}$,
sampled from a Bernoulli distribution at iteration $r$ (line 19 of Algorithm
2) is either kept (with probability $q_{\textrm{stay}}$) or perturbed to
$1-\theta_{d}^{k}$ (line $20$ of Algorithm 2). An explicit expression of the
respective kernel $\mathcal{K}_{r}^{d}$ is given by

$\mathcal{K}_{r}^{d}\left(\theta_{d,r}^{(j)}\Big{|}\theta_{d,r-1}^{(l)}\right)=\prod\limits_{k=1}^{d_{n}}\mathcal{K}_{r}^{d,k}\left(\theta_{d,r}^{k,(j)}\Big{|}\theta_{d,r-1}^{k,(l)}\right)=\prod\limits_{k=1}^{d_{n}}\left(p_{r}^{k,(l)}\right)^{\theta_{d,r}^{k,(j)}}\left(1-p_{r}^{k,(l)}\right)^{1-\theta_{d,r}^{k,(j)}},$
(12)

where

$p_{r}^{k,(l)}=\begin{cases}q_{\textrm{stay}},&\text{if
}\theta_{d,r-1}^{k,(l)}=1\\\ 1-q_{\textrm{stay}},&\text{if
}\theta_{d,r-1}^{k,(l)}=0\end{cases}.$

Throughout this work, the probability $q_{\textrm{stay}}$ is fixed across
iterations (cf. Section 4.4). Alternatively, such probability could also vary,
e.g., it may depend on some statistics of the previous population or the
number of iterations.

##### Choice of threshold levels

The initial threshold $\delta_{1}$ is obtained by a reference table
acceptance-rejection ABC pilot run. Under the given prior, we produce $10^{4}$
distances and then choose $\delta_{1}$ as their median. For $r>1$, the
threshold $\delta_{r}$ is chosen as the median of the $M$ distances computed
at the previous iteration if the acceptance rate of particles at the previous
iteration is larger than 1%. Otherwise $\delta_{r}$ is chosen as the $75$th
percentile of the $M$ distances computed at iteration $r-1$.

##### Choice of stopping criterion

The algorithm is stopped after the acceptance rate of particles at a certain
iteration has dropped below a prefixed threshold (cf. Section 4.4).

###  4.4 Implementation details

The nSMC-ABC method is coded using the statistical software R [30], combined
with the package Rcpp [15]. Algorithm 1 is coded in C++ and then integrated
into the R-Code for Algorithm 2. The for-loops of Algorithm 2 are parallelized
using the packages doParallel and foreach. All experiments are run on a
multiple core High-Performance-Cluster at the University of Klagenfurt.

The summaries (10) are computed as follows: Estimates of the spectral
densities and MSCs are obtained using smoothed periodogram estimators. We
apply the R-function crossSpectrum from the package IRISSeismic to obtain the
spectral densities $\hat{S}_{k}$ and the MSCs $\hat{Z}_{jk}$. Estimates of the
densities $\hat{f}_{k}$ and cross-correlation functions $\hat{R}_{jk}$ are
obtained using the (Gaussian) kernel density estimator density and the
R-function ccf from the stats-package, respectively. The weights $v_{l}$,
$l=1,2,3,4$, in the distance function (11) are obtained as follows: We set
$v_{1}=1$ and obtain $v_{2},v_{3},v_{4}$ by dividing the average area below
the spectral densities of the observed data by the average area below the
MSCs, cross-correlation functions and densities of the observed data,
respectively (the latter equals $1$).

The multivariate normal proposal sampler (9) is calculated using the R-package
mvnfast, which provides computationally efficient tools (e.g., through C++
code) for the multivariate normal distribution. In particular,
$\theta_{c}^{*}\sim\mathcal{N}(\theta_{c},\hat{\Sigma}_{c,r})$ is sampled
using the R-function rmvn, and the multivariate normal density is computed
with the R-function dmvn (lines $18$ and $26$ of Algorithm 2). The weighted
covariance matrix $\hat{\Sigma}_{c,r}$ is estimated with the R-function cov.wt
from the stats-package. The probability $q_{\textrm{stay}}$ of the discrete
perturbation kernel (12) is set to $0.9$ throughout. The number $M$ of kept
particles per iteration and the threshold for the stopping criterion are set
to $500$ and 0.1%, respectively, in all experiments.

Finally, sample code is provided at
https://github.com/IreneTubikanec/networkABC

##  5 Network inference from simulated data

In this section, we test the performance of the proposed nSMC-ABC algorithm on
simulated datasets. All reference datasets are generated up to time $T=20$
with time step $10^{-4}$, and then subsampled with observation time step
$\Delta=2\cdot 10^{-3}$, yielding $10^{4}$ discrete time observations of an
$N$-dimensional process.

###  5.1 Parameter vector and prior distribution

To reduce the number of continuous parameters in the model, we assume that the
coupling strength between two populations decreases with increasing distance,
where the distance is defined by the difference between their subindices. In
particular, for $j,k\in\\{1,\ldots,N\\}$ with $j\neq k$, define

$K_{jk}:=c^{|j-k|-1}L,$

where $L>0$ is a coupling strength parameter and the parameter $0\ll c<1$
describes how fast the network coupling strength decreases with increasing
distance between populations. As the model parameters $A_{k}$,
$k\in\\{1,\ldots,N\\}$, also play a central role in the (non-)activation of
neural populations (see Section 2), we also aim to infer them. The parameters
$\mu_{k}$ and $\sigma_{k}$ are fixed to $90$ and $500$, respectively, and the
remaining continuous parameters are fixed according to the standard values
reported in Table 1. Hence, applying the nSMC-ABC Algorithm 2, the goal is to
infer the $(N+2+N(N-1))$-dimensional parameter vector

$\theta=(\underbrace{A_{1},\ldots,A_{N},L,c}_{\theta_{c}},\underbrace{\textrm{vec}(\mathcal{P})}_{\theta_{d}}),$

where the discrete parameters $\theta_{d}=\textrm{vec}(\mathcal{P})$ are given
by

$\mathcal{P}=\begin{pmatrix}-&\rho_{12}&\ldots&\ldots&\rho_{1N}\\\
\rho_{21}&-&\ddots&&\vdots\\\ \vdots&\ddots&\ddots&\ddots&\vdots\\\
\vdots&&\ddots&-&\rho_{N-1N}\\\ \rho_{N1}&\ldots&\ldots&\rho_{NN-1}&-\\\
\end{pmatrix},$

with $\rho_{jk}\in\\{0,1\\},j,k=1,\ldots,N,j\neq k.$

The prior distributions for $\theta_{c}$ are independent continuous uniforms
with supports,

$A_{k}\sim\textrm{U}(2,4),\ k=1,\ldots,N,\qquad
L\sim\textrm{U}(100,2000),\qquad c\sim\textrm{U}(0.5,1).$

The prior distributions for $\theta_{d}$ are independent Bernoulli
distributions with equal probabilities,

$\rho_{jk}\sim\textrm{Bernoulli}\left(p\right),\quad p=\frac{1}{2},\quad
j,k\in\\{1,\ldots,N\\},\ j\neq k.$ (13)

###  5.2 Estimation results

We focus on $N=4$ neural populations and consider three scenarios, i.e., a
cascade structure (used to simulate the data in Figure 2), a partially and a
fully connected network. The three settings are visualised in Figure 3. We
apply Algorithm 2 to reference datasets ($4$ simultaneously observed
univariate time series, each with $10^{4}$ data points) generated under these
scenarios to estimate $\theta$.

Figure 3: Network structures of $N=4$ neural populations used in the
simulations. (a) Cascade network structure with equal coupling strengths. (b)
Partially connected network. (c) Fully connected network.

##### Setting 1: Cascade structure

Figure 4: Cascade network. (a) nSMC-ABC marginal posterior densities (blue
lines) compared to the prior densities (horizontal red lines) of the
continuous parameters in the cascade scenario, see Figure 3a. The vertical
green lines and the dotted black lines are the true parameter values and the
weighted posterior means, respectively. (b) nSMC-ABC marginal posterior
distributions of the network coupling parameters in the cascade scenario. The
vertical dashed green lines are the true parameter values.

The reference data $y$ of the cascade scenario is simulated under

$\theta=(A_{1},A_{2},A_{3},A_{4},L,\textrm{vec}(\mathcal{P}))=(3.6,3.25,3.25,3.25,700,\textrm{vec}(\mathcal{P})),$
(14)

with $\rho_{12}=\rho_{23}=\rho_{34}=1$ and all other $\rho_{jk}=0$. The
network is visualised in Figure 3a. It shows a similar behaviour as the data
in the right column of Figure 2. In this scenario, the parameter $c$ is not
present and is excluded from $\theta$.

Figure 4a shows the marginal posterior densities (blue lines) and the uniform
prior densities (horizontal red lines) of the continuous parameters. The true
parameter values used to generate the data are indicated by vertical solid
green lines, and the dotted black lines are the weighted posterior means. They
are given by

$(\hat{A}_{1},\hat{A}_{2},\hat{A}_{3},\hat{A}_{4},\hat{L})=(3.615,3.248,3.252,3.254,705.6),$

closely resembling the true values in (14). In Figure 4b, we report the
marginal posterior histograms of the discrete parameters, representing the
coupling directions. The true values (vertical dashed green lines) coincide
with the nSMC-ABC posterior modes. To summarize, Algorithm 2 yields unimodal
and narrow marginal posterior densities of all relevant continuous model
parameters, covering the true values. Moreover, the obtained nSMC-ABC marginal
posteriors for the binary parameters are Bernoulli distributed with
probabilities equal to the true parameters $\rho_{jk}$, successfully
recovering the entire network structure.

##### Setting 2: Partially connected network

The reference data $y$ of the partially connected scenario (visualised in
Figure 3b) is generated using

$\theta=(A_{1},A_{2},A_{3},A_{4},L,c,\textrm{vec}(\mathcal{P}))=(3.6,3.25,3.25,3.25,700,0.8,\textrm{vec}(\mathcal{P})),$

with $\rho_{12}=\rho_{23}=\rho_{34}=\rho_{13}=\rho_{32}=1$ and all other
$\rho_{jk}=0$. The parameter $c$ is included since there is a connection from
Population $1$ to Population $3$.

Figure 5: Partially connected network. (a) nSMC-ABC marginal posterior
densities (blue lines) compared to the prior densities (horizontal red lines)
of the continuous parameters in the partially connected scenario, see Figure
3b. The vertical green lines and the dotted black lines are the true parameter
values and the weighted posterior means, respectively. (b) nSMC-ABC marginal
posterior distributions of the network coupling parameters in the partially
connected scenario. The vertical dashed green lines are the true parameter
values.

Figure 5a shows the marginal posterior densities (blue lines) and the uniform
prior densities (horizontal red lines) of the continuous parameters
$\theta_{c}^{k},k=1,\ldots,6$. The true parameter values are shown by the
vertical solid green lines, and the dotted black lines are the weighted nSMC-
ABC posterior means, given by

$(\hat{A}_{1},\hat{A}_{2},\hat{A}_{3},\hat{A}_{4},\hat{L},\hat{c})=(3.614,3.222,3.259,3.203,693.4,0.726).$

As before, we obtain unimodal marginal posterior densities covering the true
values for $L$ and $A_{k}$, $k=1,\ldots,4$. However, the inference for $c$ is
not satisfactory. This is probably due to the low information about $c$
available in the data, as it only enters through the connection from
Population $1$ to Population $3$.

Figure 5b shows the nSMC-ABC marginal distributions of the discrete parameters
$\theta_{d}^{k}$, $k=1,\ldots,12$. As before, the posterior modes coincide
with the true values of the $\rho_{jk}$ (vertical dashed green lines). Thus,
the network structure is correctly identified.

##### Setting 3: Fully connected network

Figure 6: Fully connected network. (a) nSMC-ABC marginal posterior densities
(blue lines) compared to the prior densities (horizontal red lines) of the
continuous parameters in the fully connected scenario, see Figure 3c. The
vertical green lines and the dotted black lines are the true parameter values
and the weighted posterior means, respectively. (b) nSMC-ABC marginal
posterior distributions of the network coupling parameters in the fully
connected scenario. The vertical dashed green lines are the true parameter
values.

In the fully connected scenario (Figure 3c) the reference data $y$ is obtained
under

$\theta=(A_{1},A_{2},A_{3},A_{4},L,c,\textrm{vec}(\mathcal{P}))=(3.25,3.25,3.25,3.25,700,0.8,\textrm{vec}(\mathcal{P})),$

with $\rho_{jk}=1$ for all $j,k$.

Figure 6a depicts the marginal posterior densities (blue lines) and the
uniform prior densities (horizontal red lines) of the continuous parameters,
the true parameter values (vertical green lines), and the weighted nSMC-ABC
posterior means (dotted black lines), given by

$(\hat{A}_{1},\hat{A}_{2},\hat{A}_{3},\hat{A}_{4},\hat{L},\hat{c})=(3.235,3.247,3.228,3.246,719.7,0.791).$

As in the previous two scenarios, we obtain unimodal and narrow marginal nSMC-
ABC posterior densities, covering the true values, for $L$ and $A_{k}$,
$k=1,\ldots,4$. Compared to the partially connected case, we obtain slightly
better inference for $c$. Similar to the cascade and partially connected
scenarios, the network structure is correctly identified, see Figure 6b.

##### Acceptance rate, runtime and effective sample size

In Figure 7a, we report the acceptance rate of particles and the runtime (in
min) of Algorithm 2 as functions of the iterations $r$, for the cascade (solid
black lines), partially connected (dotted red lines) and fully connected
(dashed blue lines) network scenarios. In all settings, Algorithm 2 terminates
after iteration $43$, where the acceptance rate drops below the threshold of
$0.1\%$. The slope of both the acceptance rate and the runtime curve is
comparable in the three scenarios. Remarkably, the correct network is inferred
(i.e., the marginal posterior modes of the discrete parameters coincide with
the true values) earlier than the algorithm terminates, i.e., from iteration
$9$, $20$ and $13$ for the cascade, partially connected and fully connected
scenarios, respectively. These iteration numbers are indicated by the vertical
lines. Note, however, that while the marginal posterior modes at these
iterations are the same as at the final one, the marginal posterior
distributions itself may still differ.

Finally, we have also considered the “effective sample size” (ESS),
$1\leq\textrm{ESS}\leq M$, calculated at iteration $r$ as
$1/(\sum_{l=1}^{M}(w_{r}^{(l)})^{2})$, as a measure of the effectiveness of
the SMC sampler, that is, intuitively, how many of the $M$ particles are
“relevant”. The ESS of the three scenarios look similar, mainly oscillating
between 350 and 450 (recall that $M$ is set to $500$ here), with overall
median across the iterations equal to 411, 418 and 407, respectively (figure
not reported).

(a) Simulated data (b) Measured EEG data

Figure 7: Acceptance rate and runtime as functions of the iterations $r$ of
Algorithm 2. (a) Simulated data, time measured in minutes. The vertical lines
mark those iterations from which the correct network is inferred, i.e., the
coupling parameters $\rho_{jk}$ are correctly estimated via the posterior
modes. The solid black, dotted red and dashed blue lines correspond to the
cascade, partially connected and fully connected scenario, respectively. (b)
EEG data with $N=4$ channels, time measured in hours. The dotted red and solid
black lines correspond to the before and during seizure periods, respectively.
The vertical lines mark that iteration from which the network estimate does
not change anymore.

##  6 Network inference from EEG data with epileptic activity

After validating the proposed algorithm on simulated data, we now use it to
estimate the connectivity structure from real EEG data.

###  6.1 Description of the data

In [33], a study on $22$ patients from the Children’s Hospital Boston, who
experience epileptic seizures, was presented. The aim of the study was to
detect seizure periods in multiple hour EEG recordings of each individual. The
data are available at https://www.physionet.org/content/chbmit/1.0.0/.

Here, we analyse recordings of an $11$ year old female patient prior to and
during seizure, from the edf-file chb01_03 in the above link. Following [26,
32], we consider the four channels FP1-F7, FP1-F3, FP2-F4 and FP2-F8, where FP
refers to the frontal lobes (the first two on the left hemisphere and the
second two on the right) and F to a set of electrodes placed behind them. The
electrode locations are according to the international 10–20 system for EEG
measurements. The recordings are visualised in Figure 8, where the vertical
dotted red lines separate the data into the period before and during seizure,
lasting 40 seconds [$s$] each, with the period from $2996$ to $3036s$
classified as a seizure in [33]. The data are sampled with a frequency of
$256$ Hz, corresponding to $10240$ discrete time measurements during each of
the two $40s$ periods. To put the recordings on the same scale as the model,
we rescale the data by multiplying each data point with the factor $0.05$.

Figure 8: EEG recordings of a $11$ year old female patient. The interval
$[2996,3036]s$ (measurements to the right of the vertical dotted red lines)
has been defined as a seizure in [33].

###  6.2 Inference

We now fit the stochastic multi-population JR-NMM to the data prior to and
during seizure, aiming to infer the underlying network structure as well as
relevant continuous model parameters in these two regimes via our proposed
nSMC-ABC Algorithm 2.

##### Parameter vector and prior distribution

In the following, we identify the channels FP1-F7, FP1-F3, FP2-F4 and FP2-F8
by Population $1$, $2$, $3$, and $4$, respectively. As the distance between
channels on different hemispheres is larger than that on the same, we have
that the distance between the channels FP1-F3 and FP2-F4 is larger than that
between FP1-F7 and FP1-F3 or between FP2-F4 and FP2-F8. For this reason, we
assume the following matrix of coupling strength parameters

$K=\begin{pmatrix}-&K_{12}&K_{13}&K_{14}\\\ K_{21}&-&K_{23}&K_{24}\\\
K_{31}&K_{32}&-&K_{34}\\\
K_{41}&K_{42}&K_{43}&-\end{pmatrix}=\begin{pmatrix}-&L&c^{2}L&c^{3}L\\\
L&-&cL&c^{2}L\\\ c^{2}L&cL&-&L\\\ c^{3}L&c^{2}L&L&-\end{pmatrix},$

with unknown $L$ and $c$. Moreover, we assume the activation parameters
$A_{k}$, $k=1,\ldots,4$, the noise intensity parameters
$\sigma_{l}:=\sigma_{1}=\sigma_{2}$ (left hemisphere) and
$\sigma_{r}:=\sigma_{3}=\sigma_{4}$ (right hemisphere) as well as the input
parameters $\mu_{l}:=\mu_{1}=\mu_{2}$ (left hemisphere) and
$\mu_{r}:=\mu_{3}=\mu_{4}$ (right hemisphere) to be unknown. For the remaining
parameters, we use the standard values reported in Table 1, except for setting
$b=20$ and $C=70$ (values chosen based on pilot experiments). We thus aim to
infer the $(10+12)$-dimensional parameter vector

$\theta=(A_{1},A_{2},A_{3},A_{4},L,c,\sigma_{l},\sigma_{r},\mu_{l},\mu_{r},\textrm{vec}(\mathcal{P})).$

We choose independent continuous uniform prior distributions for $\theta_{c}$
with broad supports,

$\displaystyle A_{k}\sim\textrm{U}(1,30),\ k=1,\ldots,4,\quad
L\sim\textrm{U}(100,5000),\quad c\sim\textrm{U}(0.5,1),$
$\displaystyle\sigma_{l},\sigma_{r}\sim\textrm{U}(100,15000),\quad\mu_{l},\mu_{r}\sim\textrm{U}(1,300).$

The priors for $\theta_{d}$ are independent Bernoulli distributions with equal
probabilities, see (13).

##### Estimation results

In Figure 9a, we report the marginal posterior densities (green lines: before
seizure, blue lines: during seizure) and the uniform prior densities
(horizontal red lines) of the continuous parameters
$\theta_{c}^{k},k=1,\ldots,10$. All posteriors for the parameters during
seizure are unimodal, showing a clear update of the priors. Comparing the
individual population parameters $A_{k}$, $k=1,...,4$, for the two periods,
there is a larger activation in all $4$ neural populations during seizure.
Moreover, as expected, the posteriors of the noise intensity parameters
$\sigma_{l}$ and $\sigma_{r}$ assume larger values in both brain hemispheres
during seizure. While we obtain a clear estimate of the posterior for the
coupling strength parameter $L$ during seizure, this is not the case before
seizure. However, the posterior suggests a coupling strength far away from
zero also in this scenario.

Figure 9: Estimation results for EEG recordings. Three models are considered:
Populations 1-4, modelling all four channels; Populations 1-2, modelling the
left hemisphere; Populations 3-4, modelling the right hemisphere. (a) nSMC-ABC
marginal posterior densities of the continuous parameters. Before seizure:
solid green ($N=4$), dotted orange ($N=2$, left hemisphere), dotted brown
($N=2$, right hemisphere). During seizure: solid blue ($N=4$), dashed grey
($N=2$, left hemisphere), dashed black ($N=2$, right hemisphere). The solid
red lines are the uniform prior distributions. (b) nSMC-ABC marginal posterior
distributions of the discrete network parameters before seizure (top panels)
and during seizure (bottom panels). The vertical dashed orange and red lines
denote the estimated posterior modes for the reduced models with $N=2$
populations of the left and right hemispheres, respectively.

Figure 9b shows the marginal posterior histograms for the coupling direction
parameters $\theta_{d}^{k},k=1,\ldots,12$ (top panels: before seizure, bottom
panels: during seizure). The posteriors during the seizure period provide a
clear connectivity structure, with the inferred network visualised in Figure
10b. The inferred network prior to seizure (Figure 10a) is similar, except for
a missing connection from Population 2 (FP1-F3) to Population 1 (FP1-F7) on
the left hemisphere, present during seizure. The inferred network also differs
for $\rho_{42}$ and $\rho_{43}$, with a connection from Population 4 (FP2-F8)
to 2 (FP1-F3) and to 3 (FP2-F4) present only before seizure. However, for
these parameters the results are less evident (estimated probabilities for a
connection smaller than 0.75 in the before seizure scenario).

To further investigate the differences, if any, between the two hemispheres,
we also apply our nSMC-ABC Algorithm 2 only to the two channels from the same
hemisphere, i.e., left hemisphere (FP1-F7 and FP1-F3) and right hemisphere
(FP2-F8 and FP2-F4), both before and during seizure. The dotted orange/dashed
grey (left hemisphere) and dotted brown/dashed black (right hemisphere) lines
in Figure 9a are the marginal nSMC-ABC posteriors of the corresponding
parameters before/during seizure. The dashed orange (left hemisphere) and red
(right hemisphere) vertical lines in Figure 9b indicate the inferred
connections, estimated as posterior modes. The corresponding posterior
probabilities that $\rho_{12}$, $\rho_{21}$, $\rho_{34}$, $\rho_{43}$ are
equal to one are $1$, $0.02$, $0.11$, $0.984$ and $1$, $1$, $0.998$, $1$
before and during seizure, respectively.

When only using data from the left hemisphere, the inferred connections
(namely, a unidi-rectional connection from Population 1 to 2 before seizure,
and bidirectional connection during seizure) agree with those obtained in the
$4$-population-case. During the seizure period, the coupling strength $L$ for
the left hemisphere is slightly larger than for $N=4$ (grey versus blue lines
in the top right panel of Figure 9a). Moreover, the coupling strengths of the
left hemisphere seem to be larger than the right one, both before and during
seizure. For the period prior to seizure, the inferred connections
($\hat{\rho}_{34}=0,\hat{\rho}_{43}=1$) from the two populations on the right
hemisphere (i.e., Population $3$, FP2-F4, and Population $4$, FP2-F8) also
agree with the 4-population-case. However, this is not the case during
seizure, as the bidirectional connections detected now for $N=2$ were not
present when using data from all 4 channels, even if the coupling strength $L$
and $c$ assume smaller values than before (black vs blue lines in the top
right panel and left bottom panel of Figure 9a).

(a) Before seizure (b) During seizure

Figure 10: Inferred network from EEG recordings. Estimated network ($N=4$,
channels FP1-F7, FP1-F3, FP2-F4 and FP2-F8) for a 11 year old female patient
before (a) and during (b) epileptic seizure. The dotted orange connections in
(a) are estimated with a probability less than 0.75 (cf. top panels of Figure
9b).

##### Fitted summaries

Figure 11: Summaries from EEG recordings and the estimated traces. Odd panels:
before seizure; Even panels: during seizure. Solid black lines: Summaries
derived from the EEG datasets visualised in Figure 8. Grey areas: Range of the
summaries obtained from synthetic datasets simulated using the kept posterior
samples from the full model with four channels ($N=4)$. Dashed lines: Median
values of the summaries derived from simulated datasets with $N=2$ populations
for the left hemisphere (orange) and the right one (red).

Figure 11 provides a comparison of the summary statistics (10) derived from
the EEG recordings (solid black lines) and those obtained from datasets
simulated under the $500$ kept posterior samples of the 4-population model
(grey regions). Median values of the summaries for the 2-population models of
the left (dashed orange lines) and right (dashed red lines) hemispheres are
also reported. The left and right panels of each subfigure correspond to the
before and during seizure scenario, respectively. The match of the observed
and simulated summaries is generally good, except for the MSCs, which are
based on the cross-spectral densities. Moreover, the oscillations of the
cross-correlation functions before seizure are not captured well, suggesting
some lack of fit of the multi-population JR-NMM to the EEG data.

##### Acceptance rate and runtime

In Figure 7b, we report the acceptance rate of particles and the runtime (in
hours) of the $4$-population-experiments as functions of the iteration $r$ of
the nSMC-ABC Algorithm 2 before (dotted red lines) and during (solid black
lines) seizure. The acceptance rate dropped below the final threshold of
$0.1\%$, causing the algorithm to terminate, after 59 and 72 iterations,
respectively. The vertical lines indicate iteration $44$ and $51$, from which
the network estimates (visualised in Figure 10) do not change anymore.

##  7 Conclusion and discussion

Understanding how different regions in the brain are connected and how they
function together, e.g., during epileptic seizure, is a problem of large
interest in the field of neuroscience. We contribute to this task by
investigating a stochastic neural mass model, which describes the activity of
multiple coupled neural populations in the brain. In particular, we extend the
$1$-population stochastic JR-NMM [1, 18] to a model of $N$ coupled neural
populations. This $6N$-dimensional SDE contains both a variety of continuous
biologically motivated model parameters and $\\{0,1\\}$-valued coupling
direction parameters, which characterize the network structure among the $N$
interacting populations. The goal is to estimate both types of parameters from
$N$ simultaneously observed time series datasets, corresponding to a multi-
channel EEG recording. This task is particularly challenging, as none of the
$6N$ coordinates of the model is directly observed. Indeed, only the
$N$-dimensional output (3), a multivariate linear function of the coordinates
of the equation, is measurable with EEG recording techniques. Moreover, due to
the intricacy of the model, the underlying likelihood function is not
available.

We tackle these issues by extending the likelihood-free ABC method presented
in [8] to the case of multiple interacting components. In particular, we
consider the gold-standard SMC-ABC algorithm and adapt it to sample and
perturb the continuous and binary parameter independently. The proposed nSMC-
ABC algorithm requires two key ingredients: a reliable and computationally
efficient numerical simulation method for data generation, and suitable
summary statistics to compare the observed and simulated data, capturing the
model dynamics and features independently of the underlying intrinsic
stochasticity. First, we construct a numerical method to simulate from the
model. The derived scheme is based on the splitting approach, it exploits the
Hamiltonian structure of the model, and it is automatized for an arbitrary
number of populations $N$. A standard integrator such as the Euler-Maruyama
method would either lead to drastic fails in the estimation or make the
algorithm computationally infeasible, as it would require a significantly
smaller time discretisation step to capture the dynamics of the model (see [8]
for an illustration of this issue on the stochastic JR-NMM with $N=1$).
Second, we propose data summaries which take into account both individual
signals as well as pairwise connectivity measures among signals. In
particular, we map the $N$-dimensional time series to the $N$ corresponding
densities and spectral densities, and derive the $N(N-1)/2$ MSCs, which are
based on the cross-spectral densities, and the $N(N-1)$ cross-correlation
functions. These quantities are derived from a given dataset using standard
estimators, e.g., a kernel density estimator and a smoothed periodogram
estimator.

The resulting algorithm yields excellent estimation results when applied to
simulated data in various scenarios, ranging from a cascade connectivity
structure to a fully connected network of four populations. Applying the
algorithm to real EEG data with epileptic activity yields a clear estimate of
both the relevant continuous model parameters and the underlying network
structure. Comparing the estimation results with those obtained from the data
before seizure, we observe a larger activation in each of the individual
neural populations, a larger noise intensity in both the left and the right
brain hemisphere, a stronger connectivity in the left brain hemisphere, and a
reduced one in the right brain hemisphere during seizure. The match of summary
statistics derived from the experimental data and those generated using the
weighted kept nSMC-ABC particles is generally good, except for the MSCs and
the oscillations in the cross-correlations before seizure. This may be due to
a lack of fit of the model, possibly due to the data not being Markovian, the
chosen modelling of the coupling strength parameters or some of the underlying
fixed parameters being slightly off.

The successful inference of the network parameters together with a relatively
large number of continuous model parameters in the case of $N=4$ populations
is remarkable, considering that previously, in the $1$-population JR-NMM, only
up to four continuous model parameters were simultaneously identified via ABC
[8] or neural posterior estimation approaches [31]. Finally, we emphasise that
the proposed nSMC-ABC algorithm, along with the presented summaries, can be
directly applied to other coupled ergodic SDE models, as long as a suitable
numerical simulation method for data generation is available. Moreover, the
proposed methodology for network inference may be incorporated into more
evolved ABC algorithms such as the guided sequential schemes introduced in
[27], aiming to further reduce the overall inferential running times.

## References

  * [1] M. Ableidinger, E. Buckwar, and H. Hinterleitner.  A stochastic version of the Jansen and Rit neural mass model: analysis and numerics.  J. Math. Neurosci., 7(8), 2017.
  * [2] L. Arnold.  Stochastic Differential Equations: Theory and Applications.  Wiley, New York, 1974.
  * [3] M. A. Beaumont, J.-M. Cornuet, J.-M. Marin, and C. P. Robert.  Adaptive approximate Bayesian computation.  Biometrika, 96(4):983–990, 2009.
  * [4] M. A. Beaumont, W. Zhang, and D. J. Balding.  Approximate Bayesian computation in population genetics.  Genetics, 162(4):2025–2035, 2002.
  * [5] S. Blanes, F. Casas, and A. Murua.  Splitting and composition methods in the numerical integration of differential equations.  Bol. Soc. Esp. Mat. Apl., 45, 2009.
  * [6] S. Blanes, F. Diele, C. Marangi, and S. Ragni.  Splitting and composition methods for explicit time dependence in separable dynamical systems.  J. Comput. Appl. Math., 235(3):646–659, 2010.
  * [7] E. Buckwar, A. Samson, M. Tamborrino, and I. Tubikanec.  A splitting method for SDEs with locally Lipschitz drift: Illustration on the FitzHugh-Nagumo model.  Appl. Numer. Math., 179:191–220, 2022.
  * [8] E. Buckwar, M. Tamborrino, and I. Tubikanec.  Spectral density-based and measure-preserving ABC for partially observed diffusion processes. An illustration on Hamiltonian SDEs.  Stat. Comput., 30(3):627–648, 2020.
  * [9] Z. Chen, B. Raman, and A. Stern.  Structure-preserving numerical integrators for Hodgkin–Huxley-type systems.  SIAM J. Sci. Comput., 42(1):B273–B298, 2020.
  * [10] J. Chevallier, A. Melnykova, and I. Tubikanec.  Diffusion approximation of multi-class Hawkes processes: Theoretical and numerical analysis.  Adv. Appl. Probab., 53(3):716–756, 2021.
  * [11] D. Cohen.  On the numerical discretisation of stochastic oscillators.  Math. Comput. Simul., 82(8):1478–1495, 2012.
  * [12] J. M. Cornuet, F. Santos, M. Beaumont, C. P. Robert, J.-M. Marin, D. J. Balding, T. Guillemaud, and A. Estoup.  Inferring population history with DIY ABC: A user-friendly approach to approximate Bayesian computation.  Bioinformatics, 24(23):2713–2719, 2008.
  * [13] P. Del Moral, A. Doucet, and A. Jasra.  An adaptive sequential Monte Carlo method for approximate Bayesian computation.  Stat. Comput., 22:1009–1020, 2012.
  * [14] S. Ditlevsen and A. Samson.  Hypoelliptic diffusions: filtering and inference from complete and partial observations.  J. Royal Stat. Soc. B, 81(2):361–384, 2019.
  * [15] D. Eddelbuettel and R. François.  Rcpp: Seamless R and C++ Integration.  J. Stat. Soft., 40(8):1–18, 2011.
  * [16] S. Filippi, P. B. Barnes, J. Cornebise, and M. P. H. Stumpf.  On optimality of kernels for approximate Bayesian computation using sequential Monte Carlo.  Stat. Appl. Genet. Mol. Biol., 12(1):87–107, 2013.
  * [17] J. Foster, G. Reis, and C. Strange.  High order splitting methods for SDEs satisfying a commutativity condition.  arXiv, 2211.11884, 2022.
  * [18] B. H. Jansen and V. G. Rit.  Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns.  Biol. Cybern., 73(4):357–366, 1995.
  * [19] B. H. Jansen, G. Zouridakis, and M. E. Brandt.  A neurophysiologically-based mathematical model of flash visual evoked potentials.  Biol. Cybern., 68:275–283, 1993.
  * [20] C. Kelly and G. Lord.  Adaptive time-stepping strategies for nonlinear stochastic systems.  Ima. J. Numer. Anal., 38(3):1523–1549, 2017.
  * [21] J. R. León and A. Samson.  Hypoelliptic stochastic FitzHugh-Nagumo neuronal model: mixing, up-crossing and estimation of the spike rate.  Ann. Appl. Probab., 28(4):2243–2274, 2018.
  * [22] X. Mao.  Stochastic Differential Equations and Applications.  Horwood Publications, Chichester, 1997.
  * [23] J.-M. Marin, P. Pudlo, C. P. Robert, and R. Ryder.  Approximate Bayesian computational methods.  Stat. Comput., 22(6):1167–1180, 2012.
  * [24] R. McLachlan and G. Quispel.  Splitting methods.  Acta Numer., 11:341–434, 2002.
  * [25] A. Melnykova.  Parametric inference for hypoelliptic ergodic diffusions with full observations.  Stat. Inference Stoch. Process., 23:595–635, 2020.
  * [26] J. Østergaard, A. Rahbek, and S. Ditlevsen.  Oscillating systems with cointegrated phase processes.  J. Math. Biol., 75:845–883, 2017.
  * [27] U. Picchini and M. Tamborrino.  Guided sequential abc schemes for intractable bayesian models.  arXiv, 2206.12235, 2022.
  * [28] P. Pilipovic, A. Samson, and S. Ditlevsen.  Parameter estimation in nonlinear multivariate stochastic differential equations based on splitting schemes.  arXiv, 2211.11884, 2022.
  * [29] Y. Pokern, A. M. Stuart, and P. Wiberg.  Parameter estimation for partially observed hypoelliptic diffusions.  J. Royal Stat. Soc. B, 71(1):49–73, 2009.
  * [30] R Development Core Team.  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria, 2011.
  * [31] P. Rodrigues, T. Moreau, G. Louppe, and A. Gramfort.  Hnpe: Leveraging global parameters for neural posterior estimation.  NeurIPS, 34, 2021.
  * [32] M. G. Ruse, A. Samson, and S. Ditlevsen.  Inference for biomedical data by using diffusion models with covariates and mixed effects.  J. R. Stat. Soc., C: Appl., 69(1):167–193, 10 2019.
  * [33] A. H. Shoeb.  Application of machine learning to epileptic seizure onset detection and treatment.  Dissertation, MIT Health Sciences and Technology Division, 2009\.
  * [34] S. A. Sisson, Y. Fan, and M. Beaumont.  Handbook of Approximate Bayesian Computation.  Chapman & Hall/CRC Handbooks of Modern Statistical Methods. CRC Press, Taylor & Francis Group, 2018.
  * [35] S. A. Sisson, Y. Fan, and M. M. Tanaka.  Sequential Monte Carlo without likelihoods.  Proc. Natl. Acad. Sci. U S A, 104(6):1760–1765, 2007.
  * [36] G. Strang.  On the construction and comparison of difference schemes.  SIAM J. Numer. Anal., 5(3):506–517, 1968.
  * [37] A. H. Strømmen Melbø and D. J. Higham.  Numerical simulation of a linear stochastic oscillator with additive noise.  Appl. Numer. Math., 51(1):89–99, 2004.
  * [38] S. Tiran and P. Maurine.  SCA with Magnitude Squared Coherence.  In S. Mangard, editor, Smart Card Research and Advanced Applications, pages 234–247, Berlin, Heidelberg, 2013. Springer Berlin Heidelberg.
  * [39] T. Toni, D. Welch, N. Strelkowa, A. Ipsen, and M. P. H. Stumpf.  Approximate Bayesian computation scheme for parameter inference and model selection in dynamical systems.  J. Roy. Soc. Interface, 6(31):187–202, 2009.
  * [40] H. F. Trotter.  On the product of semi-groups of operators.  Proc. Am. Math. Soc., 10(4):545–551, 1959.
  * [41] I. Tubikanec, M. Tamborrino, P. Lansky, and E. Buckwar.  Qualitative properties of different numerical methods for the inhomogeneous geometric Brownian motion.  J. Comput. Appl. Math., 406:113951, 2022.
  * [42] F. Wendling, J. J. Bellanger, F. Bartolomei, and P. Chauvel.  Relevance of nonlinear lumped-parameter models in the analysis of depth-EEG epileptic signals.  Biol. Cybern., 83(4):367–378, 2000.