Discovering equations from data: symbolic regression in dynamical systems

TL;DR

This paper reviews symbolic regression methods for discovering equations from data, compares five state-of-the-art algorithms on various dynamical systems, and identifies PySR as the most effective tool for inferring accurate equations.

Contribution

It provides a comprehensive overview of symbolic regression techniques and benchmarks their performance, highlighting PySR's superior ability to recover governing equations.

Findings

PySR outperforms other methods in recovering equations.

Some estimates match original analytical forms closely.

Symbolic regression shows promise for modeling real-world phenomena.

Abstract

The process of discovering equations from data lies at the heart of physics and in many other areas of research, including mathematical ecology and epidemiology. Recently, machine learning methods known as symbolic regression emerged as a way to automate this task. This study presents an overview of the current literature on symbolic regression, while also comparing the efficiency of five state-of-the-art methods in recovering the governing equations from nine processes, including chaotic dynamics and epidemic models. Benchmark results demonstrate the PySR method as the most suitable for inferring equations, with some estimates being indistinguishable from the original analytical forms. These results highlight the potential of symbolic regression as a robust tool for inferring and modeling real-world phenomena.

Tables

Table 1. Table 1: All compartmental models used and their governing equations, with transition rates as indicated in Figure 3 .

SIS	SIR	SIRV	SIRS	SEIR	SEIRD
$\dot{I}=-\dot{S}=\frac{\beta }{N}IS-\gamma I$	$\dot{S}=-\frac{\beta IS}{N}$ $\dot{I}=\frac{\beta IS}{N}-\gamma I$ $\dot{R}=\gamma I$	$\dot{S}=-\frac{\beta SI}{N}-ϵS$ $\dot{I}=\frac{\beta SI}{N}-\gamma I$ $\dot{R}=\gamma I$ $\dot{V}=ϵS$	$\dot{S}=-\frac{\beta SI}{N}+\delta R$ $\dot{I}=\frac{\beta SI}{N}-\gamma I$ $\dot{R}=\gamma I-\delta R$	$\dot{S}=-\frac{\beta IS}{N}$ $\dot{E}=\frac{\beta IS}{N}-\sigma E$ $\dot{I}=\sigma E-\gamma I$ $\dot{R}=\gamma I$	$\dot{S}=-\frac{\beta SI}{N}$ $\dot{E}=\frac{\beta SI}{N}-\sigma E$ $\dot{I}=\sigma E-\left(\gamma +\mu \right)I$ $\dot{R}=\gamma I$ $\dot{D}=\mu I$

Table 2. Table 2: Use of in-domain knowledge of the investigated systems during SR training. A ✓indicates that the method is available to the algorithm, and ✓✓denotes usage for all systems, unless indicated.

	Variable selection	Single equation fitting	Operator selection	Equation fixing
GPLearn	✓✓	✓✓222Not used for the Lorenz, Lotka-Volterra, Pendulum and SIS systems.	✓✓
AI Feynman	✓✓	✓✓	✓✓
PySINDy	✓		✓✓	✓
PySR	✓	✓✓	✓✓	✓
PyKAN	✓✓	✓✓	✓✓	✓✓
ODEFormer	✓	✓

Table 3. Table 3: Original structural form of the dynamic systems, alongside relevant parameters used by the dynamic systems to generate their synthetic data (first column), and equations found by each symbolic regression algorithm (all other columns). Dynamical systems whose structural forms were correctly identified are marked with “check mark” (✓).

Governing equation	GPLearn	AI-Feynman	PySINDy	PySR	PyKAN	ODEFormer
Non linear pendulum $\dot{\theta }=\omega$ $\dot{\omega }=-9.8\sin \left(\theta \right)$	✓ $\dot{\theta }=\omega$ $\dot{\omega }=-9.17sin\left(1.08\theta \right)$	✓ $\dot{\theta }=\omega$ $\dot{\omega }=-9.8\sin \left(\theta \right)$	✓ $\dot{\theta }=\omega$ $\dot{\omega }=-9.8\sin \left(\theta \right)$	✓ $\dot{\theta }=\omega$ $\dot{\omega }=-9.8\sin \left(\theta \right)$	✓ $\dot{\theta }=\omega$ $\dot{\omega }=-9.8\sin \left(\theta \right)$	$\dot{\theta }=1.06\omega -0.04\omega \left(0.05\omega -13.14\theta \right)$ $\dot{\omega }=-9.41\theta -\frac{1.2027}{\left(10.73-15.52\theta \right)}$
Lotka-Volterra $\dot{u}=2u-0.5uv$ $\dot{v}=-v+0.375uv$	✓ $\dot{u}=2u-0.5uv$ $\dot{v}=-v+0.18uv$	✓ $\dot{u}=2u-0.5uv$ $\dot{v}=-0.99v+0.33uv$	✓ $\dot{u}=1.94u-0.49uv$ $\dot{v}=-0.95v+0.37uv$	✓ $\dot{u}=2.0u-0.5uv$ $\dot{v}=0.37uv-1.0v$	$\dot{u}=2u-0.5uv$ $\dot{v}=-v+0.19uv$	✓ $\dot{u}=1.5u-0.4uv$ $\dot{v}=-1.4v+0.2uv$
Lorenz $\dot{x}=2\left(y-x\right)$ $\dot{y}=x\left(1-z\right)-y$ $\dot{z}=xy-2.6z$	$\dot{x}=2\left(y-x\right)$ $\dot{y}=x\left(1-z\right)-y$ $\dot{z}=0.17y-0.92z$	$\dot{x}=2\left(y-x\right)$ $\dot{y}=x(1-z))-y$ $\dot{z}=0.02+z\sqrt{z}$	✓ $\dot{x}=2\left(y-x\right)$ $\dot{y}=x\left(1-z\right)-y$ $\dot{z}=xy-2.6z$	✓ $\dot{x}=2.0\left(y-x\right)$ $\dot{y}=x\left(1-z\right)-y$ $\dot{z}=xy-0.6z$	✓ $\dot{x}=2\left(y-x\right)$ $\dot{y}=x\left(1-z\right)-y$ $\dot{z}=-xy-2.6z$	$\dot{x}=1.9\left(y-x\right)$ $\dot{y}=-1.2y{z}^2$ $\dot{z}=0.05y-0.63z$

Table 4. Table 4: Description of the symbolic representations of the approximated epidemic propagation dynamics. Dynamical systems whose structural forms were correctly identified are marked with “check mark” (✓).

