D-CODE: DISCOVERING CLOSED-FORM ODES FROM OBSERVED TRAJECTORIES
Year: 2022

Abstract:
For centuries, scientists have manually designed closed-form ordinary differential equations (ODEs) to model dynamical systems. An automated tool to distill closedform 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 (ai, bi) as training examples. However, an ODE models the time derivative ẋ(t) of a dynamical system, e.g. ẋ(t) = f(x(t), t), and the “label” ẋ(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 Closedform ODE framework (D-CODE), which advances symbolic regression beyond the paradigm of supervised learning. D-CODE uses 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.

heading: 1 INTRODUCTION
text: An ordinary differential equation (ODE) links the state of a continuous-time dynamical system x(t) to its time derivative ẋ(t) via a function f , e.g. ẋ(t) = f(x(t), t). The ODE is closed-form when f has a concise and analytical expression (e.g. f(x, t) = −t log(x)). Closed-form ODEs are ubiquitous in science and engineering. The popularity is partly due to their transparency and interpretability to human experts (Petersen et 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 a tool would accelerate the quantitative modeling of dynamical systems, with applications ranging from discovering the kinetics of biochemical reactions to modeling tumor growth (Schuster, 2019; Geng et al., 2017).
Symbolic regression is an established approach to discover a closed-form function a = f(b) for predicting the label a from the feature b. It operates in the supervised learning setting and requires a dataset of label-feature pairs for training, i.e. D = {ai, bi}Ni=1 (Schmidt & Lipson, 2009). However, symbolic regression is not directly applicable to ODEs because we usually do not observe the time derivative ẋ(t), i.e. the labels for regression are not available in the data. In contrast, we only have access to the measurements of state y(t) = x(t) + 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.
A simple way to resolve this mismatch is to estimate ẋ(t) from data (e.g. by numerical differentiation (Bickley, 1941)), and use the estimated derivatives as labels for regression (Gaucel et al., 2014).
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 (Žegklitz & Pošı́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 x(t) and the ODE f while bypassing the unobservable time derivative ẋ(t). We develop a 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 to measurement artifacts than the alternative methods. Finally, D-CODE is designed as a general framework, where some of its components can be flexibly adapted based on the application.

heading: 2 BACKGROUND AND PROBLEM SETTING
text: Curating a dataset of trajectories for ODE discovery involves many decisions, e.g. What variables to include? When and for how long to 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 J ∈ N+ variables is defined as
ẋj(t) = fj(x(t)), ∀j = 1, . . . , J, ∀t ∈ [0, T ] (1)
where we use Newton’s notation ẋj(t) for the time derivative. The functions fj : RJ → R will be sometimes referred to as the ODEs directly. We denote T ∈ R+ as the maximum time horizon we have data for. We highlight the following distinction: the trajectory xj : [0, T ]→ R is a function of time, whereas the state xj(t) ∈ R, ∀t ∈ [0, T ] is a point on the trajectory1. We denote the state vector x(t) := [x1(t), . . . , xJ(t)]
> ∈ RJ and the vector-valued trajectory function x := [x1, . . . , xJ ]. Let f∗j ’s be the true but unknown ODEs to be uncovered, and xi : [0, T ]→ RJ , i ≤ N , N ∈ N+ be the true trajectories that satisfy f∗j ’s. In practice, we only measure the true trajectories at discrete times and with noise. Denote the measurement of trajectory i at time t as yi(t) ∈ RJ ; we assume
yi(t) = xi(t) + i(t), ∀i ≤ N, t ∈ T (2)
where i(t) ∈ RJ is zero-mean noise with standard deviation σ. The measurements are made at time t ∈ T = {t1, t2, . . . , T}. We denote the dataset as D = {yi(t)|i ≤ N, t ∈ T }. Closed form. The function fj : RJ → R has a closed form if it can be expressed as a finite sequence of operations (+,÷, log, . . .), input variables (x1, x2, . . .) and numeric constants (1.5, 0.8, . . .) (Borwein et al., 2013). The functional form of fj is the expression with all the numeric constants replaced by placeholders θk’s, e.g. fj(x) = θ1x · log(θ2x). To fully uncover f∗j , we need to infer its functional form and estimate the unknown constants θk’s (if any).
Variational formulation. The variational formulation provides a direct link between the trajectory x and the ODE fj without involving ẋ (Hackbusch, 2017). We start with the following definition.
Definition 1. Consider J ∈ N+, T ∈ R+, continuous functions x : [0, T ]→ RJ , f : RJ → R, and g ∈ C1[0, T ], where C1 is the set of continuously differentiable functions. We define the functionals
Cj(f,x, g) := ∫ T 0 f ( x(t) ) g(t)dt+ ∫ T 0 xj(t)ġ(t)dt; ∀j ∈ {1, 2, . . . , J} (3)
Importantly, the functional Cj depends on the testing function g(t) and its derivative ġ(t) but not ẋj . Proposition 1 below provides the variational formulation of ODE. Essentially it specifies the infinitely many constraints that x has to satisfy in order to be a solution to the ODE. (see Appendix A).
1We restrict the domain of the function xj to [0, T ] (rather than R+) because we only have data up to T . However, after learning the ODE, we can extrapolate xj beyond T (Appendix D).
Proposition 1. (Hackbusch, 2017) Consider J ∈ N+, T ∈ R+, a continuously differentiable function x : [0, T ]→ RJ , and continuous functions fj : RJ → R for j = 1, . . . , J . Then x is the solution to the system of ODEs in Equation 1 if and only if
Cj(fj ,x, g) = 0, ∀j ∈ {1, . . . , J}, ∀g ∈ C1[0, T ], g(0) = g(T ) = 0 (4)

heading: 3 RELATED WORK

heading: 3.1 SYMBOLIC REGRESSION
text: Symbolic regression attempts to uncover a closed-form prediction function a = 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 on developing optimization algorithms. Genetic programming has been one of the most successfully and widely used methods (Koza, 1994; Schmidt & Lipson, 2009). It represents f as a tree where the internal nodes are operations and the leaves are variables or constants (Figure 1 B). It then applies genetic algorithm to search for the best tree representation (Forrest, 1993). More recently, AI Feynman introduces a set of heuristics to directly prune the search space, thereby improving the optimization efficiency (Udrescu & Tegmark, 2020). Optimization methods based on reinforcement learning (Petersen et 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) as 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 equation-data mismatch and the sensitivity to noisy or infrequent observations.

heading: 3.2 DATA-DRIVEN DISCOVERY OF CLOSED-FORM ODES
text: Two-step sparse regression performs sparse regression on the estimated time derivative ̂̇xj (Brunton et al., 2016; Rudy et al., 2017). It assumes that the true ODE can be symbolically expressed in a linear form: ẋj(t) = ∑K k=1 θkhk(x(t)) where θk ∈ R are unknown constants and hk : RD → R are pre-specified candidate functions, such as monomials x, x2, . . . (Table 1). The method uses L1 regularization on θk 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. x3/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 ̂̇xj as the label but applies symbolic regression instead (Gaucel et al., 2014) (Table 1). It may use any optimization method discussed in Section 3.1 to search for the optimal function. In this way, it removes the linear assumption on the functional form and can discover closed-form functions in general.
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-art methods for estimating derivatives, this approach still often fails when the noise is relatively low. We briefly review the methods for estimating derivatives from y(t) in Appendix E and show why this is generally a very challenging problem in Section 4.2.
The function approximator approach learns the true ODE with a function approximator f̂ , such as a neural network (NN) (Chen et al., 2018; Rubanova et al., 2019) or a Gaussian process (GP) (Heinonen et al., 2018). It estimates the unknown initial condition x̂(0) and predicts the entire trajectory x̂ 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 a 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.

heading: 4 METHOD

heading: 4.1 THE D-CODE ALGORITHM
text: D-CODE consists of a preprocessing step and an optimization step. Figure 1 provides a 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 x̂i : [0, T ] → RJ , i ≤ N from the noisy and discretely-sampled data D as an approximation to the true trajectories xi. Denoising and interpolating signals is 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 is agnostic to the exact choice of 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 x̂i, we search for the function f̂j , ∀j ≤ J that is consistent with x̂i. Specifically, we solve the following optimization problems for all j ≤ J :
f̂j = arg min f N∑ i=1 S∑ s=1 Cj(f, x̂i, gs) 2, (5)
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 gs ∈ C1[0, T ], s ≤ S (to be discussed in Section 4.2). Comparison with two-step symbolic regression. As depicted in Figure 1, D-CODE and two-step symbolic regression have two main differences. (1) The estimands in the pre-processing steps are different: D-CODE estimates the true trajectory x 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.

heading: 4.2 THEORETICAL RESULTS
text: 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 f∗. 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 metrics induced by Lp norms of the function space. However, the Lp norms consider the values in the entire domain RJ . This is a disadvantage because most dynamical systems in nature operate within specific range and scale, and RJ contains areas where we do not have knowledge about or cannot collect data from. A natural solution is to restrict the comparison on the trajectory x, which can be achieved by composing f and f∗ with x. We propose the following distance function:
dx(f, f ∗) := ||f ◦ x− f∗ ◦ x||2 = ||(f − f∗) ◦ x||2 (6)
where ◦ denotes function composition and x is a true trajectory satisfying the ODE f∗. In Appendix A, we further justify and discuss the properties of dx(f, f∗). However, in practice, the dx(f, f∗) cannot be computed from data because it depends on the unknown f∗ and x. Next, we show that the objective in Eq. 5, which is computable from data, can be used to approximate this distance.
Convergence to distance. The first step of our algorithm estimates a trajectory x̂ from the measurements {y(t)|t ∈ T } as an approximation to the true trajectory x. For a 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 1 shows that given high enough sampling frequency and large enough number of testing functions our objective converges to the squared distance dx(f, f∗)2, which measures how different f is from f∗. That justifies our objective. Theorem 1. Consider J ∈ N+, j ∈ {1, . . . , J}, T ∈ R+. Let f∗ : RJ → R be a continuous function, and let x : [0, T ] → RJ be a continuously differentiable function satisfying ẋj(t) = f∗(x(t)). Consider a sequence of functions (x̂k), where x̂k : [0, T ] → RJ is a continuously differentiable function. If (x̂k) converges to x in L2 norm. Then for any Lipschitz continuous function f
lim S→∞ lim k→∞ S∑ s=1 Cj(f, x̂k, gs) 2 = dx(f, f ∗)2, (7)
where {g1, g2, . . . } is a Hilbert (orthonormal) basis for L2[0, T ] such that ∀i, gi(0) = gi(T ) = 0 and gi ∈ C1[0, T ].
We provide the proof of Theorem 1 in the Appendix A.
Comparison with other methods. Theorem 1 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 1 place constraints on the type estimation algorithm, i.e., the estimated derivative ˆ̇x needs to converge to the true derivative ẋ. 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
and differentiation reduces its differentiability class. Namely, if a true trajectory has a maximum number of m continuous derivatives then its derivative ẋ has only m− 1 continuous derivatives. Choice of testing functions. As stated in Theorem 1, we should use testing functions that are orthonormal basis for L2[0, T ] and satisfy g(0) = g(T ) = 0. Ideally, they should also have analytical time derivatives ġ(t) for efficient computation of the functional C (Eq. 3). Hence, a natural choice is the sine functions gs(t) = √ 2/T · sin(sπt/T ), 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).

heading: 5 EXPERIMENTS AND EVALUATION
text: 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 biochemistry (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 to a oscillatory limit cycle. Finally, we 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.
Measurement settings. For each dynamical system, we consider different measurement settings specified by (1) the measurement noise level σR, (2) the sampling step size ∆t, and (3) the number of trajectories N . Since the dynamical systems have different scales, we will report the noise-tosignal ratio σR = σ/std(x(t)) for ease of comparison. We sample the system at regular intervals T = {∆t, 2∆t, . . . 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 xi(0) from a uniform distribution (specified in Appendix F). Then we obtain the true xi by solving the IVP computationally. Finally, we obtain yi(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. θ1 + x1 ≡ x1 + θ1). (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̂ approximates f∗ using the distance function in Equation 6 (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), Spline-smoothed differentiation (SR-S) (Ahnert & 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-processing 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 it does not give closed-form expressions. Appendix F contains further implementation details.
2The code is available at https://github.com/ZhaozhiQIAN/D-CODE-ICLR-2022.

heading: 5.1 RESULTS
text: Gompertz and generalized logistic growth models. These two models are widely used to capture asymmetric growth with saturation. Their governing equations are as follows:
ẋ(t) = −θ1x(t) · log ( θ2x(t) ) Gompertz Model
ẋ(t) = θ1x(t) · ( 1− x(t)θ2 ) Generalized Logistic Model
where θ1, θ2 > 0 are numeric constants. It is worth highlighting that the generalized logistic model involves a constant θ2 ∈ R+ in the power of x(t). Therefore, it cannot be symbolically expressed in a sparse linear form ∑K k=1 θkhk(x(t)) for some candidate functions hk, and the two-step sparse regression method cannot discover the exact functional form of this equation (Section 3.2).
The main simulation results are 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 in the low noise, small step size settings (σR = 0.01), where the time derivative ẋ 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 in the benchmarks compared to 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 A shows a typical observed trajectory with moderate measurement noise (σR = 0.1). Here, the Gaussian process can produce a good estimate x̂(t) by de-noising the measurements. Its posterior interval covers the true trajectory x(t). Hence, x̂(t) is suitable to guide the search for the ODE, as is done in D-CODE. On the other hand, the moderate measurement noise on x(t) translates into very high noise on the time derivative ẋ(t) (panel B). Even with the state-of-the-art numerical differentiation methods used by SR-T and SR-S, the estimated derivative is still biased on certain
intervals (panel C). Thus, applying symbolic regression on the biased estimates is unlikely to perform well, as is shown in the benchmark results. The same reasoning also applies to the scenario with large sampling step size ∆t, another challenging situation for estimating the time derivatives.
Glycolytic oscillator. This biochemical oscillator consists of two ODEs (Sel’Kov, 1968):
ẋ1(t) = θ1 − θ2x1(t)− x1(t)x2(t)2 (8) ẋ2(t) = −x2(t) + θ3x1(t) + x1(t)x2(t)2 (9)
where θ1, θ2, θ3 > 0 are numeric constants. Compared with the growth models above, the glycolytic oscillator involves more variables and has more complex expressions. The two variables in the equation interact in a highly nonlinear way through the term x1(t)x2(t)2.
The simulation results for the glycolytic oscillator are shown in Table 2. Here we report the success probability in addition to the RMSE (×10−2) on the constants θ2, θ3, 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 σR = 0.2 (The benchmarks’ RMSE are NA in this setting because they never recovered the correct functional form).
Chaotic Lorenz system. The Lorenz system is a model system for chaotic dynamics, defined as: ẋ1(t) = θ1 ( x2(t)− x1(t) ) ; ẋ2(t) = x1(t) ( θ2 − x3(t) ) − x2(t); ẋ3(t) = x1(t)x2(t)− θ3x3(t)
We set the constants as θ1 = 10, θ2 = 28, θ3 = 8/3, which leads to chaotic behaviour.
The simulation results for the chaotic Lorenz system are 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 a simulated trajectory in the phase space using the true or estimated ODEs (with σR = 0.2). D-CODE generates a 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 a very short prediction horizon (0.04), it still fails to capture the non-periodic behavior of the system.

heading: 6 D-CODE IN ACTION
text: 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 to drug resistance. Although there exist many models for unintervened tumor growth, few models can capture the effect of chemotherapy or the possibility of relapsing—in many cases exponential decay (log cell-kill model) is used for a crude approximation (Geng et 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 t = 0, and the trajectory of tumor volume are measured for up to one year T = 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 310 trajectories, among which 200 is used for training and 110 for evaluation.
The following two ODEs are discovered by D-CODE and SR-T.
ẋ(t) = 4.48t2x(t) + log(t); D-CODE ẋ(t) = 4x(t) log ( tx(t) ) log ( tx(t) + 2t ) ; SR-T
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 t = 0.6 and predicts that, under no further intervention, the tumor volume would keep growing after t > 0.6 until patient death. In comparison, SR-T’s trajectory is almost constant.
We also quantitatively evaluate the prediction errors on the testing set. Given the tumor size measured at t = 0, we solve the IVP to obtain the full trajectory x and compute the RMSE between the estimated and the observed trajectory. D-CODE achieves an error of 0.220 compared with 0.637 of SR-T.

heading: 7 DISCUSSION ON FAILURE MODES AND OPEN CHALLENGES
text: 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) ∈ RL may also influence the temporal dynamics, i.e. ẋj(t) = f ( x(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 is challenging or even impossible without strong assumptions (Ljung, 1998). In this case, D-CODE can still find an ODE based on the observed variables—this 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, we 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 or step size is too large. However, as we show experimentally, D-CODE is more robust to measurement artifacts than the alternative discovery methods.

heading: ACKNOWLEDGMENTS
text: This work is supported by the US Office of Naval Research (ONR), the National Science Foundation (NSF Grant number: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.
Reproducibility Statement. The proof of theoretical results is shown in Appendix A. The technical assumptions for the proof are 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 6 and further detailed in Appendix F. The experiment code is available at https://github.com/ZhaozhiQIAN/ D-CODE-ICLR-2022.

heading: TABLE OF SUPPLEMENTARY MATERIALS
text: 1. Appendix A: proof and discussion on the theoretical results.
2. Appendix B: implementation details of D-CODE.
3. Appendix C: discussion on dataset curation.
4. Appendix D: utilities of 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 on alternative methods.

heading: A THEORETICAL RESULTS: PROOF AND DISCUSSION
text: A.1 VARIATIONAL FORMULATION OF ODE
Before we prove Proposition 1, we need the following lemma which is a particular formulation of the Fundamental lemma of calculus of variations Elsgolc (1961).
Lemma (Fundamental lemma of calculus of variations). Let h be a continuous function on a closed interval [0, T ]. h is equal to 0 everywhere if and only if ∫ T 0 h(t)g(t) dt = 0 for all g ∈ C1[0, T ] such that g(0) = g(T ) = 0.
Proof. The forward direction is trivial, so let us focus on the converse.
Assume for contradiction that there is a point t0 ∈ [0, T ] such that h(t0) 6= 0. Without loss of generality, assume that h(t0) > 0. From the continuity of h, we know that there is actually a point t1 ∈ (0, T ), such that h(t1) > 0. Continuity of h implies that there exists δ > 0 and a small neighbourhood (t1 − δ, t1 + δ) ⊂ (0, T ) such that h(t) > 0 ∀t ∈ (t1 − δ, t1 + δ). We can now define a non-negative continuous function g : [0, T ] → [0,+∞) such that g(t) = 0 outside the neighbourhood (t1− δ, t1 + δ) and g(t) > 0 for some interval inside (t1− δ, t1 + δ). See Figure 6 for an illustrative example.
Then ∫ T 0 h(t)g(t) dt = ∫ t1+δ t1−δ h(t)g(t)dt > 0 (10) That contradicts ∫ T 0 h(t)g(t) dt = 0 for all g ∈ C1[0, T ] such that g(0) = g(T ) = 0.
Proof of Proposition 1. Observe that
ẋj(t) = fj(x(t))⇐⇒ fj(x(t))− ẋj(t) = 0 (11)
Now using Fundamental lemma of calculus of variations, we get that the following two statements are equivalent
fj(x(t))− ẋj(t) = 0 ∀t ∈ [0, T ]∫ T 0 (fj(x(t))− ẋj(t))g(t)dt = 0 ∀g ∈ C1[0, T ], g(0) = g(T ) = 0 (12)
By linearity and integration by parts, we get∫ T 0 (fj(x(t))− ẋj(t))g(t)dt = ∫ T 0 fj(x(t))g(t)dt− ∫ T 0 ẋj(t)g(t)dt
= ∫ T 0 fj(x(t))g(t)dt+ ∫ T 0 xj(t)ġ(t)dt− xj(T )g(T ) + xj(0)g(0)
= ∫ T 0 fj(x(t))g(t)dt+ ∫ T 0 xj(t)ġ(t)dt
(13)
where in the last equality we used the fact that g(0) = g(T ) = 0. That means that∫ T 0 (fj(x(t))− ẋj(t))g(t)dt = Cj(fj ,x, g) (14)
This proves that ẋj(t) = fj(x(t)) ∀j ∈ {1, . . . , J} (15)
if and only if
Cj(fj ,x, g) = 0, ∀j ∈ {1, . . . , J}, ∀g ∈ C1[0, T ], g(0) = g(T ) = 0 (16)
Discussion. As we show in the proof above, the key idea of the variational formulation is to bypass the derivative of the trajectory ẋj(t) with the derivative of a testing function ġ(t) — 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 ẋj(t) (which is unobserved and challenging to estimate), but only need to know ġ(t) for a given testing function g(t). Indeed, if g(t) has an analytical derivative, we immediately have access to ġ(t).
Note that the variational formulation builds on the exact condition Cj = 0. It does not consider the situation where Cj 6= 0, no matter how close Cj is to 0. Hence, the formulation itself is inadequate to guide the search of closed-form ODEs. For instance, suppose the Cj for two candidate functions f̂1 and f̂2 are 0.01 and 1000 respectively. The variational formulation does not imply f̂1 fits the trajectory x better than f̂2; rather, it implies that f̂1 and f̂2 are equally bad as Cj 6= 0 for both. In order to bridge this gap, we need to establish a 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 is the focus of the current work. However, the variational formulation can be generalized to other types of differential equations, e.g partial differential 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 a known ODE or DDE (Brunel et al., 2014). We envision that future works may be built on the general variational formulation to discover closed-form PDE or DDE from data.

heading: A.2 DISTANCE FUNCTION
text: In this section, we justify that dx(f, f∗) is an appropriate metric for our problem.
There are two main ways of comparing functions. We can either compare their symbolic representations or measure the differences in values they take for a set of arguments. As we want this metric to be a limit of some objective computable from the data, we focus on the latter.
Commonly used distance functions are metrics induced by Lp norms. Lp norm is defined in the following way:
||g||p = (∫ b a |g(t)|pdt ) 1 p
(17)
In practice, L1 or L2 norm is used. The induced metric is then defined as the norm of the difference between the functions, namely d(f, g) = ||f − g||p. Now the question remains which difference we should take the norm of. We have three reasonable choices:
1. f − f∗
2. x̃− x, where x̃ is trajectory prescribed by f computed from an estimated initial condition 3. f − f∗ composed with some other function
The first choice poses the challenge of choosing the range [a, b] over which we integrate f − f∗. In general we do not make any assumptions on the domain of the functions. It 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 a 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. If we compose f − f∗ with the ground truth trajectory, x, not only we have a well-defined interval over which to integrate ([0, T ]) but also we focus our attention on values of f∗ where we have a chance to recover it, and where it matters for us the most.
That is why we decided to chose dx(f, f∗) = ||(f − f∗) ◦ x||2 as our distance function. Moreover, as it is based on a familiar metric, it inherits many important properties, such as
1. dx is a pseudo-metric3, in particular, f = f∗ =⇒ dx(f, f∗) = 0 2. dx(f, f∗) = 0 =⇒ f = f∗ on x([0, T ]), the image of x 3. For a sequence fl such that dx(fl, f∗)→ 0 it holds that fl ◦ x→ f∗ ◦ x in the mean4
The third property tells us that if we keep minimising dx(f, f∗), we can expect f getting closer to f∗ on x([0, T ]), the image of x.

heading: A.3 PROOF OF THE THEOREM
text: Proof of Theorem 1. Let us consider a Hilbert space L2[0, T ] where the inner product is denoted by 〈·, ·〉 and defined as:
〈a, b〉 = ∫ T 0 a(t)b(t)dt (18)
where a and b are real functions which are square-integrable on the interval [0, T ]. The L2 norm can be written as:
||a||22 = 〈a, a〉 (19) 3It can be considered to be a metric through metric identification. We would need to declare two functions being equal if they are equal on x([0, T ]) 4This can be strengthened to imply uniform convergence if the sequence fl is equicontinuous (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 are bounded by the same constant.
Using this notation we can rewrite Cj(f, x̂k, gs) as:
Cj(f, x̂k, gs) = 〈f ◦ x̂k, gs〉+ 〈x̂kj , ġs〉 (20)
We assume that x̂k converges to x in mean. For vector-valued functions that means that∫ T 0 ||x̂k(t)− x(t)||22 dt→ 0 as k →∞ (21)
Note that || · ||2 in this equation denotes the usual Euclidean norm as the arguments are vectors not functions. This is equivalent to saying that all components of x̂ converge to the corresponding components of x. Thus
||x̂kj − xj ||2 → 0 as k →∞∀j (22) As the inner product with one argument fixed can be seen as a continuous operator, we get:
〈x̂kj , ġs〉 → 〈xj , ġs〉 as k →∞ (23)
We assume that fj is λ-Lipschitz continuous. That means that
||f ◦ x̂k − f ◦ x||22 = ∫ T 0 |f(x̂k(t))− f(x(t))|2 dt
≤ ∫ T 0 λ2||x̂k(t)− x(t)||22 dt
= λ2 ∫ T 0 ||x̂k(t)− x(t)||22 dt→ 0 as k →∞
(24)
This means that 〈f ◦ x̂k, gs〉 → 〈f ◦ x, gs〉 as k →∞ (25)
Combining equations (23) and (25) with expression (20), we conclude that
lim k→∞
Cj(f, x̂k, gs) 2 = (〈f ◦ x, gs〉+ 〈xj , ġs〉)2 (26)
Now, using integration by parts we get 〈xj , ġs〉 = ∫ T 0 xj(t)ġs(t)dt = − ∫ T 0 ẋj(t)gs(t)dt+ xj(T )gs(T )− xj(0)gs(0) (27)
We use the assumption that gs(0) = gs(T ) = 0 ∀s to obtain
〈xj , ġs〉 = −〈ẋj , gs〉 (28)
We can substitute it back into Equation 26 and use the linearity of the inner product to get
lim k→∞
Cj(f, x̂k, gs) 2 = 〈f ◦ x− ẋj , gs〉2 (29)
Now we can use Parseval’s identity and the fact that the set g1, g2, ... forms a Hilbert basis for L2[0, T ] to infer that:
lim S→∞ S∑ s=1 〈f ◦ x− ẋj , gs〉2 = ||f ◦ x− ẋj ||22 (30)
Substituting Equation 29 into Equation 30 and using the fact that the inner limit does not depend on S we obtain the statement in the theorem
lim S→∞ lim k→∞ S∑ s=1 Cj(f, x̂k, gs) 2 = ||f ◦ x− ẋj ||22 = ||(f − f∗) ◦ x||22 (31)

heading: A.4 DIVERGENT DERIVATIVES
text: The main advantage of our algorithm lies in estimating the trajectory instead of estimating its time derivative. We show how the derivatives of a 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 (hk) where hk : [0, 1]→ R is given by
hk(t) = sin(kt)
k (32)
It can be shown that (hk) converges to h(t) = 0 uniformly. However, the derivative of hk, given by
h′k(t) = cos(kt), (33)
fails to converge at any point.
We use this observation to show that the square of ||ˆ̇xj−fj ◦x̂||2, where x̂j is the derivative calculated from the recovered trajectory x̂, does not converge to dx(fj , f∗j )
2 when x̂ converges to the ground truth trajectory x.
Consider the following ODE (J = 1, T = 1):
ẋ(t) = x(t) ∀t ∈ [0, 1] (34) so f∗(x) = x. We choose x(t) = et as the ground truth trajectory from which we sample. Let (x̂k) be a sequence of trajectories given by x̂k = et + sin(kt) k . It can be shown that (x̂k) converges to x uniformly. Then for any candidate function f
||ˆ̇xk − f ◦ x̂k||22 = ∫ 1 0 (ˆ̇xk(t)− f(x̂k(t)))2dt
= ∫ 1 0 ( et + cos(kt)− f ( et + sin(kt) k ))2 dt
(35)
If ||ˆ̇xk− f ◦ x̂k||22 were to converge to dx(f, f∗)2 then in particular ||ˆ̇xk− f∗ ◦ x̂k||22 should converge to 0. However, we show this is not what happens.
||ˆ̇xk − f∗ ◦ x̂k||22 = ∫ 1 0 ( et + cos(kt)− f∗ ( et + sin(kt) k ))2 dt
= ∫ 1 0 ( et + cos(kt)− ( et + sin(kt) k ))2 dt
= ∫ 1 0 ( cos(kt)− sin(kt) k )2 dt
= ∫ 1 0 cos2(kt) + ( sin(kt) k )2 − 2cos(kt) ( sin(kt) k ) dt
(36)
As ( sin(kt) k )2 and 2cos(kt) ( sin(kt) k ) converge uniformly to a 0 function, we can use the linearity of integration to obtain
lim k→∞ ||ˆ̇xk − f∗ ◦ x̂k||22 = lim k→∞ ∫ 1 0 cos2(kt)dt = lim k→∞ ( sin(2k) 4k + 1 2 ) = 1 2 6= 0 (37)
This example shows that objectives that estimate the derivative are highly dependent on the type of estimation algorithm used. In contrast, Theorem 1 tells us that the only property of the estimation algorithm we need is that the trajectories converge to the ground truth in L2 norm.
Algorithm 1 D-CODE. Input: Dataset D = {yi(t)|i ≤ N, t ∈ T } Input: Smoothing algorithm S, numerical integration algorithm I, optimization algorithm O Input: Testing functions gs, s ≤ S, initial guesses fj0, j ≤ J . x̂i = S ( yi(t1),yi(t2), . . . ) , i ≤ N {Preprocessing step: smooth input trajectories}
for j ∈ [1, J ] do f̂j = fj0 {Initialize the jth ODE} Converge = False while Not Converge do obj = ∑N i=1 ∑S s=1 Cj(fj , x̂i, gs)
2 {Calculate objective (Eq. 5).} f̂j ,Converge = O(obj, f̂j) {Run optimization step.}
end while end for Output: The discovered ODEs f̂j , j ≤ J

heading: B D-CODE FRAMEWORK

heading: B.1 PSEUDOCODE
text: The Pseudocode is presented in Algorithm 1.

heading: B.2 DISCOVER HIGHER ORDER AND NONAUTONOMOUS ODES
text: In the main text, we focused on the system of first order autonomous ODEs (Equation 1). This is 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):
ẍ1 + ẋ1 + x1 − cos(t) = 0
The ODE is second-order as it involves the acceleration term ẍ1 and it is nonautonomous because of the time-dependent term cos(t). If we introduce two new variables: x2 = ẋ1 and x3 = t, we obtain the following first order autonomous system:
ẋ1 = x2
ẋ2 = −x2 − x1 + cos(x3) ẋ3 = 1
Since the first and the third equations are implied by the definition of x2 and x3, they do not need to be discovered from data. One can apply the discovery algorithm to uncover the second equation from the measurements of x1, x2 and x3.

heading: C DISCUSSION ON DATASET CURATION
text: The success of ODE discovery depends on the availability of a well-curated dataset. However, curating a dataset is a highly nontrivial task, and it is well beyond the scope of the current work. In fact, the curation of dataset is sometimes an important scientific contribution in itself (Blasius et al., 2020). Here, we provide a 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 on 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 a 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 a challenging high-dimensional problem). We note that there is a well-developed literature on uncovering the interaction network of a high-dimensional dynamical systems (Mangan et al., 2016; Qian et al., 2020).
Measurement settings. In practice, the measurement settings are often affected by practical constraints, 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. If we are interested in studying the cyclic behavior of the system (e.g. seasonality), the measurement should cover at least one entire period (one year). On the other hand, if we are interested in the system’s response to an external stimuli (e.g. short-term weather change due to 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 a wrong stationary point (as indicated by domain knowledge), the researcher may curate an additional dataset with longer time horizon and re-train D-CODE with the additional information.

heading: D UTILITIES OF THE DISCOVERED CLOSED-FORM ODE
text: 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 a sequence of past observations y(t1),y(t2), . . . ,y(tk), we would like to predict x(t∗) for a future time t∗ > tk. This prediction problem is very well studied and has many proven algorithms (e.g. the family of Kalman filters (Welch et 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 x(0) and different external signals (if applicable) to generate a range of trajectories. The simulated trajectories may be used in downstream tasks, e.g. computing the statistical properties.
Planning and control. In 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 of research have been devoted to deriving the optimal policy from the closed-form governing equations and the utility function (Mehrmann, 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 so on. 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).

heading: E BACKGROUND ON NUMERICAL DIFFERENTIATION
text: The task of estimating time derivative ẋ from discretely sampled and noisy data is 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
ˆ̇x = y(t+ ∆t)− y(t−∆t)
2∆t
The main disadvantage of finite difference is that it is highly sensitive to measurement noise (Chartrand, 2011).
A common way to address this shortcoming is to estimate the derivative on the smoothed trajectory. For example, spline-smoothed differentiation smooths the observed trajectories with cubic spline and estimates ẋ 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
the smoothness of the trajectory and its (higher-order) derivatives (Zhou & Wolfe, 2000): they encode a 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 use explicit regularization on the finite difference estimator. As we pointed out in Appendix A, the standard L2 regularization is not suitable for derivative estimation. The authors propose to use total variation regularization, a method originally developed in computer vision to obtain sharp denoised images (Strong & Chan, 2003). This type of regularization allows rapid (or even discountinous) change in ẋ. 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.
F IMPLEMENTATION DETAILS OF THE EXPERIMENTS

heading: F.1 DETAILED SETTINGS FOR EACH EXPERIMENT
text: The detailed settings for each experiment in Section 5 is shown in Table 3. The time horizon T and the range of initial conditions [a, b] are chosen based on the scale and properties of the system. For example, we observe that the Gompertz model reaches saturation and converges to the constant x = 1 before t = 4 (hence we set the time horizon T = 4). 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).

heading: F.2 HYPER-PARAMETER SETTINGS
text: 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 log-likelihood (Bernardo et al., 1998).
Genetic programming. The hyperparameters of genetic programming is decided based on a pilot study. In the pilot study, we consider a noiseless and frequently sampled Gompertz growth model (σ = 0, ∆t = 0.01). We use randomized search on Two-step Symbolic Regression to find the hyperparameters that achieves the best RMSE loss (Table 1). The same hyperparameters are used by all algorithms on all experiments. The values are listed below. We refer the reader to the documentation of gplearn library for a detailed explanation of these hyperparameters (Stephens, 2021).
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
The set of allowed mathematical operations is
{+,−,×,÷, log, cos, power}
D-CODE and all benchmarks use 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 for numerical differentiation (Quade & Goldschmidt, 2020). The library contain implementations of the methods discussed in Appendix E. Currently, there lacks a principled way to tune the hyperparameters of numerical differentiation because the time derivative is not observed. Existing works often rely on heuristics (Chartrand, 2011). 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 as D-CODE, which is tuned by the marginal log-likelihood on the observed trajectory.
Correctness of the functional form. As discussed in Section 2, the functional form of a closedform function f is the expression with all the numeric constants replaced by placeholders θ1, θ2, . . . In practice, the uncovered function f̂ may contain constants that are close to 0. For instance f(x) = x + 0.001x2 is close to f(x) = x when x 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 to be helpful). The thresholding is applied to 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 dx. We use Python package sympy to decide whether the discovered functional form is correct (Meurer et al., 2017), i.e. success when simplify(f∗ − f̂) returns 0.

heading: F.3 CLEANING TUMOR GROWTH DATASET
text: 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.

heading: G ADDITIONAL RESULTS
text: G.1 ADDITIONAL EVALUATION METRIC: dx
Here we report the distance function dx(f̂ , f∗) between the true and estimated ODEs. This metric tracks how well f̂ 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 dx in addition to the discovery probability.
The results are presented in Figure 7 and 8. We observe that the patterns of dx 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 frequent-sampling scenarios, but D-CODE significantly outperforms the benchmarks in the more challenging measurement settings. The agreement between dx and discovery probability suggests that when the correct functional form is not recovered, the estimated f̂ tends to perform poorly in approximating the true function as well.

heading: G.2 ADDITIONAL DYNAMICAL SYSTEMS
text: Here we consider an additional dynamical system, which is an extension to the quadratic integrate and fire (QIF) model in Neuroscience (Brunel & Latham, 2003). The governing equation is
ẋ = x2
t+ θ (38)
We highlight that this ODE is time-dependent and it involves a fraction, which makes it impossible to be symbolically represented in a linear form f(x) = ∑ k θkbk(x) (as required by the two-step sparse regression method). A similar example was given in Appendix B-2 of Brunton et al. (2016) to illustrate this failure mode of two-step sparse regression method.
The simulation results are shown in Figure 9. We observe the same pattern as the other dynamical systems considered in the main text.

heading: G.3 SENSITIVITY ANALYSIS FOR TESTING FUNCTION
text: 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
conditions in 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 too few testing functions (i.e. < 40) leads to performance loss. This agrees with Theorem 1 because the objective function Eq. 5 cannot well-approximate the distance dx (Eq. 6) with too few testing functions. We also observe that the performance is 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 error 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 a smaller integration step size or a more advanced numerical integration method.
For practical applications, the user may experiment with different numbers of 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.

heading: G.4 TRAINING TIME
text: 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 on average. This can still be a 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.

heading: G.5 TABLE OF DISCOVERED ODES
text: The top-3 most common equations discovered by different methods for each dynamical system (among the 100 simulation runs) are listed in the table on the next page.
Published as a conference paperatIC L R 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-simplistic equations for the growth models (constant or zero), possibly because the estimated derivative ̂̇x is too noisy to learn useful signal. (2) SR-T and SR-S may consistently miss a variable (e.g. the x0 term in Glycolytic 1) or miss an interaction (e.g. the x0x2 term in Lorenz 2). (3) SR-T and SR-S may include unnecessary terms in the equation (e.g. Lorenz 3). The measurements settings are the default ones listed in Table 3.
D-CODE SR-T ST-S System Discovered Equation f̂ Freq. Discovered Equation f̂ Freq. Discovered Equation f̂ Freq. Gompertz −θ1x log(θ2x) 1 θ1 0.63 θ1 0.67 - - θ1x− θ2xθ3 0.04 −θ1xθ2(xθ3 + log(x)) 0.03 - - x2 log(x) log(x2) 0.04 − log(log(θ1 + x) 0.03 Logistic θ1x− xθ2 1 0 0.96 0 0.98 - - θ1 0.02 θ1 0.02 - - θ1xθ2 − xθ3 0.02 - -
Glycolytic 1 −θ1x0 + θ2 − x0x21 0.38 θ1 − x0x21 0.5 θ1 − x0x21 0.1 −θ1x0 − θ2x1 + θ3 − x0x21 0.03 −θ1x0 + θ2 − x0x21 0.22 −θ1x0x21 + θ2 0.08 θ1 − x0(θ2 + x21) 0.03 −θ1x0x21 − θ2x0 + θ3 0.04 −θ1x0x1(θ2 + x1) + θ3 0.07
Glycolytic 2 θ1x0 + x0x
2 1 − x1 0.78 θ1x0 + x0x21 − x1 0.75 −θ1x1 + θ2x21 + θ3 + x0x1 0.09
θ1x0 + θ2 + x0x 2 1 − x1 0.13 x0(θ1 + x21)− x1 0.09 θ1x0x21 + θ2x0 − θ3x1 0.08 θ1x0 + θ2x1(θ3x0 + x0x 2 1) 0.01 θ1x0x21 + θ3x0 − θ4x1 0.02 −θ1x1 − θ2x21 + θ3 + x0x1 + x12 0.06
Lorenz 1 −θ1x0 + θ2x1 1 −θ1x0 + θ2x1 1 −θ1x0 + θ2x1 0.63 - - - - −θ1x0 + θ2x1 + θ3 0.37 - - - - - -
Lorenz 2 θ1x0 − x0x2 − x1 0.33 −θ1x0 + θ2x1 + θ3x2 − θ4 0.99 x0(θ1 − x2) 0.73 x0(θ1 − x2) + x0 − x1 0.28 θ1x0(−θ2x2 + θ3) + θ4x1 − x0 0.01 θ1x0 − x0x2 0.2 x0(θ1 − x2)− x1 0.28 - - x0(θ1 − x2)− x0 0.03
Lorenz 3 −θ1x2 + x0x1 0.8 θ1x2 − θ2 + x0x1 + x0 0.46 θ1x2 + x0x1 0.51 −θ1x2 + θ2 + x0x1 0.13 −θ1x2 + x0x1 + x0 0.11 θ1x2 + x0x1 + x0 0.11 −θ1x2 − θ2 + x0x1 0.08 θ1x0x1 − θ2x2 0.1 θ1x2 + x0x1 + x1 0.09
25
G.6 VISUALIZATION OF SMOOTHED AND TRUE TRAJECTORIES
In this section we study whether D-CODE can still discover the ODE when the smoothed trajectory x̂(t) is significantly different from the true trajectory x(t); for instance, when the signal-to-noise ratio is low or the observations are sparse.
In Section 5 we provided quantitative evaluation under scenarios with different noise levels (σR) and sampling sparsity (∆t).5 For further illustration, here we visualize x̂(t) and x(t) together with the observations y(t) in Figure 12. In the most favorable setting (top left), x̂(t) is indistinguishable from x(t), and D-CODE achieves very high success probability as expected. As the observation setting becomes more challenging (middle column), x̂(t) may recover the overall trend of x(t) 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), x̂(t) differs significantly from x(t) for extended periods of time (top right t > 2, bottom right t < 2 and t > 4). Yet, D-CODE still makes successful discovery at times.

heading: H DISCUSSION ON ALTERNATIVE METHODS
text: In this section, we discuss the alternative methods that, in theory, could also be applied to discovering closed form ODEs from data. However, as we will show, these methods suffer from various drawbacks that may severely limit their practical utility.

heading: H.1 DOUBLE OPTIMIZATION METHOD
text: The double optimization method (DO) was first proposed in the literature of ODE parametric estimation, where the functional form of f is known but it contains unknown parameters θ (see Ramsay et al. (2007) Section 1.3 for a review).6 DO solves the following optimization problem (we consider the one-dimensional case to simplify notations):
f̂ = arg min f [ min x̂i(0) N∑ i=1 ∑ t∈T ( yi(t)− x̂i(t) )2] , subject to x̂i(t) = x̂i(0) + ∫ t 0 f(x̂i(s))ds (39)
This is a nested optimization problem where the outer part searches for the unknown function f and the inner part searches for the unknown initial conditions xi(0) given the current guess of f . Note
5Note that, by definition, the signal-to-noise ratio is the inverse of σR, i.e. 1/σR. 6DO is also referred to as the nonlinear least square method (NLS) in some literature.
that the inner and outer parts are coupled because the optimal x̂i(0) depends on f , e.g. the best-fitting exponential curve would have a 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 xi(0) is typically non-convex and may take many steps to converge. To further aggravate the problem, each inner optimization step requires solving the IVP numerically.
The computational challenge becomes even more pronounced in our setting, where the functional form of f is also unknown. Genetic programming needs to evaluate many candidate functions f 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 a candidate f . The results are 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 a feasible method for discovering closed-form ODEs.

heading: H.2 NEURAL ODE WITH POST-HOC DISTILLATION
text: 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 a dataset {xi(t), fN (xi(t))}i,t; (2) apply symbolic regression on this dataset. This approach can also be viewed as a two-step regression method (Section 3), where the derivatives are estimated by a NODE.
The standard approach to train a NODE for time series modeling is to minimize the error over a prediction horizon T ∗ = [t1, t2, . . . , T ∗] (Rubanova et al., 2019), i.e. to minimize
L = N∑ i=1 ∑ t∈T ∗ ||yi(t)− x̂i(t)||22, (40)
where x̂i(t) is obtained by solving the IVP of the NODE.
However, the choice of prediction horizon T ∗ is non-trivial and plays an important role in the training. For chaotic systems, the user should choose a smaller T ∗ 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).7 In practice, the prediction horizon T ∗ is a hyper-parameter that may be difficult to decide (note that the Lyapunov time of the true system is unknown). In contrast, D-CODE does not require user to specify the prediction horizon T ∗.
Furthermore, as we alluded to earlier, this approach is an instance of two-step regression, where the derivative is estimated by NODE instead of numerical differentiation. Although NODE is known to
7The limitations of minimizing short-term prediction loss for sequence modeling is well-recognized in machine learning; for example, the sequence-to-sequence models tackle this issue by minimizing the loss over a sequence of predictions rather than the one step ahead prediction.
be a universal function approximator, its estimation efficiency (e.g. √ n 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 the NODE-S is trained on short term prediction (T ∗ = 0.04). The evaluation is performed on the Gompertz ODE (non-chaotic, one dimensional) and the Lorenz system (chaotic, three dimensional).
Figure 14 shows the success probability (bigger better) and Figure 15 shows the distance dx (Equation 6, smaller better). For the non-chaotic Gompertz ODE, both NODE-L and NODE-S are only able to recover the correct functional form when the noise is very close to zero (σR = 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 dx than NODE-S. This confirms that training over longer time horizon leads to 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 dx is 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 dx monotonically increases with the noise level. This suggests that higher noise level significantly reduces NODE-S’s ability to approximate the unknown ODE, limiting its performance on noisy data.

Review summaries:
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.

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.   Authors 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. 

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).

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.