Model Order Reduction in Neuroscience

Due to the advances in recording and imaging technologies, the number of recorded signals from brain cells increased significantly in the last few years. The recorded spatio-temporal neural activity give rise to networks with complex dynamics. In neuroscience, molecular and cellular level details are incorporated in large-scale models of the brain in order to reproduce phenomena such as learning and behavior. The rapid growth of simultaneous neuronal recordings in scale and resolution brings challenges with the analysis of the neuronal population activity. New computational approaches have to be developed to analyze, visualize, and understand large-scale recordings of neural activity. While algorithmic developments and the availability of significantly more computing power have enabled analysis of larger neuronal networks, these advances cannot keep pace with increasing size and complexity of recorded activity. The activity of complex networks of neurons can often be described by relatively few distinct patterns. Model order reduction techniques enable us to identify the coherent spatial–temporal patterns.

The presence or absence of a neural mechanism can be analyzed for neuronal populations. Dimensionality reduction methods [6] are data-driven statistical techniques for forming and evaluating hypotheses about population activity structure, which are summarized in Section 2. One of the goals of neurosciences is fast and accurate predictions of the potential propagation in neurons. The differential equations describing the propagation of potential in neurons were developed and are described by Hodgkin and Huxley equations [12]. They consists of a coupled system of ordinary and partial differential equations (ODEs and PDEs). The dimension of the associated discretized systems is very large for accurately simulating neurons with realistic morphological structure and synaptic inputs. In Section 3 we present two model order reduction approaches based on POD and DEIM [5] which can predict accurately the potential propagation in large scale neuronal networks leading to important speedups [17, 16, 2]. Using the functional neuroimagining data from electroencephalography (EEG) or functional magnetic resonance imaging (fMRI), different regions of the brain can be inferred by dynamic causal modeling (DCM) [7]. Effective connectivity is parameterised in terms of coupling among unobserved brain states and neuronal activity in different regions. In Section 4 we describe a combined state and parameter reduction for parameter estimation and identification [10] to extract effective connectivity in neuronal networks from measured data, such as data from electroencephalography (EEG) or functional magnetic resonance imaging (fMRI). In Section 5 the data-driven, equation free, model order reduction method dynamic mode decomposition (DMD) is described for identifying sleep spindle networks [3]. Reduced order models with POD and DEIM and four variants of them are presented for neuronal synaptic plasticity and neuronal spiking networks in Section 6.

Coordination of responses across neurons exists only at the level of the population and not at the level of single neurons. The presence or absence of a neural mechanism can be analyzed for neuronal populations. Dimensionality reduction methods are data-driven statistical techniques for forming and evaluating hypotheses about population activity structure. They produce low-dimensional representations of high-dimensional data with the aim to extract coherent patterns which preserve or highlight some feature of interest in the data [6]. The recorded neurons of dimension $D$ are likely not independent of each other, because they belong to a common network of neuronal populations. From the high-dimensional data of neuronal recordings, a smaller number of explanatory variables $K$ ( $K<D$) are extracted with the help of dimensionality reduction methods. The explanatory variables are not directly observed, therefore they are referred to as latent variables. The latent variables define a $K$-dimensional space representing coherent patterns of the noisy neural activity of $D$ neurons.

There exists several dimensionality reduction methods which differ in the statistical interpretation of the preserved and discarded features of the neuronal populations. We summarize the commonly used statistical methods of dimensionality reduction methods in [6], where further references about the methods can be found.

Principal component and factor analysis; The most widely known technique to extract coherent patterns from high dimensional data is the modal decomposition. A particularly popular modal decomposition technique is principal component analysis (PCA), which derives modes ordered by their ability to account for energy or variance in the data. In particular, PCA is a static technique and does not model temporal dynamics of time-series data explicitly, so it often performs poorly in reproducing dynamic data, such as recordings of neural activity. The low-dimensional space identified by PCA captures variance of all types, including firing rate variability and spiking variability, whereas factor analysis (FA) discards the independent variance for each neuron. and preserves variance that is shared across neurons.