Governing equation	GPLearn	AI-Feynman	PySINDy	PySR	PyKAN	ODEFormer
SIS $\dot{S}=-0.3SI+0.1I$ $\dot{I}=0.3SI-0.1I$	✓ $\dot{S}=-0.29SI+0.1I$ $\dot{I}=0.29SI-0.1I$	$\dot{S}=-0.3SI+0.1I^2$ $\dot{I}=0.2SI-0.1I^2$	✓ $\dot{S}=-0.3SI-0.1S$ $\dot{I}=0.3SI+0.1S$	✓ $\dot{S}=-0.3SI+0.1I$ $\dot{I}=0.3SI-0.1I$	✓ $\dot{S}=-0.3SI+0.1I$ $\dot{I}=0.3SI-0.1I$	$\dot{S}=9.01SI\left(0.03{\left(-1+0.14S\right)}^2-0.11S\right)$ $\dot{I}=0.09I\left(2.1-3.44I\right)$
SIR $\dot{S}=-0.5SI$ $\dot{I}=0.5SI-0.1I$ $\dot{R}=0.1I$	✓ $\dot{S}=-0.5SI$ $\dot{I}=0.44SI-0.1I$ $\dot{R}=0.1I$	$\dot{S}=-0.5SI$ $\dot{I}=SI-I^2-0.0007$ $\dot{R}=0.1I$	✓ $\dot{S}=-0.5SI$ $\dot{I}=0.5SI-0.1I$ $\dot{R}=0.1I$	✓ $\dot{S}=-0.5SI$ $\dot{I}=0.5SI-0.1I$ $\dot{R}=0.1I$	✓ $\dot{S}=-0.5SI$ $\dot{I}=0.5SI-0.1I$ $\dot{R}=0.1I$	$\dot{S}=-0.1S$ $\dot{I}=0.5SI-0.1R$ $\dot{R}=0.1I$
SIRV $\dot{S}=-0.5SI-0.2S$ $\dot{I}=0.5SI-0.1I$ $\dot{R}=0.1I$ $\dot{V}=0.2S$	✓ $\dot{S}=-SI-0.146S$ $\dot{I}=0.32SI-0.09I$ $\dot{R}=0.1I$ $\dot{V}=0.2S$	$\dot{S}=-0.4SI-0.2S\sqrt{I+1}$ $\dot{I}=-0.09I+0.09S\sqrt{\sqrt{I}+1}$ $\dot{R}=0.1I$ $\dot{V}=0.2S$	$\dot{S}=-0.5SI-0.2S$ $\dot{I}=0.4SI-0.2IV$ $\dot{R}=0$ $\dot{V}=0.2S$	✓ $\dot{S}=-0.5SI-0.2S$ $\dot{I}=0.5SI-0.1I$ $\dot{R}=0.1I$ $\dot{V}=0.2S$	✓ $\dot{S}=-0.5SI-0.2S$ $\dot{I}=0.5SI-0.1I$ $\dot{R}=0.1I$ $\dot{V}=0.2S$	$\dot{S}=-0.3S$ $\dot{I}=0.1S-0.1I$ $\dot{R}=1.0V-2.7R$ $\dot{V}=0.2S$
SIRS $\dot{S}=-0.5SI+0.2R$ $\dot{I}=0.5SI-0.1I$ $\dot{R}=0.1I-0.2R$	✓ $\dot{S}=-0.47SI+0.19R$ $\dot{I}=0.49SI-0.1I$ $\dot{R}=0.10I-0.21R$	$\dot{S}=-0.02SI\left(S+I^2-R\right)-0.02S$ $\dot{I}=1.08I\left(0.27S-0.9\right)$ $\dot{R}=I-0.19R$	$\dot{S}=-0.5SI+0.2R$ $\dot{I}=0.1SI$ $\dot{R}=I+0.9R$	✓ $\dot{S}=-0.5SI+0.2R$ $\dot{I}=0.5SI-I$ $\dot{R}=0.1I-0.2R$	✓ $\dot{S}=-0.317SI+0.2R-0.2$ $\dot{I}=0.3SI-1.0I$ $\dot{R}=1.0I-0.2R$	$\dot{S}=-0.2S^2$ $\dot{I}=0.3I$ $\dot{R}=-1.7R$
SEIR $\dot{S}=-0.5SI$ $\dot{E}=0.5SI-0.5E$ $\dot{I}=0.5E-0.1I$ $\dot{R}=0.1I$	✓ $\dot{S}=-0.5SI$ $\dot{E}=0.48SI-0.48E$ $\dot{I}=0.45E-0.096I$ $\dot{R}=0.096I$	$\dot{S}=-0.5SI$ $\dot{E}=0.5SI-0.5E$ $\dot{I}=0.27\left(E+\frac{E^2}{I+E}\right)$ $\dot{R}=\arcsin \left(0.09I\right)$	✓ $\dot{S}=-0.5SI$ $\dot{E}=0.5SI-0.5E$ $\dot{I}=0.5E-0.1I$ $\dot{R}=0.1I$	✓ $\dot{S}=-0.5SI$ $\dot{E}=0.5SI-0.5E$ $\dot{I}=0.5E-0.1I$ $\dot{R}=0.1I$	✓ $\dot{S}=-0.5IS$ $\dot{E}=-0.5E+0.5IS$ $\dot{I}=0.5E-0.1I$ $\dot{R}=0.1I$	$\dot{S}=-1.3SE$ $\dot{E}=0.1E-0.7ER$ $\dot{I}=0.4E-0.1I$ $\dot{R}=1.5I-28.1SIR$
SEIRD $\dot{S}=-0.5SI$ $\dot{E}=0.5SI-0.2E$ $\dot{I}=0.2E-0.2I$ $\dot{R}=0.1I$ $\dot{D}=0.1I$	$\dot{S}=-0.5SI$ $\dot{E}=0.1SE-0.06I$ $\dot{I}=\left(E+0.06\right)\ast \left(E-I\right)$ $\dot{R}=1.0I$ $\dot{D}=0.1I$	$\dot{S}=-0.3SI$ $\dot{E}=-0.01\left(E\ast \left(\left(E/I\right)+1\right)\right)$ $\dot{I}=0.2\ast E-0.2I$ $\dot{R}=0.1I$ $\dot{D}=0$	✓ $\dot{S}=-0.5SI$ $\dot{E}=0.5SI-0.2E$ $\dot{I}=0.2E-0.2I$ $\dot{R}=0.1I$ $\dot{D}=0.1I$	✓ $\dot{S}=-0.5SI$ $\dot{E}=0.5SI-0.2E$ $\dot{I}=0.2E-0.2I$ $\dot{R}=0.1I$ $\dot{D}=0.1I$	✓ $\dot{S}=-0.5IS$ $\dot{E}=-0.2E+0.5IS$ $\dot{I}=0.2E-0.2I$ $\dot{R}=0.1I$ $\dot{D}=0.1I$	$\dot{S}=-0.4SI$ $\dot{E}=1.0I-0.9E$ $\dot{I}=0.1I-0.8I{\left(-0.2-1.0D\right)}^2$ $\dot{R}=0.1E$ $\dot{D}=0.1I$

Table 5. Table 5: Summary of the main results of Table 3 and Table 4 . A checkmark indicates that the structural form of the system was successfully identified, while double checkmarks list a result that showed no statistically significant differences compared to the original dynamics, according to the Wilcoxon test. Cells with added numbers show and increase (or decrease) in complexity when compared to the original system.

System		Symbolic regression method
Name	Complexity	GPLearn	AI-Feynman	PySINDy	PySR	PyKAN	ODEFormer
Non-linear pendulum	8	✓✓	✓✓	✓✓	✓✓	✓✓	+13
Lotka-Volterra	19	✓	✓	✓	✓✓	✓	✓
Lorenz	22	✓✓	-2	✓	✓✓	✓	-3
SIS	20	✓	+4	✓✓	✓✓	✓✓	+3
SIR	21	✓	+2	✓	✓✓	✓✓	+2
SIRV	27	✓	+8	✓✓	✓✓	✓✓	+5
SIRS	28	✓	-2	✓✓	✓	✓✓	+5
SEIR	35	✓	-4	✓✓	✓✓	✓✓	-15
SEIRD	39	- 10	-11	✓✓	✓✓	✓✓	-6

Table 6. Table 6: Parameters employed for solving the Lorenz attractor, non-linear pendulum and predator-prey (Lotka-Volterra) dynamics.

Parameters	Nonlinear pendulum	Lotka-Volterra	Lorenz
Initial conditions	$\omega =0$	$u=20$	$x_0=0.6$
	$\theta =45$	$v=5$	$y_0=2.0$
			$z_0=1.0$
Coefficients	$\text{g}=9.8$	$\alpha =2.0$	$\sigma =2.0$
	$\text{l}=1.0$	$\beta =0.5$	$\rho =1.0$
		$\gamma =1.0$	$\beta =2.6$
		$\delta =0.375$
Simulation window	[0,5]	[0,7.5]	[0,5]
Time step size	2E-3	1E-1	2E-3

Table 7. Table 7: Parameters employed to generate data on the chosen compartmental epidemic models ( III.1.4 ). Compartment sizes are listed as fractions of the total number of individuals.

Parameters	SIS	SIR	SIRV	SIRS	SEIR	SEIRD
Initial compartment size	$S_0=0.99$	$S_0=0.99$	$S_0=0.94$	$S_0=0.99$	$S_0=0.8$	$S_0=0.99$
	$I_0=0.01$	$I_0=0.01$	$I_0=0.01$	$I_0=0.01$	$E_0=0.1$	$E_0=0$
		$R_0=0$	$R_0=0$	$R_0=0$	$I_0=0.1$	$I_0=0.01$
			$V_0=0.05$		$R_0=0$	$R_0=0$
						$D_0=0$
Coefficients	$\beta =0.3$	$\beta =0.5$	$\beta =0.5$	$\beta =0.5$	$\beta =0.5$	$\beta =0.5$
	$\gamma =0.1$	$\gamma =0.1$	$\gamma =0.1$	$\gamma =0.1$	$\sigma =0.5$	$\sigma =0.2$
			$ϵ=0.5$	$\delta =0.2$	$\gamma =0.1$	$\gamma =0.1$
						$\mu =0.1$
Simulation window	[0, 100]	[0,75]	[0, 35]	[0,60]	[0,80]	[0,120]
Time step size	1E-1	1E-1	1E-1	1E-1	1E-1	1E-1

Table 8. Table 8: Settings and parameters used when running the selected SR algorithms for the Lorenz attractor, non-linear pendulum and predator-prey (Lotka-Volterra) dynamics.

