Generative framework for dimensionality reduction of large scale network of non-linear dynamical systems driven by external input

TL;DR

This paper introduces a generative framework for reducing the dimensionality of large-scale nonlinear dynamical networks driven by external inputs, enabling more accurate and analytically tractable representations.

Contribution

It proposes a novel mode selection method based on input distribution symmetry, relaxing heuristic approaches in network reduction.

Findings

Order parameter is conserved in original and reduced models

Framework is applicable to networks of Van der Pol oscillators and Hindmarsh-Rose neurons

Reduction technique is analytically tractable and validated on real-world scenarios

Abstract

Several studies have proposed constraints under which a low dimensional representation can be derived from large-scale real-world networks exhibiting complex nonlinear dynamics. Typically, these representations are formulated under certain assumptions, such as when solutions converge to attractor states using linear stability analysis or using projections of large-scale dynamical data into a set of lower dimensional modes that are selected heuristically. Here, we propose a generative framework for selection of lower dimensional modes onto which the entire network dynamics can be projected based on the symmetry of the input distribution for a large-scale network driven by external inputs, thus relaxing the heuristic selection of modes made in the earlier reduction approaches. The proposed mode reduction technique is tractable analytically and applied to different kinds of real-world large-scale network scenarios with nodes comprising of a) Van der Pol oscillators b) Hindmarsh-Rose neurons. These two demonstrations elucidate how order parameter is conserved at original and reduced descriptions thus validating our proposition.

Equations

$${\tau }_n{\dot{\mathbf{x}}}_i^{(n)}=$$

$${\tau }_n{\dot{\mathbf{x}}}^{(n)}(t,\mathbf{I})=$$

$$\mathbf{x}(t,\mathbf{I})=\sum _{i=1}^M{\mathbf{\xi }}_i(t)v_i(\mathbf{I})+R(t,\mathbf{I})$$

$${\tau }_n{\dot{\mathbf{\xi }}}_i^{(n)}=$$

$${\mathbf{\xi }}_i=$$

$$\mathbf{F}\left({\mathbf{\xi }}_i\right)=$$

$$\mathbf{G}\left(\{{\mathbf{\xi }}_j\},K,i\right)=$$

$$\mathbf{H}\left(\{{\mathbf{\xi }}_j^{(m)}\},W,n,i\right)=$$

$${\mathbf{II}}_i^{(n)}=$$

$$\dot{x}=$$

$$\dot{y}=$$

$$\mathbf{f}\left({\mathbf{x}}_i\right)=$$

$$\mathbf{g}\left(\{{\mathbf{x}}_j\},K,i\right)=$$

$$\mathbf{h}\left(\{{\mathbf{x}}_j^{(m)}\},W,n\right)=$$

$$\mathbf{I}=$$

$$\mathbf{F}\left({\mathbf{\alpha }}_i\right)=$$

$$\mathbf{G}\left(\{{\mathbf{\alpha }}_j\},K,i\right)=$$

$$\mathbf{H}\left(\{{\mathbf{\alpha }}_j^{(m)}\},W,n,i\right)=$$

$$\mathbf{II}=$$

$$\dot{x}=$$

$$\dot{y}=$$

$$\dot{z}=$$

$$\mathbf{f}\left({\mathbf{x}}_i\right)=$$

$$\mathbf{g}\left(\{{\mathbf{x}}_j\},K,i\right)=$$

$$\mathbf{h}\left(\{{\mathbf{x}}_j^{(m)}\},W,n\right)=$$

$$\mathbf{I}=$$

$$\mathbf{F}\left({\mathbf{\xi }}_i\right)=$$

$$\mathbf{G}\left(\{{\mathbf{\xi }}_j\},K,i\right)=$$

$$\mathbf{H}\left(\{{\mathbf{\xi }}_j^{(m)}\},W,n,i\right)=$$

$$\mathbf{II}=$$

$$W=\left[\begin{array}{ccc}0& 1& 1\\ -1& 0& 1\\ -1& -1& 0\end{array}\right]\times w$$

$$\displaystyle error(M)=\sqrt{\sum_{p=<f,w,K,\mu>}^{T}\left(\begin{array}[]{c}C_{p}-R_{p}\end{array}\right)^{2}}$$

$$A_{ij}$$

$$a_i$$

$$I{I}_i$$

$$A_{ij}$$

$$B_{ij}$$

$$C_{ij}$$

$$a{1}_i$$

$$a{2}_i$$

$$b{1}_i$$

$$b{2}_i$$

$$c{1}_i$$

$$c{2}_i$$

$$d{1}_i$$

Full text

††preprint: APS/123-QED

