On a framework of data assimilation for neuronal networks

Wenyong Zhang1, Boyu Chen1, Jianfeng Feng2,3, Wenlian Lu1,4,5,6*

1 School of Mathematical Sciences, Fudan University, No. 220 Handan Road,
Shanghai, China

2 Institute of Science and Technology for Brain-Inspired Intelligence, Fudan
University, Shanghai, China

3 Key Laboratory of Computational Neuroscience and Brain-Inspired
Intelligence, Fudan University, Ministry of Education, Shanghai, China

4 Shanghai Center for Mathematical Sciences, Fudan University, Shanghai, China

5 Shanghai Key Laboratory for Contemporary Applied Mathematics, Fudan
University, Shanghai, China

6 Key Laboratory of Mathematics for Nonlinear Science, Fudan University,
Ministry of Education, Shanghai, China

*<EMAIL_ADDRESS>
## Abstract

When handling real-world data modeled by a complex network dynamical system,
the number of the parameters is always even much more than the size of the
data. Therefore, in many cases, it is impossible to estimate these parameters
and however, the exact value of each parameter is frequently less interesting
than the distribution of the parameters may contain important information
towards understanding the system and data. In this paper, we propose this
question arising by employing a data assimilation approach to estimate the
distribution of the parameters in the leakage-integrate-fire (LIF) neuronal
network model from the experimental data, for example, the blood-oxygen-level-
dependent (BOLD) signal. Herein, we assume that the parameters of the neurons
and synapses are inhomogeneous but independently identical distributed
following certain distribution with unknown hyperparameters. Thus, we estimate
these hyperparameters of the distributions of the parameters, instead of
estimating the parameters themselves. We formulate this problem under the
framework of data assimilation and hierarchical Bayesian method, and present
an efficient method named Hierarchical Data Assimilation (HDA) to conduct the
statistical inference on the neuronal network model with the BOLD signal data
simulated by the hemodynamic model. We consider the LIF neuronal networks with
four synapses and show that the proposed algorithm can estimate the BOLD
signals and the hyperparameters with good preciseness. In addition, we discuss
the influence on the performance of the algorithm configuration and the LIF
network model setup.

## Author summary

With the great development of brain neuroimaging technology such as MRI, in
recent decades, we can understand the brain function more generally than
before and the advanced computer technology makes it possible to construct
computational neuronal network models according to the brain neuroimaging
data. However, the number of the parameters of the neuronal network may be
very huge, usually larger than the sample size of observation data. Therefore,
it is difficult to estimate all parameters accurately with avoiding the over-
fitting problem at the same time. This paper proposes a novel framework, named
Hierarchical Data Assimilation(HDA), based on an assumption that the
observation data of the network model is independent of the specific
parameters but depends on the distribution of these parameters, proves the
convergence (asymptotic) and well-posedness of the HDA estimation approach in
the view of hierarchical Bayesian inference, and illustrate the efficiency and
power of the algorithms to estimate the BOLD signal and the hyperparameters by
numerical simulating experiments. We think that this proposed algorithm is a
promising method for a statistical inference method to integrate the large-
scale full-brain computational model with the experimental brain data.

## Introduction

Over recent decades, research within neuroscience focuses on simulating the
human brain to discover the neural mechanisms underlying perceptual, cognitive
and motor functions, and the dynamical properties of the human brain that
appears to guide the development of Artificial Intelligence, Psychiatry,
Information Science, Cognition Sciences, etc. Besides, experiments in animals
brain earlier used perturbative approaches such as implant electrodes or
surgical operations, but these approaches are difficult and unethical in
humans. Thus it is necessary to construct a computational neuronal network
model with computers or silicon chips and design experiments on the virtual
brain [7]. For modeling of neural dynamics, there are also quantities of the
model such as MP neuron [34], Hodgkin–Huxley (HH) model [22], neural mass
models [47] and other models referred to in [32]. Intuitively, we should
choose the most appropriate model but it is impossible to use implanted
electrodes to record the rates and timing of action potentials in the human
brain. Fortunately, with the development of human neuroimaging research,
alterations in brain activity can be detected non-invasively and at the whole-
brain level using electroencephalogram (EEG), magnetoencephalogram(MEG) or
functional MRI (fMRI) [12, 36, 40, 33, 39]. These techniques make simulating
the brain computationally come true: a whole-brain neuronal network has been
modeled by a normal Hopf bifurcation with resting-state MEG brain activity and
BOLD(blood oxygen level dependent) signals [8, 28]. All of these experiments
used the structural connectivity (SC) to fit the functional connectivity (FC)
from neuroimaging data and it is shown that the model parameters can be
adjusted to describe different brain states [24]. Therefore, the model
estimation based on specified datasets is an essential issue in neuronal
network simulation.

Data assimilation(DA) is an efficient method aiming at simulating a dynamic
system by a data-tuned mathematical model which was widely used in extensive
fields. The earliest successful application of DA was to handle data in
atmospheric sciences and weather forecasting which was comprehensively
described in the book [4]. Then, the methodology of data assimilation
gradually extends its application areas especially in geochemistry [11],
oceanography [23], financial econometrics [31], etc. With the development of
DA, several specific algorithms have been established. The classic algorithm
named Kalman filter (KF) provided by Kalman [26] is to deal with the Gaussian
linear dynamic system. Some generalizations of the Kalman filter methodology
such as extended Kalman filter (EKF) [2], unscented Kalman filter (UKF) [25],
ensemble Kalman filter (EnKF) [14] were lately established and proved to be
effective to nonlinear systems. And from the viewpoint of monte carlo (MC)
sampling, there are also many significant algorithms including bootstrap
filter (BF) [19], unscented particle filter (UPF) [10] and sampling importance
resampling (SIS) [18] etc. Recently, Kody Law develops a unified mathematical
perspective of DA in which a Bayesian formulation of the problem provides the
bedrock for the derivation and development of these algorithms [29].

Although DA has played a fundamental role in these areas, it has not been
systemically applied in neuroscience. Path integral approach [38, 42, 27] and
self-organizing state space model [44] have been used in parameter estimation
based on HH model [1]. Article [37] compared estimation performance of EnKF,
BF, and the optimal sequential importance re-sampling (OPT-SIRS) methods in HH
model. However, since the large parameter space of HH model makes robust
analysis difficult, Leaky Integral-Fired model (LIF) is another common neuron
model [41, 6]. Some parameter estimation algorithms for LIF model have been
put forward [35], however, LIF model is seldom simulated via the framework of
DA.

On the other hand, it is worth considering the inhomogeneity of neurons, that
is, the values of the same parameter for different neurons should be
different. For example, the decay/rising time scales are different from
neurons; the synaptic conductances are different from synapses. It was proved
by experiments that the firing activities of neurons have variability with
respect to neurons [43]. So, it is natural to believe that the parameters of
neurons should be set inhomogeneous in the neuronal network model to better
mimic the dynamics. In recent years, with the increase of the neuronal network
scale and the computing power, the data or parameters of the model are more
and more huge. Thus hierarchical Bayesian model as a novel estimation
framework has been proposed and applied in many fields including computer
vision [15], seasonal hurricane simulation [13], and solving boundary value
problems in atmospheric stream function fields [46]. In neuroscience, the
Bayesian inference is commonly used to construct a hierarchical graphical
structure in modeling cortical processing [16, 30], and the size of parameters
of each layer is usually small or large but with special hierarchical
structures such as tree-structured models [20].

In this article, we establish a framework named hierarchical data assimilation
(HDA) for estimating the hyperparameter of the dynamic equation model based on
hierarchical Bayesian inference, give two algorithms as derivatives of DA
methods to estimate the hyperparameter, and show the convergence of our
algorithms from the view of maximum likelihood estimation (MLE). Because the
number of hyperparameters and the observation have the same order of magnitude
much less than model parameters, it avoids over-fitting effectively. In
addition, we describe the neuron network model [5] and hemodynamic model [17,
9] in the following and put our approaches to a test. Finally, we show it is
efficient to estimate the hyperparameter of the model and fit the real fMRI
signals.

## Methods

In this section, we give a detailed description of our HDA frameworks,
neuronal networks, and other methods to show how we can apply HDA to neuronal
networks.

Here we give a brief introduction of our assimilation task and the
experimental results are displayed in section Results. Firstly, the
connections of neurons in networks are generated according to the specific
network topology. Using Euler iterative method, we transform the neurons ODE
model and hemodynamics model into the discrete-time equations. Then
hierarchical neuron networks assimilation estimates the hyperparameters by
iterating two processes: simulation, computing the hidden states including
model parameters and states by numerical calculation, and assimilation,
adjusting the model hidden states by HDA methods such as ensemble Kalman
filter (EnKF) according to Bold observation. It should be noted that EnKF
methods used ensemble samples to calculate the variance of hidden states and
in this paper, if not otherwise specified, ensemble samples are 100 neuron
networks with the same topology structure.

### Hierarchical Data Assimilation

Hierarchical Data Assimilation (HDA) is the crucial approach proposed in this
paper based on Data Assimilation(overview in S1 Appendix.). When the number of
the parameters in the network model is much more than the size of the
observation data, it is impossible to estimate these parameters and the exact
values of these parameters are frequently of no interest but the distribution
of the parameters may contain important information about the system and data.
Here, we establish HDA framework by introducing hyper-parameter $h$ to
describe the distribution of parameters, which may contain the important
information about these parameters.

Let us consider the following general coupled dynamic systems:

$\displaystyle\frac{dx_{i}}{dt}=f_{i}(x_{i},x^{i},\theta_{i}),~{}i=1,\cdots,m.$
(1)

Here, $x_{i}\in{\mathbb{R}}^{n}$ stands for the state vector of agent $i$,
such as the collection of states in LIF model and Balloon-Windkessel model, as
mentioned in Methods, $x^{i}=[x_{j}:~{}j\neq i]$ stands for the vectors from
the other nodes, $\theta_{i}\in{\mathbb{R}}^{p}$ stands for the parameters of
the evolution law of node $i$, which is related to states of other nodes. The
observation of $x=[x_{1},\cdots,x_{m}]$ is:

$\displaystyle y_{t_{l}}=H\circ x_{t_{l}}$ (2)

where $H$ is an observation operator that highly reduces the dimensions, i.e.,
$q\ll m\times n$, and $\\{t_{l}\\}_{l=1}^{L}$ are the sampling time points.
When the parameters $\theta_{i}$ are assumed to be independently identical
distributed (I.I.D) with the density function $\psi(\theta|h)$:

$\displaystyle\theta_{i}\sim\psi(\cdot|h)$

where $h\in{\mathbb{R}}^{r}$ is the hyperparameter vector. Giving a long
recording of $y$ from a large population of nodes, HDA aims at estimating the
hyperparameter $h$. Define

$\displaystyle\theta=[\theta_{1},\cdots,\theta_{n}]^{\top}$

to be regarded as the i.i.d. samples from $\phi(\cdot|h)$.

When dealing with real dataset, the observation is discrete-time with time
step $\Delta t$, which is assumed to be sufficiently small (high sampling
frequency), and so we consider the difference version of the above equations
$(\ref{Eq:coupled dynamic},\ref{Eq:coupled observation})$, by the basic Euler
approach, as follows:

$\displaystyle x(t+1)=F(x(t),\theta,h),\quad y(t)=H(x(t)),$ (3)

with $x(t)=[x_{1}(t),\cdots,x_{m}(t)]$. The question is how to estimate the
hyperparameters from the history of observation sequence
$Y_{t}=\\{y_{t}\\}_{i\in\mathbb{T}}$, that is, how to estimate the posterior
densities $\mathbb{P}(h|Y_{t})$. An intuitive idea under this problem
formulation is we can regard the distribution of the recording dynamic time
series as a map from the distribution, in terms of the parameter of this
distribution named hyperparameter, instead of from the exact values of all
parameters of the model.

Noticed that:

$p(h|Y_{t})=\frac{p(h)}{p(Y_{t})}p(Y_{t}|h)=\frac{p(h)}{p(Y_{t})}\int
p(Y_{t}|\theta,h)p(\theta|h)d\theta.$

We assume that the evolution of the states does not depend the specific
parameters but depend the hyperparameter. That means that there exists a
infinitesimal $\varepsilon(h,m)$ such as $\frac{1}{m}$,
$\forall\theta_{1},\theta_{2}$, satisfying

$\displaystyle\int|p(Y_{t}|\theta_{1},h)-p(Y_{t}|\theta_{2},h)|p(\theta_{1}|h)d\theta_{1}<\varepsilon(h,m),$
(4)

