Data-driven dynamics reconstruction using RBF network

Constructing the governing dynamical equations of complex systems from observational data is of great interest for both theory and applications. However, it is a difficult inverse problem to explicitly construct the dynamical equations for many real complex systems based on observational data. Here, we propose to implicitly represent the dynamical equations of a complex system using a radial basis function (RBF) network trained on the observed data of the system. We show that the RBF network trained on trajectory data of the classical Lorenz and Chen system can faithfully reproduce the orbits, fixed points, and local bifurcations of the original dynamical equations. We also apply this method to electrocardiogram (ECG) data and show that the fixed points of the RBF network trained using ECG can discriminate healthy people from patients with heart disease, indicating that the method can be applied to real complex systems.


Introduction
The dynamical properties of a system can be studied using various techniques developed for nonlinear dynamical systems, provided we have the explicit governing dynamical equations of the system [1][2][3].However, for many real complex systems, the governing equations are often not available.Therefore, considerable effort has been devoted to inferring the underlying dynamics from noisy or incomplete system observations.Without requiring the dynamical equations, the powerful delay-coordinate embedding method is used to extract statistical quantities that characterise the dynamical invariant set of the original system, while preserving the predictability of the system's evolution [4,5].According to Koopman operator theory, the nonlinear system can be considered an infinite-dimensional linear operator acting on the space of all possible measurement functions of the system, providing tremendous potential for enabling the prediction and estimation of nonlinear systems [6].By combining the delay-coordinate embedding and Koopman operator theory, chaotic nonlinear systems can be decomposed into different dynamic regimes [7].The qualitative information contained in a sequence of observations is used to derive the deterministic part of the system's dynamical equations.Although this process introduces some terms not present in the original equations, the method has shown promise [8].More recently, machine learning techniques have been used to predict the spatio-temporal behaviour of complex systems.For example, reservoir computing has been applied to wildfire models [9], parallel and hybrid reservoirs for chaotic multiscale systems [10], and a nonlinear vector autoregression machine or next-generation reservoir computing (NG-RC) for chaos prediction [11], and a parallel machine learning scheme combined with NG-RC for spatiotemporal chaos [12].Although these machine learning methods have proven effective in predicting transient behaviour, they often fail to capture the asymptotic behaviour of the system.
Based on the concept that many nonlinear dynamical systems in nature and engineering are governed by smooth functions that can be approximated by series expansions, the use of sparse optimisation for the discovery of system equations has been proposed [5], including techniques such as automatic generation of symbolic equations directly from time series data [13,14], compressive sensing for coefficient estimation [15], sparse identification of nonlinear dynamical systems (SINDs) [16], and use of a custom encoder to select measurement coordinates in SINDs [17].Although these methods provide the dynamical equations of the system, their data requirements may be too demanding for many real-world systems.Recently, a general method for extracting symbolic representations from deep learning models by introducing strong inductive biases has been proposed [18], and it has been shown that it is possible to construct new physical theories using machine learning or to apply these techniques to discrete field theory [19,20].These studies focus on reproducing the governing equations of existing models in an explicit way, which is a major challenge in many real systems.
At the same time, the smooth functions can be approximated by a neural network [21][22][23].The finite time trajectory of any n-dimensional dynamical system can be approximated by the internal states of an output unit of a recurrent neural network under certain conditions [24].The radial basis function (RBF) networks, as a special type of neural network with a simple architecture, can not only approximate the continuous functions, but also show several advantages for complex dynamical systems.
The RBF network has good local properties because it achieves global approximation through local approximation by learning the different parts within the input space.The RBF network simplifies the solution of certain problems because it can form more complex decision boundaries compared to multi-layer feedforward neural networks [25].RBF networks can be trained efficiently, resulting in shorter learning times than many other neural network architectures [26].RBF networks show better generalisation to unknown data due to the existence of RBFs.Therefore, we propose an alternative approach to approximate the governing equations of dynamical systems by harnessing the power of RBF networks and exploiting their advantages.We trained an RBF network to approximate the governing equations of complex systems, such as the classical chaotic system.Based on the trained RBF network, we can numerically determine the fixed points of the system and analyse the stability of these points.Our results show that the fixed points and their stability based on the trained RBF network are consistent with those from the original equations of the chaotic systems.We applied this robust method to the cardiac dynamical system and found that the fixed point of the trained network can distinguish patients from healthy individuals, which indicates that our RBF network-based method can be applied to real complex systems and reveal dynamic features related to their fixed points.