Time series methods: The temporal dynamics of the population activity can be identified if the data comes from a time series. The most commonly used time series methods for dimensionality reduction neural recordings are: hidden Markov models (HMM) [18], kernel smoothing followed by a static dimensionality reduction method, Gaussian process factor analysis (GPFA) [35], latent linear dynamical systems (LDS) [4] and latent nonlinear dynamical systems (NLDS) [26]. They produce latent neural trajectories that capture the shared variability across neurons. The HMM is applied when a jump between discrete states of neurons exists, other methods identify smooth changes in firing rates over time.

Methods with dependent variables: On many neuronal recordings the high-dimensional firing rate space is associated with labels of one or more dependent variables, like stimulus identity, decision identity or a time index. The dimensionality reduction aims in this case to project the data such that differences in these dependent variables are preserved. The linear discriminant analysis (LDA) can be used to find a low-dimensional projection in which the $G$ groups to which the data points belong are well separated.

Nonlinear dimensionality reduction methods: All the previous methods assume a linear relationship between the latent and observed variables. When the data lies on a low-dimensional, nonlinear manifold in the high-dimensional space, a linear method may require more latent variables than the number of true dimensions of the data. The most frequently used methods to identify nonlinear manifolds are Isomap79 [31] and locally linear embedding (LLE) [28]. Because the nonlinear methods use local neighborhoods to estimate the structure of the manifold, population responses may not evenly explore the high-dimensional space. Therefore theses methods should be used with care.

One of the goals of neurosciences is fast and accurate predictions of the potential propagation in neurons. The differential equations describing propagation of potential in neurons are described by Hodgkin and Huxley (HH) cable equations [12]. They consists of a coupled system of ordinary (ODEs) and partial differential equations PDEs. Accurate simulation of morphology, kinetics and synaptic inputs of neurons requires solution of large systems of nonlinear ODEs. The complexity of the models are determined by the synapse density of the dentritic length ($1\mu$ one micron). In simulations, for one synapse per micron on a cell $5$ mm dendrite at $5,000$ compartments each with $10$ variables are needed, which results in $50,000$ coupled nonlinear system of ODEs [17, 16]. To recover complex dynamics, efficient reduced order neuronal methods are developed using the POD and DEIM from the snapshots of the in space and time discretized coupled PDEs and ODEs [17, 16, 2]. In this section we describe two of them. They differ in the formulation of the HH cable equation and of the equations for the gating variables.

In neuroscience different regions of the brain are inferred using neuroimagining data from EEG or fMRI recordings using the method od dynamic causal modeling (DCM) [7]. Effective connectivity is parameterised in terms of coupling among unobserved brain states and neuronal activity in different regions. In DCM the neuronal activity is of the observed brain region is represented as a SISO (single input single output) linear state-space system

with the parameterized connectivity $A_{\mathrm{d}yn}(\mu)$ and external input matrices $B_{\mathrm{d}yn}$.

Linearization of the nonlinear DCM hemodynamic forward sub-model (balloon model) [7] transforms the neuronal activity to the measured BOLD (blood oxygen level dependent) response. Linearization around the equilibrium results in the following single input, single output (SISO) system:

The fMRI measurements at the $i^{th}$ brain region are reflected by the output variables $y_{i}$. For the meaning of the variables and parameters in (15) and (19) we refer to [10, 9]. The linearized forward sub-models are embedded into the fMRI connectivity model

where $\delta_{ij}\in\mathbb{R}^{4\times N_{\mathrm{d}yn}}$ denotes the Kronecker matrix with non-zero elements located at the $(i,j)^{th}$ component.

The linearized state-space forward model (20) and (21) corresponds to a multiple input, multiple output (MIMO) system

where $x\in\mathbb{R}^{N}$ is the internal state, $u\in\mathbb{R}^{J}$ the external input, $y\in\mathbb{R}^{O}$ the observed output and $\mu$ are the parameters describing different conditions.

