Effects of hidden nodes on the reconstruction of bidirectional networks

TL;DR

This paper investigates how hidden nodes affect the ability to reconstruct bidirectional networks from node dynamics, providing analytical insights and numerical verification of these effects.

Contribution

It offers new analytical results explaining the impact of hidden nodes on network reconstruction, addressing a common real-world challenge.

Findings

Hidden nodes significantly distort network reconstruction accuracy.

Analytical formulas quantify the influence of hidden nodes.

Numerical studies confirm theoretical predictions.

Abstract

Much research effort has been devoted to developing methods for reconstructing the links of a network from dynamics of its nodes. Many current methods require the measurements of the dynamics of all the nodes be known. In real-world problems, it is common that either some nodes of a network of interest are unknown or the measurements of some nodes are unavailable. These nodes, either unknown or whose measurements are unavailable, are called hidden nodes. In this paper, we derive analytical results that explain the effects of hidden nodes on the reconstruction of bidirectional networks. These theoretical results and their implications are verified by numerical studies.

Tables

Table 1. Table 1: Error rates of the reconstruction results using k-means clustering of ( Σ m − 1 ) i ​ j subscript superscript subscript Σ 𝑚 1 𝑖 𝑗 (\Sigma_{m}^{-1})_{ij} for the various networks studied. N = 100 𝑁 100 N=100 for networks A, B, C and E and N = 1000 𝑁 1000 N=1000 for network D. ρ = N L / [ n ​ ( n − 1 ) / 2 ] 𝜌 subscript 𝑁 𝐿 delimited-[] 𝑛 𝑛 1 2 \rho=N_{L}/[n(n-1)/2] , where N L subscript 𝑁 𝐿 N_{L} is the number of links among the measured nodes and n = N − n h 𝑛 𝑁 subscript 𝑛 ℎ n=N-n_{h} is the number of measured nodes. Two clusters are used for all networks except network C where three clusters are used. We show also the error rates using Σ i ​ j − 1 subscript superscript Σ 1 𝑖 𝑗 \Sigma^{-1}_{ij} in parentheses.

network	dynamics	$n_h$	$\rho$	FN/$N_L$ (%)	FP/$N_L$ (%)
A	logistic	$50$	0.202	0.81 (0.81)	0.00 (0.00)
B	logistic	$10$	0.198	$0.76\left(0.88\right)$	$0.00\left(0.00\right)$
		$30$	0.197	$1.05\left(1.26\right)$	$0.00\left(0.00\right)$
		$50$	0.202	$1.62\left(1.62\right)$	$0.00\left(0.00\right)$
		$70$	0.218	$1.05\left(3.16\right)$	$0.00\left(0.00\right)$
	FHN	$10$	0.198	$0.63\left(0.51\right)$	$0.00\left(0.00\right)$
		$30$	0.197	$1.05\left(0.84\right)$	$0.00\left(0.00\right)$
		$50$	0.202	$1.62\left(1.21\right)$	$0.00\left(0.00\right)$
		$70$	0.218	$3.16\left(1.05\right)$	$1.05\left(0.00\right)$
	Rössler	$10$	0.198	$1.01\left(0.88\right)$	$0.00\left(0.00\right)$
		$30$	0.197	$1.10\left(0.44\right)$	$0.00\left(0.00\right)$
		$50$	0.202	$0.87\left(1.31\right)$	$0.00\left(0.00\right)$
		$70$	0.218	$3.57\left(0.00\right)$	$5.95\left(0.00\right)$
C	logistic	$20$	0.198	$0.96\left(0.80\right)$	$0.00\left(0.00\right)$
		$30$	0.197	$0.85\left(0.21\right)$	$0.00\left(0.00\right)$
		$40$	0.199	$0.85\left(0.00\right)$	$1.14\left(0.00\right)$
		60	0.206	2.48 (0.00)	1.86 (0.00)
D	logistic	100	0.0039	0.44 (0.50)	0.00 (0.00)
		300	0.0037	0.45 (0.45)	0.11 (0.00)
		500	0.0040	0.60 (0.60)	1.20 (0.00)
		700	0.0017	1.05 (1.05)	3.16 (0.00)
E	logistic	$30$	0.199	$0.42\left(0.42\right)$	$29.88\left(0.00\right)$

Equations

$$P_G(x_1,\dots ,x_N)=\frac{\exp \left[-(\mathbf{x}-\mathbf{m}{)}^T{\mathbf{\Sigma }}^{-1}(\mathbf{x}-\mathbf{m})/2\right]}{(2\pi {)}^{N/2}\sqrt{\det (\mathbf{\Sigma })}}$$

$$\frac{dx_{i}}{dt}=f(x_{i})+\sum_{j\neq i}g_{ij}A_{ij}h(x_{j}-x_{i})+\eta_{i}.\$$

$$\overline{{\eta }_i(t){\eta }_j(t^{\prime })}={\sigma }_n^2{\delta }_{ij}\delta (t-t^{\prime }),$$

$$L_{ij}=s_i{\delta }_{ij}-g_{ij}{A}_{ij},s_i\equiv \sum _{k=1}^N{g}_{ik}{A}_{ik}$$

$$\frac{d}{dt}\delta x_i=-\sum _j\left(L_{ij}+a{\delta }_{ij}\right)\delta x_j+{\eta }_i$$

$$\mathbf{\Sigma }\equiv \underset{t\to \infty }{\lim }\overline{[\mathbf{x}(t)-\overline{\mathbf{x}(t)}][\mathbf{x}(t)-\overline{\mathbf{x}(t)}{]}^T}$$

$${\mathbf{\Sigma }}^{-1}=\frac{2}{{\sigma }_n^2}(\mathbf{L}+a{\mathbf{I}}_N)$$

