Network reconstruction from noise-driving dynamic data in complex networks

Abstract

Most complex social, biological, and technological systems can be described by dynamic networks. Reconstructing the structure of complex networks from measurable data of some or all nodes is a challenge in many branches of science. External influences are always present and act as noises to the networks of interest, and various difficulties extensively appear in the reconstruction of real-world networks: such as complexity of network structures; strong nonlinearity of network dynamics; diverse and unknown impacts from the interiors of nodes and externals of networks, i.e., noises; many hidden nodes of which data are not measurable in networks; and different time delays of interactions. Partial or all the above mentioned difficulties are present in reconstruction of noise-driving dynamic network. Different methods are proposed to solve these difficulties, including variable-variable correlations, velocity-variable correlations, high-order correlation, time-lagged covariance of data measurements taken at different times, and so on. This review shows partial developments of a special topic in this wide field. Moreover, we expect that all the methods in this review can be applied to the reconstruction of many realistic dynamic networks.

Graphical abstract

Keywords

network reconstruction / dynamic networks / complex systems

1 Introduction

An outstanding problem in interdisciplinary science is to study nonlinear and complex systems, including their identification, prediction, and control. The ultimate goal in the study of complex systems is to devise practically implementable strategies to control the collective dynamics [1–3]. However, a great challenge is that the network structure and the nodal dynamics are often unknown, and only limited measured time series (data) are available. To control the system dynamics, it is imperative to map out the system details from data. Reconstructing complex network structure and dynamics from the available data, the inverse problem, has thus become one of the central issues in contemporary network science and engineering [4–8]. There are broad applications of the solutions of the network reconstruction problem, due to the ubiquity of complex interacting patterns arising from many systems in a variety of disciplines [9–17].

Many important methodological directions are systematically proposed. Compressive sensing is first applied to predict the catastrophes in nonlinear dynamical systems by Wang et al [18]. Then, it is widely used to detect hidden nodes in complex networks [19], reconstruct propagation networks with natural diversity and identify hidden sources [6], and infer network structures based on evolutionary-game data [20]. Machine learning-based approaches are also used to the inference of dynamics in complex systems [21–24]. The approach of investigation of a dynamical system in terms of its response to external drivings is particularly useful for studying nonlinear systems with complex interactions [25]. Since data-based reconstruction of complex networks is a broad field, we focus on reconstruction of network subjected to fast-varying white noise. More detailed recent review [7] covers many other aspects of reconstruction of nonlinear dynamics in the network.

The rich interplay between nonlinear dynamics and stochastic fluctuations caused by noises of network nodes, often coming from microscopic level and reasonably assumed as white noises, seriously affects the performance of network reconstruction. A number of methods are proposed to deal with this problem [19, 26–41]. Under the condition that the influence of noises on the evolution of infinitesimal tangent vectors in the phase space of the underlying networked dynamical system is dominant, it can be argued that the dynamical correlation matrix, which can be computed readily from the available nodal time series, approximates the network adjacency matrix, fully unveiling the network topology [7].

This review focuses on reconstruction of network subjected to white noise with high-frequency measurable data and demonstrates how to extract diverse useful information by computing rich varieties of correlations to solve various difficulties in network reconstruction. The first two chapters consider instant interactions between network nodes where signal propagation between network nodes is much faster than that of node dynamics. In Section 2, we start from the simplest noise-driving linear networks and achieve reconstruction solution in simple matrix form by analyzing correlation matrices between data of nodes and sequences of noisy outputs. In Section 3, we extend the above correlation computation method to general nonlinear networks, and obtain compact and unified matrix form with various expansion strategies of bases which can be applied to noise-driving nonlinear networks. Sections 4 and 5 consider interactions with time-delays representing important practical networks of which time for signal propagation between network nodes is comparable with that of node evolution. Section 4 studies the basic effects of definite signal propagation time, and obtains results of intensities and total time delay of both direct and indirect interactions with data of pairs of nodes only. ln Section 5 we reconstruct network with hidden nodes, a challenging and difficult task. Different methods are suggested to overcome these difficulties [19, 42–52]. Here, we apply the results of long-distance interaction reconstruction and propose a novel concept of solvable motif subnetworks with hidden nodes, definitely identifying hidden nodes and the related hidden links.

It is expected that all the above methods, strategies, and results could be effectively applied to reconstruction of practically important networks.

2 Solving the inverse problem of noise-driven dynamic networks

A large class of dynamic networks can be most generally represented by the following coupled autonomous ordinary differential equations (ODEs) driven by noise

(1)\dot{\mathit{x}}(t)=\mathit{f}(\mathit{x}(t),\eta (t)),

with variables \mathit{x}(t)=(x_1(t),x_2(t),\cdots ,x_N(t)), noise \eta (t)=({\eta }_1(t),{\eta }_2(t),\cdots ,{\eta }_M(t)) and dynamic fields \mathit{f}(\mathit{x}(t),\eta (t))=(f_1(\mathit{x}(t),\eta (t)),\cdots ,f_N(\mathit{x}(t),\eta (t))). In this paper, we adopt white noise approximation for very short noise correlation time,

(2)⟨{\eta }_i(t)⟩=0,⟨{\eta }_i(t){\eta }_j(t^{\prime })⟩=D_{ij}\delta (t-t^{\prime }),i,j=1,2,\cdots ,M.

Around a stable and stationary phase point of deterministic system in Eq. (1) and for small noise approximation, Eq. (1) can be linearized to autonomous coupled equations

(3)\dot{\mathit{y}}(t)=\hat{\mathit{A}}(\mathit{x})\mathit{y}(t)+\mathrm{\Gamma }(t),

\mathit{y}(t)=(y_1(t),y_2(t),\cdots ,y_N(t)),\mathit{y}(t)=\mathit{x}(t)-\mathit{x},

A_{ij}(\mathit{x})={\displaystyle \frac{\mathrm{\partial }{f}_i}{\mathrm{\partial }{x}_j}}{\mid }_{\mathit{x}(t)=\mathit{x},\eta (t)=0},\mathrm{\Gamma }(t)=({\mathrm{\Gamma }}_1(t),\cdots ,{\mathrm{\Gamma }}_N(t)),

(4)⟨\mathrm{\Gamma }(t)⟩=0,⟨{\mathrm{\Gamma }}_i(t){\mathrm{\Gamma }}_j(t^{\prime })⟩=Q_{ij}\delta (t-t^{\prime }).

Now, we consider how to infer network structure \hat{\mathit{A}} merely from measurable data \mathit{y}(t). Suppose we have L pairs of variable data (\mathit{x}(t_q),\mathit{x}(t_q+\mathrm{\Delta }{t}_q),q=1,2,\cdots ,L) in a small phase space region, with all \mathrm{\Delta }{t}_q, \Vert \mathit{x}(t_{q^{\prime }})-\mathit{x}(t_q)\Vert, \Vert \mathit{x}(t_q+\mathrm{\Delta }{t}_q)-\mathit{x}(t_q)\Vert \ll 1 for all q, q^{\prime }. Thus, \hat{\mathit{A}}(\mathit{x}) in Eq. (3) can be regarded as a constant matrix in our inference computation. From the available data we can extract full information of \mathit{y}(t_q)=\mathit{x}(t_q)-\mathit{x} and \mathit{x}=\frac{1}{L}\sum _{q=1}^L\mathit{x}(t_q). The velocity term \dot{\mathit{y}}(t_q) in Eq. (3) can be measured as

(5){\dot{y}}_i(t_q)={\displaystyle \frac{y_i(t_q+\mathrm{\Delta }{t}_q)-y_i(t_q)}{\mathrm{\Delta }{t}_q}}.

Multiplying {\mathit{y}}^{\mathrm{T}}(t) to both sides of Eq. (3) and averaging all the terms in the equation for the L data samples, we can reach an identity of correlation matrices as

(6)\hat{\mathit{B}}=\hat{\mathit{A}}\hat{\mathit{C}}+⟨\mathrm{\Gamma }(t){\mathit{y}}^{\mathrm{T}}(t)⟩,

where \hat{\mathit{C}}=⟨\mathit{y}{\mathit{y}}^{\mathrm{T}}⟩ and \hat{\mathit{B}}=⟨\dot{\mathit{y}}{\mathit{y}}^{\mathrm{T}}⟩ are the variable-variable and velocity-variable correlation matrices, respectively,

(7)C_{ij}={\displaystyle \frac{1}{L}}{\displaystyle \sum _{q=1}^L}{y}_i(t_q)y_j(t_q),B_{ij}={\displaystyle \frac{1}{L}}{\displaystyle \sum _{q=1}^L}{\dot{y}}_i(t_q)y_j(t_q),

and

Since noise \mathrm{\Gamma }(t_q^{\prime }) contributing to the velocity \dot{\mathit{y}}(t_q) in Eq. (5) (and also contributing to matrix \hat{\mathit{B}}) is taken during (t_q,t_q+\mathrm{\Delta }{t}_q) and \mathrm{\Delta }{t}_q is much longer than the noise correlation time, \mathrm{\Gamma }(t_q^{\prime }) is thus not correlated with {\mathit{y}}^{\mathrm{T}}(t_q) and the second term in Eq. (6) is vanishing. We have thus identified noise-variable decorrelation

(8)⟨\mathrm{\Gamma }(t){\mathit{y}}^{\mathrm{T}}(t)⟩=0

and reach a simple algorithm \hat{\mathit{B}}=\hat{\mathit{A}}\hat{\mathit{C}}, leading to the solution

(9)\hat{\mathit{A}}=\hat{\mathit{B}}{\hat{\mathit{C}}}^{-1}.

We can also multiply the both sides of Eq. (3) by {\mathit{y}}^{+\mathrm{T}}(t)={\mathit{y}}^{\mathrm{T}}(t+\mathrm{\Delta }t), and compute the corresponding correlations, and arrive at

(10){\hat{\mathit{B}}}^{+}=\hat{\mathit{A}}\hat{\mathit{C}}+⟨\mathrm{\Gamma }(t){\mathit{y}}^{+\mathrm{T}}(t)⟩,

\begin{array}{rl}{\hat{\mathit{B}}}^{+}=& ⟨\dot{\mathit{y}}(t){\mathit{y}}^{+\mathrm{T}}(t)⟩,⟨\mathrm{\Gamma }(t){\mathit{y}}^{+\mathrm{T}}(t)⟩\\ \approx & \frac{1}{L}\sum _{q=1}^L\mathrm{\Gamma }(t_q^{\prime }){\mathit{y}}^{+\mathrm{T}}(t_q),\end{array}

where {\mathit{y}}^{+}(t_q)=\mathit{y}(t_q+\mathrm{\Delta }{t}_q) is taken at time instant t_q+\mathrm{\Delta }{t}_q, and thus is definitely correlated with \mathrm{\Gamma }(t_q^{\prime }) which is taken during (t_q,t_q+\mathrm{\Delta }{t}_q). This correlation can be calculated in the continuous limit by solving Eq. (3) as \mathit{y}(t)=\mathit{y}(t_0)+{\mathrm{e}}^{\hat{\mathit{A}}t}{\int }_{t_0}^t{\mathrm{e}}^{-\hat{\mathit{A}}{t}^{\prime }}\mathrm{\Gamma }(t^{\prime })\mathrm{d}{t}^{\prime }, Due to the delta function of noise correlation, the integrals are singular at time t, with {\mathit{y}}^{+}(t) (\mathit{y}(t)) including (not including) the noise impact at time t. And thus we have

(11a)\begin{array}{c}⟨\mathrm{\Gamma }(t){\mathit{y}}^{\mathrm{T}}(t)⟩=0,\end{array}

(11b)\begin{array}{c}⟨\mathrm{\Gamma }(t){\mathit{y}}^{+\mathrm{T}}(t)⟩=\hat{\mathit{Q}}.\end{array}

Substituting Eq. (11a) to Eq. (6) we obtain Eq. (9), and substituting Eq. (11b) to Eq. (10) we have {\hat{\mathit{B}}}^{+}=\hat{\mathit{A}}\hat{\mathit{C}}+\hat{\mathit{Q}}, leading to,