For large number of $M:=N^{2}$ parameters, the computational cost for inferring the parameters and states is very high. In [10, 8] a combined state and parameter model order reduction is developed for parameter estimation and identification to extract effective connectivity. The inversion procedure consists of two phases, the offline and online phases. In the offline phase, the underlying parameterized model is reduced jointly in states and parameters. In online phase, the reduced order model’s parameters are estimated to fit the observed experimental data. Using state and parameter reduction, the computational cost is reduced in the offline phase. The simultaneous reduction of state and parameter space is based on Galerkin projections with the orthogonal matrices for the state $V\in\mathbb{R}^{N\times n}$ and for the parameters $P\in\mathbb{R}^{M\times m}$. The reduced model is of lower order $n<<N,\;m<<M$ than the original full order model. The reduced states $x_{r}(t)\in\mathbb{R}^{n}$ and the reduced parameters $\mu\in\mathbb{R}^{m}$ are computed as

with a reduced initial condition $x_{r,0}=V^{T}x_{0}$ and the reduced components

In the online phase, the optimization based inverse problem is combined with the reduction of state and parameter space. The inversion is based on generalized data-driven optimization approach to construct the ROMs in [23] and enhanced with the Monte-Carlo method to speed up the computations. The state projection $V\in\mathbb{R}^{N\times n}$ and parameter projection $P\in\mathbb{R}^{m\times m}$ are determined iteratively based on a greedy algorithm by maximizing the error between the high-fidelity original and the low-dimensional reduced model in the Bayesian setting.

Numerical experiments with the DCM model in [23] show highly dimensional neuronal network system is efficiently inverted due to the short offline durations. In the offline phase, Monte-Carlo enhanced methods require more than an order of magnitude less offline time compared to the original and data-misfit enhanced methods. In the online phase the reduced order model has a speedup factor about an order of magnitude more compared to the full-order inversion. The output error of the data-misfit enhanced method is close to full order method. The output-errors of Monte-Carlo decrease with increasing number of simulation but does not reach the output error of the full order system. The source code is available in MATLAB [8].

Dynamic mode decomposition (DMD) is a data-driven, equation free ROM technique [20]. It was initially developed to reduce the high dimensional dynamic data obtained from experiments and simulations in the fluid mechanics into a small number of coupled spatial–temporal modes [29, 30]. DMD was applied to explore spatial–temporal patterns in large-scale neuronal recordings in [3]. DMD can be interpreted as combination of discrete Fourier transform (DFT) in time and principal component analysis (PCA) [15] in space. Both PCA and independent component analyses (ICA) [13] are static techniques, which perform poorly in reproducing dynamic data, such as recordings of neural activity.

The data is taken from electrocorticography (ECoG) recordings. Voltages from $n$ channels of an electrode array sampled every $\Delta t$. These measurements are arranged at snapshot $k$ to the column vectors ${\mathbf{x}}_{k}$. The $m$ snapshots in time construct to data matrices shifted in time with $\Delta t$

These matrices are assumed to be related linearly in time

The DMD of the data matrix pair ${\mathbf{X}}$ and ${\mathbf{X}}^{\prime}$ is given by the eigendecomposition of ${\mathbf{A}}$ using the singular value decomposition (SVD) of the data matrix ${\mathbf{X}}=U\Sigma V^{*}$ by computing the pseudoinverse ${\mathbf{A}}\approx{\mathbf{X}}^{\prime}{\mathbf{X}}^{\dagger}\equiv{\mathbf{X}}^{\prime}{\mathbf{V}}{\mathbf{\Sigma}}^{-1}{\mathbf{U}}^{*}.$ The spatio-temporal modes are computed by the exact DMD algorithm [32].

Because DMD does not contain explicit spatial relationship between neighboring measurements, traveling waves occurring in neuronal networks can not be captured well with a few coherent modes. DMD was also used as a windowed technique with a temporal window size constrained by lower bound as for the discrete Fourier transformation (DFT). In contrast to the fluid dynamics where $n>>m$, in neuroscience the electrode arrays that have tens of channels $n$ in the recordings with $m$ number of snapshots in the windows data per second, so that $n<m$. The number of singular values $v$ of ${\mathbf{X}}$ are less than $n$ and $m-1$, which restricts the maximum number of DMD modes and eigenvalues to $n$. Because of this the dynamics can be captured over $m$ snapshots. The rank mismatch is resolved by appending to the snapshot measurements with $h-1$ time-shifted versions of the data matrices. The augmented data matrix ${\mathbf{X}}_{\mathrm{a}ug}$ is given as