$${\mathbf{\Sigma }}_{ij}^{-1}=-\frac{2}{{\sigma }_n^2}{g}_{ij}{A}_{ij},i\ne j$$

$${\Sigma }_{ij}\approx ⟨[x_i(t)-⟨x_i(t)⟩][x_j(t)-⟨x_j(t)⟩]⟩$$

$${\mathbf{\Sigma }}_m\equiv \underset{t\to \infty }{\lim }\overline{[\mathbf{y}(t)-\overline{\mathbf{y}(t)}][\mathbf{y}(t)-\overline{\mathbf{y}(t)}{]}^T}$$

$$({\Sigma }_m{)}_{ij}\approx ⟨[y_i(t)-⟨y_i(t)⟩][y_j(t)-⟨y_j(t)⟩]⟩$$

$$\mathbf{\Sigma }=\left(\begin{array}{cc}{\mathbf{\Sigma }}_m& \mathbf{U}\\ {\mathbf{U}}^{\mathbf{T}}& {\mathbf{\Sigma }}_h\end{array}\right)$$

$$({\Sigma }_h{)}_{\mu \nu }=\underset{t\to \infty }{\lim }\overline{[x_{\mu +n}(t)-\overline{x_{\mu +n}(t)}][x_{\nu +n}(t)-\overline{x_{\nu +n}(t)}]}$$

$$U_{i\mu }=\underset{t\to \infty }{\lim }\overline{[y_i(t)-\overline{y_i(t)}][x_{\mu +n}(t)-\overline{x_{\mu +n}(t)}]}$$

$$\mathbf{L}=\left(\begin{array}{cc}{\mathbf{L}}_m& \mathbf{E}\\ {\mathbf{E}}^T& {\mathbf{L}}_h\end{array}\right)$$

$$E_{i\mu }=-g_{i,\mu +n}{A}_{i,\mu +n}$$

$${\mathbf{I}}_n$$

$$0$$

$${\mathbf{\Sigma }}_m^{-1}=\frac{2}{{\sigma }_n^2}\left[{\mathbf{L}}_m+a{\mathbf{I}}_n-\mathbf{E}{\left({\mathbf{L}}_h+a{\mathbf{I}}_{n_h}\right)}^{-1}{\mathbf{E}}^T\right]$$

$$\mathbf{C}\equiv \mathbf{E}{\left({\mathbf{L}}_h+a{\mathbf{I}}_{n_h}\right)}^{-1}{\mathbf{E}}^T$$

$${\left({\Sigma }_m^{-1}\right)}_{ij}=-\frac{2}{{\sigma }_n^2}\left(g_{ij}{A}_{ij}+C_{ij}\right),i\ne j$$

$${\left({\Sigma }_m^{-1}\right)}_{ij}={\Sigma }_{ij}^{-1}-\frac{2}{{\sigma }_n^2}{C}_{ij}i,j=1,\dots ,n$$

$$S_{\mu \nu }$$

$$W_{\mu \nu }$$

$$C_{ij}=\sum _{{\mu }_1=1}^{n_h}{F}_0{W}_{i{\mu }_1}{W}_{{\mu }_{1}j}+\sum _{{\mu }_1,{\mu }_2=1}^{n_h}{F}_1{W}_{i{\mu }_1}{W}_{{\mu }_1{\mu }_2}{W}_{{\mu }_{2}j}$$

$$+\sum _{k=3}^{n_h-1}\sum _{{\mu }_1,\cdots ,{\mu }_{k+1}=1}^{n_h}{F}_k{W}_{i{\mu }_1}\left(\prod _{j=1}^k{W}_{{\mu }_j{\mu }_{j+1}}\right)W_{{\mu }_{k+1}j}$$

$$P({\Sigma }_{ij}^{-1}=x)=(1-r)P_0(x)+r{P}_1(x)$$

$$m_1=-\frac{2}{{\sigma }_n^2}⟨g⟩;{\sigma }_1=\frac{2}{{\sigma }_n^2}{\sigma }_g$$

$$P(({\Sigma }_m^{-1}{)}_{ij}=x)=(1-r){\tilde{P}}_0(x)+r{\tilde{P}}_1(x),i\ne j$$

$${\tilde{m}}_0$$

$${\tilde{\sigma }}_0^2$$

$${\tilde{m}}_1$$

$${\tilde{\sigma }}_1^2$$

Full text

I Introduction

Many systems of interest in physics and biology are represented by complex networks of a large number of elementary units or nodes that interact or link with each other Strogatz . A substantial amount of data have been obtained for various networks especially biological networks, and a grand challenge is to reveal the structure of these networks, namely the links, their direction and relative coupling strength, from the measured data. It is expected that Timme0 the structure of a network controls its dynamics and thus one might be able to uncover information about the structure of a network from its dynamics. Much research effort has been devoted to developing methods for reconstructing a network from the dynamics of the nodes (see e.g. TimmeReview ; PhysReportYCLai2016 for review). Counterintuitively, it has been demonstrated that the presence of noise acting on the network can be beneficial for network reconstruction as the noise induces a relation between measurable quantities from dynamics and the network structure noise . Making use of different relations of this kind, a number of methods PRE ; PRErapid ; SolvingPRE ; PRErapid2 ; LaiPRE ; ChenEPL ; correlated have been proposed for reconstructing networks solely from the dynamics of the nodes. In all these methods, in order to calculate the quantities that are related to the network structure, the measurements of the dynamics of all the nodes are required.

