Neural Symbolic Regression of Complex Network Dynamics

## Abstract

Abstract Complex networks describe important structures in nature and society, composed of nodes and the edges that connect them. The evolution of these networks is typically described by dynamics, which are labor-intensive and require expert knowledge to derive.
However, because the complex network involves noisy observations from multiple trajectories of nodes, existing symbolic regression methods are either not applicable or ineffective on its dynamics.
In this paper, we propose Physically Inspired Neural Dynamics Symbolic Regression (PI-NDSR), a method based on neural networks and genetic programming to automatically learn the symbolic expression of dynamics. Our method consists of two key components: a Physically Inspired Neural Dynamics (PIND) to augment and denoise trajectories through observed trajectory interpolation; and a coordinated genetic search algorithm to derive symbolic expressions. This algorithm leverages references of node dynamics and edge dynamics from neural dynamics to avoid overfitted expressions in symbolic space. We evaluate our method on synthetic datasets generated by various dynamics and real datasets on disease spreading. The results demonstrate that PI-NDSR outperforms the existing method in terms of both recovery probability and error.

## 1 Introduction

Complex networks (Gerstner et al., ; Gao et al., ; Bashan et al., ; Newman et al., ) describe important structures in nature and society,
which is composed of a set of nodes and a set of edges that connect them.
Complex networks can model various real-world systems, including social networks (Kitsak et al., ), epidemic networks (Pastor-Satorras & Vespignani, ), brain networks (Laurence et al., ; Wilson & Cowan, ), and transportation networks (Kaluza et al., ).
Extensive works (Zang & Wang, ; Murphy et al., ; Gao & Yan, ) have been conducted to analyze the dynamics of complex networks (Pastor-Satorras et al., ; MacArthur, ; Kuramoto & Kuramoto, ), which is crucial for understanding the underlying mechanisms of complex systems and predicting their future behaviors.
Among them, the symbolic complex network dynamics is of great importance as it offers a clear and concise depiction of the internal mechanisms and their impact on the overall system behavior.

Obtaining the symbolic expressions for complex network dynamics is challenging and is usually done by mathematicians or physicists.
The process of symbolic regression of complex network dynamics typically requires repetitive trial-and-error attempts and profound expertise in the relevant field.
Therefore, it is often unlikely to find the symbolic expressions without any prior knowledge of the system (Virgolin & Pissis, ).
The goal of this paper is to develop a machine-learning method to automatically learn the symbolic expressions of complex network dynamics from data.

Symbolic regression aims to discover the underlying symbolic formula from data (Schmidt & Lipson, ; Petersen et al., ; Cranmer et al., ; Shi et al., ).
Various techniques have been applied to symbolic regression for dynamics, including genetic programming (Gaucel et al., ; Kronberger et al., ), sparse regression (Brunton et al., ), deep-learning-based genetic programming (Qian et al., ), and Transformer (d’Ascoli et al., ).
However, all these methods are designed for symbolic regression of dynamics for a single trajectory.
Directly applying them to complex network dynamics ignores common patterns between trajectories and leads to inefficient learning and poor performance (Gao & Yan, ), as the complex network is an evolving system with multiple observed trajectories of nodes sharing the same underlying dynamics.
To deal with multiple trajectories, recent work (Gao & Yan, ) proposes a Two-Phase Sparse Identification of Nonlinear Dynamics (TP-SINDy) for symbolic regression of complex network dynamics.
The two-phase optimization in TP-SINDy effectively narrows down the symbolic search space and improves its performance.
However, TP-SINDy relies on the estimated time derivatives of the node states, which are noisy and inaccurate.
Additionally, since it is built on SINDy (Brunton et al., ), there exists a pre-defined function library for regression. This can make it challenging to recover the symbolic expression beyond the confines of the function library.
In this paper,
we propose Physically Inspired Neural Dynamics Symbolic Regression (PI-NDSR) for the complex network dynamics based on the observed trajectory,
which is more accurate and
robust than the previous methods.

Figure: Figure 1: Our method is designed for symbolic regression of complex network dynamics. The proposed method interpolates and denoises complex network observations with neural dynamics to avoid inaccurate estimation of time derivatives and applies a coordinated genetic search algorithm to derive the symbolic expressions of complex network dynamics.
Refer to caption: x1.png

The proposed method (Fig. 1) consists of two parts: (a) A physically inspired neural dynamics (neural dynamics refers to the dynamics represented by neural network) disentangles the node dynamics and edge dynamics to augment and denoise trajectories with the observed trajectory interpolation.
Neural dynamics provides high-quality supervision signals for calculating fitness in the genetic search algorithm and gives references for symbolic expressions to enhance search efficiency.
(b) A coordinated genetic search algorithm regresses the symbolic expression of network dynamics that uses references of node dynamics and edge dynamics from neural dynamics to constrain the symbolic search space, which helps avoid overfitted expressions.
Compared with the previous methods which relies on the estimated time derivative, PI-NDSR applies to more general complex network observations with larger noise and sampling intervals and does not rely on a pre-defined function library. The proposed method is evaluated on synthetic datasets generated with various dynamics and a real dataset from disease spreading. The results demonstrate that our method surpasses the previous approach in terms of recovery probability and error.

#### Notations

Matrix, vector, and scalar are denoted as bold capital letters $\mathbf{X}$, bold lowercase letters $\mathbf{x}$, and lowercase letters $x$, respectively.
The element in $i$-th row and $j$-th column of matrix $\mathbf{X}$ is denoted as $\mathbf{X}_{ij}$.
The $v$-th row of matrix $\mathbf{X}$ is denoted as $\mathbf{x}_{v}$.
A complex network is denoted as $\mathbf{G}=(\mathsf{G},\mathbf{X}(t),t\in\mathcal{T})$. $\mathsf{G}=(V,E)$ is the structure of complex network where $V$ is the set of nodes and $E$ is the set of edges.
$\mathbf{X}(t)=[\mathbf{x}_{1}(t)^{\top},\cdots,\mathbf{x}_{N}(t)^{\top}]^{\top
}\!\!\in\mathbb{R}^{N\times d}$ is a $d$-dimensional node states of $N$ nodes at the timestamp $t$, and $\mathcal{T}=\{t_{0},t_{1},\cdots,t_{K-1}\}$ is the set of $K$ timestamps of complex network observations. $\dot{\mathbf{x}}(t)$ represents the time derivatives of $\mathbf{x}(t)$.

## 2 Related work

### 2.1 Symbolic regression

Throughout the history of physics, extracting elegant symbolic expressions from extensive experimental data has been a fundamental approach to uncovering new formulas and validating hypotheses. Symbolic Regression (SR) is a notable topic in this context (Schmidt & Lipson, ; Petersen et al., ; Cranmer et al., ; Kamienny et al., ), aiming to mimic the process of deriving an explicit symbolic model that accurately maps input $X$ to output $y$ while ensuring the model remains concise. Traditional methods for deriving formulas from data have predominantly relied on genetic programming (GP) (Schmidt & Lipson, ; Koza, ; Worm & Chiu, ), a technique inspired by biological evolution that iteratively evolves populations of candidate solutions to discover the most effective mathematical representations.

More recently, due to the remarkable accomplishments of neural networks across diverse domains, there has been an increasing interest in leveraging deep learning for symbolic regression. Specifically, some recent works (Cranmer et al., ; Chen et al., ; Qian et al., ; Udrescu & Tegmark, ; Martius & Lampert, ; Mundhenk et al., ; Shi et al., ) have explored guiding genetic programming with the output of neural networks to improve the efficiency and accuracy of symbolic regression. This approach takes advantage of the powerful pattern recognition and generalization capabilities of neural networks to inform the evolutionary processes of genetic programming, resulting in more effective and efficient discovery of symbolic expressions.
Another line of works (Kamienny et al., ; Biggio et al., ; Valipour et al., ; d’Ascoli et al., ) applies Transformer to symbolic regression and achieves comparable performance to GP-based methods.

### 2.2 Complex network dynamics learning

With the development of machine learning in recent decades, many works have been proposed to learn the dynamics of complex networks from data.
The work for learning dynamics of complex networks can be divided into two categories: dynamics learning with neural networks and dynamics learning with symbolic regression.

Dynamics learning with neural networks utilizes deep neural networks to learn dynamics. Neural dynamics usually follows the encode-process-decode paradigm (Hamrick et al., ; Zang & Wang, ) with an encoding neural network
to pre-process the initial node states by mapping them to latent states.
Then the latent states are processed by a graph neural network to capture the evolution and interaction of the complex network.
Murphy et al. () propose to use GNN to model the evolution of complex networks with regularly sampled
observation.
NDCN (Zang & Wang, ) models the continuity of the dynamics with graph neural ODE (Chen et al., ) and the interaction of the dynamics with GNN. Even though neural dynamics learning can capture the complex dynamics of complex networks, it is difficult to interpret the learned dynamics due to the black-box nature of neural networks.