Large-scale dynamical systems are useful tools to explain a wide variety of complex phenomena in nature e.g. financial markets HSIEH , jamming transitions Charbonneau et al. (2017), human mobility dynamics Song et al. (2010), weather patterns Antoulas et al. (2010) and brain dynamics Deco et al. (2008). While increase in scale or dimensions may increase the predictive power of the model system, nonetheless a reduction to simpler descriptions at lower dimensions is critical for having relevance to empirical observations and analytical tractability of underlying mechanisms governing empirical observations. One robust approach of reducing dimensions is defining modes on which the original system can be projected Antoulas et al. (2010); Stefanescu and Jirsa (2008). The selection of a mode is often heuristically motivated, and the mode can also be an order parmeter from the perspective of the paradigmatic framework of Synergetics Haken (1983). In Neuroscience, reduction of dynamical systems with respect to modes constructed from distribution of external input has been performed earlier on small-scale network of linearly coupled excitable systems Stefanescu and Jirsa (2008). Since this reduction retain important network dynamics, large-scale networks were conceptualized by coupling these reduced systems with long-range coupling Becker et al. (2015); Sanz-Leon et al. (2015), the later being heuristically argued from symmetry properties. In present work, we perform reduction on a large-scale network where connection among nodes involve global and local coupling mimicking a real-world system. Subsequently, long-range coupling term between modes in the reduced system is derived analytically as part of the reduction process. Global coherence is an order parameter that can be computed both at the level of original dynamical system as well as from the mode dynamics in the reduced system. Conservation of global coherence at both levels is used to validate the generality of our approach in two distinct networks. First, we simulate a large-scale network where each node is a Van der Pol oscilator having 2-dimensional dynamics and coupled using local and global parameters. Second each node is a Hindmarsh-Rose neuron, a 3-dimensional dynamical system which can exhibit different time scales of oscillations resulting in bursting along with tonic spiking behavior.

Each $i^{th}$ excitable system in node $n$ is represented by a vector of state variables (state vector) $\bm{x}^{(n)}_{i}$. Input to $i^{th}$ excitable system from interaction with other excitable systems within $n^{th}$ node is given by a vector function $\bm{g}\left(\{\bm{x}^{(n)}_{j}\},K,i\right)$ where $K$ is coupling constant and $\{\bm{x}^{(n)}_{j}\}$ is a set of $J$ state-vectors ($j=1,...,J$) in $n^{th}$ node. Input to every excitable system in $n^{th}$ node from other nodes is given by another vector function $\bm{h}\left(\{\bm{x}_{j}^{(m)}\},W,n\right)$ where elements of matrix $W$ ($W\in\Re^{J\times N}$) are the weights of connections between nodes and $\{\bm{x}_{j}^{(m)}\}$ is a set of state-vectors in $N$ nodes ($m=1,...,N$) (Figure 1A). Vector function $\bm{f}(\bm{x}_{i}^{(n)})$ contributes to the local dynamics of $i^{th}$ excitable system in $n^{th}$ node. Then, the dynamics of the entire system is described by the following set of equations

[TABLE]

where $\dot{\bm{x}}$ is time derivative of state vectors, $\bm{\phi}$ is a function of the external input ($\bm{I}$) and $\tau_{n}$ is the time-constant of $n^{th}$ node which is a differentiating factor between nodes. However, within $n^{th}$ node, external current to $i^{th}$ excitable system ($\bm{I}^{(n)}_{i}$) differentiates it from the rest. For a large scale network comprising of individual excitable nodes, (1) can be written as

[TABLE]

where $I$ is a continuous variable having normal distribution $\mathbb{N}(\mu,\sigma)$.

Now, we can represent $\bm{x}(t,I)$ as a superposition of $M$ bi-orthogonal modes $\{v_{i}\}$

[TABLE]

where $R(t,\bm{I})$ is the residual and $M\lll J$ (Figure 1B). The nature of reduction is such that dynamical system given in (1) is reduced to solving for the mode coefficients $\bm{\xi}_{i}$ as described in the following set of equations

[TABLE]

where

[TABLE]

and $\{v_{i}^{+}\}$ are the adjoint basis for the biorthogonal modes $\{v_{i}\}$.

Large-scale network of Van der Pol oscillators:

A Van der Pol oscillator van der Pol (1920) has two state variables $x$ and $y$ which follows the following equations

[TABLE]

A large-scale network where individual node is essentially a Van der Pol oscillator can be cast into equation (1) with the following relations

[TABLE]

where $a$ is a constant and $\mathbb{E}[\{x_{j}\}]$ is the expected value of the state-variables $\{x_{j}\}$. For the reduced system described in (4) with $\xi_{i}=[\alpha_{i},\beta_{i}]^{T}$, we derive

[TABLE]

where the constants $a_{i}$ and $A_{ij}$ are computed by applying bi-orthogonal assumption of modes as stated in Appendix A.

Large-scale network of Hindmarsh-Rose neurons:

In Neuroscience, Hindmarsh-Rose (HMR) neuron is a three dimensional model of single neuron firing dynamics having three state variables Hindmarsh and Rose (1984)

[TABLE]

where $a$, $b$, $c$, $d$, $r$, $s$ and $x0$ are the constants. A network of excitatory and inhibitory HMR neurons have been used to describe a small-scale network of neurons Stefanescu and Jirsa (2008). Thus, a node in the brain (Fig 1a) can be expressed as a six dimensional state space with 3 excitatory and 3 inhibitory variables represented as $\bm{x}=\left[x_{1},y_{1},z_{1},x_{2},y_{2},z_{2}\right]^{T}$. Extending this architecture, for the large-scale system in (1), we obtain