In real-world problems, it is common that either some nodes of a network of interest are unknown or the measurements of some nodes are unavailable. These nodes, either unknown or whose measurements are unavailable, are called hidden nodes. It is thus important to study and understand the effects of hidden nodes on the reconstruction of networks DunnPRE2013 ; YCLai2014 ; ECC2014 ; WXWangPRL2015 ; HuangJPhyA2015 ; HuGang2017 ; Marc2017 . This task is highly challenging and, as of today, there is not yet a general and analytical understanding of the effects of hidden nodes.

A usual practice infers links from the correlation of the measurements, with a larger correlation coefficient interpreted as a higher probability of link StuartScience2003 ; BrainPRL ; frontiers . However, correlation between measurements of two nodes cannot be equated with direct interactions between the two nodes. In fact, it has been clearly shown that for networks of neurons, the spiking dynamics of neurons can have weak pairwise correlations even though they are strongly coupled weak . This study further shows that the spiking dynamics of neurons are quantitatively captured by the probability distribution $P_{\rm Ising}(\hat{\sigma}_{1},\ldots,\hat{\sigma}_{N})\propto\exp[\sum_{i<j}^{N}J_{ij}\hat{\sigma}_{i}\hat{\sigma}_{j}+\sum_{i}h_{i}\hat{\sigma}_{i}]$ of an Ising model, which is the maximum entropy distribution of a system of $N$ nodes with binary state variables $\hat{\sigma}_{i}=\pm 1$ that is consistent with the measured averages and covariances. This leads to extensive studies of the inverse Ising problem: the inference of the couplings $J_{ij}$, which are interpreted as effective interactions between nodes $i$ and $j$, from measured averages and covariances of the data (see e.g. Advance for review).

In general, systems of interest have nodes with continuous state variables. In this case, the maximum entropy distribution consistent with the measured averages and covariances is the multivariate Gaussian distribution Conrad

[TABLE]

where $\bm{x}=(x_{1},x_{2},\ldots,x_{N})^{T}$, $\bm{m}=(m_{1},\ldots,m_{N})^{T}$ with $m_{i}$ being the measured average of $x_{i}$, $\bm{\Sigma}$ is the measured covariance matrix, and the superscript $T$ denotes the transpose. Equation (1) suggests that information of links in bidirectional networks is contained in $\bm{\Sigma}^{-1}$, the inverse of the covariance matrix. A mathematical relation between the weighted adjacency matrix of the network and the inverse of the covariance matrix has indeed been derived PRErapid ; correlated for a model class of bidirectional networks with diffusive-like coupling and subjected to Gaussian white noise. Based on this theoretical result, a method that reconstructs bidirectional networks from $\bm{\Sigma}^{-1}$ has been developed correlated .

In this paper, we address the question of how hidden nodes affect the reconstruction of bidirectional networks. Using the same model class of bidirectional networks studied in PRErapid ; correlated , we derive analytical results relating quantities involving the hidden nodes to the inverse of the covariance matrix of the measured data and use these results to explain the effects of hidden nodes on the reconstruction results. Then we carry out numerical studies to verify these theoretical results and their implications.

II Formulation of the problem

We consider weighted bidirectional networks of $N$ nodes with nonlinear dynamics and diffusive-like coupling. The dynamics of each node is described by a variable $x_{i}(t)$, $i=1,2,\ldots,N$, and the time evolution of $x_{i}(t)$ is given by

[TABLE]

The adjacency matrix element $A_{ij}$ is $1$ when node $j$ is linked to node $i$ by the diffusive-like coupling function $h$ with coupling strength $g_{ij}$; otherwise $A_{ij}=g_{ij}=0$. The coupling is bidirectional so $A_{ij}=A_{ji}$ and $g_{ij}=g_{ji}$, and the graph of the networks has no self-loops, that is, $A_{ii}\equiv 0$. As discussed PRErapid ; correlated , $f$ describes the intrinsic dynamics that is generally nonlinear and identical for all the nodes, and the diffusive-like coupling function $h$ satisfies $h(-z)=-h(z)$ and $h^{\prime}(0)>0$. Thus excitatory or activating links have with $g_{ij}>0$ whereas inhibitory links have $g_{ij}<0$. Here we take $h^{\prime}(0)=1$. External influences are modelled by a Gaussian white noise $\eta$ with zero mean and variance $\sigma_{n}^{2}$:

[TABLE]

where the overbar is an ensemble average over different realizations of the noise.

The weighted Laplacian matrix of the network, ${\bm{L}}$, is given by

[TABLE]

and contains connectivity information of the network. Here $s_{i}$ is the weighted degree or the strength of node $i$. For these networks, $x_{i}(t)$’s approach $X_{0}$ with $f^{\prime}(X_{0})<0$ in the absence of noise. Let $\delta x_{i}(t)=x_{i}(t)-X_{0}$, then the linearized system around the noise-free steady state is given by

[TABLE]

where $a\equiv-{f^{\prime}(X_{0})}>0$. We consider systems that have stationary dynamics and this implies $\bm{L}+a{\bm{I}}$ is positive definite Arnold . Using Eq. (5), it has been derived correlated that the covariance matrix $\bm{\Sigma}$, defined by

[TABLE]

is related to $\bm{L}$ by

[TABLE]

where ${\bm{I}}_{N}$ is the $N\times N$ identity matrix. Equation (7) is an exact result for the linearized system [Eq. (5)] and a good approximation for the original nonlinear network [governed by Eq. (2)] when the noise is weak. All the theoretical results presented in this paper should be understood in this manner. An important consequence of Eq. (7) is

[TABLE]

which indicates that the off-diagonal elements of $\Sigma^{-1}_{ij}$ would separate into two groups according to $A_{ij}=0$ or 1. Making use of this result, a reconstruction method of the adjacency matrix and thus the links of the network by performing clustering analysis of the off-diagonal elements of $\bm{\Sigma}^{-1}$ has been developed correlated . In this method, ${\bm{\Sigma}}$ is evaluated by approximating the ensemble average by time average:

