Dynamical and Coupling Structure of Pulse-Coupled Networks in Maximum Entropy Analysis

Abstract

Maximum entropy principle (MEP) analysis with few non-zero effective interactions successfully characterizes the distribution of dynamical states of pulse-coupled networks in many fields, e.g., in neuroscience. To better understand the underlying mechanism, we found a relation between the dynamical structure, i.e., effective interactions in MEP analysis, and the anatomical coupling structure of pulse-coupled networks and it helps to understand how a sparse coupling structure could lead to a sparse coding by effective interactions. This relation quantitatively displays how the dynamical structure is closely related to the anatomical coupling structure.

1. Introduction

Binary-state networks—each node in one sampling time bin is binary-state—arise from many research fields, e.g., gene regulatory modeling and neural dynamics [1,2,3]. Statistical distributions of network states are essential to encode information [4,5,6,7,8]. For example, with statistical distributions of network states, experimental studies show that rats can perform awake replays of remote experiences in the hippocampus [9]. The number of network states of n binary-state nodes, $2^n$, exponentially grows as the node number, n, increases, which creates a challenge of characterizing the probability distribution of network states. Many works effectively characterize the distribution of network states in various systems, e.g., a network of ∼100 neurons [10], with a low-order maximum entropy principle (MEP) analysis [10,11,12,13,14,15,16,17,18]—a method with few (far less than $2^n$) non-zero effective interactions (see a precise definition in Equation (1)) constrained by low-order statistics. In the MEP analysis, the dynamical structure of the network, i.e., how nodes interact with each other in the recording of dynamical data, is characterized by effective interactions. This dynamical structure has been used to study the functional connectivity of networks [16,19]. For example, experimental studies show that the second-order effective interaction map of the retina is sparse and dominated by local overlapping effective interaction modules [19]. In this sense, those effective interactions can be regarded as a sparse coding of the information that encoded in the state distribution. Here, the sparseness is defined as the ratio of the number of non-zero effective interactions to $2^n$ (the number of total effective interactions). Several studies show that high-order effective interactions could be important for characterizing observed distributions of network states [6,10,20]. Although high-order effective interactions are required, the number of required non-zero high-order effective interactions are very small in these experimental studies [10]. What leads to the sparsity of effective interactions remains to be clarified.

We address how a sparse anatomical coupling structure (In the following, for simplicity, we use “coupling structure” instead of “anatomical coupling structure”.) could lead to a sparse coding by effective interactions. Network dynamical structure often closely relates to the underlying coupling structure [21]. For example, when the input of each node is independent from others, (i) high-order (≥2) effective interactions are zero in a network of no connections, and (ii) high-order effective interactions are large in a dense and strong connected excitatory network. To efficiently encode information, a realistic system often incorporates a coupling structure with certain features [22,23], e.g., sparsity, small-world, or scale-free. However, it is still unclear how the coupling structure affects the dynamical structure of effective interactions.

In this work, we consider a general class of pulse-coupled networks. The state of each node is binary-state, i.e., active when the node sends pulses to its child nodes, otherwise, silent. We first establish the connection between the coupling structure and the number of non-zero effective interactions in the full-order MEP analysis (constrained by all moments) through an observed Fact—which is independent from node dynamics. We then examine the observed Fact by numerical simulations. Through our analysis, we can estimate the number of non-zero effective interactions for a given coupling structure when the external input of each node is independent from one another. Our results show that a sparse network could lead to a lot of vanishing high-order effective interactions. For illustration, we verify our results by estimating the number of non-zero effective interactions for each order in a network with Erdos–Renyi connection structure, in which our estimated number is much smaller than $C_n^k$, the number of all possible kth-order effective interactions. Our results establish a connection between the dynamical structure and the network coupling structure. This connection provides an insight into how a sparse coupling structure can lead to a sparse coding scheme. In this work, we would mainly use neural networks as examples for illustration, while our results apply to general binary-state networks.

2. Results

