D-CODE: Discovering Closed-form ODEs from Observed Trajectories
Authors: Zhaozhi Qian
Krzysztof Kacprzyk
Mihaela van der Schaar
Venue: ICLR 2022 Spotlight (2022)
Keywords: Symbolic Regression
Ordinary Differential Equation

Abstract: For centuries, scientists have manually designed closed-form ordinary differential equations (ODEs) to model dynamical systems. An automated tool to distill closed-form ODEs from observed trajectories would accelerate the modeling process. Traditionally, symbolic regression is used to uncover a closed-form prediction function $a=f(b)$ with label-feature pairs $(a_i, b_i)$ as training examples. However, an ODE models the time derivative $\dot{x}(t)$ of a dynamical system, e.g. $\dot{x}(t) = f(x(t),t)$, and the "label" $\dot{x}(t)$ is usually *not* observed. The existing ways to bridge this gap only perform well for a narrow range of settings with low measurement noise, frequent sampling, and non-chaotic dynamics. In this work, we propose the Discovery of Closed-form ODE framework (D-CODE), which advances symbolic regression beyond the paradigm of supervised learning. D-CODE leverages a novel objective function based on the variational formulation of ODEs to bypass the unobserved time derivative. For formal justification, we prove that this objective is a valid proxy for the estimation error of the true (but unknown) ODE. In the experiments, D-CODE successfully discovered the governing equations of a diverse range of dynamical systems under challenging measurement settings with high noise and infrequent sampling.