[TABLE]

where $\langle\cdots\rangle$ denotes a time average. To evaluate $\bm{\Sigma}^{-1}$, the measurements of $\bm{x}(t)$ from all the $N$ nodes are required.

We study the problem of network reconstruction when there are $n_{h}$ hidden nodes and only measurements from $n=N-n_{h}<N$ nodes are available. We call these $n$ nodes the measured nodes and, for clarity, denote their measured dynamics by $y_{i}(t)$, $i=1,2,\ldots,n$ and the corresponding covariance matrix of the measured data by ${\bm{\Sigma}}_{m}$:

[TABLE]

where $\bm{y}(t)\equiv(y_{1},y_{2},\cdots,y_{n})^{T}$. Similarly,

[TABLE]

We would like to answer the following questions: How would the hidden nodes affect the reconstruction results based on ${\bm{\Sigma}}_{m}^{-1}$? Whether and when the links among the $n$ measured nodes can be reconstructed from the measured $y_{i}(t)$’s?

III Theoretical relation for ${\bm{\Sigma}}_{m}^{-1}$

Without loss of generality, we let $y_{i}(t)=x_{i}(t)$, $i=1,2,\ldots,n$. Then we partition ${\bm{\Sigma}}$ into four block matrices

[TABLE]

where the $n\times n$ block matrix is ${\bm{\Sigma}}_{m}$, the $n_{h}\times n_{h}$ block matrix ${\bm{\Sigma}}_{h}$ is the covariance matrix of the hidden nodes given by

[TABLE]

and the $n\times n_{h}$ block matrix ${\bm{U}}$ measures the covariance between the measured and hidden nodes with

[TABLE]

For clarity, we use Roman subscripts $i,j,\ldots$ for the measured nodes and Greek subscripts $\mu,\nu,\ldots$ for the hidden nodes. We partition the weighted Laplacian matrix in a similar fashion:

[TABLE]

The $n\times n$ block matrix ${\bm{L}}_{m}$ and the $n_{h}\times n_{h}$ block matrix ${\bm{L}}_{h}$, with elements $(L_{m})_{ij}=L_{ij}$ and $(L_{h})_{\mu\nu}=L_{\mu+n,\nu+n}$, contain information of the connectivity among the measured nodes and among the hidden nodes respectively while the $n\times n_{h}$ block matrix ${\bm{E}}$ contains information of the connectivity between the measured and hidden nodes with elements

[TABLE]

Using Eq. (7), we obtain

[TABLE]

which imply

[TABLE]

As ${\bm{L}}+a{\bm{I}}_{N}$ is positive definite, ${\bm{L}}_{h}+a{\bm{I}}_{n_{h}}$ is also positive definite and is thus invertible. We define

[TABLE]

then the off-diagonal elements of Eq. (19) can be written as

[TABLE]

which shows that the hidden nodes affect the reconstruction results based on ${\bm{\Sigma}}_{m}^{-1}$ by introducing corrections given by $\bm{C}$. This can be seen directly by using Eq. (7) to rewrite Eq. (19) as

[TABLE]

Equations (19) and (22) are our major theoretical results and we shall use them to answer the questions of interest.

IV Corrections due to hidden nodes

From Eq. (20) and using Eq. (16), we immediately see that $C_{ij}=0$ when $g_{i,\mu+n}=0$ or $g_{j,\mu+n}=0$ for all $\mu=1,2,\ldots,n_{h}$, that is when at least one of the measured nodes $i$ and $j$ is not connected to any hidden node. We let $\bm{M}\equiv\bm{L}_{h}+a\bm{I}_{n_{h}}\equiv{\bm{S}}-{\bm{W}}$ where ${\bm{S}}$ and ${\bm{W}}$ are defined by

[TABLE]

for $\mu,\nu=1,2,\ldots,n_{h}$. $\bm{S}$ is a diagonal matrix with the diagonal elements related to the strength of the hidden nodes and ${\bm{W}}$ is the weighted adjacency matrix of the hidden nodes. Then we obtain (see Appendix)

[TABLE]

where $F_{k}$ generally depends on $s_{\mu_{l}+n}+a$, $l=2,\ldots,k$ and the eigenvalues of $\bm{M}$. We have also similarly defined $W_{i\mu}\equiv g_{i,\mu+n}A_{i,\mu+n}$ for $i=1,2,\ldots,n$, $\mu=1,2,\ldots,n_{h}$ and $W_{\mu i}=W_{i\mu}$. The product $W_{i\mu_{1}}W_{\mu_{1}\mu_{2}}\cdots W_{\mu_{m-1}\mu_{m}}W_{\mu_{m}j}$ is nonzero only if there is a path connecting the measured node $i$ and the measured node $j$ via hidden nodes $\mu_{1},\mu_{2},\ldots,\mu_{m-1},\mu_{m}$ (see Fig. 1). Thus $C_{ij}$ would be zero if there does not exist any path connecting the measured nodes $i$ and $j$ via hidden nodes only. In general, we expect $C_{ij}$ to be nonzero for some pairs of measured nodes $i$ and $j$.

As indicated by Eq. (8), the distribution of the off-diagonal elements of $\Sigma^{-1}_{ij}$ for $i,j=1,2,\ldots,n$ can be written as

[TABLE]