In the following analysis, we use binary vector $V\left(l\right)=({\sigma }_1,\cdots ,{\sigma }_n)\in {\{0,1\}}^n$ to represent the state of n nodes within the sampling time bin labeled by l. To obtain correlations up to the mth-order requires to evaluate all ${〈{\sigma }_{i_1}\cdots {\sigma }_{i_M}〉}_E$, where $1\le i_1<i_2<\cdots <i_M\le n$, $1\le M\le m$, and ${〈·〉}_E$ is defined by ${〈g\left(l\right)〉}_E={\sum }_{l=1}^{N_T}g\left(l\right)/N_T$ for any function $g\left(l\right)$ and $N_T$ is the total number of sampling time bins in the recording. The mth-order MEP analysis is to find the desired probability distribution $P\left(V\right)$ for n nodes by maximizing the entropy $S\equiv -{\sum }_{V}P\left(V\right)\log P\left(V\right)$ subject to correlations up to the mth-order ($m\le n$). To solve this optimization problem [24], one introduces Lagrangian multipliers for each constraint, that is, $-J_{i_1\cdots i_k}$ for the constraint of ${〈{\sigma }_{i_1}\cdots {\sigma }_{i_k}〉}_E$ and $(Z-1)$ for the constraint of ${\sum }_{V}P\left(V\right)$ = 1. The optimization problem is

$$P_m\left(V\right)=\arg \underset{P\left(V\right)}{\max }-\sum _{V}P\left(V\right)\log P\left(V\right)-\sum _{k=1}^m\sum _{i_1<\cdots <i_k}^n{F}_{i_1\cdots i_k}-F_0,$$

where

$$F_0=(Z-1)(\sum _{V}P\left(V\right)-1),$$

$$F_{i_1\cdots i_k}=-J_{i_1\cdots i_k}\left(\sum _{V}P\left(V\right)\prod _{j=1}^k{\sigma }_{i_j}\left(V\right)-{〈{\sigma }_{i_1}\cdots {\sigma }_{i_k}〉}_E\right),$$

and ${\sigma }_{i_j}\left(V\right)$ is the state of the $i_j$th node in the network state V. Since the entropy S is a convex function of $P\left(V\right)$, the unique distribution can be solved by taking a derivative of the above objective function with respect to each $P\left(V\right)$ as follows:

$$P_m\left(V\right)=\frac{1}{Z}\exp (\sum _{k=1}^m\sum _{i_1<\cdots <i_k}^n{J}_{i_1\cdots i_k}{\sigma }_{i_1}\cdots {\sigma }_{i_k}),$$

(1)

where, following the terminology of statistical physics, we call $J_{i_1\cdots i_k}$ a kth-order effective interaction ($1\le k\le m$), the partition function Z is the normalization factor. Equation (1) is referred to as the mth-order MEP distribution. When $m<n$, the constants $J_{i_1\cdots i_k}$ can be determined by all constraints through a commonly used iteration method [13]. When $m=n$, all order moments of the experimentally observed distribution are used as constraints; thus, the above optimization problem has only one feasible solution, that is, the experimentally observed distribution [25,26].

First, we discuss the relationship between effective interactions and the statistical distribution of network states. By taking the logarithm of both sides of Equation (1) for $P_n\left(V\right)$, we can get a set of linear equations of all-order effective interactions for all states V. Since $P_n$ is the same as the experimentally observed distribution [25], we can obtain the effective interactions in $P_n$ in terms of the experimentally observed distribution [27]. For example, $n=3$, we can obtain $J_1=\log (P_{100}/P_{000})$ and $J_{12}=\log (P_{110}/P_{010})-J_1$, where $P_{{\sigma }_1{\sigma }_2{\sigma }_3}$ represents the probability of the network state $({\sigma }_1,{\sigma }_2,{\sigma }_3)$. By applying $P({\sigma }_1,{\sigma }_2,{\sigma }_3)=P\left({\sigma }_1\right|{\sigma }_2,{\sigma }_3)P({\sigma }_2,{\sigma }_3)$, we have $J_1=\log \frac{P({\sigma }_1=1|{\sigma }_2=0,{\sigma }_3=0)}{P({\sigma }_1=0|{\sigma }_2=0,{\sigma }_3=0)}$ and $J_{12}=\log \frac{P({\sigma }_1=1|{\sigma }_2=1,{\sigma }_3=0)}{P({\sigma }_1=0|{\sigma }_2=1,{\sigma }_3=0)}-J_1\triangleq J_1^1-J_1$. Our earlier study has shown a recursive structure among effective interactions, that is, the $(k+1)$st-order effective interaction $J_{123\cdots (k+1)}$ can be obtained as follows [27]: first, we switch the state of the $(k+1)$st node in $J_{123\cdots k}$ from silent to active to obtain a new term $J_{123\cdots k}^1$, e.g., from $J_1$ to $J_1^1$; then, we subtract $J_{123\cdots k}$ from the new term to obtain $J_{123\cdots (k+1)}$, i.e.,