and $\varepsilon(h,m)$ converges to zero almost everywhere when $m$ goes to
infinity. This assumption is not trivial but depends on the property of the
observation function(see the Remark).

Then we have the following approximation using the I.I.D. samples $\theta$

$\displaystyle p(h|Y_{t})\approx\frac{p(h)}{p(Y_{t})}p(Y_{t}|\theta,h),a.e.$
(5)

According to the evolution of the states by the transition probability density
and Markov property, we have

$\displaystyle p(Y_{t}|\theta,h)=$ $\displaystyle\int
p(Y_{t}|X_{t})p(X_{t}|\theta,h)dX_{t}$ $\displaystyle=$ $\displaystyle\int
p(Y_{t}|X_{t})\prod_{t=1}^{T}p(x_{t}|x_{t-1},\theta,h)p(x_{0})dx_{0}\cdots
dx_{t}$ $\displaystyle=$ $\displaystyle\int
p(Y_{t}|X_{t})\prod_{t=1}^{T}p(x_{t}|x_{t-1},\theta)p(\theta|h)p(x_{0})dx_{0}\cdots
dx_{t}$

From the perspective of DA, the hyperparameters $h$ are considered as states
following stochastic differential equation $\frac{d}{dt}h(t)=\sigma dB_{t}$,
which evolve discretely via random walks:

$\displaystyle h(t+1)=h(t)+\zeta(t)$

where $\zeta_{t}$ is a white Gaussian noises with components independent and
small variances $\sigma_{\zeta}^{2}$ which are known preliminarily. The
evolution of the parameters $\theta$ is tricky. Here, we also consider the
following law:

$\displaystyle\theta(t+1)=\Phi(\theta(t),h(t+1),h(t),\eta(t)),$

where $\eta(t)$ are independent (w.r.t time and component) and distributed so
that $\theta(t+1)$ follows the distribution with parameters $h(t+1)$. The
specific realization methods and more details are discussed in next
subsection.

Comparing with the approaches in DA, we consider the following two methods:
ensemble Kalman filter and particle filter. Let write $[x(t),h(t)]$ by $X(t)$
for simplicity. Here we show the EnKF algorithm combined with our HDA as
Algorithm 1 and other details are referred to in S2 Appendix..

Algorithm 1 EnKF combined with HDA

0:
$N,[\mu,C_{0}],F_{0:T}(\cdot),\Sigma_{h},\Sigma_{x},\eta_{0:T},[H,\Gamma],y_{0:T}$.

0: $\\{X_{0:T}^{(n)}\\}_{n=1}^{N}$.

1: draw
$X_{0}^{(n)}\sim\mathcal{N}(\mu_{0},C_{0})),\theta_{0}^{(n)}\sim\mathbb{P}(\cdot|h_{0}^{(n)}),\hat{h}_{0}^{(n)}=h_{0}^{(n)},\qquad\forall
n=1:N$

2: for $t=1:T$ do

3: draw
$\hat{h}^{(n)}_{t}\sim\mathcal{N}(h_{t-1}^{(n)},\Sigma_{h}),\qquad\forall
n=1:N$.

4:
$\theta^{(n)}_{t}=\Phi(\theta^{(n)}_{t-1},\hat{h}_{t}^{(n)},\hat{h}_{t-1}^{(n)},\eta(t)),\qquad\forall
n=1:N$.

5: draw
$\hat{x}^{(n)}_{t}\sim\mathcal{N}(F_{t-1}(x_{t-1}^{(n)},\theta_{t}^{(n)}),\Sigma_{x}),\qquad\forall
n=1:N$.

6: $\hat{\mu}_{t}\leftarrow\frac{1}{N}\Sigma_{n=1}^{N}\hat{X}_{t}^{(n)}$.

7:
$\hat{C}_{t}\leftarrow\frac{1}{N-1}\Sigma_{n=1}^{N}(\hat{X}_{t}^{(n)}-\hat{\mu}_{t})(\hat{X}_{t}^{(n)}-\hat{\mu}_{t})^{T}$.

8: $S_{t}\leftarrow H\hat{C}_{t}H^{T}+\Gamma$.

9: $K_{t}\leftarrow\hat{C}_{t}H^{T}S_{t}^{-1}$.

10: draw $y_{t}^{(n)}\sim\mathcal{N}(y_{t},\Gamma),\qquad\forall n=1:N$.

11: $\delta_{t}^{(n)}\leftarrow y_{t}^{(n)}-H\hat{x}_{t}^{(n)},\qquad\forall
n=1:N$.

12:
$X_{t}^{(n)}\leftarrow\hat{X}_{t}^{(n)}+K_{j}\delta_{t}^{(n)},\qquad\forall
n=1:N$.

13: end for

14: return $\\{X_{0:T}^{(n)}\\}_{n=1}^{N}$.

#### Remark

In our HDA framework, a hierarchical structure is used to replace estimating
the distribution of parameters by estimating hyperparameters. Thus we have not
changed the topology structure of the network. Since we assume that each group
of parameters sampled from the same distribution should have the same system
observation, the parameters initialization of each ensemble sample is
independent in Algorithm 1. And we improve the accuracy of hyperparameters
estimation by increasing the number of parameters rather than the times of
parameters resampled. So it allows us not only to expand the network scale but
also to simulate quickly by using limited computing resources.

### Computational Neuronal network model

The computational neuron model is generally a nonlinear operator from a set of
input synaptic spike trains to an output axon spike train, described by three
components: The subthreshold equation of membrane potential that describes the
transformation from the synaptic currents of diverse synapses; the synaptic
current equation describes the transformation from the input spike trains to
the corresponding synaptic currents; the threshold scheme gives setup of the
condition for triggering a spike by membrane potential.

We, following [5], consider the leakage integrate-and-fire (LIF) model as
neuron. A capacitance-voltage (CV) equation describes the membrane potential
of neuron $i$, $V_{i}$, when it is less than a given voltage threshold
$V_{th,i}$:

$\displaystyle
C_{i}\dot{V}_{i}=-g_{L,i}(V_{i}-V_{L})+\sum_{u}I_{u,i}+I_{ext,i},\qquad
V_{i}<V_{th,i}$ (6)

where $C_{i}$ is the capacitance of the neuron membrane, $g_{L,i}$ is the
leakage conductance, $V_{L}$ is leakage voltage, $I_{u,i}(t)$ stands for the
synaptic current at the synapse type $u$ of neuron $i$ and $I_{ext,i}$ is the
external stimulus from the environment and we set $I_{ext,i}=0~{}nA$ in this
paper. Here, we consider at least four synapse types: $AMPA$, $NMDA$,
$GABA_{A}$ and $GABA_{B}$ [48].

As far the synapse model, We consider an exponentially temporal convolution as
this map:

$\displaystyle\begin{array}[]{ll}I_{u,i}&=g_{u,i}(V_{u,i}-V_{i})J_{u,i}\\\
\dot{J}_{u,i}&=-\frac{J_{u,i}}{\tau_{u,i}}+\sum_{k,j}w_{ij}^{u}\delta(t-t_{k}^{j})\end{array}$
(9)

Here, $g_{u,i}$ is the conductance of synapse type $u$ of neuron $i$, $V_{u}$
is the voltage of synapse type $u$, $\tau_{u,i}$ is the time-scale constants
of synapse type $u$ of neuron $i$, $w_{ij}^{u}$ is the connection weights from
neuron $j$ to $I$ of synapse type $u$, $\delta(\cdot)$ is the Dirac-delta
function and $t_{k}^{j}$ is the time point of the $k-th$ spike of neuron $j$.

When $V_{i}=V_{th,i}$ at $t=t_{n}^{i}$, neuron registers a spike at time point
$t_{n}^{i}$ and the membrane potential is reset at $V_{rest}$ during a
refractory period:

$\displaystyle V_{i}(t)=V_{rest},t\in[t_{k}^{i},t_{k}^{i}+T_{ref}].$ (10)

After then, $V_{i}$ is governed by CV equation again.

Based on the above neuron model, the neuronal network is composed of
$N_{E}=4/5N$ excitatory neurons, denoted by $\mathcal{N}_{E}$ and $N_{I}=1/5N$
inhibitory neurons, denoted by $\mathcal{N}_{I}$. Let
$W_{u}=[w^{u}_{ij}]_{i,j=1}^{N}$ be the adjacent matrix of the synapse type
$u$ where $w^{u}_{ij}\neq 0$ means there is a link from neuron $j$ to $i$ and
the conductance of this link is $w^{u}_{ij}$, otherwise $w^{u}_{ij}=0$. Each
neuron has $D$ neighbors and the network graph topology is generated by the
following rule: each neuron randomly chooses a synaptic junction from $4/5D$
excitatory neurons as excitatory synapses and $1/5D$ inhibitory neurons as
inhibitory synapses. Thus the edge of each $u$-synapse graph is $N\times D$.
For convenience, we transform these ODE into the iterated map by the Euler
method of time step 1 msec and calculate the mean of neuron activity in the
neuronal network as output:

$\displaystyle z(t)=\frac{1}{N}\sum_{i\in
N,k\in\mathbb{N}}\delta(t-t^{i}_{k}).$ (11)

### BOLD signal and hemodynamical model

FMRI is the widely used towards human objects as BOLD signal that indirectly
measures the neural activity in a low spatial and temporal resolution. A
useful mathematical model is the hemodynamical model that takes the neural
activity quantified by the spike rate of a pool of neurons and outputs bold
signal. Let $z_{t}$ be the time series of neural activity from Eq (11) and
hemodynamical model can be generally written as $\dot{g}=G(g,z)$. Here $g$
denotes the states related to blood volume and blood oxygen consumption. And
the bold signal is read as a function of $g$ i.e., $y(t)=h(g(t))$. There are
diverse hemodynamical models and herein, we introduce the Balloon-Windkessel
model [17] as shown:

$\displaystyle\left\\{\begin{array}[]{lr}\dot{s}_{i}=\varepsilon z_{i}-\kappa
s_{i}-\gamma(f_{i}-1)\\\ \dot{f}_{i}=s_{i}\\\
\tau\dot{v}_{i}=f_{i}-v_{i}^{\frac{1}{\alpha}}\\\
\tau\dot{q}_{i}=\frac{f_{i}(1-(1-\rho)^{\frac{1}{f_{i}}})}{\rho}-\frac{v_{i}^{\frac{1}{\alpha}}q_{i}}{v_{i}}\\\
y_{i}=V_{0}[k_{1}(1-q_{i})+k_{2}(1-\frac{q_{i}}{v_{i}})+k_{3}(1-v_{i})]\\\
\end{array}\right.$ (17)

Here, $i$ is the label of ROI/voxel, $s_{i}$ is vasodilatory signal, $f_{i}$
is blood inflow, $v_{i}$ is blood volume and $q_{i}$ is deoxyhemoglobin
content. $\varepsilon,\kappa,\tau,\gamma,\alpha,\rho,V_{0},k_{1},k_{2},k_{3}$
are the constant parameter of the Balloon-Windkessel model. All parameters of
constant term in neuronal network model and hemodynamic model are set in Table
1.

Table 1:  Parameters in Neuronal Network Model and Hemodynamic Model.

Neuronal Network Model  | Hemodynamic Model
---|---
$\bf symbol$ |  $\bf value$ | $\bf symbol$ | $\bf value$
$C_{i}$ | 1$\mu f$ | $\varepsilon$ | 200
$g_{L,i}$ | 0.03$mS$ | $\kappa$ | 1.25
$V_{L}$ | -75$mV$ | $\gamma$ | 2.5
$V_{th}$ | -50$mV$ | $\alpha$ | 0.2
$V_{rest}$ | -65$mV$ | $\rho$ | 0.8
$[V_{AMPA,i},V_{NMDA,i},V_{GABA_{a},i},V_{GABA_{b},i}]$ | [0,0,-70,-100]$mV$ | $V_{0}$ | 0.02
$[\tau_{AMPA,i},\tau_{NMDA,i},\tau_{GABA_{a},i},\tau_{GABA_{b},i}]$ | [2,40,10,50]  | $\tau$ | 1
$w_{i,j}^{u}$ |  $\sim$ rand(0,1)  | $k_{1}$ | 5.6
$T_{ref}$ | 5$ms$ | $k_{2}$ | 2
|  | $k_{3}$ | 1.4

Table shows parameters of constant term in Eq (6), Eq (9), Eq (10) and Eq
(17). $C_{i}$ is the capacitance of the neuron membrane, $g_{L,i}$ is the
leakage conductance, $V_{L}$ is leakage voltage, $V_{th}$ is threshold
voltage, $V_{rest}$ is rest voltage after spike, $V_{u}$ is the voltage of
synapse type $u$, $\tau_{u,i}$ is the time-scale constants of of synapse type
$u$, $u\in\\{AMPA,NMDA,GABA_{a},GABA_{b}\\}$, $w_{i,j}$ is I.I.D from a
uniform distrition over $(0,1)$. Parameters in hemodynamic model are referred
in detail in [17].

Furthermore, to mimic the time series of BOLD signals that were measured by
the fMRI, we conduct a down-sampling process of the recording data:

$\displaystyle TS(l)=y(t_{l})+\eta(l),~{}l=1,2,\cdots,L$

where $\eta(l)$ is a white Gaussian noise with
$cov(\eta(l)\eta(s))=\sigma_{n}\delta(l-s)$ with $\sigma_{n}=10^{-8}$,
$t_{l}$, $l=1,2,\cdots,L$ are sampling time points that are periodic with
frequency $1.25$ HZ, as same as frequency in fMRI.

### Sampling parameters from the hyperparameter

Since our algorithm is based on the framework of Hierarchical Bayesian
Inference, it is necessary to introduce how hyperparameter determines
parameters in the algorithm. As mentioned in section Hierarchical Data
Assimilation, the evolution of the parameters $\theta$ should be tricky
because it should ensure that the parameter follows the distribution with
hyperparameter $h$. Here, we give two methods to adjust the parameter.

The first one is simple and effective by using the cumulative density function
(cdf) and its inverse function. As we can see, the update rule of parameter in
the previous subsection is according to the following equation:

$\displaystyle\tilde{g}_{i}(t+1)=\mathcal{F}_{h(t+1)}^{-1}\circ\mathcal{F}_{h(t)}\left(g_{i}(t)\right),$
(18)

where $\mathcal{F}_{c}(\cdot)$ denotes the cdf by exponential distribution
with mean $c$. Therefore if the parameter follows the specific distribution,
we just use the cdf to map $\theta$ to $\vartheta=\mathcal{F}_{h(t)}(\theta)$,
following uniform distribution and hold random variables $\vartheta$ in this
probability space. Then the map
$\vartheta\rightarrow\hat{\theta}=\mathcal{F}_{h(t+1)}^{-1}(\vartheta)$ ensure
$\hat{\theta}$ following the distribution with new hyperparameter $h(t+1)$.

And the other method, finding a transform rule from the $\theta_{t-1}$ to
$\theta_{t}$ based on the old probability density function(pdf) $p_{0}(\cdot)$
and the new one $q(\cdot)$ according to MCMC method, is to deal with the
situation when the cdf is difficult to be calculated numerically, such as
Gussian distribution or Gamma distribution. Denote the parameters of the last
iteration as $\theta=(\theta_{1},\cdots,\theta_{N})$ and the new parameters as
$\hat{\theta}$. We aim to make the distribution of new parameter as close as
possible to new distribution $\mu_{q}$. By probability theory, the pdf of
$\hat{\theta}$ becomes
$p_{\hat{\theta}}(x)=\Sigma_{n=1}^{N}\frac{1}{N}\mathbb{I}_{\\{x=\hat{\theta}\\}}$.
The algorithm is as follow Algorithm 2.

Algorithm 2 Transformation based on MCMC method

0:
$\theta_{0}=(\theta_{0}^{1},\cdots,\theta_{0}^{N}),p_{0}(\cdot),q(\cdot),T,\Sigma_{0}$

0: $\hat{\theta}$

1: for $t=0:T-1$ do

2: for $n=1:N$ do

3: draw $\theta^{n}_{t+1}\sim\mathcal{N}(\theta^{n}_{t},\Sigma_{0}),$

4: calculate $q_{(\theta^{n}_{t})}$ and $q_{(\theta^{n}_{t+1})}$

5: draw $\mu\sim Union(0,1)$

6: let $\theta^{n}_{t}=\theta^{n}_{t+1},\qquad$ if
$\mu\leqslant\frac{q_{(\theta^{n}_{t+1})}}{q_{(\theta^{n}_{t})}}$

7: end for

8: denote $\theta_{t+1}=(\theta^{1}_{t+1},\cdots,\theta^{N}_{t+1})$

9: end for

10: return $\hat{\theta}=\theta_{T}$

Fig 1(a) illustrates performance of Algorithm 2 at $N=1000$ with $T=1000$ to
generate Gamma distribution $\Gamma(2,1)$. According to this method, it allows
us to choose different distribution of parameters to model the brain network
(See more details in Discussion).

Fig 1: Sampling parameters from the hyperparameter via MCMC Method.

a: These four pictures show the change from exponential distribution to Gamma
distribution during MCMC method iteration and the iteration are 0, 1, 10, 100
respectively. b Based on MCMC method, HDA can also assimilite arbitrarily
distribution of parameters effectively.

## Results

### Asymptotic and well-posedness analyses

In this subsection, we give mathematical proof of the convergence
(asymptotics) of the HDA estimation process in the view of MLE.

Let $Y$ be an $\mathbb{R}^{m}$ -valued random variable, whose probability
distribution is defined by a series of random variable
$\Theta_{m}=[\theta_{1},\cdots,\theta_{m}]$. Meanwhile,
$[\theta_{1},\cdots,\theta_{m}]$ called parameters are i.i.d. following a
specific distribution $\psi(\cdot|h)$ with hyperparameter $h$. Thus the pdf of
$Y$ can be represented by conditional probability:

$\mathbb{P}(y=Y|\Theta_{m},h).$ (19)

It is assumed that $(Y_{1},\cdots,Y_{n})$ called available data is a set of
independent realizations of $Y$, for example coming from observations or
numerical simulations. To simplify the presentation, it is assumed that there
is no noise in the data. We try to prove hyperparameter $h$ can be estimated
as the size of parameters and the times of observation of data goes to
infinity. Since we only have information that parameters are sampled from a
given distribution, it can not get post distribution accurately. so it is
necessary to construct a likelihood function that depends on an unknown
hyperparameter $h$ firstly.

Let $\hat{h}_{n,m}$ be the solution of the following optimization problem,

$\displaystyle\hat{h}_{n,m}=\arg\min L(h,Y_{1},\cdots,Y_{n},\Theta_{m}),$ (20)

in which $Y_{i}$ is sample of state $Y$ and $L(\cdot)$ is the log-likelihood
function defined by

$\displaystyle L(h,Y_{1},\cdots,Y_{n},\Theta_{m})$
$\displaystyle=-\frac{1}{N}\sum_{i=1}^{N}L(h,Y_{i},\Theta_{m}^{i})$
$\displaystyle=-\frac{1}{N}\sum_{i=1}^{N}\log\mathbb{P}(Y_{i}|h,\Theta_{m}^{i})$
(21)

If the log-likelihood function is differentiable and $\hat{h}_{n,m}$ is a
global maximum of the optimization problem, we have

$\frac{\partial}{\partial
h}L(h,Y_{1},\cdots,Y_{n},\Theta_{m})|_{h=\hat{h}_{n,m}}=0.$ (22)

For the equation 22, we have a Taylor expansion for $h$ in the neighborhood
$\mathcal{N}(\tilde{h})$ about the point $h=\tilde{h}$, we assume
$h\in{\mathbb{R}}$ for simplify, as follow

$\displaystyle\frac{\partial}{\partial h}L(h,Y_{1},\cdots,Y_{n},\Theta_{m})$
$\displaystyle=\frac{\partial}{\partial
h}L(h,Y_{1},\cdots,Y_{n},\Theta_{m})|_{h=\tilde{h}}$ (23)
$\displaystyle+(h-\tilde{h})\frac{\partial^{2}}{\partial
h^{2}}L(h,Y_{1},\cdots,Y_{n},\Theta_{m})|_{h=\tilde{h}}$
$\displaystyle+\frac{1}{2}\zeta(h-\tilde{h})^{2}H(Y_{1},\cdots,Y_{n},\tilde{h}),$

where $\zeta\in(0,1)$ and for all $Y$,

$\displaystyle|\frac{\partial^{3}}{\partial
h^{3}}L(h,Y_{1},\cdots,Y_{n},\Theta_{m})|\leqslant H(Y_{1},\cdots,Y_{n},h).$

We assume:

  * R1.

For each $h\in\mathcal{N}(\tilde{h})$,

$\displaystyle\frac{\partial}{\partial
h}L(h,Y,\Theta_{m}),\frac{\partial^{2}}{\partial
h^{2}}L(h,Y,\Theta_{m}),\frac{\partial^{3}}{\partial h^{3}}L(h,Y,\Theta_{m}),$

exist, for all $Y$.

  * R2.

For each $h\in\mathcal{N}(\tilde{h})$, there exist functions $f(Y),g(Y),H(Y)$
(possibly depending on h) such that the relations

$\displaystyle|\frac{\partial}{\partial h}\mathbb{P}(Y|h,\Theta_{m})|\leqslant
f(Y),|\frac{\partial^{2}}{\partial h^{2}}\mathbb{P}(Y|h,\Theta_{m})|\leqslant
g(Y),|\frac{\partial^{3}}{\partial h^{3}}L(h,Y,\Theta_{m})|\leqslant H(Y)$

a.s. hold, for all $Y$, and

$\displaystyle\int f(Y)dY<\infty,\int
g(Y)dY<\infty,\mathbb{E}_{Y}[H(Y)]<\infty.$

  * R3.

For all $\Theta_{m}$ and $\Phi_{m}$ following the distribution $\psi(\cdot|h)$
with same hyperparameter $h$, then there exist $k,K,K_{\infty}$ satisfying

$\displaystyle 0<k\leqslant\int[\frac{\partial}{\partial
h}L(h,Y,\Theta_{m})]^{2}\mathbb{P}(Y|h,\Phi_{m})dY\leqslant K<\infty,$
$\displaystyle a.s.$ $\displaystyle 0<k\leqslant\int[\frac{\partial}{\partial
h}L(h,Y,\Theta_{m})]^{4}\mathbb{P}(Y|h,\Phi_{m})dY\leqslant K<\infty,$
$\displaystyle a.s.$
$\displaystyle\lim_{m\rightarrow\infty}\int[\frac{\partial}{\partial
h}L(h,Y,\Theta_{m})]^{2}\mathbb{P}(Y|h,\Theta_{m})dY=K_{\infty},$
$\displaystyle a.s.$

  * R4.

The Hellinger distance between two distributions is defined by

$\displaystyle
d_{Hell}(\mathbb{P}_{1},\mathbb{P}_{2})=\big{(}\frac{1}{2}\int(\sqrt{p_{1}(y)}-\sqrt{p_{2}(y)})^{2}dy\big{)}^{\frac{1}{2}}.$

We assume there exists a distribution $\mathbb{P}(Y|h)$ satisfies that for
arbitrary distribution $\mathbb{P}(Y|h,\Theta_{m})$,

$\displaystyle
d_{Hell}(\mathbb{P}(Y|h,\Theta_{m}),\mathbb{P}(Y|h))\leqslant\varepsilon(m,h)\xrightarrow[m\rightarrow\infty]{a.s.}0,$

i.e. $\varepsilon(m,h)$ converges to zero almost everywhere when m goes to
infinity. Thus, for arbitrary two distributions $\mathbb{P}(Y|h,\Theta_{m})$
and $\mathbb{P}(Y|h,\Phi_{m})$, we have

$\displaystyle
d_{Hell}(\mathbb{P}(Y|h,\Theta_{m}),\mathbb{P}(Y|h,\Phi_{m}))\leqslant
2\varepsilon(m)\xrightarrow[m\rightarrow\infty]{a.s.}0.$

When $\\{Y_{i}\\}$ drawn from pdf $\mathbb{P}(\cdot|\tilde{\theta},\tilde{h})$
with hyperparameter $\tilde{h}$ and parameter
$\tilde{\theta}=[\tilde{\theta}_{1},\cdots,\tilde{\theta}_{m}]$ which are the
real hyper-parameter and parameters of the specific dynamic system, Theorem 1
shows $\tilde{h}$ and $\hat{h}_{n,m}$ are very close to each other as
$m,n\rightarrow\infty$ and the prove is in S3 Appendix..

###### Theorem 1

Assume conditions (R1-R4), the solution of equation (22) $\hat{h}_{n,m}$ and
the real hyperparameter $\tilde{h}$ satisfy that $|\tilde{h}-\hat{h}_{n,m}|$
converges to zero as $m,n$ goes to infinity.

On the other hand, if we have two different simulation of the data,
$Y^{{}^{\prime}}$ and $Y^{{}^{\prime\prime}}$ respectively, the difference of
hyperparameter estimator based on different observations,
$|\tilde{h}^{{}^{\prime}}-\tilde{h}^{{}^{\prime\prime}}|$, can be controlled
by the $|Y^{{}^{\prime}}-Y^{{}^{\prime\prime}}|$ [29]. It means our estimator
given by MLE depends continuously on the observation data, shown in Theorem 2
and the prove is in S4 Appendix..

###### Theorem 2

Let $x$ denote all model states including hyperparameter and Bold model
states, $\mu$ and $\mu^{\prime}$ be the posterior probability distributions of
model states associated with two different data sets $Y=y_{1:\mathbb{T}}$ and
$Y^{{}^{\prime}}=y^{\prime}_{1:\mathbb{T}}$ respectively, $\mu_{0}$ be the
prior measure on $x$. Assume that
$\mathbb{E}_{\mu_{0}}w(x_{1:\mathbb{T}})<\infty$, where $w(x_{1:\mathbb{T}})$
is given by

$\displaystyle
w(x_{1:\mathbb{T}})=\Sigma_{t=1}^{\mathbb{T}}(1+|H(x_{t})|^{2}).$ (24)

Then there exists $c=c(r)$ such that for all
$|y_{1:\mathbb{T}}|,|y^{\prime}_{1:\mathbb{T}}|<r$,

$\displaystyle d_{Hell}(\mu,\mu^{\prime})\leqslant
c|y_{1:\mathbb{T}}-y^{\prime}_{1:\mathbb{T}}|.$ (25)

#### Remark

In this paper, we use HDA framework to simulate the hyperparameters of
ordinary differential equations with assumption 4, and here we will show the
assumption can be achieved by a specific observation function of the system.
Assumption 4 essentially represents not each In this paper, we use HDA
framework to simulate the hyperparameters of ordinary differential equations
with assumption 4, and here we will show the assumption can be achieved by a
specific observation function of the system. Assumption 4 essentially
represents not each parameters but the hyperparameters are vital to the
observation series $y_{t}$. And a necessary condition is that as the number of
parameters $m$ goes to infinite, each parameter $\theta$ in ODE has little
effect to observation, such as $\lim_{m\rightarrow\infty}\frac{\partial
y_{t}}{\partial\theta}=0$. It is obvious that choosing the mean value function
of all states as an observation function satisfies the condition and here we
give a simple model.

Let us consider one of the simplest case of discrete time. $x_{i}(t)$ be an
i.i.d. Gaussian process: $x_{i}(t)=\mu_{i}+\varepsilon_{i,t}$, $i=1,\cdots,m$,
and $t=1,\cdots,T$, where $\mu_{i}\in{\mathbb{R}}$, $\sigma_{i}>0$ and
$\varepsilon_{i,t}$ is an standard white Gaussian process which are
independent with respect to $i$ and $t$. Suppose that all $\mu_{i}$ are
independently identically Gaussian distributed with $N(\mu_{0},1)$. The
observation $y(t)=\frac{1}{m}\sum_{i=1}^{m}x_{i}(t)$ is the average of all
$x_{i}$. After a simple calculation, we have
$y(t)=\frac{1}{m}\sum_{i=1}^{m}(\mu_{i}+\varepsilon_{i,t})\sim
N(\mu_{0},\frac{2}{m})$, $\lim_{m\rightarrow\infty}\frac{\partial
y(t)}{\partial\mu_{i}}=0$, $\mathbb{E}y(t)=\mu_{0}$. And if $\mu$ and
$\mu^{\prime}$ are both i.i.d. with $N(\mu_{0},1)$,
$\lim_{m\rightarrow\infty}|P(y^{\prime}|\mu^{\prime},\mu_{0})-P(y|\mu,\mu_{0})|=0$.
Thus the observation satisfies the condition4. From this simple model, we can
also find that if the task is to estimate hyperparameter $\mu_{0}$, larger $m$
makes $y$ and $\mu_{0}$ more closer. In other words, the estimation statistics
of the hyperparameter converges to the real value (with respect to the
invariant measure) as $m$ (number of equations/states) goes to infinity. The
hyperparameters are vital to the observation series $y_{t}$. And a necessary
condition is that as the number of parameters $m$ goes to infinite, each
parameter $\theta$ in ODE has little effect to observation, such as
$\lim_{m\rightarrow\infty}\frac{\partial y_{t}}{\partial\theta}=0$. It is
obvious that choosing the mean value function of all states as an observation
function satisfies the condition and here we give a simple model.

Let us consider one of the simplest case of discrete time. $x_{i}(t)$ be an
i.i.d. Gaussian process: $x_{i}(t)=\mu_{i}+\varepsilon_{i,t}$, $i=1,\cdots,m$,
and $t=1,\cdots,T$, where $\mu_{i}\in{\mathbb{R}}$, $\sigma_{i}>0$ and
$\varepsilon_{i,t}$ is an standard white Gaussian process which are
independent with respect to $i$ and $t$. Suppose that all $\mu_{i}$ are
independently identically Gaussian distributed with $N(\mu_{0},1)$. The
observation $y(t)=\frac{1}{m}\sum_{i=1}^{m}x_{i}(t)$ is the average of all
$x_{i}$. After a simple calculation, we have
$y(t)=\frac{1}{m}\sum_{i=1}^{m}(\mu_{i}+\varepsilon_{i,t})\sim
N(\mu_{0},\frac{2}{m})$, $\lim_{m\rightarrow\infty}\frac{\partial
y(t)}{\partial\mu_{i}}=0$, $\mathbb{E}y(t)=\mu_{0}$. And if $\mu$ and
$\mu^{\prime}$ are both i.i.d. with $N(\mu_{0},1)$,
$\lim_{m\rightarrow\infty}|P(y^{\prime}|\mu^{\prime},\mu_{0})-P(y|\mu,\mu_{0})|=0$.
Thus the observation satisfies the condition4. From this simple model, we can
also find that if the task is to estimate hyperparameter $\mu_{0}$, larger $m$
makes $y$ and $\mu_{0}$ more closer. In other words, the estimation statistics
of the hyperparameter converges to the real value (with respect to the
invariant measure) as $m$ (number of equations/states) goes to infinity.

### Simulation of LIF neuronal network

We use our HDA Algrothms to simulate LIF neuron network mentioned in Methods
of $N=1000$ with degree $D=20$ up to $T=40000~{}msec$ by the Euler method of
time step $1~{}msec$. Thus, we have a BOLD-signal time series of $50$ time
points. The brain region activity is recorded as the sum of neuron activity
for all $i\in N$:

$\displaystyle z(t)=\frac{1}{N}\sum_{i\in
N,k\in\mathbb{N}}\delta(t-t_{i}^{k})+\eta(t).$

where $\eta(t)$ is a white Gaussian noise with
$cov(\eta(t)\eta(s))=\sigma_{n}\delta(t-s)$ with $\sigma_{n}=10^{-8}$. After
getting brain region activity, we use Balloon-Windkessel model mentioned above
to generate the BOLD signal. Then we focus on how to estimate the
hyperparameters $h$ of synaptic conductance $g_{ui}$ from the observation
series $\\{y(t_{l})\\}_{l=1}^{T}$. To estimate these hyperparameters,
parameters of synaptic conductance is sampled from the exponential
distribution with the given hyperparameter independently and the procedure(see
in S2 Appendix.) is similar to Algorithm 1. When we choose one of the
hyperparameters to assimilate, such as $g_{AMPA,i}$ $i\in\mathcal{N}$, all
other parameters including $g_{u,i}(u\neq AMPA)$ are assumed to be known. All
other neuron network model parameters are recapitulated in Table 1 and the
details of algorithm configuration are in Discussion.

Fig 2(a, b) illustrates BOLD error and hyperparameter error when we choose
different combinations of hyperparameters to assimilate. The main criterion is
how close the BOLD signal generated by toy-model with specific hyperparameters
$h$ is to the BOLD signal simulation with hyperparameters simulation under
specify configuration of algorithm parameters, as follow:

$\displaystyle\varepsilon_{BOLD}$
$\displaystyle=\frac{1}{T}\sum_{t=1}^{T}\frac{|BOLD_{simulation}(t)-BOLD_{real}(t)|^{2}}{|BOLD_{real}(t)|^{2}}.$
$\displaystyle\varepsilon_{h}$
$\displaystyle=\frac{1}{T}\sum_{t=1}^{T}\frac{|h_{simulation}(t)-h_{real}(t)|^{2}}{|h_{real}(t)|^{2}}.$

As shown, the best performance in single hyperparameter assimilation is
$u=AMPA$, and in the conbination of excitatory and inhibitory synaptic
conductance hyperparameters assimilation is $u=NMDA,GABA_{b}$.

Fig 2: BOLD error and hypermeter error under different combinations of
hyperparameter.

a, b: The difference of BOLD singles and hyperparameters between simulation
and real are shown respectively on single synaptic conductance hyperparameter
assimilation or the conbination of excitatory and inhibitory synaptic
conductance hyperparameters assimilation. The error funcition is defined in
above. c: Sliding window is used to calculate the fire rate of sampling
neurons and 40 percent of neurons have a firing rate of less than 5 Hz. d:
Although the observation is at the level of second, the spike of sample
neurons can be recorded at each millisecond based on our model. e: Bold signal
and hyperparmeter serie of real data are both close to simulation under AMPA
synaptic conductance assimulation.

### Performance of HDA

Here we show the basic behavior of the network during a particular
assimilation trial, $AMPA$ synaptic conductance assimilation. We record the
spikes of all neurons at each millisecond, use a sliding window to average the
spike every 800 msec, and calculate the spike frequency distribution shown in
Fig 2(c). It is easy to see that spike sequence and frequency distribution are
similar to the record from the human brain, which shows the availability of
our neuron network and method from the side. Besides, we choose 200 neurons to
record their activities over 1000 milliseconds and we can see there is no
synchronization or excessive activation in sample neurons in Fig 2(d).
Besides, in Fig 2(a, b), it is shown that $AMPA$ synaptic conductance
assimilation has the best performance which means the simulation signal is
most close to the real Bold signal. We show Bold signal and hyperparameter
series of HDA simulation and real data under $AMPA$ synaptic conductance
assimilation in Fig 2(e).

### Using MCMC to sample parameters from the hpyerparameter

The Algorithm 2 (in Methods) makes it be realized to assimilate arbitrary
parameters distribution with known probability density function. Without
losing generality, we choose Gamma distribution instead of exponential
distribution as a sample whose cumulative distribution function is hard to
calculate. It is assumed that the parameters of AMPA synaptic conductance is
i.i.d Gamma distribution with $\alpha=2$ and the original distribution in HDA
is Gamma distribution with $\alpha=1$. As shown in Fig 1b, with more iteration
times, HDA method still assimilates the real hyperparameter and thus is still
valid.

## Discussion

### Compare the setup of the algorithms

In this section, we explore that how different learning parameters setup (in
Table 2) in our method and different structures of network in out model
influence the simulation result by experimental data. The main criterion is
how close the BOLD signal generated by toy-model is to the BOLD signal
simulation under specify configuration of different parameters and how close
the estimation of hyperparameter is to the real hyperparameter.

Table 2:  Compare the different configuration of algorithm.

Different configuration | Standard deviation of | Standard deviation of | Number of | Parameter of
---|---|---|---|---
of HDA Algorithm 1 | hyperparameter noise | observation noise | the samples | WS network
| configuration | 0.01 | $10^{-4}$ | 5 | 0
Config1 | hp error $(\times 10^{-3})$ | $4.66\pm 2.85$ | $683\pm 265$ | $112\pm 195$ | $55.2\pm 65.3$
| bold error $(\times 10^{-3})$ | $6.25\pm 2.87$ | $184\pm 102$ | $29.8\pm 31.5$ | $3.18\pm 0.02$
| configuration | 0.1 | $5\times 10^{-4}$ | 10 | 0.1
Config2 | hp error $(\times 10^{-3})$ | $4.89\pm 3.27$ | $5.94\pm 1.18$ | $26.7\pm 58.4$ | $63.7\pm 71.5$
| bold error $(\times 10^{-3})$ | $5.70\pm 2.58$ | $5.46\pm 0.69$ | $9.92\pm 11.5$ | $3.23\pm 0.02$
| configuration | 0.2 | $\mathbf{10^{-3}}$ | 50 | 0.2
Config3 | hp error $(\times 10^{-3})$ | $4.90\pm 1.97$ | $\mathbf{5.29\pm 0.79}$ | $5.71\pm 1.95$ | $67.4\pm 65.2$
| bold error $(\times 10^{-3})$ | $4.87\pm 1.07$ | $\mathbf{4.39\pm 0.21}$ | $4.59\pm 0.65$ | $3.21\pm 0.02$
| configuration | $\mathbf{0.5}$ | $5\times 10^{-3}$ | $\mathbf{100}$ | 0.3
Config4 | hp error $(\times 10^{-3})$ | $\mathbf{5.29\pm 0.79}$ | $4.18\pm 0.61$ | $\mathbf{5.29\pm 0.79}$ | $67.9\pm 55.5$
| bold error $(\times 10^{-3})$ | $\mathbf{4.39\pm 0.21}$ | $4.10\pm 0.26$ | $\mathbf{4.39\pm 0.21}$ | $3.20\pm 0.03$
| configuration | 1 | $10^{-2}$ | 200 | 0.5
Config5 | hp error $(\times 10^{-3})$ | $5.85\pm 0.85$ | $2.97\pm 0.64$ | $5.19\pm 0.87$ | $74.5\pm 58.5$
| bold error $(\times 10^{-3})$ | $4.78\pm 0.35$ | $5.04\pm 0.61$ | $4.32\pm 0.11$ | $3.18\pm 0.02$
| configuration | 2 | $5\times 10^{-2}$ | 300 | 0.8
Config6 | hp error $(\times 10^{-3})$ | $8.24\pm 1.51$ | $7.72\pm 1.87$ | $5.47\pm 0.45$ | $87.3\pm 98.7$
| bold error $(\times 10^{-3})$ | $6.65\pm 0.72$ | $23.4\pm 4.18$ | $4.25\pm 0.07$ | $3.03\pm 0.01$
| configuration | 5 | $10^{-1}$ | 500 | 1
Config7 | hp error $(\times 10^{-3})$ | $109\pm 142$ | $13.2\pm 2.34$ | $5.36\pm 0.37$ | $69.7\pm 102$
| bold error $(\times 10^{-3})$ | $70.6\pm 87.8$ | $54.9\pm 8.75$ | $4.28\pm 0.09$ | $3.12\pm 0.02$

Table shows different configurations of algorithm parameters discussed in this
section to exhibit the error simply and intuitively. The highlight of the
table is the final configurations used in Results.

#### Standard deviation of noise

When we use Enkf to estimate the state of the real system or ground truth
model, we regard the hyperparameters as system states as mentioned before, and
therefore it is necessary to add Gaussian noise in the model samples and
observation dataset. It is obvious that bigger covariance of noise causes more
randomness in the learning process, thus we should choose proper covariance to
make sure the learning efficiency. Table 2 shows the hyperparameter error and
bold error under different covariance of model noise $\sigma_{h}$ and
different covariance of observation noise $\sigma_{o}$. In the previous
section, we set $\sigma_{h}=1,\sigma_{o}=10^{-6}$ in Fig 2(e) and this
configuration is best according to the Table LABEL:Table:2, considering
hyperparameters and BOLD error both.

#### Number of the samples

In the traditional Enkf method, it is necessary to use ensemble samples to
calculate the covariance of states. Thus in our algorithm, ensemble samples
are also used and it is obvious that more samples cause more accuracy during
simulating bold signals and estimating the hyperparameters. Table 2 shows the
hyperparameters and BOLD error under the different sizes of samples. In our
trial, we use a NVIDIA Telsa V100 GPU to simulate 20 samples of 1000 neurons
with network degree 20, and 5 GPU to realize 100 samples as reasonable
ensemble number defaults.

### Effect of network topologies

In the previous part of this paper, we assumed the neuron network is a random
network with a specific degree. For one thing, the random network can improve
randomness to make assumption 4 easily realized, for another, the random
network can reduce the difference between different samples in EnKF. Here we
try other network structures such as WS-small-world network [45] and BA-scale-
free network [3]. In WS-small-world network, the parameter of the network
means the probability of rewiring each edge, thus larger parameter means the
network has a more random connection. In BA-scale-free network, the number of
edges to attach from a new node to existing nodes is set to 10. As shown in
Table 2 and Fig 3(a), BOLD can be simulated with small error but
hyperparameters simulation is worse than random network.

Fig 3: BOLD error and hypermeter error under different combinations of
hyperparameter.

a: When generating BA-scale-free network, the new nodes are regard as
postsynaptic neurons so that most neurons have same number of parent nodes. In
this case, the observation is greatly affected by hub nodes thus
hyperparameter is hard to estimate accurately. b: We replace ensemble Kalman
filter with particle filter to assimalate the AMPA synaptic conductance. The
other configurations of Algorithm remain same and the numbers of samples in
two trials are both 100.

### HDA by particle filter

As we showed in Methods, except Kalman filter, HDA framework can also combine
with other filters such as particle filter. Different from Kalman filter
trying to minimize mean-squared posterior error, particle filter seeks to
maximize a posterior probability and the specified procedure is introduced in
detail in S2 Appendix.. Because of the complexity of the neuron model and the
scale of network, particle filter usually needs more samples to calculate
posterior probability and thus needs more storage space and computation. Here
we apply particle filter with 100 samples to make a comparison to EnKF filter.
As shown in Fig 3(b), particle filter also has a good effect in both
hyperparameters estimation and BOLD signal simulation. But when the dimension
of the parameter is high, the number of the samples increases exponentially to
ensure the simulation accuracy.

### Transformation function of hyperparameters

As shown in S1 Appendix., HDA framework adds hyperparameters into ODE states
by following the equation, $\dot{h}=0$. However, if we use $h^{\prime}=f(h)$
as new hyperparameters, $h^{\prime}$ should also satisfy
$\dot{h^{\prime}}=\frac{\partial f}{\partial h}\dot{h}=0$. Thus,
theoretically, we can choose proper hyperparameters representation by adding
constraints to get better simulation results. In this paper, we replace $h$
with $h^{\prime}$ to have a random walk, where

$\displaystyle h=h_{min}+\frac{(h_{max}-h_{min})\exp(\lambda
h^{\prime})}{1+\exp(\lambda h^{\prime})},$ (26)

$\lambda=0.1$, $h_{min}$ and $h_{max}$ are one quarter and double of ground
truth of synaptic conductance respectively. In this situation, we make sure
$h$ belong to $(h_{min},h_{max})$, when random walk of $h^{\prime}$ is on
${\mathbb{R}}$.

## Conclusion

When focusing on simulation or estimation of states in large-scale network
systems, the number of model parameters is frequently much larger than the
number of the observed samples. Especially when investigating the human brain,
it is much more difficult to measure each single neuron activity but able to
collect the BOLD signals of voxels or brain regions. In this paper, we give an
example (in Discussion) that the observation data does not depend on the
specific parameters but depend on the distribution of the parameters. Thus we
claim that, in many cases, the distribution of the parameters in the model
rather than the exact values of these parameters are physiologically
interesting in a very large-scale complex network dynamical system. By
describing the distribution of random variables by introducing a
hyperparameterized distribution, we proposed an HDA framework by estimating
hyperparameter based on hierarchical Bayesian inference. According to the
trails in previous sections, the HDA framework can estimate the
hyperparameters with the true values well with small errors, and the error of
the sampled BOLD signals with the true simulated data is small, too. This
result is shown to be stable and robust for different parameter resamples from
the estimated hyper-parameter.

In detail, we design a tricky evolution of parameters to guarantee that the
parameters are independent and identically distributed throughout. When the
cumulative distribution function (cdf) of the parameters is easy to be
numerically calculated, as we have shown in Methods, the evolution is achieved
by the cdf and its inverse function. However, when the cdf is complicated to
be numerical calculated, Algorithm 2 is an effective method based on the MCMC
method. We arise two HDA algorithms according to the particle filter and the
ensemble Kalman filter (EnKF) approaches in data assimilation. These methods
are both proved efficient for a large-scale neuron network model with the
simulated BOLD signals. Furthermore, we compare the BOLD signals simulation
brought out by different configurations of our algorithm and different
parameters in our neuron model. Estimating synaptic conductance of $AMPA$ has
the best performance in single hyperparameter assimilation, and estimating
both $NMDA$ and $GABA_{b}$ are the best in the combination of excitatory and
inhibitory synaptic conductance assimilation. For network topology, BOLD
signals can be simulated well in all of WS-small-world network, BA-scale-free
network and random network, while hyperparameters estimation is more accurate
in a random network. Finally, we prove that the HDA algorithm can
asymptotically estimate the hyperparameter feasibly under certain assumptions.

We conclude that, the proposed HDA framework present an efficient method to
statistically infer the large-scale full-brain neuronal network with
inhomogeneous parameters, since the number of hyperparameter is relatively
small that can avoid over-fitting results. Last but not least, the HDA is an
extensible framework to equip with diverse biophysical neuronal network model
other than LIF network, and diverse modules of brain data other than MRI. This
could be one of our future research orients.

## Supporting information

##### S1 Appendix.

Overview of Data Assimilation.
Let $F(\cdot)\in C^{1}(\mathbb{R}^{n},\mathbb{R}^{n})$, $H(\cdot)\in
C^{1}(\mathbb{R}^{n},\mathbb{R}^{q})$ and consider the following general
ordinary differential equation:

$\displaystyle\frac{d}{dt}x(t)=$ $\displaystyle F(x(t)),$ (27) $\displaystyle
y(t)=$ $\displaystyle H(x(t)),$ (28)

where $x\in\mathbb{R}^{n}$ stands for the state vectors of the dynamic system
and the observation of $x$ is $y\in\mathbb{R}^{q}$. Here $H(\cdot)$ is called
a observation map.

When $F$ is so complicated that we fail to get accurate formula solution, it
is effective to numerical simulation by transforming the ODE into the iterated
map, or discrete-time dynamical system. Due to numerical errors and modeling
errors, we are also describe the physical system by the following discrete-
time stochastic dynamical model of the form:

$\displaystyle x(t)$ $\displaystyle=F(x(t-1))+\xi_{t},$ (29) $\displaystyle
y(t)$ $\displaystyle=H(x(t))+\eta_{t},\quad t\in\mathbb{T}=\\{1,\cdots,T\\}.$
(30)

Here $\xi=\\{\xi_{t}\\}_{t\in\mathbb{T}}$ is an i.i.d.sequence, independent of
$x$ and $\eta=\\{\eta_{t}\\}_{t\in\mathbb{T}}$ is an i.i.d.sequence,
independent of $(x,y,\xi)$.

If observational data $Y_{T}=\\{y(t)\\}_{t\in\mathbb{T}}$ is acquired, DA is a
kind of iterative process by correcting the states of the dynamical model
based on the observation and model information and can be derived using
Bayesian statistics and recursive Bayesian estimation. There are generally two
classes of DA methods. The first one is filtering methods which update the
posterior filtering distribution $\mathbb{P}(x_{t}|Y_{T})$ at each time and
the second one is smoothing methods which aim at approximating the posterior
smoothing distribution $\mathbb{P}(X_{T}|Y_{T})$. Generally speaking, DA
methods is optimization to seek the posterior distribution.DA methods is
optimization to seek the posterior distribution and the criteria for measuring
the optimality can be referred in [book Guzzi]:

  1. 1.

The Minimum Mean-Squared Error (MMSE) that is defined in terms of prediction
error

$\displaystyle\mathbb{E}[||x_{t}-\hat{x}_{t}||^{2}|y_{\mathbb{T}}]=\int||x_{t}-\hat{x}_{t}||^{2}p(x_{t}|y_{\mathbb{T}})dx_{k},$

where the conditional mean $\hat{x}_{t}=\mathbb{E}[x_{t}|y_{\mathbb{T}}]=\int
x_{t}\ p(x_{t}|y_{\mathbb{T}})dx_{t}$.

  2. 2.

Maximum A Posteriori (MAP) that is aimed to find the mode of posterior
probability $p(x_{t}|y_{\mathbb{T}})$ which is equal to minimize a cost
function.

  3. 3.

The Maximum Likelihood (ML) neglect the prior, that is to find maximum
likelihood estimator by maximize likelihood function
$L=\mathbb{P}(y_{\mathbb{T}}|x_{\mathbb{T}})$.

  4. 4.

The Minimax that is to minimize the maximum of cost function of the posterior
probability.

In discrete-time stochastic dynamical model, if we assume that $F(\cdot)$ and
$H(\cdot)$ are determined by the parameter vector $\theta$ and $\theta$ is
regarded as a state variable with dynamic equation $\dot{\theta}=0$, DA has
proved to be an efficient framework in parameter estimation.

##### S2 Appendix.

EnKF and Particle filter.
Ensemble Kalman filter as a kind of Kalman filter derivatives tries to
minimize mean-squared error. EnKF is a approximate Guassian filter proposed by
G. Evensen, that has a natural generalization to non-Gaussian problems.
Compared with Kalman Filter approach, EnKF draws some samples of the dynamic
system and uses the covariance of samples to substitute for the covariance of
posterior distribution. The complete EnKF procedure is as follows.

  1. 1.

Generate $N$ independent realisations of all random terms in the model,
including the Gaussian noised/error in measurements and random walks, denoted
by $\alpha^{n}(t)$, $n=1,\cdots,N$;

  2. 2.

Generate $N$ realisations of the predictive density $p(X(t)|Y_{t-1})$:
$\hat{X}^{n}(t)$, $j=1,\cdots,N$, by the map
$\hat{X}^{n}(t)=F(\tilde{X}^{n}(t-1),\alpha^{n}(t))$;

  3. 3.

Calculate the mean $\mu(t)=\frac{1}{N}\sum^{N}_{n=1}\hat{X}^{n}(t)$

  4. 4.

Calculate the covariance matrix
$C(t)=\frac{1}{N-1}\sum^{N}_{n=1}(\hat{X}^{n}(t)-\mu(t))(\hat{X}^{n}(t)-\mu(t))^{T}$;

  5. 5.

Calculate the Kalman gain matrix $K(t)=C(t)H^{T}(HC(t)H^{T}+\Gamma)^{-1}$;

  6. 6.

Generate $N$ realizations of the filter density $p(X(t)|Y_{t})$:
$\tilde{X}^{n}(t)=\hat{X}^{n}(t)+K(t)(Y_{t}-HX^{n}(t))$, $n=1,\cdots,N$.

From the view of minimization principle [29], we note the posterior filter
distribution proportional to $\exp(-I_{filter})$. While the Kalman filter is
restricted to linear Gaussian problem, the mean in Kalman filter can be
written as:

$\displaystyle m_{j+1}$ $\displaystyle=\arg\min_{v}I_{KF}(v)$ $\displaystyle
I_{KF}(v)$
$\displaystyle=\frac{1}{2}|y_{j+1}-Hv|_{\Gamma}^{2}+\frac{1}{2}|v-Fm_{j}|_{\hat{C}_{j+1}}^{2},$

where $F\in R^{n\times n},H\in R^{m\times n}$ are linear operator. When we
consider the nonlinear case, the mean in EnKF is given by

$\displaystyle m_{j+1}$ $\displaystyle=\arg\min_{v}I_{EnKF}(v)$ (31)
$\displaystyle\hat{\mu}_{j+1}$ $\displaystyle=F(m_{j})$ $\displaystyle
I_{EnKF}(v)$
$\displaystyle=\frac{1}{2}|y_{j+1}-H(v)|_{\Gamma}^{2}+\frac{1}{2}|v-\hat{\mu}_{j+1}|_{\hat{C}_{j+1}}^{2}.$

Particle filter approach seeks maximum a posteriori. By Bayesian and total
probability formulas, we have

$\displaystyle p(X(t)|Y_{t})$
$\displaystyle=p[X(t)|(y(t),Y_{t-1})]=\frac{p(y(t)|X(t)|Y_{t-1})p(X(t)|Y_{t-1})}{p(y(t)|Y_{t-1})}$
$\displaystyle=\frac{p(y(t)|X(t))}{p(y(t)|Y_{t-1})}p(X(t)|Y_{t-1})$ (32)

noting that $p(y(t)|X(t)|Y_{t-1})=p(y(t)|X(t))$ because $y(t)$ is an
observation of $X(t)$, with

$\displaystyle p(y(t)|Y_{t-1})=\int p(y(t)|X(t))p(X(t)|Y_{t-1})dX(t)$ (33)

Monto-Carlo filter is deployed here to construct the posterior density. In
this approach, each density is approximated by the sampling data. For example,
for a probability density $p(\cdot)$, let $z_{1},\cdots,z_{N}$ be $N$
independent samples from this probability distribution. Then, we have

$\displaystyle p(z)\approx\frac{1}{N}\sum_{i=1}^{N}g(z-z_{i})$

where $g(\cdot)$ is a high peak function centered at $0$, for example, the
Gaussian-type function: $g(z)\propto\exp(-z^{2}/(2\alpha_{g}^{2}))$ with a
sufficiently small $\alpha_{g}$.

Note

$\displaystyle p(X(t)|Y_{t-1})=\int p(X(t)|X(t-1))p(X(t-1)|Y_{t-1})dX(t-1).$
(34)

Beside the estimation problem $p(X(t)|Y_{t})$, which is also named filter, the
prediction problem is also concerned, i.e., $p(x_{t}|Y_{t-1})$, and related,
as indicated in (34).

Giving the size of sampling as $N$, which is very large, the following
bootstrap algorithm is used. Initiated with the realizations from the prior
distribution, the iteration from $p(X(t-1)|Y_{t-1})$ to $p(X(t)|Y_{t})$ is
conducted as follows:

  1. 1.

Generate $N$ independent realisations of all random terms in the model,
including the Gaussian noised/error in measurements and random walks, and the
Poisson point processes at synapses, denoted by $\hat{\alpha}^{j}(t)$,
$j=1,\cdots,N$;

  2. 2.

Generate $N$ realisations of the predictive density $p(X(t)|Y_{t-1})$:
$\hat{X}^{j}(t)$, $j=1,\cdots,N$, by the map
$\hat{X}^{j}(t)=F(\tilde{X}^{j}(t-1),\hat{\alpha}^{j}(t))$;

  3. 3.

Calculate the density weights $\omega^{j}(t)=p(y(t)|X(t)=\hat{X}^{j}(t))$;

  4. 4.

Generate $N$ realizations of the filter density $p(X(t)|Y_{t})$:
$\tilde{X}^{j}(t)$, $j=1,\cdots,N$, by re-sampling (with replacement) of
$\hat{X}^{j}(t)$ with sampling probabilities propositional to $\omega^{j}(t)$.

Thus, the filter density can be constructed by the sampling $\tilde{X}^{j}(t)$
and the prediction density by $\hat{X}^{j}(t)$.Thus, the estimative
expectation of any function $\mathbb{E}(f(X(t))|Y_{t})$ can be calculated as
follows:

$\displaystyle\mathbb{E}(f(X(t))|Y_{t})\approx\frac{1}{N}\sum_{j=1}^{n}f(\tilde{X}^{j}(t)).$

and the predictive expectation of any function $\mathbb{E}(f(X(t))|Y_{t-1})$
can be calculated as follows:

$\displaystyle\mathbb{E}(f(X(t))|Y_{t-1})\approx\frac{1}{N}\sum_{j=1}^{n}f(\hat{X}^{j}(t)).$

##### S3 Appendix.

Prove of the Theorem 1.
Before we proof the theorem, let proof the following lemma:

###### Lemma 1

Let $p_{1},p_{2}$ be the probably distribution function and
$f:{\mathbb{R}}^{l}\rightarrow{\mathbb{R}}$ be such that

$\displaystyle\int|f(y)|^{2}p_{1}(y)dy+\int|f(y)|^{2}p_{2}(y)dy<\infty$

Then

$\displaystyle|\int f(y)(p_{1}(y)-p_{2}(y))dy|\leqslant
2\big{(}\int|f(y)|^{2}p_{1}(y)dy+\int|f(y)|^{2}p_{2}(y)dy\big{)}^{\frac{1}{2}}\
d_{Hell}(p_{1},p_{2})$

Proof

$\displaystyle|\int f(y)(p_{1}(y)-p_{2}(y))dy|$
$\displaystyle\leqslant\int|f(y)||p_{1}(y)-p_{2}(y)|dy$
$\displaystyle=\int\sqrt{2}|f(y)|\sqrt{p_{1}}+\sqrt{p_{2}}|\cdot\frac{1}{\sqrt{2}}|\sqrt{p_{1}}-\sqrt{p_{2}}|dy$
$\displaystyle\leqslant\big{(}2\int|f(y)|^{2}|\sqrt{p_{1}}+\sqrt{p_{2}}|^{2})dy\big{)}^{\frac{1}{2}}({\frac{1}{2}}\int|\sqrt{p_{1}}-\sqrt{p_{2}}|^{2}dy)^{\frac{1}{2}}$
$\displaystyle\leqslant
2\big{(}\int|f(y)|^{2}p_{1}(y)dy+\int|f(y)|^{2}p_{2}(y)dy\big{)}^{\frac{1}{2}}\
d_{Hell}(p_{1},p_{2})$