(12)\begin{array}{c}\hat{\mathit{A}}=({\hat{\mathit{B}}}^{+}-\hat{\mathit{Q}}){\hat{\mathit{C}}}^{-1},\end{array}

and to solve further noise correlation

(13)\begin{array}{c}\hat{\mathit{Q}}={\hat{\mathit{B}}}^{+}-\hat{\mathit{B}}.\end{array}

Due to the identity

\begin{array}{rl}{B}_{ij}^{+}=& \frac{1}{L}\sum _{q=1}^L{\dot{y}}_i(t_q)y_j^{+}(t_q)=\frac{1}{L}\sum _{q=1}^L\frac{y_i^{+}(t_q)-y_i(t_q)}{\mathrm{\Delta }{t}_q}{y}_j^{+}(t_q)\\ =& \frac{1}{L}\sum _{q=1}^L\frac{y_i^{+}(t_q)y_j^{+}(t_q)-y_i(t_q)y_j^{+}(t_q)}{\mathrm{\Delta }{t}_q}\\ \approx & \frac{1}{L}\sum _{q=1}^L\frac{y_i(t_q)y_j(t_q)-y_i(t_q)y_j^{+}(t_q)}{\mathrm{\Delta }{t}_q}\\ =& \frac{1}{L}\sum _{q=1}^L\frac{y_j(t_q)-y_j^{+}(t_q)}{\mathrm{\Delta }{t}_q}{y}_i(t_q)=-\frac{1}{L}\sum _{q=1}^L{\dot{y}}_j(t_q)y_i(t_q)\\ =& -B_{ji},\end{array}

we have

{\hat{\mathit{B}}}^{+}=-{\hat{\mathit{B}}}^{\mathrm{T}},

and

(14)\begin{array}{c}\hat{\mathit{Q}}={\hat{\mathit{B}}}^{+}-\hat{\mathit{B}}=-(\hat{\mathit{B}}+{\hat{\mathit{B}}}^{\mathrm{T}})=-2{\hat{\mathit{B}}}_s,\end{array}

where {\hat{\mathit{B}}}_s and {\hat{\mathit{B}}}^{\mathrm{T}} are the symmetric part and the transposition of \hat{\mathit{B}}, respectively.

Now, a novel double correlation matrices and noise-decorrelation (DCMND) method is proposed to generally solve the inverse problem of Eq. (3) and also explicitly depict noise statistical correlation matrix \hat{\mathit{Q}}, by using the simple and unified algorithms (9) and (14) [30].

Actually, Eq. (3) can be generally derived for any phase space point where output data are available and \hat{\mathit{A}} is thus \mathit{x} dependent. Around a stable steady state of noise-free system, we can linearize Eq. (1) directly around the fixed point {\mathit{x}}_0 and set \mathit{x}={\mathit{x}}_0 and \mathit{y}(t)=\mathit{x}(t)-{\mathit{x}}_0 for computing (5) and (7). Now, our task is to reveal both Jacobian matrix \hat{\mathit{A}} in (3) and noise correlation \hat{\mathit{Q}} in (4) from measurable variable data \mathit{y}(t).

In Fig. 1 we consider two examples in Eq. (3) for numerical simulations. The network size in Figs. 1(a) and (b) (Case 1) is N=100. Among {10}^3 links there are 500 active links A_{ij}(a)=1 and 500 repressive ones A_{ij}(r)=-1 arbitrarily chosen, A_{ij}(0)=0 for all else off-diagonal matrix elements, and the diagonal terms are set to A_{ii}=-3 for keeping the network evolution bounded. Moreover, we take Q_{ij}={\sigma }_i{\delta }_{ij},{\sigma }_i=0.01. Running system (3) from y_i(t=0)=0 we produce L=5\times {10}^5 sets of data. Assume we know nothing about network structure \hat{\mathit{A}} and noise statistics \hat{\mathit{Q}} but only the data sequences of \mathit{y}(t), among which a projective trajectory is plotted in a 2D (y_1(t), y_2(t)) phase plane in Fig. 1(a), which seems fully disordered. All the elements of matrices \hat{\mathit{B}} and \hat{\mathit{C}} can be computed from the data \mathit{y}(t) with Eqs. (5) and (7), and then matrix \hat{\mathit{A}} can be retrieved by Eq. (9). The results are presented in Fig. 1(b) where interactions depicted by \hat{\mathit{B}}{\hat{\mathit{C}}}^{-1} agree well with the actual interactions \hat{\mathit{A}}. In Figs. 1(c) and (d) (Case 2), we do the same as Figs. 1(a) and (b) with the noise statistics changed to Q_{ij}={\sigma }_i{\delta }_{ij} with {\sigma }_i\in (0.005,0.015), and A_{ij} set to A_{ij}(a)\in (1.0,2.0), A_{ij}(r)\in (-2.0,-1.0) (all randomly chosen with uniform distributions in their ranges) and the diagonal terms are set to A_{ii}=-5. Though the trajectory behaviors of Figs. 1(a) and (c) deviate from each other substantially due to different noise correlations \hat{\mathit{Q}} and Jacobian matrix \hat{\mathit{A}}, and these different data sets can surely yield considerably different matrices \hat{\mathit{B}} and \hat{\mathit{C}}, it is remarkable that in Figs. 1(b) and (d) the DCMND method can correctly deduce both interaction matrices \hat{\mathit{A}} by applying the same algorithm \hat{\mathit{A}}=\hat{\mathit{B}}{\hat{\mathit{C}}}^{-1}. For confirming the conclusions of Eq. (14), we plot -2{\hat{\mathit{B}}}_s computed from the variable data against actual noise statistics \hat{\mathit{Q}} in Fig. 1(e) for the two data sets of Figs. 1(a) and (c). All these dots locate very closely around the diagonal line, convincingly justifying the prediction of Eq. (14). The whole theoretical derivations and results of statistical computations are robust against different noise intensities. In Figs. 1(f) and (g), we do exactly the same as Figs. 1(c) and (d), respectively, with incomparably stronger noise used (100 times of {\sigma }_i in (c)(d)). Matrix \hat{\mathit{A}} can be retrieved equally well.

The results of Eqs. (9)(14) are exact in the limits of L\to \mathrm{\infty }, \mathrm{\Delta }{t}_q\to 0 and \mathrm{\Delta }{t}_q\gg \tau with \tau being the correlation time of noise (now \tau \equiv 0 with white noise approximation). In Fig. 1(h) we use the systems of Figs. 1(a) and (c) to numerically compute the standard deviation (SD) of inferred values of interactions defined as

(15)\begin{array}{c}SD(L)=\sqrt{{\displaystyle \frac{{{\displaystyle \sum }}_{i=1}^N{{\displaystyle \sum }}_{j=1}^N([B{C}^{-1}{]}_{ij}-A_{ij}{)}^2}{N^2}}},\end{array}

where summations of i and j run over all matrix elements and L is the total number of computational samples. It is clear that SD\propto \frac{1}{\sqrt{L}}, agreeing with the conclusion of exact inference solutions for L\to \mathrm{\infty }. The theoretical conclusion that algorithms (9) and (14) can reveal both interaction structure \hat{\mathit{A}} and noise statistics \hat{\mathit{Q}} are remarkably verified by numerical simulations based on which we can reconstruct the stochastic dynamic networks of Eq. (3) from their variable outputs only. And this capacity of the DCMND method is unique in all known inference methods. We observe, of course, certain errors in all Figs. 1(b), (d), (e) and (g) with small (large) while finite \mathrm{\Delta }{t}_q (L). Decreasing \mathrm{\Delta }{t}_q and increasing L by keeping \tau =0 can surely increase the precision of depiction.

3 Reconstruction of noise-driven nonlinear networks from node outputs by using correlations of nonlinear functions

Let us consider a general noise-driven nonlinear network of Eq. (1) as

(16)\begin{array}{c}\begin{array}{rl}& \dot{\mathit{x}}(t)=\mathit{f}(\mathit{x}(t))+\mathbf{\Gamma }(t),\\ & \mathit{x}=(x_1,x_2,\cdots ,x_N{)}^{\mathrm{T}},\\ & \mathit{f}=(f_1,f_2,\cdots ,f_N{)}^{\mathrm{T}},\\ & \mathbf{\Gamma }=({\mathrm{\Gamma }}_1,{\mathrm{\Gamma }}_2,\cdots ,{\mathrm{\Gamma }}_N{)}^{\mathrm{T}},\end{array}\end{array}

where T represents the operation of transposition. Dynamical noises {\mathrm{\Gamma }}_i(t), i=1,2,\cdots ,N, represent the impacts from microscopic world, and they are expected to have very short correlation time {\tau }_d\ll 1, much smaller than the characteristic time of deterministic dynamics assumed to be of order 1. Then noises are approximated as white ones,

(17)\begin{array}{c}⟨{\mathrm{\Gamma }}_i(t)⟩=0,⟨{\mathrm{\Gamma }}_i(t){\mathrm{\Gamma }}_j(t^{\mathrm{\prime }})⟩=Q_{ij}(\mathit{x})\delta (t-t^{\mathrm{\prime }}),\end{array}

with i,j=1,2,\cdots ,N. It is emphasized that multiplicative noises have been seldom considered so far in the study of network reconstruction, though this type of noises exist extensively in some practical circumstances. In Eq. (16) we have all measurable data in our hand, namely, we measure

(18)\begin{array}{c}\begin{array}{rl}& \mathit{x}(t_1),\mathit{x}(t_2),\cdots ,\mathit{x}(t_k),\cdots ,\mathit{x}(t_L),\\ & \mathrm{\Delta }t=t_{k+1}-t_k\ll 1;k=1,2,\cdots ,L;L\gg 1.\end{array}\end{array}

With \mathrm{\Delta }t\ll 1 we can compute velocities of \mathit{x} in Eq. (16) and with L\gg 1 we have sufficiently large samples to perform statistical analysis. These conditions are not always available, but they may be available in many important practical experiments, or can be realized on purpose in many realistic measurements.

In Eq. (16) all functions f_i, i=1,2,\cdots ,N, are unknown. The noise statistic matrix \hat{\mathit{Q}}(\mathit{x})=(Q_{ij}(\mathit{x})) is also unknown. Only the output variables (18) are available for analysis. The task is to specify dynamic functions f_i, the interactions between different nodes, and noise statistics Q_{ij}, i,j=1,\cdots ,N.

First we assume f_is can be generally expanded by a certain basis set [34, 42]

(19)\begin{array}{c}{f}_i(\mathit{x})={\displaystyle \sum _{\mu =1}^{\mathrm{\infty }}{A}_{i,\mu }{Y}_{i,\mu }(\mathit{x}),i=1,2,\cdots ,N,}\end{array}

where all constant coefficients A_{i,\mu }, \mu =1,2,\cdots ,M_i,\cdots ,\mathrm{\infty } are unknown, while all functions Y_{i,\mu }(\mathit{x}) called as bases are known. For treating nonlinearities in \mathit{f}(\mathit{x}), the chosen basis set should be nonlinear and complete for expanding field \mathit{f}(\mathit{x}). In Eq. (16) we give a freedom to use different basis sets for expanding different field \mathit{f}(\mathit{x}). It seems that the expansions of Eq. (19) should include infinite terms for arbitrary functions \mathit{f}(\mathit{x}). One has to truncate the expansion to finite terms. There is a systematical and self-consistent method, described in the following part to make such truncation. At present, we just assume a truncation at M_i for f_i(\mathit{x}) expansion. Then Eq. (19) can be simplified as

(20)\begin{array}{c}\begin{array}{rl}& f_i(\mathit{x}(t))={\hat{\mathit{A}}}_i^{\mathrm{T}}{\mathit{Y}}_i(t),\\ & {\hat{\mathit{A}}}_i^{\mathrm{T}}=(A_{i,1},A_{i,2},\cdots ,A_{i,M_i}),\\ & {\mathit{Y}}_i=(Y_{i,1}(t),Y_{i,2}(t),\cdots ,Y_{i,M_i}(t){)}^{\mathrm{T}}.\end{array}\end{array}