The augmented matrices ${\mathbf{X}}_{{\mathrm{a}ug}}$ and ${\mathbf{X}}^{\prime}_{{\mathrm{a}ug}}$ are Hankel matrices, which are constant along the skew diagonal, as in the Eigenvalue Realization Algorithm (ERA) [14]. The number of the stacks $h$ is chosen such that $hn>2m$. A measure to determined the optimal number of stacks $h$ is the approximation error

where $||\cdot||_{F}$ is the Frobenius norm. The approximation error $E$ is decreasing with increasing number of stacks $h$ and reaches a plateau, so that the DMD accuracy does not significantly increases.
DMD is applied in [3] as an automated approach to detect and analyze reliably spatial localization and frequencies of sleep spindle networks from human ECoG recordings. A MATLAB implementation is available at github.com/bwbrunton/dmd-neuro/.

Recently reduced order models for ODEs

are constructed using POD and DEIM to investigate input-output behavior of the neuronal networks in the brain [22, 21], where $x(t)$ are state, and $u(t)$ are input variables, respectively.

The model in [22] is based on the chemical reactions of molecules in synapses, that are the intercellular information transfer points of neurons. The signaling pathways in striatal synaptic plasticity are modeled in [19]. This model describes how certain molecules, which are a prerequisite for learning in the brain, act in synapses. The stoichiometric equations obey the law of mass action, which leads to a deterministic system of $44$ ODEs of the form (27) . The state $x(t)$ of the control system (27) is a collection of ions, molecules, and proteins that act in neuronal synapses. The linear part of (27) is sparse, the nonlinearities are quadratic. The time dependent stimulus $u(t)$ consists of molecules that are important for neuronal excitability and plasticity, calcium and glutamate.

In [21], a nonlinear biophysical network model is considered, describing synchronized population bursting behavior of heterogeneous pyramidal neurons in the brain [27]. Neurons communicate by changing their membrane voltage to create action potentials (spikes), propagating from cell to cell. Spiking is the fundamental method of sensory information processing in the brain, and synchronized spiking is an emergent property of biological neuronal networks. The ODE system (27) in [21] consists of the states $x(t)$ as a collection of $50$ neurons, each modeled with $10$ ODEs, leading totally to a system of ODEs of dimension $500$. Each cell is modeled with Hodgkin-Huxley equations, where each cell has only two compartments (dendrites and soma) and these compartments have different ion channels. The state variables $x(t)$ include the voltages of somatic and dendritic compartments, dendritic calcium concentration, synaptic and ion channel gating variables and the input $u(t)$ is an injected current. Additionally, the soma compartment voltages are coupled to dentritic compartments of randomly chosen cells. This networking of the output of cells as input to other cells is key for producing synchronized population behavior. The nonlinearities are Hodgkin-Huxley type, i.e. exponential functions as well as cubic and quartic polynomials.

In [22], POD+DEIM was applied to a data-driven biological model of plasticity in the brain (27). The ROMs with POD-DEIM reduce significantly the simulation time and error between the original model and reduced order solutions can be tuned by adjusting the number POD and DEIM bases independently. When the ROMs are trained in a matching time interval of $10000$ seconds, accurate results are obtained. However, generalizing the reduced model to longer time intervals is challenging, which is characteristic for all nonlinear models.

In [21], the network model (27) is reduced with localized DEIM (LDEIM) [24], discrete adaptive POD (DAPOD) [33, 34], and adaptive DEIM [25]. DEIM and the variations are used here in combination with POD. ROMs require large number POD and DEIM bases, to accurately simulate the input-output behavior in the ROMs. In this model, every cell is heterogeneous in parameters and there are also jump/reset conditions, which are factors that pose additional challenges to the reduced order methods. However, the ROMs in [21] were able to replicate the emergent synchronized population activity in the original network model. DAPOD and ADEIM perform best in preserving the spiking activity of the original network model. ADEIM is too slow and does not allow low enough dimensions to offset the computational costs of online adaptivity. DAPOD is able to find a lower dimensional POD basis online than the other methods find offline, but the runtime is close to the original model.