On the other hand, dynamics learning with symbolic regression aims to learn symbolic expressions of complex network dynamics from data.
Two-phase sparse identification of nonlinear dynamical systems (TP-SINDy) (Gao & Yan, ) is proposed to obtain the dynamics of complex networks from data.
The proposed method includes a function library including a wider range of elementary functions than the orthogonal basis functions in SINDy (Brunton et al., ) and a two-phase regression method to learn the dynamics of complex networks.
However, TP-SINDy requires an accurate estimation of the time derivative of the node states and can only learn the symbolic expressions from the pre-defined function library.
Our work uses the supervision signal and dynamics references provided by neural dynamics, based on genetic programming to learn the symbolic expressions of complex network dynamics.

## 3 Problem setup

The complex network dynamics is defined by the following differential ordinary equation:

$$ $\dot{\mathbf{x}}_{v}(t)=F\left(\mathbf{x}_{v}(t)\right)+\sum\nolimits_{u\in N_ {v}}\!\!\!a_{vu}G\left(\mathbf{x}_{v}(t),\mathbf{x}_{u}(t)\right),$ (1) $$

In (1), $\dot{\mathbf{x}}_{v}(t)$ denotes the time derivative of $\mathbf{x}_{v}(t)$, $F(\mathbf{x}_{v}(t))$ denotes the node dynamics term of node $v$, which includes processes like influx, degradation, or reproduction. $G(\mathbf{x}_{v}(t),\mathbf{x}_{u}(t))$ is the edge dynamics describing the interactions between node $v$ and node $u$, accounting for processes such as spreading and competition.
$G$ is shared across all edges in the network because of the universality in network dynamics (Barzel & Barabási, ; Gao et al., ).
$a_{vu}$ is the weight of the edge between node $v$ and node $u$, and $N_{v}$ is the set of neighbors of node $v$.

Given observations of the node states $\left\{\mathbf{X}(t)|t\in\mathcal{T}\right\}$,
the symbolic regression of complex network dynamics (Gao & Yan, ) aims to find the symbolic expressions of function $F$ and $G$ in (1).
Compared with traditional symbolic regression, we identify two main challenges of symbolic regression for network dynamics:
(a) Directly regressing symbolic function on time derivative is inaccurate and unstable because the estimation of derivative $\dot{\mathbf{X}}(t)$ is noisy and its accuracy relies on the regular sampling intervals of the observation.
(b) The network dynamics consist of two parts, namely node dynamics $F$ and edge dynamics $G$.
In symbolic regression, it is possible that $F$ and $G$ jointly perform well on the observed (interpolation) trajectory but exhibit significant discrepancies in predicting node dynamics and edge dynamics compared to the real dynamics, which results in bad extrapolation performance (see examples in Section B.4). We refer to this phenomenon as the overfitting problem in symbolic space.
In this paper, we propose a physically inspired symbolic regression algorithm,
addressing the first challenge with interpolated network observations and the second challenge with a coordinated genetic optimization algorithm.

## 4 Method

To address the challenge of inaccurate estimation of $\dot{\mathbf{X}}(t)$ in symbolic regression, we instead utilize the state of the nodes as supervision.
We train a neural network $f_{\theta}$ with parameter $\theta$ to denoise and interpolate the observed trajectory as the approximation for the true observation $\mathbf{X}(t)$. To address the challenges in symbolic regression of node dynamics and edge dynamics in (1), we design a coordinated genetic optimization algorithm based on the algorithmic alignment between neural network $f_{\theta}$ and network dynamics (1).

### 4.1 Trajectory denoising and augmentation with interpolated observations

In symbolic regression, it is beneficial to regress symbolic expressions with many low-noise network observations.
However, the network dynamics observation is usually sparse (i.e., sample time interval is large) and noisy.
To address this issue, we propose Physically Inspired Neural Dynamics (PIND) to interpolate and denoise the network observations, obtaining denser and less noisy trajectories. Except for the interpolated trajectory, the neural dynamics aligns with network dynamics in (1) and provides an estimation of the node dynamics and edge dynamics, which is used as references for coordinated genetic optimization in deriving the symbolic $F$ and $G$.

Following existing works (Zang & Wang, ; Murphy et al., ),
the neural dynamics is designed based on encode-process-decode paradigm (Battaglia et al., ).
To obtain the denser and less noisy trajectory and the high-quality estimation for node dynamics and edge dynamics for genetic search, we incorporate the algorithmic inductive bias in (1) to design a physically-inspired neural dynamics and train the neural dynamics from the observed trajectory.
In neural dynamics, the raw observation $\mathbf{x}_{v}(t)$ is first encoded by an encoding network $\mathsf{Enc}$ to obtain the latent state $\mathbf{h}_{v}(t)$, and the graph neural network is applied to calculate the time derivative of the latent state $\dot{\mathbf{h}}_{v}(t)$, which is later decoded by a decoding network $\mathsf{Dec}$ to obtain the time derivatives $\dot{\mathbf{x}}_{v}^{\circ}(t)$. Finally, differential ODESolver (Chen et al., ) is applied to $\dot{\mathbf{x}}_{v}^{\circ}(t)$ to evolve the neural dynamics to obtain more network observations.

Based on (1), the time derivative $\dot{\mathbf{x}}_{v}^{\circ}(t)$ for node $v$ is designed with inductive bias from (1):

$$ $\displaystyle\dot{\mathbf{x}}_{v}^{\circ}(t)=\mathsf{Dec}(\dot{\mathbf{h}}_{v} (t)),\quad\text{s.t.}\quad\begin{cases}\dot{\mathbf{h}}_{v}(t)=\phi^{n}\big{(} \mathbf{h}_{v}(t),t\big{)}+\sum\nolimits_{u\in{N}_{v}}\phi^{e}(\mathbf{h}_{v}( t),\mathbf{h}_{u}(t),t),\\ \mathbf{h}_{v}(t)=\mathsf{Enc}(\mathbf{x}_{v}(t)),\mathbf{h}_{u}(t)=\mathsf{ Enc}(\mathbf{x}_{u}(t)),u\in N_{v},\end{cases}$ (2) $$

where $\phi^{n}$ and $\phi^{e}$ are two MLPs aligning with the node dynamics and edge dynamics in (1), respectively, and $\mathbf{h}_{v}(t)\in\mathbb{R}^{d^{\prime}}$ is the latent state of node $v$.
In (2), the neural node dynamics $\phi^{n}$ captures the evolution of nodes influenced by their properties, and the neural edge dynamics $\phi^{e}$ captures the interactions between two end nodes of an edge. Therefore, PIND of node $v$ is written as

$$ $f_{\theta}(\mathsf{G},\mathbf{X}(t_{0}),t)_{v}=\text{ODESolver}(\dot{\mathbf{x }}_{v}^{\circ}(t),\mathbf{X}(t_{0}),t_{0},t).$ (3) $$

The alignment between neural dynamics (3) and dynamics formulation (1) enables better learning of complex network dynamics.
To train the neural dynamics, we minimize the error between $f_{\theta}(\mathsf{G},\mathbf{X}(t_{0}),t)_{v},\forall v\in V$ and the observed trajectories $\{\mathbf{X}(t)|t\in\mathcal{T}\}$ with standard deep learning optimization techniques.
After the training, we use the interpolated trajectory $\hat{\mathbf{X}}(t)$ from $f_{\theta}(\mathsf{G},\mathbf{X}(t_{0}),t)_{v},\forall v\in V$ as the supervision signal for symbolic regression and the estimated node dynamics and edge dynamics as references for the genetic search in deriving the symbolic $F$ and $G$ in Section 4.2.

### 4.2 Coordinated genetic search for symbolic regression

To regress the symbolic formulas $F$ and $G$, we propose an algorithm called coordinated genetic search to obtain the symbolic expressions of complex network dynamics.
Based on the physically-inspired design of (2), we can obtain the references for symbolic expressions by disentangling the neural dynamics into node dynamics and edge dynamics.
Then, we use the neural references to coordinate the search process of symbolic expressions of $F$ and $G$ to avoid the overfitted expressions.

Because of the inductive bias and algorithmic alignment of (1), the neural dynamics often have better generalization performance (Xu et al., ; Veličković & Blundell, ).
Therefore, the neural node dynamics and edge dynamics can be used as estimations of the node dynamics and edge dynamics in (1), which can be used as references for coordinating the search process of symbolic expressions of $F$ and $G$.
The neural node dynamics and edge dynamics can be computed with the following equations:

$$ $\hat{F}\left(\mathbf{x}_{v}(t)\right)=\mathsf{Dec}\left(\phi^{n}\left(\mathsf{ Enc}(\mathbf{x}_{v}(t)),t\right)\right),\quad\hat{G}\left(\mathbf{x}_{v}(t), \mathbf{x}_{u}(t)\right)=\mathsf{Dec}\left(\phi^{e}\left(\mathsf{Enc}(\mathbf{ X}(t)),t\right)\right).$ (4) $$

Unlike Cranmer et al. (), we do not conduct the genetic search on (4) directly because the neural dynamics are not accurate enough for the fitness calculation. Instead, coordinated genetic search uses the neural dynamics for reference and denoised interpolated trajectories for fitness calculation to guide the search process.
In the coordinated genetic search, we have two populations containing symbolic node dynamics and edge dynamics, i.e., $\mathcal{F}$ and $\mathcal{G}$.
However, in our problem setup, the populations $\mathcal{F}$ and $\mathcal{G}$ often have different distances to their corresponding ground-truth dynamics.
Jointly evolving $\mathcal{F}$ and $\mathcal{G}$ may cause the population closer to ground truth to be overfitting and another population to be underfitting. Therefore, we select population to coordinate their evolution process by comparing the distances of $\mathcal{F}$ and $\mathcal{G}$ to the references $\hat{F}$ and $\hat{G}$ in (4).
In the search process, we calculate the distances between populations $\mathcal{F}$ and $\mathcal{G}$ to the references $\hat{F}$ and $\hat{G}$, respectively with

$$ $d(\mathcal{F})=\sum\nolimits_{F\in\mathcal{F}}\|F-\hat{F}\|^{2}\quad\text{and} \quad d(\mathcal{G})=\sum\nolimits_{G\in\mathcal{G}}\|G-\hat{G}\|^{2},$ (5) $$

where $\|\cdot\|$ represents the distance between functions and is the average of absolute error between two functions on randomly sampled points in our paper.
Then, we select to evolve the population with a larger distance to the reference, i.e., if $d(\mathcal{F})>d(\mathcal{G})$, we evolve the node dynamics population $\mathcal{F}$; otherwise, we evolve the edge dynamics population $\mathcal{G}$.
The coordinated strategy in genetic search helps to avoid the overfitted expressions in symbolic space and improve the quality of the symbolic expressions of $F$ and $G$.

To evolve the selected population, we assess the fitness of each symbolic expression in the selected populations by comparing the error between the integral of combined symbolic node dynamics and edge dynamics with the interpolated trajectory. For the node dynamics population $\mathcal{F}$ and edge dynamics population $\mathcal{G}$, the fitness of a symbolic node dynamics $F\in\mathcal{F}$ and a symbolic edge dynamics $G\in\mathcal{G}$ are respectively calculated with

$$ $\displaystyle\!\!\!\!\!\mathsf{f}_{F}\!=\!\text{Mean}\!\circ\!\text{BigK}\! \Big{\{}\sum\nolimits_{v\in V,t\in T}\!-e\!\Big{(}\!\int_{0}^{t}\!\!\!F( \mathbf{x}_{v})+\!\!\!\sum\nolimits_{u\in N_{v}}\!\!\!G(\mathbf{x}_{v},\mathbf {x}_{u})dt,f_{\theta}\left(G,\mathbf{X}(t_{0}),t\right)_{v}\!\!\Big{)}\!\Big{| }G\in\mathcal{G}\!\Big{\}},$ (6) $\displaystyle\!\!\!\!\!\mathsf{f}_{G}\!=\!\text{Mean}\!\circ\!\text{BigK}\! \Big{\{}\sum\nolimits_{v\in V,t\in T}\!-e\!\Big{(}\!\int_{0}^{t}\!\!\!F( \mathbf{x}_{v})+\!\!\!\sum\nolimits_{u\in N_{v}}\!\!\!G(\mathbf{x}_{v},\mathbf {x}_{u})dt,f_{\theta}\left(G,\mathbf{X}(t_{0}),t\right)_{v}\!\!\Big{)}\!\Big{| }F\in\mathcal{F}\!\Big{\}},$ (7) $$

where $e(\cdot,\cdot)$ is the error function between the symbolic integral value and interpolated trajectory, $T$ is the set of timestamps interpolated from training time ranges, BigK is the function that selects the $K$ largest negative errors from a set, Mean is the mean operator and $\circ$ is the composition operator.
In (6) and (7), a large fitness indicates that the corresponding symbolic expression is close to the ground-truth dynamics because the integral of the symbolic expression of dynamics is close to the interpolated trajectory.
Then, the expressions in the selected population are selected, crossed, and mutated by the genetic algorithm to generate the next population.

The coordinated genetic search algorithm is presented in Algorithm 1. In the algorithm, we first initialize the node dynamics and edge dynamics populations $\mathcal{F}^{(0)}$ and $\mathcal{G}^{(0)}$ with random symbolic expressions in line 1. In the evolution process, the populations $\mathcal{F}^{(i-1)}$ and $\mathcal{G}^{(i-1)}$ are evolved to generate the next populations $\mathcal{F}^{(i)}$ and $\mathcal{G}^{(i)}$ in lines 3-14. The process is repeated until the fitness from (6) and (7) is below a threshold (lines 6 & 11) or reaching the maximum number $M$ of iterations. In the end, the algorithm selects and outputs the symbolic node dynamics and edge dynamics with the lowest fitness in lines 16-18.

Figure: Algorithm 1 Coordinated genetic search for symbolic regression

### 4.3 Comparison

We compare SymDL (Cranmer et al., ), NASSymDL (Shi et al., ), D-CODE (Qian et al., ), TP-SINDy (Gao & Yan, ) and PI-NDSR in Table 1.
These methods cater to different problem settings, utilizing distinct forms of input and output.
SymDL and NASSymDL perform general symbolic regression, finding a function $y=f(x)$ from input-output pairs ${(x_{i},y_{i})}_{i=1}^{N}$.
D-CODE focuses on symbolic regression of dynamics, taking a single trajectory ${x(t)|t\in\mathcal{T}}$ to output the governing ODE $\dot{x}=\nicefrac{{dx}}{{dt}}$.
TP-SINDy and PI-NDSR target symbolic regression of complex network dynamics, using multiple trajectories $\left\{\mathbf{X}(t)|t\in\mathcal{T}\right\}$ to output symbolic network dynamics $F$ and $G$.

We compare these algorithms in two aspects: proxy models and formula regression. ‘Proxy models’ are trained to fit data and serve as a basis for deriving symbolic expressions. ‘Formula regression’ directly extract symbolic expressions from raw data or proxy models.
For proxy model, SymDL uses graph networks (GN) with inductive bias, and NASSymDL employs neural architecture search (NAS) for skeleton search. D-CODE can incorporate any suitable regressor. PI-NDSR fits multiple trajectories using PIND, a graph neural ODE aligned with network dynamics for better generalization.
SymDL and NASSymDL rely on potentially noisy estimated derivatives for dynamics fitting, whereas D-CODE and PI-NDSR train directly on raw observations for improved accuracy. TP-SINDy operates without a proxy model.

In formula regression, methods using genetic search employ elementary operations (e.g., $+,-,\times,\div,\sin,\exp$) to represent formula, offering more flexibility and requiring less prior knowledge than TP-SINDy’s linear combination of functions in predefined function library.
SymDL and NASSymDL use the internal functions (Cranmer et al., ) from the proxy models as supervision to compute fitness in genetic search. D-CODE uses the interpolated trajectories as supervision. TP-SINDy is based on sparse regression and uses estimated derivatives for symbolic regression, which can be noisy and inaccurate over large time intervals. PI-NDSR design a coordinated genetic search, using (4) from proxy model as references and interpolated trajectories as supervision.

**Table 1: Comparison with different methods for symbolic regression. (GN: graph network, NAS: neural architecture search, GS: genetic search)**
| Category | Design | SymDL | NASSymDL | D-CODE | TP-SINDy | PI-NDSR |
| --- | --- | --- | --- | --- | --- | --- |
| Input | input-output<br>pairs | input-output<br>pairs | single<br>trajectory | multiple<br>trajectories | multiple<br>trajectories |  |
| Proxy<br>model | Model design | GN w/<br>inductive bias | NAS | any regressor | – | PIND |
| Dynamics<br>fitting data | estimated<br>derivatives | estimated<br>derivatives | raw<br>observations | – | raw<br>observations |  |
| Formula<br>regression | Prior<br>knowledge | elementary<br>operation | elementary<br>operation | elementary<br>operation | function<br>library | elementary<br>operation |
| Method | GS | GS | GS | sparse<br>regression | coordinated<br>GS |  |
| Supervision | internal<br>functions | internal<br>functions | interpolated<br>trajectory | estimated<br>derivatives | network ref &<br>interp. trajectories |  |
| Output | input-output<br>mapping | input-output<br>mapping | ODE | Graph<br>ODE | Graph<br>ODE |  |

