Inferring Network Topology from Complex Dynamics

Inferring network topology from dynamical observations is a fundamental problem pervading research on complex systems. Here, we present a simple, direct method to infer the structural connection topology of a network, given an observation of one collective dynamical trajectory. The general theoretical framework is applicable to arbitrary network dynamical systems described by ordinary differential equations. No interference (external driving) is required and the type of dynamics is not restricted in any way. In particular, the observed dynamics may be arbitrarily complex; stationary, invariant or transient; synchronous or asynchronous and chaotic or periodic. Presupposing a knowledge of the functional form of the dynamical units and of the coupling functions between them, we present an analytical solution to the inverse problem of finding the network topology. Robust reconstruction is achieved in any sufficiently long generic observation of the system. We extend our method to simultaneously reconstruct both the entire network topology and all parameters appearing linear in the system's equations of motion. Reconstruction of network topology and system parameters is viable even in the presence of substantial external noise.

Abstract. Inferring network topology from dynamical observations is a fundamental problem pervading research on complex systems. Here, we present a simple, direct method to infer the structural connection topology of a network, given an observation of one collective dynamical trajectory. The general theoretical framework is applicable to arbitrary network dynamical systems described by ordinary differential equations. No interference (external driving) is required and the type of dynamics is not restricted in any way. In particular, the observed dynamics may be arbitrarily complex; stationary, invariant or transient; synchronous or asynchronous and chaotic or periodic. Presupposing a knowledge of the functional form of the dynamical units and of the coupling functions between them, we present an analytical solution to the inverse problem of finding the network topology. Robust reconstruction is achieved in any sufficiently long generic observation of the system. We extend our method to simultaneously reconstruct both the entire network topology and all parameters appearing linear in the system's equations of motion. Reconstruction of network topology and system parameters is viable even in the presence of substantial external noise.

Background
Understanding the relations between network topology and its collective dynamics is at the heart of current interdisciplinary research on networked systems [1]. Often it is possible to observe the dynamics of the individual units of the network, whereas the coupling strengths between them and the underlying network topology cannot be directly measured. Hence, various methods have been proposed to solve the inverse problem of inferring network structure from observation and control of dynamics.
Perturbing a fixed point of a network dynamical system constitutes the simplest controlled intervention of a system. The method of Tegnér et al. [2] perturbs the steady state expression levels of selected genes. Their iterative algorithm can reveal the structure of an underlying gene regulatory network by analysing resultant dynamical changes in the pattern of gene expression levels. A similar iterative method based on multiple regression, coupled with transcriptional perturbations to the fixed points of a genetic network, has been used to successfully identify a nine-gene sub-network [3]. A method introduced by Timme [4] extends reconstruction to networks of smoothly coupled limit-cycle oscillators with periodic collective dynamics. The underlying idea is that the asymptotic response dynamics of a network to different externally induced driving conditions is a function of its topology and of the (external) driving signals. Thus measuring the response to suitable driving signals in different experiments restricts the set of network topologies that are consistent with the driving-response pairs, yielding the network's topology for sufficiently many experiments.
Can we reconstruct a network displaying collective dynamics richer than simple fixed points or limit-cycles? Yu et al. [5] introduced a synchronization method to identify networks of chaotic Lorenz oscillators up to N =17 units. Assuming full knowledge of all model parameters, the network topology of a clone of the system is varied progressively via error minimisation until it synchronizes with the original system. The topology of the clone is then recognised as that of the original network. An extension of this synchronisation method [6] involves additional 'control signals' to externally drive the system to steady states, allowing the inference of interaction topology for sparse, symmetric networks.
Here, we introduce a simple, direct and intervention-free method to reconstruct networks of arbitrary topology from the mere observation of generic collective dynamics. Given the functional form of the intrinsic and interaction dynamics, we show that all other factors of the network dynamical system, such as network topology and coupling strengths (and even typical parameters), can be reliably and efficiently reconstructed.