Back to the prove of Theorem 1.

Proof Let

$\displaystyle A_{n,m}$ $\displaystyle=\frac{\partial}{\partial
h}L(h,Y_{1},\cdots,Y_{n},\Theta_{m})|_{h=\tilde{h}}$
$\displaystyle=-\frac{1}{n}\sum_{i=1}^{n}\frac{\partial}{\partial
h}\log\mathbb{P}(Y_{i}|h,\Theta_{m}^{i})|_{h=\tilde{h}}$

By(R1,R2),

$\displaystyle\int\frac{\partial}{\partial
h}\mathbb{P}(Y|h,\Theta_{m})dY=\frac{\partial}{\partial
h}\int\mathbb{P}(Y|h,\Theta_{m})dY=\frac{\partial}{\partial h}1=0$
$\displaystyle\int\frac{\partial^{2}}{\partial
h^{2}}\mathbb{P}(Y|h,\Theta_{m})dY=\frac{\partial}{\partial
h}\int\frac{\partial}{\partial h}\mathbb{P}(Y|h,\Theta_{m})dY=0$

and it follows that

$\displaystyle\int\frac{\partial}{\partial h}\log\mathbb{P}(Y|h,\Theta_{m})\
\mathbb{P}(Y|h,\Theta_{m})dY=\int\frac{\partial}{\partial
h}\mathbb{P}(Y|h,\Theta_{m})dY=0$