Without noise {\mathrm{\Gamma }}_i(t), all the unknown coefficients can be solved by algebraic equations if sufficient data are accumulated [18, 53]. With strong noises the inference computations are much more difficult. Here for network reconstruction we use a method to compute two-time correlations to filter noise effect, together with using correlations of the chosen nonlinear basis sets to reconstruct interactions, nonlinearities of the network nodes and multiplicative noises [30, 33, 34].

For arbitrary node i in the network, multiplying Eq. (16) from the right side by a functional vector {\mathit{Z}}_i^{\mathrm{T}}(\mathit{x}(t-\tau ))

(21)\begin{array}{c}{\mathit{Z}}_i(\mathit{x})=(Z_{i,1}(\mathit{x}),Z_{i,2}(\mathit{x}),\cdots ,Z_{i,M_i}(\mathit{x}){)}^{\mathrm{T}},\end{array}

and computing all related correlations, we obtain a linear matrix algebraic equation

(22)\begin{array}{c}{\hat{\mathit{B}}}_i^{\mathrm{T}}(-\tau )={\hat{\mathit{A}}}_i^{\mathrm{T}}{\hat{\mathit{C}}}_i(-\tau )+⟨{\mathrm{\Gamma }}_i(t){\mathit{Z}}_i^{\mathrm{T}}(t-\tau )⟩\end{array}

with

(23)\begin{array}{c}\begin{array}{rl}& {\hat{\mathit{B}}}_i(-\tau )=(B_{i,1}(-\tau ),B_{i,2}(-\tau ),\cdots ,B_{i,M_i}(-\tau ){)}^{\mathrm{T}},\\ & B_{i,\mu }(-\tau )=⟨{\dot{x}}_i(t)Z_{i,\mu }(t-\tau )⟩=\frac{1}{L-p}\sum _{k=p+1}^L{\dot{x}}_i(t_k)Z_{i,\mu }(t_{k-p}),\\ & {\dot{x}}_i(t_k)=\frac{x_i(t_{k+1})-x_i(t_k)}{\mathrm{\Delta }t},\tau =p\mathrm{\Delta }t,\\ & {\hat{\mathit{C}}}_i(-\tau )=(C_{i,\mu \nu })=(\frac{1}{L-p}\sum _{k=p+1}^L{Y}_{i,\mu }(t_k)Z_{i,\nu }(t_{k-p})).\end{array}\end{array}

{\mathit{Z}}_1(\mathit{x}(t)),{\mathit{Z}}_2(\mathit{x}(t)),\cdots ,{\mathit{Z}}_N(\mathit{x}(t)) are called correlators. They can be arbitrarily chosen for computing the correlation matrix {\hat{\mathit{C}}}_i, provided that the entries of each correlator vector Z_i(x) are linearly independent. This condition guarantees that the correlation matrix {\hat{\mathit{C}}}_i is nonsingular and thus invertible. In Eq. (23) we should have 0\approx {\tau }_d\ll \tau \ll 1 with \tau being larger than the correlation time of dynamical noises {\tau }_d, and much smaller than the characteristic times of deterministic network dynamics, previously assumed to be of order 1. Therefore, noises and correlators must be decorrelated,

(24)\begin{array}{c}⟨{\mathrm{\Gamma }}_i(t){\mathit{Z}}_i(t-\tau {)}^{\mathrm{T}}⟩\approx 0,\mathrm{f}\mathrm{o}\mathrm{r}\tau >{\tau }_{\mathrm{d}},\end{array}

since the fast-varying noises of Eq. (17) have no correlation with any variable data of earlier times, disregarding any forms of multiplicative noises Q_{ij}(\mathit{x}). Now with the noise-decorrelation of Eq. (24), Eq. (22) can be reduced to

(25)\begin{array}{c}{\hat{\mathit{B}}}_i^{\mathrm{T}}(-\tau )={\hat{\mathit{A}}}_i^{\mathrm{T}}{\hat{\mathit{C}}}_i,\end{array}

leading to

(26)\begin{array}{c}{\hat{\mathit{A}}}_i^{\mathrm{T}}={\hat{\mathit{B}}}_i^{\mathrm{T}}(-\tau ){\hat{\mathit{C}}}_i^{-1},\end{array}

where {\hat{\mathit{B}}}_i(-\tau ), {\hat{\mathit{C}}}_i and {\hat{\mathit{A}}}_i are given in Eqs. (23) and (20), respectively. We delete the notion (-\tau ) in {\hat{\mathit{C}}}_i because \tau \ll 1 does not considerably change the values of C_{i,\mu \nu }. All elements of vector {\hat{\mathit{B}}}_i(-\tau ) and matrix {\hat{\mathit{C}}}_i can be computed with known output variables \mathit{x}(t), and thus all the unknown linear and nonlinear coefficients in Eq. (16) can be inferred with a simple and unified matrix equation (26), though the noise matrix \hat{\mathit{Q}} in Eq. (17) is unknown.

We consider an example of a three-node nonlinear network, the Lorenz system, one of the most famous models in chaos study,

(27)\begin{array}{c}\begin{array}{rl}& \dot{x}=f_1(x,y,z)+{\mathrm{\Gamma }}_1(t),f_1=\sigma (y-x),\\ & \dot{y}=f_2(x,y,z)+{\mathrm{\Gamma }}_2(t),f_2=x(\rho -z)-y,\\ & \dot{z}=f_3(x,y,z)+{\mathrm{\Gamma }}_3(t),f_3=xy-\beta z.\end{array}\end{array}

Though network (27) looks very simple and low-dimensional, it serves as one of the most prevalent prototypes for data analysis. We consider noises

(28)\begin{array}{c}\hat{\mathit{Q}}=\left(\begin{array}{ccc}{a}_1+a_2|z|& 0& 0\\ 0& b_1+b_2|x|& 0\\ 0& 0& c_1+c_2|y|\end{array}\right),\end{array}

where we adopt absolute values |x|, |y|, |z| in Q_{ii} functions for Q_{ii} should be positive definite or zero from the definition of Eq. (17). In Eq. (27) both noise and nonlinearity essentially effect and noises with different multiplicative factors can produce considerably different data sets, as shown in Fig. 2 where trajectories for additive noises [Figs. 2(a) and (b)] and multiplicative noises [Figs. 2(c) and (d)] are presented. In Fig. 2, the time sequences and trajectories are obtained by directly integrating the noise-driven Lorenz system in Eq. (27) with the additive and multiplicative noises specified above. We then apply the general algorithm in Eq. (26) to reconstruct the Lorenz system from these data. For the reconstruction calculation, we choose the correlator vector to be the same as the basis vector for the nonlinear field expansion with truncation M. Specifically, without particular prior information, a power-series basis is adopted, and we set

(29)\begin{array}{c}\begin{array}{rl}& {\mathit{Y}}_i={\mathit{Z}}_i=\mathit{Y},i=1,2,3;\\ & \mathit{Y}=(Y_1,Y_2,\cdots ,Y_M{)}^{\mathrm{T}}=(1;x,y,z;x^2,xy,xz,y^2,yz,z^2;\cdots {)}^{\mathrm{T}}.\end{array}\end{array}

Here we introduce an idea of self-consistent checking of the truncation. At first we take M_i^0 bases in Eq. (29), compute all elements of corresponding \widehat{{\mathit{C}}_i}, {\hat{\mathit{B}}}_i and obtain {\hat{\mathit{A}}}_i(M_i^0). Then we further consider some more bases, i.e., M_i^1>M_i^0 bases in Eq. (29), and obtain results of \{{\hat{\mathit{A}}}_i(M_i^1)|M_i^1>M_i^0\}. M_i^0 is concluded as a suitable truncation if the results satisfy two conditions. Condition (i): A_{i,\mu }(M_i^1)\approx A_{i,\mu }(M_i^0) for all coefficients obtained with M_i^0 truncation; Condition (ii): A_{i,\mu }(M_i^1)\approx 0 for all coefficients of bases not included by M_i^0 truncation. If any of the above two conditions is not satisfied we should go on to include more bases in Eq. (29), M_i^2>M_i^1, M_i^3>M_i^2 and so on, until the two conditions are satisfied at M_i^k truncation. We then conclude M_i^{k-1} is the suitable truncation. Now we use the above method of truncations, successively from low orders to high orders of Eq. (29).

A self-consistent justification on correct truncation for Eq. (27) is illuminated in Fig. 3. In Fig. 3(a), it is apparent that M_i=4 (bases: 1,x,y,z) is not a proper truncation, since the reconstruction of M_i=4 is considerably different from that of M_i=10, where all bases of powers of the second order are taken into account. In Fig. 3(b) we compare the reconstruction results truncated by the second order (M_i=10) with those of the third order (M=20), and the above self-consistent checking conditions (i) and (ii) are both fulfilled. Therefore we can conclude that it is enough to truncate the expansion at M_i=10 with properly chosen bases in Eq. (29). In Fig. 3(c) we compare the reconstructed results for M_i=10 with the actual A_{i,\mu }, and the nonlinear and interactive structures of Eq. (27) are satisfactorily recovered indeed.

In this review, we consider the Lorentz equation with small number of nodes and finite power expansion terms. The method can be easily applied to reconstruction of nonlinear networks with large number of nodes and with infinitely many expansion terms [34]. There exist many methods to treat network reconstruction problems. Similar approach of self-consistent truncation checking can be also used for specifying multiplicative white noise matrix, for more details see Ref. [34]. For example, the compressive sensing method can exclude many zero coefficients from the inference computations and effectively reduce parameter number and computational difficulties [54, 55].

4 Reconstructing direct and distant interactions of multiple paths between pairs of perceptible nodes in dark networks

We consider again a complex network of N nodes, which is described by a set of differential equations where all functions f_i,i=1,2,\cdots ,N are prior unknown,

(30a)\begin{array}{c}{\dot{x}}_i(t)=f_i\left({\mathit{x}}_i(t)\right)+{\mathrm{\Gamma }}_i(t),i=1,2,\cdots ,N,\end{array}

where

{\mathit{x}}_i(t)=\left(\begin{array}{c}{x}_1(t-{\tau }_{i1}),\cdots ,x_i(t),\cdots ,x_N(t-{\tau }_{iN})\end{array}\right).

f_i is a general form of the dynamical function of node i, including the settling local dynamics and the interactions from other nodes in the network (\frac{\mathrm{\partial }{f}_i}{\mathrm{\partial }{x}_j}=0 if interaction from j is absent). A positive constant {\tau }_{ij} denotes time delay signal propagation from j to i. A dynamical model which considers the transmission time delay is close to practical systems and is capable of elucidating much richer behaviors in real world. We adopt a Gaussian white noise {\mathrm{\Gamma }}_i(t) for inherent noise disturbances satisfying

(30b)\begin{array}{c}⟨{\mathrm{\Gamma }}_i(t)⟩=0,⟨{\mathrm{\Gamma }}_i(t){\mathrm{\Gamma }}_i(t+t^{\prime })⟩=Q_i\delta (t^{\prime }),\end{array}

where Q_i is the noise intensity, a positive constant. The delta function \delta reflects that correlation time of the noise is much shorter than the characteristic time of the system dynamics. Here, we take Q_{ij}=0,i\ne j, because noises coming from microscopic level of different nodes fluctuate independently.

We can measure the variable x_i(t) of node i that changes continuously over time if it is perceptible. Our goal is to reconstruct the interactions between two perceptible nodes based on the measurable data of these two nodes only. The other information, such as, the dynamics of all nodes, the connection structure of entire network, the transmission time delay between adjacent nodes, background noise, etc., are all unknown.

For reconstructing interaction from one node (say, node A) to another node (node B), we inject a varying external driving D_A(t) at node A as

(31)\begin{array}{c}{\dot{x}}_A(t)=f_A\left({\mathit{x}}_A(t)\right)+{\mathrm{\Gamma }}_A(t)+D_A(t),\end{array}

where the requirement for the driving term D_A(t) will be induced below. Assuming that there is a direct connection from A to B, for node B we have

(32)\begin{array}{c}{\dot{x}}_B(t)=f_B\left({\mathit{x}}_B(t)\right)+{\mathrm{\Gamma }}_B(t),\end{array}