## 5 Experiments

In this section, we conduct experiments on both synthetic and real-world datasets to evaluate the PI-NDSR.
All experiments are implemented with PyTorch (Paszke et al., ), PyTorch Geometric (Fey & Lenssen, ), and gplearn (Stephens, ) in NVIDIA GeForce RTX 4090 GPUs and AMD EPYC 7763 Processors.

### 5.1 Experiments on synthetic dataset

#### Baseline

We compare our method with baselines SymDL (Cranmer et al., ), SINDy (Brunton et al., ) and Two-Phase SINDy (TP-SINDy)(Gao & Yan, ). Cranmer et al. () can not be aplied directly to symbolic regression of dynamics, we adopt it as a baseline by first using numerical methods (the same metods as TP-SINDy) to estimate derivatives and then applying their method. SINDy (Brunton et al., ) is a sparse regression methods to find symbolic dynamics. SINDy first numerically estimates the derivative of each node’s activity through the five-point approximation (Sauer, ) and then optimizes the coefficients of the linear combination of predefined orthogonal basis functions. TP-SINDy is an improved version of it, which contains more elementary functions and a extra finetuning phase to remove terms with small coefficients. NASSymDL (Shi et al., ) is not compared because it cannot find a better network architecture than PIND. We also do not compare against D-CODE (Qian et al., ) because it does not have a natural extension to the dynamics regression of multiple trajectories.

#### Dataset

We investigate the following four network dynamics in experiments, i.e., Susceptible-Infected-Susceptible (SIS) Epidemics Dynamics (Pastor-Satorras et al., ), Lotka-Volterra (LV) Population Dynamics (MacArthur, ), Wilson-Cowan Neural Firing Dynamic (Laurence et al., ; Wilson & Cowan, ) and Kuramoto Oscillators(KUR) model (Kuramoto & Kuramoto, ). Their dynamics are shown in Table 2.
We conduct experiments on two complex network structures, i.e., Erdős-Rényi (ER) graph (Erdos & Renyi, ) and Barabási-Albert (BA) graph (Barabási & Albert, ) with 200 nodes.

We randomly initialize the state of all nodes and regularly sample $K$ timestamps ${t_{0},t_{1},\cdots,t_{K-1}}$ from the range $[0,T]$ because all other baselines rely on the equal time interval to estimated time derivatives. Then we simulate the whole dynamics to get the node states $[\mathbf{X}(t_{0}),\mathbf{X}(t_{1}),...,\mathbf{X}(t_{K-1})]$. The edge weight $a_{vu}$ is set to binary values, i.e., $a_{vu}=1$ if there is an edge between node $v$ and node $u$, otherwise $a_{vu}=0$.

**Table 2: Dynamics for generating synthetic dataset.**
|  | node dynamics | edge dynamics |
| --- | --- | --- |
| SIS | $-\delta x_{i}(t)$ | $(1-x_{i}(t))x_{j}(t)$ |
| LV | $x_{i}(t)(\alpha-\theta x_{i}(t))$ | $-x_{i}(t)x_{j}(t)$ |
| WC | $-x_{i}(t)$ | $(1+\exp(-\tau(x_{j}(t)-\mu)))^{-1}$ |
| KUR | $\omega$ | $\sin(x_{i}(t)-x_{j}(t))$ |

#### Evaluation metrics

The performance is evaluated by two metrics. (a) The recovery probability (Rec. Prob.) of formulas with correct skeletons. (See Appendix A for computation details). (b) The mean squared error (MSE) between the simulated trajectories using the recovered symbolic expression of dynamics and the ground truth observations.
For the fair evaluation, we only compute MSE for the symbolic expressions with correct skeletons. So the MSE can measure how accurate the constants in the formulas are.

#### Results

The comparison results are shown in Table 3.
The proposed PI-NDSR generally has a higher recovery probability. For SIS and LV dynamics, TP-SINDy is not stable enough to recover the formula with the correct skeleton.
This may result from the instability in the numerical estimation step for derivatives and
the failure in narrowing down model space because normalized data may cause the overfitting of candidate functions.
For the WC dynamics, the TP-SINDy always fails to regress the correct skeleton, this is because the edge dynamics evolves a parametric function that can not be represented by a linear combination of predefined functions. For the KUR dynamics, both TP-SINDy and PI-NDSR succeed with recovery probability $1$.
SINDy does not contain the finetuning phase which TP-SINDy has. As a result, it exhibits a lower recovery probability compared to TP-SINDy. SymDL also relies on the unstable numerical estimation of derivatives, and its symbolic regression process relies solely on the proxy model while PI-NDSR utilizes additional information from the denoised and augmented trajectories. As a result, PI-NDSR achieves the highest recovery probability.

Even in experiments where both methods successfully produce the correct formula skeleton, PI-NDSR consistently achieves better performance. This is because other baselines include a numerical estimation step for derivatives, which can introduce additional errors and result in inaccuracies in the constants of the formula.

**Table 3: Performance comparison on synthetic datasets. MSE values are scaled by $10^{-2}$ and multiply by $10^{-2}$ to obtain the actual values. (TP: TP-SINDy, PI: PI-NDSR, NA: MSE is not applicable because of failure of the correct skeleton recovery.)**
| Graphs | Dynamics | Rec. Prob.$\uparrow$ | MSE$\downarrow$ ($10^{-2}$) |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| SymDL | SINDy | TP | PI | SymDL | SINDy | TP | PI |  |  |
| BA | SIS | 0.35 | 0.11 | 0.15 | 1 | 0.979±0.173 | 0.484±0.056 | 0.434±0.052 | 0.312±0.012 |
| LV | 0.16 | 0.12 | 0.20 | 1 | 2.075±0.303 | 1.170±0.049 | 0.875±0.057 | 0.136±0.008 |  |
| KUR | 0.80 | 0.87 | 1 | 1 | 0.064±0.018 | 0.175±0.016 | 0.040±0.003 | 0.007±0.001 |  |
| WC | 0.56 | 0 | 0 | 1 | 0.362±0.057 | NA | NA | 0.092±0.004 |  |
| ER | SIS | 0.31 | 0.11 | 0.17 | 1 | 1.173±0.095 | 0.468±0.059 | 0.386±0.051 | 0.119±0.025 |
| LV | 0.15 | 0.09 | 0.19 | 1 | 1.784±0.236 | 0.941±0.041 | 0.763±0.077 | 0.251±0.007 |  |
| KUR | 0.87 | 0.78 | 1 | 1 | 0.071±0.018 | 0.087±0.022 | 0.069±0.019 | 0.017±0.001 |  |
| WC | 0.40 | 0 | 0 | 1 | 0.266±0.047 | NA | NA | 0.044±0.003 |  |

### 5.2 Experiments on real dataset

#### Dataset

We demonstrate the effectiveness of PI-NDSR on the real epidemic network.
We use the same influenza A (H1N1) spreading dataset as (Gao & Yan, ).
In this dataset, each node represents a country or region, with the daily counts of newly reported cases serving as the state of these nodes.
The edges of the complex network are defined by the global aviation routes, depicting human mobility between regions.
Our goal is to uncover the dynamics that govern the spread of the disease.
To ensure a fair comparison, we employed the same data preprocessing procedures as (Gao & Yan, ),
such as constructing the adjacency matrix and cleaning the data.

#### Results

We use PI-NDSR and TP-SINDy to regress the symbolic expression of dynamics for the spreading of influenza A.
The symbolic expressions of the spreading of influenza A regressed by PI-NDSR is

$$ $\dot{\mathbf{x}}_{v}(t)=a\mathbf{x}_{v}(t)+\sum\nolimits_{u\in N_{v}}\frac{b}{ 1+\exp{-\left(m\mathbf{x}_{v}(t)+c\right)}}\mathbf{x}_{u}(t),$ (8) $$