Theory of Direct Reconstruction from Dynamical Trajectories
Given an observation of the collective trajectory of a dynamical network, how can we infer its underlying topology? We consider a dynamical system consisting of N units, where the dynamics of each unit is specified by an arbitrary set of dynamical equations, x where i, j ∈ {1, 2, . . . N }, and i , · · · , x (D) i ∈ R D describes the state of the i -th unit, and the functions f i , g i : R D → R D mediate intrinsic and interaction dynamics of the D-dimensional unit, and are known. Methods based on copy-synchronisation [5,6] rely on the construction of a new network, with dynamics governed by Eq. 1 and network parameters J ij that are tuned to that of the real network by an error minimisation procedure. Here, we reduce the same reconstruction problem by evaluating the states and their derivatives directly, recognising that the only remaining unknowns in Eq. 1 are the coupling strengths, which are to be determined. The dynamics of the d -th dimension of the i -th unit is given bẏ where t m ∈ R are the times we evaluate Equation 1 and we now writeż for the rate of change d dt z of a scalar variable z. If there are M such times, m ∈ {1, . . . , M }, we have M equations of the forṁ separately. As the equations (3) are uncoupled for any two different dimensions d and d , we treat these separately and drop the index (d) from now on.
Repeated evaluations of the equations of motion (2) of the system at different times t m thus comprise a simple and implicit restriction on the network topology J ij as follows: writing X i,m =ẋ i,m − f i,m , these equations constitute the matrix equation Here the elements of the i -th row of J are given by J i and comprise the sequence (J ij ) j∈{1,...,N } of all input coupling strengths to unit i. Can we rewrite this equation explicitly for J i ? Generically, M > N , and we wish to solve this overdetermined problem by minimising the error function given by for the best (in Euclidean ( 2 ) norm) solutionĴ i . HereĴ ik represents the reconstructed value of the real coupling strength J ik . Equating to zero the derivatives of the error function with respect to the matrix elements, ∂ ∂J ik E i Ĵ i ! = 0, yields an analytical solution to 2 error-minimisation given bŷ and thus the set of input coupling strengths (and input connectivity) of unit i. Evaluating such equations for all i ∈ {1, . . . , N } yields the complete reconstructed networkĴ. This mathematical form of minimum 2 -norm solution is implemented in many mathematical packages (e.g. as the mrdivide function in MATLAB [7]).

Performance for Different Collective Dynamics
How does this theoretical method perform in applications on data? To illustrate the performance of the method and its insensitivity to the type of dynamics, we apply it to four distinct collective dynamical states ranging from simple periodic synchronous dynamics to very complex, highly chaotic asynchronous states. We choose networks of Rössler oscillators that can exhibit a rich repertoire of collective dynamics from multidimensional chaos to periodicity and from global synchrony to asynchrony (see Fig. 1), depending on local unit parameters and coupling functions.

Successful reconstruction from different dynamics
The dynamics of each Rössler oscillator [8] is given by the three ordinary differential equationsẋ where a i , b i , c i are local parameters. The coupling function was set to f (x i , x j ) = x j − x i to induce synchronisation and to f (x i , x j ) = sin (x j ) to prevent it. The local unit parameters a i , b i and c i of the Rössler oscillators are chosen to induce either chaotic or periodic dynamics. The parameters were treated as unknown and are not needed for the reconstruction of the network topology. This is the case in this example, since the equations where the parameters of network topology occur (in the x -dimension) do not contain any other (unknown) parameter. We now demonstrate a reconstruction of the network in all four dynamical paradigms illustrated in Fig. 1: periodic synchronous, chaotic synchronous, periodic asynchronous and chaotic asynchronous collective states.
Reconstruction in praxis works as follows: for networks exhibiting synchronised dynamics, the coupling function f (x i , x j ) = x j − x i is uniformly zero for all time for all units, revealing no information about the network topology. Nevertheless, the network can still be reconstructed from its transient dynamics towards the synchronous state. In general, by substituting X i,m =ẋ i (t m )+y i (t m )+z i (t m ) and G ij,m = f (x i (t m ), x j (t m )) in (4), we find the least-squares reconstructed network according to (6). Since each unit has at most N − 1 incoming links of unknown weights, the diagonals in the reconstructed J are zero. A reliable reconstruction of the network from trajectories (as shown in Figure 1) is illustrated in Figure 2, for all four dynamical paradigms considered.