and

$\displaystyle\int\frac{\partial^{2}}{\partial
h^{2}}\log\mathbb{P}(Y|h,\Theta_{m})\ \mathbb{P}(Y|h,\Theta_{m})dY$ (35)
$\displaystyle=\frac{\partial}{\partial h}\int\frac{\partial}{\partial
h}\log\mathbb{P}(Y|h,\Theta_{m})\
\mathbb{P}(Y|h,\Theta_{m})dY-\int\frac{\partial}{\partial
h}\log\mathbb{P}(Y|h,\Theta_{m})\frac{\partial}{\partial
h}\mathbb{P}(Y|h,\Theta_{m})dY$
$\displaystyle=-\int\biggl{(}\frac{\partial}{\partial
h}\log\mathbb{P}(Y|h,\Theta_{m})\biggr{)}^{2}\mathbb{P}(Y|h,\Theta_{m})dY$

By the lemma and R4,

$\displaystyle|$ $\displaystyle\int\frac{\partial}{\partial
h}\log\mathbb{P}(Y|h,\Theta_{m})\ \mathbb{P}(Y|h,\tilde{\Theta}_{m})dY|$ (36)
$\displaystyle=|\int\frac{\partial}{\partial h}\log\mathbb{P}(Y|h,\Theta_{m})\
(\mathbb{P}(Y|h,\tilde{\Theta}_{m})-\mathbb{P}(Y|h,\Theta_{m}))dY|$
$\displaystyle\leqslant 2\big{(}\int|\frac{\partial}{\partial
h}\log\mathbb{P}(Y|h,\Theta_{m})|^{2}(\mathbb{P}(Y|h,\tilde{\Theta}_{m})+\mathbb{P}(Y|h,\Theta_{m}))dY\big{)}^{\frac{1}{2}}d_{Hell}(\mathbb{P}(Y|h,\tilde{\Theta}_{m}),\mathbb{P}(Y|h,\Theta_{m}))$
$\displaystyle\leqslant 4\sqrt{2K}\varepsilon(m),$