where $a=0.0740$, $b=0.0015$, $m=-0.0041$ and $c=9.9643$.
Node dynamics in (8) is a linear function, which aligns with the exponential growth of the epidemic.
Edge dynamics in (8) is proportional to the neighboring region’s state, which is consistent with the fact that the epidemic spreads increases with the number of infected cases in neighboring regions.
The other factor of edge dynamics consists of a composition of a linear transformation followed by a sigmoid activation.
This indicates that the infection from neighboring regions may gradually decrease as the number of infected cases in the current region increases.
This trend may be caused by the reduction of the willingness of people to travel to epidemic areas or the decrease of basic reproduction number $(R_{0})$ when the density of infected people is high.

The symbolic expression of the spreading of influenza A regressed by TP-SINDy is

$$ $\dot{\mathbf{x}}_{v}(t)=a^{\prime}\mathbf{x}_{v}(t)+\sum\nolimits_{u\in N_{v}} \frac{b^{\prime}}{1+\exp{-\left(\mathbf{x}_{v}(t)-\mathbf{x}_{u}(t)\right)}},$ (9) $$

where $a^{\prime}=0.074$ and $b^{\prime}=7.130$. (9) fails to capture the spreading pattern of the epidemic.
When there are no infected cases in the complex network, the edge dynamics of (9) will still result in a non-zero growth rate, which is not reasonable from a physical perspective.
In the same scenario, (8) yields zero growth rate, which is consistent with the fact that the epidemic will not spread when there are no infected cases.

Figure: Figure 2: Visualizing the predicted number of newly reported cases in two regions using symbolic expressions from TP-SINDy and PI-NDSR.
Refer to caption: x2.png

We compare the trajectories of symbolic expressions in (8) and (9) with the real infected cases.
Fig. 2 visualizes the simulation results of inferred dynamics in two regions, i.e., “Finland” and “Saint Pierre and Miquelon”.
The trajectories of the infected cases in the two regions inferred by PI-NDSR are consistent with the ground truth, while the trajectories inferred by TP-SINDy deviate from the ground truth.

We quantitatively evaluate the errors of regressed dynamics. As the scale of infected cases varies across different regions, we normalize the infected cases to the range of $[0,1]$ by the maximal value of each region. The mean square error (MSE) of TP-SINDy is $0.9028$, while the MSE of PI-NDSR is $0.8261$. The results show that PI-NDSR fits the neural dynamics better.

### 5.3 Robustness on PI-NDSR and TP-SINDy

Figure: (a) Recovery probability versus noise level
Refer to caption: x3.png

We adopt the evaluation settings in
(Qian et al., ) to test the robustness of our method. We evaluate the performances when observations are noisy or the time interval is large.
We select KUR dynamics to evaluate the robustness of PI-NDSR and TP-SINDy.

Firstly, we add Gaussian noise to all node states to evaluate its performance under noise.
The magnitude of noise is measured by the signal-to-noise ratio (SNR). The results of recovery probability are presented in Fig. 3(a). In the experiments, as the signal-to-noise ratio (SNR) decreases from 70 dB to 25 dB, our method consistently maintains a 100% recovery rate. In contrast, TP-SINDy’s recovery rate drops to 0% when the SNR reaches 30 dB.
In Fig. 3(b), the proposed method consistently produces more accurate symbolic expressions that have lower MSE.
The reason is that TP-SINDy relies on numerically estimating time derivatives that are noisy and inaccurate. On the other hand, our method uses neural dynamics to denoise and interpolate observations directly.
Deep neural networks can handle noisy observations effectively because they have the ability to learn and extract meaningful patterns from large amounts of data, even when the data contains significant noise.
Using the accurately denoised observations, PI-NDSR can predict constants in the formula better and produce a more accurate trajectory when noise exists.

TP-SINDy relies on the equal time interval to estimate time derivatives.
So we increase the size of time intervals to evaluate the performances of both TP-SINDy and PI-NDSR.
We sample the timestamps regularly from $[0,100]$ with different intervals $\Delta t$. As shown in Fig. 3(c),
our model achieved a $100\%$ recovery probability across all values of $\Delta t$, whereas TP-SINDy consistently failed to recover the correct formula skeleton when the time interval is large.
Fig. 3(d) shows that PI-NDSR always produces more accurate results when both methods produce the correct skeleton of dynamics. This advantage arises because the interpolated observations in PI-NDSR are better suited when the time interval is large compared to the estimated time derivatives used by TP-SINDy.
Visualization of interpolated trajectories and estimated time derivatives are shown in Fig. 5 of the appendix.

### 5.4 Ablation study

We conduct ablation studies to demonstrate the necessity of the components in PI-NDSR. We test two variants of PI-NDSR: (a) Using the original observations instead of interpolated and denoised trajectory when doing the coordinated genetic search. (without Interp.) (b) Using the interpolated trajectory but evolving both populations $\mathcal{F}$ and $\mathcal{G}$ every time. (without Coord.)

Table 4 shows the results of SIS and LV dynamics in the BA graph. Without the interpolated and denoised observations,
both the recovery probability and the accuracy of PI-NDSR drop.
This indicates that the interpolated and denoised trajectories can provide high-quality supervision data for symbolic regression. Without the coordinated genetic search, the recovery probability of PI-NDSR drops significantly, and when the skeleton can be successfully recovered the MSE increases slightly, indicating that the search strategy can avoid overfitted expressions in symbolic space and help find correct dynamics expressions.

**Table 4: Ablation study with experiment results on SIS and LV dynamics in BA graph.**
| Model | SIS | LV |  |  |
| --- | --- | --- | --- | --- |
| Rec. Prob.$\uparrow$ | MSE$\downarrow$ ($10^{-2}$) | Rec. Prob.$\uparrow$ | MSE$\downarrow$ ($10^{-2}$) |  |
| PI-NDSR | 1 | 0.312±0.012 | 1 | 0.136±0.008 |
| PI-NDSR(w/o Interp.) | 0.81 | 0.408±0.027 | 0.86 | 0.588±0.028 |
| PI-NDSR(w/o Coord.) | 0.31 | 0.326±0.015 | 0.47 | 0.142±0.014 |

## 6 Conclusion

In this paper, we propose physically inspired symbolic regression to learn symbolic expressions of complex network dynamics from data.
Our method is based on neural dynamics to augment and denoise the trajectory of complex network, and then apply the coordinated genetic search to infer the symbolic expressions based on the dynamics reference from neural dynamics.
In the future, we plan to extend our method to learn the symbolic expression of more complex systems and apply it to more real-world applications. See limitations in Appendix C.

## Appendix A Details on experiments

### A.1 Dataset statistics

The BA graph is generated with 200 nodes and the initial degree of each node is set to 3. The ER graph is generated with 200 nodes and the probability for edge creation is set to 0.02.
The initial states of SIS, LV, and WC dynamics are generated by randomly sampling from $[0,1]$. For KUR dynamics, the initial states are generated by randomly sampling from $[0,2\pi]$.
For SIS dynamics, we set $\delta=0.5$. For LV dynamics, we set $\alpha=0.75$, $\theta=0.5$,. For KUR dynamics, we set $\omega=0.75$. For WC dynamics, we set $\tau=0.75$, $\mu=0.5$. We regularly sample 100 timestamps from $[0,1]$ and simulate the dynamics to generate the observation data.

### A.2 Details for network training

We split the timestamps randomly into training, validation, and testing sets with a ratio of 0.8, 0.2, 0.1 to train the NerualODE. We train the neural dynamics for 1000 epochs using optimizer AdamW. The learning rate is searched in the range of $[1e-3,1e-2]$, the weight decay is set to $0.001$.
We use MLPs as the encoder and decoder of neural dynamics. The hidden dimension of the neural dynamics is set to 10. The details of the network structure are shown in Table 5.

**Table 5: Details of network structure for different dynamics.**
|  | SIS | LV | KUR | WC | real dataset |
| --- | --- | --- | --- | --- | --- |
| Hidden dimension | 10 | 10 | 10 | 10 | 10 |
| Activation of $\phi^{n}$ | ReLU | ReLU | ReLU | ReLU | Sigmoid |
| Activation of $\phi^{e}$ | ReLU | Tanh | Tanh | Sigmoid | Sigmoid |
| Activation of Encoder | ReLU | Tanh | Tanh | ReLU | Tanh |
| Activation of Decoder | ReLU | Tanh | Tanh | ReLU | Tanh |
| Layer of $\phi^{n}$ | 2 | 2 | 1 | 1 | 2 |
| Layer of $\phi^{e}$ | 2 | 2 | 3 | 2 | 3 |

### A.3 Details of genetic search