$$J_{123\cdots (k+1)}=J_{123\cdots k}^1-J_{123\cdots k}.$$

(2)

Without loss of generality, we randomly select two nodes labeled by 1 and 2. By the recursive relation, any kth-order effective interaction that includes node 1 and 2 can be expressed as the summation of terms with the following basic form:

$$\begin{array}{cc}{J}_{12}^b({\sigma }_3,\cdots ,{\sigma }_n)\hfill & =\log \frac{P({\sigma }_1=1|{\sigma }_2=1,{\sigma }_3,\cdots ,{\sigma }_n)}{P({\sigma }_1=0|{\sigma }_2=1,{\sigma }_3,\cdots ,{\sigma }_n)}\hfill \\ & -\log \frac{P({\sigma }_1=1|{\sigma }_2=0,{\sigma }_3,\cdots ,{\sigma }_n)}{P({\sigma }_1=0|{\sigma }_2=0,{\sigma }_3,\cdots ,{\sigma }_n)}.\hfill \end{array}$$

(3)

Note that we use superscript b to denote “basic”. For any i and j ($i\ne j$), the states of node l ($l\ne i,j$) in $J_{ij}$ is silent. $J_{ij}^b$ is a function of states of nodes ${\sigma }_l$’s ($l\ne i,j$) by replacing the silent state of every node l ($l\ne i,j$) in $J_{ij}$ with variable ${\sigma }_l\in \{0,1\}$. For example, $J_{123}=J_{12}^b(1,0,\cdots ,0)-J_{12}^b(0,0,\cdots ,0)$ and $J_{1234}=[J_{12}^b(1,1,0,\cdots ,0)-J_{12}^b(0,1,0,\cdots ,0)]-J_{123}$. We can observe that if nodes 1 and 2 are independent conditioned on all other nodes, i.e., $P\left({\sigma }_1\right|{\sigma }_2=1,{\sigma }_3,\cdots ,{\sigma }_n)=P\left({\sigma }_1\right|{\sigma }_2=0,{\sigma }_3,\cdots ,{\sigma }_n)$, any effective interaction containing these two nodes is zero.

Next, we would show what kind of coupling structure could entail the conditional independence of two nodes. Here, we define some notations. In any sampling time bin $[0,\Delta )$ with state $V=({\sigma }_1,\cdots ,{\sigma }_n)$, $\forall t\in [0,\Delta )$, we denote $I_{i,t}$ as node i’s input from the outside of the network, denote $w_{ij}\left(t\right)$ as the input from the node i to node j, denote $\mathrm{C}\left(i\right)$ as the set of all child nodes of node i, denote $U_i=\mathrm{C}\left(i\right)\cup \left\{i\right\}$, denote $P\left(e\right)$ as the probability of event e, denote $U_0=\{1,2,\cdots ,n\}$.

Fact 1.

For n pulse-coupled nodes with binary-state dynamics on a network with a coupling structure $G_0$, in any sampling time bin $[0,\Delta )$, $\forall t\in [0,\Delta )$, $\forall i_1,j_1\in U_0$, we assume that: (a) the external inputs of each node are independent from one another, i.e., $P(I_{i_1,t},I_{j_1,t})=P\left(I_{i_1,t}\right)P\left(I_{j_1,t}\right)$; (b) whether a parent node sends spikes to its child nodes only depends on its state, i.e., $P(w_{i_1{j}_1}\left(t\right),V)=W({\sigma }_{i_1},i_1,j_1,t)$, where $W(·,·,·,·)$ is a real function. $\forall i,j\in U_0$, if they neither are connected nor share any common child node, i.e., $U_i\cap U_j=\varphi$, then, node i and j are independent conditioned on the state of all other nodes, i.e.,