where $r$ is the fraction of connected pairs among the measured nodes and is equal to the number of connected pairs of measured nodes divided by $n(n-1)/2$, $P_{0}$ and $P_{1}$ are the distributions of $\Sigma^{-1}_{ij}$ with $A_{ij}=0$ and $A_{ij}=1$ respectively. If the positive and negative coupling strength of the links are described by two different distributions, $P_{1}$ can further be a mixture of two distributions $P_{1+}$ and $P_{1-}$, which correspond to $g_{ij}>0$ and $g_{ij}<0$ respectively [c.f. Eq. (48)]. In the limit of infinite number of data, $P_{0}(x)$ approaches $\delta(x)$. For a finite number of data points, numerical studies correlated show that $P_{0}(x)$ is well-approximated by a Gaussian distribution of mean $m_{0}=0$ and standard deviation $\sigma_{0}$, and $\sigma_{0}$ decreases when the number of data points increases. Equation (8) implies that $P_{1}(x)$ would have a mean $m_{1}$ and standard deviation $\sigma_{1}$ given by

[TABLE]

where $\langle g\rangle$ and $\sigma_{g}$ are the average and standard deviation of the coupling strength $g_{ij}$ of the links. In the presence of hidden nodes, the corrections $C_{ij}$ will modify the distributions $P_{0}$ and $P_{1}$ to ${\tilde{P}}_{0}$ and ${\tilde{P}}_{1}$:

[TABLE]

The mean and standard deviation of ${\tilde{P}}_{i}$, denoted by $\tilde{m}_{i}$, $\tilde{\sigma}_{i}$, $i=0$ and $1$, would be modified. Using Eq. (21) and the results for $P_{0}(x)$, we obtain

[TABLE]

where

[TABLE]

Here $\mu_{C,i}$ and $\sigma_{C,i}$ for $i=1,2$ are the conditional average and conditional standard deviation of $C_{ij}$ for $A_{ij}=0$ and $A_{ij}=1$ respectively, and $K(g_{ij},C_{ij})$ measures the correlation of $C_{ij}$ with the coupling strength $g_{ij}$ for measured nodes $i$ and $j$ that are connected. Hence, the corrections $C_{ij}$ would shift the means, broaden and distort the distributions of $P_{0}$ and $P_{1}$ to $\tilde{P}_{0}$ and $\tilde{P}_{1}$.

Suppose $P_{0}$ and $P_{1}$ are distinguishable. If ${\tilde{P}}_{0}$ and ${\tilde{P}}_{1}$ remain distinguishable even though with a larger extent of overlap, then the links among the measured nodes can still be reconstructed amid with a larger error rate. For a mixture of two general distributions, there is no simple criterion on when the component distributions are distinguishable. Nonetheless, the component distributions are likely to be distinguishable when the absolute value of the difference between their means are larger than a certain multiple of the sum of their standard deviations. Let $m_{1}-m_{0}=\gamma(\sigma_{0}+\sigma_{1})$ with $|\gamma|>1$. Using Eqs. (29)-(32), we obtain

[TABLE]

Thus ${\tilde{P}}_{0}$ and ${\tilde{P}}_{1}$ are likely to remain distinguishable (with $|{\tilde{m}}_{1}-{\tilde{m}}_{0}|/({\tilde{\sigma}}_{0}+{\tilde{\sigma}}_{1})>1$) if $|\mu_{C,1}-\mu_{C,0}|$, $\sigma_{C,0}$, $\sigma_{C,1}$ and $|K(g_{ij},C_{ij})|$ are sufficiently small.

To shed further light on this, we use Eq. (20) to obtain a crude estimate of $C_{ij}$. Substitute ${\bm{L}}_{h}+a{\bm{I}}_{n_{h}}={\bm{S}}-{\bm{W}}={\bm{S}}({\bm{I}}_{n_{h}}-{\bm{S}}^{-1}{\bm{W}})$ into Eq. (20), we have

[TABLE]

If the Neumann series $\sum_{k=0}^{\infty}({\bm{S}}^{-1}{\bm{W}})^{k}$ converges, it converges to $({\bm{I}}_{n_{h}}-{\bm{S}}^{-1}{\bm{W}})^{-1}$. The necessary and sufficient condition for the Neumann series to converge is the spectral radius of ${\bm{S}}^{-1}{\bm{W}}$, denoted by $\rho({\bm{S}}^{-1}{\bm{W}})$, is less than 1 note . For networks with $g_{ij}\geq 0$, one can easily show that the infinity norm $||{\bm{S}}^{-1}{\bm{W}}||_{\infty}\equiv\max_{\mu}\{\sum_{\nu}|({\bm{S}}^{-1}{\bm{W}})_{\mu\nu}|\}<1$, and thus $\rho({\bm{S}}^{-1}{\bm{W}})\leq||{\bm{S}}^{-1}{\bm{W}}||_{\infty}<1$. For networks with both positive and negative $g_{ij}$, we have checked numerically that $\rho({\bm{S}}^{-1}{\bm{W}})<1$ for all the cases studied. When the Neumann series converges, we keep two terms in the series, namely $({\bm{I}}_{n_{h}}-{\bm{S}}^{-1}{\bm{W}})^{-1}\sim{\bm{I}}_{n_{h}}+{\bm{S}}^{-1}{\bm{W}}$, to obtain a crude estimate of $C_{ij}$:

[TABLE]

Using Eq. (40), one sees that the magnitude of $C_{ij}$ depends on three factors: (1) the number of paths connecting the measured nodes $i$ and $j$ via the hidden nodes which determines the number of nonzero terms in the sums, (2) the strength of the hidden nodes and (3) the coupling strength of the links in these paths. Regarding to the second factor, we note that hidden nodes with larger strength actually give rise to smaller corrections in contrary to what one might have guessed. If these factors do not differ much between the two groups of unconnected or connected measured nodes, then $\mu_{C,1}\sim\mu_{C,0}$; if these factors do not vary much among the measured nodes in each group, then $\sigma_{C,0}$ and $\sigma_{C,1}$ would be small and if these factors do not correlate with the magnitude of $g_{ij}$ for connected measured nodes, then $|K(g_{ij},C_{ij})|$ would be small even when the magnitudes of the corrections $C_{ij}$’s themselves might be large. Such situations are expected when the the hidden nodes are not preferentially linked to the measured nodes in any manner. In this case, ${\tilde{P}}_{0}$ and ${\tilde{P}}_{1}$ remain distinguishable and it is possible to reconstruct the links among the measured nodes from $(\Sigma_{m}^{-1})_{ij}$, $i\neq j$.