Quality of Reconstruction
How accurate is such a reconstruction? The quality of reconstruction is defined as the fraction of coupling strengths that are considered correct. Here α ≤ 1 is the required accuracy of the coupling strengths and H is the Heaviside step function, H (x) = 1 for x ≥ 0 and is H(x) = 0 otherwise. The normalized element-wise difference between the reconstructed and real network is where J max = max i ,j |J i ,j | , Ĵ i ,j . Typically, the quality of reconstruction increases with both the sampling rate and the length of trajectory observed, becoming close to 1 even at lower sampling rates for longer times (see Fig. 3).

Required observation time
What is the minimum length of observation required to reconstruct a network? We define to be the minimum length of time of observation at a sampling rate ω required for accurate reconstruction at quality level q at accuracy α. A general observation is that with increasing sampling rate ω or increasing observation time T , the quality increases, due to more accurate information that is obtained about the system's states. In general, however, it is not only the total number T × ω of restricting equations (per node and per dimension) that controls the quality. For instance, at a given observation time, high quality close to Q 0.95 = 1 is reached even at lower sampling frequency if the collective dynamics is more irregular, cf., Figure 3a,b vs. Figure 3c,d. We ascribe this to the lower degree of correlation among observations at different times that is required for accurate reconstruction numerics in (6), e.g. for irregular dynamics, sample points more distant in time provide more relevant information increase about the system as they are less correlated than closer-by points. How does the minimum required observation time scale with system size? Figure 4 shows T 0.98,0.95,200 , the minimum length of time of observation required to have at least q = 98% of the links accurate in strength to an accuracy of at least α = 0.95, sampled at a rate of ω = 200, as a function of N. The numerics suggest that, at fixed ω, T q,α,ω scales sublinearly with network size N for reasonably small 0 < 1 − α 1 and 0 < 1 − q 1 (reasonably large α and q). This implies that the cost of observation does not grow prohibitively quickly, and that even large networks can be reconstructed by a single observation of the collective dynamics.

Robustness: Substantial noise and unknown parameters
Is reconstruction still feasible in the presence of substantial noise? Is it possible to find unknown parameters of the local unit systems? In the preceding examples, none of the unknown parameters appeared in the dimension of coupling, making the problem of reconstructing network topology distinct from that of inferring dynamical parameters. In the following example we illustrate reconstructing networks with intrinsic unit dynamics that are governed by arbitrary functions with K unknown parameters, and where the dynamics are influenced by substantial additive noise ξ have a finite variance, we note that an observation of the dynamics of the system yields a system of equations linear in the K +N unknowns, which we can solve as before. A simple example is a network of Lorenz oscillators [9], where the dynamics The fit suggests that T q,α,ω ∝ N γ , where the exponent of scaling γ ≈ 0.5. (b) An algebraic scaling (black) fits the data (K =100, ) best. Data plotted on a bilogarithmic scale. Linear best-fits (green dashes) overestimate and logarithmic fits (red dots) underestimate reconstruction time. All fits to data from N =100 to N =2000.
of each oscillator is given bẏ where the parameters σ i , ρ i and β i are unknown, and chosen randomly from intervals where the Lorenz system is known to be chaotic: σ i ∈ [9,11] , ρ i ∈ [20, 35] , β i ∈ [2,3]. The noise terms ξ (d) i , d ∈ {x, y, z}, i ∈ {1, . . . , N } are independent identically distributed random variables (on the discrete simulation time scale of ∆t = 0.001) drawn from the standard normal distribution. The network topology and all parameters of the system can be reconstructed and the reconstruction method works despite substantial interference from noise. Figure 5 shows a successful reconstruction of the network and of all dynamical parameters for a network of heterogeneous Lorenz oscillators, where the noise amplitude η = 0.5 is chosen such that it drastically alters the dynamics from its deterministic counterpart (black vs. blue curve in Figure 5a). This illustrates by example that the theory is insensitive to additive noise and capable of successful reconstruction, as desired for generic real-world systems.

