Inferring Network Connectivity by Delayed Feedback Control

We suggest a control based approach to topology estimation of networks with elements. This method first drives the network to steady states by a delayed feedback control; then performs structural perturbations for shifting the steady states times; and finally infers the connection topology from the steady states' shifts by matrix inverse algorithm () or -norm convex optimization strategy applicable to estimate the topology of sparse networks from perturbations. We discuss as well some aspects important for applications, such as the topology reconstruction quality and error sources, advantages and disadvantages of the suggested method, and the influence of (control) perturbations, inhomegenity, sparsity, coupling functions, and measurement noise. Some examples of networks with Chua's oscillators are presented to illustrate the reliability of the suggested technique.


Introduction
The research on complex networks [1][2][3][4] pervades almost all biological sciences, from gene network [5,6] to system biology [7], from physiology [8][9][10] to psychology [11], to name just a few. Recent developments [12] in the quantitative analysis of complex networks, based largely on graph theory, have been rapidly translated to studies of brain networks. Mathematically, brain networks [12] can be described as graphs that are composed of nodes (vertices) denoting neural elements (neurons or brain regions) that are linked by edges representing physical connections (synapses or axonal projections) or functional ones based on imaging data. Current studies of brain networks focus on understanding the relation between network connectivity and function [12]. It turns out that small perturbations of structural and functional connectivity may dramatically change the function of networks and even lead to the occurrence of cognitive dysfunctions. In the context of brain functional networks based on imaging data [12], for example, one may quantify the functional connectivity between brain regions by analyzing the topological parameters (such as clustering coefficient, connectivity distribution, and average network distance) of the functional network, and the change of the topological properties has been considered as the pathophysiological mechanism of cognitive dysfunctions. In order to infer the emergent function of a real network, one first has to identify the underlying (functional and structural) connection topology.
The Pearson's correlation method [13][14][15] is based on the following assumption: If the value of Pearson's correlation coefficient between two brain imaging time-series, representing the activities of two brain regions of interest, exceeds a threshold, then there exists a linkage between the two brain regions; otherwise, there is no connection between them. However, how to determine suitable thresholds is still an open problem and the assumption that correlation implies connections (or causality) is logically not sound [30,31]. This problem also occurs with the phase synchronization approach [16] that depends on the following assumption: If the phase synchronization degree (or index) between two brain imaging time-series is above a threshold, then there exists a linkage between the two brain regions; otherwise, there is no connection between them. Again, the determination of the threshold remains a nontrivial problem. Furthermore, how to define the phase of complex systems still remains an open problem.
Bayesian estimation methods [17,18] have been used to evaluate the connectivity between brain regions of interest with imaging data, but their efficiency and feasibility depend on the validity of the priors and the model adopted.
Network topology estimation using identical synchronization (which is conceptually equivalent to adaptive observer) was first developed in Ref. [19]. However, synchronization of networks may become an obstacle of topology estimation because synchronization leads to a situation where network connectivity information is hidden. Therefore, one has to complete the estimation processing as soon as possible (before the network is synchronous), otherwise one requires proper external perturbations to shift the network out of synchronous state.
Perturbation based method [20,21] transforms the topology estimation problem into a matrix inversion task. It has been shown [21] that for sparsely connected networks, this matrix-inverse problem can be solved effectively using an ' 1 -norm optimization strategy in combination with the well-known singular value decomposition technique. The perturbation method [20,21], however, depends on the steady-state assumption (more precisely, it is assumed that the network to be analyzed always reaches a stable stationary state automatically) which is a restriction for some network systems with complex dynamical behaviors (including chaos). When the external perturbation matrix is unknown, a recursive strategy [22] can be used to estimate both the perturbation and connectivity matrices.
Some authors [23,24] recently developed a so-called compressive sensing method that first formulates the dynamical system of interest as the following equation with F(X) being a polynomial function and P being parameter vector to be estimated, then obtains two data matrices Y1 _ X X(t 1 ), _ X X(t 2 ), . . . , _ X X(t m ) and Z~½F(X(t 1 )),F(X(t 2 )), . . . ,F(X(t m )), satisfying Y~ZP, and finally estimate P by an ' 1 -norm convex optimization processing. They showed that their method is effective to reconstruct dynamical systems [23] and network topology [24]. However, such a method requires a differential estimator that may be sensitive to measurement noise. Furthermore, complex dynamical systems usually cannot be described by Eq. (1), more precisely, their dynamical equation in general is non-polynomial and does not linearly depend on the parameters. It should be remarked that the performance of ' 1 -norm convex optimization strategy usually becomes bad when the sparsity of networks decreases, as will be shown below (cf. Fig. 8).
Timme's recent work [26] analyzed the possibility of direct topology reconstruction from dynamical trajectories. Remarkably, the question how the parameters (e.g. sampling rate, observation time, and external noise) influence the performance of topology reconstruction is discussed in detail. The reliability of his method [26] has been demonstrated clearly. As a minor drawback, his method requires some prior knowledge about local dynamics of each node, and a differential estimator that may be sensitive to measurement noise.
To use the perturbation method also for networks with complex dynamics, a linear state feedback control based method [27][28][29] was suggested very recently, and can be used to estimate topology by exploiting information obtained from the observed steady-state responses of each node. However this method has some limitations. For instance, one generally has to assume that all state variables of each node are completely measurable and all state components of each node admit an external input. Furthermore, a high-gain feedback control will be involved in some cases.
In brief, most of developed topology estimation methods have their advantages and disadvantages, and thus far the topology estimation issue remains an open problem. Here we make an effort to remove some drawbacks of previous methods, and show that the connection topology of complex dynamical networks can be identified by exploiting information obtained from shifted steady states that are stabilized by means of multiple delay feedback control (MDFC) [32]. This control approach is combined with some methods [21] for detecting connectivity of networks under the assumption that a stable stationary state exists (also called steady state assumption). However, in contrast to that work, our topology detection method is applicable to dynamical networks with complex dynamical behaviors (far from stationarity) and does not depend on the steady state assumption. Furthermore, our method is possible to be applied in a challenging scenario where only one state variables of each node are measurable and accessible, and does require only a little structure information about the networks under study.