Three steps to reconstruct the dynamical system
We use an RBF network to implicitly represent the governing dynamical equations of the system.There are three steps to train an RBF network using trajectory data of a dynamical system (figure 1).The first step is data preparation.For an n-dimensional continuous dynamical system: We can numerically integrate the system using Runge-Kutta method and obtain the trajectory Γ as a time series X i (i = 1, 2, . ..).We then numerically obtain the differential of the time series, Ẋi (i = 1, 2, . ..), using the fourth-order five-point differentiation schemes as follows [27]: where f ′ (t 0 ) is numerical differential of time series at t 0 .For a real complex system, if we have the time series of the observation X i , we can numerically obtain the derivates using the fourth-order five-point differentiation schemes.The second step is to train the RBF network.We use the time series X i from the real system or dynamic equations as the input and the derivate Ẋ as the output to train the RBF network.We apply the MATLAB function newrb and newrbe to accomplish the training process, by varing the centre vector µ, variance of activation function Σ, the connection weight from hidden neurons to output neurons W, and the number of hidden neurons.Then we obtain the well-trained RBF network, which is an approximation of function F(X) and the governing function of the original system can be implicitly written as: The third step is to use numerical method analyse the approximated system which implicitly represents the original governing equations.The trajectory of the approximated system, Γ N , can be obtained applying the traditional Runge-Kutta schemes to the implicit dynamical system (equation ( 3)).The trajectory predicts the evolution of the system.Furthermore, the fixed points of the implicit dynamical system, XN (Network( XN ) = 0) , can be found out using the Particle Swarm Optimization (PSO) algorithm [28].By introduce the small perturbation δe i to the fixed point, the Jacobian matrix of the fixed point can be numerically obtained as: Then the eigenvalue of the Jacobian matrix can be worked out and the type and stability of the fixed point can be determined accordingly.

RBF network overview
We recall the basic elements of the RBF network, which consists of the hidden and the output layers of neurons.A hidden neuron receives all inputs and sends information to all output neurons.The activation of hidden neurons follows RBF as: for l = 1, 2, 3, . . .L, where X is input vector with dimension N, L is the number of hidden neurons, µ l is the centre vector of the lth hidden neuron, Σ is the covariance matrix of the Gaussian activation function.The output neuron weights the output of each hidden neuron as follows: for k = 1, 2, . . ., K, where w l,k is connection weight from lth hidden neuron to kth output neuron.Then we can obtain the error of the network: By varying the covariance Σ, centre vector µ, connection weights W and the number of hidden neurons, the RBF network can approximate arbitrary continuous function at desire accuracy.

Data preparation and processing
In this study, we utilized the Lorenz system, the Chen system, and electrocardiogram (ECG) data as real-world example to demonstrate the process of reconstructing dynamical systems through a data-driven approach.
The dynamic equations of Chen system is described as [30] where parameters a and b were fixed at 35 and 3, respectively, while parameter c was varied from 14 to 40.Data collection for Chen system was similar to that for the Lorenz system, but with trajectory lengths of 500 time steps.
For real complex systems, we analysed an open research dataset for 12-lead ECG [31].The dataset includes 10 s, 500 Hz, 12-dimension ECGs, along with labels for rhythms and other conditions for each patient.To increase the number of data points, we upsampled the ECG data from 500 Hz to 3000 Hz and then applied the Butterworth bandpass filter to smooth the data at lower and higher cut-off frequencies of 1 Hz and 40 Hz, respectively.Since the 12-dimensional ECG signal can be linearly transformed into a 3-dimensional vectorcardiography (VCG) without loss of information about cardiac dynamics, we simplified the data by transforming it into a 3-dimensional VCG using specified equations [32]: In this way, we obtained a three-dimensional VCG sequence.The data were then divided into five groups, with the first 5000 data points from each participant reserved for the construction of the dataset.