Published as a conference paper at ICLR 2022
D-CODE: DISCOVERING CLOSED-FORM ODES FROM OBSERVED TRAJECTORIES
Zhaozhi Qian University of Cambridge zq22 cam accam.ac.uu uckuk
Krzysztof Kacprzyk University of Cambridge kk @cam. ac. ukk
Mihaela van der Schaar University of Cambridge, UCLA,
The Alan Turing Institute mv 72 cam. ac uk
ABSTRACT
For centuries, scientists have manually designed closed-form ordinary differential equations (ODEs) to model dynamical systems. An automated tool to distill closed- form ODEs from observed trajectories would accelerate the modeling process. Traditionally, symbolic regression is used to uncover closed-form prediction function = f(b) with label-feature pairs (a;, b,) as training examples. However, ODE models the time derivative (tt) of dynamical system, e.g. .i(t) fax(t), t), and the "label" x(t) is usually not observed. The existing ways to bridge this gap only perform well for a narrow range of settings with low measurement noise and frequent sampling. In this work, we propose the Discovery of Closed- form ODE framework (D-CODE), which advances symbolic regression beyond the paradigm of supervised learning. D-CODE uses novel objective functioo based on the variational formulation of ODEs to bypass the unobserved time derivative. For formal justification, we prove that this objective is valid proxy for the estimation error of the true (but unknown) ODE. In the experiments, D-CODE successfully discovered the governing equations of a diverse range of dynamical systems under challenging measurement settings with high noise and infrequent sampling.
INTRODUCTION
An ordinary differential equation (ODE) links the state of a continuous-time dynamical system ax() to its time derivative (t)) via function f. e.g. f+(t) = f(x=(t). t). The ODE closed-form when f has concise and analytical expression (e.g. [ar, t) -f log( Closed-form ODEs ubiquitous in science and engineering. The popularity is partly due to their transparency and interpretability to human experts (Petersen al., 2019). Their concise functional form also facilitates the analysis of many key properties of the dynamical system, e.g. asymptotic stability (Lyapunov, 1992).
However, the discovery of closed-form ODEs has been laborious and time-consuming as it heavily relies on human experts (Simmons, 1972). The goal of this work is to discover closed-form ODEs in an automated and data-driven way. We envision that such tool would accelerate the quantitative modeling of dynamical systems, with applications ranging from discovering the kinetics of bio- chemical reactions to modeling tumor growth (Schuster, 2019; Geng et al., 2017).
Symbolic regression is an established approach to discover closed-fom ffom function = {bb for predicting the label a from the feature b. It operates the supervised learning setting and requires a dataset of label-feature pairs for training, i.e. D = {a;; be J:W (Schmidt & Lipson, 2009).
However, symbolic regression is not directly applicable to ODEs because we usually do not observe the time derivative i(t), i.c. the labels for regression are not available in the data. In contrast, only have access to the measurements of state y(t) = z(t) + e at some discrete time points and with noise. For instance, the datasets for tumor growth modeling typically contain tumor volumes measured at different clinical visits with substantial noise ( (Wilkerson et al., 2017; Mazaheri et al., 2009). We refer to this problem as the equation data mismatch.
simple way to resolve this mismatch is to estimate (t) from data (c.g. by numerical differentiation (Bickley, 1941) )), and use the estimated derivatives labels for regression (Gaucel et al., 2014).

---

Published as a conference paper at ICLR 2022
However, it is very challenging to recover the time derivative under high measurement noise or infrequent sampling (Cullum, 1971), and symbolic regression tends to perform poorly with inaccurate labels (Zegklitz Pošík, 2021). This is an substantial drawback because many applications, such as healthcare, may involve noisy or infrequently sampled data (Jensen et al., 2014).
In this work, we develop the Discovery of Closed-form ODE framework (D-CODE), which extends symbolic regression beyond the supervised learning setting. The key insight behind D-CODE is the variational formulation of ODE (Hackbusch, 2017), which establishes a direct link between the trajectory xxt) and the ODE while bypassing the unobservable time derivative (t)). We develop novel objective function based on this insight, and prove that it is a valid proxy for the estimation error of the true (but unknown) ODE. We demonstrate via extensive experiments that D-CODE can uncover the governing equations for a diverse range of dynamical systems while being substantially more robust measurement artifacts than the alternative methods, Finally, D-CODE designed as general framework, where some of its components can be flexibly adapted based on the application.
2 BACKGROUND AND PROBLEM SETTING
Curating a dataset of trajectories for ODE discovery involves many decisions, e.g. What variables to include? When and for how long take measurements? (we discuss dataset curation in Appendix C.) In this work, we assume the dataset is given and the variables can be modeled by a system of first-order autonomous ODEs (Eq. 1). We note that higher-order or time-dependent ODEs can be represented in this form by curating a dataset with additional variables (Simmons, 1972) (Appendix B). The system with N+ variables is defined as
tt({) = (x(())) Vj I....,, J, Vt E [0,
where we use Newton's 's notation (t) for the time derivative. The functions ʃ RJ - RR will be sometimes referred to as the ODEs directly. We denote T E R as the maximum time horizon we have data for. We highlight the following distinction: the trajectory æj [0, + R is function of time, whereas the state xz (() E IR, vl T) is a point on the trajectory'. We denote the state vector at) = (t) (1, (t)
R and the vector-valued trajectory function æ [r,)
Let 's be the true but unknown ODEs to be uncovered, and z; : [0, 7] - R i N, N N+ be the true trajectories that satisfy fj's In practice, only measure the true trajectories at discrete times and with noise. Denote the measurement of trajectory at time as (t) E TR:'; we assume
y = (t) + (t), Vi N, TET
(2)
where (t) E R is zero-mean noise with standard deviation σ.. The measurements are made at time t (t1, .2........ T). We denote the dataset as D (t)l < N,t T).
Closed form. The function fy R + R has a closed if can be expressed as a finite sequence of operations (+, +, log....), input variables (I1, x2, and numeric constants (1.5, ..8....) (Bor- wein et al., 2013). The functional form of f, is the expression with all the numeric constants replaced by placeholders 02 's, e.g. (x) = 0, log(02x). To fully uncover fj,, we need to infer its functional form and estimate the unknown constants 0g 's (if any).
Variational formulation. The variational formulation provides a direct link between the trajectory and the ODE f without involving å (Hackbusch, 2017). We start with the following definition.
Definition 1. Consider J N+ T e R+, continuous functions [0, T]
f RJ R. and (o, T), where C¹ is the set of continuously differentiable functions. We define the functionals
C,(f,2,9)= 1'' (((()))(()dt+ + S." (,(()((t)dt; VjE (1,2,
Importantly, the functional C₁ depends on the testing function g(t) and its derivative g(t) but not i;, Proposition below provides the variational formulation of ODE, Essentially specifies the infinitely many constraints that x has to satisfy in order to be a solution to the ODE. (see Appendix A).
We restrict the domain of the function zy to [0, T] (rather than R because only have data up to T. However, after learning the ODE, we can extrapolate 2, beyond (Appendix D).
2

---

Published as a conference paper at ICLR 2022
Table 1: Comparison of related works. "Data": the observed variables. "Allowed f the space of discoverable functions. "Est.": the quantities estimated in the intermediate step. Free": is the method not reliant on à?? (=( Free": is the method not reliant on initial condition a00}? "Objective": the objective function. References: [1] Schmidt & Lipson (2009), [2] Brunton et al. (2016), [3] Gaucel et al. (2014), [4] Chen et al. (2018).
Method
Data
Allowed f
Est.
t Free
az a=0) Free
Objective
Symbolic Reg [1]
a, b
Closed-form
None
||a f(b)]
2-step Sparse [2]
y(t)
Σ 0ahn(x))
Σ, (ā()) - (()(t))
2-step Symbolic [3]
y(t)
Closed-form
[|(()) -
( (()(1))
ODE Approx [4]
v(L)
Neural nets
2(0)
x
L. [1y((() - #(01]22
D-CODE
(()
Closed-form
Equation
Proposition 1. (Hackbusch, 2017) Consider J N+ T RRt, continuously differentiable function
a : [0, T) + R. and continuous functions Jj : TR. + R. j = =1,...... I.....J. Then x is the solution to
the system of ODEs in Equation if and only if
Cy(J,, = 0. e (l..... Vg E Co g(0) -(T))= =0
RELATED WORK
3.1 SYMBOLIC REGRESSION
Symbolic regression attempts to uncover a closed-form prediction function f(b using supervised learning (Table 1). The main challenge is optimization: searching for the optimal f is thought to be NP-hard because the space of closed-form functions is vast and complex it is combinatorial in the functional form and continuous in the numeric constants (Lu et al., 2016). Therefore, most existing work focuses developing optimization algorithms. Genetic programming has been one of the most successfully and widely used methods (Koza, 1994; Schmidt & Lipson, 2009). It represents tree where the internal nodes operations and the leaves are variables constants (Figure 1 B). It then applies genetic algorithm to search for the best tree representation (Forrest, 1993). More recently, AI Feynman introduces set of heuristics directly prune the search space, thereby improving the optimization efficiency (Udrescu Tegmark, 2020). Optimization methods based on reinforcement learning (Petersen al., 2019), pre-trained neural networks (Biggio et al., 2021), and Meijer G-functions (Alaa & van der Schaar, 2019) have also been proposed.
All methods above use prediction error (e.g. RMSE) the objective, which is not applicable to ODE discovery because the label (time derivative) is not observed. D-CODE advances symbolic regression beyond the supervised learning setting and addresses the unique challenges for ODE discovery: the cquation-data mismatch and the sensitivity to noisy or infrequent observations.
3.2 DATA-DRIVEN DISCOVERY OF CLOSED-FORM ODES
Two-step sparse regression performs sparse regression on the estimated time derivative ū, (Brunton et al., 2016; Rudy et al., 2017). It assumes that the true ODE can be symbolically expressed in linear form: (t) Ek_ hyh (x(t)) where f R are unknown constants and h2 R R pre-specified candidate functions, such monomials 2, x² .... (Table 1). The method uses L1 regularization on 0ş to ensure the learned function only contains a few terms. We review the methods to estimate time derivative in Appendix E.
However, the assumed linear form is very restrictive (e.g. 3/2 or 1/(x + 2) are not allowed). In fact, it cannot symbolically represent many well-known ODEs (e.g. generalized logistic model). Furthermore, the choice of candidate functions heavily depends on the human experts, which deviates from our goal of automated and data-driven ODE discovery.
Two-step symbolic regression also uses the estimated time derivatives ②, as the label but applies symbolic regression instead (Gaucel et al., 2014) (Table 1). It may any optimization method discussed Section 3.1 to search for the optimal function. In this way, removes the linear assumption the functional form and can discover closed-form functions in general.
3

---

Published as a conference paper at ICLR 2022
(A)
y(t)
D-CODE
Estimated Trajectory &
(B)
01 a log(02 =)
Observed Trajectory
(1)
Equation
Two-step Symbolic Regression
Estimated Derivative xit]
RMSE (Table 1)
0
1.5
0.8
Figure 1: (A) Illustration of D-CODE (compared with two-step symbolic regression). Both algorithms involve two steps, but they (1) estimate different variables in step one and (2) optimize different objective functions in step two. (B) An example of closed-form function and its tree representation.
The main difficulty is that the true label (the derivative) may not be accurately recoverable due to measurement and discretization error, but symbolic regression is sensitive to inaccurate labels because of the large search space (Agapitos et al., 2012). As we show experimentally in Section 5, even with the state-of-the- relatively low. We briefly review the methods for estimating derivatives from y(() in Appendix E and show why this is generally a very challenging problem in Section 4.2.
methods for estimating derivatives, this approach still often fails when the noise is
The function approximator approach learns the true ODE with a function approximator f. such as a neural network (NN) (Chen al., 2018; Rubanova al., 2019) Gaussian process (GP) (Heinonen et al., 2018). It estimates the unknown initial condition @(0) and predicts the entire trajectory z by solving the initial value problem (IVP) of the approximated ODE f. The model is trained by minimizing the prediction error between predicted and the observed trajectories (Table 1).
Importantly, the function approximator approach does not give concise closed-form expression to describe the dynamics, which is the main objective of this work. Furthermore, its performance strongly depends on the prediction horizon used for training. Using an inappropriate horizon may fail with certain type of systems, e.g. chaotic systems (Ditto & Munakata, 1995). We further discuss this issue and provide experimental results in Appendix H.2.
4
METHOD
4.1 THE D-CODE ALGORITHM
D-CODE consists of a preprocessing step and optimization step. Figure provides schematic illustration of D-CODE and Appendix B provides the pseudocode. The two steps are detailed below.
Preprocessing. D-CODE starts by estimating the trajectories G; [0, + R < N from the noisy and discretely-ammpled data D as an approximation to the true trajectories æ;, Denoising and interpolating signals a well-studied problem in statistics and signal processing, with many proven solutions including Gaussian process and spline regression (Bernardo et al., 1998; Marsh & Cormier, 2001). The D-CODE framework agnostic to the exact choice the smoothing algorithm. The user should choose a suitable algorithm based on the application and adopt best practices (e.g. Is the noise distribution Gaussian? Are the measurements made at regular intervals? etc.).
Optimization, After estimating the trajectories R;, we search for the function Fy. Vj < that is consistent with &. Specifically, we solve the following optimization problems for all j J:
N ΣV
f, = arg min
M C.
(5)
4

---

Published as a conference paper at ICLR 2022
where Cj is the functional defined in Equation 3 and computed using numerical integration (Davis & Rabinowitz, 2007). The search space of f is the set of closed- form functions, and the user may instantiate D-CODE with any optimization algorithm proposed for symbolic regression (Section 3). We use a set of pre-defined testing functions E C² T). 8 S $ (to be discussed in Section 4.2).
Comparison with two-step symbolic regression. As depicted in Figure D-CODE and two-step symbolic regression have two main differences. (1) The estimands in the pre-preee pre-prccsssing steps are different: D-CODE estimates the true trajectory while the two-step symbolic regression estimates the time derivative. (2) The optimization objectives are different: D-CODE optimizes the loss function in Equation 5 while the two-step symbolic regression optimizes the prediction loss, e.g. RMSE.
4.2 THEORETICAL RESULTS
In this section we provide a formal justification of the objective in Eq. 5 and describe the shortcomings of the objectives used in other methods (Table 1). For ease of exposition, we assume that our dataset contains only one trajectory, i.e. N 1. The results can be easily extended to multiple trajectories.
Distance between f and fr. A reasonable objective should measure the distance between a candidate function f and the target function f (i.e. the true ODE). By minimizing the distance, the candidate f would better approximate f". A common way to measure the distance between functions is to use the metries induced by L" norms of the function space. However, the I," norms consider the values in the entire domain R. This is disadvantage because most dynamical systems in nature operate within specific range and scale, and R. contains areas where we do have knowledge about or cannot collect data from. A natural solution is to restrict the comparison on the trajectory z, which can be achieved by composing f and f with æ. We propose the following distance function:
(f.ff)) == II o - oxl12- |K - f*) æ112
(6)
where o denotes function composition and æ is a true trajectory satisfying the ODE f". In Appendix A, further justify and discuss the properties of (J, However, in practice, the dz (S, cannot be computed from data because it depends on the unknown f and æ. Next, we show that the objective in Eq. 5, which computable from data, can be used approximate this distance.
Convergence to distance. The first step of our algorithm estimates a trajectory 2 from the measure- ments (t)l E 7 as an approximation to the true trajectory æ. For suitable smoothing algorithm, if we increase the number of samples on the trajectory, the recovered trajectory should converge to the ground truth. This is achieved, for instance, by spline regression (Stone, 1985; 1994) or Gaussian processes Choi & Schervish (2007). If this is the case, Theorem shows that given high enough sampling frequency and large enough number of testing functions our objective converges to the squared distance ff which measures how different ff is from f*. That justifies our objective.
Theorem 1. Consider J E N+,
.J)......].. E RR., Let f R + R be a continuous functton, and let x: : [0, + I be a continuously differentiable function satisfying Lt(t) = f (t))).
Consider a sequence of functions (āx)) where et : [0. 2] + IR a continuously differentiable function. If (mx) converges to x in L norm. Then for any Lipschitz continuous function f
lim lim Σ Ec,(F.n.9.)) = d)).
8=1
where f91.92,.... .2..... and gi E
} is Hilbert (orthonormal) basis for I.² T such that Vi, 9:(0) (T) 0
We provide the proof of Theorem in the Appendix A.
Comparison with other methods. Theorem does not place any additional constraints regarding the smoothing algorithm apart from the convergence of the regression function. In general, the convergence of functions does not imply the convergence of the derivatives (example in Appendix A). That means that in contrast to D-CODE, other objectives in Table place constraints on the type estimation algorithm, i.e., the estimated derivative d needs to converge to the true derivative f. Even if the estimation algorithm satisfies this constraint (for instance, spline regression (Zhou & Wolfe, 2000)), the convergence rate for derivatives is slower than for the functions themselves (Stone, 1982). Convergence rate depends on the smoothness (differentiability class) of the function
5

---

Published as a conference paper at ICLR 2022
and differentiation reduces its differentiability class. Namely, if true trajectory has maximum number of m continuous derivatives then its derivative è has only m continuous derivatives.
Choice of testing functions. As stated in Theorem 1, should use testing functions that are orthonormal basis for [0 I"] and satisfy g(0) = g(T) = 0. Ideally, they should also have analytical time derivatives gt) for efficient computation of the functional C' (Eq. 3). Hence, natural choice the sine functions (t) = 22/T sin(stt/ which we will use in all the main experiments. Another possibility is the cubic spline functions (investigated in Appendix G). Although Theorem 1 uses infinitely many testing functions to establish convergence, we observe in the experiments that typically 40-60 sine functions are sufficient (a sensitivity analysis is conducted in Appendix G).
5 EXPERIMENTS AND EVALUATION
In this section, we perform a series of simulations to evaluate whether the algorithms can discover the underlying closed-form ODEs that govern the observed trajectories. 2
Choice of dynamical systems. In this study, we select five dynamical systems governed by nine closed-form ODEs (one of them is discussed and examined in Appendix G due to space limit). The selected systems have varying complexity and different properties of the temporal dynamics. We start with two common growth models, the Gompertz model and the generalized logistic model (Gompertz, 1825; Richards, 1959). Both models involve one variable governed by a nonlinear ODE and converging to a global fixed point. Next, we consider the glycolytic oscillator in bio- chemistry (Sel' Kov, 1968), which is a standard benchmark problem for dynamical system prediction and inference (Daniels & Nemenman, 2015a;b). The system involves two nonlinearly interacting variables that converge oscillatory limit cycle. Finally, consider the chaotic Lorenz system, which involves three variables forming a strange attractor (Lorenz, 1963). Together, these systems represent a range of temporal dynamics and application scenarios.
Mcasurement settings. For each dynamical system, we consider different measurement settings specified by (1) the measurement noise level an, (2) the sampling step size At, and (3) the number of trajectories N. Since the dynamical systems have different scales, we will report the noise-to- signal ratio OR //std( (x (t)) for ease of comparison. We sample the system at regular intervals T = (t 2 ... T to show the effect of changing sampling frequency. For each setting, we perform 100 independent simulation runs to compute the evaluation metric and its confidence interval.
Data generation. For each trajectory, we first sample the initial condition x, (0) from a uniform distribution (specified in Appendix F). Then we obtain the true x; by solving the IVP computationally, Finally, we obtain y: (t) by adding independent Gaussian noise and sampling at discrete time steps.
Evaluation metrics. We use three metrics to study various aspects of the algorithms. (1) The probability of successfully recovering the functional form (Success Prob.). We use a computer algebra system to judge whether the true and discovered functional forms are equivalent (Meurer et al., 2017) (e.g. 0, + xi æ + 0,). (2) When the functional form is correct, we evaluate the accuracy of the estimated constants ê using root mean square error (RMSE). (3) We evaluate how well the estimated f apprrximmaes f* using the distance function in Equation (Dist.). Note that it is possible that f has a wrong functional form but still approximates f" well.
Algorithm and Benchmarks. We instantiate D-CODE with Gaussian process for smoothing and genetic programming for optimization due to their success and popularity in the literature (Gramacy, 2020; Schmidt & Lipson, 2009). Our main benchmark is the two-step symbolic regression method. We consider three variants with different methods for estimating the derivative: Total variation regularized differentiation (SR-T) (Chartrand, 2011), 1), & Abel, 2007), and Gaussian process smoothed differentiation (SR-G) (detailed in Appendix E). Note that SR-G and D-CODE both fit Gaussian process in the pre-proce processessing step but they use different objective functions. For fairness of comparison, we use genetic programming for optimization in all these benchmarks. We also compare with Neural ODEs in some cases, although does not give closed-form expressions. Appendix F contains further implementation details.
smoothed differentiation (SR-S) (Ahnert
Thhe code is available at wtthhttt////g/ http:///ggthuub..oo///hhhhhh QI AN/ COCOR- ICLR- 2022.

---

Published as a conference paper at ICLR 2022
Gompertz Model
Success Prob. Success Prob.
Varying noise level dR
Varying step size At
Varying sample size
0.5
0.5
0.5
0.0
0.0
10
10
10
50
100
Generalized Logistic Model
1.0
1.0
Methods
0.5-
0.5
0.5-
D-CODE
SR.T
SR-S
10
10-3
10a 10
100
50
75
100
SR-G
OR
Figure 2:
The success probability of the two growth models under different settings. The three columns correspond to various noise levels dR, step sizes At, and numbers of trajectories N. The rows correspond the Gompertz model and generalized logistic model. The shaded area is the 95% confidence interval.
(A) Trajectory
(B) Derivative
(C) Estimated Derivative
0.6
0.4
0.4
X(U
Xit
High noise
0.4
0.2
SR-5
Blas
0.0
0.0
3
3
2 Time
Time
Time
Figure 3: Illustration of performance gain using trajectory with moderate noise oR 0.1. Panel A: the trajectory #(t) estimated from the measurements y(t) very close to the true trajectory a+t); hence, the first step of D-CODE tends to have small crror. Panels B and C: even with moderate noise on x (f), the estimated derivatives (t)) may suffer from very high variance (B: finite difference) or systematic bias (C: SR-T, SR-S). Hence, two-step symbolic regression tends to under-perform due to high emror in the first step.
5.1 RESULTS
Gompertz and generalized logistic growth models. These two models are widely used to capture asymmetric growth with saturation. Their governing equations are as follows:
a(t) = =00x(()) log (0zr(1) Gompertz Model a(t) = (tx xt))
Generalized Logistic Model
where 0,,0; are numeric constants. It is worth highlighting that the generalized logistic model involves a constant 0g E R+ in the power of rff). Therefore, it canmot be symbolically expressed in sparse linear form Σ k= hy (x(t)) for some candidate functions hg, and the two-step sparse regression method cannot discover the exact functional form of this equation (Section 3.2).
The main simulation results presented in Figure 2 (the RMSE and distance metrics follow the same pattern and are shown in Appendix G. The discovered ODEs are shown in Appendix G Table 4). D-CODE achieves comparable performance as the benchmarks the low noise, small step size settings (aR = 0.01), where the time derivative i can be accurately estimated from the observed trajectories. However, D-CODE consistently out-performs the benchmarks under the more challenging measurement settings. With increased noise or step size, the success probabilities decrease much more rapidly the benchmarks compared D-CODE. Note that although SR-G and D-CODE both use Gaussian process in the pre-processing step, D-CODE achieves much stronger performance. This suggests that D-CODE's performance gain is mainly due to its objective function.
Figure 3 contextualizes the reason for performance gain with an illustrative example. The panel shows a typical observed trajectory with moderate measurement noise (aR = 0.1). Here, the Gaussian process can produce good estimate @(() by de-noising the measurements. Its posterior interval covers the true trajectory (t). Hence, ((t) is suitable to guide the search for the ODE, is done in D-CODE. On the other hand, the moderate measurement noise xtt) translates into very high noise on the time derivative (t)) (panel B). Even with the state-of ofate-of-the-ar numerical differentiation methods used by SR-T and SR-S, the estimated derivative is still biased on certain
7

---

Published as a conference paper at ICLR 2022
Table 2: Simulation results of the glyeolytic oscillator under different noise levels. The success prohability and the RMSE on 0,, 0g are reported for the two equations. Standard deviations are shown the brackets.
Success Prob.
RMSE (10-3)
Equation
Method
OR 0.01
0.1
0.2
OR 0.01
0.1
0.2
SR-T
0.45 (.05)
0.27 (.05)
0.00 (.00)
1.19 (.15)
2.37 (.38)
NA
SR-S
0.44 (.06)
0.11 (.04)
0.00 (.00)
1.67 (.24)
2.07 (.47)
NA
Eq. 8
SR-G
0.44 (.05)
0.05 (.02)
0.00 (.00)
1.87 (.22)
2.18 (.46)
NA
D-CODE
0.58 (.05)
0.51 (.05)
0.26 (.05)
1.01 (.12)
1.55 (.28)
2.00 (.27)
SR-T
0.99 (.03)
0.88 (.03)
0.00 (.00)
0.14 (.03)
0.43 (.03)
NA
SR-S
0.95 (.02)
0.04 (.02)
0.00 (.00)
(C10)
0.41 (.18)
NA
Eq.
SR-G
(.01)
0.94 (.02)
0.25 (.04)
0.22 (.02)
1.10 (.07)
2.85 (.02)
D-CODE
1.00 (.00)
0.91 (.03)
0.65 (.05)
0.07 (.01)
0.40 (.07)
0.61 (.05)
Success Prob. (t(t)
Success Prob. Xy(t)
Success Prob. Xy(t)
1.0
1.0
1.0
0.5
0.5
0.5
0.0
0.10
0.15
0.20
0.25
0.30
0.10
0.15
0.20
0.25
0.30
0.10
0.15
0.20
0.25
0.30
Noise
Noise
Noise
Ground truth
D CODE
SR-T
Neural ODE
40
40
20
20
20
20 0
20
20
20 20
20
0
20
0
Figure 4: Simulation results of chaotic Lorenz system. First row: success probabilities the three equations under different noise levels. Second row: simulated trajectories using true and estimated equations.
intervals (panel C). Thus, applying symbolic regression on the biased estimates is unlikely to perform well, is shown in the benchmark results. The same reasoning also applies to the scenario with large sampling step size At, another challenging situation for estimating the time derivatives.
Glycolytic oscillator. This biochemical oscillator consists of two ODEs
1968):
c (1) 0, 02 (1)=x=1111))))))))))))))
2b (1) (1) 0321
(1)²
where 01, 22 03 > O are numeric constants. Compared with the growth models above, the glycolytic oscillator involves more variables and has more complex expressions. The variables in the equation interact in a highly nonlinear way through the term .r (t)
The simulation results for the glycolytic oscillator shown Table Here we report the success probability in addition to the RMSE (x 10-2) on the constants 02, 03. whose true values are set to 0.1 (other metrics are reported in Appendix G). D-CODE performs better or equally well compared to the benchmarks in all metrics. Note that D-CODE is the only model that successfully recovered the functional form under nn 0.2 (The benchmarks' RMSE are NA this setting because they never recovered the correct functional form).
Chaotic Lorenz system. The Lorenz system is model system for chaotic dynamics, defined as: (t( = 0, (x2(t) - =t))): itt() = .z (t) (02 - xz((t) x2(();; (t( z (t) (yxx(()
We set the constants 0 10,02 28, 0g 8/3, which leads chaotic behaviour.
8

---

Published as a conference paper at ICLR 2022
The simulation results for the chaotic Lorenz system shown in Figure 4. D-CODE's robustness to measurement noise is still evident when the underlying system is chaotic and non-periodic. The benchmarks achieve high success probability only when the measurement noise is low. The second row of Figure 4 shows simulated trajectory in the phase space using the true or estimated ODEs (with n = 0.2). D-CODE generates trajectory that is closest to the true one. SR-T fails because it is unable to recover the true ODE under this noise level. Although the Neural ODE is trained on very short prediction horizon (0.04), it still fails to capture the non-periodic behavior of the system.
D-CODE IN ACTION
We apply D-CODE to model the temporal effect of chemotherapy on the tumor volume. For most patients, the tumor volume decreases initially in response to the treatment but it eventually grows again due drug resistance. Although there exist many models for unintervened tumor growth, few models can capture the effect of chemotherapy or the possibility of relapsing--nn many cases exponential decay (log cell-kill model) is used for crude approximation (Geng al., 2017).
We used the dataset collected by Wilkerson et al. (2017) based on eight clinical trials on cancer patients. The patients started chemotherapy at time = 0, and the trajectory of tumor volume are measured for up to one year 7 = 1. In practice, the tumor volumes are inferred from medical scans and are subject to significant measurement noise (Mazaheri et al., 2009). After data cleaning, we obtained 3 310 trajectories, among which 200 is used for training and 1 110 for evaluation.
The following two ODEs are discovered by D-CODE and SR-T.
(tt) = 4.48t (t) ) + log(t); D-CODE
Observed
<(t) = (t) log (t)) log (Lx(!) + 21); SR-T
0.-
The equation discovered by SR-T is much more complex than the one found by D-CODE while still being a poor representation of the data. Figure 5 shows an observed trajectory and the trajectories generated by the two ODEs. D-CODE successfully captures the relapse at = 0.6 and predicts that, under no further intervention, the tumor volume would keep growing after > 0.6 until patient death. In comparison, SR-T's trajectory is almost constant.
0.0
0.0
0.2
0.6
1.0
1.0
D-CODE
0.0
02
0.4
0.5
0.8
1.0
5R-T
0.6
We also quantitatively evaluate the prediction errors on the testing set. Given the tumor size measured t 0, solve the IVP obtain the full trajectory x and compute the RMSE between the est- imated and the observed trajectory. D-CODE achieves error of 0.220 compared with 0.637 of SR-T.
Figure 5: An observed trajectory and the trajectories generated by the two discovered ODEs. Values are indexed at O.
7 DISCUSSION ON FAILURE MODES AND OPEN CHALLENGES
Scientific discovery is a challenging process with no success guarantee. The same holds for ODE discovery with D-CODE. Here we discuss the failure modes and the challenges for future work.
Unobserved variables. In some cases, certain external variables u(t) E R may also influence the temporal dynamics, i.e. (t) = (=((), u(t)). When u(t) is observed or known, D-CODE can include these variables to discover the expanded equation. However, when they are unobserved, the identification of system dynamics challenging or even impossible without strong assumptions (Ljung. 1998). In this case, D-CODE can still find an ODE based on the observed variables- ODE may still be a useful approximation but it will leave some dynamics unexplained.
Complex equations. The discovery algorithm may fail to discover highly complex ODEs. The complexity can arise from two aspects: (1) x(t) may include many variables (high dimensionality) or (2) the equation f may include many mathematical operations (long expression). However, highlight that these settings are difficult in general, even for human experts with domain knowledge.
Extreme measurement settings. D-CODE may still fail when the measurement noise step size is too large. However, as we show experimentally, D-CODE is more robust to measurement artifacts than the alternative discovery methods.

---

Published as a conference paper at ICLR 2022
ACKNOWLEDGMENTS
This work is supported by the US Office of Naval Research (ONR), the National Science Foundation (NSF Grant number: 1iter:1722511 1722516), and Roche. We thank the anonymous reviewers, Professor William Zame and the members of vanderschaar-lab for many insightful comments and suggestions.
Ethics Statement. In this work, we use a tumor growth dataset that is anonymized and licensed for research purpose (Section 6). The dataset has also been used in prior research in the literature (Wilkerson et al., 2017). We discuss the limitations and open challenges for D-CODE in Section 7. We stress that D-CODE is not intended to or capable of fully replacing human experts in modeling dynamical systems. We envision that D-CODE would facilitate the modeling process by automatically extracting closed-form ODEs from observed trajectories allowing human experts to easily inspect, analyze and refine the distilled ODEs.
The proof of theoretical results is shown in Appendix A. The technical assumptions for the proof listed in Proposition 1 and Theorem 1 in the main text and further discussed in Appendix A. The experiment settings are discussed in Section 5 and and further detailed in Appendix F. The experiment code is available at https: /ett://gittttt /ggtttu.......///gg com/ zhaczhiQIAn D- I-CODE-ICLR--0222
Reproducibility Statement.
REFERENCES
Alexandros Agapitos, Anthony Brabazon, and Michael O' O'Neill. Controlling overfitting in symbolic regression based on bias/variance error decomposition. In International Conference on Parallel Problem Solving from Nature, pp. 438-447 Springer, 2012.
Karsten Ahnert and Markus Abel. Numerical differentiation of experimental data: local versus global methods. Computer Physics Communications, 177(10):764-7777 774, 2007.
Ahmed M Alaa and Mihaela van der Schaar. Demystifying black-box models with symbolic metamodels. Advances Neural Information Processing Systems, 32:1 14111304-11114 2019.
J Bernardo, Berger, APAFMS Dawid, A Smith, et al. Regression and classification using gaussian process priors. Bayesian statisties, 6:475, 1998.
WG Bickley. Formulae for numerical differentiation. The Mathematical Gazette 25(26)):19-27, 1941.
Luca Biggio, Tommaso Bendinelli, Alexander Neitz, Aurelien Lucchi, and Giambattista Parascandolo. Neural symbolic regression that scales. arXiv preprint urXiv:2106.06422 2021.
Bernd Blasius, Lars Rudolf, Guntram Weithoff, Ursula Gaedke, and Gregor F Fussmann. Long-term cyclic persistence an experimental predator-prey system. Nature, 52267789) 2020.
Jonathan M Borwein, Richard E Crandall, et al. Closed forms: what they are and why we care.
Notices of the AMS, 60(1):00-65, 2013.
Felix E Browder. Fixed point theory and nonlinear problems. Proc. Sym. Pure. Math, 39:49-88, 1983.
Nicolas Brunel and Peter E Latham. Firing rate of the noisy quadratic integrate-and-fin neuron. Neural computation, ((10))2281-2306 2003.
Nicolas JB Brunel, Quentin Clairon, and Florence d' Alché Buc. Parametric estimation of ordinary dif- ferential equations with orthogonality conditions. Journal of the American Statistical Association, 109(505):733-88, 2014.
Steven L Brunton, Joshua L Proctor, and J Nathan Kutz. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the national academy of sciences, 113(15):332-3393) 2016.
Rick Chartrand. Numerical differentiation of noisy, nonsmooth data. International Scholarly Research Notices, 2011, 201
10

---

Published as a conference paper at ICLR 2022
Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud. Neural ordinary differential equations. arXiv preprint arXiv:1886.07366, 2018.
Taeryon Choi and Mark J. Schervish. On posterior consistency in nonparametrie regression problems. Journal of Multivariate Analysis, 98(10):199- 1987, November 2007.
S-N Chow and Jack K Hale. Methods of bifurcation theory, volume 251. Springer Science & Business Media, 2012.
Jane Cullum, Numerical differentiation and regularization. SIAM Journal on numerical analysis, (2):254 265, 1971.
Bryan C Daniels and Ilya Nemenman. Automated adaptive inference of phenomenological dynamical models. Nature communications, 6(1):1-8, 201
Bryan C Daniels and llya Nemenman. Efficient inference of parsimonious phenomenological models of cellular dynamies using s-systems and alternating regression. PloS one, 10(3):e01 19821, 201 5b.
Philip Davis and Philip Rabinowitz. Methods of numerical integration. Courier Corporation, 2007,
Carl De Boor. A practical guide to splines, volume 27. springer-verlag New York, 1978.
William Ditto and Toshinori Munakata. Principles and applications of chaotic systems. Communica- tions of the ACM, 38(11);96-100)) 1995.
Arnaud Doucet, Simon Godsill, and Christophe Andrieu. On sequential monte carlo sampling methods for bayesian filtering. Statistics and computing, 10(3):97-2088, 2000.
L. E. Elsgolc. Calculus of variations / L. E. Elsgolc. International series of monographs in pure and applied mathematies Volume 19, Oxford, New York, Pergamon Press 1961, 1961. ISBN 1-4831-3756-2.
William F. Ford and James A. Pennline. When does convergence in the mean imply uniform convergence? The American Mathematical Monthly, 114(1):58-66 60, 2007.
Stephanie Forrest. Genetic algorithms: principles of natural selection applied to computation. Science, 261(5123):872-11) 878. 1993.
Sébastien Gaucel, Maarten Keijzer, Evelyne Lutton, and Alberto Tonda. Learning dynamical systems using standard symbolic regression. In European Conference Genetic Programming, pp. 25-36. Springer, 2014.
Changran Geng, Harald Paganetti, and Clemens Grassberger. Prediction of treatment response for combined chemo-and radiation therapy for non-small cell lung cancer patients using a bio- mathematical model. Scientific reports, 7(1):1-12, 2017.
Nitin S Gokhale. Practical finite element analysis. Finite infinite, 2008
Benjamin Gompertz. Xxiv. on the nature of the function expressive of the law of human mortality, and new mode of determining the value of life contingencies. in a letter to francis baily, esq. frs &c. Philosophical transactions of the Royal Society of London, pp. 513-583, 1825.
Robert B Gramacy. Surrogates: Gaussian process modeling, design, and optimization for the applied sciences. Chapman and Hall/CRC, 2020.
Shota Gugushvili and Chris AJ Klaassen.
Vri-consistent parameter estimation for systems of ordinary differential equations: bypassing numerical integration via smoothing. Bernoulli, 18(3):1061-1088 2012.
Wolfgang Hackbusch. Variational Formulation, pp. 159-180. Springer Berlin Heidelberg, Berlin, Heidelberg, 2017. ISBN 978-3-666-554661-2 doi: 10.
Wassim Haddad and VijaySekhar Chellaboina. Nonlinear dynamical systems and control. Princeton university press, 2011.
11

---

Published as a conference paper at ICLR 2022
Markus Heinonen, Cagatay Yildiz, Henrik Mannerström, Jukka Intosalmi, and Harri Lähdesmäki. Learning unknown ode models with gaussian processes. In International Conference on Machine Learning, pp. 1959-1968. PMLR, 2018.
Anders Boeck Jensen, Pope L Moseley, Tudor Oprea, Sabrina Gade Ellesøe, Robert Eriksson, Henriette Schmock, Peter Bjødstrup Jensen, Lars Juhl Jensen, and Søren Brunak. Temporal disease trajectories condensed from population-wis widd registry data covering 6.2 million patients. Nature communications, 5-1):1-10, 2014.
Holger Kantz and Thomas Schreiber. Nonlinear time series analysis, volume 7. Cambridge university press, 2004.
John R Koza. Genetic programming means for programming computers by natural selection. Statistics and computing, 4(2):87-112, 1994.
Lennart Ljung. System identification. In Signal analysis and prediction, pp. 163-173. Springer, 1998.
Edward N Lorenz. Deterministic nonperiodic flow. Journal of atmospheric sciences, 2002):130-141,
1963.
Qiang Lu, Jun Ren, and Zhiguang Wang. Using genetic programming with prior formula knowledge to solve symbolic regression problem. Computational intelligence and neuroscience, 2016, 2016.
Aleksandr Mikhailovich Lyapunov. The general problem of the stability of motion. International journal of control, 55(3):531-534 1992.
Niall M Mangan, Steven L Brunton, Joshua L Proctor, and J Nathan Kutz. Inferring biological net- works by sparse identification of nonlinear dynamics. IEEE Transactions on Molecular, Biological and Multi-Scale Communications, 2(1):52 63, 2016.
Lawrence C Marsh and David R Cormier. Spline regression models. Sage, 2001.
Yousef Mazaheri, Hedvig Hricak, Samson W Fine, Oguz Akin, Amita Shukla-Dave, Nicole M Ishill, Chaya S Moskowitz, Joanna E Grater, Victor E Reuter, Kristen Zakian, al. Prostate tumor volume measurement with combined t2- t2-weithted imaging and diffusion liffusion-weighhte mr: correlation with pathologic tumor volume. Radiology, 252(2):449-45 457, 2009.
Volker Ludwig Mehrmann. The autonomous linear quadratic control problem: theory and numerical solution, volume 163. Springer, 1991.
Aaron Meurer, Christopher P Smith, Mateusz Paprocki, Ondfej Čertík, Sergey Kirpichev, Matthew Rocklin, AMiT Kumar, Sergiu Ivanov, Jason Moore, Sartaj Singh, et al. Sympy: symbolic computing python. PeerJ Computer Science, 3:3103, 2017.
Brenden K Petersen, Mikel Landajuela Larma, T Nathan Mundhenk, Claudio Santiago, Soo K Kim, and Joanne T Kim. Deep symbolic regression: Recovering mathematical expressions from data via risk-sceking policy gradients. arXiv preprint arXiv: 191..0v9112047777 2019.
Zhaozhi Qian, Ahmed Alaa, Alexis Bellot, Mihaela Schaar, and Jem Rashbass. Learning dynamic and personalized comorbidity networks from event data using deep diffusion processes. In International Conference on Artificial Intelligence and Statistics, pp. 3295-3305. PMLR, 2020.
Markus Quade and Andy Goldschmidt. Numerical differentiation of noisy time series data in python.
https: //ppi//yp..a////////////e. wwwgooooooooooooooo hwtgg/proecc..oooo prroject/deriveee derivattiiiii viriiiittiiiii... 2020.
Jim O Ramsay, Giles Hooker, David Campbell, and Jiguo Cao. Parameter estimation for differential equations:
generalized smoothing approach. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 69(5):741-796, 2007.
FJ Richards. A flexible growth function for empirical use. Journal of experimental Botany, 10(2): 290-301, 1959.
Yulia Rubanova, Ricky TQ Chen, and David K Duvenaud. Latent ordinary differential equations for pregularly-samphe time series. Advances in Neural Information Processing Systems, 32: 5320-5330, 2019.
12

