Compressed sensing

The general problem that CST addresses is to reconstruct a vector XεRN from linear measurements Y about X in the form

where YεRM and Φ is an M × N matrix. The striking feature of CS is that the number of measurements can be much less than the number of components of the unknown vector, that is, M<

where is the L1 norm of X and the matrix Φ satisfies the restricted isometry property. Solutions to the convex optimization are now standard353637383940. (More details of the CST can be found in Supplementary Note 1.) Our goal is to develop a framework to cast the problem of reconstructing propagation networks into form (1).

Reconstruction framework

To present our framework in a transparent manner, we first consider the relatively simple case where there is no hidden source. Further, we assume that the disease starts to propagate from a fraction of the infected nodes. As we will see, based on this framework, it is feasible to locate any hidden source based solely on time series after outbreak of infection. The state of an arbitrary node i is denoted as Si, where

Owing to the characteristic difference between the SIS dynamics and CP, we treat them separately (see Methods).

For the SIS dynamics, the probability of an arbitrary node i being infected by its neighbors at time t is

where λi is the infection rate of i, aij stands for the elements of the adjacency matrix (aij=1 if i connects to j and aij=0 otherwise), Sj(t) is the state of node j at t, and the superscript 01 denotes the change from susceptible state (0) to infected state (1). At the same time, the recovery probability of i is , where δi is the recovery rate of node i and the superscript 10 denotes the transition from infected state to susceptible state. Equation (4) can be rewritten as

If measurements at different times t=t1, t2, ···, tm are available, equation (5) can be written in the matrix form Ym × 1=Φm × (N−1)·X(N−1) × 1, where Y contains at different t, Φ is determined by the state Sj(t) of nodes except i, and X comprising the links and infection rates of i is sparse for a general network (see Methods). The main challenge here is that the infection probabilities at different times are not given directly from the time series of the nodal state.

To develop a method to estimate the probability from the nodal states, we set a threshold Δ pertaining to the normalized Hamming distance between strings composed of Sj(t) (j≠i) at different t to identify a base string at and a set of strings subject to the base. According to the law of large numbers, the probability can be estimated by the average over the state Si(t+1) at all proper time. By setting another threshold Θ associated with the normalized Hamming distance, we can identify a set of base strings. This process finally gives rise to a set of reconstruction equations in the matrix form:

where correspond to the time associated with m base strings and ‹·› denote the average over all satisfied t (see Methods). The vector Ym × 1 and the matrix Φm × (N−1) can then be obtained based solely on time series of nodal states and the vector X(N−1) × 1 to be reconstructed is sparse, rendering applicable the CS framework. As a result, we can achieve exact reconstruction of all neighbours of node i from relatively small amounts of observation. In a similar manner the neighbouring vectors of all other nodes can be uncovered from time series, enabling a full reconstruction of the whole network by matching the neighboring sets of all nodes.

For the CP dynamics, the infection probability of an arbitrary node i is given by

where ki is the degree of the node i, and the recovery probability is (see Methods). In close analogy to the SIS dynamics, we have

We then choose a series of base strings using a proper threshold Θ to establish a set of equations, expressed in the matrix form Ym × 1=Φm × (N−1)·X(N−1) × 1 (see Supplementary Note 2), where Φ has the same form as in equation (6), but Y and X are given by

Our reconstruction framework based on establishing the vector Y and the matrix Φ is schematically illustrated in Fig. 1. It is noteworthy that our framework can be extended to directed networks in a straightforward manner due to the feature that the neighbouring set of each node can be independently reconstructed. For instance, the neighbouring vector X can be defined to represent a unique link direction, for example, incoming links. Inference of the directed links of all nodes yields the full topology of the entire directed network.

Reconstructing networks and infection and recovery rates

To quantify the performance of our method in terms of the number of base strings (equations) for a variety of diffusion dynamics and network structures, we study the success rates for existent links (SRELs) and null connections (SRNCs), corresponding to non-zero and zero element values in the adjacency matrix, respectively. We impose the strict criterion that the network is regarded to have been fully reconstructed if and only if both success rates reach 100%. The sparsity of links makes it necessary to define SREL and SRNC separately. As the reconstruction method is implemented for each node in the network, we define SREL and SRNC on the basis of each individual node and the two success rates for the entire network are the respective averaged values over all nodes. We also consider the issue of trade-off in terms of the true positive rate (for correctly inferred links) and the false positive rate (for incorrectly inferred links).