Training the RBF network
A supervised training approach is used, where the training set consists of time sequences and corresponding derivative sequences.These sequences are randomly divided into a training set and a validation set, with 70% of the data used for training and 30% for validation.We use the newrb function available in MATLAB to train the network, using a backpropagation algorithm in conjunction with the minimum mean squared error criterion.This algorithm optimises weights and biases to minimise the mean squared error between the network output and the target.Gradient descent methods are used during training to update the weights and biases.
The network has a single hidden layer.The number of neurons in the hidden layer is automatically determined by the newrb function, which selects the number and distribution of RBFs in the data space based on the distribution of the training data and the target error requirements.
The initialisation of connection weights is done automatically by the newrb function, using a random initialisation method.This means that the initialised weights are random numbers following a uniform distribution, which helps prevent the network from falling into local optima and increases network diversity.The initial number of hidden layer units is adaptively determined by the function based on the training data and the target error requirements.
When training the network for the Chen and Lorenz systems, we set two hyperparameters of the newrb function to: spread = 2 and goal = 10 −8 .The spread parameter controls the operating range of the neurons and must be adjusted according to the characteristics of the data set.Smaller spread values result in a more localised response range of the neurons, while larger values result in a wider response range.The goal parameter controls the stopping conditions of the training process and represents the error threshold between the network output and the target output.Smaller goal values bring the network closer to the training data but risk overfitting, while larger goal values can lead to underfitting.These values need to be determined based on the specific task and data set.For the cardiac dynamical system we choose 25 000 hidden units and set Spread to 1 in the newrb function, other settings are identical to those of the previous network for classical chaotic systems.
The MATLAB code for this work is available on GitHub at https://github.com/congcongdu/RBF.

Reconstruction of Lorenz system
We use the RBF network to approximate the Lorenz system, including its chaotic behaviour.We first numerically integrate the system using the fourth-order Runge-Kutta method with a time step ∆t = 0.001 over a period T, given the parameter   We note that the orbit is very similar to the orbit integrated from the original Lorenz system starting from the same initial point.This suggests that the RBF network can effectively capture and reconstruct the chaotic behaviour of the system.In addition to the orbit, we identified three fixed points of the approximated system using the PSO algorithm: (0.004, 0.007, −0.026), (−8.485, −8.486, 26.999), (8.485, 8.485, 27.002).The numerical calculation shows that the eigenvalues of the fixed points are −22.854,11.826, −2.643 for the origin and −13.855, 0.0960 ± 10.189i for the symmetric fixed points, respectively.The analytical results of the fixed points and their eigenvalues from the original system are (0, 0, 0),(±6 √ 2, ±6 √ 2, 27) and their eigenvalues are −22.328,12.328, −8/3 for the origin and −13.855, 0.094 ± 10.195i for the symmetric fixed points, respectively.It's obvious that the reconstructed fixed points and their eigenvalues agree with the original system.To further explore the local bifurcation of the approximated system, we change the parameter ρ from 0 to 35 and train the RBF network.Panels (a) and (b) in figure 3 show the input structure of these networks and their network structure.For the Lorenz system, the theoretical values of the equilibrium points can be derived directly from its dynamic equations.If ρ ⩽ 1, there is a single fixed point at the origin (0, 0, 0), and if ρ > 1, there are three fixed points: the origin (0, 0, 0) and two symmetric fixed points . Both local and global bifurcations have been extensively discussed in previous literature, e.g.[33].
We then compute the equilibrium of the trained networks.Using the PSO algorithm, we identify the fixed points of each trained RBF network and derive the eigenvalues for each fixed point.As shown in figure 3(c), the RBF network reconstructed system exhibits a supercritical pitchfork bifurcation at ρ = 1 and a Hopf bifurcation at ρ = 24.7.Before the pitchfork bifurcation there is a single stable fixed point near the origin.If 1 < ρ < 24.7, there is a saddle point with a positive eigenvalue near the origin and two stable foci.If ρ > 24.7, there is one saddle point and two unstable foci with a negative real eigenvalue.Therefore, the local bifurcation of the system reconstructed by the RBF network corresponds to that of the original Lorenz system.
We observe that the number of hidden neurons in the trained RBF network varies with the bifurcation parameter ρ.When the parameter is close to the supercritical pitchfork bifurcation point (ρ = 1), the trajectory of the original Lorenz system is more 'regular' , requiring more hidden neurons in the network to extract the contained information and approximate the original system.Conversely, when the parameter ρ is significantly greater than 1, the trajectory of the original Lorenz system exhibits more randomness.In this case, the RBF network can more easily capture and represent the properties of the original system, resulting in a trained network with fewer hidden neurons, as shown in figure 3(d).

Reconstruction of Chen system
For the Chen system, we set the parameters a = 35, b = 3 and change the parameter c from 14 to 40, using the system reconstruction method proposed in this paper.Figure 4(a) shows the evolution trajectories of the network and the real system starting from the same point when parameter c = 28.Although the trajectories do not completely overlap, their evolution patterns are remarkably similar.We therefore examine their equilibrium points and bifurcation conditions.It's easy to see that the fixed points and their eigenvalues of the RBF network reconstructed system match those of the original Chen system, as shown in figure 4(b).This indicates that the local bifurcation is consistent with that of the original Chen system.