Theory
We consider a network of interacting dynamical systems, given by where . . ,f in T : R n ?R n describes the dynamics of the ith element. For simplicity we assume that only the first components of each element are connected to each other (a more general case will be treated elsewhere). Here h ij : R?R is a coupling function and C~½1,0, . . . ,0 T . The topology of the network connections is determined by the adjacency matrix A~(a ij ): a ij~1 if there exists a connection from the jth node to the ith node; and a ij~0 otherwise. We shall show that MDFC [32] is very efficient to shift the steady states and the steady states' shifts enable us to uncover the connection topology in terms of an estimation of the elements of the matrix A~(a ij ).
We restrict ourselves to the case that only the coupling variables, namely x i , can be measured (or monitored) and we add the control term to only the first equation of each element, where delay times t 1 and t 2 and control gains k 1 and k 2 are uniform for all elements. For D i~0 , the control signal (3) becomes the original MDFC [32].
If system (2) without any perturbation has at least one equilibrium, which usually is satisfied for most of networks, then equation F(X)~0 has at least one real root. By using the continuity of function F (because f i is continuous for all i), it follows that when constants i1 and i2 are close to zero, equation F(X)~D has at least one real root. This indicates that Assumption 1 is not really a restriction in practice.
The following theorem is the foundation of topology identification and provides conditions under which the system (4) is locally asymptotical stable at a stationary state. Detailed discussion about Theorem 1 can be found in Discussion Part.
Theorem 1: The system (4) (with D i [ ½ i1 , i2 for all i) is locally asymptotical stable at a stationary state S, satisfying where Eq. (5) has been used. Locally linearizing the above system around the origin results in Therefore, in terms of the standard linear system theory, the stability of the error system (7) determines by the characteristics equation det½lI{DF(S)zk 1 0. If all roots of the characteristics equation have negative real parts, then the asymptotic stability of the error system (7) is satisfied. This completes the proof.
Steady-state shifts. If proper k i , t i , and D i are chosen such that Theorem 1 is fulfilled (see Discussion Part for further information), then one can stabilize the steady state (x 1s , . . . ,x ns ), satisfying Vi where . . ,g in T is nonsingular, then one can conclude from the implicit function theory [33] that there exists a mapping w i : R?R n{1 such that Substituted into the first equation of Eq. (8) this yields where g i (x) :~f i1 (x,w i (x)).
As will be shown below, Eq. (10) is the foundation of the topology estimation method to be suggested, and has reduced the original n-dimensional problem to an 1-dimensional (1D) one. It should be remarked that Eq. (10) is satisfied, provided that (i) Equation (8) has at least one real solution; (ii) the steady state satisfying Eq. (8) can be stabilized by MDFC; and (iii) We now show that shifting the steady states of the network system M times by M structural perturbations enables us to uncover the network connectivity (where M depends on the network size N).
For the mth perturbation, we replace the control constant D i by D i zdD m i for each node i such that the steady states of the coupling variables are shifted from g i1 to g i1 zdg m i1 for all i. Then the resulting steady state response equations of the coupling variables read For sufficiently small perturbations dD m i , we approximate Then the set of equations (12) can be rewritten in a compact form which contains N equations that restrict the N 2 elements b ij . Perturbing the steady state of the network system M times, we achieve where To summarize the above analysis, Eq. (15) is fulfilled if and only if: (i) Equation (8) has at least one real solution; (ii) the steady states satisfying Eqs. (10) and (11) can be stabilized by MDFC; . . ,g in T is nonsingular; and (iv) perturbations dD m i are sufficiently small for all m.
Topology estimation using matrix inverse algorithm. Equation (15) actually contains NM conditions that restrict the N 2 elements b ij . Hence, after performing M~N perturbations, all elements b ij can be estimated byb b ij , given bŷ if the inverse of L N exists. It follows that if all elements b ij can be estimated with high accuracy (more precisely, there exists a sufficiently small such that jb ij {b b ij jƒ) , then all off-diagonal elements a ij can be identified from Eq.(13): a ij~0 when jb b ij jƒ; and a ij~1 otherwise. In practice, one may follow the SDTIA algorithm [28] and divide all values jb b ij j into two sets: I 0 containing all elements jb b ij j corresponding to a ij~0 and I 1 containing all elements jb b ij j corresponding to a ij~1 , by the following steps: Step 1. Calculate elements jb b ij j for all i,j.
Step 2. Order (or arrange) all elements jb b ij j in an ascending sequence and obtain a new sequence fs i g.
Step 3. The critical point sequence number i c of set I 0 is determined by the rule: As clearly shown in Fig. 1 that when w 1 w3 with w 1~m in i,j,i=j jb ij j and w 2~m ax i,j,i=j jb ij j, the distance between sets I 0 and I 1 is larger than the length of set I 0 , and thereby one can distinguish the sets I 0 and I 1 by the above steps (SDTIA algorithm [28]) and reconstruct the network topology in terms of an estimation of all elements of the matrix A~(a ij ), where the distance between two point sets is equal to the minimal distance between any two points which are taken from different sets, and the length of a point set is the difference between the maximal and minimal values in the set. Therefore, the smaller the value of and the bigger the value of w 1 , the higher the possibility of successful topology reconstruction.
Topology estimation using ' 1 -norm optimization strategy. Topology estimation using Eq. (16) requires N perturbations and becomes ''costly'' and less effective when the network size N is very large. However, for sparsely connected networks, it turns out that by using a ' 1 -norm convex optimization strategy to be shown below, we can accurately and efficiently approximate all elements b ij from Eq. (15) with M%N.
We transpose Eq. (15) and rewrite it as where p j and q j are the jth column vector of matrices B T and V T M , respectively.
The estimated value of each p j , referred to asp p j , can be determined by solving the following convex optimization problem where is the tolerance (in the following simulations, ~10 {5 ), ExE 1~Pi jx i j is the ' 1 -norm of vector x and ExE 2 2~Pi x 2 i . The advantage of choosing the formulation (18) is that one can determine the network with a minimal number of connections (each vectorp p j will have a minimal number of nonzero elements) and it can be solved in polynomial time with some standard scientific softwares (e.g., Matlab toolbox CVX Ver1.1 [34]). By this ' 1 -norm convex optimization strategy, we can determine the matrix with minimal driving connections for each node; hence we can effectively estimate all b ij for sparsely connected networks when M%N perturbations are performed, as will be illustrated below. Again, one may follow the SDTIA algorithm [28] (shown above) for an effective topology reconstruction.
Topology estimation quality. Following Timme's work [26], we define the normalized estimation error e ij of each element b ij by whereb b ij is an estimation of b ij , and b max~m ax i;j fjb b ij j; jb ij jg.
We further define the estimation accuracy a [ ½0, 1 such that b ij can be identified correctly if This implies from Fig. 1 that ~2b max (1{a) and thereby the topology can be estimated correctly if where jmax i,j,i=j jb ij j{w 2 jƒ2b max (1{a) and jw 2 jƒb max are used. Therefore, the bigger the values of g and a, the higher the topology estimation accuracy. Based on the condition (21), the minimal value of g being supported for a successful topology reconstruction is determined by the maximal value of a satisfying the condition (20). The estimation accuracy of b ij is crucial for topology reconstruction, so it is of importance to quantify the estimation quality of values b ij . Here we qualify the estimation accuracy of all non-diagonal elements b ij as a whole by the variable Q a , given by where H is the Herviside step function, i.e., H(x)~1 for x §0 and H(x)~0 otherwise. This definition is a little bit different from Timme's work [26] that considered the estimation of all elements b ij . It is clear that the bigger the values of a and Q a , the higher the estimation accuracy of all non-diagonal elements b ij . Based on this observation, we restrict ourselves and assume that an effective network topology reconstruction is said to occur when Q 0:98 §0:99.

Simulation
To illustrate the above topology estimation methods, we use a network of coupled Chua's circuits, given by  Figure 1. The condition to ensure a successful topology reconstruction using the SDTIA algorithm [28]. Sets I 0 and I 1 contain all elements jb b ij j corresponding to a ij~0 and that corresponding to a ij~1 , respectively, where w 1~m in i,j,i=j jb ij j and w 2~m ax i,j,i=j jb ij j. doi:10.1371/journal.pone.0024333.g001 Fig. 2. In the following, we first illustrate the steady state stabilization and shifts numerically. Then, based on steady state shifts and measurement, we show two methods for topology estimation, i.e., matrix inverse and ' 1 -norm convex optimization strategy, with estimation accuracy being quantified by Q a . Following Ref. [32], we can determine suitable control parameter values k i and t i by a search strategy. We numerically found that there is a big window for the control parameters k i and t i such that system (2) can be driven to a steady state by the MDFC (3), as illustrated in Fig. 3 as a typical example. It is clear from Fig. 3D that MDFC is very effective for steady state stabilization. Furthermore, when MDFCs with distinct D i are used, the steady-state response shift phenomenon can be observed (cf. Fig. 4 for a representative result). Those steady state shifts are the foundation of topology estimation, as shown above (cf. Theory Part).
When system (23) is driven to a steady state (x 1s , . . . ,x ns ) with x is~½ g i1 ,g i2 ,g 3n T being the steady state of the ith element, then one can easily confirm that is nonsingular and where g i2~n (g i1 ) is the unique solution of equation g i2 z l(g i2 )~g i1 . Therefore, Eq. (10) is fulfilled. This implies that shifting and measuring the steady state response of the first state of each node becomes possible for a successful topology reconstruction. In the following, we show two methods for topology estimation, i.e., matrix inverse and ' 1 -norm convex optimization strategy.
As a representative result using the matrix inverse algorithm (16) for topology estimation, Fig. 5A shows the estimation error e ij of elements b ij for a random directed network of interacting Chua's oscillators. It is clear that all elements b ij have been reconstructed effectively with Q 0:98~1 (due to e ij v0:02 for all i,j). With this high (normalized) estimation accuracy of b ij , one may identify all parameters a ij correctly by the SDTIA algorithm [28] (also shown above), as illustrated in Figs. 5C-5D where the estimated jb ij j (with %) corresponding to a ij~1 are bigger than that (with D) corresponding to a ij~0 .
The matrix inverse method for topology reconstruction requires N perturbations and becomes ''costly'' and less effective when the network size N is very large. However, such a drawback for sparsely connected networks may be relaxed by the ' 1 -norm convex optimization strategy described in Eq. (18). As typical examples, Fig. 5B and Fig. 6B shows that an acceptable topology estimation accuracy (i.e., Q 0:98~0 :9975, Q 0:98~0 :9954) can be obtained when only M%N perturbations are performed.
Furthermore, the matrix inverse method may lead to wrong conclusion in some cases due to the ill-condition problem of the matrix inverse, as illustrated in Fig. 6A where Q 0:98~0 :2433, implying a bad estimation result, is achieved. However, for sparse networks, such a drawback may be removed by the ' 1 -norm convex optimization strategy, as shown in Fig. 6B where Q 0:98~0 :9954.
The question how to choose control parameters becomes crucial for steady state shifts which are the foundation for topology reconstruction. For simplicity, in the above simulation, we let D i~0 and choose parameters dD i [ ½{0:3,0:3 randomly. We now analyze the influence of perturbation parameters dD i on topology estimation. Figure 7 summarizes our results and shows that the estimation accuracy Q 0:98 using the ' 1 -norm convex optimization strategy changes with the node-pair connection possibility p [ f0:1,0:2,0:3, . . . ,1g for two cases, i.e, undirected (cf. yellow bars) and directed (cf. red bars) networks. There, each bar represents the result of averaging over 30 random perturbations (with dD i being uniformly distributed in the range ½{0:3, 0:3) and the standard square error is given as well. From Fig. 7 we may draw the following conclusions: (i) the performance of topology reconstruction using the ' 1 -norm optimization strategy becomes bad when p increases; (ii) The estimation accuracy Q 0:98 is not sensitive to the choice of perturbation parameters dD i when Q 0:98 is close to one; (iii) There is no distinct difference between undirected (cf. yellow bars) and directed (cf. red bars) networks. This indicates that the performance of topology reconstruction using the ' 1 -norm convex optimization strategy is not sensitive to the inhomegeneity but sparsity (cf. Fig. 8).
Note that the ' 1 -norm convex optimization strategy is very effective for sparsely connected networks only. Hence, for nonsparsely connected networks, this optimization method usually has to require that almost all nodes are perturbed, as illustrated in Fig. 9A as a representative result. In this case, the ' 1 -norm convex optimization strategy has no any clear advantage compared to the matrix inverse algorithm (cf. Fig. 9B) that uses an ' 2 -norm optimization processing.
As mentioned above, we restrict ourselves and assume that an effective network topology reconstruction is said to occur when Q 0:98 §0:99. Based on this rule, we now analyze numerically the relation between the minimal number of perturbations, referred to as M min , that are required for a successful topology reconstruction satisfying Q 0:98 §0:99, and the network size N. Figure 10 summarizes our results and shows the logarithmic-linear plot of the relation of N and M min for two cases, i.e., 4-nearest-neighbor coupled network and directed network of nodes randomly connected with possibility p~0:1. There is a clear logarithmic-linear relation between N and M min . This result is consistent with Timme's work [21], and implies that we need less control applications (perturbations) than the size of the networks under study.
Measurement noise cannot be avoided in some cases and usually deteriorates the control performance of high-gain control methods because measurement noise is largely amplified. Fortunately, the MDFC method does not belong to high-gain control [28] and can stabilize stationary states with very small gain (indeed k 1~k2~1 was used in all simulation results in this paper). This implies that our topology estimation method is applicable to network systems in the presence of measurement noise, as illustrated in Fig. 11A where results are shown obtained from observed signals contaminated with 5% measurement noise. We found that more perturbations are generally required in the presence of measurement noise (cf. Fig. 11A where M~75, Q 0:98~0 :9935) compared to the case in the absence of measurement noise (cf. Fig. 11B where M~70, Q 0:98~0 :9967).
Finally, we analyze the influence of g on topology estimation, and revisit the network (23) but assume h ij (x)~w ij sin(x) with w ij being uniformly distributed in range ½w 1 , w 2 such that the value of g can be changed with the choice of parameters w 1 and w 2 . Figures 12 and 13 summarize our results and show in both cases (i.e., w 1~0 :01,w 2~2 and w 1~0 :001,w 2~2 ) that the minimal value of estimated jb ij j corresponding to a ij~1 is more than twice the maximal value of that corresponding to a ij~0 , and thereby one may identify all parameters a ij correctly by the SDTIA algorithm [28]. Furthermore, the ratio of the distance between sets I 0 and I 1 to the maximal value of set I 0  roughly increases with the value of g where the definition of sets I 0 and I 1 is illustrated in Fig. 1. Therefore, there exists a critical value g c such that if gwg c is fulfilled, then one may identify all elements a ij correctly. On the other hand, when gvg c , the boundary between sets I 0 and I 1 will become unclear and some elements a ij cannot be identified correctly. Even under such a circumstance, it is still possible to estimate partial elements a ij correctly if a suitable strategy is used to delete those elements jb b ij j contaminating the boundary between sets I 0 and I 1 . Detailed analysis is now under our investigation and will be reported elsewhere. and that corresponding to a ij~0 after being sorted with ascending order, respectively. It is clear from Panels (C)-(D) that one may identify all parameters a ij correctly by the SDTIA algorithm [28]. doi:10.1371/journal.pone.0024333.g005

Delayed feedback control
It has been shown experimentally [35][36][37][38] that Pyragas's delayed feedback control method [39], which feeds the amplified difference of a monitor (or measurable) variable and its delayed component back to the controlled system, is applicable and very effective to stabilize unstable period orbits as well as unstable equilibrium points. Some advantages of Pyragas's delayed feedback control method include: (i) it feeds the amplified difference of a monitor (or measurable) variable and its delayed component back to the controlled system but does not use any structure information about the controlled system; (ii) it is noninvasive, that is, the control signal approaches to zero after a unstable period orbit or a unstable equilibrium point is stabilized; and (iii) it can easily be realized using analog or digital devices. Some extended versions using more delayed components have also been developed for improving further the control performance, such as extended time delay auto synchronization [40,41] and multiple delay feedback control [32,42,43] methods.
Although the reliability of all delayed feedback control methods for stabilizing unstable period orbits and unstable equilibrium points has been illustrated by various experiments, the theoretical analysis and mechanism of delayed feedback control is still far from strictness and completeness [44][45][46][47]. Fortunately, we found that the steady state stabilization based on MDFC is always possible for a large class of dynamical networks. In practice, one can usually determine suitable control parameter values by a search strategy, as illustrated in previous work [32].
Thus far the research on delayed feedback control focused on stabilizing unstable period orbits and unstable equilibrium points of chaotic systems. In this paper, we show a potential application of using delayed feedback control for topology reconstruction. Compared to previous linear state feedback control method [27][28][29] which in general requires high-gain control and full state feedback (i.e., all state components of each node are measurable and accessible), the suggested delayed feedback control method is applicable even in a challenging scenario where only one state variables of each node are measurable and accessible.

Extension to more general coupling functions
Our method can also be extended to networks with more general coupling functions but does not limit to those with only the state-difference form h ij (x j {x i ). To demonstrate this point more clearly, we consider the following network Figure 8. The influence of the sparsity on topology reconstruction. The estimation error Q 0:98 changes with the sparsity of directed random networks. There, the sparsity is defined as the ratio of the number of zero non-diagonal elements a ij to N(N{1), and Q 0:98 is calculated using the ' 1 -norm convex optimization strategy (with N~64,M~60). Furthermore, each black point represents the result of averaging over 30 random perturbations (with D i~0 and dD i being uniformly distributed in the range ½{0:3, 0:3) and the standard square error is given as well. doi:10.1371/journal.pone.0024333.g008 where all variables follow the same definition in system (2) except the coupling functions h ij . Here h ij : R|R?R. Again, we assume that system (26) can be driven to a steady state by the control signal (3). In this case, following similar steps developed for the state-difference form, one can easily see that Eq. (12) now reads where the first order approximation h ij (x j , This implies that Eq. (14) is again fulfilled but the matrix a ij (Lh ij =Lg i1 ), for i~j: Therefore, our methods using matrix inverse algorithm and ' 1norm convex optimization strategy can be extended to topology reconstruction of network (26) with more general coupling form, as illustrated in Fig. 14 where h ij (x j ,x i )~sin(x j ){sin(x i ) and the network topology can be estimated effectively.

Implementation and error sources
We briefly outline our method for topology estimation: i.
Drive the network (with N nodes) to a steady state by control signal (3) with D i (usually D i~0 ), and measure the resulting steady state response g i1 for all i; ii. Perturb the control signal (3) (i.e., replace D i by D i zdD m i where dD m i is randomly chosen from the range [-v, v]) M times, and measure the resulting steady state response g i1 zdg m i1 for all i; iii. Estimate all non-diagonal elements b ij using the matrix inverse algorithm (M~N) or the ' 1 -norm convex optimization strategy (MƒN); iv. Infer all non-diagonal elements a ij from estimated b ij by the SDTIA algorithm [28].
One may see from the above steps that the topology estimation error may come from different sources: (i) Steady state control; (ii) Steady state measurement; (iii) The first order approximation concerning functions h ij (x); (iv) The matrix inverse operation error (for the matrix inverse algorithm) or the optimization error (for the ' 1 -norm convex optimization strategy); and (v) The value of g.
As described above, delayed feedback control methods [35][36][37][38] are very efficient for stabilizing stationary states in various real systems such as optics, semiconductors, networks of chemical oscillators, and reaction-diffusion systems. Therefore, steady state control usually cannot be considered as an error source, as illustrated in Figs. 5, 6B, 9B, 10, 11, and 12.
Measurement of steady states also cannot be taken as a major error source, as illustrated in Fig. 11 where acceptable results are  For the matrix inverse algorithm, a major error source may come from the inverse operation itself, as illustrated in Fig. 6A where a bad estimation result (with Q 0:98~0 :2433) is achieved due to the ill-condition problem of the matrix inverse operation.
For the ' 1 -norm convex optimization strategy, a major error source may come from the sparsity of networks under study, as illustrated in Fig. 8. This is consistent with the fact that the ' 1 -norm convex optimization strategy is effective for sparsely connected networks only.
The influence of g on topology estimation has been illustrated in Figs. 12-13. It is clear that the ratio of the distance between sets I 0 and I 1 to the maximal value of set I 0 roughly increases with the value of g where the definition of sets I 0 and I 1 is illustrated in Fig. 1. Therefore, there exists a critical value g c such that if gwg c is fulfilled, then one may identify all elements a ij correctly. It should be remarked that the value of g c is determined by the control signal (3), the coupling functions, the equilibria of network (2), and the initial states. If the network under study has more than one equilibrium, then it is still possible to change the value of g c by choosing the proper time to perform the steady state control to shift the equilibrium of the network dramatically. However, such a strategy in principle has to Figure 13. Topology reconstruction in the case of w 1~1 |10 {3 , w 2~2 , and g~6:776|10 {4 . (A) The estimation error surface of a directed network with N~100 and node-pair connection probability p~0:1 is calculated using ' 1 -norm convex optimization strategy with M~80. With the normalized estimation errors e ij shown in Panel (A), Panels (B)-(C) plot the estimated jb ij j corresponding to a ij~1 and that corresponding to a ij~0 after being sorted with ascending order, respectively. Insert in Panel (b) shows a local augment. It is clear that the minimal value of estimated jb ij j shown in Panel (b) is more than twice the maximal value of estimated jb ij j shown in Panel (c), and thereby one may identify all parameters a ij correctly by the SDTIA algorithm [28]. doi:10.1371/journal.pone.0024333.g013 The estimation error surface of a directed network with N~100 and node-pair connection probability p~0:1 is calculated using ' 1 -norm convex optimization strategy with M~80. With the normalized estimation errors e ij shown in Panel (A), Panels (B)-(C) plot the estimated jb ij j corresponding to a ij~1 and that corresponding to a ij~0 after being sorted with ascending order, respectively. Insert in Panel (B) shows a local augment. It is clear that the minimal value of estimated jb ij j shown in Panel (b) is more than twice the maximal value of estimated jb ij j shown in Panel (C), and thereby one may identify all parameters a ij correctly by the SDTIA algorithm [28]. doi:10.1371/journal.pone.0024333.g012 require some prior knowledge about the the equilibria of the network, and thereby has its restriction in some applications.

Advantages and disadvantages of our method
Some advantages of our method include: i.
If network synchronization occurs and leads to vanishing coupling terms, then the network connectivity information is hidden and cannot be recovered with time-series analysis methods [19,24,26]. However, our topology reconstruction method is applicable to synchronous networks; ii. Previous topology reconstruction method [27,28] based on steady-state stabilization generally has to assume that all state variables of each node are completely measurable and all state components of each node admit an external input. However, our method is applicable even in a challenging scenario where only one state variables of each node are measurable and accessible; iii. Our method requires only small control injection and does not belong to a kind of high-gain control [27,28]. Hence it is not sensitive to measurement noise and may achieve better performance than high-gain control method [27,28] and the methods using differential estimator in the presence of measurement noise; iv. Previous time-series methods [19,26] require a lot of information about the local dynamics of each node and coupling functions. This is really a restriction in some applications. However, our method does require only a little structure information about the controlled networks, and provides a promising solution for topology reconstruction if the required control perturbations are allowed.
On the other hand, our method also possesses some disadvantages: i.
Our method is applicable to topology estimation of sparsely connected networks with size N when M%N perturbations are performed, but in general one has to measure the steady state response of all nodes and the measurement ''cost'' increases linearly with the size of networks, even when only partial connections of interest require to be estimated. Such a drawback also exists for most of previous methods except the high-gain control method [28]; ii. Steady state stabilization and shifts are the foundation of our method. However, such a kind of steady state control (or perturbation) will influence the dynamical behavior of systems, so our method may fail for systems that do not support the required steady state control. In this case, previous time-series methods [19,24,26] might be considered as a potential strategy for topology reconstruction. iii. Our method may in principle fail when time-varying topology is required to be reconstructed. In such a circumstance, previous time-series methods [19,24,26] might be applicable for correct estimation.

Potential applications
Previous works have shown the importance of topology connections on spatiotemporal pattern of networks of coupled chemical oscillators [48][49][50][51]. Furthermore, delayed feedback control has effectively been applied to stabilize (unstable) steady states of chemical oscillators (cf. Ref. [52] for a representative result). Therefore, our method is possible to be used to reconstruct the connection topology of interacting chemical oscillators. Another possible application is to reconstruct topology of gene networks [22] by delayed feedback control, provided online measurement and injection techniques are feasible. Generally, the suggested technique enables us to identify the connection topology of real networks (including circuit networks and interacting coupled chemical oscillators [48][49][50][51]) which allow the required control applications (perturbations). Some possible experimental research is now under our investigation.