We implement the coordinated genetic search based on gplearn (Stephens, ). gplearn (Stephens, ) represent the symbolic expressions as a syntax tree, where the functions are interior nodes, and the variables and constants make up the leaves. Evolution such as crossover, mutation, and reproduction are performed on the syntax tree.
The population size of $\mathcal{F}$ and $\mathcal{G}$ are set to 200. The maximum generation of the genetic search $M$ is set to 50 and the stopping threshold $\epsilon=10^{-5}$. The $K$ in Algorithm 1 equals to 20. The function set includes addition, subtraction, multiplication, division, sine, cosine, and exponential. The constants are constrained in the range $[-1,1]$.
Other hyperparameters of gplearn are set as: p_crossover=0.6, p_subtree_mutation=0.1, p_hoist_mutation=0.05, p_point_mutation=0.1, parsimony_coefficient=0.01. We conduct the genetic search in 256 parallel threads to speed up the search process. Our CPUs are two AMD EPYC 7763 Processors.

### A.4 Computational details of Rec. Prob.

The recovery probability is calculated as the ratio of the number of successful recovery of formula skeletons to the total number of experiments. We automatically check the correctness of the recovered formula skeletons using the method for verifying skeletons provided in Qian et al. (). Basically, we replace the constants in the formulas with placeholders and use the $\mathsf{simplify}(f^{\prime}-f)==0$ criterion from the Sympy package to determine if the skeleton is correct.

### A.5 Discussions on the choice of metrics

#### Compute MSE between trajectories instead of the constants.

We do not directly compute the MSE between predicted and true constants. This is because our goal is to evaluate how well the obtained symbolic expressions predict trajectories, which is crucial for real-world scenarios like epidemic forecasting. Directly computing constant errors is insufficient, as different constants impact the trajectory differently. Some constants require high precision, with small deviations causing significant errors, while others are less critical and can tolerate some errors.

#### Compute MSE for formulas with correct skeletons.

For simulated datasets, we choose MSE restricted to correctly recovered skeletons because the baseline methods often exhibit large MSE when recovering incorrect skeletons. Filtering out these formulas allows the baselines to achieve comparable performance. For real datasets, since the true dynamics skeleton is unknown, we directly compare the MSE of the trajectories without filtering by the skeletons.

## Appendix B Additional results

### B.1 Ablation study

We also experiment on the robustness of the ablation variants. Table 6 shows the results of KUR dynamics in the BA graph when the observations are noisy or the time interval is large. Different from the results in Table 4, the success Prob. significantly drop when removing the interpolation part. This proves the effectiveness of neural dynamics in denoising and augmenting trajectories.

**Table 6: The robustness of two variants compared with full method on KUR dynamics in BA graph.**
| Models | Noise (SNR=35dB) | Time interval ($\Delta t=1.28$) |  |  |
| --- | --- | --- | --- | --- |
| Rec. Pro.$\uparrow$ | MSE$\downarrow$ ($10^{-2}$) | Rec. Pro.$\uparrow$ | MSE$\downarrow$ ($10^{-2}$) |  |
| PI-NDSR | 1 | 0.478±0.103 | 1 | 0.454±0.213 |
| PI-NDSR(w/o Interp) | 0.84 | 6.970±1.870 | 0.78 | 2.645±1.138 |
| PI-NDSR(w/o Alter) | 0.76 | 0.504±0.094 | 0.66 | 0.570±0.241 |

### B.2 Runtime

In Table 7, PI-NDSR saves 30.1% running time on SIS dynamics and 39.0% running time on LV dynamics compared with PI-NDSR(w/o Alter). The results show that the coordinated genetic search can significantly reduce the search space and improve the efficiency of the search process.

**Table 7: The runtime (minutes) of PI-NDSR and PI-NDSR(w/o Coord.).**
| Model | SIS | LV |
| --- | --- | --- |
| PI-NDSR | 61.5 | 50.9 |
| PI-NDSR(w/o Coord.) | 88.0 | 83.4 |

### B.3 Examples of expressions from symbolic expression

In this section, we provide examples of symbolic expressions of PI-NDSR, TP-SINDy (Rec.), and TP-SINDy (Fail) on SIS, LV, KUR, and WC dynamics in the BA graph.
TP-SINDy (Rec.) represents the symbolic expressions of TP-SINDy when the skeleton of the dynamics is successfully recovered, while TP-SINDy (Fail) represents the symbolic expressions of TP-SINDy when the skeleton of the dynamics is not successfully recovered.
The expressions are shown in Table 8.

**Table 8: Symbolic regressions of PI-NDSR, TP-SINDy (Rec.), and TP-SINDy (Fail) on SIS, LV, KUR, and WC dynamics in the BA graph.**
| Dynamics | Models | Node dynamics | Edge dynamics |
| --- | --- | --- | --- |
| SIS | GT | $-0.5x_{i}(t)$ | $(1-x_{i}(t))x_{j}(t)$ |
| PI-NDSR | $-0.48540x_{i}(t)$ | $(1-x_{i}(t))x_{j}(t)$ |  |
| TP-SINDy (Rec.) | $-0.46640x_{i}(t)$ | $(0.99119-1.09637x_{i}(t))x_{j}(t)$ |  |
| TP-SINDy (Fail) | $-0.47256x^{2}_{i}(t)-0.12596$ | $0.10860\text{sigmoid}(x_{j}(t)-x_{i}(t))+$<br>$0.18835(x_{j}(t)-x^{2}_{i}(t))$<br>$0.19917(x_{j}(t)-x_{i}(t))+0.35416\sin(x_{j}(t))$ |  |
| LV | GT | $x_{i}(t)(0.75-0.5x_{i}(t))$ | $-x_{i}(t)x_{j}(t)$ |
| PI-NDSR | $x_{i}(t)(0.75034-0.48812x_{i}(t))$ | $-0.99428x_{i}(t)x_{j}(t)$ |  |
| TP-SINDy (Rec.) | $x_{i}(t)(0.69882-0.41853x_{i}(t))$ | $-0.91701x_{i}(t)x_{j}(t)$ |  |
| TP-SINDy (Fail) | $0.03984+0.36330*\sin(x_{i}(t))$ | $-0.945810x_{i}(t)x_{j}(t)-0.11895x_{i}(t)x^{2}_{j}(t)$ |  |
| KUR | GT | $0.75$ | $\sin(x_{i}(t)-x_{j}(t))$ |
| PI-NDSR | $0.75002$ | $\sin(1.0001x_{i}(t)-x_{j}(t))$ |  |
| TP-SINDy (Rec.) | $0.75014$ | $0.99899\sin(x_{i}(t)-x_{j}(t))$ |  |
| TP-SINDy (Fail) | NA | NA |  |
| WC | GT | $-x_{i}(t)$ | $\text{sigmoid}(-0.75(x_{j}(t)-0.5))$ |
| PI-NDSR | $-x_{i}(t)$ | $\text{sigmoid}(-0.74503(x_{j}(t)-0.49128))$ |  |
| TP-SINDy (Rec.) | NA | NA |  |
| TP-SINDy (Fail) | $-0.82267x_{i}(t)$ | $0.08513\text{sigmoid}(x_{j}(t)-x_{i}(t))$<br>$+0.68484\text{sigmoid}(x_{j}(t))$ |  |

### B.4 Example of overfitting

Take the SIS dynamics in the BA graph as an example. The symbolic expressions of PI-NDSR, TP-SINDy (Rec.), and TP-SINDy (Fail) are shown in Table 8. We compute the MSE of the predicted trajectories under interpolated and extrapolated settings. The results are shown in Table 9. Although the symbolic expressions of TP-SINDy (Rec.) and TP-SINDy (Fail) have relatively low MSE values under the interpolated setting, their MSE values are much higher under the extrapolated setting. This indicates that the symbolic expressions of TP-SINDy are overfitted and cannot generalize well to extrapolated setting. We double the time range from $[0,1]$ to $[0,2]$ to evaluate the extrapolation performance.

**Table 9: The MSE of PI-NDSR, TP-SINDy (Rec.), and TP-SINDy (Fail.) under interpolated and extrapolated settings on SIS dynamics in the BA graph.**
|  | PI-NDSR | TP-SINDy (Rec.) | TP-SINDy (Fail.) |
| --- | --- | --- | --- |
| Interpolation | $3.2\times 10^{-3}$ | $5.2\times 10^{-3}$ | $147.2\times 10^{-3}$ |
| Extrapolation | $3.9\times 10^{-3}$ | $46.9\times 10^{-3}$ | $697.0\times 10^{-3}$ |

Note that all failure cases in Table 8 can also be viewed as examples of overfitting. The symbolic expressions of PI-NDSR are more interpretable and simpler while the overfitted symbolic expressions of TP-SINDy are more complex and contain more terms.

### B.5 Visualization of neural dynamics

#### Node/edge dynamics estimation