where \frac{\mathrm{\partial }{f}_B}{\mathrm{\partial }{x}_A}\ne 0, called B is distance 1 (d=1) from A.

To analyze output response in B, first we take the derivative of two sides of Eq. (32) with respect to t,, which gives

(33)\begin{array}{c}\begin{array}{rl}{\ddot{x}}_B(t)=& \sum _{i=1}^N\frac{\mathrm{\partial }{f}_B\left({\mathit{x}}_B(t)\right)}{\mathrm{\partial }{x}_i(t-{\tau }_{Bi})}\frac{d{x}_i(t-{\tau }_{Bi})}{dt}+{\dot{\mathrm{\Gamma }}}_B(t)\\ =& \sum _{i=1}^N\frac{\mathrm{\partial }{f}_B\left({\mathit{x}}_B(t)\right)}{\mathrm{\partial }{x}_i(t-{\tau }_{Bi})}{f}_i\left({\mathit{x}}_i(t-{\tau }_{Bi})\right)\\ & +\sum _{i=1}^N\frac{\mathrm{\partial }{f}_B\left({\mathit{x}}_B(t)\right)}{\mathrm{\partial }{x}_i(t-{\tau }_{Bi})}{\mathrm{\Gamma }}_i(t-{\tau }_{Bi})\\ & +\frac{\mathrm{\partial }{f}_B\left({\mathit{x}}_B(t)\right)}{\mathrm{\partial }{x}_A(t-{\tau }_{BA})}{D}_A(t-{\tau }_{BA})+{\dot{\mathrm{\Gamma }}}_B(t),\end{array}\end{array}

where {\ddot{x}}_B is the second derivative in time of x_B and {\dot{\mathrm{\Gamma }}}_B denotes an internal noise-induced term which has no correlation with the injected noise signal D_A(t). It can be concluded that there is only a single driving-induced term which builds a direct quantitative relationship with coupling function, provided D_A(t) is independent of system,

(34)\begin{array}{c}{\displaystyle \frac{\mathrm{\partial }{f}_B\left({\mathit{x}}_B(t)\right)}{\mathrm{\partial }{x}_A(t-{\tau }_{BA})}}{D}_A(t-{\tau }_{BA}).\end{array}

We collect the response data of node B with a fixed time interval (\mathrm{\Delta }t\ll 1), denoted as

(35)\begin{array}{c}[x_B(t_1),x_B(t_2),\cdots ,x_B(t_l),\cdots ,x_B(t_L)],\end{array}

where t_{l+1}-t_l=\mathrm{\Delta }t, T=L\mathrm{\Delta }t. Note, the information included in these data not only the L data values but also various differential relationships of the data in sequences, from which we can produce large independent correlation quantities. From data sequence, the first and the second derivative of x_B are given by

{\dot{x}}_B(t_l)\approx x_B^{(1)}(t_l)={\displaystyle \frac{x_B(t_{l+1})-x_B(t_l)}{\mathrm{\Delta }t}},

{\ddot{x}}_B(t_l)\approx x_B^{(2)}(t_l)={\displaystyle \frac{x_B^{(1)}(t_{l+1})-x_B^{(1)}(t_l)}{\mathrm{\Delta }t}},

and for \mu \ge 1,

(36)\begin{array}{c}{\displaystyle \frac{{\mathrm{d}}^{\mu +1}}{\mathrm{d}{t}^{\mu +1}}}{x}_B(t_l)\approx x_B^{(\mu +1)}(t_l)={\displaystyle \frac{x_B^{(\mu )}(t_{l+1})-x_B^{(\mu )}(t_l)}{\mathrm{\Delta }t}}.\end{array}

To detect the effect of the driving component D_A included in the data of x_B(t), we calculate the correlation function between the driving signal and x_B^{(2)} where t^{\prime }\in [0,T]

(37)\begin{array}{c}\begin{array}{rl}{E}_{BA}^{2,1}(t^{\prime })=& \frac{1}{T-t^{\prime }}{\int }_0^{T-t^{\prime }}{D}_A(t)x_B^{(2)}(t+t^{\prime })\mathrm{d}t\\ =& \frac{1}{T-t^{\prime }}{\int }_0^{T-t^{\prime }}\frac{\mathrm{\partial }{f}_B\left({\mathit{x}}_B(t)\right)}{\mathrm{\partial }{x}_A(t-{\tau }_{BA})}{D}_A(t)D_A(t+t^{\prime }-{\tau }_{BA})\mathrm{d}t\\ & +E_{{\mathrm{\Gamma }}_B}^1+E_{f_B}^1.\end{array}\end{array}

This yields an intensity spectrum by varying t^{\prime }, the first term of which is directly related to the autocorrelation property of driving D_A(t). There are also two other terms E_{{\mathrm{\Gamma }}_B}^1 and E_{f_B}^1 which are evidently unrelated to D_A(t).

For the term of E_{{\mathrm{\Gamma }}_B}^1, the emergence of noise is bound to cause a random deviation. From Eqs. (32, 33), it is known that the leading deviation component is given by

(38a)\begin{array}{c}{E}_{{\mathrm{\Gamma }}_B}^1=⟨D_A{\mathrm{\Gamma }}_B^{(1)}⟩=⟨D_A(t){\displaystyle \frac{{\mathrm{\Gamma }}_B(t+\mathrm{\Delta }t)-{\mathrm{\Gamma }}_B(t)}{\mathrm{\Delta }t}}⟩.\end{array}

Additionally, for the term of E_{f_B}^1, the system dynamics can induce a considerable bias,

(38b)\begin{array}{c}{E}_{f_B}^1=⟨D_A{\displaystyle \sum _{i=1}^N\left[{\displaystyle \frac{\mathrm{\partial }{f}_B\left({\mathit{x}}_B\right)}{\mathrm{\partial }{x}_i}}{f}_i\left({\mathit{x}}_i\right)\right]}⟩.\end{array}

According to Eq. (37), the correlation function E_{BA}^{2,1}(t^{\prime }) can be calculated based on the measured data. It is expected that the driving signal can be well designed so that both terms, E_{{\mathrm{\Gamma }}_B}^1 and E_{f_B}^1, approach to negligible quantities. For doing so the driving signal should have the following characteristics:

(i) The driving signal is sampled from a series of quantified stochastic numbers D_A(t_1), D_A(t_2)\cdots ,D_A(t_l),\cdots ,D_A(t_L), i.e.,

D_A(t)=D_A(t_l),t\in (t_{l-1},t_l].

where the tunable time interval t_l-t_{l-1} should be much smaller than the system characteristic time (so D_A(t) is called fast-varying driving).

(ii) The statistical mean of the driving signal is zero, i.e.,

⟨D_A(t)⟩=\underset{t_L\to \mathrm{\infty }}{lim}{\displaystyle \frac{1}{t_L}}{\displaystyle \sum _{l=1}^L{D}_A(t_l)(t_l-t_{l-1})\to 0.}

The zero-mean property can help to weaken the effect of the driving signals on the system and reduce the bias in the correlation function. E_{f_B} is proved to be equivalently small as (t_l-t_{l-1}) by theoretical derivation and can be ignored,

E_{f_B}\to O(t_l-t_{l-1}).

In addition, the dynamics of the systems depend not only on the current state but also on the previous states due to the effects of time delay on signal propagation between two nodes, which may induce complicated or practically important behaviors of system dynamics (e.g., excitability or chaotic dynamics). Therefore, the reconstruction of time delay {\tau }_{BA} is essential. According to Eq. (37), to distinguish time delay in the autocorrelation the driving signal should have the following characteristics.

(iii) The driving signal is nonperiodic.

(iv) The driving signal has short correlation time, i.e.,

\begin{array}{rl}& ⟨D_A(t_m)D_A(t_n)⟩\ne 0,m=n,\\ & ⟨D_A(t_m)D_A(t_n)⟩=0,m\ne n.\end{array}

It is interesting that white noise, which is common and inherent in many practical systems, satisfies all the above requirements. Many previous studies on network reconstruction show that the network can generate rich distinctive data in the presence of white noises. Thus, we designed a stochastic sequence to simulate noise-like driving for network reconstruction:

(39a)\begin{array}{c}{D}_A(t)=D_A(t_l),t\in (l\mathrm{\Delta }{t}_D-\mathrm{\Delta }{t}_D,l\mathrm{\Delta }{t}_D),\end{array}

where the signal changes with a fixed time interval \mathrm{\Delta }{t}_D. The value of D_A(t_l) is randomly sampling from a Gaussian distribution

(39b)\begin{array}{c}{D}_A(t_l)\sim N(0,Q_{DA}/\mathrm{\Delta }{t}_D).\end{array}

Q_{DA} is the noise intensity. For simplicity, we set \mathrm{\Delta }{t}_D=\mathrm{\Delta }t in our study, and hence,

\begin{array}{rl}& \underset{L\to \mathrm{\infty }}{lim}\frac{1}{L}{\sum }_{l=1}^L{D}_A(t_l)\to 0,\\ & \underset{L\to \mathrm{\infty }}{lim}\frac{1}{L}{\sum }_{l=1}^L{D}_A(t_l)D_A(t_l)\mathrm{\Delta }t\to Q_{DA},\end{array}

and for any k>0

\underset{L\to \mathrm{\infty }}{lim}{\displaystyle \frac{1}{L}}{\displaystyle \sum _{l=1}^{L-k}{D}_A(t_l)D_A(t_l+k\mathrm{\Delta }t)\to 0.}

Based on the above well-designed stochastic driving signal for direct interaction from A to B (d=1), Eq. (37) can be rewritten as

(40)\begin{array}{c}{E}_{BA}^{2,1}(t^{\prime })=\{\begin{array}{l}O(\mathrm{\Delta }t),t^{\prime }\ne {\tau }_{BA},\\ I_{BA}(A\to B)Q_{DA}+O(\mathrm{\Delta }t),t^{\prime }={\tau }_{BA},\end{array}\end{array}

leading to the solution

(41)\begin{array}{c}I(A\to B)=⟨{\displaystyle \frac{\mathrm{\partial }{f}_B\left({\mathit{x}}_B(t)\right)}{\mathrm{\partial }{x}_A(t-{\tau }_{BA})}}⟩,\end{array}

where I(A\to B) is the effective intensity of the interaction from A to B. We can reconstruct the fundamental information of this adjacent link, including the accurate connection strength and the time delay {\tau }_{BA}. However, the interaction between two direct-connected nodes can only provide limited information compared to the overall system, and it has been a important and challenging task to reconstruct indirect and even long-distant connections in the network. Intuitively, an indirect interaction from A to B cannot be inferred with only the measurable data of the destination node B, since the dynamics of node is also strongly influenced by other hidden nodes.

Considering a path P:A\to j_1\to B with one intermediate node j_1, i.e., the distance from A to B is d=2. The driving-induced term disappears in {\ddot{x}}_B due to \mathrm{\partial }{f}_B\left({\mathit{x}}_B(t)\right)/\mathrm{\partial }{x}_A(t-{\tau }_{BA})\equiv 0. We should go further to calculate the higher-order derivatives of x_B(t), which is given by

(42)\begin{array}{c}\begin{array}{rl}{x}_B^{(3)}(t)\approx {\left({\ddot{x}}_B(t)\right)}^{(1)}=& \frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{\mathrm{\partial }{f}_B\left({\mathit{x}}_B(t)\right)}{\mathrm{\partial }{x}_{j_1}(t-{\tau }_{B{j}_1})}\right)f_{j_1}\left({\mathit{x}}_{j_1}(t-{\tau }_{B{j}_1})\right)\\ & +\frac{\mathrm{\partial }{f}_B\left({\mathit{x}}_{j_1}(t)\right)}{\mathrm{\partial }{x}_{j_1}(t-{\tau }_{B{j}_1})}\frac{\mathrm{\partial }{f}_{j_1}\left({\mathit{x}}_{j_1}(t-{\tau }_{B{j}_1})\right)}{\mathrm{\partial }{x}_A(t-{\tau }_{B{j}_1}-{\tau }_{j_{1}A})}(f_A\left({\mathit{x}}_A(t-{\tau }_{B{j}_1}-{\tau }_{j_{1}A})\right)\\ & +{\mathrm{\Gamma }}_A(t-{\tau }_{B{j}_1}-{\tau }_{j_{1}A})+D_A(t-{\tau }_{B{j}_1}-{\tau }_{j_{1}A})).\end{array}\end{array}