Here we assume that there is no hidden source and the spreading process starts from a fraction of infected nodes, and record the binary time series. Figure 2a shows the reconstructed values of the components of the neighbouring vector X of all nodes. Let be the number of base strings normalized by the network size N. For small values of , for example, , the values of elements associated with links and that associated with null connections (actual zeros in the adjacency matrix) overlap, leading to ambiguities in the identification of links. In contrast, for larger values of , for example, , an apparent gap emerges between the two groups of element values, enabling us to correctly identify all links by simply setting a cut-off within the gap (see Supplementary Fig. 1a and Supplementary Note 4 for a method to set the cut-off). The success rates (SREL and SRNC) as a function of for SIS and CP on both homogeneous and heterogeneous networks are shown in Fig. 2b,c, where we observe nearly perfect reconstruction of links insofar as exceeds a relatively small value—an advantage of compressed sensing. The exact reconstruction is robust in the sense that a wide range of values can yield nearly 100% success rates. Our reconstruction method is then effective for tackling real networks in the absence of any a priori knowledge about its topology. In particular, the existence of a clear gap in the reconstructed vector X represents a successful reconstruction for a real network.

Note that a network is reconstructed through the union of all neighbourhoods, which may encounter ‘conflicts’ with respect to the presence/absence of a link between two nodes as generated by reconstruction centred at the two nodes. The conflicts would reduce the accuracy in the reconstruction of the entire network. To characterize the effects of edge conflicts, we study the consistency of mutual assessment of the presence or absence of link between each pair of nodes, as shown in Fig. 2b,c. We see that inconsistency arises for small values of but vanishes completely when the success rates reach 100%, indicating complete consistency among the mutual inferences of nodes and consequently guaranteeing accurate reconstruction of the entire network. Detailed results of success rates and trade-off measures with respect to a variety of model and real networks are displayed in Table 1, Supplementary Figs 2 and 3. and Supplementary Note 5.

Although the number of base strings is relatively small compared with the network size, we need a set of strings at different time with respect to a base string to formulate the mathematical framework for reconstruction. We study how the length of time series affects the accuracy of reconstruction. Figure 3a,b show the success rate as a function of the relative length nt of time series for SIS and CP dynamics on both homogeneous and heterogeneous networks, where nt is the total length of time series from the beginning of the spreading process divided by the network size N. The results demonstrate that even for very small values of nt, most links can already be identified, as reflected by the high values of the success rate shown. Figure 3c,d show the minimum length required to achieve at least 95% success rate for different network sizes. For both SIS and CP dynamics on different networks, decreases considerably as N is increased. This seemingly counterintuitive result is due to the fact that different base strings can share strings at different times to enable reconstruction. In general, as N is increased, will increase accordingly. However, a particular string can belong to different base strings with respect to the threshold Δ, accounting for the slight increase in the absolute length of the time series (see Supplementary Fig. 4 and Supplementary Note 5) and the reduction in (see Supplementary Note 3 on the method to choose base and subordinate strings). The dependence of the success rate on the average node degree ‹k› for SIS and CP on different networks has been investigated as well (see Supplementary Fig. 5 and Supplementary Note 5). The results in Figs 2 and 3, Supplementary Figs 2–5 and Table 1 demonstrate the high accuracy and efficiency of our reconstruction method based on small amounts of data.

In practice, noise is present and it is also common for time series from certain nodes to be missing, and it is necessary to test the applicability of our method in more realistic situations. Figure 4a,b show the dependence of the success rate on the fraction nf of states in the time series that flip due to noise for SIS and CP dynamics on two types of networks. We observe that the success rates are hardly affected, providing strong evidence for the applicability of our reconstruction method. For example, even when 25% of the nodal states flip, we can still achieve about 80% success rates for both dynamical processes and different network topologies. Figure 4c,d present the success rate versus the fraction nm of unobservable nodes, the states of which are externally inaccessible. We find that the high success rate remains mostly unchanged as nm is increased from 0 to 25%, a somewhat counterintuitive but striking result. The high degree of robustness against the limit to accessing nodal states is elaborated further in Supplementary Fig. 6 and Supplementary Note 5. We find that, in general, missing information can affect the reconstruction of the neighbouring vector, as reflected by the reduction of the gap between the reconstructed values associated with actual links and null connections. However, even for high values of nm, for example, nm=0.3, there is still a clear gap, indicating that a full recovery of all links is achievable. We have also found that our method is robust against inaccurately specified diffusion processes with fluctuation in infection rates (see Supplementary Fig. 7 and Supplementary Note 5). Taken together, the high accuracy, efficiency and robustness against noise, missing information and inaccurately modelling of real dynamical processes provide strong credence for the validity and power of our framework for binary time-series-based network reconstruction.