and $\lim_{m\rightarrow\infty}4\sqrt{2K}\varepsilon(m)=0$.

Thus $A_{n,m}$ is the mean of $n$ independent random vectors and the expected
value of each random vectors is $0$ as $m$ goes to infinity. Therefore, by the
Weak Law of Large Numbers(WLLN), $\lim_{n,m\rightarrow\infty}A_{n,m}=0$ in
probability.

Let

$\displaystyle B_{n,m}$ $\displaystyle=\frac{\partial^{2}}{\partial
h^{2}}L(h,Y_{1},\cdots,Y_{n},\Theta_{m})|_{h=\tilde{h}}$
$\displaystyle=-\frac{1}{n}\sum_{i=1}^{n}\frac{\partial^{2}}{\partial
h^{2}}\log\mathbb{P}(Y_{i}|h,\Theta_{m}^{i})|_{h=\tilde{h}}$ $\displaystyle
C_{n,m}$ $\displaystyle=H(Y_{1},\cdots,Y_{n},\tilde{h})$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}H(Y,\tilde{h})$

Combined with equation 23, we have

$\displaystyle 0$ $\displaystyle=\frac{\partial}{\partial
h}L(h,Y_{1},\cdots,Y_{n},\Theta_{m})|_{h=\hat{h}}$ (37)
$\displaystyle=\frac{\partial}{\partial
h}L(h,Y_{1},\cdots,Y_{n},\Theta_{m})|_{h=\tilde{h}}+(\hat{h}-\tilde{h})\frac{\partial^{2}}{\partial
h^{2}}L(h,Y_{1},\cdots,Y_{n},\Theta_{m})|_{h=\tilde{h}}$
$\displaystyle+\frac{1}{2}\zeta(\hat{h}-\tilde{h})^{2}H(Y_{1},\cdots,Y_{n},\tilde{h})$
$\displaystyle=A_{n,m}+(\hat{h}-\tilde{h})\biggl{(}B_{n,m}+\frac{1}{2}\zeta(\hat{h}-\tilde{h})C_{n,m}\biggr{)}$