$$P({\sigma }_i,{\sigma }_j|H)=P\left({\sigma }_i\right|H)P\left({\sigma }_j\right|H),$$

(4)

where H is a possible state of nodes in $U_0\setminus \{i,j\}$.

We justify our two assumptions as follows: to avoid the influence of correlation in external inputs when we are studying the relation between the dynamical structure and the coupling structure, we assume that the external input of each node is independent from one another, i.e., assumption (a). The second assumption implicates a Markov-like property; that is, for a connected pair of pulse-coupled nodes in an equilibrium state, the pulse from the parent node to the child node only depends on the state of the parent node but is independent from inputs imposed on the parent node. For example, in neural networks, a neuron sends out spikes only when this neuron is active, regardless of what inputs are imposed on the neuron.

The argument for the conclusion in Equation (4) is as follows. By assumption (a), node i and node j can be dependent only through the coupling structure $G_0$. When we are considering how node i and node j affect each other by changing their states through the coupling structure $G_0$, we can consider a simplified coupling structure, $G_1$, which ignores those connections that are independent from states of node i and node j, i.e., ${\sigma }_i$ and ${\sigma }_j$. $\forall k\in U_o\setminus \{i,j\}$, i.e., any other node k, its state ${\sigma }_k$ is fixed when we are considering the conditional probability in Equation (4). By assumption (b), for node k’s any child node l, the input from node k to node l is independent from ${\sigma }_i$ and ${\sigma }_j$. Thus, the connections started from those nodes in $U_o\setminus \{i,j\}$ are fixed for different states of ${\sigma }_i$ and ${\sigma }_j$. Therefore, $G_1$ is a simplified coupling structure that only keeps those connections that originated from node i and node j in $G_0$. In $G_1$, any connection only exists in either sub-network $U_i$ or sub-network $U_j$. Under the condition $U_i\cap U_j=\varphi$, i.e., they neither are connected nor share any common child node, sub-network $U_i$ and sub-network $U_j$ are two isolated sub-networks. ${\sigma }_i$ and ${\sigma }_j$ cannot affect each other by changing their states through the coupling structure $G_1$, that is, node i and j are independent conditioned on the states of all other nodes.

Figure 1 displays an example to illustrate our observed Fact. The coupling structure $G_0$ is shown in Figure 1a. We focus on node 1 and node 2, where they neither are connected nor share any child node. When the state of other nodes (black) are fixed, all outputs from black nodes can be ignored in the simplified coupling structure $G_1$, as shown in Figure 1b. Node 1 and node 2 respectively belong to two separate sub-networks. Therefore, nodes 1 and node 2 are independent conditioned on the state of all other nodes.

Based on the recursive structure of effective interactions and the observed Fact, we reach the following conclusion: with the two assumptions in the observed Fact, for a group of nodes $\{i_1,i_2,\cdots ,i_k\}$, if there exists at least one pair of nodes that neither are connected nor share any child node, effective interaction $J_{i_1,i_2,\cdots ,i_k}$ is zero.

The system we would use to examine our conclusion is an integrate-and-fire (I&F) network, a general pulse-coupled network, with both excitatory and inhibitory nodes [21]. For the ith node, the dynamics of its state variable $x_i$ with time scales $\tau$ are governed by

$${\dot{x}}_i=-\frac{x_i}{\tau }-(g_i^{\mathrm{bg}}+g^{\mathrm{ex}})(x_i-x_{\mathrm{ex}})-g_i^{\mathrm{in}}(x_i-x_{\mathrm{in}}),$$

(5)