Having reconstructed the network structure, we can estimate the infection and recovery rates of individuals to uncover their diversity in immunity. This is an essential step to implement target vaccination strategy in a population or on a computer network to effectively suppress/prevent the spreading of virus at low cost, as a large body of literature indicates that knowledge about the network structure and individual characteristics is sufficient for controlling the spreading dynamics47484950. Here we offer an effective method to infer the individuals’ infection rates λi based solely on the binary time series of the nodal states after an outbreak of contamination. (To our knowledge, there was no prior work addressing this critical issue.) In particular, after all links have been successfully predicted, λi can be deduced from the infection probabilities that can be approximated by the corresponding infection frequencies (see Methods). These probabilities depend on both λi and the number of infected neighbours. The reproduced infection rates λi of individuals for both SIS and CP dynamics on different networks are in quite good agreement with the true values with small prediction errors (see Supplementary Fig. 8 and Supplementary Note 6). Results from a comprehensive error analysis are listed in Table 1, where the uniformly high accuracy validates our method. The inhomogeneous recovery rates δi of nodes can be predicted from the binary time series in a more straightforward way, because δi do not depend on the nodal connections (see Supplementary Fig. 9 and Supplementary Note 6). Thus, our framework is capable of predicting characteristics of nodal diversity in terms of degrees, and infection and recovery rates based solely on binary time series of nodal states.

Locating the hidden source of propagation

We assume that a hidden source exists outside the network but there are connections between it and some nodes in the network. In practice, the source can be modelled as a special node that is always infected. Starting from the neighbourhood of the source, the infection originates from the source and spreads all over the network. We collect a set of time series of the nodal states, except the hidden source (see Methods). The basic idea of ascertaining and locating the hidden source is based on missing information from the hidden source when attempting to reconstruct the network. In particular, to reconstruct the connections belonging to the immediate neighbourhood of the source accurately, time series from the source are needed to generate the matrix Φ and the vector Y. However, as the source is hidden, no time series from it are available, leading to reconstruction inaccuracy and, consequently, anomalies in the predicted link patterns of the neighbouring nodes. It is then possible to detect the neighbourhood of the hidden source by identifying any abnormal connection patterns51, which can be accomplished by using different data segments. If the inferred links of a node are stable with respect to different data segments, the node can be deemed to have no connection with the hidden source; otherwise, if the result of inferring a node’s links varies significantly with respect to different data segments, the node is likely to be connected to the hidden source. The s.d. of the predicted results with respect to different data segments can be used as a quantitative criterion for the anomaly. Once the neighbouring set of the source is determined, the source is then precisely located topologically.

Figure 5 presents an example, where a hidden source is connected with four nodes in the network (Fig. 5a), as reflected in the network adjacency matrix (Fig. 5b). We implement our reconstruction framework on each accessible node by using different sets of data in the time series. For each data set, we predict the neighbours of all nodes and generate an adjacency matrix. Averaging over the elements corresponding to each location in all the reconstructed adjacency matrices, we obtain Fig. 5c, in which each row corresponds to the mean number of links in a node’s neighbourhood. The inferred links of the immediate neighbours of the hidden source exhibit anomalies. To quantify the anomalies, we calculate the structural s.d. σ from different data segments, where σ associated with node i is defined through the ith row in the adjacency matrix as

where j denotes the column, represents the element value in the adjacency matrix inferred from the kth group of the data, is the mean value of aij and g is the number of data segments. Applying equation (10) to the reconstructed adjacency matrices gives the results in Fig. 5d, where the values of σ associated with the immediate neighbouring nodes of the hidden source are much larger than those from others (which are essentially zero). A cut-off value can be set in the distribution of σi to identify the immediate neighbours of the hidden source (see Supplementary Fig. 1b and Supplementary Note 4). The performance of locating hidden source by means of the trade-off measures (true positive rate versus false positive rate) are displayed in Table 1.

Discussion

We have developed a general framework to reconstruct complex propagation networks on which epidemic spreading takes place from binary time series. Our paradigm is based on compressed sensing, is completely data-driven and practically significant for controlling epidemic spreading through targeted vaccination. Both theoretically and practically, our framework can be used to address the extremely challenging problem of reconstructing the intrinsic interacting patterns of complex stochastic systems based on small amounts of polarized time series. The key to success of our method lies in our development of a novel class of transformation, allowing the network inference problem to be converted to the problem of sparse signal reconstruction, which can then be solved by the standard compressed-sensing algorithm. The accuracy and efficiency of our framework in uncovering the network structure, the natural diversity in the nodal characteristics, and any hidden source are guaranteed by the CST with rigorous proof for low-data requirement and convergence to optimal solution. The feasibility of our framework has been demonstrated using a large number of combinations of epidemic processes and network structures, where in all cases highly accurate reconstruction is achieved. Our approach opens up a new avenue towards fully addressing the inverse problem in complex stochastic systems in a highly efficient manner, a fundamental stepping stone towards understanding and controlling complex dynamical systems in general.