Thus

$\lim_{m,n\rightarrow\infty}(\hat{h}-\tilde{h})\biggl{(}B_{n,m}+\frac{1}{2}\zeta(\hat{h}-\tilde{h})C_{n,m}\biggr{)}=0\quad\text{in
probability}.$ (38)

Analyzing as same as the $A_{n,m}$, we have

$\displaystyle\int\frac{\partial^{2}}{\partial
h^{2}}\log\mathbb{P}(Y|h,\Theta_{m})\
\mathbb{P}(Y|h,\tilde{\Theta}_{m})dY-\int\biggl{(}\frac{\partial}{\partial
h}\log\mathbb{P}(Y|h,\Theta_{m})\biggr{)}^{2}\mathbb{P}(Y|h,\Theta_{m})dY$
(39) $\displaystyle=\int\frac{\partial^{2}}{\partial
h^{2}}\log\mathbb{P}(Y|h,\Theta_{m})\
\mathbb{P}(Y|h,\tilde{\Theta}_{m})dY+\int\frac{\partial^{2}}{\partial
h^{2}}\log\mathbb{P}(Y|h,\Theta_{m})\ \mathbb{P}(Y|h,\Theta_{m})dY$
$\displaystyle\leqslant 4\sqrt{2K}\varepsilon(m).$

According to $\lim_{m\rightarrow\infty}\varepsilon(m)=0$ and assumption R3, we
have

$\displaystyle\lim_{n,m\rightarrow\infty}B_{n,m}=K_{\infty}\quad\text{in
probability},$
$\displaystyle\lim_{n,m\rightarrow\infty}C_{n,m}=\mathbb{E}_{Y}[H(Y)]<\infty\quad\text{in
probability}.$

According to the Slutsky Theorem,

$\displaystyle|\hat{h}-\tilde{h}|=\lim_{n,m\rightarrow\infty}\frac{|A_{n,m}|}{|B_{n,m}+\frac{1}{2}C_{n,m}\zeta(\hat{h}-\tilde{h})|}=0\quad\text{in
distribution}.$ (40)

Furthermore, based on central limit theorem,

$\displaystyle\begin{aligned}
\sqrt{N}|\hat{h}-\tilde{h}|&=\frac{\sqrt{N}|A_{n,m}|}{|B_{n,m}+\frac{1}{2}C_{n,m}\zeta(\hat{h}-\tilde{h})|}\\\
\end{aligned}$ (41)

converges to Gaussian distribution with variance of $1/{K_{\infty}}$.

##### S4 Appendix.

Prove of the Theorem 2.
This proof refers to the book [29].
Proof Let $\rho$, $\rho^{\prime}$ denote the Lebesgue densities on $\mu$,
$\mu^{{}^{\prime}}$, respectively. Then

$\displaystyle\rho(x)$
$\displaystyle=\mathbb{P}(x|y_{1:\mathbb{T}})\propto\mathbb{P}(y_{1:\mathbb{T}},x)\propto\mathbb{P}(y_{1:\mathbb{T}}|x)\mathbb{P}(x)$
$\displaystyle\rho^{\prime}(x)$
$\displaystyle=\mathbb{P}(x|y^{\prime}_{1:\mathbb{T}})\propto\mathbb{P}(y^{\prime}_{1:\mathbb{T}},x)\propto\mathbb{P}(y^{\prime}_{1:\mathbb{T}}|x)\mathbb{P}(x)$

As the equation 31 shows,

$\displaystyle\mathbb{P}(y_{1:\mathbb{T}}|x)$
$\displaystyle=\prod_{t=1}^{\mathbb{T}}\mathbb{P}(y_{t}|x)\propto\prod_{t=1}^{\mathbb{T}}\exp(-I_{filter,t})$
$\displaystyle\propto\prod_{t=1}^{\mathbb{T}}\exp(-\frac{1}{2}|y_{t}-H(x)|_{\Gamma}^{2}-\frac{1}{2}|x-\hat{\mu}_{t}|_{\hat{C}_{t}}^{2})$
$\displaystyle\propto\exp(-\sum_{t=1}^{\mathbb{T}}\frac{1}{2}|y_{t}-H(x)|^{2}_{\Gamma})$
$\displaystyle\mathbb{P}(y^{\prime}_{1:\mathbb{T}}|x)$
$\displaystyle\propto\prod_{t=1}^{\mathbb{T}}\exp(-I_{filter,t})\propto\exp(-\sum_{t=1}^{\mathbb{T}}\frac{1}{2}|y^{\prime}_{t}-H(x)|^{2}_{\Gamma})$

Let $\rho_{0}(x)=\mathbb{P}(x)$ denote the Lebesgue densities on $\mu_{0}$,

$\displaystyle\varrho(x)=\sum_{t=1}^{\mathbb{T}}\frac{1}{2}|y_{t}-H(x)|^{2}_{\Gamma},$
$\displaystyle\varrho^{\prime}(x)=\sum_{t=1}^{\mathbb{T}}\frac{1}{2}|y^{\prime}_{t}-H(x)|^{2}_{\Gamma}.$

Then

$\displaystyle\rho(x)=\frac{1}{Z}\exp(-\sum_{t=1}^{\mathbb{T}}\frac{1}{2}|y_{t}-H(x)|^{2}_{\Gamma})\rho_{0}(x)=\frac{1}{Z}\exp(-\varrho(x))\rho_{0}(x),$
$\displaystyle\rho^{\prime}(x)=\frac{1}{Z^{\prime}}\exp(-\sum_{t=1}^{\mathbb{T}}\frac{1}{2}|y^{\prime}_{t}-H(x)|^{2}_{\Gamma})\rho_{0}(x)=\frac{1}{Z^{\prime}}\exp(-\varrho^{\prime}(x))\rho_{0}(x),$

where

$\displaystyle
Z=\int\exp(-\sum_{t=1}^{\mathbb{T}}\frac{1}{2}|y_{t}-H(x)|^{2}_{\Gamma})\rho_{0}(x)dx,$
$\displaystyle
Z^{\prime}=\int\exp(-\sum_{t=1}^{\mathbb{T}}\frac{1}{2}|y^{\prime}_{t}-H(x)|^{2}_{\Gamma})\rho_{0}(x)dx.$

Thus,

$\displaystyle d_{Hell}(\mu,\mu^{\prime})^{2}=$
$\displaystyle~{}\frac{1}{2}\int|\sqrt{\rho(x)}-\sqrt{\rho^{\prime}(x)}|^{2}dx,$
$\displaystyle=$
$\displaystyle~{}\frac{1}{2}\int|\frac{\sqrt{\exp(-\varrho(x))}}{\sqrt{Z}}-\frac{\sqrt{\exp(-\varrho^{\prime}(x))}}{\sqrt{Z^{\prime}}}|^{2}\rho_{0}(x)dx,$
$\displaystyle\leqslant$
$\displaystyle~{}\frac{1}{2}\int\frac{1}{Z}|\sqrt{\exp(-\varrho(x))}-\sqrt{\exp(-\varrho^{\prime}(x))}|^{2}\rho_{0}(x)dx$
$\displaystyle+\frac{1}{2}|\frac{1}{\sqrt{Z}}-\frac{1}{\sqrt{Z^{\prime}}}|^{2}\int\exp(-\varrho^{\prime}(x))\rho_{0}(x)dx,$
$\displaystyle=$
$\displaystyle~{}\frac{1}{2Z}\int|\exp(-\varrho(x)/2)-\exp(-\varrho^{\prime}(x)/2)|^{2}\rho_{0}(x)dx$
$\displaystyle+\frac{1}{2Z}|\sqrt{Z^{\prime}}-\sqrt{Z}|^{2}.$

And we assume that $y,y^{\prime}$ are both contained in a ball of radius $r$
in the Euclidean norm on $\mathbb{R}^{\mathbb{T}\times q}$. Thus, $\rho_{0}$
is bounded independently of $r$, $\rho$ and $\rho^{\prime}$ are bounded
dependently of $r$. Thus there exists $K_{1}$ depended on $r$ satisfy:

$\displaystyle\frac{1}{2Z}|\sqrt{Z^{\prime}}-\sqrt{Z}|^{2}$
$\displaystyle=\frac{1}{2Z}\frac{|Z-Z^{\prime}|^{2}}{|\sqrt{Z}+\sqrt{Z^{\prime}}|^{2}}$
(42) $\displaystyle\leqslant K_{1}|Z-Z^{\prime}|^{2}$
$\displaystyle=K_{1}|\int\exp(-\varrho(x))\rho_{0}(x)dx-\int\exp(-\varrho^{\prime}(x))\rho_{0}(x)dx|^{2}$
$\displaystyle\leqslant
K_{1}\int|\exp(-\varrho(x))-\exp(-\varrho^{\prime}(x))|^{2}\rho_{0}(x)dx$

Since $\varrho(x)>0$ and $\varrho^{\prime}(x)>0$, using the equation 42 and
the fact that $\exp(-x)$ and $\exp(-x/2)$ are both Lipschitz on
${\mathbb{R}}^{+}$, there exists $K_{2}$ depended on $r$ satisfy

$\displaystyle d_{Hell}(\mu,\mu^{\prime})^{2}\leq$
$\displaystyle~{}\frac{1}{2Z}\int|\exp(-\varrho(x)/2)-\exp(-\varrho^{\prime}(x)/2)|^{2}\rho_{0}(x)dx$
(43)
$\displaystyle+K_{1}\int|\exp(-\varrho(x))-\exp(-\varrho^{\prime}(x))|^{2}\rho_{0}(x)dx$
$\displaystyle\leq$
$\displaystyle~{}K_{2}\int|\varrho(x)-\varrho^{\prime}(x)|^{2}\rho_{0}(x)dx$
$\displaystyle\leq$
$\displaystyle~{}K_{2}\int|\frac{1}{2}\sum_{t=1}^{\mathbb{T}}\left(|y_{t}-H(x_{t})|^{2}_{\Gamma}-|y^{\prime}_{t}-H(x_{t})|^{2}_{\Gamma}\right)|^{2}\rho_{0}(x)dx$
$\displaystyle\leq$
$\displaystyle~{}\frac{K_{2}}{4}\int\biggl{(}\sum_{t=1}^{\mathbb{T}}|y_{t}-y^{\prime}_{t}|_{\Gamma}\times|y_{t}-H(x_{t})+y_{t}-H(x_{t})|_{\Gamma}\biggr{)}^{2}\rho_{0}(x)dx$
$\displaystyle\leq$
$\displaystyle~{}K_{2}\int\biggl{(}\sum_{t=1}^{\mathbb{T}}|y_{t}-y^{\prime}_{t}|^{2}_{\Gamma}\biggr{)}\times\biggl{(}\sum_{t=1}^{\mathbb{T}}|\frac{y_{t}+y^{\prime}_{t}}{2}-H(x_{t})|^{2}_{\Gamma}\biggr{)}\rho_{0}(x)dx$