With the inductive bias, the neural node dynamics and edge dynamics computed with (4) are accurate enough to be used as the neural references. We visualize the result for LV dynamics in the BA graph in Fig. 4. It can be seen that the neural estimations of dynamics are very close to the ground truth.

Figure: (a) Node dynamics. $\hat{F}$ is in red, ground truth is in black.
Refer to caption: x7.png

#### Observation denoising and interpolating

When observations are noisy or time interval is large, the neural dynamics can denoise and interpolate the observations to provide high-quality supervision data for symbolic regression. On the other hand, the numerical estimation is sensitive to noise and needs the sample interval to be small enough. We visualize the interpolate trajectories and the estimated time derivatives in Fig. 5, which is consistent with our contributions.

Figure: (a) Interpolated trajectory with noisy observation.
Refer to caption: extracted/5922682/img/regressed_vs_gt_vs_noise.png

## Appendix C Limitations

Our method has several limitations. Firstly, PI-NDSR cannot successfully recover symbolic expressions when the formulations are highly complex or the dimension of node states is high.
Secondly, our method cannot deal with the complex network dynamics when some variables are missing in the observations.

## Peer Reviews

confidence: 3
contribution: 2
initial_rating: 6
presentation: 2
questions: Did I understand correctly that to train PIND you use four networks? An Encoder, a Decoder, $\phi^e$, $\phi^n$? I believe section 4.1 could be slightly expanded to show how the neural dynamics is trained, maybe supported with a figure.

I also found these typos:
In section 5.1, in "baseline": metods should be methods. "can not be aplied directly" should be "cannot be applied directly".
In Appendix A2, "NerualODE" should be "NeuralODE"
soundness: 3
strengths: The results on these datasets show the good performance of the method, with some proof of robustness too. This applies to the real world dataset too, where the extracted dynamics of this method don't have "un-physical" properties. The ablation studies show that both parts of the algorithm are necessary for the best performance
summary: The paper introduces Physically Inspired Neural Dynamics Symbolic Regression (PI-NDSR), a new method for symbolic regression of complex network dynamics. The approach has two main components: neural dynamics, which denoises and augments observed data, and a coordinated genetic search, which leverages the previously build neural references to optimize the search for symbolic expressions. Evaluated on synthetic an real-world datasets, PI-NDSR outperforms previous methods, offering more accurate symbolic representations and robustness against noise and sampling irregularities.
weaknesses: It is not clear if this method could scale to real world datasets that have more complex dynamics, since there is only a single real-world dataset.

---

confidence: 2
contribution: 3
initial_rating: 6
presentation: 3
questions: (1)  Have the authors attempted to explore solutions for the missing data problem?

(2) The paper mentions that the computational cost of PI-NDSR is relatively high. Could the authors provide some suggestions for reducing this computational cost?
soundness: 3
strengths: **Novel Approach** The use of genetic programming for symbolic regression of complex network dynamics is innovative as it constrains the symbolic search space, thereby helping to mitigate overfitting. Moreover, this approach does not rely on a predefined function library.

**Improved Performance** PI-NDSR outperforms existing methods in terms of recovery probability and error, as demonstrated by experiments on various datasets.

**Robustness** By incorporating neural ODEs, the method demonstrates strong robustness to noise and can effectively handle large time intervals between observations.
summary: This paper presents PI-NDSR, a novel method for symbolic regression of complex network dynamics. It leverages neural networks for trajectory interpolation and denoising, and combines it with a coordinated genetic search algorithm that utilizes references from the neural dynamics.
weaknesses: **Limited Novelty in Applying NeuralODE for Interpolation** The use of neural networks to model differential equations for interpolation is not particularly novel and has been addressed in prior research. With respect to the interpolation method, the paper does not introduce any methodological innovations.

---

confidence: 3
contribution: 2
initial_rating: 6
presentation: 3
questions: 1. **Can Interpolation Truly Address Data Sparsity?**
   - The paper uses neural network-based interpolation and denoising to handle issues related to data sparsity and noise in complex network dynamics. However, interpolation fundamentally relies on existing observed data, meaning that any missing information between observation points could still lead to deviations from the true dynamics. Particularly in cases of high data sparsity, interpolation methods may produce overly smooth results, potentially missing critical dynamic variations. Thus, interpolation alone may not fully resolve the information loss caused by sparsity, as its effectiveness remains dependent on the quality and density of the original data. This could be further addressed by exploring additional complementary data sources or methods that can better capture system dynamics under sparse conditions.

2. **Unclear Motivation in the Paper’s Methodology**
   - The paper highlights challenges in symbolic regression for complex network dynamics, such as noise, data sparsity, and the need for decoupling node and edge dynamics. However, it does not fully explain the specific reasons and impacts of these challenges. In particular, the motivation for decoupling node and edge dynamics lacks sufficient detail, making it difficult for readers to understand why this decoupling is critical to obtaining high-quality symbolic expressions. Simply mentioning that "traditional methods tend to lead to overfitting in symbolic space" is somewhat vague and lacks specific examples or theoretical support. Clarifying the distinct roles of node and edge dynamics in symbolic expressions, and why decoupling is necessary for accurate representation, would enhance the paper's persuasiveness.
soundness: 2
strengths: **Addresses a Novel Challenge of Data Sparsity and Noise in Symbolic Regression**: The PI-NDSR method effectively tackles the issue of sparse and noisy data in symbolic regression by introducing a neural dynamics module that performs interpolation and denoising, making it suitable for complex network data with inconsistent quality.

 **Improves Genetic Algorithm for Symbolic Regression**: By incorporating a coordinated genetic search that uses neural dynamics as a reference, PI-NDSR refines the traditional genetic algorithm approach, allowing for more accurate and efficient search of symbolic expressions while avoiding overfitting.

**Demonstrates Promising Performance**: The method shows competitive performance in both synthetic and real-world datasets, achieving higher recovery accuracy and lower error rates than baseline methods, indicating its potential as a robust solution for complex network dynamics modeling.
summary: This paper introduces a Physically Inspired Neural Dynamics Symbolic Regression (PI-NDSR) method for symbolic regression of complex network dynamics. The motivation behind this research is to address the challenges of symbolic regression in complex networks, particularly the issues of noisy, sparse observations and the need to decouple node dynamics from edge dynamics. PI-NDSR features two main innovations: (1) a neural dynamics module (PIND) that interpolates and denoises the observed data to generate dense, low-noise trajectories, thereby avoiding the instability of directly modeling derivatives; and (2) a coordinated genetic search algorithm that uses neural dynamics as a reference to guide symbolic regression of node and edge dynamics separately, preventing overfitting and improving search efficiency. Experiments on synthetic and real-world epidemic data demonstrate that PI-NDSR outperforms traditional methods in recovery accuracy and error rates, showing enhanced robustness in high-noise and sparse-sampling scenarios.
weaknesses: **Lack of Novelty in Using Neural Networks for Differential Equation Interpolation**: Modeling differential equations with neural networks to achieve interpolation is not particularly novel and has been explored in previous works. This may limit the perceived innovation in the proposed approach’s foundational methodology.

**Potential Oversmoothing from Interpolation without Fully Addressing Data Sparsity**: While the method employs interpolation to handle data sparsity, it does not adequately address how to mitigate potential oversmoothing effects in the interpolated results. This could mean that the sparsity issue remains partially unresolved, as the interpolation might overlook important variations in sparse data.

**Unclear Motivation for Decoupling Node and Edge Dynamics**: The rationale for decoupling node and edge dynamics is not clearly explained. The claim that "traditional symbolic regression methods struggle to find high-quality symbolic expressions without decoupling node and edge dynamics, leading to overfitting in symbolic space" is somewhat vague and lacks a concrete justification, potentially weakening the motivation behind this design choice.

---

confidence: 3
contribution: 2
initial_rating: 6
presentation: 2
questions: - What could be the potential issues when using the model for highly complex formulations?

- For simulated datasets, would it be better to report both filtered and unfiltered MSE, or to use a separate metric for skeleton recovery?

- It would be beneficial to expand the discussion on the limitations of the model.
soundness: 3
strengths: - The model introduced in this paper is groundbreaking, offering a method designed to overcome limitations of other approaches, such as high sensitivity to noise.

- Comprehensive experiments are provided, comparing the proposed model with established methods across synthetic and real-world datasets.