We have focused on two types of spreading dynamics, SIS and CP, where an infected individual can recover and becomes susceptible again. In this regard, even if an outbreak occurs, control strategy such as targeted vaccination or quarantine can be helpful to eliminate the virus eventually. A main purpose of our work is to identify the key individuals in the network to implement target control and to locate the source of infection to isolate it so as to prevent recurrent infection in the future. Although for any spreading dynamics, the most effective way to prevent a large-scale outbreak is to implement control during the early stage, this may be impractical in many situations. If we miss the early stage, which is possible especially in complex networks where the epidemic threshold can be near zero, to be able to reconstruct the spreading network is of tremendous value. Besides disease spreading, our framework is applicable to rumour or information spreading. In this case, identifying the source of rumour is important, a problem that our framework is capable of solving.

Our work raises a number of questions to further and perfect the theoretical and algorithmic development of reconstructing complex dynamical systems. For example, if partial knowledge about the network structure is available, the information can be incorporated into our framework to further reduce the required data amount. Moreover, for non-Markovian spreading processes, our current reconstruction framework may fail. This raises the need to develop new and more general approaches. Nevertheless, our theory, due to its generality and applicability to various types of inhomogeneous interactions, can be applied to networks of networks or interdependent networks, in which there may be different spreading patterns associated with distinct layers or components. Taken together, our results provide strong credence to the proposition that complex networks can be fully decrypted from measurements, even when stochastic disturbance and hidden sources are present. This can offer a deeper understanding of complex systems in general and significantly enhance our ability to control them based on, for example, the recently developed controllability theory of complex networks52535455565758.

Spreading processes

The SIS model is a classic epidemic model that has been used frequently to study a variety of spreading behaviours in social and computer networks. Each node of the network represents an individual and links are connections along which infection can propagate to others with certain probability. At each time step, a susceptible node i in state 0 is infected with rate λi if it is connected to an infected node in state 1. If i connects to more than one infected neighbour, the infection probability P01 is given by equation (4). At the same time, infected nodes are continuously recovered to be susceptible at the rates δi. The CP model has been used extensively to describe, for example, the spreading of infection and competition of animals over a territory, where λi is determined by equation (7). The main difference between SIS and CP dynamics lies in the influence on a node’s state from its vicinity. In both SIS and CP dynamics, λi and δi depend on the individuals’ immune systems and are selected from a uniform distribution characterizing the natural diversity (see Supplementary Note 7 for details of numerical simulations). Moreover, a hidden source is regarded as infected for all time.

Mathematical formulation of reconstruction based on CST

For SIS dynamics, suppose measurements at a sequence of times t=t1, t2, ···, tm are available. equation (5) leads to the following matrix form Ym × 1=Φm × (N−1)·X(N−1) × 1:

where the vector X(N−1) × 1 contains all possible connections between node i and all other nodes, and it is sparse for a general complex network. We see that if the vector Ym × 1 and the matrix Φm × (N−1) can be constructed from time series, X(N−1) × 1 can then be solved by using CST. The main challenge here is that the infection probabilities at different times are not given directly by the time series of the nodal states. To devise a heuristic method to estimate the probabilities, we assume that the neighbouring set Γi of the node i is known. The number of such neighbouring nodes is given by ki, the degree of node i and their states at time t can be denoted as

To approximate the infection probability, we use Si(t)=0 so that at t+1, the node i can be infected with certain probability. In contrast, if Si(t)=1, Si(t+1) is only related with the recovery probability δi. Hence, we focus on the Si(t)=0 case to derive . If we can find two time instants: t1, t2εT (T is the length of time series), such that Si(t1)=0 and Si(t2)=0, we can then calculate the normalized Hamming distance between and , where the normalized Hamming distance between two strings of equal length is defined as ratio of the number of positions with different symbols between them and the length of string. If , we can regard the states at the next time step, Si(t1+1) and Si(t2+1), as i.i.d Bernoulli trials. In this case, using the law of large numbers, we have

A more intuitive understanding of (equation (12)) is that if the states of i’s neighbours are unchanged, the fraction of times of i being infected by its neighbours over the entire time period will approach the actual infection probability . Note, however, that the neighbouring set of i is unknown and to be inferred. A strategy is then to artificially enlarge the neighbouring set to include all nodes in the network, except i. In particular, we denote