V Numerical results and discussions

We check our theoretical results using data from numerical simulations. We study five different networks, four of $N=100$ each and one of $N=1000$.

• (1)

Network A: it consists of two random networks, each of 50 nodes and a connection probability of 0.2, connected to each other by one link and $g_{ij}$’s, taken from a Gaussian distribution $\mathcal{N}(10,2^{2})$ of mean 10 and standard deviation 2, are all positive. We take all the 50 nodes of one of the random network as hidden nodes.

• (2)

Network B: it is a random network of connection probability 0.2 and $g_{ij}$’s also taken from $\mathcal{N}(10,2^{2})$ and are all positive. We choose the hidden nodes randomly from the network with $n_{h}\leq 70$ such that the number of links among the measured nodes is at least of the order of 100.

• (3)

Network C: it is similar to network B except that $g_{ij}$ of 80% of the links taken from $\mathcal{N}(10,2^{2})$ and the remaining 20% taken from $\mathcal{N}(-10,2^{2})$. As a result, about 80% of the $g_{ij}$’s are positive and about 20% are negative. The hidden nodes are chosen randomly from the network.

• (4)

Network D: it is a scale-free network of $N=1000$ BA with degree distribution obeying a power law and $g_{ij}$’s taken from $\mathcal{N}(10,2^{2})$ are all positive. The hidden nodes are chosen randomly from the network.

• (5)

Network E: it is constructed by linking 30 additional nodes to a random network of $70$ nodes and connection probability 0.2 with the restriction that every one of the additional nodes is only commonly connected to randomly selected pairs of unconnected nodes in the random network; and the additional nodes are randomly connected among themselves with the same connection probability 0.2. $g_{ij}$’s are taken from $\mathcal{N}(10,2^{2})$ and are all positive. We take the 30 additional nodes as hidden nodes.

For the dynamics, we mainly study nonlinear logistic function

[TABLE]

and diffusive coupling function

[TABLE]

and take $\sigma_{n}=1$ for the noise. To explore how general our theoretical results are, we go beyond the description by Eq. (2) and study two additional cases. In the first additional case, the nodes of network B have two-dimensional state variables $(x_{i}(t),y_{i}(t))$ with nonlinear FitzHugh-Nagumo (FHN) dynamics FHN

[TABLE]

where $\epsilon=0.01$ and $\alpha=0.95$. In the second additional case, the nodes of network B have three-dimensional state variables $(x_{i}(t),y_{i}(t),z_{i}(t))$ with nonlinear Rössler dynamics Rossler and nonlinear coupling PRErapid :

[TABLE]

where $c_{1}=c_{2}=0.2$ and $c_{3}=9$. In these two additional cases, the system does not approach a steady state in the absence of noise, and has chaotic dynamics when the nodes are decoupled in the second case with Rössler dynamics. We integrate the equations of motion using the Euler-Maruyama method and record the time series $x_{i}(t)$ with a sampling interval $\delta t=5\times 10^{-4}$. For all the cases studied, including the cases with FHN and Rössler dynamics, we calculate ${\bm{\Sigma}}$ using $x_{i}(t)$’s with a time average over $N_{\rm data}=2\times 10^{6}$ data points.

For network A, since there is only one link connecting the hidden nodes and the measured nodes, there is no path connecting any pair of measured nodes via the hidden nodes thus $C_{ij}=0$ for all $i\neq j$. As a result, Eq. (22) implies that $(\Sigma_{m}^{-1})_{ij}=\Sigma_{ij}^{-1}$ for $i\neq j$. We show the distributions of $P(\Sigma_{ij}^{-1})$ and $P((\Sigma_{m}^{-1})_{ij})$ for $i\neq j=1,2,\ldots,n=50$ in Fig. 2. As expected, the two distributions coincide with each other. Moreover, $P(\Sigma_{ij}^{-1})$ is bimodal with the peak around zero corresponding to $P_{0}$ for unconnected nodes and the peak around $x_{m}\approx-20$ corresponding to $P_{1}$ for connected nodes in accord with Eq. (26). Furthermore, the value of $x_{m}$ is in excellent agreement with the theoretical value of $\mu_{1}=-2\langle g\rangle/\sigma_{n}^{2}$ [see Eq. (27)]. Hence in this case, the links among the measured nodes can be reconstructed from $(\Sigma_{m}^{-1})_{ij}$ with $i\neq j$, which can be calculated using the dynamics $y_{i}(t)$ of the measured nodes only.

