{'fraction_non_alphanumeric': 0.02346576388197418, 'fraction_numerical': 0.009326562448139666, 'mean_word_length': 4.999601752289925, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Backpropagation on Dynamical Networks

Dynamical networks are versatile models that can describe a variety of behaviours such as synchronisation and feedback. However, applying these models in real world contexts is difficult as prior information pertaining to the connectivity structure or local dynamics is often unknown and must be inferred from time series observations of network states. Additionally, the influence of coupling interactions between nodes further complicates the isolation of local node dynamics. Given the architectural similarities between dynamical networks and recurrent neural networks (RNN), we propose a network inference method based on the backpropagation through time (BPTT) algorithm commonly used to train recurrent neural networks. This method aims to simultaneously infer both the connectivity structure and local node dynamics purely from observation of node states. An approximation of local node dynamics is first constructed using a neural network. This is alternated with an adapted BPTT algorithm to regress corresponding network weights by minimising prediction errors of the dynamical network based on the previously constructed local models until convergence is achieved. This method was found to be succesful in identifying the connectivity structure for coupled networks of Lorenz, Chua and FitzHugh-Nagumo oscillators. Freerun prediction performance with the resulting local models and weights was found to be comparable to the true system with noisy initial conditions. The method is also extended to non-conventional network couplings such as asymmetric negative coupling.

1. Introduction. Dynamical networks are a common occurrence when studying complex systems where the goal is to model large systems of multiple interacting components. In its simplest form, a dynamical network may be defined as having three main components: local node dynamics, coupling dynamics and connectivity structure (see Equation 1.1). One key feature of dynamical networks is its ability to recreate and describe rich and interesting dynamics such as chimera states, synchronisation [1,2,37] and cascading failure commonly encountered in real world systems [17]. They have also been applied to other systems such as neuron networks [10,17], power grids [34,14], epidemic spread [9], and cardiac arrhythmia [23].
In Eq 1.1, nodes are indexed with i, c represents the coupling structure, f and g represents the local and coupling dynamics respectively. Despite its relatively simple form, data-driven applications of this framework to model real world systems can be difficult if information regarding the coupling structure and dynamics is not accessible. In these cases, these components need to be inferred before a dynamical networks framework can be applied. One simple example would be the data driven analysis of a network of chaotic oscillators. The inverse problem of inferring network characteristics from data has been well studied under specific assumptions and restrictions. However, this problem is difficult due to the complicated interplay between local models and coupling behaviour. Solution attempts often assume that either the local or coupling dynamics are fully known and are used to infer the remaining components or restrict inferences to statistical arguments with uncertainty.
Numerous approximation and inference methods for inferring network connectivity have been developed each with varying success. These methods may be classified into three main approaches: 1. Direct Approach -For systems where both local f and coupling dynamics g is explicitly known, a direct algebraic approach using derivatives have been used to recover connectivity structure [25]. 2. Perturbation Methods -These methods try and recreate a proxy dynamical network that performs similarly to the target system. Small perturbations are applied to individual nodes in the network and observed to see how each signal propagates. These propagations are then used to infer the network connectivity structure [27,4,3] 3. Signal Correlation -This approach focuses on quantifying the correlation between node signals [18]. Node pairs with highly correlated signals are used to infer the connectivity structure. Maximum route entropy has also been used as an alternative to achieve the same result [32]. An extension of this approach is the employment of some notion of causality such as Granger causality [19,13] and short term causal dependence [27]. The inference of connectivity structure also poses an additional challenge of false positive couplings and remains an open research area. Granger causality and similar statistical based approaches have been found to particularly struggle with this problem, even for simple network structures such as rings and chains [13,35] . Additionally, dynamical networks may also have the task of identifying coupling strength in the case of weighted networks.
The second problem of constructing models to describe the local dynamics f in multi-node networks has been studied for pseudo-periodic systems [23,19]. Outside these, many approaches focus on the case of a single node without external interactions and utilise a time-delay embedding with appropriately chosen delay and embedding dimension [28]. Proposed modelling methods for long term prediction include neural network predictors [36], radial basis predictors [26] and reservoir computing [11] among several others. Assessing the quality of these local models rely on the comparison of dynamical invariants such as Lyapunov exponents or correlation dimension [21,16,7] and attractor homology [29].
Despite these topics being well studied, the literature on inferring both connectivity structure and node dynamics from time series observations is not rich. However, we note that extensive work has been done in the inference of coupling functions of dynamical networks where node dynamics are restricted to pseudo-periodic behaviour [23,20]. These assumptions allow a filtration of the phase-dynamics and calculates coupling effects based on the remaining phase defects in the signal. Another proposed method utilises a mean field approach to construct an 'effective network' that acts as an approximate proxy for the real dynamical network, reproducing many of its properties and behaviours. Under specific assumptions such as weak coupling, network heterogeneity and limited links, this method has been shown to perform well in identifying critical transitions and the connectivity structure of a cat cortex [6].
An alternative approach based on the increasingly popular method of compressive sensing [?] has also been proposed. The general approach for compressive sensing assumes that system equations may be approximated by a power series with a small number of terms. Compressive sensing tries to find the sparsest combination of coefficients with respect to a power series basis via convex optimisation with respect to the L 1 -norm. Link and coupling weight identification can also be included into this approach as well [30,31]. However, we note that this method may be limited to simple systems with analytical or well behaved equations.
There is also a third related problem of identifying coupling equations g from signals. Solutions have been proposed using phase dynamics for networks exhibiting pseudo-periodic behaviour, but not necessarily for aperiodic behaviour [23]. For simplicity, the work in this paper assumes that the coupling equations are known.
In this paper, we propose a data-driven approach from simultaneously inferring both the local node dynamics and connectivity structure purely from the observation of node signals. Recognising the structural similarities between dynamical systems and RNNs, we propose an adaptation of the backpropagation through time (BPTT) algorithm commonly used to train RNNs [33], as a potential method of regressing the connectivity structure of a dynamical network. This regression is then used to decouple the observed node signals and simultaneously construct a local dynamics model.
We test the performance of the proposed method on a dynamical network of Lorenz [15] and Chua oscillators [12] with a randomly initialised weighted connectivity structure. We find that an adapted backpropagation regression method was able to successfully reproduce the network coupling structure and local dynamics for a 64 node network. Additional cases with asymmetric and negatively weighted coupling are also explored, yielding similar results. To test the performance on a more realistic application, the proposed backpropagation regression method was applied to a network or driven FitzHugh-Nagumo oscillators operating in the chaotic regime to simulate a simple biological neuron network. Similar to the previous cases, the backpropagation method was able to successfully reproduce both the network structure and local dynamics.

Background.
2.1. Dynamical Networks. Dynamical networks can be defined as a network G = (f, g, C) consisting of three main components, local dynamics f , coupling dynamics g and connectivity structure C. Node states evolve according to Equation 1.1 with local behaviour governed by f and additional external influences from neighbours due to coupling components g and C.
The effect of coupling between nodes results in a high dimensional system with multiple interacting sub-systems. For small or highly connected networks, sufficiently large coupling can cause network-wide synchronisation of node states. As coupling is decreased, chimera states and group synchronisation emerges [1,2,37]. When no coupling is present, all nodes behave independent of each other.
Here, we focus on dynamical networks based on two different 3-dimensional chaotic systems, Lorenz and Chua, with difference coupling in the first component. For the latter, we will focus on the Chua system with cubic nonlinearity φ operating in the single scroll mode [12] to analyse dynamics contrasting with the Lorenz system.

Backpropagation and Recurrent Neural
Networks. Recurrent neural networks (RNNs) is a modelling architecture developed in the field of a machine learning that is widely used for sequence and time series prediction [22,5,24]. RNNs consist of an input layer with weights C (in) , a fully connected hidden layer of nodes h(t) with feedback weights C (h) , and an output layer with readout weights C (out) . For prediction, an input sequence x(t) is fed into the RNN via the input layer with input weights . These inputs are then sequentially fed into the recurrent hidden layer whose nodes feedback to each other according to weights before applying a nonlinear transformation σ. The use of the latter allows the network to approximate nonlinear dependencies between states with nodes evolving according to Equation 2.1.
Predictionsx (t) are calculated using the readout weights and hidden node states as given in Equation 2.2, The training of an RNN involves the small adjustment of network weights C (in) , C (h) , C (out) and biases b. This is usually achieved via the backpropagation through time (BPTT) algorithm [19]. The BPTT algorithm can be interpreted as an adaptation of the widely used gradient descent training algorithm for neural networks that accounts for the propagation of network weight effects through the history of the RNN's hidden node states.
The BPTT algorithm works by unfolding the network predictions at each time step. For a given RNN model, an input sequence x in = {x 0 , x 1 , ..., x τ } is fed in to produce a sequence of outputs x out = {x τ +1 , x τ +2 , ..., x τ +k }. Prediction losses are calculated with respect to some metric d and the real output The gradient of the loss taken with respect to node weights is calculated by propagating the error backward through time steps and used to update node weights (see. Figure 1).
From observing its architecture, the hidden layer of an RNN functions identically to a dynamical network. Hence, we aim to investigate if the BPTT approach can be adapted to regress the weights of a dynamical network given a set of input node states.
3. Methods. This paper proposes a novel method of simultaneously inferring the local dynamics f and connectivity structure C of a dynamical network with a known coupling function g by only observing node states. The proposed method consists of three main stages, namely initialisation, backpropagation and decoupling. An algorithm and flowchart outlining this process is provided (see Figure 2).

Backpropagation Regression Algorithm.
Initialisation. The initialisation stage consists of constructing a local modelf that approximates the true local dynamics f . Due to the complex interplay between the local f and coupling g dynamics, getting an accurate coupling-free estimation of f directly from node signals is not possible.
The estimated local modelf is constructed as a mapping from state space x(t) to vector fieldẋ(t). The time series in numerically differentiated in order to retrievė For this stage of the method, we make the assumption that the magnitude of coupling effects from g is relatively small compared to local dynamics f . This allows a mean field approach to be taken when constructing an approximate local modelf . Here we calculate the estimate true vector field at each point in state space as an average over its K nearest neighbours (see Figure 3), δt . The estimated local dynamicsf is defined as a mappingf : x(t) →x(t). A feedforward network with 2 hidden layers was used to learn the mapping due to its simplicity and ease of implementation. Each hidden layer is composed of 128 neurons with a rectified linear activation function. In addition to the local dynamicsf , values for the estimated coupling adjacency matrixĈ must also be initialised. However, there is no restriction on what values ofĈ may be taken.

Algorithm 3.1 Backpropagation Regression Algorithm
Require: Dynamical network time series x i (t), coupling function g(x i , x j ) Ensure: Local model f (x), Coupling matrix C Step 1: Initialisation Construct imperfect local modelf : x →ẋ for iteration do Step 2: Backpropagation for epoch do Step 3: Decoupling Decouple time series withĈ Retrain local modelf : x →ẋ on decoupled time series

Return Estimates
Backpropagation. The main aim of the backpropagation stage is to provide an improved estimate for the coupling adjacency matrixĈ. Similar to the BPTT algorithm, a forward pass consisting of freerun predictions is first calculated using the estimatedf andĈ with randomly selected initial values from the observed time series in order to properly sample across the whole attractor.
Freerun predictions for all nodesx i (t) are calculated by recursively evaluating Equations 3.3. The loss in freerun predictions due to errors in the local modelf and network weightsĈ can be calculated using some distance function. Here, n is the length of the freerun prediction, N is the number of nodes. This loss at each node, E i , is calculated using the L 2 -norm across the total freerun trajectory {x i (t 0 ), ...,x i (t n )}, For regressing network weights, the loss gradient with respect toĈ must be calculated and used to update the estimate ofĈ, This requires the estimation of partial derivatives with respect to the coupling weights c ij . The expression for these derivatives are given recursively: To provide more stability in the regression, the gradient for each step of backpropagation is averaged over the multiple randomly chosen initial inputs x(t 0 ). An overview of the modified backpropagation is given in Figure 4.
The convergence criteria of the adapted backpropagation algorithm can be calculated for the simple 2 node dynamical network and is given by where E[g 12 ] is the average expected coupling effect between the two node trajectories and α LR is the learning rate hyperparameter. However, convergence for the arbitrary N node dynamical network may not be guaranteed as it simplifies to a variant of the N -body problem. Detailed derivations on the convergence criteria are provided in the Appendix. Notable Hyperparameters. The backpropagation algorithm requires the selection of various hyperparameters that govern the regression behaviour.
• Momentum (M ) -This quantity includes a notion of momentum in the gradient update by allowing previous calculated update steps to propagate into future steps with decaying effect. This technique is also commonly used in RNN backpropagation to improve convergence, • Model Learning Rate (η) -This term governs the weight step size of the feedforward network during the initial construction of the model and retraining. • Node Learning Rate (ᾱ) -The average learning rate that will be applied to each coupling weight if all real links weights are equal in the calculated gradient and exist with probability p. Its relationship to the real learning rate is given by, • Input History Length (t in ) -The amount of steps over which to unfold the backpropagation regression algorithm. Longer history results in the accumulation of coupling weight effects over a longer period and provides better convergence at the cost of slower computation.
Decoupling and Retraining. The use of a mean field approximation to construct an initial local modelf inherently results in significant errors in the coupled dimensions. These errors limit the ability for the algorithm to estimate the true coupling weights as accurate regressions assume that f =f . By using estimated coupling valuesĈ to partially decouple calculated values of the vector field, the output space of f changes to better reflect the local dynamics f . To adapt to this change, the model must be retrained on data of the new partially decoupled output space.
After attaining an improved estimate of the weighted adjacency matrixĈ, the calculated coupling weights can be used to partially decouple the observed node time series according to Equation 3.11.
The new calculated value of the vector fieldf can then be used to improve the previously constructed local dynamics model. The final decoupling step is then alternated with the backpropagation step to gradually improve both the local model and estimate weighted adjacency matrixĈ ij . A much lower model learning rate η must be used to ensure that the structure of the initial model is preserved. The usage of a small learning rate also has the added benefit of providing additional stability to the algorithm as it minimises drastic changes in the loss landscape that can affect convergence.

Performance Assessment.
We need to introduce several measures to assess performance for the proposed backpropagation method. Firstly, prediction horizons were measured for the regressed dynamical network to assess the overall quality of the combined model inclusive of local dynamics and coupling weights. The prediction horizon for node i was defined as the time t p,i where prediction errors exceeded a set threshold averaged across all nodes.
Log weight errors are used to assess the overall error in the regressed weighted adjacency matrix and is normalised the by the true weight values, Average total mutual information I(τ ) was calculated to assess the improvements in the local model without the influence of coupling weights. This is done by using the trained local model to produce freerun predictions of separate system with only a single node (i.e. no external interactions). The average mutual information between freerun prediction trajectories and real model for a given prediction length τ is given as, (3.14) The cumulative area under the mutual information curve (AUC) is taken to track the collapse of mutual information. Larger values indicate more persistent mutual information and a better predictive local model,

Results.
Three main tests were done to assess the performance of the proposed adapted backpropagation regression algorithm. These were the cases of the fully observable network, and non-positive coupling weights. All tests were conducted on an either a 16 or 64 node random network with connection weight probabilities of p = log(N )/N to minimise likelihood disjoint networks. Node coupling strengths were normally distributed N (µ, σ 2 ) with mean and variance dependent on the type of dynamical system to account for different ease of synchronisation between the Lorenz and Chua systems. The values for (µ, σ) were set to be (0.25, 0.05) and (0.15, 0.02) for the Lorenz and Chua systems respectively. Mean coupling strength µ assumed to be non zero to ensure that regressed values clearly distinguish between the absence and presence of edges.
Input data of 25000 steps with an additional 2000 steps washout were generated with timestep dt = 0.02 using the standard RK4 algorithm. Washout steps were discarded to filter out any transient dynamics. An initial learning rate of η = 0.001 was used for training the feedforward model during initialisation. A lower learning rate of η = 0.0002 was used for all subsequent retraining iterations with each training period consisting of 40 epochs.

Fully Observable Network.
The fully observable case assumes that all node time series are accessible for regression. Gradient losses L are calculated at each backpropagation step and used to regress the network adjacency matrixĈ. The regression algorithm was run for a total of 30 refit iterations each consisting of 300 backpropagation steps followed by decoupling and retraining of the model. A total of 8 randomly generated 16 node networks were tested. The quality of the regressed weights and local model for each network was assessed by running freerun predictions with 16 randomly generated initial conditions.
The backpropagation method was able to accurately regress the weights of both the Lorenz ( Figure 5) and Chua ( Figure 6) dynamical networks. The resulting normalised error in the regressed weights for these systems at the end of 30 iterations was 17% and 7% for the Lorenz and Chua networks respectively. The constructed trained local model for both systems also showed a large increase in the AUC for the mutual information with 117% increase for Lorenz and a 54% increase for Chua when compared to the initial local model based purely on the mean field approach (i.e. Iteration 0).
Prediction horizon for both systems was found to increase with each additional iteration. To assess network prediction performance, the predictions using the constructed local model and regressed weights were compared against the control case where the real local model and coupling weights are known, but initial conditions are perturbed with Gaussian noise with a mean of 0 and standard deviation of ξ = 0.001. This method was chosen over the calculation of dynamical invariants such as Lyapunov exponents due to its ease of computation, especially in high dimensional systems such as dynamical network. For the Chua system, the resulting combined model and regressed weights was able to outperform the control resulting in a significant increase in the prediction horizon. This was not the case in Lorenz. However, this is attributed to the natural size scales of Lorenz and Chua. Due to the inherently larger scales of the Lorenz system, a perturbation of ξ = 0.001 is much smaller compared to the Chua system whose oscillations have a much smaller spatial variation.

4.2.
Large 64-node Network. A single larger system of 64 random coupled Lorenz oscillators was used to test the scalability of the system. Weights were gener-ated with µ = 0.15 and σ = 0.02 respectively to ensure there was no synchronisation. The algorithm was run for 80 iterations to account for slightly slower convergence. Performance metrics identical to the regular 16 node case was used with the network tested against 16 different initial conditions. The results of this case is given in Figure  7.
We find that the backpropagation method was also able to accurately regress the desired weights of the the coupled system with a final normalised weight error of 9%. The mutual information of the local model was also found to increase by 58% by the end of the iterations.

Negative Asymmetric
Coupling. Many dynamical systems assume the case of positive coupling resulting in the attraction of nearby trajectories. In addition to this common case, we consider dynamical systems with both negative and positive coupling where repulsive behaviour between nodes is possible. Additionally, we also consider the effects where the adjacency matrix is not symmetrical. To do this, a 16 node adjacency matrix was generated using the same parameters for p. However each node has a probability 0.25 of being changed to a negative value. The time series was then generated similarly to the original test case with the Lorenz equations. The results of this case is given in Figure 8.
The backpropagation method was found to also work equally well in regressing cases with negative weights and asymmetric coupling. The results for the network prediction horizon, log weight errors and mutual information measures show a similar trend to their symmetric, positively weighted counterparts with a final normalised weight error of 6.7% and a 113% relative increase in mutual information for the constructed local model.

Conclusion.
Dynamical networks are an interesting and useful framework to analyse complex systems with multiple interacting components. However, application of this framework is usually difficult due to lack of information pertaining to features such as the connectivity structure and local dynamics. As a result, there is a need for a method to extract these features purely from observed node states.
In this paper, we propose an adapted backpropagation method for inferring the local model and connectivity structure of dynamical networks purely from node time series. This method draws inspiration from the backpropagation through time (BPTT) algorithm commonly used to train RNNs and adapts it for the purpose of regressing couplings weights in a dynamical network. Regressed values for the network coupling is then used to decouple the observed node time series and improve the construction of the local model without coupling effects. The two steps of backpropagation weight regression and decoupling alternate until sufficient convergence is achieved yielding a final local model and regressed coupling weights.
The method was tested on two 16-node networks of Lorenz and Chua oscillators with successful results. The backpropagation regression method was able to accurately identify the coupling weights down to a margin of 7%. Additionally, the constructed local model was compared against the control case and found to show significant improvements. Combined, the regressed weights and constructed local model was able to achieve respectable prediction horizons for the 16 node dynamical network compared to a control case with perturbed initial conditions where the exact model and coupling weights are known. Despite its poor scalability, the method was also able to accurately reconstruct the local model and coupling weights of a larger 64 node network with similar results to the smaller 16 node counterparts. Additional testing was also done on dynamical networks with negative and asymmetric coupling with similar performance.
Acknowledgments. D.C.C. is supported by the Australian Research Council through the Centre for. Transforming Maintenance through Data Science (grant number IC180100030), funded by the Australian Government.
Appendix A. FitzHugh-Nagumo Neuron Network. An additional system consisting of a network of FitzHugh-Nagumo neuron oscillators [8] was used to test the capabilities of the backpropagation regression method. The FitzHugh-Nagumo system is a 2 dimensional driven oscillator modelling the dynamics of neuron voltage blocking (w) and depolarisation (V ) in response to an input driving current (I). This biologically inspired system is able to reproduce stereotypical neuron behaviours such as periodic spike trains and bursting dynamics. This model also exhibits chaotic behaviour when driven by a sinusoidal current with a well-tuned driving frequency ω.
The neuron network presents an additional challenge when constructing a datadriven model due to the presence of disparate time scales in the dynamics (fast spiking depolarisation and slow repolarisation). It also demonstrates the performance of the backpropagation method in a more realistic context, in this case, the inference of neuronal networks. We focus on the FitzHugh-Nagumo system operating under a chaotic regime as given by Hong [8] with equations, Forward Pasŝ For simplicity, we introduce the intermediate variables ∆c and g 12 corresponding to the coupling error and directional coupling respectively. ∆c = c * − c, Hence, we may express the error gradients following expressions, with corresponding coupling weight updates δc given in terms of the gradient and some selected learning rate α LR .
For convergence, we require that errors in the coupling weights decrease at each backpropagation step, If at each step of backpropagation, N initial starting times are used to calculate multiple weight updates δc i . the final weight update taken as the average over all samples can be written as where E[G ij ] and E[X i ] are the expected values of g(x i (t), x j (t)) and x i (t) at any given time t. Assuming that local model errors are independent of position in state space ( f ⊥ X i ) and errors have zero mean, the convergence criterion simplifies to The value of E[G ij ] varies depending on the degree of synchronisation and can be split into three different cases. Case 1: No/Weak Synchronisation X 1 and X 2 are uncorrelated. Therefore X 1 ⊥ X 2 . Convergence is dependent on the spatial spread of the attractor.
Case 2: Strong/Full Synchronisation X 1 ≈ X 2 and system is highly correlated =⇒ . Therefore, the interaction term G 12 = X 2 − X 1 ≈ 0. Therefore the limiting behaviour is divergent and results in the following expression of the convergence criterion, Case 3: Partial/Intermittent Synchronisation Convergence criterion is given by,