---

Published as a conference paper at ICLR 2022
Samuel H Rudy, Steven L Brunton, Joshua L Proctor, and Nathan Kutz. Data-drives discovery of partial differential equations. Science Advances, 3(4):e 16066144 2017.
Dominik Schmidt, Georgia Koppe, Zahra Monfared, Max Beutelspacher, and Daniel Durstewitz. Identifying nonlinear dynamical systems with multiple time scales and long-range dependencies. In International Conference on Learning Representations, 2020.
Michael Schmidt and Hod Lipson. Distilling free-form natural laws from experimental data. science, 324(5923):81 -888, 2009.
Peter Schuster. What is special about autocatalysis? Monatshefte für Chemie-Chemical Monthly, 150 (5):663-775, 2019.
EE Sel'Kov. Self-oscillations in glycolysis 1. simple kinetic model. European Journal of Biochem- istry, 4(1):79-86, 1968.
Xiaotong Shen, DA Wolfe, and Zhou. Local asymptotics for regression splines and confidence regions. The annals of statistics, 2665):170--1782, 1998.
George Finlay Simmons. Differential equations, with applications and historical notes. New York, McGraw-Hill. 1972.
T Stephens. gpleam: Genetic programming in python, with scikit-learn inspired and compatible api. https: http//ppppaaaaaaeeee/eee/.e readthedocs io/ io/en///eee stable/ index. html, 2021.
Charles J Stone. Optimal global rates of convergence for nonparametric regression. The annals of statistics, pp. 1040-1053, 1982. Publisher: ISTOR.
Charles J. Stone. Additive regression and other nonparametric models. The Annals of Statistics, 13 (2):689-705, 1985.
Charles J. Stone. The use of polynomial splines and their tensor products in multivariate function estimation. The Annals of Statistics, 22(1):1 18-171, 1994.
David Strong and Tony Chan. Edge-preserving and scale-dependent properties of total variation regularization. Inverse problems, 19(6):S165, 2003.
Silviu-Marian Udrescu and Max Tegmark. Ai feynman: physics-inspired method for symbolic regression. Science Advances, 6(16):caay2631, 2020.
Greg Welch, Gary Bishop, al. An introduction to the Kalman filter. Chapel Hill, NC, USA, 1995.
Julia Wilkerson, Kald Abdallah, Charles Hugh-Jones, Greg Curt, Mace Rothenberg. Ronit Simantov, Martin Murphy, Joseph Morrell, Joel Beetsch, Daniel Sargent, et al. Estimation of tumour regres- sion and growth rates during treatment patients with advanced prostate cancer: retrospective analysis. The Lancet Oncology, 18(1):143-155) 2017.
Jan Žegklitz and Petr Pošík. Benchmarking state-of the-art symbolic regression algorithms. Genetic programming and evolvable machines, 22(1):5-33, 2021.
Shanggang Zhou and Douglas A Wolfe. On derivative estimation in spline regression. Statistica Sinica, pp. 93-108, 2000.
Olgierd Cecil Zienkiewicz, Robert Leroy Taylor, Perumal Nithiarasu, and JZ Zhu. The finite element method, volume 3. McGraw-hin hill London, 1977.
13