Reconstruction of cardiac dynamical system
In the previous sections we have shown that the RBF network can reconstruct a classical chaotic system.Our method allows us to predict the evolution of the system, an observable process, as well as to locate the fixed points of the system, which are usually not observable in an evolving real system.To test the validity of our method, we choose the ECG as a representative example and use the RBF network to reconstruct the dynamics of the ECG and identify its fixed point.Previous research has shown that the RBF network can detect ECG patterns [34].We postulate that the heart, as a typical dynamical complex system, should contain specific fixed points.Furthermore, we believe that the fixed point of a healthy individual's heart may be different from that of a patient with heart disease.
The trained RBF network can generate time series of vectorcardiography (VCG) and its derivatives.As shown in figure 5, the derivatives produced by the trained RBF network for each group of participants are consistent with the numerical differentiation of the VCG data from the first participant in each group.The right panel of figure 5 shows the difference in derivatives between the VCG data and the trained RBF network, indicating that the dynamics of the VCG can be effectively captured by the trained RBF network.
Internal medicine physicians can differentiate between healthy individuals and heart disease patients based on the ECG waveform.We propose that diseases not only change the ECG waveform, but also change the fixed points of the heart's dynamical system.Therefore, we should discriminate between healthy individuals and patients using the fixed points of the dynamical system reconstructed by the RBF network.We used the PSO algorithm to identify the fixed point of each reconstructed system.Despite the possibility of multiple fixed points in the reconstructed system, we focused on a single fixed point for each network in this study.Given our methodology of using a 'leave one out' approach and introducing new subject data to train a new RBF network, the zero point of the new RBF network should not change significantly.We start from the fixed point fp of the trained network with single fixed point and calculate the distance between fp and each fixed point of other trained networks in the same group.We kept the fixed point fp * that was closest to fp.As shown in figure 6, the retained fixed points of each trained network can be easily classified into healthy individuals and patients by a hyperplane obtained by iterative optimization using centroid-guided initialization and random perturbations.The confusion matrix in table 1 shows the performance of the classifier,where the accuracy is (95 + 138)/(95 + 3 + 13 + 138) ∼ 0.936, the True Predict Rate is 95/98 ∼ 0.969, the False Predict Rate is 13/151 ∼ 0.086, suggesting a good performance of the classification [35,36].
Figure 6 shows that the fixed points of the trained RBF networks can be classified into two groups: one representing healthy individuals and the other representing heart disease patients.Thus, the trained RBF network not only successfully reconstructs the dynamics of the cardiac system, but also reveals the underlying fixed points of the system, thereby capturing the critical factors of the original dynamical system.

Conclusion and discussion
In this study, we used the RBF network trained on trajectory data from classical Lorenz and Chen systems, as well as from real complex systems, to reconstruct their dynamics.The reconstructed system not only predicts system evolution by reproducing trajectories, but also reveals the fixed points behind the observable data of complex systems.Our approach represents a promising method for data-driven modelling of dynamical systems with known or unknown equations and for reconstructing the dynamics of complex systems based on observational data, thereby improving our overall understanding of the dynamical mechanisms underlying the evolution of complex systems.Several points are worth noting: First, the RBF network remains robust to noise in the training data.In fact, we replicated the Lorenz system experiment with fixed parameters of [10, 28, 8/3] thirty times, with each instance using a different training set, resulting in different trained RBF networks.However, the zeros of these RBF networks showed very little variation (0.04%), indicating that the differences are insignificant.Therefore, the RBF network robustly reconstructed the original dynamical system.
Second, we considered the trained RBF network as a function and used the Particle Swarm Optimisation (PSO) algorithm [28] to identify the zeros of the RBF network, which represent the fixed points of the original system.In dynamical systems, fixed points determine the local and possibly global properties of complex systems [1][2][3].We therefore propose a method for identifying the fixed points of some real complex systems, revealing many of the key properties of the system.
Third, as we noted in the introduction, considerable effort has been devoted to reconstructing the dynamics of complex systems from observational data.The delay-coordinate embedding method can extract the geometric properties of the system's invariant sets [4], but it cannot predict the system's evolution.Koopman operators have been used to decompose complex systems into different dynamic states [7].NG-RC and other machine learning methods can effectively predict the transient rather than the asymptotic behaviour of complex systems [9][10][11].SINDs [17] or series expansion approximations [5,15] have been used for explicit reconstruction of specific dynamical system equations.These approaches require high-quality observational data, and the terms of the reconstructed equations are sparse.In the field of Recurrent Neural Networks (RNN), studies have understood the mechanisms of RNN operation through fixed points and the linearised dynamics around these points [37].There are also Tensorflow toolboxes for identifying and characterising RNN fixed points [38].Our focus differs from these methods in that we prioritise the use of neural networks to approximate and analyse real-world dynamical systems.Based on the universal Approximation Theorem, we used the flat RBF networks to approximate the governing equations of dynamical systems.The method is simple but versatile.
Fourth, we also considered several other network architectures, such as reservoir networks [9] and graph neural networks (GNNs).Because the dynamics of the units in reservoir networks and their output always change with time, they cannot be considered as constant functions, making them unsuitable for the task at hand.GNNs, on the other hand, generally require a known network structure to function properly.For example, in the study by Zhang et al [39], they used the Gumbel softmax trick and GNNs to reconstruct the network structure and predict system dynamics from time series data.However, our input is only a time series data without any information about the network structure, so it is completely data-driven.If we consider the simultaneous use of time series data and network structure, the choice of neural network would favour GNNs over the currently used RBF neural networks.
However, our approach has some limitations.While the universal Approximation Theorem exists [21,22], the architecture of our network is overly simplistic, making it difficult to identify suitable parameters to approximate the governing equation in some complex systems.One possible solution is to use deeper or more complex networks to reconstruct the dynamics of complex systems.Indeed, recent work has used a deep neural network architecture to approximate the nonlinear operator [40].In addition, certain dimensions of observed data from real complex systems may depend on other dimensions of the same data.Therefore, we need to decouple the correlated dimensions before using our method to reconstruct the dynamics of the system.