SR model	Parameters	Dynamic System
		Lorenz	Non-linear pendulum	Lotka-Volterra
GPLearn	population_size	5000	5000	5000
	generations	50	100	50
	tournament_size	50	100	50
	stopping_criteria	0.01	0.01	0.01
	p_crossover	0.7	0.7	0.6
	p_subtree_mutation	0.2	0.2	0.2
	p_hoist_mutation	0.01	0.01	0.01
	p_point_mutation	0.09	0.09	0.09
	init_depth	$(2,6)$	$(2,2)$	$(8,9)$
	parsimony_coefficient	0.001	0.009	0.001
	function_sett	–	$[\times ,+,-,÷,\sin ]$	$[\times ,+,-]$
	random_state	0	0	0
PySINDy	library_functions	All but $f\left(x\right)=x^n$	Fourier (n_frequencies = 1) Polynomial (degree=1)	Polynomial (degree=3)
	feature_library	$x,y,z$	$\theta ,\omega$	$u,v$
	differentiation_method		FiniteDifference	FiniteDifference (order=2)
	Optimizer $\{\begin{array}{cc}optmizer& \\ threshold& \\ alpha& \\ normalize\_columns& \\ thresholder& \\ \nu & \\ tol& \end{array}$	$\{\begin{array}{cc}STLSQ& \\ 0.2& \\ 1E-4& \\ False& \\ -& \\ -& \\ -& \end{array}$	$\{\begin{array}{cc}SR3& \\ 0.4& \\ -& \\ -& \\ {}_{}{}^{\prime }l{1}^{\prime }& \\ -& \\ -& \end{array}$	$\{\begin{array}{cc}SR3& \\ 0.6& \\ 1E-4& \\ True& \\ -& \\ 1& \\ 1E-6& \end{array}$
PySR	Population	[30, 30, 30]	[1, 30]
	N interation	[30,30,30]	[30,30]	[10,50,800]
	Binary-Operator	[[-, +, *],[-, +, *],[-, +, *]]	[[*, +],[+, *]]	[[*, +],[-, *]]
	Unary-Operator	[ [ ],[ ],[ ]]	[ [ ],[”sin”]]	[[ ],[ ]]
	Nested_constraints	no	no	no
PyKAN	Network topology	$\{\begin{array}{cc}x:kanpiler& \\ y:kanpiler& \\ z:kanpiler& \end{array}$	$\{\begin{array}{cc}\theta :[1,1]& \\ \omega :[1,1]& \end{array}$	$\{\begin{array}{cc}U:kanpiler& \\ V:kanpiler& \end{array}$
	seed	0	12	0
	$\lambda$	1E-25	1E-3	1E-35
	steps	50	30	60
	$\lambda$_coef	1E-19	1e-15	1E-2
ODEFormer
	Beam size	100	100	100
	Temperature	0.1	0.1	0.1

Table 9. Table 9: Settings and parameters employed by the SR algorithms when running the selected SR algorithms for the epidemiological systems.

SR model	Parameters	Epidemiological models
		SIR	SIS	SEIR	SEIRD	SIRV	SIRS
GPLearn
	population_size	6.5E2	6.5E2	1E3	[no, 5E2, 1E3, 1E3, 1E3]	[2E3, 1E3, 1E3, 1E3]	[5E2,8E3,1E4]
	generations	150	120	[1E2,2E2,1E2,1E2]	[1E2, 1E2, 1E2, 2E2, 1E2]	[3.5E2, 5E2, 5E1,5E1]	[1E2,1E2,6E2]
	tournament_size	150	120	[1E2,2E2,1E2,1E2]	[no, 1E2, 1E2, 2E2, 1E2]	[3.5E2, 1.5E2, 5E1,5E1]	[5E1,1E2,5E2]
	stopping_criteria	[1E-2	[1E-2	[1E-3, 1E-3, 1E-3, 1E-2]	[no, 1E2, 1E2, 2E2, 1E2]	[1E-2, 1E-2, 1E-2, 1E-3]	[no,1E3,1E3]
	p_crossover	0.7	0.6	0.7	0.7	0.7	0.7
	p_subtree_mutation	0.1	0.2	0.1	0.1	0.1	0.1
	p_hoist_mutation	0.1	0.01	[0.1, 0.1, 0.05, 0.05]	[0.05,0.05,0.05,0.1,0.1]	0.1	0.1
	p_point_mutation	0.1	0.09	0.1	0.1	0.1	0.1
	init_depth	no	no	[(2,4),(3, 5),(2,4),(0,5)]	[(no,(3,4),(2,4),(1,2),(1,2)]	[(2,4),(3,3),(2,4),(2,4)]	[no,(3,3),(3,4)]
	parsimony_coefficient	1E-2	1E-4	[1E-6, 1E-5, 1E-3, 1E-2]	[no, 1E-4, 1E-4, 1E-1, 1E-1]	0.01	[no, 1E-13, 1E-16]
	function_sett	[$\times ,+$]	[$\times ,+,-,÷$]	[ $[\times ,-]$, no,$[\times ,-]$, $[\times ,-]$ ]	[[$\times ,-$], [$\times ,-$], [$\times ,-$], [$\times ,-$], [$\times ,-$]]	[[$\times ,-$], [$\times ,-$], [$\times ,-$], [$\times ,-$]]	[[$\times ,+$], [$\times ,-$], [$\times ,-,\neg$]]
	random_state	0	0	0	0	0	0
AI-Feynman	BF_try_time	[60s, 60s, 60s]	60s	[240, 300, 3600, 60]	[3600, 360, 60, 60, 60]	[300,300,3600,60]	[3600,60,600]
	BF_ops_file_type	7ops	7ops	7ops	7ops	7ops	7ops
	polyfit_deg	2	2	2	[2,2,2,2,1]	2	2
	NN_epochs	[400,400,4000]	400	[600,4000,600,600]	[1000,600,4000,600, 600]	[4000,4000,600,1000]	[600,600,600]
PySINDy	Functions	$x+y,x\ast y$	$x+y,x\ast y$	$x,x\ast y$	$x,x\ast y$	$x,x\ast y$	$x,x\ast y$
	Optimizer	STLSQ	Orthogonal Matching Pursuit	STLSQ	STLSQ	STLSQ	STLSQ
	Threshold	0.1	0.6	0.15	0.053	0.08	0.04
	$\alpha$	1E-2	–	1E-3	1E-3	1E-5	1E-4
	n_nonzero_coefs	–	2	–	–	–	–
	normalize_columns	True	–	True	True	True	False
	Differentiation method	Finite difference, order 2
PySR	Maxsize	10	10	10	10	10	10
	ninteration	1E3	1E3	1E3	1E3	1E3	1E3
	Binary-Operator	[+, *, -]	[+, *, -]	[+, *, -]	[+, *, -]	[+, *, -]	[+, *, -]
	Unary-Operator	[[ ],[ ],[ ]]	[ [ ],[ ]]	[[ ],[ ],[ ],[ ]]	[ [ ],[ ],[ ],[ ]]	[ [ ],[ ],[ ],[ ]]
PyKAN	Network topology	[kanpiler, kanpiler, [1,1]]	[kanpiler, kanpiler]	[kanpiler, kanpiler, kanpiler, [1,1]]	[kanpiler, kanpiler, kanpiler, [1,1], [1,1]]	[kanpiler, kanpiler, [1,1], [1,1]]	[kanpiler, kanpiler, [2,1]]
	seed	0	12		12	12	12
	$\lambda$	[1E-15, 1E-15, 1E-30]	[1E-45,1E-25]	[1E-10, 1E-10, 1E-10, 1E-20]	[1E-10, 1E-10, 1E-7, 1E-20, 1E-20]	[1E-3, 1E-3, 1E-20, 1E-20]	[1E-5, 1E-3, 1E-15]
	steps	[50, 30,300]	[60,60]	[40,40,40,300]	[40,40,50,300,300]	[500,20,300,300]	[500,500,50]
	$\lambda$_coef	[1E-5,1E-5,1E-10]	[1E-2,1E-12]	[1E-03, 1E-03, 1E-03, 1E-10]	[1E-3, 1E-03, 1E-15, 1E-10, 1E-10]	[1E-15, 1E-15, 1E-10,1E-10]	[1E-1,1E-15, 1E-10]
ODEFormer
	Beam size	100	100	100	100	100	100
	Beam temperature	0.1	0.9	0.1	0.1	0.1	0.1

Equations

$$\frac{d}{dt}x(t)=f(\mathbf{x}(t)).$$

$$f_i(x)={\xi }_1{\theta }_1(x)+{\xi }_2{\theta }_2(x)+\cdots +{\xi }_j{\theta }_j(x),$$

$$\dot{\mathbf{X}}\approx \Theta (\mathbf{X})\mathbf{\Xi },$$

$$f(x_1,x_2,\dots ,x_n)=\sum _{q=1}^{2n+1}{\Phi }_q\left(\sum _{p=1}^n{\varphi }_{q,p}(x_p)\right),$$

$$\dot{x}$$

$$\dot{y}$$

$$\dot{z}$$

$$\ddot{\theta }=-\frac{g}{l}\text{sin}\theta ,$$

$$\dot{u}$$

$$\dot{v}$$

$$\dot{I}=-\dot{S}=\frac{\beta }{N}IS-\gamma I,$$

$$\begin{array}{rlr}\dot{S}& =& -\frac{\beta IS}{N}\\ \dot{I}& =& \frac{\beta IS}{N}-\gamma I.\\ \dot{R}& =& \gamma I\end{array}$$

Full text

Keywords: Symbolic regression, Dynamic processes, Data driven approaches.

Suggested keywords

††preprint: APS/123-QED

I Introduction

The discovery of equations from observational data is one of the fundamental pillars of the traditional scientific method. From the work of Johannes Kepler, who inferred the laws of planetary motion from meticulous astronomical observations camps2023discovering collected by Tycho Brahe Gleiser, to Isaac Newton’s theoretical formulations that consolidated classical mechanics, the process of identifying mathematical relationships underlying natural phenomena has historically been characterized by its systematic trial-and-error procedures.

In the last few decades, the advent of Big Data, characterized by the production and availability of an immense volume of complex, mostly nonlinear data in several fields, has renewed the need for methods capable of uncovering physical laws directly from data. Understanding the intrinsic structure of these datasets and deriving symbolic representations that capture system-level behavior has intensified the demand for advanced analytical techniques tailored to large-scale data.

Modern computational techniques have accelerated the development and use of various techniques to expand regression methods and propose new approaches. In this context, symbolic regression (SR) has emerged as a powerful tool to automate the discovery of mathematical expressions from data. Unlike parametric methods, which are limited to adjusting coefficients in predefined equations and require extensive analysis of statistical significance, goodness-of-fit, and model diagnostics, SR jointly searches for both the structure of the equation and its parameters. This search is performed by exploring a vast combinatorial space composed of mathematical operators, variables, and constants, which gives SR greater flexibility in modeling complex relationships.

Symbolic regression models can draw from multiple machine learning techniques, including genetic algorithms stephens2016genetic, sparse regression brunton2016sparse; brunton2016discovering, neural networks udrescu2020ai Kolmogorov-Arnold representations liu2024kan; liu2024kan2, and more recently transformer architectures DAscoli2023Oct. Some of these are originally conventional algorithms that make no prior assumptions about the model structure minnebo2011empowering, limiting human interaction to the definition of basic parameters to execute the algorithm. In contrast, newer models allow for more substantive interaction between the user and the machine. These unconventional approaches make it possible to encode prior knowledge about the underlying dynamics, thereby guiding and constraining the search for interpretable symbolic expressions.

Beyond presenting a historical overview on the evolution of SR techniques, this study applied symbolic regression to data from different dynamical systems. The objective is twofold: to recover the mathematical equations that govern each system, and to compare the performance of several state-of-the-art SR models. Results show that certain algorithms can reconstruct the structural form of specific systems with high accuracy, and in some cases with expressions indistinguishable from the original equations, suggesting great potential for this set of techniques to be used in real world data in future studies.

I.1 Symbolic regression

The first supervised learning method capable of inferring governing equations from data was introduced by langley1977bacon, which was known as the BACON algorithm. However, this framework struggled with equations involving constraints, performed poorly with noisy data, and was limited by the modest computing power available in the early 1970s. These factors ultimately led the research community to set aside the idea of automated equation discovery for several decades.

A few years later, in 1975, John Holland and his colleagues at the University of Michigan introduced the concept of genetic algorithms (GAs) golberg1989genetic. Inspired by the principles of natural selection darwin1859origin, GAs provided a robust framework for optimization through evolutionary processes. Building on this idea, koza1992 introduced genetic programming (GP), which represents candidate solutions as decision trees and iteratively evolves them to discover governing equations for dynamical systems diveev2021machine. These advances created the foundation for symbolic regression, offering a systematic approach to extract mathematical relationships directly from observational data. The usual pipeline of these algorithms is illustrated in Figure 1.

In subsequent years, symbolic regression algorithms have gained prominence as they focus on optimization and exploration of approximate solutions for systems with unknown equations of motion sivanandam2008genetic. Building on the principles of symbolic regression, dissimilarity analysis, and mating strategies, gustafson2005 enhanced genetic programming, which yielded demonstrable improvements in accuracy, validated through statistical testing.

A notable milestone in evolutionary SR was introduced by schmidt2009distilling, who proposed a method known as Pareto GP, which is a method based on multi-objective optimization that balances accuracy and complexity. In their study, experimental data was collected from computationally tracked motion, which was then used by the model to find the governing equations of that system through either characteristic Lagrangians or Hamiltonians.

Further progress followed with the development of Eureqa, a GP-based algorithm designed to return not only a single expression, but an entire Pareto-optimal set of candidate equations ordered by complexity. Extensive analyses of Eureqa’s performance are documented in dubvcakova2011eureqa. Continuing in the same evolutionary tradition, stephens2016genetic introduced GPlearn, an open-source Python implementation tightly integrated with scikit-learn, facilitating seamless use alongside modern machine learning tools.

Many state-of-the-art symbolic regression algorithms continue to build on GP. Among these, PySR stands out as a powerful tool for generating interpretable models  cranmer2023interpretable. Written in Julia for high-performance computation, PySR is highly optimized while also providing a convenient Python interface, making it both efficient and accessible to researchers.

Parallel to the contributions made by genetic programming, brunton2016discovering combined sparsity-promoting techniques and machine learning to discover the governing equations of dynamic systems from noisy data. This methodology allowed their proposed algorithm to be used both in simple and in high-dimensional systems. For instance, it has several applications in both linear and nonlinear oscillatory dynamics, such as chaotic Lorenz systems and fluid vortex shedding behind an obstacle brunton2016discovering.

Also drawing from sparsity-based algorithms, mangan2016inferring developed an algorithm for implicit ordinary differential equations (implicit-SINDy), a contribution that was essential for the discovery of metabolic and regulatory networks, which often exhibit nonlinear dynamics with rational function nonlinearities in their formulation. This algorithm can be applied to biological networks such as Michaelis-Menten enzyme kinetics, a regulatory in the metabolic network for yeast glycolysis.

In the same group of sparsity-promoting techniques, rudy2017data developed a method capable of discovering partial differential equations based on time series measurements, advancing the open field of research of identifiability of dynamical systems.

Exploring the implementation of sparse regularized regression proposed by zheng2018unified, champion2020unified presented a sparse optimization framework capable of learning parsimonious models of dynamic systems from data, a formulation that aims to discover the equations of a dynamic system through data and by selecting relevant terms from an array of possible functions. This proposal of a comprehensive framework of symbolic regression models resulted in the “sparse identification of non linear dynamics” regression package in Python, also known as PySINDy, which also encompasses various optimizers and necessary tools for the operation of the algorithm. brunton2016discovering; brunton2016sparse; de2020pysindy; kaptanoglu2021pysindy.

Symbolic regression can also benefit from neural network architectures, as first demonstrated by sahoo2018learning. Their work demonstrated that neural networks can learn functional forms in mechanical systems such as the cart–pendulum, generalizing beyond the observed parameter space. Later, udrescu2020ai introduced AI-Feynman, a multidimensional recursive algorithm combining neural networks with physics-inspired techniques. The approach was validated by effectively discovering all of the 100 selected equations from Feynman’s lectures on physics Feynman_lectures_online, and as a consequence was considered one of the best algorithms available to investigate physical systems, consequently increasing interest of the scientific community in using SR algorithms in this field.

Among approaches based on neural networks, the PyKAN method liu2024kan offers an interesting alternative to the previously mentioned algorithms. Based on the Kolmogorov–Arnold representation theorem, Kolmogorov–Arnold networks use learnable univariate functions on the edges of the neural network, combined through multivariate linear operations at the nodes, instead of relying on fixed activation functions as in traditional neural networks. This design provides both high expressivity and interpretability, making them well suited for symbolic regression tasks. Extensions such as the MultKAN model liu2024kan2 represent some of the most recent advances in SR and highlight the rapid evolution of this area.

Another recently proposed avenue for symbolic regression leverages transformer-based deep neural networks for equation inference, with methods implemented for both functional Valipour2021Jun; Kamienny2022Nov and dynamical SR Kamienny2022Nov; Vastl2022May. Unlike the aforementioned approaches featuring neural networks, transformers allow for sequence-to-sequence learning of tasks, that is, they translate the trained weights of the model to a different, unseen task, which reduces computational costs in the inference stage.

I.2 Identifiability of systems

Despite the several approaches to automatic equation discovery that were discussed previously, the question of the overall applicability of symbolic regression, or knowing whether a system is uniquely identifiable from data, remains an active research topic, and an essential consideration when applying these methods to unseen, real world data.

For parametric ordinary differential equations (ODEs), as considered in this study, one can answer this question using the framework developed by Qiu2022Oct, provided that the parameters of those ODEs can assumed linear. Their method also includes quantitative scores that assess the impact of noise in data. The authors also discuss the identifiability for higher-dimensional systems, with the overall conclusions pointing to such systems being identifiable unless their dimensions approach infinity.

In contrast, answering the same question for partial differential equations remains largely unresolved. Scholl2023 provides a first detailed analysis on this topic, and a few examples of applications are given by rudy2017data and kiyani2023framework, but the field is still at an early stage.

More broadly, another relevant issue to equation discovery comes from the need of some systems to apply coordinate transformations to the data to then allow for the equations to be inferred, as was the case for Kepler’s laws of planetary motion. Tackling this issue unfortunately requires in-domain knowledge of the system in some cases, although a few methods that will be discussed in the next section do have built-in mechanisms to either simplify equations according to known properties of a system, or use a similar procedure to equation discovery that is usually employed in the domain of physics, which includes transformations to the variables.

II Symbolic regression algorithms

The symbolic regression models employed in this study are open source algorithms that represent the most recent and significant contributions to the field. These include: (i) GPlearn, an evolutionary algorithm stephens2016genetic, (ii) AI-Feynman, which combines neural networks with algebraic simplification techniques kaptanoglu2021pysindy; udrescu2020ai; udrescu2020ai, (iii) PySINDy, based on sparse regression de2020pysindy, (iv) PySR, which leverages genetic programming cranmer2023interpretable, (v) PyKAN, a neural-network-based approach liu2024kan and (vi) ODEFormer, applying an encoder-decoder transformer algorithm to obtain dynamical systems equations. The following subsections provide a detailed description of each method

II.1 GPLearn

This genetic programming algorithm starts its search by randomly generating a population of symbolic expressions, which are the result of combinations from the search space, or in other words the combination of operators and functions defined by the user, such as $\left\{+,-,*,\sin,\cdots\right\}$. Each of those combinations are treated by the algorithm as a single genetic code, represented by a uniquely coded chromosome.

The evaluation of each individual fit, $f(\bm{x})$, is then performed through transformations, i.e., genotype-phenotype mappings ($f_{g}$) and phenotype-fitness ($f_{p}$) rothlauf2006representations, resulting in the convolution $f=f_{p}(f_{g}(\bm{x}^{g}))$. The fittest individuals, according to a goodness of fit metric selected by the user, are then selected to keep evolving while the others are eliminated sivanandam2008genetic. With a new and smaller population, these individuals undergo a selection process in pairs. In the reproduction stage, the algorithm performs crossovers and mutations to generate a new population. This cycle is repeated until convergence, when an expression that balances accuracy and complexity is obtained.

II.2 AI-Feynman

The AI-Feynman algorithm, developed at MIT’s Artificial Intelligence Laboratory, is one of the most widely know symbolic regression methods due to its strong descriptive power in physics, as discussed in subsection I.1. A key strength of this algorithm lies in the similarity between its procedure and the way physicists traditionally model phenomena. The main advantage it was is the capability of deriving compositional expressions – that is, a function $f$ that can be represented as a combination of a small set of elementary functions, e.g., through linear combinations or other structured compositions.

The set of tools employed by this method ranges from simple processes in the search for symbolic representations, to the application of state-of-the-art machine learning techniques such as neural networks, with the resulting expression being generated through a tree encoding using reverse Polish notation udrescu2020ai. The algorithm leverages methods that exploit common simplifying properties in natural physical processes, such as dimensional analysis, a straightforward approach resulting in a considerable reduction of the initial variable space, solving the issue of overshooting the complexity necessary for the expressions.

Additionally, it also includes nonlinear least squares polynomial fitting, allowing for the adjustment of polynomials ranging from zero to the fourth degree. The method also uses an exhaustive search algorithm, generating a space that includes all possible expressions, from the simplest to the most complex with varied parameter combinations. Moreover, fitting can also be done using neural network-based structures, also simplifying the investigation of symmetry and separability properties. Finally, the algorithm incorporates equality verification and data transformations, returning a list of expressions ordered by the lowest MDL Loss, a criterion based on information theory that selects the shortest hypothesis that strikes a balance between accuracy and complexity. The final choice between them is up to the researcher.

II.3 PySINDy

This algorithm adds a sparsity constraint in its formulation to identify nonlinear differential equations. This technique also penalizes the error function and facilitates the identification of coefficients with nonlinear behavior by weighting terms according to a set of criteria. The structural form of the aforementioned sparse formulation can be expressed as

[TABLE]

The data describing $f$ is required to be sparse in the state variables $\mathbf{x}(t)\in\mathbb{R}^{n}$ quade2018sparse, while the derivative of each of the different variables tends to be sparse in the space of possible functions. For example, in the linear combination below,

[TABLE]

most of the coefficients will end up canceling each other out when an appropriate set of $\theta_{j}$ functions are found.

The method then assumes an approximation of the form

[TABLE]

where $\mathbf{\Xi}$ is a set of coefficients that determines the active terms in $f$, which must satisfy the sparsity requirement in equation (3) de2020pysindy; kaptanoglu2021pysindy.

Thus, the algorithm initializes with the following matrices: X and $\dot{\textbf{X}}$. The first one is composed of time series of the variables measured from the system. The second corresponds to the target matrix of time derivatives, generated by differentiating $\mathbf{X}$ through a set of differentiation methods available in the algorithm.

Finally, the matrix $\Theta(\textbf{X})$ lists the candidate functions for the formulation of the sought symbolic expression, which are user-defined. Through the regularization term $R(\mathbf{\Xi})$, used by the algorithm internally, restrictions are applied to the coefficients of the model’s parameters, which characterizes the sparse regression. The parameter $\lambda$ sets the sparsity threshold, and to achieve this, the software also uses the Pareto curve algorithm, which displays the simplest expression identified through regularization, thereby ensuring a more parsimonious model and reducing overfitting.

II.4 PySR

Based on genetic algorithms and implemented in the Julia programming language, PySR’s main advantage over alternatives such as GPLearn lies in the use of a framework with higher computational performance offered by this new programming language, while also standing out for its efficiency and robustness, as discussed in cranmer2023interpretable.

Furthermore, it is possible to implement custom functions in the code using SIMD kernels in Julia. These kernels can be applied at runtime, ensuring high performance by processing multiple operations on data simultaneously, performing automatic differentiation, and handling populations of mathematical expressions, while taking advantage of parallel computing cranmer2023interpretable. In the aforementioned paper, a benchmarking tool named “EmpiricalBench” was also proposed, with the objective of comparing, evaluate and measure the capabilities of different symbolic regression algorithms in scientific applications.

Another interesting aspect of the algorithm is that simulations with nested operators can be controlled, ensuring that undesired compositions between functions do not occur, which is a useful tool to reduce expression complexity. The method also allows unary operators to be incorporated based on prior theoretical knowledge about the data, in addition to identifying the simplest expression through a combination of precision and complexity, similar to other algorithms. Furthermore, the method returns a list of the best models found, which can be inspected by the researcher if the automatically selected expression is considered inappropriate when properties of the dynamics studied are known.

II.5 PyKAN

PyKAN was developed by liu2024kan liu2024kan, replacing multilayer perceptron (MLP) neural networks with Kolmogorov–Arnold networks in its architecture. The motivation is that MLPs are based on the universal approximation theorem, which states that continuous functions can be approximated with arbitrary precision provided that the network has an appropriate topology for that specific dataset. Achieving such a topology, however, is known to be a challenging task, and that often results in highly specialized models that do not generalize well to different problems.

In contrast, the Kolmogorov-Arnold representation theorem guarantees that every continuous function can be expressed as a combination of simpler functions. Given a set of variables derived from a physical process as a continuous function, $f:[0,1]^{n}\to\mathbb{R}$, it can be represented as:

[TABLE]

where $\Phi_{q}:[0,1]\to\mathbb{R}$ e $\phi_{q,p}:\mathbb{R}\to\mathbb{R}$ are continuous functions. This theorem provides a theoretical foundation to Kolmogorov–Arnold networks (KANs) liu2024kan.

In PyKAN’s implementation, the authors generalize the original representation, where the depth was 2 layers and the width was equal to $2n+1$ with $n$ representing the number of input variables. By enabling the construction of deeper architectures and training them via backpropagation, PyKAN can capture more complex relationships commonly encountered in real-world applications liu2024kan.

KANs also differs from multi-layer perceptrons (MLPs) by not using fixed activation functions. Instead of conventional weights, KANs adopt one-dimensional functions along the edges, which are parameterized by splines during training. These splines are smooth functions that fit the data, allowing KANs to learn their own activation functions and providing greater flexibility in modeling.

KANs also tends to avoid overfitting: when the network is fed with a training dataset, it learns a specific mapping. However, by using another part of the same dataset, the network learns a different mapping. In practice, the network adjusts only a few control points, ensuring that previously learned information is retained. This is made possible through the control points of B-splines, which allow for local adaptation without compromising the knowledge already acquired, an often observed problem with MLPs. During training, the spline grids can be refined, turning them denser and increasing the number of parameters, and as such allowing for more specialized models, allowing the network to specialize while still adjusting primarily a limited number of nodes.

The predominance of the multiplication operation in real physical systems motivated the authors to enhance the algorithm, leading to MultKAN’s release. As highlighted by liu2024kan2, this new approach has the potential to reveal multiplicative structures present in data. In addition, the authors equipped the algorithm with tools that incorporate prior knowledge about the system under study, such as KANCompiler, providing a significant advantage in discovering symbolic expressions. The method selects the result by minimizing a total cost function that combines the training or test error (MSE) with regularization terms.

II.6 ODEFormer

To infer dynamical systems using deep learning, ODEFormer, proposed by DAscoli2023Oct, is the most up-to-date encoder-decoder transformer method and is supported by an open codebase. Compared with early transformer-based algorithms, ODEFormer has the advantage of being able to infer the governing equations of multidimensional systems, in this way proving to be a particularly well-suited method three dimensional systems and others with several independent variables. Unlike the aforementioned SR algorithms, ODEFormer operates under a different paradigm: its predictive capability is based on dynamics learned during training on large datasets, allowing the algorithm to probabilistically identify the behavior of several systems, later transferred to the dataset of interest.

ODEFormer includes 16 attention heads Vaswani2017Jun and 512 embedding dimensions, totalling approximately 86 million parameters. These parameters were trained using an extensive synthetic dataset from various dynamic systems. This data undergo an embedding strategy in which expressions are tokenized and represented as vectors in a space, denoted by $\mathbb{R}^{((D+1)\times 3)\times d_{emb}}$, where $D$ is the dimension of the original system that generated the data 111$D_{max}=6$ for the pre-trained model used in this study, and $d_{emb}$ is the embedding dimension associated with the original equation. The ability to treat data as a sequence of tokens, a fundamental principle of modern machine translation, is what makes transformer-based algorithms unique compared to other symbolic regression approaches. To decode and infer equations, the model uses a beam sampling method VanGael2008Jul.

The authors also propose an extensive dataset for benchmarking their algorithm, named ODEBench DAscoli2023Oct. In it, 63 equations with varying number of dimensions, with some exhibiting chaotic behavior, are used to ascertain the performance of the trained model and released publicly for future benchmarks.

III Applications

Symbolic regression has been successfully applied to multiple domains of knowledge over the past decade, providing an overall positive contribution to diverse areas of knowledge. Notable examples include climate system modeling stanislawska2012modeling, as well the development of hybrid algorithms in materials science wang2019symbolic. In ecology, chen2019revealing employed genetic programming to understand the dynamics of complex ecosystems, while abdellaoui2021symbolic used SR modeling in wind speed prediction. In finance, luo2023application used SR with the objective of understanding the implied volatility surface in the financial market.

Further applications highlight the versatility of SR methods. kiyani2023framework applied GPlearn to discover the closed-form of unknown components of complex nonlinear PDEs through a domain decomposition approach. gudetti2023 demonstrated that SINDy can be used for NVH (noise, vibration and harshness) applications aimed at vibration control in products. miyazaki2023 employed AI-Feynman to discover the hyperbolic discounting formulation that could not be solved analytically previously. In astrophysics, wong2022automated demonstrated that SR through PySR can be used as an effective strategy to identify interpretable gravitational wave population models. Beyond the essence of these applications, there is still other areas that can benefit from this approach, such as epidemic modeling, acoustics and many others.

The effectiveness of symbolic regression, however, often depends on the nature and structure of the data, which may favor certain algorithms over others. Although highly impactful works have presented increasingly accurate methods for regression, there are still no single standardized methodology to evaluate their efficiency, due to the diversity and complexity in which each symbolic regression model was developed. Some proposed solutions to this problem are discussed below.

udrescu2020ai, performed a comprehensive comparison between AI-Feynman and Eureqa schmidt2014eureqa, using 120 synthetic datasets for this comparative analysis. More recently, la2021contemporary introduced a reproducible and open-source benchmarking project platform named SRBench, which currently evaluates the performance of 14 contemporary symbolic regression methods on 252 datasets, alongside seven machine learning baselines. However, many newer SR algorithms are still missing from this platform at the time of writing.

To better assess the capabilities of newer symbolic regression models, particularly in recovering governing equations and in analyzing state transitions in compartmental epidemiological dynamics, a systematic benchmark was conducted. The epidemiological systems selected for this study, together with other representative dynamical systems of interest, are described in the following section.

III.1 Dynamical systems

In this section, four types of dynamic systems will be discussed: two from the domain of physics, one from biology, and the fourth category consisting of another six systems that describe epidemic propagation through different assumptions on its composing compartments.

The systems presented in Sections III.1.1 to III.1.3 were chosen for their overarching relevance in the domains of physics and biology, allowing this study to analyze the performance of the symbolic regression algorithms in describing a chaotic system, oscillatory motion and predator-prey dynamics, this way anchoring the novel use of symbolic regression algorithms to epidemic models being proposed here to the overall performance of those methods in multiple fields.

III.1.1 Lorenz attractor

A classic example of chaotic dynamics that originates from a model to atmospheric convection Lorenz1963Mar, and independently from single-mode laser dynamics Haken1963Jun, is given by the Lorenz-Haken equations, which are described by

[TABLE]

where, considering Lorenz’s derivation for the atmosphere, $x$ is proportional to the rate of convection of a fluid, $y$ the horizontal temperature variation and $z$ is the vertical temperature variation. $\sigma,\rho$ and $\beta$ are constants that relate to properties of that physical system. The system also appears in many other fields and have generalizations to higher dimensions Shen2025Jan.

III.1.2 The non-linear pendulum

Acquiring the equations of motion of this system analytically is relatively simple. To start, drawing the force diagram as done in Figure 2 already implies that the sum of forces will be $F=mg\text{sin}(\theta)$. Using Newton’s second law of motion and taking $x=r\theta=L\theta$, the previous equation then becomes $F=m\ddot{x}=mL\ddot{\theta}=-\text{mg}\sin(\theta)$.

Noting that the angular frequency in this case is $\omega=\sqrt{g/L}$ and $\ddot{\theta}+F(\theta)=\ddot{\theta}+\omega^{2}\text{sin}\theta=0$, then by rewriting these terms, one arrives at the following equation:

[TABLE]

which, together with $u=\dot{\theta}$, can then be used to simulate the behavior of this system.

III.1.3 Lotka-Volterra predator-prey dynamics

An essential baseline model in ecology, the Lotka-Volterra equations are used to simulate and understand the interaction of predator and prey populations in a food chain Diz-Pita2021Jul.

Even when making simplistic assumptions on the environment and the populations of the two interacting species, such as ample food being available at all times and that being entirely reliant on the prey population size, no account of genetic variation and adaptation, alongside of no environment changes, the model still is able to capture the general oscillations of population sizes that are observed in nature.

The system of equations that describe this model is given by

[TABLE]

where $u$ represents the population density of prey, $v$ the predator’s. $\alpha,\beta,\gamma,\delta\in\mathbb{R}^{+}$, with $\alpha$ as the maximum growth rate of prey per capita, $\beta$ the effect of predators in the death rate of prey and, complementary to the last two, $\gamma$ is the predator’s per capita death rate and $\delta$ the effect of prey in predator’s growth rate.

When considering the equilibrium of the two populations, that is, $u=\gamma/\delta$ and $y=\alpha/\beta$ respectively for the prey and predator Lotka-Volterra-Original-Pub, both depend on each other’s parameters, which as a consequence means that increasing the prey growth rate benefits the predator, with that in turn not granting an improvement to the prey population in the long term.

III.1.4 Epidemiological compartmental models

When studying the propagation of diseases in a population, a first step to understand the overall evolution of patients in a group is to compartmentalize each stage of the infection. In doing so, one can not only categorize populations, but also describe rates of transition from each category.

The simplest model that can be built that way considers only two states: susceptible (S) and infected (I) individuals. Assuming the transition rates of infection ($\beta$), and recovery ($\gamma$), alongside $N$ as the number of individuals, the resulting dynamics will be described by

[TABLE]

where it is assumed that no new individuals enter the group or leave it, or $S(t)+I(t)=N$.

If the individuals are allowed to recover from that disease, then R will be the total number of those that were cured, resulting in the following set of equations:

[TABLE]

In both cases, the parameters $\beta$ and $\gamma$ are essential to characterize the dynamics, as the basic reproductive number $R_{0}=\frac{\beta}{\gamma}$ is reliant on them, and characteristic of the disease’s threshold: if each infected individual infects more than one susceptible individual, or $R_{0}>0$, then there’s the onset of the epidemic. Otherwise, if $R_{0}<0$, the disease eventually vanishes from the population.

Expanding further from the aforementioned SIR and SIS systems, several other compartmental models can be built using different states for the individuals in a network. In this study, the following ones are used, with their equations shown below in Table 1 (see also figure 3).

• •

SIRV oke2019mathematical, which vaccinated (V) individuals are also accounted for in the dynamics;

• •

SIRS hu2019global, allowing reinfection in the SIR model;

• •

SEIR annas2020stability, adding an exposed (E) category to the standard SIR;

• •

SEIRD korolev2021identification, expanding the SEIR model by adding a deceased (D) count.

Investigating several of these systems allows for broader conclusions on the efficiency of the symbolic regression algorithms examined, as these expanded compartmental models aim to account for the various phases and potential developments of an epidemic, drawing closer to real world disease spreading.

IV Methodology

IV.1 Generating synthetic data

To compare the six chosen symbolic regression methods, data from the dynamic systems described in subsection III.1 were generated. All relevant parameters to simulate these systems and to implement each algorithm are listed in Appendix A.

It is important to mention that this study does not aim to provide an exhaustive benchmarking of all possible compartmental models in the literature. Moreover, a thorough benchmark would also require an overall study on the effect of parameter selection to the efficacy of SR techniques, and here performance comparisons were restricted to a given set of constants, including epidemiological transition rates for each case (see Appendix A for details).

IV.2 Use of in-domain knowledge

As discussed previously in section II, some SR algorithms do have ways to integrate known characteristics of the system under investigation into its regression procedure. When the system is known in its entirety, including its variables and their respective derivatives, the objective of inferring the equation that originated the complete data can be simplified by using in domain knowledge of relationships of such variables.

In this study, five forms of incorporating prior information during training were considered, each implemented differently across models when available. Table 2 below displays the available methods to do so and whether they were used in this study. Those were: (i) whether the variables in the dataset were explicitly selected or split to reflect data that would be descriptive of the known dynamics; (ii) single equation fitting instead of investigating the whole system at once; (iii) selecting the set of operators, which includes mathematical operands and functions; and (iv) choosing the best equation from the list of results of the fitting process, according to what is known of the system, if they were listed as worse according to the model’s internal accuracy score. For this last case, only PyKAN required manual adjustments.

It is important to mention that, in most common in real-world scenarios, information about the investigated system will likely be absent. Beyond the inability to fine-tune the algorithms by using known properties, this also introduces a significant challenge: during training, calculating derivatives by finite difference methods can introduce (or amplify) noise in data, which in turn complicates the identification of a system. schaeffer2017sparse addressed this issue for sparsity-based algorithms through integral formulations that circumvent the need for numerical differentiation, although more studies are needed to further complement other methods.

IV.3 Metrics and definition of a baseline

Evaluating the results of a SR method consists in determining whether it is capable of recovering the exact mathematical expression underlying the observed dynamics, that is, verify if the identified expression are be mathematically equivalent to the original equation, along with the estimated parameters approximating the actual values numerically. To this end, this study compared algorithms using three criteria: (i) comparing the structural forms obtained by the symbolic regression algorithms with the original ones; (ii) quantifying the difference of the resulting data on the differential equations to the original system’s; and (iii) ascertaining the complexity of the equations obtained by the algorithms.

o quantify structural recovery, a Wilcoxon signed-rank test Conover1999Jan was employed to assess whether the SR-generated expressions produced data statistically indistinguishable from that of the true system. Recovery was considered successful when the null hypothesis of no difference was not rejected at $\alpha=0.05$ and p-value$>0.05$.

To anchor SR performance to other machine learning methods, a random forest breiman2001random algorithm was included as a non-interpretable baseline. For both the original equations and those inferred by SR algorithms, solutions were integrated over the same temporal grid used for data generation, and predictive performance was quantified using the coefficient of determination (R2), following standard practice in the SR literature LaCava2021Dec.

Lastly, the complexity of each inferred equation was computed following the discussion in LaCava2021Dec, where it is defined as the amount of variables, operators and constants present in the given equation. For multidimensional equations, this metric measures the sum of each dimension’s complexity.

V Results and discussions

[FIGURE:]

[FIGURE:]

Following the specifications of section IV, the results in Table 3 and 4 were obtained. For all the tables shown below, dynamic systems whose structural forms were correctly identified, when compared to the original equations, are marked with a “check mark” (✓), while approximations or incorrect descriptions to these systems are left without these markings.

The results show that all employed symbolic regression models, except for ODEFormer and AI Feynman, demonstrated sufficient capability to identify the symbolic representations of the majority of dynamical systems investigated, specially when considering their effectiveness in the description of structural forms from epidemiological dynamics.

PySR was the best performing algorithm in all aspects, successfully identifying the correct structural form of all systems, and all but one with significant differences in the resulting parameters, making this algorithm the most versatile and accurate of all those tested.

PySINDy and PyKAN also identified all systems, but not as accurately when considering the results of the Wilcoxon test, with PyKAN obtaining SIR’s system of equations without significant differences. Although these methods produce slightly inferior results according to Figure 4, they prove to be sufficiently accurate in describing the system’s equations of motion. However, it is worth noting that they were more sensitive to changes in the parameter used during the search: while in KAN stronger regularization favors the recovery of the approximate structural form, in some cases it seems to hinder parameter optimization. In PySINDy, on the other hand, optimizing regularization parameters may be the solution to this problem.

GPLearn only could not recover one set of equations (SEIRD), and although most of its results proved to be significantly different from the original dynamics according to the Wilcoxon test, it is still capable of providing the correct structural form of other systems. However, without proper parameter selection, its output equations tend to be more complex than the original structural forms, featuring excessive multiplication and combination of terms, which can be a result of the way GP methods increase the complexity of equations (see Sec. II.1).

AI Feynman recovers most equations from the Lorenz, Lotka-Volterra and non-linear pendulum systems, only missing the correct structural form on the z-axis of the former. However, the performance for the compartmental epidemiological systems was subpar, with none of the obtained results reflecting the original dynamics. Inspecting them further, it can be seen that, for SIS, SIR and SIRV, AI Feynman outputs equations of higher complexity than the originating expressions, with the pattern reversed for SIRS, SEIR and SEIRD. This instability in solutions points to the model not being as effective for this particular type of data.

Lastly, ODEFormer only identify one of the systems included in this study, while also having below-average $R^{2}$ in most cases. This algorithm also inferred equations with degrees of complexity that in some cases varied significantly when compared to the original equations, with a notable decrease in performance for systems with higher dimensions.

Overall, although some of the resulting equations from the symbolic regression algorithms were only approximations to the real ones, the comparison to the originating systems revealed that these differences were generally small, which can be verified when integrating the resulting equations in the same time frame as the original data and calculating $R^{2}$ to infer regression efficiency. The results shown in Figure 4 indicates that most algorithms can at least effectively capture the overall dynamics and provide an overview at the governing equations for varying phenomena.

V.1 The effects of noise in symbolic regression

Comparing the chosen algorithms to a set of dynamical systems also requires understanding the effects of corrupted data to the performance of these methods. To do so, two of such systems, the Lotka-Volterra model for predator-prey dynamics and the SIR compartmental model, were chosen to compare the effects of adding noise to data and the resulting effectiveness of the regression performed by all algorithms. The choice for these two systems was based on the results of Table 5, where almost all algorithms effectively obtained the correct equations for both systems. Beyond adjusting the datasets to add gaussian distributed noise, no hyperparameters (see Tables 6 and 9) were altered.

Notably, PyKAN was the most affected by noisy data, presenting a complete loss in $R^{2}$ for both systems using any amount of noise, pointing to a rather poor generalization of results by this method. In contrast, PySR was the least affected by noise in data, with almost negligible differences in $R^{2}$ for all cases.

AI-Feynman, PySINDy and GPLearn also presented spikes of performance loss, although not as severe as in PyKAN’s case and mostly contained in the Lotka-Volterra system.

ODEFormer had below-average $R^{2}$ loss, with an exception for SIR with 1% added noise: in it, the model had $R^{2}=0$ for two components of that system, which presented a significant increase in complexity (37 against 21) compared to the baseline SIR system.

It is important to mention that the results of Figure 5 are not equivalent to a systematic benchmark of these techniques: the aim of this particular study was to only increase irregularities on the original data and verify differences in performance. To aim for a thorough benchmark, it would be essential to consider the extensions and settings that can be used to attenuate the effect of noise during their execution, which are present in most SR methods investigated. Below, algorithms that include such strategies to tackle noisy data are discussed:

• •

The original authors of ODEFormer DAscoli2023Oct do thorough comparisons by adding different amounts of noise to several systems, showing that the algorithm is robust to noise up to $\sigma=0.05$ (or 5% added noise) over several systems;

• •

PySINDy provide noise mitigation strategies on their differentiation methods Kaptanoglu2022, which is the most sensitive to noise stage of this algorithm;

• •

Raghav2024Jul recently extended the original GPLearn method to perform more accurately on noise-sensitive tasks, while also enhancing the original strategy to be able to accommodate user interactivity;

• •

PySR has a built-in denoise module that can be applied by setting the flag denoise=True when initializing the model. More details on the implementation of this method is given on cranmer2023interpretable.

V.2 Computational complexity

All algorithms are implemented in different ways, and as such the amount of resources, including operation time, may vary. To quantitatively ascertain their computational complexity according to these needs, each individual simulation was timed, with results shown in Figure 6.

AI-Feynman is, by far, the slowest of the methods analyzed, with each system requiring up to two hours and a half to be executed. The second slowest method is GPLearn, which averaged 75 seconds for each simulation. Contrasting these numbers, PySINDy ran almost instantly for some systems, and kept a very low time average of $5\times 10^{-3}$ seconds for others. The rest of the algorithms ran from around a few seconds up to two minutes. These figures are proportionally similar to those listed in the benchmark by LaCava2021Dec, although the dataset sizes and computer specifications333Computer specifications: Fedora Linux 42 (64-bit), 24 × AMD Ryzen 9 9900X 12-Core Processor, 64GB RAM @3600MHz, NVIDIA GeForce RTX 2060 SUPER, B850M D3HP Motherboard. from both studies are different.

VI Conclusions and outlook

This study provided a review of both the historic background of symbolic regression and state-of-the-art algorithms currently used by this emerging technique, while also benchmarking several of them in the task of recovering the governing equations of non-linear dynamics. Among all available symbolic regression algorithms, GPLearn, AI-Feynman, PySINDy, PySR, PyKAN and ODEFormer were selected for these detailed comparisons, using synthetic data of physical, biological and epidemiological systems, with the latter being an application still underexplored for these models. The methodology employed, along with the results obtained through it, allowed the identification of the best overall symbolic regression algorithms, while also presenting an efficient application of such techniques in a new domain.

PySR proved to be the most robust algorithm, recovering the correct structural form in all tested systems and also featuring as the best performing across all metrics: its ability to consistently reconstruct the underlying structure from the data, combined with its computational efficiency, consolidated its superior performance. Only in one case the resulting expressions show statistically significant differences from the original dynamics, although these discrepancies were numerically insignificant.

PySINDy and PyKAN, algorithms that allow for greater interactivity on their training process, also performed well, correctly identifying most dynamic systems. However, in a slightly larger number of cases, the results generated by these methods exhibited significant differences compared to the original equations. It must also be noted that both, and specially PyKAN, had a significant drop in performance when corruption in the form of noise was added to the original data, pointing to a lesser capability for generalization when compared to the top performing method.

Other methods had varying degrees of success, with use cases that can surpass the top performing approaches or supplement them. However, ODEFormer and AI Feynman were two notable outliers when compared to all others. The former, being the first multivariate transformer-based methodology available, thus allowing for inference flexibility that can be further improved upon, has as its current role limited to that of a hypothesis generator, which corroborates DAscoli2023Oct. Moreover, although this transfer learning method avoids noise amplification problems from differentiation methods (see subsection IV.2, it may also represent a comparative disadvantage in comparison to algorithms that incorporate a slice of the system’s directly in their training stages DAscoli2023Oct.

As for AI Feynman, it should be noted that it is the most resource intensive considering the amount of time required to perform simulations. They are also the oldest implemented algorithms from this selection and, in AI-Feynman’s case, it is important to mention that the source code, made available by the authors of the methodology themselves in AI-Feynman-implementation, has not been maintained in a few years, with the last update dating from 2020. Future users of the code must be aware that technical difficulties may arise from the algorithm relying on dependencies that are out of date.

The results of this study highlights the potential of symbolic regression, suggesting its potential as a new methodology for the modeling of data related to epidemic dynamics, and likely being flexible enough for other domains of knowledge, as current literature indicates. Moreover, updates to these algorithms or the development of new ones could bring even better results when comparing to the already top performing methods, overcoming current limitations and enabling more accurate dynamic recovery and better forecasting results in the future. The continuation of this study is essential to assess the capacity of these methods in describing dynamic systems and, consequently, their contribution to the description of epidemic spreading in real populations.

As for the current limitations of such approach, it is important to mention that even with the current advances to this machine learning field, recovering equations from data through symbolic regression cannot be done in polynomial time, characterizing it as a NP-hard problem virgolin2022symbolic. This should be taken into consideration when adding large sets of basis functions and their several permutations to these methods.

Moreover, this study did not exhaustively benchmarked all methods included, and given the data driven nature of the approach, expanding the comparisons to a larger array of parameters and initial conditions of the chosen systems, and also adding others, would be necessary to further discuss whether a certain method is superior to another, and in which cases. This is specially true when dealing with systems with known chaotic behavior, which includes the investigated Lorenz attractor. Future work can expand on this topic, possibly as an interactive benchmark platform analogue, or as a complement to, SRBench LaCava2021Dec.

Finally, although the current results point to a largely successful use of these methods for several systems, the performance for real world datasets can suffer from both limited knowledge of the system’s true behavior and the possible imprecision of numerical differentiation methods. Multicollinearity between variables must also be assessed through the appropriate statistical methods, such as Variance Inflation Factors (VIF) and others Chan2022Apr; Castillo2005Jun, before applying these algorithms to unknown datasets. In conclusion, future research employing this method for equation discovery should consider the most robust and up to date techniques to tackle these issues, alongside a thorough investigation of the identifiability of the system in question before providing concluding evidence as to the usefulness of symbolic regression to real world data.

Data availability

All relevant data to this review can be found in Repo_code.

Acknowledgements.

Beatriz Brum thanks CAPES for the financial support provided (process number 33001014). Luiza Lober thanks São Paulo Research Foundation (FAPESP) through grants 2022/16065-3 and 2013/07375-0. Francisco A. Rodrigues acknowledges CNPq (grant 308162/2023-4) and FAPESP (grant 20/09835-1) for the financial support given for this research. This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior – Brasil (CAPES) – Finance Code 001.

Appendix A Parameter tables

Tables 6 describe the parameters, initial conditions and time intervals used to generate synthetic data. The solve_ivp function from the SciPy library was used to numerically integrate these equations over time. The selection of these parameters was made in a case-by-case basis, with the intent of representing a realistic scenario for each system.

As for the SR algorithm’s settings, Table 8 and Table 9 list the employed parameters for the non-linear pendulum, Lotka-Volterra and Lorenz systems, and the epidemiological models, in that order. Any other available parameters that were not mentioned were kept as default. These settings were chosen to allow for the best performance possible, and kept the same when investigating the impact of noise in these algorithms (see subsection V.1).

[FIGURE:]