---

Published as a conference paper at ICLR 2022
TABLE OF SUPPLEMENTARY MATERIALS
1. Appendix A: proof and discussion on the theoretical results.
2. Appendix B: implementation details of D-CODE.
3. Appendix C: discussion dataset curation.
4. Appendix D: utilities the discovered closed-form ODE.
5. Appendix E: background on numerical differentiation.
6. Appendix F: implementation details of the experiments.
7. Appendix G: additional simulation results and sensitivity analysis.
8. Appendix H: discussion alternative methods.
A THEORETICAL RESULTS: PROOF AND DISCUSSION
A. VARIATIONAL. FORMULATION OF ODE
Before we prove Proposition 1, we need the following lemma which a particular formulation the Fundamental lemma of calculus of variations Elsgoic (1961).
Lemma (Fundamental lemma of calculus of variations). Let h be continuous function closed interval [0, T h is equal everywhere if and only if S h(t)g(t) dt 0 for all C¹ [O T such that g()) (T) 0.
Proof. The forward direction is trivial let focus the converse.
Assume for contradiction that there is a point tu E [0, T such that h(tp) * 0. Without loss of generality, assume that h(tu) > 0. From the continuity of h, know that there is actually point 4 E (), such that 1(4)) > 0. Continuity of à implies that there exists & > 0 and small neighbourhood
8, + 6) C (0, T) such that h(() VI E (41 5,11 + 5)-
We can now define a non-negative continuous function g [0, 7] - [0, +00) such that g(t) = outside the neighbourhood (t1 8,61 8) and g(t) for some interval inside (t1 6, ti + 8). See Figure for an illustrative example.
4--5
1 4 6
Figure 6: A visual illustration of the proof of Fundamental lemma of calculus of variations
Then
S 1(D0(1) di
11+8 J h(1)g(t)d
/
(10)
That contradicts %f' h())(()) dt 0 for all g E (0, such that g(0) (TT)
Proof of Proposition I. Observe that
at(1) [z((((() + fz(x(1) 4,(1)
(11)
14

---

Published as a conference paper at ICLR 2022
Now using Fundamental lemma of calculus of variations, get that the following two statements are equivalent
f(((6)) - -,(0) |0,
(12)
(((0))))(())(()) =1l0))- c (o, g(0) g(T)
By linearity and integration by parts, we get
I/" (ff(((t)) (((a(t)))+((t)))9()) (t)))g(t)dt== L'" f)(a(t))((t)dt
L'" 2z(t)g(t)dt
l [)(x(t))((tddt L [)(m(())(()d) + L/ x,(()(()dl
+ L.
- -2((T)(()++
(13)
where in the last equality we used the fact that g(0) = g(T)) = 0. That means that
[ (,(((00) -(()(())((()))())
C,(fs,z,9) =
(14)
This proves that if and only if
2,(t) f)(x(1)) E [I....,.]
(15)
C,(,, 0, Vj {I.....,J} J} Ag Co g(0) 9(T) 0
(16)
Discussion. As we show in the proof above, the key idea of the variational formulation is bypass the derivative of the trajectory æ (t) with the derivative of a testing function j(4) - the swapping of derivatives is achieved by integration by parts (Equations 13). An important consequence is that we no longer need to know or estimate (t)t) (which is unobserved and challenging to estimate), but only need to know (t)) for a given testing function g()). Indeed, if gt) has an analytical derivative, we immediately have access to g()).
Note that the variational formulation builds on the exact condition C, = 0. It does not consider the situation where Cj o 0, no matter how close C₁ is to O. Hence, the formulation itself is inadequate to guide the search of closed-form ODEs. For instance, suppose the C, for two candidate functions fi and f2 are 0.01 and 1000 respectively. The variational formulation does not imply fi fits the trajectory æ better than fz; rather, it implies that fi and f2 are equally bad as C * 0 for both. In order to bridge this gap, we need to establish notion of distance between functions and link the functional Cj to this distance they are addressed in Appendix A.2 and A.3 respectively.
The formulation above is tailored for ODEs, which the focus of the current work. However, the variational formulation can be generalized to other types of differential equations, e.g partial differen- tial equations (PDE) and delay differential equations (DDE). In fact, the variational formulation is the theoretical foundation of the Finite Element Method (FEM) the industrial standard for solving PDEs in complex engineering systems and fluid dynamics (Zienkiewicz et al., 1977; Gokhale, 2008). It has also been used to estimate the parameters of known ODE DDE (Brunel et al., 2014). We envision that future works may be built on the general variational formulation to discover closed-form PDE DDE from data.
A.2 DISTANCE FUNCTION
In this section, we justify that df, f*) is appropriate metric for our problem.
15

---

Published as a conference paper at ICLR 2022
There two main ways of comparing functions. We either compare their symbolic representa- tions or measure the differences in values they take for set of arguments. As we want this metric to be a limit of some objective computable from the data, focus the latter.
Commonly used distance functions are metrics induced by L norms. LP norm is defined in the following way:
|| (f (0(0)))))))))
(17)
In practice, L or L norm is used. The induced metric is then defined as the norm of the difference between the functions, namely d(f.g) Ilf -l1pe
Now the question remains which difference we should take the norm of. We have three reasonable choices:
1. f f
2. Ø- where is trajectory prescribed by f computed from estimated initial condition
3. f f composed with some other function
The first choice poses the challenge of choosing the range (a, bj over which we integrate f - f _ In general we do not make any assumptions on the domain of the functions. also seems unreasonable to compare the functions far from trajectory, where we have no information about their behaviour.
The second choice is similar to what is done in the Neural ODE approach. As we mentioned before, this involves solving the Initial Value Problem (IVP). Even if the IVP has unique solution, it might be very sensitive to the initial condition, making it unfit for chaotic systems. We verify this claim in the experiments (Figure 4)
The third proposition can solve the main issue of the first option. we compose f- with the ground truth trajectory, æ; not only have well-defined interval over which to integrate (JO, T)) but also we focus our attention on values of f where we have chance to recover it, and where matters for the most.
That is why we decided to chose dz. = ||0 - /") o a as our distance function. Moreover, it is based on familiar metric, it inherits many important properties, such
1. de is a pseudo-metric', in particular, f f = dx(f.f") =
2. (f) 0 f= f x (100, TT), the image of æ
3. For a sequence f/ such that (f. f*) 0 it holds that fi x - f x in the mean'
The third property tells us that if we keep minimising d=()"") we can expect f getting closer to f* on x(0,T)) the image of æ.
A.3 PROOF OF THE THEOREM
Proof of Theorem I. Let us consider a Hilbert space L 2'J where the inner product denoted by ( and defined as:
(a, (a.b) 6) = J' a(t)b(t)dt
(18)
where and are real functions which are square-integrable on the interval [O, The L2 norm can be written as:
||a||2 = (a, a)
(19)
It can be considered to be metric through metric identification. We would need declare functions heing equal they equal (=00, T])
'This can be strengthened to imply uniform convergence if the sequence fi is equicontinouus (Ford & Pennline, 2007). This can be satisfied if, for instance, they all share the same Lipschitz constant. In particular, this is true if their derivatives bounded by the same constant.
16

---

Published as a conference paper at ICLR 2022
Using this notation rewrite Cy @k. go) as: Cj(f, kk. 9s) (f 3x 9n) (Fkj.
(20)
We assume that Ty converges in mean. For vector-valued functions that means that
[x((t) xt1) dt as k +
(21)
Note that || 1|2 in this equation denotes the usual Euclidean norm the arguments are vectors not functions. This is equivalent to saying that all components of w converge to the corresponding components of Thus
k. 2y 2|112 + 0 as k
(22)
As the inner product with one argument fixed can be seen as continuous operator, we get:
(azj. 96) + (r)z 96) as h + oc We assume that fj is X-Lipschitz continuous, That means that
(23)
[| o 3t fo (())x(() (a((1))) di
< (1) x=(1] de
(24)
X2 (" (|āx(1) =(1) dt + 0 as & - 00
1"
This means that
(f Cx. go) + (f o x. go) as k + x Combining equations (23) and (25) with expression (20), conclude that
(25)
lim C2 gs) 2 ((5 x. gs) (x), k0000
(26)
Now, using integration by parts we get
(rjs (5. S' az(O),())dt ())00(()))= S/² (t(()g + xr, (T)gg) - I. (0)g0)=)
(27)
We use the assumption that (0) = (T) = V to obtain
(x Oz) = - 3 9-)
(28)
We substitute it back into Equation 26 and use the lincarity of the inner product to get
lim Cf (f 2z)
(29)
Now we can use Parseval's identity and the fact that the set 91, 92. forms Hilbert basis for L, 7] to infer that:
lim 500
1ffox-1,.90 1lf 20112
(30)
Substituting Equation 29 into Equation 30 and using the fact that the inner limit does not depend on we obtain the statement in the theorem
lim lim Σ C, Ilf 2,11 II(f 1r)
(31)
S-
17

---

Published as a conference paper at ICLR 2022
A.4 DIVERGENT DERIVATIVES
The main advantage of our algorithm lies in estimating the trajectory instead of estimating its time derivative. We show how the derivatives of uniformly convergent sequence of functions may fail to converge, and we illustrate how it impacts objectives based on estimating the derivative.
We highlight the following mathematical fact:
Remark. Uniform convergence of functions does not guarantee convergence of the derivatives
This is in contrast to definite integrals where such implication holds. Intuitively, this follows from the fact that the derivative is a very local property of the function. To see an example of such a sequence, consider a sequence of functions (hx) where h2 [0, + given by
h₂(t) sin(kt)
(32)
k
It can be shown that (hx) converges to h(t) 0 uniformly. However, the derivative of h2, given by
hig (t) cos(k:),
(33)
fails to converge at any point.
We use this observation to show that the square of [| - f where ăy is the derivative calculated from the recovered trajectory é, does not converge to dz (fe, fz)² when converges to the ground truth trajectory æ.
Consider the following ODE (J 1,T
=(() xx) 4 E [0.1]
(34)
so fx)x = x. We choose x(t) e as the ground truth trajectory from which we sample. Let (&x)
be a sequence of trajectories given by d = + sin(kt) uniformly. Then for any candidate function f
It can be shown that (fx) converges to o
110x- f 2x S (BaCC) (fa(t)- ((t)(t)))
(35)
= L/' (d''==cosii cos((kt) - f ((t sin(kt) (d==in(h0)) dt
k
If [1 - fo a l were to converge to dz f", then in particular l| - f 2 should converge to O. However, we show this is not what happens.
11.b fr o x = 1/-
c² + cos(kt) frr (c' sin(kt)
siin(ho)
k )
dt
e' + cos(kt) (ee sin(kit) )) dt
(36)
(cos(k))
sin(kt)
dl.
J(¹ cos² (kil) +
ssin(k)) 2
ssn((k))
2cos(kt)
dl
k
As sin(kt)
and (kt)(kt)
sin(k()
converge uniformly function, we use the linearity
of integration to obtain
[ coss (kttddt = lim
sin(2k)
1
lim | - f 1 lim k
#
(37)
k-
k0000
4k
2
2
This example shows that objectives that estimate the derivative highly dependent on the type of estimation algorithm used. In contrast, Theorem tells us that the only property of the estimation algorithm need is that the trajectories converge the ground truth in 7,2 norm.
18

---

Published as a conference paper at ICLR 2022
Algorithm D-CODE.
Input: Dataset D (t) & NT
Input: Smoothing algorithm S, numerical integration algorithm Z, optimization algorithm O
Input: Testing functions Js. S S, initial guesses fςo, j S
T. = ((y) ....))....)))) for j E [1 do
(Preprocessing step: smooth input trajectories)
f f=
(Initialize the juh ODE)
Converge False
while Not Converge do
obj = >i E C2 9a)
{Calculate objective (Eq. 5).} (Run optimization step.)
J. Convergee = ((bb {z)
end while
end for Output: The discovered ODEs F,, j S
B D-CODE FRAMEWORK
B. PSEUDOCODE
The Pseudocode is presented in Algorithm 1.
B.2 DISCOVER HIGHER ORDER AND NONAUTONOMOUS ODES
In the main text, focused the system of first order autonomous ODEs (Equation 1). This because higher order or nonautonomous ODEs can be transformed into first order autonomous systems by including more variables. For instance, consider the damped, driven harmonic oscillator, a model system for classical mechanics (Simmons, 1972):
5 + 1 + #1 cos(1)
The ODE second-order as it involves the acceleration term 31 and it is nonautonomous because of the time-dependent term cos(t). If we introduce two new variables: I2 d1 and 23 = t, we obtain the following first order autonomous system:
Φ 42
22 -2 21 + cos(Ts)
Since the first and the third equations implied by the definition of 22 and 33, they do not need to be discovered from data. One can apply the discovery algorithm to uncover the second equation from the measurements of 11, 22 and Z3.
C DISCUSSION ON DATASET CURATION
The success of ODE discovery depends the availability of well-curated dataset. However, curating dataset is highly nontrivial task, and it is well beyond the scope of the current work. In fact, the curation of dataset sometimes an important scientific contribution in itself (Blasius et al., 2020). Here, we provide brief discussion on this topic.
Selection of variables. The researcher needs to decide which variables to include in the dataset this defines the scope of problem. In general, variable selection depends the scientific problem being addressed and the background knowledge about the interactions between variables. For instance, one may use the known chemical reaction network to choose small set of compounds to include in the data. When such knowledge is unavailable, the researcher may need to perform pilot studies to find out the interaction network (otherwise they may face challenging high- dineennnnnii problem). We note that there a well-developed literature on uncovering the interaction network of a high-dimensional dynamical systems (Mangan al., 2016; Qian al., 2020).
19

---

Published as a conference paper at ICLR 2022
Measurement settings. In practice, the measurement settings often affected by practical con- straints, e.g. the precision of measurement device, the measurement cost, the time limit, etc. However, the researcher should ensure that the measurement horizon [0. T] covers the pattern of scientific interest. To illustrate this point, consider an application in climate modeling. we are interested in studying the cyclic behavior of the system (e.g. seasonality), the measurement should cover least one entire period (one year), On the other hand, if we are interested in the system's 's response to an external stimuli (e.g. short-term weather change due volcanic activity), the measurement should cover the period before and after the stimuli (volcano eruption).
Iterative nature of discovery. Scientific discovery is by no means a single-step monolithic process. Rather, it is an iterative process involving many components: hypothesis generation, data curation, modeling, validation, etc. D-CODE focuses on the modeling component, but it can facilitate the entire loop of discovery because it distills an interpretable close-form ODE, For instance, if the discovered ODE converges to wrong stationary point (as indicated by domain knowledge), the rescarcher may curate an additional dataset with longer time horizon and re-train D-CODE with the additional information.
D UTILITIES OF THE DISCOVERED CLOSED-FORM ODE
In the main text, we focus on the problem of distilling closed-form ODE from data. Here we illustrate how the learned closed-form ODE can facilitate various downstream tasks and applications.
Prediction. The learned closed-form ODE can be used to make predictions about the future states of a dynamical system. Formally speaking, given sequence of past observations y(t)), (t2))..... y (tx).. we would like to predict x(t") for a future time t* > tg. This prediction problem is very well studied and has many proven algorithms (e.g. the family of Kalman filters (Welch al., 1995) and the Bayesian methods (Doucet et al., 2000)).
Simulation and synthetic data. The learned closed-form ODE can also be used as a simulation model to generate synthetic datasets. The user can specify different initial conditions a=0) and different external signals (if applicable) to generate a range trajectories. The simulated trajectories may be used in downstream tasks, c.g. computing the statistical properties.
Planning and control. some cases, the system contains variables that can be intervened on (e.g. external forces), and we would like to manipulate these variables to maximize certain utility function. This is the canonical problem setting for Control Theory, where decades research have been devoted to deriving the optimal policy from the closed-form governing equations and the utility function (Mchrmann, 1991; Haddad & Chellaboina, 2011).
Understanding the properties. Often we are interested in understanding the high-level properties of the dynamical system. Examples include asymptotic stability, fixed points, periodicity, bifurcation, and The closed-form ODE would facilitate the study of these properties because the analytical and the computational tools developed in the dynamical system literature are both applicable (Browder, 1983; Lyapunov, 1992; Chow & Hale, 2012).
E BACKGROUND ON NUMERICAL DIFFERENTIATION
The task of estimating time derivative from discretely sampled and noisy data commonly referred to as numerical differentiation. Finite difference is a family of the most basic numerical differentiation methods (Bickley, 1941). For example, the central difference estimator calculates
y(1 At) v(1 At)
2At
The main disadvantage of finite difference is that is highly sensitive to measurement noise (Char- trand, 2011). 1).
common way to address this shortcoming is to estimate the derivative on the smoothed trajectory. For example, opline-smoothed differentiation smooths the observed trajectories with cubic spline and estimates i from the fitted spline functions (Shen et al., 1998; Zhou & Wolfe, 2000). (This method is used in the benchmarks SR-S). However, this family of methods is usually based on assumptions on
20

---

Published as a conference paper at ICLR 2022
Table 3: The detailed settings for each simulation: noise level øn, time step size At, number of trajectories N, time horizon T, range of initial conditions (a, b), numeric constants 0. Bold values are defaults.
System
OR
At
N
a, bj
Gompertz
1.3]
0.1, 0.2, 0.5,
20, 50,
4
00 0.0.01]
01
1.5
Logistic
1.3
0.1, 0.2, 0.5,
5, 20, 50, 100
10
0,0.1
01 = 1,02 = 0.5
Glycolytic
0.01, 0.1,
0.1. 0.2, 0.5,
5, 20, 20, 100
15
0,0.1
6: 0.75, 82 - 0.1
Lorenz
00.1,.03)
0.01, 0.02, 0.04, 0.1
2,, 50, 100
10
[0, 10]
10,28 28, 8/3
the smoothness of the trajectory and its (higher-order) derivatives (Zhou & Wolfe, 2000); they encode bias towards smooth and slow-changing derivatives. Such bias may the hamper the estimation of dynamical systems with abrupt changes (Chartrand, 2011).
To address this challenge, Chartrand (2011) proposes to explicit regularization the finite difference estimator. As we pointed out in Appendix A, the standard 1.2 regularization is not suitable for derivative estimation. The authors propose developed in computer vision to obtain sharp denoised images (Strong & Chan, 2003). This type of regularization allows rapid (or even discountinous) change æ. The method has been adopted in the two-step sparse regression literature (Brunton et al., 2016), and it is used in the benchmark SR-T.
total variation regularization, method originally
IMPLEMENTATION DETAILS OF THE EXPERIMENTS
F.1 DETAILED SETTINGS FOR EACH EXPERIMENT
The detailed settings for each experiment in Section 5 is shown in Table 3. The time horizon and the range of initial conditions (a, example, we observe that the Gompertz model reaches saturation and converges to the constant x = 1 before = 4 (hence set the time horizon T In practical applications, the time horizon is decided by the problem scope and the availability of experimental or observational data (see the discussion on dataset curation in Appendix C).
chosen based on the scale and properties of the system. For
F.2 HYPER-PARAMETEE SETTINGS
Smoothing algorithm. In the experiments, D-CODE uses Gaussian process in the pre-processing step. We use the radial basis kernel with the hyper-parameters tuned based on marginal likeliielihood (Bernardo et al., 1998).
Genetic programming. The hyperparameters of genetic programming is decided based on a pilot study. In the pilot study, we consider noiseless and frequently sampled Gompertz growth model (a = 0, Al = 0.01). We use randomized search on Two-step Symbolic Regression to find the hyperparamettrs that achieves the best RMSE loss (Table 1). The same hyperparameters by all algorithms on all experiments. The values are listed below. We refer the reader to the documentation of gplearn library for detailed explanation of these hyperparameters (Stephens, 2021).
used
1. population size: 15000
2. tournament size: 20
3. p crossover: 0.6903
4. p subtree mutation: 0.133
5. p hoist mutation: 0.0361
6. P point mutation: 0.0905
7. generations: 20
8. parsimony coefficient: 0.01
21

---

Published as a conference paper at ICLR 2022
The set of allowed mathematical operations is
{+, - x, log, cos, power)
D-CODE and all benchmarks the same set of mathematical operations for fairness of comparison. To speed up computation, we run a pilot study for each dynamical system in the noiseless setting described above to select a subset of these operations for the main experiments.
Numerical differentiation. We use the python package derivative veriiiiiin ive for numerical differentiation (Quade Goldschmidt, 2020). The library contain implementations of the methods discussed in Appendix E. Currently. there lacks principled way to tune the hyperparameters of numerical differentiation because the time derivative is not observed. Existing works often rely on heuristics (Chartrand, 2011) 1). Here, we tune the hyperparameters based on the true derivative using grid search over (0.001, 0.01, 0.1, 1, 10, 100, 1000). Note that this procedure may allow the time derivative to be estimated more accurately than what is possible (because the true derivative is not available in real data). Hence, it will give an advantage to SR-T and SR-S because D-CODE does not use the true derivative in any way. SR-G uses the same Gaussian process D-CODE, which is tuned by the marginal log-likeling ligllikeiihood on the observed trajectory.
Correctness of the functional form. As discussed in Section 2, the functional form of closed- form function f is the expression with all the numeric constants replaced by placeholders 01, Oz.... In practice, the uncovered function f may contain constants that are close to O. For instance f(x) = x + 0.001x is close to (x) x when is small. Hence, we apply thresholding on the estimated constants and set all constants between [-0.05, 0.05] as 0 (only for logistic model and Glycolytic oscillator, where we found thresholding be helpful). The thresholding is applied the results of all methods to evaluate success probability of recovering the functional form and it is not used when calculating the RMSE of parameter estimation or the distance dz We use Python package sympy to decide whether the discovered functional form is correct (Meurer al., 2017). i.c. success when simplify(f
returns 0.
F.3 CLEANING TUMOR GROWTH DATASET
The original tumor growth dataset (Section 6) contains heterogeneous growth patterns and outliers. However, the dataset does not include any patient covariate that could help to model or explain the heterogeneity, i.e. some variables that would drive the tumor growth dynamics are not recorded in the data, e.g. age and comorbidity (see the discussion in Section 7). To address this issue, we performed k-means clustering on the observed trajectories (with 10 clusters), such that the growth pattern within each cluster is more or less homogeneous. We selected the largest cluster with 310 trajectories for the analysis.
G ADDITIONAL RESULTS
G.1 ADDITIONAL EVALUATION METRIC: dee
Here we report the distance function de (f) between the true and estimated ODEs. This metric tracks how well ff approximates f*. The rationale for reporting this metric is that it may be possible that f has the wrong functional form but it still well approximates f". Hence, we need the distance dn in addition to the discovery probability.
The results are presented in Figure 7 and 8. We observe that the patterns of d₂ and discovery probability (reported in the main text Section 5) are very similar - the performance gap between the methods is minor in the low-noise frequeent-sampling scenarios, but D-CODE significantly out- performs the benchmarks in the more challenging measurement settings. The agreement between de and discovery probability suggests that when the correct functional form is not recovered, the estimated / tends to perform poorly in approximating the true function as well.
22

---

Published as a conference paper at ICLR 2022
Gompertz Model
Varying noise level
Varying step size At
Varying sample size
Objective dx Objective d,
0.10
0.05
0.05
10
10
10
Generalized Logistic Model
0.10
Methods
0.02
0.05
0.02
SR-T
0.00
SR-S
10
10
10º
2.0
100
SR-G
OR
At
Figure The metric dæ (Eq. 6) two growth models under different settings (smaller better). The three columns correspond to various noise levels d, step sizes At, and numbers of trajectories N. The two rows correspond to the Gompertz model and the generalized logistic model respectively.
Objective d, for kx())
Objective d, for Xy(f)
Objective d, for Xy(()
20
10
364
10
5R-5
56-5
0.10
0.15
0.20
0.25
0.30
0.10
0.15
0.20
0.25
0.30
0.10
0.15
0.20
0.25
0.30
On
0a
Figure 8: The metric du (Bq. 6) for the chaotic Lorenz system under different noise levels (smaller better)..
Discover Prob.
Objective d,
1.00
DOCDD
0.75
0.50
0.00
0.0
0.2
0.6
0.8
0.0
0.2
0.4
0.6
0.8
Figure 9: Simulation results for the dynamical system in Eq 38 under various noise levels. Left: the success probability (higher better) and Right: distance dz (lower better).
G.2 ADDITIONAL DYNAMICAL SYSTEMS
Here consider additional dynamical system, which is extension the quadratic integrate and fire (QIF) model in Neuroscience (Brunel Latham, 2003). The governing equation is
x2
(38)
t t
We highlight that this ODE is fime-dependen and it involves a fraction, which makes it impossible to be symbolically represented in a linear form f(x) = Ex 02 (as required by the two-step sparse regression method). similar example was given in Appendix B-2 Brunton et al. (2016) to
illustrate this failure mode of two-step sparse regression method.
The simulation results are shown Figure 9. We observe the same pattern as the other dynamical
systems considered in the main text.
G.3 SENSITIVITY ANALYSIS FOR TESTING FUNCTION
Here we study the sensitivity of D-CODE to the choice and number of testing functions. We consider two families of testing function: sine functions and cubic spline functions, both of which satisfy the
23

---

Published as a conference paper at ICLR 2022
GompertZODE sine
GompertZODE cubic
100
100
Lorenz sine
Lorenz cubic
1.0
0.5
0.5
0.0
0.0
40 Number of basis
100
Number of basis
Figure 10: The success probability with different families and different numbers of testing functions y (sine and cubic spline).
Computation time Seconds
1750
LHCOOL
S6-11
1250
750
250
Competix
Logiskic:
Glycu,
Lorenz
Figure 11: Comparison average training time.
conditions Section 4.2. We refer the readers to De Boor (1978) for the exact recursive definition of cubic spline functions.
Figure 10 reports the discovery probability for Gompertz model (first row) and Lorenz system (second row). We first observe that using few testing functions (i.e. < 40) leads performance loss. This agrees with Theorem because the objective function Eq. 5 cannot well- rell-approximate the distance dæ (Eq. with few testing functions. We also observe that the performance reasonably stable when enough testing functions are included (e.g. 40 80 sine functions, and 80 - 100 cubic splines). Finally, we see a performance drop for using too many (e-g. 100) sine functions. After investigation, we find that the drop is due to the numerical eтror in the integration procedure - sine functions with very high frequency are oscillating too rapidly for accurate numerical integration. One may resolve this issue by using smaller integration step size more advanced numerical integration method.
For practical applications, the user may experiment with different numbers testing functions and examine the discovered ODEs. A key advantage of discovering close-form ODEs is that it allows the user to examine and refine the ODEs in an iterative way.
G.4 TRAINING TIME
Figure 11 compares the average training time of D-CODE and two-step symbolic regression (SR-T). We expect D-CODE to have longer training time due to the numerical integration step in evaluating the functional C (Eq. 3). However, even for the longest-running task, D-CODE finishes training within 30 minutes average. This can still be very significant speed up compared with manual discovery by human experts. Experiments are run on a computer with Intel Xeon E3- 12xx v2 CPU (16 cores) and 60 GB memory.
G.5 TABLE OF DISCOVERED ODES
The top-3 most common equations discovered by different methods for each dynamical system (among the 100 simulation runs) are listed in the table the next page.
24

---

Published as conference paper at ICLR 2022
Table 4: The top-3 most common equations discovered by different methods for each dynamical system (among the 100 simulation runs). "Freq" denotes the chance of discovering the equation. We make the following observations: (1) SR-T and SR-S produce over simplisii equations for the growth models (constant zero), possibly because the estimated derivative is too noisy to learn useful signal. (2) SR-T and SR-S may consistently miss variable (c.g. the zu term Glycolytic 1) or miss an interaction (c.g. the Zu term in Lorenz 2). (3) SR-T and SR-S may include unnecessrry terms in the equation (c.g. Lorenz 3). The measurements settings are the default ones listed in Table 3.
D-CODE
SR-T
ST-S Discovered Equation
System
Discovered Equation f
Freq.
Discovered Equation f
Freq.
Freq.
Fr log((227) (f2z))
0.63
0.67
Gompertz
01
0g x² log (a)) log(x**)
0.04
- 01 (z + log(z))
0.03
0.04
log(log( (0: + x)
01 02
0.96
0.98
Logistic
01
0,02
0
0.02
0,
2
0.02
0130 + 02 2p2l
0.38
81 xu
0.5
0.1
35
Glycolytic
01 0221 + Og- Foo
0.03
-0,000 - 0₂ F02
0.22
-0 =01 + 02
0.08
81 z0 (02+ + x
0.03
01002 02200
0.04
-0, (02 + 31 + Os
0.07
ra zozî 31
0.78
0130 + Fuz 21
0.75
f22f 00 + Z021
0.09
Glycolytic 2
0100 + 02 + xoti
0.13
ro(0) + xx) 21
0.09
+ 0200 0zz1
0.08
0150 O2æi (Ogaro =ox1)
0.01
01 0300 04w
0.02
02a7
0s + 2022
m1?
0.06
8130 8221
+ fgw
0.63
Lorenz
8120
0222 + 0s
0.37
Fu Tu22 21
0.33
-100 0231 0g22 0g
0.99
zu x2))
0.73
Lorenz
[0 x2)
20 31
0.28
01r0( 02322 + 03) 0421 20
0.01
0,20 J072
0.2
00 x2)
F1
0.28
zo (f) x2) 32
0.03
0112 + To
0.8
0102
02 + 00021 + x0
0.46
0,22 + x001
0,51
I Lorenz 3
01 + 0g + T0
0.13
-0177
- Jg + 310
0.11
0,72 + To
+ 300
0.11
0132 02 + 7071
0.08
01 0222
0.1
0,72 + T071
0.09

---

Published as a conference paper at ICLR 2022
Gompert2ODE:
GompertZODE: Success Prob-0.98
GompertzDDE:
Succ2ss
Xit
10
8D
05
0.5
0.5
0.0
10.0
GompertzODE: de 0.7. At 0.1 Succass Prab-D.99
CompErZODE de - 0.7. At 0.2 Success Prob-0.98
GompertZODE: dv D.7. A0-10 Success Prob-0.97
1.5
1.5
1.5
ytD xtD
1.0
1.0
1.0
0.5
0.0
0.0
Visualization of the observations yt). the true trajectory 2(t) and the smoothed trajectory (()). First row: increasing noise level ø- Second row: increasing sampling sparsity At. The success probability of D-CODE is shown in the headings.
Figure 12:
G.6 VISUALIZATION OF SMOOTHEE AND TRUE TRAJECTORIES
In this section we study whether D-CODE can still discover the ODE when the smoothed trajectory (tt) is significantly different from the true trajectory (t):: for instance, when the signal-to-noiss ratio is low or the observations are sparse.
In Section 5 we provided quantitative evaluation under scenarios with different noise levels (a)) and sampling sparsity (At).5 For further illustration, here we visualize #(t) and æxt)) together with the observations ( in Figure 12. In the most favorable setting (top left), #tt) is indistinguishable from att), and D-CODE achieves very high success probability as expected. As the observation setting becomes more challenging (middle column), (()) may recover the overall trend of (() but it may also contain some estimation artifacts (e.g. the unnecessary zigzag pattern). However, D-CODE is still able to recover the true ODE with very high probability. In the most challenging settings (right column), (t) differs significantly from a=(() for extended periods of time (top right /. 2, bottom right / < and > 4). Yet, D-CODE still makes successful discovery at times.
H DISCUSSION ON ALTERNATIVE METHODS
In this section, diseuss the alternative methods that, in theory, could also be applied discovering closed form ODEs from data. However, as we will show, these methods suffer from various drawbacks that may severely limit their practical utility.
H.1 DOUBLE OPTIMIZATION METHOD
The double optimization method (DO) was first proposed in the literature of ODE parametric estimation, where the functional form of / known but it contains unknown parameters 0 (see Ramsay et al. (2007) Section 1.3 for a review). DO solves the following optimization problem (we consider the one-dimensional case to simplify notations):
N
Σ (1=(1));;1)) =. subject to a.(() = (0) + S f(ix(s))ds (39)
f = arg min
min (0)
1=1 LET
This is nested optimization problem where the outer part searches for the unknown function / and the inner part searches for the unknown initial conditions z;(0) given the current guess of f. Note
Note that, by definition, the signal-to-noise ratio is the inverse of dR. i.c. 1/08- 'DO also referred the nonlinear least square method (NI.S) some literature.
26

---

However, the choice of prediction horizon 7 non-trivial and plays an important role in the training. For chaotic systems, the user should choose smaller 7* because the system is unpredictable in the long term. On the other hand, making T too small may lead to a myopic model that fails to capture long-term patterns (Schmidt et al., 2020; Kantz & Schreiber, 2004).² In practice, the prediction horizon 7 is hyper-parameter that may be difficult to decide (note that the Lyapunov time of the specify the prediction horizon Published as a conference paper at ICLR 2022
Gompertz ODE
Logistic ODE
Runnira time (Seconda))
2 uanaaaa time (Socco ada
+
Mean 3.09
Meen 0.07
Mear
Mean 0:15
DO
D.CODE
DC
D-CODE
seconds) evaluate candidate function f. The data are generated from Gompertz and Logistic ODEs respectively. The results are calculated from 100 independent runs. The summary statistics (mean, median and quartiles) are shown in the boxplot. DO is 30 50 times slower than D-CODE.
Figure 13: The computation
that the inner and outer parts are coupled because the optimal (0) depends on f. e.g. the best-fitting exponential curve would have different initial condition from the best linear fit.
A well-known shortcoming of DO is its computational complexity (Gugushvili & Klaassen, 2012).
The inner optimization over (0) is typically non-convex and may take many steps converge. To further aggravate the problem, each inner optimization step requires solving the IVP numerically.
The computational challenge becomes even more pronounced our setting, where the functional form of / is also unknown. Genetic programming needs to evaluate many candidate functions / but DO requires solving the inner optimization for each evaluation. In contrast, D-CODE (Equation 5) does not involve any nested optimization.
To illustrate this point, we have implemented DO and measured its computation time to evaluate candidate f. The results shown in Figure 13. We found that DO is typically 30 50 times slower than D-CODE. Its computational cost is too high to be considered as feasible method for discovering closed-form ODEs.
H.2 NEURAL ODE WITH POST-IIOC DISTILLATION
Although Neural ODE (NODE) does not directly produce closed-form ODEs, one may attempt to distill a closed-form equation from the trained network. The method takes two steps: (1) train a NODE fN using observed trajectories and query the trained NODE to produce dataset (t), JN ( (2) apply symbolic regression on this dataset. This approach can also be viewed as two-step regression method (Section 3), where the derivatives are estimated by NODE,
The standard approach to train a NODE for time series modeling is to minimize the error over prediction horizon T [t1, 12, ... T**] (Rubanova al., 2019), i.c. minimize
C
(1 - (1)
(40)
EET
where é; (t) is obtained by solving the IVP of the NODE.
true system is unknown). In contrast, D-CODE does not require T""
Furthermore, alluded earlier, this approach an instance of two-step regression, where the derivative is estimated by NODE instead of numerical differentiation. Although NODE is known to
'The limitations of minimizing short-term prediction loss for sequence modeling is well-recognized in machine learning; for example, the eequence-to-eqqenne models tackle this issue by minimizing the loss over sequence of predictions rather than the step ahead prediction.
27

---

Published as a conference paper at ICLR 2022
Gompertz Model
Lorenz X1(8)
Lorenz Xx(8)
Lorenz Xy(t)
1.00
1.00
1.00 00
1.00
Success ProD.
Methods
0.75
0.75
0.75
0.75
D-CODE
NODE-L
0.50
0.50
0.50
0.50
NODE-S S
0.25
0.00
0.00
0.00
0.00
0.0
0.5
1.0
0.1
0.2
0.3
0.1
0.2
0.3
0.1
0.2
0.3
da
Ga
On
On
Figure 14: The probability of successfully recovering the unknown ODEs. Two versions of NODE are compared with D-CODE on the Gompertz ODE (first panel) and the Lorenz system (second to fourth panels), In the rightmost panel, the results of NODE-L and NODE-S overlap and only one line displayed.
Gompertz Model
Lorenz (t)
Lorenz x2(()
Lorenz xx(()
80
Methods
Objective d.
100
100
D-CODE
NODE-L
NODE-S
50-
01
0
0.0
0.5
1.0
0.1
0.2
0.3
0.1
0.2
0.3
0.1
0.2
0.3
Noise
Noise
Noise
Noise de
Figure 15: The distance da between the true and the learned ODE.
be a universal function approximator, its estimation efficiency (e.g. vn consistency), finite sample performance, and robustness to noise are still largely under investigation. Hence, we should remind the readers that NODE is not guaranteed to outperform numerical differentiation in estimating derivatives.
To illustrate these points, we evaluate the performance of NODE with post-hoc distillation. We consider two versions of NODE: the NODE-L is trained on long term prediction (T" = 25) and NODE-S trained on short term prediction (T 0.04). The evaluation is performed on the Gompertz ODE (non-chaotie, one dimensional) and the Lorenz system (chaotic, three dimensional).
Figure 14 shows the success probability (bigger better) and Figure shows the distance de (Equation 6, smaller better). For the non-chaotic Gompertz ODE, both NODE-L and NODE-S only able to recover the correct functional form when the noise is very close to zero (Ø 0.01). Note that this behavior is very similar to other two-step regression methods (SR-T, SR-S and SR-G, Figure 3). We also observe that the NODE-L achieves consistently better da than NODE-S. This confirms that training over longer time horizon leads better estimate of the ODE when the system is non-chaotic.
Turning to the chaotic Lorenz system, we observe that NODE-L always fails to discover the correct ODE and its de higher than the other two methods. This verifies our expectation that training with long-term prediction is not feasible for chaotic systems. As before, NODE-S performs well only when the noise level is relatively small. We also observe that its d, monotonically increases with the noise level. This suggests that higher noise level significantly reduces NODE-S's ability to approximate the unknown ODE, limiting performance on noisy data.
28


Reviews:
[{"review_id": "Vk8787RNRuh", "reviewer": "Reviewer_79Ft", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The method is clearly described, and the potential benefits are clearly understandable (see summary). However, I find the baseline comparison somewhat incomplete.\n\nMy main questions are about the NeuralODE comparison. I do not agree that learning a chaotic dynamics directly with a NODE is not possible. The author argue that we need to estimate the initial condition what makes it impossible due to the chaotic instability. Only thing we actually need to train a NODE is to minimize the error in a shorter time horizon. Specifically, if we have samples more often than the Lyapunov time we should be able to learn the mapping. We integrate the dynamics forward inside a smaller than Lyapunov time window, and use this error to learn. \nPlease clarify this point, and train the NODE with a loss appropriate, and show that result on the rightmost plot of Figure 4. A simulation should stay on the attractor, of course the concrete trajectory is impossible to reproduce.\n\nThe valid argument still remains that NODE is not interpretable while closed form ODE is. I have no doubt on this point. However if NODE training is actually possible, it can be used similarly as GP and others in the two-step methods, and the gradient networks can be distilled to a closed form equation.\n", "summary_of_the_paper": "The paper addresses the problem of learning closed form ODEs from observational data, when the observation can be noisy and not frequent enough that instantaneous derivatives can be estimated with low enough error. The work uses a variational criterion on the solution of the ODE, circumventing the need to evaluate instantaneous derivatives. It uses existing symbolic regression techniques but instead of the regression loss, it minimizes a loss term quantifying the violation of the variational constraints. It is a work combining existing results in a novel way, but offers new perspective and there are synergies between the different parts applied.", "summary_of_the_review": "The work is derivative in methodology, but provide enough insight to be interesting. I have some questions on the Neural ODE comparison what would need some clarification. If this concern is addressed it would make it easier to support acceptance.\n", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works."}, "extra_text": {}}, {"review_id": "IAM0iK0rJFW", "reviewer": "Reviewer_EBvJ", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper is well-motivated and easy to follow. The proposed variational approach is sound for by-passing time-derivatives of noisy observations. It is also a model-agnostic framework, thus might be widely applicable regardless own’s specific setting.\n\nEmpirical validation on the proposed framework (D-CODE) shows that D-CODE is the right way to recover ODEs under some undesirable but unavoidable measurement errors.\n\nThe tumor growth experiment is very interesting because its closed-form ODE is indeed unknown. D-CODE can capture an interesting behavior (tumor volume increase due to drug resistance?) while the baseline cannot.\n\nI have not found any major drawbacks to this paper. For these reasons, I think this paper is worthy of publication for ICLR. The followings are some (relatively minor) questions:\n\n* Q1. Treating second-order non-autonomous systems as a set of first-order autonomous ODEs requires non-trivial effort since the number of variables increases to 2J + 1. Is there any suggestion about when one should consider such a situation (without prior knowledge of the target system)? It seems that the authors have found non-autonomous ODE for the tumor growth task, thus it will be nice if they can explain their modeling procedure briefly.\n\n* Q2. Is the search space for mathematical operations used in this paper enough for arbitrary systems? If not, and at the same time if the variational objective does not converge to 0 well for an unknown target system, then what would be the best option for the next behavior, e.g., adding some additional operations, or considering higher-order non-autonomous systems such as Q1? It may depend on the target problem, but it will be nice if the authors can discuss about it.\n\n* Q3. The authors state that neural ODEs should estimate the unknown initial condition (thus they are not very effective for chaotic systems), but it seems that one can rather directly use y(0) = x(0) + eps(0) as an initial condition for neural ODEs. Are neural ODEs still bad for such a case, though it might depend on the noise level and Lyapunov exponent of the target dynamics? The current evaluation of neural ODEs in Figure 4 seems to be extremely poor (and I felt the evaluation is slightly unfair), thus I wonder it is truly a limitation of neural ODEs or due to a bad estimation of the initial state.\n", "summary_of_the_paper": "This paper suggests a variational formulation for discovering ODE systems in a closed-form. The main advantage of such a formulation is that it can find the system based on the observed trajectories (and some analytical testing functions) rather than estimating (possibly noisy) derivatives of them. The variational objective is equal to the distance between estimated and true velocity vector fields on the manifold when the estimated trajectory converges to the true one in L2 sense and there are infinitely many testing functions (while a finite number of functions are okay empirically). The proposed method outperforms its counterparts for various ODE systems and tumor volume dataset (whose true dynamics is unknown).", "summary_of_the_review": "As I mentioned above, I think the paper is interesting and thus recommend the acceptance of this paper.\n\n***\nPost-rebuttal: I appreciate the authors’ clarifications on my questions. I'd like to keep my score.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "DK4TZf1rje4", "reviewer": "Reviewer_ZddY", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper is well written and addressing an important problem for scientific discovery and interpretability of ODE methods. \n\n\n- My main comment is that I think another baseline to investigate the importance of the proposed objective function would be the following : you optimize over an initial value and the function f and you minimize the reconstruction error. Is this a baseline you could compare against ? As it does not involve computing the derivative either, it could also be competitive and really assess the added value of the proposed objective function. I also believe this setup would fit Theorem 1. If this is not a good idea, can authors explain why ? As I don't see a clear advantage of the proposed method over the one I just described, it currently undermines my evaluation of the significance/novelty of this work. However, I'd be very happy to change my mind on this topic.\n\n- In Table 2 it seems that the performance differ massively between both dimensions (8) and (9). Yet, at a glance, their functional form look really similar. It would be nice to provide some extra information or insights in that regard in the text.\n\n- Because all these symbolic regression approaches make use of non-differentiable optimizers, the computational cost can be quite high. I would appreciate some discussion of the computational complexity of the different methods and to have an order of magnitude of the computation time.\n\n- It's also not exactly clear to me how the interpolation steps work in the multi-dimensional setup. In practice do you use single dimension GPs for each dimension or do you use multi-task GP to model the correlations ? \n\n- In Equation 3, the first integral should have $f_j$, as $f$ is not a scalar function.", "summary_of_the_paper": "This paper proposes a new methodology to infer symbolic ODE representation from observed time series. In contrast to previous methods, they bypass the inference of the time derivative and rather proceed by first estimating a continuous approximation of the trajectory and then optmizing a novel objective function. The first step can be performed with the interpolation method of choice (GP or splines). The second step can use any optimization scheme that does not require derivative (such as genetic algorithm). The authors then evaluate their approach. on a series of dynamical system and show improved performance compared to the baselines under consideration.", "summary_of_the_review": "I think this is a nice paper but to fully motivate the approach, it should compare against another obvious baseline. I would like the authors to motivate why the setup I described above is not considered in the experiments. ", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "DF2jinkPevJ", "reviewer": "Reviewer_vaG3", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Using symbolic regression to approximate closed-form ODEs is novel and very interesting. The results presented in the paper are compelling but I would want to see more comparisons to other methods, as well as results on more datasets. \n\nThe method section is very well described, with useful and precise details found in the appendix.\n\nThe weakness of this paper lies in the experiments and results section. \n\nAuthors have selected five dynamical systems governed by a wide range of nine closed-form ODEs. However, all these system have independent Gaussian noise (with different noise level). I assume that if the noise is properly handled by the smoother, this will not impact the model, but I would like to see experiments on datasets where the smoother does not necessary perform very well (low signal to noise ratio). In that case, will D-CODE be robust enough to recover correct parameters of the ODE? \nIf this is not the case, I don't really see the advantage of D-CODE in addition to \"converting\" an already inferred trajectory into an ODE for the sake of interpretability. \n\nSimilarly, how will it handle sparse observations? One advantage of methods inferring ODE/SDE parameters for dynamical systems are the ability to approximate trajectories in regions with few observations. Here, won't the method suffer from the smoother if not able to capture the trajectory correctly? \n\nTo understand better the advantages of D-CODE over other methods, I think that showing results on more than one real dataset (I understand the dataset presented here is very challenging), with comparisons to more than one method (for example, Neural ODE, discussed here in the simulated data experiments). \n\nAlso, is D-CODE restricted to one dimensional observations? (would it be possible to perform symbolic regression for higher dimensional latent variables?) If it is possible, then an extension of this framework to higher-dimensional observations and latent processes would be very compelling and allow for comparisons with many models from the Dynamical systems literature that infer latent ODEs.\n\nMinor comment: the figures (especially Figure 2) could be made bigger and clearer (by removing white space between the panels).\n\nPost rebuttal : I increased my score from 6 to 8 given that D-CODE allows to approximate the latent ODE even when the smoother does not perform well. ", "summary_of_the_paper": "Authors propose D-CODE, a framework to approximate closed-form ODEs using symbolic regression. The method first estimate the trajectory by smoothing the observations (using any convenient smoother) before inferring a closed form (finite sequence of operations) transition function. This allows for better interpretability of the system. \n\nAuthors test D-Code on 5 simulated datasets, where they show high performance of the method, and on one real dataset where they use D-Code to model the temporal effect of chemotherapy on the tumor volume, and compare it to Total variation regularized differentiation (SR-T), the closest method in the literature. ", "summary_of_the_review": "I would lean towards acceptance of the paper. Using symbolic regression to approximate closed-form ODEs is novel and very interesting, and the results presented in the paper are compelling. However, I would want to see more comparisons to other methods, as well as results on more datasets to gain a better understanding of when this method will be useful and outperform existing methods. I am hoping that the authors will be able to to answer my concerns by presenting more results. ", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper introduces a new technique for discovering closed-form functional forms (ordinary differential equations) that explain noisy observed trajectories x(t) where the "label" x'(t) = f(x(t), t) is not observed, but without trying to approximate it. The method first tries to approximate a smoother trajector x^hat(t), then relies on a variational formulation using a loss function over functionals {C_j}_j, defined in terms of an orthonormal basis {g_1, …, g_S} of sampling functions such that the sum of squares of all the C_j approximates the theoretical distance between f(x) and the solution f*(x). These sampling functions are typically chosen to be a basis of sine functions. The method is evaluated on several canonical ODEs (growth model, glycolitic oscillator, Lorenz chaotic attractor) and compared to gaussian processes-based differentiation, to spline-based differentiation, regularised differentiation, and applied to model the temporal effect of chemotherapy on tumor volume. 

Reviewers found that the paper was well-motivated and easy to follow (EBvJ), well evaluated (EBvJ), offering new perspectives to symbolic regression (79Ft). Reviewer vaG3 had their concerns addressed. Reviewer ZddY had concerns about the running time (a misunderstanding that was clarified) and the lack of comparison to a simple baseline consisting in double optimisation over f and x^hat(0) using Neural ODEs (the authors have added a Neural ODE baseline but were in disagreement with ZddY and 79Ft about their limitations).

Reviewers engaged in a discussion with the authors, and the scores are 6, 6, 8, 8. I believe that the paper definitely meets the conference acceptance bar and would advocate for its inclusion as a spotlight in the conference.
Decision: Accept (Spotlight)