Finally, we note that since all norms are equivalent on finite-dimensional
spaces, there is constant $K_{3}$ such that

$\sum_{t=1}^{\mathbb{T}}|y_{t}-y^{\prime}_{t}|^{2}_{\Gamma}\leq
K_{3}|y_{1:\mathbb{T}}-y^{\prime}_{1:\mathbb{T}}|^{2}$ (44)

Thus,

$\displaystyle d_{Hell}(\mu,\mu^{\prime})^{2}\leq$ $\displaystyle~{}K_{2}\int
K_{3}|y_{1:\mathbb{T}}-y^{\prime}_{1:\mathbb{T}}|^{2}\times\sum_{t=1}^{\mathbb{T}}|\frac{y_{t}+y^{\prime}_{t}}{2}-H(x_{t})|^{2}_{\Gamma}\rho_{0}(x)dx$
(45) $\displaystyle\leq$ $\displaystyle~{}K_{2}\int
K_{3}|y_{1:\mathbb{T}}-y^{\prime}_{1:\mathbb{T}}|^{2}\times\sum_{t=1}^{\mathbb{T}}\big{(}r^{2}+|H(x_{t})|^{2}\big{)}\rho_{0}(x)dx$
$\displaystyle\leq$
$\displaystyle~{}K_{2}K_{3}r^{2}|y_{1:\mathbb{T}}-y^{\prime}_{1:\mathbb{T}}|^{2}\int\sum_{t=1}^{\mathbb{T}}\big{(}1+|H(x_{t})|^{2}\big{)}\rho_{0}dx$
$\displaystyle=$
$\displaystyle~{}K_{2}K_{3}r^{2}|y_{1:\mathbb{T}}-y^{\prime}_{1:\mathbb{T}}|^{2}\mathbb{E}_{\mu_{0}}w(x_{1:\mathbb{T}})$

Since $\mathbb{E}_{\mu_{0}}w(x_{1:\mathbb{T}})<\infty$, there exists $c=c(r)$
such that for all $|y_{1:\mathbb{T}}|,|y^{\prime}_{1:\mathbb{T}}|<r$,

$\displaystyle d_{Hell}(\mu,\mu^{\prime})\leqslant
c|y_{1:\mathbb{T}}-y^{\prime}_{1:\mathbb{T}}|.$ (46)

## Acknowledgments

This work is jointly supported by the National Key R & D Program of China (No.
2018AAA0100303), the National Natural Science Foundation of China (No.
62072111), the Shanghai Municipal Science and Technology Major Project (No.
2018SHZDZX01) and the ZHANGJIANG LAB.

## References

  * 1. H. Abarbanel.  Predicting the Future.  Predicting the Future, 2013.
  * 2. Bdo Anderson, J. B. Moore, and M. Eslami.  Optimal filtering.  IEEE Transactions on Systems Man & Cybernetics, 12(2):235–236, 2007.
  * 3. A. L. Barabasi and R. Albert.  Albert, r.: Emergence of scaling in random networks. science 286, 509-512.  Science, 286(5439):509–512, 1999.
  * 4. L. Bengtsson, M. Ghil, and E Källén.  Dynamic Meteorology: Data Assimilation Methods.  Dynamic Meteorology: Data Assimilation Methods, 1981.
  * 5. N. Brunel and X. J. Wang.  Effects of neuromodulation in a cortical network model of object working memory dominated by recurrent inhibition.  Journal of Computational Neuroscience, 37(1):63–85, 2001.
  * 6. A. N. Burkitt.  A review of the integrate-and-fire neuron model: I. homogeneous synaptic input.  Biological Cybernetics, 95(1):1, 2006.
  * 7. G. Deco, J. Cabral, V. M. Saenger, M. Boly, and M. L. Kringelbach.  Perturbation of whole-brain dynamics in silico reveals mechanistic differences between brain states.  Neuroimage, 169:46–56, 2018.
  * 8. G. Deco, J. Cabral, M. W. Woolrich, A. Stevner, Tj Van Hartevelt, and M. L. Kringelbach.  Single or multiple frequency generators in on-going brain activity: A mechanistic whole-brain model of empirical meg data.  NeuroImage, 152:538, 2017.
  * 9. G. Deco, M. L. Kringelbach, A. Arnatkeviciute, S. Oldham, and A. Fornito.  Dynamical consequences of regional heterogeneity in the brain’s transcriptional landscape.  Science Advances, 7(29):eabf4752, 2021.
  * 10. A. Doucet, N. D. Freitas, and E. Wan.  The unscented particle filter.  Proc Nips, 13:584–590, 2001.
  * 11. H. Eibern and H. Schmidt.  A four-dimensional variational chemistry data assimilation scheme for eulerian chemistry transport modeling.  JOURNAL OF GEOPHYSICAL RESEARCH, 1999.
  * 12. Elena, Boto, Niall, Holmes, James, Leggett, Gillian, Roberts, Vishal, and Shah.  Moving magnetoencephalography towards real-world applications with a wearable system.  Nature, 2018.
  * 13. J. B. Elsner and T. H. Jagger.  A hierarchical bayesian approach to seasonal hurricane modeling.  Journal of Climate, 17(14):2813–2827, 2004.
  * 14. G. Evensen.  Data assimilation. The ensemble Kalman filter.  false, 2006.
  * 15. L. Fei-Fei and P. Perona.  A bayesian hierarchical model for learning natural scene categories.  In 2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’05), 2005.
  * 16. Friston and Karl.  The free-energy principle: a unified brain theory?  Nature Reviews Neuroscience, 11(2):127, 2010.
  * 17. K. J. Friston, A. Mechelli, R. Turner, and C. J. Price.  Nonlinear responses in fmri: The balloon model, volterra kernels, and other hemodynamics.  Neuroimage, 12(4):466–477, 2000.
  * 18. Afmse Gelfand.  Bayesian statistics without tears: A sampling-resampling perspective.  American Statistician, 1992.
  * 19. N. J. Gordon, D. J. Salmond, and Afm Smith.  Novel approach to nonlinear/non-gaussian bayesian state estimation.  IEE Proceedings F - Radar and Signal Processing, 140(2):107–113, 2002.
  * 20. S. Guo, Z. Yu, D. Fei, X. Hu, and C. Feng.  Hierarchical bayesian inference and learning in spiking neural networks.  IEEE Transactions on Cybernetics, PP(99):1–13, 2017.
  * 21. Rodolfo Guzzi.  Data Assimilation: Mathematical Concepts and Instructive Examples.  09 2015.
  * 22. A. L. Hodgkin and A. F. Huxley.  A quantitative description of membrane current and its application to conduction and excitation in nerve.  The Journal of Physiology, 1952.
  * 23. I. Hoteit, D. T. Pham, G. Triantafyllou, and G. Korres.  A new approximate solution of the optimal nonlinear filter for data assimilation in meteorology and oceanography.  Monthly Weather Review, 136(1):317–334, 2008.
  * 24. B. M. Jobst, R. Hindriks, H. Laufs, E. Tagliazucchi, G. Hahn, A. Ponce-Alvarez, Aba Stevner, M. L. Kringelbach, and G. Deco.  Increased stability and breakdown of brain effective connectivity during slow-wave sleep: Mechanistic insights from whole-brain computational modelling.  Rep, 7(1):4634, 2017.
  * 25. S. J. Julier and J. K. Uhlmann.  A new extension of the kalman filter to nonlinear systems.  Proc.aerosense Int.symp.aerospace/defence Sensing Simulation & Controls, 1997.
  * 26. R.E. Kalman.  A new approach to linear filtering and prediction problems.  Trans of the Asme Journal of Basic Engineering, 1960.
  * 27. Mark Kostuk, Bryan Toth, C Daniel Meliza, Daniel Margoliash, and Henry Abarbanel.  Dynamical estimation of neuron and network properties ii: Path integral monte carlo methods.  Biological cybernetics, 106:155–67, 04 2012.
  * 28. M. L. Kringelbach, A. R. Mcintosh, P. Ritter, V. K. Jirsa, and G. Deco.  The rediscovery of slowness: Exploring the timing of cognition.  Trends in Cognitive Sciences, 19(10):616–628, 2015.
  * 29. Kjh Law, A. M. Stuart, and K. C. Zygalakis.  Data Assimilation: A Mathematical Introduction.  Data Assimilation: A Mathematical Introduction, 2015.
  * 30. T. S. Lee and D. Mumford.  Hierarchical bayesian inference in the visual cortex.  Journal of the Optical Society of America, 20(7):1434–48, 2003\.
  * 31. H. F. Lopes and R. S. Tsay.  Particle filters and bayesian inference in financial econometrics.  Journal of Forecasting, 30(1):168–209, 2011.
  * 32. C. W. Lynn and D. S. Bassett.  The physics of brain network structure, function, and control.  American Physical Society, 2018.
  * 33. Marcus, E., and Raichle.  Colloquium paper: Behind the scenes of functional brain imaging: A historical and physiological perspective.  Proceedings of the National Academy of Sciences of the United States of America, 1998.
  * 34. Warren S. Mcculloch and Walter Pitts.  A logical calculus of the ideas immanent in nervous activity.  biol math biophys, 1943.
  * 35. P. Mullowney and S. Iyengar.  Parameter estimation for a leaky integrate-and-fire neuronal model from isi data.  Journal of Computational Neuroscience, 24(2):179–194, 2008.
  * 36. A. E. Pereda.  Electrical synapses and their functional interactions with chemical synapses.  Nature Reviews Neuroscience.
  * 37. N. Politi, J. Feng, and W. Lu.  Comparing data assimilation filters for parameter estimation in a neuron model.  In International Joint Conference on Neural Networks, 2016.
  * 38. J. C. Quinn and Hdi Abarbanel.  State and parameter estimation using monte carlo evaluation of path integrals.  Quarterly Journal of the Royal Meteorological Society, 136(652), 2010.
  * 39. Adeel Razi, Mohamed Seghier, Yuan Zhou, Peter McColgan, Peter Zeidman, Hae-Jeong Park, Olaf Sporns, Geraint Rees, and Karl Friston.  Large-scale dcms for resting state fmri.  Network Neuroscience, 1:1–41, 05 2017.
  * 40. Stam and J. Cornelis.  Modern network science of neurological disorders.  Nature Reviews Neuroscience, 15(10):683, 2014.
  * 41. Stein and B Richard.  A theoretical analysis of neuronal variability.  Biophysical Journal, 5(2):173–194, 1965.
  * 42. Bryan Toth, Mark Kostuk, C Daniel Meliza, Daniel Margoliash, and Henry Abarbanel.  Dynamical estimation of neuron and network properties i: Variational methods.  Biological cybernetics, 105:217–37, 10 2011.
  * 43. C Vanvreeswijk and H. Sompolinsky.  Chaos in neuronal networks with balanced excitatory and inhibitory activity.  Science, 274(5293):p.1724–1726, 1996.
  * 44. D. V. Vavoulis, V. A. Straub, Jad Aston, J. Feng, and L. J. Graham.  A self-organizing state-space-model approach for parameter estimation in hodgkin-huxley-type models of single neurons.  Plos Computational Biology, 8(3):e1002401, 2012.
  * 45. D. J. Watts and S. H. Strogatz.  Collective dynamics of ’small-world’ networks.  Nature, 1998.
  * 46. Christopher K. Wikle, L. Mark Berliner, and Ralph F. Milliff.  Hierarchical bayesian approach to boundary value problems with stochastic boundary conditions.  Monthly Weather Review, 131(6):1051–1062, 2003.
  * 47. H. R. Wilson and J. D. Cowan.  Excitatory and inhibitory interactions in localized populations of model neurons.  Biophysical Journal, 12(1):1–24, 1972.
  * 48. G. R. Yang, J. D. Murray, and X. J. Wang.  A dendritic disinhibitory circuit mechanism for pathway-specific gating.  Nature Communications, 7(1):12815, 2016.