For network B with the hidden nodes randomly chosen, there are nonzero $C_{ij}$’s for some pairs of measured nodes $i$ and $j$. We first consider the case with logistic dynamics and calculate $C_{ij}$ using Eq. (20) and together with ${\Sigma}^{-1}_{ij}$, we obtain the theoretical results for $(\Sigma_{m}^{-1})_{ij}$ using Eq. (22). We compare these theoretical results with $(\Sigma_{m}^{-1})_{ij}$ directly calculated from the measured dynamics $y_{i}(t)$’s in Fig. 3 and perfect agreement is found for all the values of $n_{h}$ studied. For FHN and Rössler dynamics, the system is not described by Eq. (2) thus $a=-f^{\prime}(X_{0})$ is not defined. We put $a=0$ in Eq. (20) and obtain the theoretical estimate for the off-diagonal $(\Sigma_{m}^{-1})_{ij}$ as $\Sigma^{-1}_{ij}-(2/\sigma_{n}^{2})(EL_{h}^{-1}E^{T})_{ij}$. Interestingly, these theoretical estimates are in good agreement with the directly calculated $(\Sigma_{m}^{-1})_{ij}$’s in most cases, as shown in Figs. 4 and 5. For FHN dynamics with larger $n_{h}$, an improved theoretical estimate is obtained by $\Sigma^{-1}_{ij}-(2/\sigma_{n}^{2})[(EL_{h}^{-1}E^{T})_{ij}+b]$, where $b$ is a constant. This indicates the general applicability of our theoretical results beyond the model class studied.

As shown in Figs. 3-5, the distribution $P((\Sigma_{m}^{-1})_{ij})$ is a mixture of the modified distributions $\tilde{P}_{0}$ and $\tilde{P}_{1}$, in accord with Eq. (28), which remain distinguishable as expected since the hidden nodes are chosen randomly. Thus it is possible to reconstruct the links among the measured nodes from $(\Sigma_{m}^{-1})_{ij}$ with $i\neq j$. We note that this is true for all the three kinds of dynamics studied and even when the hidden nodes outnumber the measured nodes. In Table 1, we compare the error rates of the reconstruction results obtained using k-means clustering of $(\Sigma_{m}^{-1})_{ij}$ from the measured dynamics only and of ${\Sigma}^{-1}_{ij}$ from the dynamics of the whole network. We measure the error rates by the ratios of false negatives (FN) and false positives (FP) over the number of actual links $N_{L}$ among the measured nodes. These error rates are related to the sensitivity and specificity usually used for a predictive test: sensitivity is given by $1-{\rm FN}/N_{L}$ and specificity is given by $1-({\rm FP}/N_{L})\rho/(1-\rho)$, where $\rho=N_{L}/[n(n-1)/2]$ is the link density of the measured nodes. For networks with low link density $\rho$, the error rates can be rather high even when specificity is close to 1 so the error rates are better measures of the accuracy of the reconstruction results PRErapid2 . As can be seen, the accuracy of the reconstruction results using $(\Sigma_{m}^{-1})_{ij}$ from the measured nodes only is comparable to that obtained using $(\Sigma^{-1})_{ij}$ from all the nodes.

For network C, the positive and negative $g_{ij}$’s follow two different distributions so $P_{1}$ can be further decomposed into a weighted sum of $P_{1+}$ and $P_{1-}$, which correspond to $g_{ij}>0$ and $g_{ij}<0$ respectively. Therefore, Eq. (26) can be rewritten as

[TABLE]

where $\beta$ is the fraction of the links among the measured nodes with positive $g_{ij}$’s. Similarly, Eq. (28) is also rewritten as

[TABLE]

As $g_{ij}$ can now be either positive or negative, the terms in the summations contributing to $C_{ij}$ [see Eq. (40)] could cancel one another. Thus we expect the magnitudes of $\mu_{C,0}$, $\mu_{C,1+}\equiv\langle C_{ij}|A_{ij}=1,g_{ij}>0\rangle$ and $\mu_{C,1-}\equiv\langle C_{ij}|A_{ij}=1,g_{ij}<0\rangle$ to be smaller than the magnitudes of $\mu_{C,0}$ and $\mu_{C,1}$ for network B with logistic dynamics. It can indeed be clearly seen that the shifts of $\tilde{P}_{0}$, $\tilde{P}_{1+}$ and $\tilde{P}_{1-}$ from $P_{0}$, $P_{1+}$ and $P_{1-}$ in Fig. 6 are smaller than the shifts of $\tilde{P}_{0}$ and $\tilde{P}_{1}$ from $P_{0}$ and $P_{1}$ in Fig. 3. Moreover, $\tilde{P}_{0}$, $\tilde{P}_{1+}$ and $\tilde{P}_{1-}$ are again only slightly broadened as compared with $P_{0}$, $P_{1+}$ and $P_{1-}$ because the hidden nodes are randomly chosen. Thus for network C, the links among the measured nodes can also be accurately reconstructed from clustering of $(\Sigma_{m}^{-1})_{ij}$ with $i\neq j$ (see Table 1).

For the scale-free network D, as the link density $\rho$ is very small, most of the measured nodes are not linked via a path of hidden nodes and thus most $C_{ij}$’s vanish. But as the nodes have a power-law degree distribution, both the strength of the hidden nodes and the number of paths connecting measured nodes via hidden nodes could have a large variation leading to a large variation in the magnitude of $C_{ij}$’s. This implies a large $\sigma_{C,0}/\mu_{C,0}$ and $\sigma_{C,1}/\mu_{C,1}$ as compared to the case of random network B. In particular, this results in a larger distortion from $P_{0}$ to ${\tilde{P}}_{0}$ as seen in Fig. 7. The effect is more evident for $P_{0}$ because there are far more unconnected measured nodes than connected measured nodes due to the small $\rho$. Nonetheless, ${\tilde{P}}_{0}$ and ${\tilde{P}}_{1}$ remain distinguishable and the error rates of the reconstruction of the links among the measured nodes remain low (see Table 1).