If H[S−i(t1), S−i(t2)]=0, the condition will be ensured. Consequently, due to the nature of i.i.d Bernoulli trials, from the law of large numbers, we have

Hence, the infection probability of a node at can be evaluated by averaging over its states associated with zero-normalized Hamming distance between the strings of other nodes at some time associated with . In practice, to find two strings with absolute zero-normalized Hamming distance is unlikely. We thus set a threshold Δ so as to pick the suitable strings to approximate the law of large numbers, that is

where serves as a base for comparison with S−i(t) at all other times and . As is not exactly zero, there is a small difference between and . We thus consider the average of for all tν to obtain , leading to the right-hand side of equation (14). We denote and . To reduce the error in the estimation, we implement the average on S−i(t) over all selected strings through equation (14). The averaging process is with respect to the nodal states Sj, j,≠i(t) on the right-hand side of the modified dynamical equation (5). Specifically, averaging over time t restricted by equation (14) on both sides of equation (5), we obtain . If λi is small with insignificant fluctuations, we can approximately have (see Supplementary Fig. 10 and Supplementary Note 8), which leads to . Substituting by , we finally get

Although the above procedure yields an equation that bridges the links of an arbitrary node i with the observable states of the nodes, a single equation does not contain sufficient structural information about the network. Our second step is then to derive a sufficient number of linearly independent equations required by CST to reconstruct the local connection structure. To achieve this, we choose a series of base strings at a number of time instants from a set denoted by Tbase, in which each pair of strings satisfy

where and correspond to the time instants of two base strings in the time series and Θ is a threshold. For each string, we repeat the process of establishing the relationship between the nodal states and connections, leading to a set of equations at different values of in equation (15), as described in the matrix form (equation (6)). See Supplementary Fig. 11, 12 and Supplementary Note 8 for the dependence of success rate on threshold Δ and Θ for SIS and CP dynamics in combination with four types of networks.

Inferring inhomogeneous infection rates

The values of the infection rate λi of nodes can be inferred after the neighbourhood of each node has been successfully reconstructed. The idea roots in the fact that the infection probability of a node approximated by the frequency of being infected calculated from time series is determined both by its infection rate and by the number of infected nodes in its neighbourhood. To provide an intuitive picture, we consider the following simple scenario in which the number of infected neighbours of node i does not change with time. In this case, the probability of i being infected at each time step is fixed. We can thus count the frequency of the 01 and 00 pairs embedded in the time series of i. The ratio of the number of 01 pairs over the total number of 01 and 00 pairs gives approximately the infection probability. The infection rate can then be calculated by using equations (4) and (7) for the SIS and CP dynamics, respectively. In a real-world situation, however, the number of infected neighbours varies with time. The time-varying factor can be taken into account by sorting out the time instants corresponding to different numbers of the infected neighbours, and the infection probability can be obtained at the corresponding time instants, leading to a set of values for the infection rate whose average represents an accurate estimate of the true infection rate for each node.

To be concrete, considering all the time instants tν associated with kI infected neighbors, we denote , ∀ tν, and Si(tν)=0, where Γi is the neighbouring set of node i, kI is the number of infected neighbours and represents the average infected fraction of node i with kI infected neighbours. Given , we can rewrite equation (4) by substituting for and for λi, which yields . To reduce the estimation error, we average with respect to different values of kI, as follows:

where Λi denotes the set of all possible infected neighbours during the epidemic process and denotes the number of different values of kI in the set. Analogously, for CP, we can evaluate from equation (7) by

where is the node degree of i. Insofar as all the links of i have been successfully reconstructed, can be obtained from the time series in terms of the satisfied Si(tν+1), allowing us to infer via equations (17) and (18).

Note that the method is applicable to any type of networks insofar as the network structure has been successfully reconstructed.

Networks analysed

Model networks and real networks we used are described in Supplementary Note 10 and Table 1.

Author contributions

W.-X.W., Z.S., Y.F., Z.D. and Y.-C.L. designed research; Z.S. and W.-X.W. performed research; Y.F. and Z.D. contributed analytic tools; Z.S., W.-X.W., Y.F., Z.D. and Y.-C.L. analysed data; and W.-X.W. and Y.-C.L. wrote the paper.

Additional information

How to cite this article: Shen, Z. et al. Reconstructing propagation networks with natural diversity and identifying hidden sources. Nat. Commun. 5:4323 doi: 10.1038/ncomms5323 (2014).