[TABLE]

where $K_{11}$, $K_{12}$ and $K_{21}$ are coupling constants between excitatory and inhibitory nodes. For the reduced system in (4), $\xi_{i}=[\alpha_{1},\beta_{1},\gamma_{1},\alpha_{2},\beta_{2},\gamma_{2}]^{T}$, we derive

[TABLE]

where the constants ${a_{1}}_{i}$, ${a_{2}}_{i}$, ${b_{1}}_{i}$, ${b_{2}}_{i}$, ${c_{1}}_{i}$, ${c_{2}}_{i}$, ${d_{1}}_{i}$, ${d_{2}}_{i}$, ${p_{1}}_{i}$, ${p_{2}}_{i}$, $A_{ij}$, $B_{ij}$ and $C_{ij}$ are defined in Appendix.

We simulated a network of three nodes with each functional unit in a node is governed by 1) Van der Pol (VDP) oscillator or 2) HMR neurons. We select the node connection matrix of the given form

[TABLE]

where $w$ is a scalar in the range [math] to $1$ representing negligible node interaction to strong node interaction. We also considered different values of $K$ for VDP (and Inhibition-excitation ratio, $IER=\frac{K_{12}}{K_{11}}$ for HMR) and for different values of mean current $\mu$. For validating reduced system with $M$ modes, we compare the Global Coherence (GC) evaluated using the complete system ($C_{p}$) with GC evaluated using the reduced system ($R_{p}$) using the following equation.

[TABLE]

where $f$ is frequency and $T$ is total number of pairs $<f,w,K,\mu>$. GC is computed between mean activities of each node. For the reduced system, the mean activity of each node is the mean of individual excitable systems’ activities estimated using (3) without the residual.

As $M$ increases, error in proposed reduction process decreases more rapidly as compared to previous approaches (Figure 2A & B) for both large-scale networks using Van der-Pol and Hindmarsh-Rose models as nodes, thus, validating our reduction approach. For Van der-Pol model, reduction proposed in Sanz-Leon et al. (2015) generated numerical instability during simulation. Global coherence calculated from proposed reduction for an exemplar parameter space matches closely with the original system unlike heuristic approaches (Figure 2C). We further validate our framework by showing the reproduction of time-series of mean-field activity of each HMR and VDP node (Figure 3) for several different parameters spaces.

In Stefanescu and Jirsa (2008), local-interaction between state-variables was facilitated via the mean field activity of the node which by itself is a small-scale network. In the reduced model this interaction was preserved as the input to state variable of mode $i$ was characterized via mode-specific output of the node ($\sum_{j=1}^{M}A_{ij}\alpha_{j}$) governed by matrix $A$ which is obtained as a part of reduction. However, in the case of the large-scale network of these reduced nodes the input to state-variable of mode $i$ of node $m$ was either the activity of state-variable of mode $i$ from other nodes Becker et al. (2015) or it was the sum of activities of state-variables of all modes from other nodes Sanz-Leon et al. (2015). In this paper, the input to mode $i$ of node $m$ is derived by projecting the long-range interaction term of complete large-scale network to modes of the external input. Thus, our proposed scheme preserves the original long range interaction (Figure 2 and Figure 3).

In summary, we propose a generalized scheme for reduction of the dynamics of a large-scale network into lower dimensional mode description based on properties of the external input. Obviously, any such reduction lowers the computational complexity. However an important point to note is a model’s benefit is not necessarily limited to mimicking the complex dynamics of real-world system. For example, a detailed model of the cortical layer will be highly informative Markram et al. (2015), but not necessarily insightful for explaining the cortical interactions during a specific behavioral task. How do we develop task-specific insights from a high dimensional data with minimal assumptions will be critical question for future studies and how they could profit from our proposed formalism. Our approach is best-placed to application where task inevitably involves reduction from high-dimensional state space to functionally well defined lower dimensional modes. Such problems are also not necessarily exclusively brain specific, but also pertinent for climate dynamics, traffic problems, or anywhere dynamical systems are driven by external time varying and distributed input and the approach presented in this article may be insightful for a wide range of scientific disciplines.

This research was funded by NBRC core and the grants Ramalingaswami fellowship (BT/RLF/Re-entry/31/2011) and Innovative Young Bio-technologist Award (IYBA) (BT/07/IYBA/2013) from the Department of Biotechnology (DBT), Ministry of Science & Technology, Government of India to AB. DR is supported by the Ramalingaswami fellowship (BT/RLF/Re-entry/07/2014) and DST extramural grant (SR/CSRI/21/2016). We also thank Prof. Michael Breakspear for helpful comments and discussions.

Appendix A: Reduction coefficients for Van der Pol oscillator network

[TABLE]

$g(I)$ is the pdf of external input.

Appendix B: Reduction coefficients for Hindmarsh-Rose neuronal network

[TABLE]

$g1(I)$ and $g2(I)$ are the pdfs of external input to excitatory and inhibitory sub-populations respectively.