where $x_{\mathrm{ex}}$ and $x_{\mathrm{in}}$ are the reversal values of excitation (ex) and inhibition (in), respectively. $g_i^{\mathrm{bg}}=f{\sum }_k\mathbb{H}(t-T_{i,k}^F)\exp [-(t-T_{i,k}^F)/{\sigma }^{\mathrm{ex}}]$ is the background input with magnitude f and time scale ${\sigma }^{\mathrm{ex}}$, $T_{i,k}^F$ is a Poisson process with rate $\mu$, $\mathbb{H}(·)$ is the Heaviside function, $g_i^{\mathrm{ex}}={\sum }_j{\sum }_k{S}_{ij}^{\mathrm{ex}}\mathbb{H}(t-T_{j,k}^{\mathrm{ex}})\exp [-(t-T_{j,k}^{\mathrm{ex}})/{\sigma }^{\mathrm{ex}}]$ is the excitatory input from other jth excitatory nodes, and $g_i^{\mathrm{in}}={\sum }_j{\sum }_k{S}_{ij}^{\mathrm{in}}\mathbb{H}(t-T_{j,k}^{\mathrm{in}})\exp [-(t-T_{j,k}^{\mathrm{in}})/{\sigma }^{\mathrm{in}}]$ is the inhibitory input from other jth inhibitory nodes. The jth excitatory (inhibitory) node $x_j$ evolves continuously according to Equation (5) until it reaches a firing threshold $x_{\mathrm{th}}$. That moment in time is referred to as a firing event (say, the kth spike) and denoted by $T_{j,k}^{\mathrm{ex}}$ ($T_{j,k}^{\mathrm{in}}).$ Then, $x_j$ is reset to the reset value $x_r$ ($x_{\mathrm{in}}<x_r<x_{\mathrm{th}}<x_{\mathrm{ex}}$) and held $x_r$ for an absolute refractory period of ${\tau }_{\mathrm{ref}}.$ Each spike emerging from the jth excitatory (inhibitory) node causes an instantaneous increase $S_{ij}^{\mathrm{ex}}$ ($S_{ij}^{\mathrm{in}}$) in $g_i^{\mathrm{ex}}$ ($g_i^{\mathrm{in}}$), where $S_{ij}^{\mathrm{ex}}$ and $S_{ij}^{\mathrm{in}}$ are the excitatory and inhibitory coupling strengths, respectively. The model (5) describes a general class of physical networks [1,3,21,28,29]. To be intuitive, the I&F model described in Equation (5) can be understood through a resistor–capacitor circuit. Each neuron is a leaky circuit, which consists of a capacitor with dimensionless capacitance as 1 in parallel with two resistors. $x_i$ is the voltage. ${\dot{x}}_i$ on the left-hand side is the current which passes through the capacitor. On the right-hand side, $-x_i/\tau$ is the leaky current. The first resistor has a reversal potential $x_{ex}$ with conductance $(g_i^{\mathrm{bg}}+g^{\mathrm{ex}})$. Therefore, the second term on the right-hand side obtained by the Ohm’s law is the current that passes through the first resistor. Note that the conductance is affected by inputs. It is similar for the third term on the right-hand side, i.e., the current that passes through the second resistor.

In the first example, two excitatory and two inhibitory I&F nodes form a ring coupling structure (Figure 2a). For any pair of nodes, say, node i and j, we compute ${\Delta }_{ij}\left(H\right)=|P({\sigma }_i=1|{\sigma }_j=1,H)-P({\sigma }_i=1|{\sigma }_j=0,H)|$, where H is a state vector of other two nodes. By our observed Fact, the conditional independent pairs are $(\mathrm{neuron}1,\mathrm{neuron}3)$ and $(\mathrm{neuron}2,\mathrm{neuron}4)$, and other pairs are categorized as dependent pairs. In Figure 2b, the strengths of ${\Delta }_{ij}\left(H\right)$ of independent pairs (green) are almost two orders of magnitude smaller than those of dependent pairs (red). We then shuffle spike trains of each node. We similarly compute ${\Delta }_{ij}\left(H\right)$ for 10 different pieces of shuffled data. Blue dots and cyan dots in Figure 2b are results of all shuffled data of dependent pairs and independent pairs, respectively. The strength of ${\Delta }_{ij}\left(H\right)$ of independent pairs (green)—computed from the observed data—are within the statistical error of shuffled data. We then solve effective interactions in the full-order MEP analysis $P_n$ for this ring network. As shown in Figure 2c, the effective interaction strengths of independent pairs ($J_{24}$ and $J_{13})$ are within the statistical error of shuffled results (red). Since every high-order (≥3) effective interaction includes at least one independent pair of nodes, as predicted, the strengths of all high-order effective interactions are within the statistical error of shuffled results as shown in Figure 2d.