And the driving signal D_A can propagate along the path A\to j_1\to B via the following driving-induced term,

{\displaystyle \frac{\mathrm{\partial }{f}_B\left({\mathit{x}}_{j_1}(t)\right)}{\mathrm{\partial }{x}_{j_1}(t-{\tau }_{B{j}_1})}}{\displaystyle \frac{\mathrm{\partial }{f}_{j_1}\left({\mathit{x}}_{j_1}(t-{\tau }_{B{j}_1})\right)}{\mathrm{\partial }{x}_A(t-{\tau }_{B{j}_1}-{\tau }_{j_{1}A})}}{D}_A(t-{\tau }_{B{j}_1}-{\tau }_{j_{1}A}).

To detect the driving force in x_B(t), we calculate the correlation between D_A(t) and x_B^{(3)}(t), which is given by

(43)\begin{array}{c}\begin{array}{rl}{E}_{BA}^{3,1}(t^{\prime })=& \frac{1}{T-t^{\prime }}{\int }_0^{T-t^{\prime }}{D}_A(t)x_B^{(3)}(t+t^{\prime })\mathrm{d}t\\ =& \{\begin{array}{ll}O(\mathrm{\Delta }t),& t^{\prime }\ne {\tau }_{B{j}_1}+{\tau }_{j_{1}A},\\ I(A\to j_1\to B)Q_{DA}+O(\mathrm{\Delta }t),& t^{\prime }={\tau }_{B{j}_1}+{\tau }_{j_{1}A},\end{array}\end{array}\end{array}

where

I(A\to j_1\to B)=⟨{\displaystyle \frac{\mathrm{\partial }{f}_B\left({\mathit{x}}_{j_1}(t)\right)}{\mathrm{\partial }{x}_{j_1}(t-{\tau }_{B{j}_1})}}\frac{\mathrm{\partial }{f}_{j_1}\left({\mathit{x}}_{j_1}(t-{\tau }_{B{j}_1})\right)}{\mathrm{\partial }{x}_A(t-{\tau }_{B{j}_1}-{\tau }_{j_{1}A})}⟩.

Finally, we consider the case d>2, where there are d-1 different hidden nodes along the path from A to B. Denote number \nu as the distance of intermediate node j_{\nu } away from A, we have the path from A to B as

(44a)\begin{array}{c}P:A\to j_1\to \cdots j_{\nu }\to \cdots \to j_{d-1}\to B.\end{array}

The structural effective interaction intensity from A to B along this path P is computed as

(44b)\begin{array}{c}{I}_{BA}(P)=⟨{\displaystyle \frac{\mathrm{\partial }{f}_B}{\mathrm{\partial }{x}_{j_{d-1}}}}\cdot \prod _{\nu =2}^{d-1}{\displaystyle \frac{\mathrm{\partial }{f}_{j_{\nu }}}{\mathrm{\partial }{x}_{j_{\nu -1}}}}\cdot {\displaystyle \frac{f_{j_1}}{\mathrm{\partial }{x}_A}}⟩,\end{array}

where for i,j=j_1,j_2,..,j_{d-1},B, we have

{\displaystyle \frac{\mathrm{\partial }{f}_i}{\mathrm{\partial }{x}_j}}={\displaystyle \frac{\mathrm{\partial }{f}_i\left({\mathit{x}}_i(t+{\tau }_{iA})\right)}{\mathrm{\partial }{x}_j(t+{\tau }_{jA}(P))}},

and time delay from A to B through P reads

(44c)\begin{array}{c}{\tau }_{BA}(P)={\tau }_{j_{1}A}+{\displaystyle \sum _{\mu =2}^{d-1}{\tau }_{j_{\mu }{j}_{\mu -1}}+{\tau }_{B{j}_{d-1}}.}\end{array}

To detect the driving force component, we calculate the correlation between the driving signal D_A(t) and the d+1th derivative of x_B(t), which is given by

(45)\begin{array}{c}{E}_{BA}^{d+1,1}(t^{\prime })\approx \{\begin{array}{ll}O(\mathrm{\Delta }t),& t^{\prime }\ne {\tau }_{BA}(P),\\ I_{BA}(P)Q_{DA}+O(\mathrm{\Delta }t),& t^{\prime }={\tau }_{BA}(P),\end{array}\end{array}

where I_{BA}(P) and {\tau }_{BA}(P) are given in Eqs. (44b) and (44c), respectively. The calculation of E_{BA}^{d+1,1}(t^{\prime }) is only dependent on the data of x_B(t) and D_A(t). According to the theoretical deduction, there is a discontinuity in Eq. (45) at t^{\prime }={\tau }_{BA}, which is called the characteristic discontinuity of the certain path P.

Note that there are two more terms E_{{\mathrm{\Gamma }}_B}^d and E_{f_B}^d in the correlation function. The leading terms of the random deviation caused by noises and the bias caused by the system dynamics are given by

(46a)\begin{array}{c}{E}_{{\mathrm{\Gamma }}_B}^d=⟨D_A(t){\mathrm{\Gamma }}_B^{(d)}⟩,\end{array}

(46b)\begin{array}{c}{E}_{f_B}^d=⟨D_A(t){{\displaystyle \sum }}_{i,j=1}^N{F}_f^d⟩,\end{array}

where F_f^d is a simplified function given as

F_f^d\left(f_i,{\displaystyle \frac{\mathrm{\partial }{f}_i}{\mathrm{\partial }{x}_j}},{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left({\displaystyle \frac{\mathrm{\partial }{f}_i}{\mathrm{\partial }{x}_j}}\right),\cdots ,{\displaystyle \frac{{\mathrm{d}}^d}{\mathrm{d}{t}^d}}\left({\displaystyle \frac{\mathrm{\partial }{f}_i}{\mathrm{\partial }{x}_j}}\right)\right).

Under the condition of the well-designed stochastic driving signal, these two quantities should be small and negligible for long time averages.

Due to the discretizations of measurement time, the above correlation can be rewritten as

(47a)\begin{array}{c}{E}_{BA}^{d+1,1}(k\mathrm{\Delta }t)={\displaystyle \frac{\mathrm{\Delta }t}{T-k\mathrm{\Delta }t}}{\displaystyle \sum _{l=1}^{L-k}{D}_A(t_l)x_B^{d+1}(t_l+k\mathrm{\Delta }t),}\end{array}

where \mathrm{\Delta }t is the time interval in the recording process. The reconstructed interaction intensity can be represented by

(47b)\begin{array}{c}{\tilde{I}}_{BA}(P)=E_{BA}^{d+1,1}({\tilde{\tau }}_{BA}(P))/Q_{DA},\end{array}

where {\tilde{\tau }}_{BA}(P) is the experimental time delay (ceil(x) means rounding up number x)

(47c)\begin{array}{c}{\tilde{\tau }}_{BA}(P)=\mathrm{\Delta }t\cdot ceil({\tau }_{BA}(P)/\mathrm{\Delta }t).\end{array}

We first consider a linear system without the inherent noise where an external driving signal is loaded on node A (A=1 as an example)

(48a)\begin{array}{c}{\dot{x}}_{A=1}(t)=-A_{11}{x}_1(t)+{\displaystyle \sum _{j\ne 1}^{}{A}_{1j}{x}_j(t-{\tau }_{1j})+D_1(t),}\end{array}

and for the other nodes i\ne A

(48b)\begin{array}{c}{\dot{x}}_i(t)=-A_{ii}{x}_i(t)+{\displaystyle \sum _{j\ne i}^{}{A}_{ij}{x}_j(t-{\tau }_{ij}),}\end{array}

where A_{ii}=4 and \hat{A} is the connection matrix. Node B (B=2) is perceptible and the interaction paths from A to B are tested for deep reconstruction. Both the values of measurement intervals (\mathrm{\Delta }t) and the loading interval of driving signal (\mathrm{\Delta }{t}_D) are set to be the same, \mathrm{\Delta }{t}_D=\mathrm{\Delta }t={10}^{-2}.

In the following part, we demonstrate that our method works in the condition of multiple paths with different distances. We generate a network with N=5 named as Net,

Net:\hat{\mathit{A}}=\left(\begin{array}{ccccc}0& 0& 0& 0& 0\\ 0& 0& 1.2& 0& 1.5\\ 1.2& 0& 0& 0& 0\\ 1.5& 0& 0& 0& 0\\ 0& 0& 0& 1.5& 0\end{array}\right).

The parameters of time delay are as follows:

\hat{\tau }=\left(\begin{array}{ccccc}-& -& -& -& -\\ -& -& 0.835& -& 0.720\\ 0.739& -& -& -& -\\ 0.604& -& -& -& -\\ -& -& -& 0.945& -\end{array}\right).

There are two paths from A (A=1) to B (B=2) with different distances, i.e., P_1:1\to 3\to 2 and P_2:1\to 4\to 5\to 2, as shown in Fig. 4(a). To reconstruct the multiple paths of the Net, the correlations E_{BA}^{n,1}(t^{\prime }) are calculated based on data of x_A(t) and x_B(t). We find that a singularity emerges in E_{BA}^{3,1}(t^{\prime }) at t^{\prime }=1.57 firstly as shown in Fig. 4(b), which means there is a shortest path with d=2, and another singularity emerges in E_{BA}^{4,1}(t^{\prime }) at t^{\prime }=2.27 as shown in Fig. 4(c), informing there is a second shortest path with d=3. For convenience, we color the two characteristic singular peaks in (b) and (c) to distinguish the shortest path and the second shortest path (red for the shortest path of d=2 and blue for the second shortest path of d=3). The hollow red and blue points represent the positions of theoretical value I_{BA} and {\tau }_{BA}, also red for the shortest and blue for the second shortest path. Note that the value of E_{BA}^{4,1}(t^{\prime }) at t^{\prime }=1.57 is much larger than that in E_{BA}^{3,1}(t^{\prime }). This is because correlation function comes from the same data of x_B(t), and correlation in E_{BA}^{3,1}(t^{\prime }) is enlarged by 1/\mathrm{\Delta }t in E_{BA}^{4,1}(t^{\prime }), not considered as a new singularity. Similarly, the same driving signal exist in multiple paths with different distances and time delays, which can cause random background perturbation in E_{BA}^{4,1}(t^{\prime }). Thus we consider only the amplitude of singularity at t^{\prime }=2.27 as the effective intensity of the second shortest path. The reconstructed results are shown in Table 1, and the accuracy of intensities for both paths is very close to 1. According to the above numerical simulation, a sketch map of the reconstructed network is shown in Fig. 4(d).

A considerable amount of works in the field of biological, chemical, and physical systems have been dedicated to excitable media [56–59]. The interest in this topic is motivated by the fact that wave propagation in these media provides an efficient mechanism for communication between distant locations. Seminal examples include the conduction of electrical impulses along nerve axons. A modified version of FitzHugh-Nagumo model, called the Bär Model is investigated in this paper [60–62], described as follows ({\delta }_{1,i\ne 1}=0, {\delta }_{1,1}=1):

(49a)\begin{array}{c}{\displaystyle \frac{\mathrm{d}{u}_i}{\mathrm{d}t}}=-{\displaystyle \frac{1}{e}}{u}_i(u_i-1)(u_i-{\displaystyle \frac{v_i+b}{a}})+W_i+D_1{\delta }_{1,i},\end{array}

(49b)\begin{array}{c}{\displaystyle \frac{\mathrm{d}{v}_i}{\mathrm{d}t}}=f(u_i)-v_i,\end{array}

(49c)

where u_i and v_i are the membrane potential and inhibitory currents of the ith node (i=1,2,\cdots ,N). W_i is the coupling effect from other nodes, W_i=\sum _{j=1}^N{M}_{ij}(u_j(t-{\tau }_{ij})-u_i(t)), (i\ne j). M_{ij} is the connection intensity and {\tau }_{ij} is the transmission time delay between adjacent nodes. The function {\delta }_{1,i} means that the external signal is loaded on node 1. Parameters of identical local dynamics for all nodes are a=0.84;b=0.07;e=0.04.

A randomly generated network of size N=30 is considered (more details are shown in caption of Fig. 5). The schematic diagram of the network is shown in Fig. 5(a), where all the nodes are arranged according to the shortest distance from A. For example, the shortest distances from node A=1 to nodes 19 and 26 are d(1\to 19)=1 and d(1\to 26)=4. The driving intensity on node A=1 is set Q_{DA}=0.1. According to the reconstruction process, we detect the distance {\tilde{d}}_{BA}B=2,3,\cdots ,N based on the collected data from this excitable network. If the distance of the shortest paths from A to B are accurately reconstructed, the node B is marked with white hollow circle in Fig. 5(a), otherwise it is marked in grey. We can further infer and calculate effective intensity and time delay for paths with correct distance [Eq. (47)]. Figure 5(b) presents the reconstructed intensities {\tilde{I}}_{BA}(P) against the corresponding structural values I_{BA}(P) for all reconstructed paths together, where the red dots for the shortest paths and the blue ones for the second shortest paths. The results are mostly quite accurate for the shortest paths, however there are some relatively large errors for the second shortest ones.

To optimize the experimental results, we increase the amount of data and recollect N_{run}=20 times different dataset. In each run the system evolves from different initial states and the correlation function in the ith dataset is denoted as E_{BA,i}^{d+1,1}(t^{\prime }). The averaged correlation is thus given by

(50)\begin{array}{c}{\overline{E}}_{BA}^{d+1,1}(t^{\prime })={\displaystyle \frac{1}{N_{run}}}{\displaystyle \sum _{i=1}^{Times}{E}_{BA,i}^{d+1,1}(t^{\prime }).}\end{array}

The new reconstruction results based on the averaged correlation function {\overline{E}}_{BA}^{d+1,1}(t^{\prime }) are shown in Fig. 5(c), which is greatly improved where all values of both the shortest and the second shortest paths are all distributed near the diagonal line, i.e., {\tilde{I}}_{BA}\approx I_{BA}. In addition, the transmission time delays of these paths are perfectly inferred, shown in Fig. 5(d).

5 Uncovering hidden nodes and hidden link structure by applying noise injection

Neural networks have attracted huge attention in practice, where a number of units are often not accessible while important for realistic functions, e.g., various neural systems [11–15]. Let us consider again the neural network with the units described by the B\ddot{a}r model of Eq. (49).

Figure 6(a) illustrates an example of an artificial network with N=8, where M=6 red nodes are accessible, that is, their data can be measured, certain signals can be injected, and N-M=2 black nodes and all black links are hidden. The task is to detect the hidden nodes and infer all links with available data of accessible nodes only. To do so, white noise is injected into an accessible node, for example, node A, and the data of two accessible nodes A and B (B\ne A) are measured, then the differential equation for A is changed to

(51a)\begin{array}{c}\begin{array}{rl}\frac{\mathrm{d}{x}_A(t)}{\mathrm{d}t}=& -\frac{1}{\epsilon }{x}_A(t)(x_A(t)-1)(x_A(t)-\frac{y_A(t)+b}{a})\\ & +\sum _{j=1,j\ne A}^N{w}_{Aj}(x_j(t-{\tau }_{Aj})-x_A\left(t\right))+{\mathrm{\Gamma }}_A(t),\end{array}\end{array}

where noise {\mathrm{\Gamma }}_A(t) injected to node A is the signal that provides rich information on the hidden structure in the network. The noise has a correlation time much smaller than the characteristic time of network dynamics, and is approximated by white noise, satisfying

(51b)\begin{array}{c}⟨{\mathrm{\Gamma }}_A(t)⟩=0;⟨{\mathrm{\Gamma }}_A(t){\mathrm{\Gamma }}_A(t+t^{\prime })⟩=Q_A\delta (t^{\prime }).\end{array}

Note, here {\mathrm{\Gamma }}_A(t) is randomly generated and its time series is unknown, unlike the completely known noise sequences D_A(t) applied in Eq. (31). We can measure the output data from accessible nodes A and B as

(51c)\begin{array}{c}\begin{array}{rl}& x_A(t_1),x_A(t_2),\cdots ,x_A(t_k),\cdots ,x_A(t_L),\\ & x_B(t_1),x_B(t_2),\cdots ,x_B(t_k),\cdots ,x_B(t_L),\end{array}\end{array}

with the sampling interval \mathrm{\Delta }t

A and B can be selected among all the accessible nodes of Fig. 6(a).

Rich information about the network structure can be extracted by diverse correlation calculations with data of x_A(t) and x_B(t) as

(52)\begin{array}{c}{C}_{BA}^{\nu ,1}(t^{\prime }=k^{\prime }\mathrm{\Delta }t)={\displaystyle \frac{1}{L}}{\displaystyle \sum _{k=1}^L{x}_B^{(\nu )}(t_{k+k^{\prime }})x_A^{(1)}(t_k),}\end{array}

where x_i^{(\nu )} is the \nuth-order “derivative” of x_i, defined by

x_i^{(\nu )}(t_k)={\displaystyle \frac{x_i^{(\nu -1)}(t_{k+1})-x_i^{(\nu -1)}(t_k)}{\mathrm{\Delta }t}},x_i^{(0)}(t_k)=x_i(t_k),

where we take x_i(t_k)=0 for all k\le 0 and k>L. In Eq. (52) correlations for different \nu {}^{\prime }s include various and often mutually independent information. If there is an interaction path from A to B with distance d, we have

(53a)\begin{array}{c}{I}_{BA}(d)=w_{j_{1}A}{\displaystyle \prod _{\mu =2}^{d-1}{w}_{j_{\mu }{j}_{\mu -1}}{w}_{B{j}_{d-1}},}\end{array}

(53b)\begin{array}{c}{\tau }_{BA}(d)={\tau }_{j_{1}A}+{\displaystyle \sum _{\mu =2}^{d-1}{\tau }_{j_{\mu }{j}_{\mu -1}}+{\tau }_{B{j}_{d-1}},}\end{array}

where I_{BA}(d), {\tau }_{BA}(d) indicates the actual interaction intensity and total time delay along the path of length d

(53c)\begin{array}{c}A\to j_1\to \cdots \to j_{\mu }\to \cdots \to j_{d-1}\to B.\end{array}

Based on Eq. (53) and for an actual interaction path from node A to node B with distance d, we have

x_B^{(d+1)}(t)={\tilde{I}}_{BA}(d)x_A^{(1)}(t-{\tilde{\tau }}_{BA}(d))+\text{other terms}.

Equation (52) can be further derived to

(54)\begin{array}{c}\begin{array}{rl}{E}_{BA}^{\nu =d+1,1}(t^{\prime })=& C_{BA}^{d+1,1}(t^{\prime })\mathrm{\Delta }t\\ =& \{\begin{array}{cc}0(\mathrm{\Delta }t),& t_k^{\prime }\ne {\tau }_{BA}(d);\\ Q_A{\tilde{I}}_{BA}(d)+0(\mathrm{\Delta }t),& t_k^{\prime }={\tau }_{BA}(d).\end{array}\end{array}\end{array}

Hence, the interaction from A to B of I_{BA}(d) and {\tau }_{BA}(d) can be inferred from accessible data of nodes A and B by computation of the peak height I_{BA}(d)\approx {\tilde{I}}_{BA}(d) and peak position {\tau }_{BA}(d)\approx {\tilde{\tau }}_{BA}(d) in Eq. (54), where {\tilde{I}}_{BA}(d) and {\tilde{\tau }}_{BA}(d) refer to the reconstructed interaction intensity and total time delay along the path, respectively.

By analysis of the output data of the accessible nodes of Eq. (51c), rich information on paths can be extracted with different distances using Eqs. (53) and (54) from d=1 to, in principle, any large d if data of sufficiently large length L\gg 1 and small enough sampling interval \mathrm{\Delta }t\ll 1 are available. In addition, this information can be employed to infer interaction structure in networks not only between accessible nodes (d=1), but also between accessible and hidden nodes (d=2), and even between multiple hidden nodes (d\ge 3).

We demonstrate the results of the above analyses by computing a simple network structure in Fig. 6(a) with data of accessible nodes only. Figures 6(e)−(j) present some cases of E_{BA}^{d+1,1}(t^{\prime }) calculated from Eq. (54) for several d=2 and d=3 pairs with nonzero peaks. In Figs. 6(b, e, f), it can be concluded that there is definitely one hidden node along each interaction path 1\to 3,4 and 2\to 3,4. Without further analysis, the four peaks in Figs. 6(e, f) may be associated with any of the four cases from one to four hidden nodes. However, by integrating the correlations of all the four peaks in addition to appropriate computation and derivation, we conclude definitely that all the four paths pass through a single hidden node (say node 7). For instance, by comparing the peak heights (identifying Q_A{\tilde{I}}_{BA}(d)) and peak positions (showing {\tilde{\tau }}_{BA}(d)), we can convincingly reach the above conclusion, that is, if equalities

(55a)\begin{array}{c}\begin{array}{rl}& {\tilde{\tau }}_{B_1{A}_1}(d=2)-{\tilde{\tau }}_{B_1{A}_2}(d=2)\\ & ={\tilde{\tau }}_{B_2{A}_1}(d^{\prime })-{\tilde{\tau }}_{B_2{A}_2}(d^{\prime })\\ & ={\tilde{\tau }}_{h_1{A}_1}(1)-{\tilde{\tau }}_{h_1{A}_2}(1)\\ & \approx {\tau }_{h_1{A}_1}(1)-{\tau }_{h_1{A}_2}(1),\end{array}\end{array}

(55b)\begin{array}{c}\begin{array}{rl}& {\tilde{I}}_{B_1{A}_1}(d=2)/{\tilde{I}}_{B_1{A}_2}(d=2)\\ & ={\tilde{I}}_{B_2{A}_1}(d^{\prime })/{\tilde{I}}_{B_2{A}_2}(d^{\prime })\\ & ={\tilde{I}}_{h_1{A}_1}(1)/{\tilde{I}}_{h_1{A}_2}(1)\\ & \approx I_{h_1{A}_1}(1)/I_{h_1{A}_2}(1),\end{array}\end{array}

are justified, a hidden node h_1 and its inputs from A_1 and A_2 nodes can be determined; and if equations

(55c)\begin{array}{c}\begin{array}{rl}& {\tilde{\tau }}_{B_1{A}_1}(d=2)-{\tilde{\tau }}_{B_2{A}_1}(d=2)\\ & ={\tilde{\tau }}_{B_1{A}_2}(d^{\prime })-{\tilde{\tau }}_{B_2{A}_2}(d^{\prime })\\ & ={\tilde{\tau }}_{B_1{h}_2}(1)-{\tilde{\tau }}_{B_2{h}_2}(1)\\ & \approx {\tau }_{B_1{h}_2}(1)-{\tau }_{B_2{h}_2}(1),\end{array}\end{array}

(55d)\begin{array}{c}\begin{array}{rl}& {\tilde{I}}_{B_1{A}_1}(d=2)/{\tilde{I}}_{B_2{A}_1}(d=2)\\ & ={\tilde{I}}_{B_1{A}_2}(d^{\prime })/{\tilde{I}}_{B_2{A}_2}(d^{\prime })\\ & ={\tilde{I}}_{B_1{h}_2}(1)/{\tilde{I}}_{B_2{h}_2}(1)\\ & \approx I_{B_1{h}_2}(1)/I_{B_2{h}_2}(1),\end{array}\end{array}

are valid, a single hidden node h_2 and its outputs to nodes B_1 and B_2 can be inferred. In Eq. (55), all these d-1 and d^{\prime }-1 nodes are hidden in the considered paths. It can be easily verified that the four peaks of Figs. 6(e, f) satisfy all Eqs. (55a)−(55d) with d=d^{\prime }=2, A_1=1, A_2=2, B_1=3, B_2=4, and identify h_1=h_2=7 and its inputs and outputs by blue color shown in Fig. 6(b). The four peaks of Figs. 6(g, h) satisfy Eqs. (55a) and (55b) with d=2, d^{\prime }=3, A_1=1, A_2=2, B_1=3, B_2=6, and determine h_1=7 and its inputs and output by purple color as shown in Fig. 6(c). The four peaks of Figs. 6(i, j) justify Eqs. (55c) and (55d) with d=2, d^{\prime }=3, A_1=4, A_2=2, B_1=6, B_2=5, and identify h_2=8 and its input and outputs by purple color as shown in Fig. 6(d). With substructure, the motif in Fig. 6(b), hidden node 7 can uniquely be identified, that is, identify node 7 and all its inputs from observable nodes and outputs to observable nodes from all the rest of the hidden nodes and other links. The motif in Fig. 6(c) [Fig. 6(d)] determines hidden node 7 (8) and identifies the hidden node and its inputs (outputs) links with observable nodes from any other hidden nodes and hidden links. Hence, the simplest detectable structure is distinguished as a reconstruction motif, with which a hidden node and some of its links with accessible nodes and even the link form hidden node 7 to hidden node 8 can be recognized according to the analysis of measured accessible node data in the motifs. In this work, three different motifs are shown in Figs. 6(b−d), with which we can evidently distinguish all hidden nodes and hidden links in the network of Fig. 6(a).

We demonstrate the above in a simple network how to use noise injection and analyze output data of the accessible nodes to uncover hidden nodes and links. This method can be directly extended to more complex and larger networks by utilizing the network motifs in Figs. 6(b−d). A complex network of size N=20 is illustrated in Fig. 7(a). In this network, ten red nodes [Fig. 7(b)] are accessible, while all the other black nodes in Fig. 7(a) are hidden, and the interaction structure of the whole network and network dynamics are all unknown. The time series [e.g., Eq. (51c)] of the accessible nodes are generated using Eqs. (49, 51) with noise being injected to one of the accessible nodes [i.e., node A in Fig. 7(b)]. The aim is network reconstruction, that is, to explore additional hidden nodes in Fig. 7(a) and infer hidden links associated with these uncovered nodes.

Based on the correlation results of d=1, all direct interactions between accessible nodes [red arrows in Fig. 7(c)] can be easily predicted. In addition, hidden nodes and links can be inferred by using E_{BA}^{d+1,1}(t{}^{\prime }) with d=2 [blue nodes and links in Fig. 7(d)] through motif in Fig. 6(b) and by using E_{BA}^{d+1,1}(t{}^{\prime }) with d=2,3 [purple nodes and links in Fig. 7(e)] based on motifs in Figs. 6(c, d). By integrating the results and systematically applying Eq. (55) with network motifs in Figs. 6(b)−(d), we obtain the structures of Figs. 7(d) and (e) step by step, fully reconstructing network of Fig. 7(a). In Figs. 7(f, g), we plot ratio of intensities \tilde{R}=\frac{{\tilde{I}}_{h_{\mu }{A}_1}}{{\tilde{I}}_{h_{\mu }{A}_2}} and time delay differences \mathrm{\Delta }\tilde{\tau }={\tilde{\tau }}_{h_{\mu }{A}_1}-{\tilde{\tau }}_{h_{\mu }{A}_2} between any pairs of inputs and also \tilde{R}=\frac{{\tilde{I}}_{B_1{h}_{\mu }}}{{\tilde{I}}_{B_2{h}_{\mu }}} and \mathrm{\Delta }\tilde{\tau }={\tilde{\tau }}_{B_1{h}_{\mu }}-{\tilde{\tau }}_{B_2{h}_{\mu }} for outputs with \mu =1,2,\cdots ,N-M for all hidden nodes, and compare these reconstructed values with actual ones (R=\frac{I_{h_{\mu }{A}_1}}{I_{h_{\mu }{A}_2}} or R=\frac{I_{B_1{h}_{\mu }}}{I_{B_2{h}_{\mu }}} and \mathrm{\Delta }\tau ={\tau }_{h_{\mu }{A}_1}-{\tau }_{h_{\mu }{A}_2} or \mathrm{\Delta }\tau ={\tau }_{B_1{h}_{\mu }}-{\tau }_{B_2{h}_{\mu }}). All dots distribute along the diagonal lines, confirming acceptable reconstruction of these network quantities. It should be indicated that for each hidden node, two arbitrary factors, one for time delays and the other for interaction intensity, are not determined. Note that network Fig. 7(e) has 50\mathrm{\%} hidden nodes and more than 85\mathrm{\%} links output from or input to hidden nodes. It is a fascinating success for us to conclude the complete reconstruction. Moreover, all the approaches used in the section can be applied to any networks where signal injection to nodes is difficult while internal fast-varying noises exist in network nodes. This advantage is important in reconstructions of many practical networks.

In Fig. 7(e), all the hidden nodes and hidden links of network Fig. 7(a) have been reconstructed by employing the motifs in Figs. 6(b−d) only. With a given network, a larger ratio of hidden nodes relative to observable nodes can certainly render the task of uncovering the hidden structure more difficult. In these cases, the calculation of additional correlations of data of measurable nodes and identification of additional motifs with larger d^{\prime }s (e.g., d>3) may aid in complete reconstruction. Nevertheless, the conditions for complete reconstruction are still unknown and significantly warrants further investigation [41].

6 Discussion and future perspectives

More than two decades of intensive researches on complex networks have resulted in a large body of knowledge about network structures and their effects on various dynamical processes [63, 64]. A typical approach in the field is to implement a particular dynamical process of interest on networks whose connecting topologies are completely specified [65]. In addition to this line of research that is necessary to discover and understand various fundamental phenomena in complex networks, the “inverse” or “reverse engineering” problem of predicting network structure and dynamics from data becomes extremely important. The basic hypothesis underlying the inverse problem is that the detailed structure of the network and the nodal dynamics are unknown, but only a limited set of signals or time series measured from some nodes of the network are available [5]. The question is whether the intrinsic structure of the network and the dynamical processes can be inferred solely from the set of measured time series and how to reconstruct the network structure if the answer is yes. Data-based reverse engineering of complex networks, with great application potential, is important not only to advance network science and engineering at a fundamental level but also to meet the need to address an array of applications in complex practical networks [14, 18, 25, 66–76]. For example, inference of some practical neural networks with experimentally measured data has been conducted [77–80].

In our review, two key points, existence of fast-varying noises and time-delay for signal propagation through interactions in the network, are considered. It is well known that these two factors are pervasive in realistic systems [81–85] and they can bring difficulties in analyses of network behaviors by producing random fluctuations and making network dynamics more complex. However, the applicability of the presented methods is limited by some strong assumptions (e.g., availability of high-frequency data measurement, feasibility of signal injection, favorable noise conditions, and so on), which may be not satisfied in many real-world systems. Nevertheless, we should emphasize that although these conditions are not always available, they may be available in many important practical experiments, or can be realized on purpose in many realistic measurements.

In this review, we have described an aspect of this field in detail how noise-induced correlations can be used for uncovering the full topology and nodal dynamical processes of complex, nonlinear oscillator networks, such as, for revealing the interaction patterns and dynamic structures of both linear and nonlinear networks in a compact and unified matrix form; for exploring both direct interaction and indirect distant links by computing multiple orders of correlations and considering effects of time-delays of signals propagation between network nodes; for detecting hidden nodes and the related hidden interaction structures with data of partial accessible nodes only.

Data-based reconstruction of complex networks is a broad field. A number of tasks of future development may be important: to identify and find data of high quality; to develop effective methods to extract as much as possible information from limited data for overcoming reconstruction difficulties; and to apply the developed methods to better understand, predict and control realistic social, biological, and technological complex networks.

References

[1] Y. Y. Liu , J. J. Slotine , and A. L. Barabási , Controllability of complex networks, Nature 473(7346), 167 (2011)

[2] G. Yan , J. Ren , Y. C. Lai , C. H. Lai , and B. Li , Controlling complex networks: How much energy is needed, Phys. Rev. Lett. 108(21), 218703 (2012)

[3] Y. Y. Liu and A. L. Barabási , Control principles of complex systems, Rev. Mod. Phys. 88(3), 35006 (2016)

[4] D. Marbach , J. Costello , R. Küffner , et al. Wisdom of crowds for robust gene network inference, Nat. Methods 9(8), 796 (2012)

[5] X. Han , Z. Shen , W. X. Wang , and Z. Di , Robust reconstruction of complex networks from sparse data, Phys. Rev. Lett. 114(2), 28701 (2015)

[6] Z. Shen , W. X. Wang , Y. Fan , Z. Di , and Y. C. Lai , Reconstructing propagation networks with natural diversity and identifying hidden sources, Nat. Commun. 5(1), 4323 (2014)

[7] W. X. Wang , Y. C. Lai , and C. Grebogi , Data based identification and prediction of nonlinear and complex dynamical systems, Phys. Rep. 644, 1 (2016)

[8] M. Nitzan , J. Casadiego , and M. Timme , Revealing physical interaction networks from statistics of collective dynamics, Sci. Adv. 3(2), e1600396 (2017)

[9] D. Q. Nykamp , Pinpointing connectivity despite hidden nodes within stimulus-driven networks, Phys. Rev. E 78(2), 21902 (2008)

[10] H. Haehne , J. Casadiego , J. Peinke , and M. Timme , Detecting hidden units and network size from perceptible dynamics, Phys. Rev. Lett. 122(15), 158301 (2019)

[11] D. Soudry , S. Keshri , P. Stinson , M. Oh , G. Iyengar , and L. Paninski , Efficient “Shotgun” inference of neural connectivity from highly sub-sampled activity data, PLOS Comput. Biol. 11(10), e1004464 (2015)

[12] B. A. W. Brinkman , F. Rieke , E. Shea-Brown , and M. A. Buice , Predicting how and when hidden neurons skew measured synaptic interactions, PLOS Comput. Biol. 14(10), e1006490 (2018)

[13] A. Das and I. R. Fiete , Systematic errors in connectivity inferred from activity in strongly recurrent networks, Nat. Neurosci. 23(10), 1286 (2020)

[14] W. Gilpin , Y. Huang , and D. B. Forger , Learning dynamics from large biological data sets: Machine learning meets systems biology, Curr. Opin. Syst. Biol. 22, 1 (2020)

[15] S. Panzeri , M. Moroni , H. Safaai , and C. D. Harvey , The structures and functions of correlations in neural population codes, Nat. Rev. Neurosci. 23(9), 551 (2022)

[16] K. A. Blaha , A. Pikovsky , M. Rosenblum , M. T. Clark , C. G. Rusin , and J. L. Hudson , Reconstruction of two-dimensional phase dynamics from experiments on coupled oscillators, Phys. Rev. E 84(4), 046201 (2011)

[17] R. Yang , Y. Wang , X. Wang , Z. Lei , Z. Zheng , and Y. Qian , The extended master stability function approach to alternating synchronization modes on networked oscillator systems, New J. Phys. 28(1), 013901 (2026)

[18] W. X. Wang , R. Yang , Y. C. Lai , V. Kovanis , and C. Grebogi , Predicting catastrophes in nonlinear dynamical systems by compressive sensing, Phys. Rev. Lett. 106(15), 154101 (2011)

[19] R. Q. Su , W. X. Wang , and Y. C. Lai , Detecting hidden nodes in complex networks from time series, Phys. Rev. E 85(6), 065201 (2012)

[20] W. X. Wang , Y. C. Lai , C. Grebogi , and J. Ye , Network reconstruction based on evolutionary-game data via compressive sensing, Phys. Rev. X 1(2), 021021 (2011)

[21] H. Zhao , Inferring the dynamics of “black-box” systems using a learning machine, Sci. China Phys. Mech. Astron. 64(7), 270511 (2021)

[22] Y. Du , Q. Li , H. Fan , M. Zhan , J. Xiao , and X. Wang , Inferring attracting basins of power system with machine learning, Phys. Rev. Res. 6(1), 013181 (2024)

[23] H. Luo , Y. Du , H. Fan , X. Wang , J. Guo , and X. Wang , Reconstructing bifurcation diagrams of chaotic circuits with reservoir computing, Phys. Rev. E 109(2), 024210 (2024)

[24] T. T. Gao and G. Yan , Autonomous inference of complex network dynamics from incomplete and noisy data, Nat. Comput. Sci. 2(3), 160 (2022)

[25] M. Timme , Revealing network connectivity from response dynamics, Phys. Rev. Lett. 98(22), 224101 (2007)

[26] W. X. Wang , Q. Chen , L. Huang , Y. C. Lai , and M. A. F. Harrison , Scaling of noisy fluctuations in complex networks and applications to network prediction, Phys. Rev. E 80(1), 016116 (2009)

[27] J. Ren , W. X. Wang , B. Li , and Y. C. Lai , Noise bridges dynamical correlation and topology in coupled oscillator networks, Phys. Rev. Lett. 104(5), 058701 (2010)

[28] W. X. Wang , J. Ren , Y. C. Lai , and B. Li , Reverse engineering of complex dynamical networks in the presence of time-delayed interactions based on noisy time series, Chaos 22(3), 033131 (2012)

[29] E. S. C. Ching , P. Y. Lai , and C. Y. Leung , Extracting connectivity from dynamics of networks with uniform bidirectional coupling, Phys. Rev. E 88(4), 042817 (2013)

[30] Z. Zhang , Z. Zheng , H. Niu , Y. Mi , S. Wu , and G. Hu , Solving the inverse problem of noise-driven dynamic networks, Phys. Rev. E 91(1), 012814 (2015)

[31] E. S. C. Ching , P. Y. Lai , and C. Y. Leung , Reconstructing weighted networks from dynamics, Phys. Rev. E 91(3), 30801 (2015)

[32] Y. Chen , S. Wang , Z. Zheng , Z. Zhang , and G. Hu , Depicting network structures from variable data produced by unknown colored-noise driven dynamics, Europhys. Lett. 113(1), 18005 (2016)

[33] E. S. C. Ching and H. C. Tam , Reconstructing links in directed networks from noisy dynamics, Phys. Rev. E 95(1), 10301 (2017)

[34] Y. Chen , Z. Zhang , T. Chen , S. Wang , and G. Hu , Reconstruction of noise-driven nonlinear networks from node outputs by using high-order correlations, Sci. Rep. 7(1), 44639 (2017)

[35] H. C. Tam , E. S. C. Ching , and P. Y. Lai , Reconstructing networks from dynamics with correlated noise, Physica A 502, 106 (2018)

[36] E. S. C. Ching and P. H. Tam , Effects of hidden nodes on the reconstruction of bidirectional networks, Phys. Rev. E 98(6), 062318 (2018)

[37] Z. Zhang , Y. Chen , Y. Mi , and G. Hu , Reconstruction of dynamic networks with time-delayed interactions in the presence of fast-varying noises, Phys. Rev. E 99(4), 042311 (2019)

[38] C. H. Cheng and P. Y. Lai , Efficient reconstruction of directed networks from noisy dynamics using stochastic force inference, Phys. Rev. E 106(3), 034302 (2022)

[39] X. Wang , Y. Mi , Z. Zhang , Y. Chen , G. Hu , and H. Li , Reconstructing distant interactions of multiple paths between perceptible nodes in dark networks, Phys. Rev. E 106(1), 014302 (2022)

[40] C. Sun , K. C. Lin , C. Y. Yeung , E. S. C. Ching , Y. T. Huang , P. Y. Lai , and C. K. Chan , Revealing directed effective connectivity of cortical neuronal networks from measurements, Phys. Rev. E 105(4), 044406 (2022)

[41] Z. Zhang , X. Wang , H. Li , Y. Chen , Z. Qu , Y. Mi , and G. Hu , Uncovering hidden nodes and hidden links in complex dynamic networks, Sci. China Phys. Mech. Astron. 67(4), 240511 (2024)

[42] H. Risken , Fokker–Planck Equation, pp 63–95, Springer, Berlin, Heidelberg, (1996)

[43] B. Gemao and P. Y. Lai , Effects of hidden nodes on noisy network dynamics, Phys. Rev. E 103(6), 062302 (2021)

[44] H. Huang , Effects of hidden nodes on network structure inference, J. Phys. A Math. Theor. 48(35), 355002 (2015)

[45] R. Schmidt , H. Haehne , L. Hillmann , J. Casadiego , D. Witthaut , B. Schäfer , and M. Timme , Inferring topology of networks with hidden dynamic variables, IEEE Access 10, 76682 (2022)

[46] Z. Yan , L. Gui , K. Xu , and Y. Lan , Reconstructing dynamics of complex systems from noisy time series with hidden variables, New J. Phys. 25(8), 083011 (2023)

[47] Y. Zhang , H. Li , Z. Zhang , Y. Qian , and V. Pandey , Network reconstruction from binary-state time series in presence of time delay and hidden nodes, Acta. Phys. Sin. 67, 203 (2020)

[48] R. Shi , W. Jiang , and S. Wang , Detecting network structures from measurable data produced by dynamics with hidden variables, Chaos 30(1), 013138 (2020)

[49] Y. Chen , C. Zhang , T. Chen , S. Wang , and G. Hu , Reconstruction of noise-driven nonlinear dynamic networks with some hidden nodes, Sci. China Phys. Mech. Astron. 60(7), 070511 (2017)

[50] R. Shi , C. Deng , and S. Wang , Detecting directed interactions of networks by random variable resetting, Europhys. Lett. 124(1), 18002 (2018)

[51] C. Zhang , Y. Chen , and G. Hu , Network reconstructions with partially available data, Front. Phys. (Beijing) 12(3), 128906 (2017)

[52] C. Zhang , Y. Chen , and G. Hu , Inference of targeted interactions of networks with data of driving and driven nodes only by applying fast-varying noise signals, Phys. Lett. A 381(31), 2502 (2017)

[53] S. G. Shandilya and M. Timme , Inferring network topology from complex dynamics, New J. Phys. 13(1), 013004 (2011)

[54] E. J. Candes , J. Romberg , and T. Tao , Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information, IEEE Trans. Inf. Theory 52(2), 489 (2006)

[55] E. J. Candès , J. K. Romberg , and T. Tao , Stable signal recovery from incomplete and inaccurate measurements, Commun. Pure Appl. Math. 59(8), 1207 (2006)

[56] Z. Qu , G. Hu , A. Garfinkel , and J. N. Weiss , Nonlinear and stochastic dynamics in the heart, Phys. Rep. 543(2), 61 (2014)

[57] Z. Zhang , G. Hu , Y. Zhang , and Z. Qu , Kramers rate theory of pacemaker dynamics in noisy excitable media, Phys. Rev. Lett. 129(4), 048101 (2022)

[58] Z. Zhang and Z. Qu , Bistable nerve conduction, Biophys. J. 121(18), 3499 (2022)

[59] Z. Zhang , Y. Zhang , and Z. Qu , Bistable spiral wave dynamics in electrically excitable media, Phys. Rev. E 108(6), 064405 (2023)

[60] M. Bär and M. Eiswirth , Turbulence due to spiral breakup in a continuous excitable medium, Phys. Rev. E 48, 1635 (1993)

[61] X. Liao , Q. Xia , Y. Qian , L. Zhang , G. Hu , and Y. Mi , Pattern formation in oscillatory complex networks consisting of excitable nodes, Phys. Rev. E 83(5), 056204 (2011)

[62] Y. Qian , X. Huang , G. Hu , and X. Liao , Structure and control of self-sustained target waves in excitable small-world networks, Phys. Rev. E 81(3), 036101 (2010)

[63] A. L. Barabási and R. Albert , Emergence of scaling in random networks, Science 286(5439), 509 (1999)

[64] D. J. Watts and S. H. Strogatz , Collective dynamics of mall-world’ networks, Nature 393(6684), 440 (1998)

[65] A. Arenas , A. Díaz-Guilera , J. Kurths , Y. Moreno , and C. Zhou , Synchronization in complex networks, Phys. Rep. 469(3), 93 (2008)

[66] T. S. Gardner , D. Bernardo , D. Lorenz , and J. J. Collins , Inferring genetic networks and identifying compound mode of action via expression profiling, Science 301(5629), 102 (2003)

[67] T. R. Lezon , J. R. Banavar , M. Cieplak , A. Maritan , and N. V. Fedoroff , Using the principle of entropy maximization to infer genetic interaction networks from gene expression patterns, Proc. Natl. Acad. Sci. USA 103(50), 19033 (2006)

[68] E. Schneidman , M. J. II Berry , R. Segev , and W. Bialek , Weak pairwise correlations imply strongly correlated network states in a neural population, Nature 440(7087), 1007 (2006)

[69] E. Bullmore and O. Sporns , Complex brain networks: graph theoretical analysis of structural and functional systems, Nat. Rev. Neurosci. 10(3), 186 (2009)

[70] S. Cocco , S. Leibler , and R. Monasson , Neuronal couplings between retinal ganglion cells inferred by efficient inverse statistical physics methods, Proc. Natl. Acad. Sci. USA 106(33), 14058 (2009)

[71] S. Hempel , A. Koseska , J. Kurths , and Z. Nikoloski , Inner composition alignment for inferring directed networks from short time series, Phys. Rev. Lett. 107(5), 054101 (2011)

[72] Z. Levnajić and A. Pikovsky , Network reconstruction from random phase resetting, Phys. Rev. Lett. 107(3), 034101 (2011)

[73] Z. Bar-Joseph , A. Gitter , and I. Simon , Studying and modelling dynamic biological processes using time-series gene expression data, Nat. Rev. Genet. 13(8), 552 (2012)

[74] Y. Roudi and G. Taylor , Learning with hidden variables, Curr. Opin. Neurobiol. 35, 110 (2015)

[75] J. Casadiego , M. Nitzan , S. Hallerberg , and M. Timme , Model-free inference of direct network interactions from nonlinear collective dynamics, Nat. Commun. 8(1), 2192 (2017)

[76] J. Casadiego , D. Maoutsa , and M. Timme , Inferring network connectivity from event timing patterns, Phys. Rev. Lett. 121(5), 054101 (2018)

[77] A. K. Fletcher and S. Rangan, Scalable inference for neuronal connectivity from calcium imaging, in: Z. Ghahramani, M. Welling, C. Cortes, N. Lawrence, and K. Q. Weinberger (Eds.), Advances in Neural Information Processing Systems, Vol. 27, Curran Associates, Inc., 2014

[78] L. Aitchison , L. Russell , A. M. Packer , J. Yan , P. Castonguay , M. Hausser , and S. C. Turaga , Model-based Bayesian inference of neural activity and connectivity from all-optical interrogation of a neural circuit, in: I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett (Eds.), Advances in Neural Information Processing Systems, Vol. 30, Curran Associates, Inc., (2017)

[79] I. Magrans de Abril , J. Yoshimoto , and K. Doya , Connectivity inference from neural recording data: Challenges, mathematical bases and research directions, Neural Netw. 102, 120 (2018)

[80] A. Banerjee , S. Chandra , and E. Ott , Network inference from short, noisy, low time-resolution, partial measurements: Application to elegans neuronal calcium dynamics, Proc. Natl. Acad. Sci. USA 120(12), e2216030120 (2023)

[81] B. Lindner , J. Garca-Ojalvo , A. Neiman , and L. Schimansky-Geier , Effects of noise in excitable systems, Phys. Rep. 392(6), 321 (2004)

[82] J. A. White , J. T. Rubinstein , and A. R. Kay , Channel noise in neurons, Trends in Neurosci. 23(3), 131 (2000)

[83] A. A. Faisal , L. P. J. Selen , and D. M. Wolpert , Noise in the nervous system, Nat. Rev. Neurosci. 9(4), 292 (2008)

[84] C. A. A. de Carvalho and H. M. Nussenzveig , Time delay, Phys. Rep. 364(2), 83 (2002)

[85] K. K. Sreenivasan and M. D’Esposito , The what, where and how of delay activity, Nat. Rev. Neurosci. 20(8), 466 (2019)

RIGHTS & PERMISSIONS

Higher Education Press