- The paper is generally well-written (except for the conclusion) with a clear and linear structure that makes it easy to follow the proposed methodology, experiments, and results. The explanations of the Physically Inspired Neural Dynamics Symbolic Regression framework and the challenges it addresses are presented in a logical sequence, enhancing reader comprehension.
summary: The paper introduces Physically Inspired Neural Dynamics Symbolic Regression, a method for automatically discovering symbolic expressions that describe complex network dynamics. The authors present a two-part solution: a Physically Inspired Neural Dynamics component that augments and denoises trajectory data through interpolation, and a coordinated genetic search algorithm that derives symbolic expressions using references from the neural dynamics. The method addresses key limitations of existing approaches by handling noisy observations from multiple trajectories without requiring pre-defined function libraries or estimated time derivatives. The authors demonstrate PI-NDSR's effectiveness through evaluation on both synthetic datasets with various dynamics and real datasets on disease spreading, showing improvements in recovery probability and error rates compared to existing methods.
weaknesses: - The evaluation of the model is somewhat lacking, particularly in the choice of metrics for synthetic and real datasets, which is not fully explained. For example, filtering out models that don’t fully recover the skeleton in synthetic datasets may overlook valuable insights into how different methods handle structural errors. Often, a method’s ability to approximate the skeleton, even imperfectly, can still yield good performance on downstream tasks.

- The conclusion is brief, as it mainly summarizes the process without emphasizing key findings. Improvements over existing methods are not highlighted, and it lacks an impact statement or discussion of broader implications. The future work section is vague; it’s unclear what is meant by “more complex systems” and “more real-world applications.” Although limitations are mentioned, it’s not evident why the method fails for highly complex formulations.

---

confidence: 4
contribution: 2
initial_rating: 3
presentation: 2
questions: -Please clarify the design rationale for the encoder and decoder, and explain how to ensure the decoder maintains consistency between training and inference. I suggest that the authors provide empirical evidence or theoretical justification for the effectiveness of the encoder-decoder design.

-Can the method proposed in the paper recover the symbolic expressions of network dynamics when the state of nodes has multiple features? Please either provide results for multi-dimensional node states, or discuss the potential challenges and necessary modifications for extending the method to such cases.

-As described in line 749 of the appendix, Rec.Prob is computed by “simplify,” which suggests that Rec.Prob takes the value 0 or 1. Why is there a float number in Table 3? Please clarify how to calculate the Rec.Prob values in Table 3, and explain why some values are not binary if the calculation method described in the appendix is correct.

-What is the specific meaning of the "error" in line 215 and how is it calculated?

-Lack of discussion on highly relevant work. DNND [1] employs the GNN proxy model to decompose the dynamics and then uses the SINDy to obtain the symbolic expression. Please explain the difference between the method in the paper and DNND and suggest adding comparisons. Please include a comparison with DNND in experiments, or at minimum, to provide a detailed discussion of how the proposed method differs from and potentially improves upon DNND.

-Due to the method's design in this paper being similar to the graph ODE models, the proposed method appears to be able to handle irregularly-sampled data, which could have been considered an advantage but was not mentioned and validated in this paper. I suggest that the authors conduct additional experiments with irregularly-sampled data to validate this capability, or discuss why this potential advantage was not explored in the current work.

-There are some typos and unclear points as follows: 
    The inputs in Eq. (2) and Eq. (4) are inconsistent. 
    In Eq. (6) and (7), too many G's appear, and these G's seem to have different meanings.

[1] Liu, B., Luo, W., Li, G., Huang, J., Yang, B.: Do we need an encoder-decoder to model dynamical systems on networks? In: Proceedings of the Thirty-Second International Joint Conference on Artificial Intelligence, pp. 2178–2186 (2023)
soundness: 2
strengths: (1)The proposed method can handle observation data with noise for inferring interpretable network dynamics.
(2)The proposed method circumvents the need for expert knowledge elaboration, necessitating only the basic operator regression symbolic model.
summary: This paper proposes a physically-inspired neural dynamics symbolic regression (PI-NDSR) to attempt to address the limitations of sparse regression methods for nonlinear network dynamics that require finite differences and expert design function libraries. It comprises two principal components: a physically-inspired neural dynamic model that enhances trajectories and reduces noise by observing trajectory interpolation; and a coordinated genetic search algorithm to derive symbolic expressions.
weaknesses: (1)The rationality of the method design needs to be further clarified. The design and assumptions of the encoder and decoder especially need to be explained. When training the neural network, the input of the Dec in Eq. 2 is \phi^n + \phi^e, but when inferring equations, the input of the Dec becomes \phi^n or \phi^e. This means that the Dec needs to satisfy Dec(\phi^n + \phi^e)=Dec(\phi^n)+Dec(\phi^e), which is a strong assumption for the design of the decoder, thus affecting the corresponding inverse operation, i.e. the encoder, which also satisfies the corresponding assumption. There is reason to doubt the effectiveness and rationality of the encoder and decoder.
(2)The contribution of the proposed method lacks clear clarification.
(3)The experimental analysis is not comprehensive. For example, there is a lack of time cost analysis while time consumption is an important metric for algorithm performance evaluation and the states of all network dynamics evaluated in the paper are one-dimensional.

---

confidence: 4
contribution: 2
initial_rating: 5
presentation: 3
questions: 1. How does the computational efficiency of GA compare to TP-SINDy in the symbolic regression component? Please include an appropriate discussion on this part. If the method in this paper takes too long, using a predefined library, as in TP-SINDY, remains a good option. Determining which approach is better requires evaluating the time consumption of each method.
soundness: 2
strengths: 1. The author proposes a network dynamics symbolic regression method that uses neural networks to complete data, which can help mitigate the problem of noisy data trajectories.  
2. By applying genetic programming for symbolic regression, this method avoids the need for a pre-designed knowledge base.
summary: This paper introduces Physically Inspired Neural Dynamics Symbolic Regression (PI-NDSR), a technique that combines neural networks with genetic programming to autonomously discover symbolic representations of dynamics. PI-NDSR includes two main components: a Physically Inspired Neural Dynamics (PIND) module, which enhances and denoises trajectories 
through observed trajectory interpolation; and a coordinated genetic search algorithm that extracts symbolic expressions. This algorithm uses references from both node and edge dynamics provided by neural dynamics to prevent overfitting in symbolic space.
weaknesses: 1. The innovation of the method in the article is limited. First, the use of neural networks to fit data in order to enhance symbolic regression performance has been widely adopted in many works. For instance, in AI-Feynman[1], a neural network is first used to fit the data, and then the fitted data is analyzed to infer the nature of the equation, which facilitates subsequent symbolic regression. Similarly, this paper also uses a neural network to fit data to avoid inaccuracies in derivative estimation caused by sparse or noisy observational data. Essentially, both approaches leverage the generalization capability of neural networks to mitigate the negative impact of sparse data points. Second, using surrogate models for symbolic regression is also a common approach. For example, SymDL [2], like this paper, uses an MLP locally as a surrogate model for node dynamics G and edge dynamics F.  
2. The experiments in the article are not sufficiently comprehensive; they are limited to relatively simple dynamical scenarios and have not been applied to more complex scenarios such as neuronal, genetic, social, and coupled oscillator dynamics. Moreover, the experiments are restricted to one-dimensional data, while many other works, like TP-SINDy, can handle higher-dimensional data. For example, the FHN dynamics is a higher-order network equation, and this paper does not include experimental discussion on this equation.  


References:  
[1] Silviu-Marian Udrescu, Max Tegmark ,AI Feynman: A physics-inspired method for symbolic regression.Sci. Adv.6,eaay2631(2020).DOI:10.1126/sciadv.aay2631.   
[2] Cranmer M, Sanchez Gonzalez A, Battaglia P, et al. Discovering symbolic models from deep learning with inductive biases[J]. Advances in neural information processing systems, 2020, 33: 17429-17442.  
[3] Gao, TT., Yan, G. Autonomous inference of complex network dynamics from incomplete and noisy data. Nat Comput Sci 2, 160–168 (2022). https://doi.org/10.1038/s43588-022-00217-0.

Meta-review: The paper proposes a framework for symbolic regression of network dynamics using graph neural networks and genetic algorithms (for extracting a symbolic representation).

The paper had a large amount of reviews (6), which scores ranging from 3 to 6. After the rebuttal phase, the most negative score improved to 5, with 2 reviewer negative on the current status of the paper. They both raise a few concerns that should be addressed. In general, there was limited interaction in the rebuttal stage, and as a result some concerns remained open.

There is a consensus that the paper's novelty is unclear, and that the design of the method requires more details (especially concerning the assumptions that are made). The authors have commented that the use of a genetic algorithm conditioned on the specific neural formulation is novel enough. All reviewers who have interacted during the rebuttal have remained mostly negative.

Due to the unaddressed concerns and the generally unsatisfactory rebuttal, I have leaned towards a negative evaluation of the paper, although I mention this is a borderline case which could have benefitted from a stronger discussion.

Decision: Reject