In the second example in the second row in Figure 2, results are similar in that dependent pairs and independent pairs can be identified through our observed Fact, and the strength of any effective interaction that includes the independent pair of nodes (node 1 and node 3) is within the statistical error of shuffled data. In this example, as shown in Figure 2h, $J_{124}$ is very small, i.e., within the statistical error of shuffled results. However, $J_{124}$ does not include the independent pair of nodes (node 1 and node 3); thus, the theoretical estimation of zero-strength effective interactions misses $J_{124}$. This example indicates that our theoretical estimation result may give rise to an upper bound of the number of non-zero effective interactions. In contrast, for a network of all excitatory nodes with the same coupling structure as the one in Figure 1e, $J_{124}$ is significantly larger than zero (not shown). Since the estimation of the strength of high-order effective interactions involves the estimation of high dimensional joint probability distribution, a very long recording due to the curse of dimensionality constrains us from examining ${\Delta }_{ij}\left(H\right)$ for a large network.

Based on the relation between the coupling structure and effective interactions, the number of non-zero high-order effective interactions can be small in a sparsely connected network compared with $C_n^k$, which is the number of all possible kth-order interactions. For example, we estimate the number of each-order non-zero effective interactions in a network with an Erdos–Renyi connection structure. We generate 1000 Erdos-Renyi random networks of 100 nodes (the same connection probability but different random samples). In each network, for any pair of node i and node j, we assign the connection from node i to node j by the following rule: we generate a random number from the uniform distribution on $[0,1]$; if the number is smaller than $0.05$, we assign a connection from node i to node j. As shown in Figure 3, the number of non-zero kth-order ($k>1$) effective interactions is much smaller than $C_{100}^k$ (too large to be shown). The number of high-order effective interactions (order higher than 11th) almost vanishes (order higher than 20th not shown). In this example, the sparseness, i.e., the ratio of the number of non-zero effective interactions to $2^n$, is less than ${10}^{-24}$.

3. Conclusions and Discussion

In summary, we have established a relation between effective interactions in MEP analysis and the anatomical coupling structure (“coupling structure” for simplicity in the following) of pulse-coupled networks to understand how a sparse coupling structure could lead to a sparse coding by effective interactions. For example, the sparseness of the case in Figure 3 is less than ${10}^{-24}$. Since effective interactions characterize how nodes interact with each other in the recording of dynamical data, i.e., network dynamical structure, our study quantitatively displays how the dynamical structure closely relates to the coupling structure.

Even though high-order effective interactions are often much smaller compared with low-order ones [27], it is still unclear why small high-order effective interactions do not accumulate to have a significant effect in a large network [10,30]. For example, it has been shown that MEP distribution with sparse low-order effective interactions—non-zero effective interactions are sparse and vanish when the order is higher than the eighth-order—can well capture the state distribution of 99 ganglion cells in the salamander retina responding to a natural movie clip or natural pixel [10]. In this study, we show that a large amount of effective interactions vanish in a sparse coupling structure, thus rationalizing the absence of the accumulation of high-order interactions for a large network.

Finally, we point out that some important issues remain to be elucidated in the future. First, we have ignored correlations in external inputs when estimating the number of non-zero effective interactions. Correlated inputs can induce non-zero high-order effective interactions [31]. It still needs to be considered how the statistics of inputs affect the sparsity of effective interactions. Second, current algorithms for estimating non-zero effective interactions (not limited to the second-order) for a large network (e.g., ∼100 nodes) are very slow, e.g., Monte Carlo based methods [30,32]. Our undergoing work is developing a fast algorithm that exploits the sparsity of effective interactions. We have seen an indication that the algorithm can work well for an I&F network with sparse coupling structure; however, that work has yet to be fully verified to be conclusive.