In network E, every one of the $n_{h}=30$ hidden nodes is only commonly linked to randomly selected pairs of unconnected measured nodes. We first randomly choose a pair of unconnected measured nodes and link all the $n_{h}$ hidden nodes to both of them. Then we link the hidden nodes to a second pair of unconnected measured nodes with the restriction that no hidden nodes are commonly linked to connected measured nodes that might exist among the first and second pairs of unconnected nodes. Therefore, the number of hidden nodes commonly linked to the second pair can be less than $n_{h}$. We repeat the process for all the remaining pairs of unconnected measured nodes. In this way, the number of hidden nodes commonly linked to a given pair of measured nodes $i$ and $j$ or the number of nonzero terms in the first sum in Eq. (40) is identically zero for $i$ and $j$ that are connected, and varies among $i$ and $j$ that are unconnected. This preferential connection of the hidden nodes to unconnected measured nodes gives rise to $\mu_{C,0}>\mu_{C,1}$ and $\sigma_{C,0}>\sigma_{C,1}$. Thus the distortion of $P_{0}$ is large leading to a larger overlap of $\tilde{P}_{0}$ and $\tilde{P}_{1}$ as shown in Fig. 8. As expected, the error rates of the reconstruction results are larger with FP/$N_{L}\approx 30$% (see Table 1). However, we note that even with this error rate, the specificity is above 90%. Moreover, the error rate FN/$N_{L}$ remains less than 1% and thus more than 99% of the actual links among the measured nodes are correctly reconstructed.

VI Conclusions

We have addressed the interesting question of how hidden nodes affect reconstruction of bidirectional networks. By using a model class of bidirectional networks with nonlinear dynamics and diffusive-like coupling and subjected to a Gaussian white noise, as described by Eq. (2), we have derived analytical results, Eqs. (19) and Eq. (22), that allow us to answer this question precisely. Hidden nodes affect the reconstruction results by introducing corrections $C_{ij}$. These corrections $C_{ij}$ are nonzero when the measured nodes $i$ and $j$ are connected via a path of hidden nodes as depicted in Fig. 1. Our estimate of $C_{ij}$, as shown in Eq. (40), shows that three factors determine $C_{ij}$’s, namely the number of paths of hidden nodes connecting the two measured nodes $i$ and $j$, the coupling strength of the links and the strength of the hidden nodes in these paths. Interestingly, hidden nodes with larger strength give rise to smaller corrections when the other two factors remain the same. When the hidden nodes are not preferentially linked to the measured nodes in any manner, these three factors would not differ much between or among the two groups of connected and unconnected measured nodes and, as a result, the hidden nodes would have little effects on the reconstruction of the links among the measured nodes. This is true even when the hidden nodes outnumber the measured nodes. In the event that the hidden nodes are preferentially linked to the measured nodes such that one or more of the above three factors vary significantly either between or among the two groups, the accuracy of the reconstruction results would deteriorate. Yet useful information can still be uncovered. We have verified our theoretical results and their implications using numerical simulations and our numerical results indicate the applicability of our results and analytical understanding beyond the model class of bidirectional networks studied.

Hence our work shows that the method based on the inverse of covariance is useful for reconstructing bidirectional networks even when there are hidden nodes. Most networks of interest in real-world problems are directed. It is highly challenging to derive analytical results for the effects of hidden nodes on the reconstruction of general directed networks. It would thus be interesting to investigate whether and how the present results and understanding for bidirectional networks can be extended to general directed networks.

Acknowledgements.

The work of ESCC and PHT has been supported by the Hong Kong Research Grants Council under grant no. CUHK 14304017.

Appendix A Derivation of Eq. (25)

Denote the eigenvalues of $\bm{M}\equiv\bm{L}_{h}+a{\bm{I}}_{n_{h}}$ by $\lambda_{k}$, $k=1,\ldots,n_{h}$. By Cayley-Hamilton theorem, $\bm{M}$ satisfies its own characteristic equation. Therefore

[TABLE]

where ${\bm{M}}^{0}=\bm{I}_{n_{h}}$ and $e_{m}(\lambda_{1},\ldots,\lambda_{k})$, $1\leq m\leq n_{h}$, are the elementary symmetric polynomials of $\lambda_{k}$’s. For example, $e_{1}=\sum_{k=1}^{n_{h}}\lambda_{k}$ and $e_{n_{h}}=\prod_{k=1}^{n_{h}}\lambda_{k}$. Multiplying ${\bm{M}}^{-1}$ to Eq. (50) and rearranging terms, we express ${\bm{M}}^{-1}$ as a finite polynomial of $\bm{M}$:

[TABLE]

Using $\bm{M}\equiv\bm{S}-\bm{W}$ with $\bm{S}$ and $\bm{W}$ defined in Eqs. (23) and (24), we obtain the elements of ${\bm{M}}^{m}$ in terms of the elements of $\bm{S}$ and $\bm{W}$. For $m=1$ and $2$:

[TABLE]

and for $3\leq m\leq n_{h}-1$:

[TABLE]

Here $f^{(m)}_{0}=(s_{\mu+n}+a)^{m}$ and $f^{(m)}_{m}=(-1)^{m}$ for $m=1,\ldots,n_{h}-1$. For $m\geq 2$, $f^{(m)}_{k}$, $1\leq k\leq m-1$, generally depends on $\hat{s}_{\mu+n}\equiv{s}_{\mu+n}+a$, $\hat{s}_{\nu+n}$ and $\hat{s}_{\alpha_{i}+n}$, $i=1,\ldots,k-1$. Explicit results for $f^{(m)}_{k}$, $1\leq k\leq m-1$, for $m=2$ and 3 are

[TABLE]

Substituting Eqs. (52)-(54) into Eq. (51), we obtain

[TABLE]

where

[TABLE]

depends generally on $\hat{s}_{\mu+n}$, $\hat{s}_{\nu+n}$, $\hat{s}_{\alpha_{i}+n}$, $i=1,\ldots,k-1$ and $\lambda_{k}$. Putting Eq. (58) into Eq. (20), we thus obtain

[TABLE]

which is just Eq. (25).