Conclusion
In summary, we have introduced a simple robust method to infer connection topology from observations of deterministic and noisy network dynamical systems, where the functional form of the evolution equations is known.
Our method is unique in the following ways. A simple sufficiently long observation of the system dynamics suffices to reconstruct the network topology and coupling strengths. In many cases, experimental access to the system, say, to introduce 'control signals', as in [4,6] may not be possible, and the method proposed here does not require any form of system intervention. Furthermore, the type of collective dynamics is not restricted to, e.g., fixed points, periodic orbits, synchronous states, or any other specific type of motion, cf. [3,4]. Using entire dynamical trajectories, including transients, for the reconstruction, we demonstrate robust reconstruction for a wide range of observed dynamical states, from asynchronous chaotic states to transient states towards global synchrony. Moreover, we treated the system as a 'gray box', where we know some general principles of the system (like the coupling functions) but lack the details, like network structure or intrinsic parameters. Given (e.g. experimental) dynamical observation data, our theory provides an explicit analytic solution to the inverse problem of finding the network structure. This solution is a direct restatement of the differential equations governing the dynamics of the system, and is thus conceptually the simplest possible. This simplicity suggests highest attainable quality for such inverse problems.
The method scales sublinearly with network size, seems robust against substantial addition of noise, and thus provides a promising complement to existing reconstruction methods. Thus, our method offers a conceptual simplification over other methods that make the same assumptions we do, but rely on more complex techniques like copysynchronisation (called auto-synchronization in [6]) or the use of topology estimating clone models or control signals [5]. Further, the method suggested here is capable of reconstructing network structure from a simple observation of the system's dynamics, without resorting to any external intervention to drive the system into or from some canonical state, as in [4].
Efforts to understand the general interplay of network structure and dynamics have yielded several promising approaches, mainly applicable to smaller systems. Notable among the forward methods, i.e., methods that predict dynamical features from knowledge of network topology are those that study the propagation of a harmonic perturbation through a network of coupled phase oscillators [10,11], and methods to predict disordered dynamics from structures of strongly connected components in the network [12]. Cimponeriu et al. [13] introduce two methods to estimate the interaction delay in weak coupling between two self-sustained oscillators from observed dynamical time series. Arenas et al. [14] show that times of synchronisation can reveal the hierarchical structure of a network, revealing a connection between synchronisation dynamics and topological clustering. In an alternative approach, Memmesheimer and Timme [15,16] present an analytical method to design networks of spiking neurons that display a required spike pattern. Other inverse methods have relied on stochastic optimisation [17] to fit a model of a network of spiking neurons to an observation of a real network to infer its topological parameters.
Several avenues for further research present themselves. As suggested by previous work [4,2], minimising the 1 norm, instead of the 2 norm as used here may result in more efficient reconstruction of sparse networks [18]. In addition, our preliminary studies suggest that this highly overdetermined inverse problem can be reduced to an exactly-determined problem by selectively choosing points on the time series to ensure that the resultant system of equations is maximally linearly independent. This reduction significantly reduces the cost of computation to reconstruct large networks. Furthermore, it is straightforward to extend this method to coupled map networks and to systems with delay. Finally, recent studies show that a method analogous to the one presented here for smoothly coupled systems is capable of reconstructing networks of pulse-coupled systems such as integrate-and-fire neurons [19].