Figure 1 .
Figure 1.Three steps to reconstruct the dynamical system using RBF network.Numerical derivate is the first step, training the RBF network is the second step, and numerical analysis is the last step.
[σ, ρ, β] = [10, 28, 8/3] [29].The resulting time series is represented by the matrix [r x ; r y ; r z ] 3× T ∆t .Using fourth-order five-point differentiation schemes, we compute the derivative matrix [de x ; de y ; de z ] 3× T ∆t from the time series {r i }.The RBF network is then trained with these data using the MATLAB function newrb, yielding a trained network and an approximation of the Lorenz system: [de x ; de y ; de z ] 3× T ∆t = Network([r x ; r y ; r z ] 3× T ∆t ).Panels (a)-(c) in figure 2 illustrate a short segment of the time series from the training set, the derivative series, and the structure of the RBF network, respectively.Panel (d) shows the trajectory of the autonomous evolution of the network after training and its equilibrium points.

Figure 2 .
Figure 2. Using the RBF network to approximate Lorenz system.(a) Time series of chaotic Lorenz attractor.(b) The numerical derivatives of the time series in (a) using the fourth-order five-point numerical differentiation schemes.(c) Schematic representation of the RBF network structure.(d) Example trajectory from the Lorenz equations and from the approximated RBF network.

Figure 3 .
Figure 3. Reconstruction of the local bifurcation of the Lorenz system.(a) The time series of the Lorenz system with different parameters ρ.(b) Schematic of the RBF networks.(c) The local bifurcation of the approximated system and the original Lorenz system.The line for the original Lorenz system and points for the RBF approximation.Blue circles represent saddle points, while red and blue asterisks indicate stable and unstable foci, respectively.(d) The number of hidden neurons in the trained RBF network for different ρ.

C-C Du et alFigure 4 .
Figure 4. Reconstruction of the Chen system.(a) Example trajectory from the original Chen system and the system reconstructed by the RBF network.(b) The local bifurcation with respect to the parameter c.Lines for the original system, points and stars for the reconstructed system.Blue circles represent saddle points, while red and blue stars indicate stable and unstable foci, respectively.

Figure 5 .
Figure 5.The derivatives of the VCG from the data and from the results of the trained RBF network for subject 1 in each subject type.Left panel: the derivative of each group, red lines are the numerical derivative of the VCG data and green lines are the series reconstructed by the trained RBF network.Right panel: the error between the VCG data and the reconstructed series, blue, green and red represent the dimensions of v1, v2 and v3 respectively.

C-C Du et alFigure 6 .
Figure 6.The retained fixed points of the reconstructed dynamical system approximated by the trained RBF network are shown in v1 − v2 − v3 3D space, v1 − v3 plane and v1 − v2 plane, respectively.Green dots represent the reconstructed network of healthy subjects and red dots represent the patients.A classification plane was iteratively optimised using centroid-guided initialisation and random perturbations.The final equation providing optimal separation is 0.7416v1 + 0.6701v2 −0.0319v3 + 0.0091 = 0.

Table 1 .
Confusion matrix of classification on the retained fixed points of the RBF network.