No Equations Needed: Learning System Dynamics Without Relying on Closed-Form ODEs
Authors: Krzysztof Kacprzyk
Mihaela van der Schaar
Venue: ICLR 2025 Poster (2025)
Keywords: dynamical systems
interpretability
ordinary differential equations
symbolic regression
Primary area: applications to physical sciences (physics, chemistry, biology, etc.)

Abstract: Data-driven modeling of dynamical systems is a crucial area of machine learning. In many scenarios, a thorough understanding of the model’s behavior becomes essential for practical applications. For instance, understanding the behavior of a pharmacokinetic model, constructed as part of drug development, may allow us to both verify its biological plausibility (e.g., the drug concentration curve is non-negative and decays to zero in the long term) and to design dosing guidelines (e.g., by looking at the peak concentration and its timing). Discovery of closed-form ordinary differential equations (ODEs) can be employed to obtain such insights by finding a compact mathematical equation and then analyzing it (a two-step approach). However, its widespread use is currently hindered because the analysis process may be time-consuming, requiring substantial mathematical expertise, or even impossible if the equation is too complex. Moreover, if the found equation's behavior does not satisfy the requirements, editing it or influencing the discovery algorithms to rectify it is challenging as the link between the symbolic form of an ODE and its behavior can be elusive. This paper proposes a conceptual shift to modeling low-dimensional dynamical systems by departing from the traditional two-step modeling process. Instead of first discovering a closed-form equation and then analyzing it, our approach, direct semantic modeling, predicts the semantic representation of the dynamical system (i.e., description of its behavior) directly from data, bypassing the need for complex post-hoc analysis. This direct approach also allows the incorporation of intuitive inductive biases into the optimization algorithm and editing the model's behavior directly, ensuring that the model meets the desired specifications. Our approach not only simplifies the modeling pipeline but also enhances the transparency and flexibility of the resulting models compared to traditional closed-form ODEs.

Published as a conference paper at ICLR 2025
No EQUATIONS NEEDED: LEARNING SYSTEM DYNAM- ICS WITHOUT RELYING ON FORM ODEs
Krzysztof Kacprzyk University of Cambridge 2k7510cas ac.uk uk
Mihaela der Schaar University of Cambridge mv4 720 cam. ac uk
ABSTRACT
Data-driven modeling of dynamical systems is crucial area of machine learning. In many scenarios, thorough understanding of the model's behavior becomes essential for practical applications. For instance, understanding the behavior of a pharmacokinetic model, constructed as part of drug development, may allow us to both verify its biological plausibility (e.g., the drug concentration curve is non- negative and decays to zero in the long term) and design dosing guidelines (e.g., by looking at the peak concentration and its timing). Discovery of closed- ordinary differential equations (ODEs) can be employed to obtain such insights by finding compact mathematical equation and then analyzing it (a two-step approach). However, its widespread use is currently hindered because the analysis process may be coneccon,,,,ng, requiring substantial mathematical expertise, even impossible if the equation is too complex. Moreover, if the found equation's 's behavior does not satisfy the requirements, editing it influencing the discovery algorithms to rectify it is challenging as the link between the symbolic form of an ODE and its behavior can be clusive. This paper proposes a conceptual shift modeling oow-dimensional dynamical systems by departing from the traditional two-step modeling process. Instead of first discovering a dloeed-form equation and then analyzing it, our approach, direct semantic modeling, predicts the semantic representation of the dynamical system (i.e., description of its behavior) directly from data, bypassing the need for complex post-hoc analysis. This direct approach also allows the incorporation of intuitive inductive biases into the optimization algorithm and editing the model's behavior directly, ensuring that the model meets the desired specifications. Our approach not only simplifies the modeling pipeline but also enhances the transparency and flexibility of the resulting models compared to traditional closed-form ODEs.
INTRODUCTION
Background: data-driven modeling of dynamical systems through ODE discovery. Modeling dynamical systems is a pivotal aspect of machine learning (ML), with significant applications across various domains such as physics (Raissi et al. 2019), biology (Neftci & Averbeek 2211), engincering (Brunton & Kutz 2022), and medicine (Lee et al. 2020). In real-world applications, understanding the model's 's behavior is crucial for verification and other domain-specific tasks. For instance, in drug development, it is important to ensure the pharmacokinetic model (Mould & Upton! 2012) is biologically plausible (e.g., the drug concentration is non-negative and decays to zero), and the dosing guidelines may be set up based on the peak concentration and its timing (Han al. 2018). One effective approach to gain such insights is the discovery of closed-form ordinary differential equations (ODEs) (Bongard & Lipson 2007 Schmidt & Lipson 2009, Brunton al. 2016a). where a concise mathematical representation is first found by an algorithm and then analyzed by human.
Motivation: the primary goal of discovering a closed-form ODE its semantic representation. We assume that the primary objective of discovering closed-form ODE, as opposed to using black-box model, is to have a model representation that can be analyzed by humans to understand the model's behavior (Qian et al., (022). Under this assumption, the specific form of the equation, its syntactic representation, is just a medium that allows one to obtain the description of the model's behavior, its semantic representation, through post-hoe mathematical analysis. We call the process of

---

Published as a conference paper at ICLR 2025
discovering equation and then analyzing it two-step modeling approach. An illustrative example showing the difference between a syntactic and semantic representation of the same ODE (logistic growth model (Verhulst 1845)) can be seen in Figure I
Limitations of the traditional two-step modeling. The traditional two-step modeling pipeline, where an ODE is first discovered and then analyzed to understand its behavior, presents several limitations. The analysis process can be time-consuming, and requiring substantial mathematical expertise. It may even be impossible if the discovered equation is too complex. Furthermore, as the link between syntactic and semantic representation may not be straightforward. modifying the discovered equation to adjust the model's 's behavior may pose significant challenges. This complicates the refinement process and limits the ability to ensure that the model meets specific requirements.
Proposed approach: direct semantic modeling. To overcome these limitations, we propose novel approach, called direct semantic modeling. that shifts away from the traditional two-step pipeline. Instead first discovering a closed-form ODE and then analyzing it, our approach generates the semantic representation of the dynamical system directly from data, eliminating the need for post-hoc mathematical analysis. By working directly with the semantic representation, our method allows for intuitive adjustments and the incorporation of constraints that reflect the system's 's behavior. This direct approach also facilitates more flexible modeling and improved performance, as it does not rely a compact closed-form equation.
Contributions and outline. In [Section 3) of ODEs, discuss the limitations of the traditional two-step modeling pipeline and introduce direct semantic modeling as an alternative. We formalize semantic representation (Section 4) and then use to introduce Semantic ODE in Section 5] a concrete instantiation of our approach for modeling ID systems. Finally, we illustrate its practical usability and flexibility in (Section
define the syntactic and semantic representation
Syntactic representation
Semantic representation
intlcction point (p(re)
0
S.4% < 1.4
1445% kg < 2.8
8g 28
288 <
Mid point
1.4)
x(c)
x(X)
x(O
x(0)
(o(ta).No/2 1.9)
x() x(1) (1
x(t)
2B8
((5)
280
o(xa)
In(2)
0.0
0.0
14
0.0
0.0 2.0
0.0
0.0
Xa
lims (tt))
Iim o(%)
0
Figure 1: Syntactic representation of a logistic growth model refers to its symbolic form, whereas semantic representation describes its behavior for different initial conditions.
2 FORECASTING MODELS AND DISCOVERY OF CLOSED-FORN ODES
In this section, we formulate the task of discovering closed-form ODEs from data and show how it can be reinterpreted as a more general problem of fitting forecasting model.
Let f R+1 R,M, and let T = (lo, +ox). A system of M ODEs is described as
(+t) f(x(() Vt ET
(I)
where æ: -> M is called a trajectory and d.m di We also assume each Fm E CT),), i.e., is twice continuously differentiable on T We denote the dataset of observed trajectories as D
drm. the derivative of am with respect to t.
{(thd) yn ( )n=i)a=1>
(d)
1D where each y&" represents the noisy
measurement of some ground truth trajectory (d) governed by f at time point ((())
A closed-form equation QQan al. 2022 is mathematical expression consisting of a finite number of variables, constants, binary arithmetic operations (+, x, +)) and some well-known functions such exponential or trigonometric functions. A system of ODEs is called closed-fomm when each function f' is closed-form. The task is to find a closed-form f given D.
Traditionally (Bongard & Lipson 2007 Schmidt & Lipson
2009)
discovery of governing equations
has been performed using genetic programming (Koza) 1992). In a seminal paper, Brunton et al.
We assume Im E (T) instead of C² (T),, so that we discuss curvature and inflection points.
2

---

Published as a conference paper at ICLR 2025
(2016a) proposed to represent ODE linear combination of terms from prespecified library. This was followed by numerous extensions, including implicit equations (Kaheman et al. 2020).
equations with control (Brunton et al.) 2016b), and partial differential equations (Rudy et al. 2017) Approaches based on weak formulation of ODEs that allow to circumvent derivative estimation have also been proposed (Msssenger Bortz, 2021a Qian al. 2022). The extended related works
section can be found in Appendix F
Each system of ODEs P defines forecasting mode through the initial value problem (IVP), i.e., for each initial condition (#66) Zo E IRM,, F maps 2o to a trajectory governed by f satisfying this initial condition. Therefore, ODE discovery can be treated as special case of fitting forecasting model F : R M + (T).
3 FROM DISCOVERY AND ANALYSIS TO DIRECT SEMANTIC MODELING
In this section, we define the syntactic and semantic representations, describe the traditional two-step modeling and its limitations, and introduce approach, direct semantic modeling.
3.1 SYNTAX VS. SEMANTICS.
ODEs usually represented symbolically as closed-form equations. For instance, =(t) (tt))((t). We refer to this kind of representation as syntactic.
Syntactic representation of a closed-form ODE refers to its symbolic form, i.c., the arrangement of variables, arithmetic operations, numerical constants, and some well-known functions.
The output of the current ODE discovery algorithms is in the form of syntactic representation. We assume that the primary objective of discovering a closed-form ODE, as opposed to using black-box model, is to have model representation that can be analyzed by humans to understand its behavior. Such understanding is necessary to ensure that the model behaves as expected; for instance, operates within the range of values and exhibits trends consistent with domain knowledge. We call the description of the dynamical system's 's behavior its semantic representation.
Semantic representation describes the behavior of dynamical system. Semantic representation of single trajectory may include its shape, properties, and asymptotic behavior, whereas semantic representation a forecasting model, including a system of ODEs, describes how they change under different conditions, e.g., for different initial conditions.
Comparison between the syntactic and semantic representation of the same ODE is shown in Figure
3.2 TWO-STEP MODELING AND ITS LIMITATIONS
Semantic representation of a dynamical system usually obtained by first discovering an equation (e.g., using ODE discovery algorithm) and then analyzing it. This two-step modeling approach has several limitations (depicted in FFiuure
Analysis of closed-form ODE may be time-consuming, and requiring mathematical expertise. It may be impossible if the discovered equation is too complex. As a result, it may introduce a trade-off between better fitting the data and being simple enough be analyzed by humans.
Obtained insights may be nonactionable. As the link between syntactic and semantic represen- tations often far from trivial, it is difficult to edit the syntactic representation of the model to cause specific change in its semantic representation and to provide feedback to the optimization algorithm solicit a model with different behavior.
Incorporation of prior knowledge. Often the prior knowledge about the dynamical system concerns its semantic representation rather than its syntax. For instance, we may know what shape the trajectory should have (e.g., decreasing and approaching a horizontal asymptote) rather than what kind of terms or arithmetic operations are present in the best-fittingg equation.
3.3 DIRECT SEMANTIC MODELING
To address the limitations of two-step modeling, we propose conceptual shift in modeling low- dimensional dynamical systems. Instead of discovering an equation from data and then analyzing
"Wit some regularity conditions to ensure uniqueness of solutions.
In our work refer forecasting model as any model that outputs trajectory.
3

---

Published as a conference paper at ICLR 2025
to obtain its semantic representation, our approach, direct semantic modeling, generates the semantic representation directly from data, eliminating the need for post-hoc mathematical analysis.
Forecasting model determined by semantic representation A major difference between our approach and traditional two-step modeling is how the model ultimately predicts the values of the trajectory. Given a system of closed-form ODEs f. a forecasting model F is directly given by the equation. We just need to solve the initial value problem (IVP) for the given initial condition. There are plenty of algorithms to do so numerically, the forward Euler method being the simplest (Butcher 2016). In contrast, the result of direct semantic modeling is a semantic predictor Fsom (that corresponds to the semantic representation of the model) that predicts the semantic representation of the trajectory. Then it passes it to a trajectory predictor Fing whose role to find trajectory in a given hypothesis space that has a matching semantic representation. The matching does not need to be unique but Fraj needs to be deterministic. Defining F as Funj o F'sm has multiple advantages. No post- hoc mathematical analysis is required as the semantic representation of F is directly accessed through Fsem- The model can be easily edited to enforce a specific change in the semantic representation because we can directly edit F sem- Incorporating prior knowledge and feedback into the optimization algorithm is also streamlined and more intuitive. Finally, as the resulting model does not need to be further analyzed, it does not need have compact symbolic representation, increasing its fiexibility. Figure 2 compares two-step modeling and direct semantic modeling.
Semantic ODE as a concrete instantiation We have outlined the core principles of direct semantic modeling above. In the following sections, we propose concrete machine learning model that realizes these principles. It is a forecasting model that takes the initial condition 2g E IR and predicts a 1-dimensional trajectory, x : + R. We call Semantic ODE because maps an initial condition to trajectory (like ODEs implicitly do). Although Semantic ODE can only model 1-dimensional trajectories, we believe direct semantic modeling can be successfully applied to multi-dimensional systems. We describe proposed roadmap for future research to achieve that goal in Apppeniix G.2 Before we describe the building blocks of Semantic ODE in Section 5 semantic representation.
we need a formal definition of
step modeling
Direct semantic modeling
Prior knooladgo
Prior knowladge
Difficutt to inject
Challenging edit X
Easy inject
Detaset
Fitting
Predictive model
analysis (may complex)
Datoset
Fitting
Determines
Prodictive modill
P
Somandio nppresonsoriony
Bes.
Cannot provide feedback!
Provide feedbback
Figure 2: Comparison between two-step modeling and direct semantic modeling. Left: The discovery of closed-fone ODEs often allows for human analysis, but editing the equation or providing feedback to the optimization algorithm is challenging. Right: We propose to predict the semantic representation directly from data, which allows for editing the model and steering the optimization algorithm.
4 FORMALIZING SEMANTIC REPRESENTATION
To propose concrete instantiation of direct semantic modeling in Section 5 called Semantic ODE, we need to formalize the definition of semantic representation in [Section 3 make it operational. We consider setting where F= R - (T) is 1D forecasting model (any ODE can be treated as such a model). We first define a semantic representation of a trajectory x E C²(T) and then use it to define a semantic representation of F' itself.
Semantic representation as composition and properties Our definition of semantic representation is motivated by the framework proposed by Kacprzyk et al.) (2024b). Following that work, each trajectory æ can be assigned a composition (denoted c₂) that describes the general shape of the trajectory and the set of properties (denoted pz) which is set of numbers that describes this shape quantitatively. The composition of the trajectory depends on the chosen set of motifs. Each motif describes the shape of the trajectory on particular interval. For instance, "increasing and strictly convex", Given a set of motifs, we can then subdivide T into shorter intervals such that x is described by single motif on each of them. This results in a motif sequence and the shortest such sequence is called a composition. The points between two motifs and on the boundaries are called transition points. An example of trajectory, its composition, and its transition points is shown in Figure 3a
4

---

Published as a conference paper at ICLR 2025
84-6
8 8-40 4
8--b
12
12
A2
1/2
Motifs:
Fe-t incressing concave
So-t decreasing conceva
Sust decressing convex
Se+8 incrussing convex
fend
hocal
Composition:
(S.-b. Su--b S-+ (
Transition points:
5 and
B B-
start
Inflection
Fend
S
maximum
point minimum
bend
(b)
Figure 3: (a) Composition and transition points of ==()= sin(t) on [0, 2nr.. (b) Motif set used in the proposed formalization of semantic representation.
Extending dynamical motifs The set of motifs we choose is inspired by the original set of dynamical motifs (Kaeprzyk et 2024b) but we adjusted and extended it to cover unbounded time domains and different asymptotic behaviors. We define set of ten motifs, four bounded motifs and six unbounded motifs. Each motif is of the form $1.1., i.c. is described by two symbols (each + or -) and one letter (b/ b//u/h). hh. The symbols refer to its first and second derivatives. The letter b signifies the motif is for bounded time domains (e.g., for interval (t,, t2)). Both h and u refer to unbounded time domains. These motifs are always the last motif of the composition, describing the shape on (lend: + where fend is the (-coordinate of the last transition point. / specifically describes motifs with horizontal asymptotes. For instance, +h is an unbounded motif that describes a function that is decreasing (-), strictly convex (+) and with horizontal asymptote (h). All motifs are visualized in Figure 3b. Note that we excluded the three original motifs describing straight lines to simplify the modeling process. If necessary, they can be approximated by other motifs with infinitesimal curvature. We denote the set of all compositions constructed from these motifs as C.
Properties Apart from the composition, the semantic representation of a trajectory also involves a set of properties. Ideally, the properties should be sufficient to visualize what each of the motifs looks like and to constrain the space of trajectories with the corresponding semantic representation. Following the original work, we include the coordinates of the transition points as they characterize bounded motifs well. In contrast to their bounded counterparts, the unbounded motifs are not described by their right transition point but by a set of motif properties. These, in turn, depend on how we describe the unbounded motif. For instance, we could parameterize 84++ as sr()) ar(eea) ((( B, where (fend: :ffend)) is the position of the last transition point. In that setting, B is the property of 8444 that describes the doubling time of x (z(t + B) = (x)t).. ). In reality, choosing good parametrization with meaningful properties is challenging, and we discuss it in more detail in Appenddix D.2 The set of properties also includes the first derivative at the first transition point (tv) and the first and the second derivative at the last transition point (fend). They are needed for the trajectory predictor described in /Section 5.2] Each composition E C may require different set of properties that we denote Pe. For instance, a trajectory x with C₂ = (8++b, 84- will have Pz = (fa, f1, (ffo), (((())))((()))))))))))) ). (((o), .it(t)). #(tz), ((tz),, h, (122), where each (1z, ((:)) is a transition point, and (h, (1/2) are the properties of the unbounded motif (see Figure 4). We denote all possible sets of properties as P. where P = Ucec P...
We are finally ready provide formal definition of the semantic representation of trajectory x E (T) and a forecasting model F : R - C² (T). Given this formal definition of semantic representation, we introduce our model, Semantic ODE, in the next section.
Definition 1. The semantic representation of a trajectory x E (²) is a pair (cz, Px), where Ca E C is the composition of z and Pz E Pec, is the set of properties as specified by Cz.
Definition 2. The semantic representation of F : R (T) is pair (CF, (F): R + C x P defined as follows. CF R + C is called a composition map and it maps any initial condition 2g E R to composition of the trajectory determined by its initial condition. Formally, CF (zo) CF(z)) Pr : RR + P is called a property map, and it maps any initial condition To E R to the properties of the predicted trajectory Cr(z Formally, Pr( (=o))
5

---

Published as a conference paper at ICLR 2025
5 SEMANTIC ODE
ODE discovery methods aim discover the governing ODE f and thus forecasting model F (defined by IVP), which is later analyzed to infer its semantic representation (Cp, Pp) in a two-step modeling process. In this section, introduce novel forecasting model, called Semantic ODE, that allows for direct semantic modeling. As described in [Section 3.3 our model F consists of two submodels, Fsem and Finj (where F Funj o F. We can now define formally Fsem function that maps an initial condition E I to the semantic representation of a trajectory, i.c., Fsem R. such representation, i.e., Fini : C.
P. Funj then takes this semantic representation and matches it trajectory with
(T) such that if = [eij p) then (Ca,Pz) = (c.p). This is visualized in Figure 4] Crucially, by definition, the semantic predictor is the semantic representation of F. Indeed, let F(xu) Fugi (Fse (xo)). Then by the definition of Fug above, (Cp, (zo) = (F(a))))(()))))))))))))))))))))) (Cz. = (em) Unlike the two-step modeling approach, there no need for post-hoc mathematical analysis. The semantic representation of the model can be directly inspected through the semantic predictor. In [Section 5.1] we propose an implementation for we describe the trajectory predictor.
the semantic predictor, and in Section 5.2,
Initial condition
Semantic predictor
Semantic representation of
Trsjectory predictor
Trajectory
Composition:
(f-)))
x(0) 0.4
Frem Section 5.1
Transition points:
(18.44))
Firaj
Derivatives:
x(G) 0.3.k(G,) 0.7.X(t.)
Section 5.2
(18.144
(00,0.0)
Properties of
h 28, Eyp 2.9
Figure 4: Semantic ODE (F) maps the initial condition to the semantic representation of the trajectory (using Fsem and then uses it to predict the actual trajectory (through Finj). Formally, F == Fimj Fsem-
5.1 SEMANTIC PREDICTOR
Semantic predictor F'sm consists of two models. One that predicts the composition denoted Foom (corresponding the composition map CF), and one that predicts the properties, denoted F peop (corresponding to the property map Pr).
Foom classification algorithm from to C' C C, where C' our chosen composition library. We model it as a partition of IR into intervals (here called branches), each mapped to different composition. The maximum number of branches I E N selected by the user. For instance, the logistic growth model example in Figure 1 would have 3 branches: Feom (zo) = (8++ 8+-1) for S Zo < 1.4, From (zp) = (8+ for 1.4 S xo S 2.8, and Feom (2o) = (s_+h) for 2.8 < Z0.
As mentioned earlier, although zo 2.8 should be straight line, approximate it with 84-
Prop is modeled as a set univariate functions. Each of them describes one of the coordinates of the transition point, the value of the first or second derivative, or the properties of the unbounded motif. These functions are different for different compositions. Continuing our logistic growth example, Fprop is a piecewise function consisting of three composition-speciin property sub- sub-maps P/pp p($=+>,8--b, on +0 that describes the properties of ($4+4) -4) The list of properties includes the coordinates of both transition points, first and second derivatives, and two properties of the unbounded motif (8) This is visualized in Figure 5
ppppp- p F that correspond to the three branches described above. Let us focus
Property map composition 31 Trarokonn puints teccodivates)
Trans tios points oionnss
Derivatives
Properris of unborrded mott
0.6-
atuol
555
125
0-
0.
3.50
8.75
110
325
323
050
375
3.00
123
0.5
355
1.0
18
325
0.50
1.00
1.5
Figure 5: Property sub-map of the logistic growth model (for composition (84+6, 84-))) describes how various properties depend on the initial condition. See Figure 9 for full semantic representation.

---

Published as a conference paper at ICLR 2025
Training of Fsem performed in two steps. First, train the composition map From- Then, the dataset divided into separate subsets, each for a different compostion-accordingg to the composition map and separate property sub-map trained each of the subsets. A block diagram presented in Figure 6 and the pseudocode of the training procedure can be found in Appendix
each (He, f. For each c E C: return mio P2 (=,)) 1
Dynamic
programming
Room
Feom
D2
property submap
Pirop
Divide according to Fooms
Fooap
D2
Fit property aubmep
Krop
Y, F2 (44"),
(Haatt" Mautttl" -apgettt y
Figure 6: Block diagram showing main elements of training the semantic predictor.
5.2 TRAJECTORY PREDICTOR
The goal of the trajectory predictor (Fing) map the semantic representation (c, to a trajectory E X' such that (Cz, Px) (c. In that case, say that conforms to (c,p) and write it as x (c,p). This mapping requires specifying the hypothesis space X' for the predicted trajectories. As mentioned previously, we would also like each trajectory to be twice continuously differentiable, ie., X C C We choose X' to be a set of piecewise functions defined as
ar (1)
1 85 lend
x(t)
2
(1z(t)
8 lend
where lend is the last transition point (before the unbounded motif), 2 E X called the bounded part of x. and | E 2 is called the unbounded part of x. Finding z and jæ: is done separately and we discuss respectively in Sections 5.2.I and 5.2.2
5.2.1 BOUNDED PART OF THE TRAJECTORY
In this section, define X and describe how find the bounded part of the trajectory (17) given semantic representation (c, p) without the unbounded motif and its properties, denoted (().
Cubic splines We decide to define X as a set of cubic splines. They are piecewise functions where each piece is defined as cubic polynomial. The places where two cubics are joined are called knots. Cubic splines require that the first and second derivatives at the knots be the same for neighboring cubies so that the cubic spline guaranteed to be twice continuously differentiable. Cubic splines are promising because they are flexible, and for fixed set of knots, the equations for their values and derivatives are linear in their parameters. We come up with two different ways finding a cubic spline that conforms to a particular semantic representation that leads us to develop two trajectory predictors: Fi C x P C and F
C2. We use F00 during training of Fsem as is fast and differentiable, but the found trajectory may not be in C2 (only continuous). At inference, we use Pily that is slower and differentiable but ensures that the trajectory is in C2, We describe both approaches briefly below. More details available in Appendix D.I
C trajectory predictor FO describes each motif as a single cubic. Each cubic is found by solving a set four linear equations. Two are for the positions of the transition points, and two are for the derivatives, one at each transition point. As each transition point (apart from to) is either local extremum or an inflection point, we know that either the first or the second derivative vanishes. The first derivative at to is specified directly in pl. During training, also make that this derivative always in the range described in Table 9] We prove why this is sufficient to guarantee that the predicted trajectory conforms (c) 22) in Appendix D.1.2
7

---

Published as a conference paper at ICLR 2025
Ca trajectory predictor Fig describes each motif two cubics with an additional knot between every two transition points. Given K cubics (and K + 1 knots), the traditional way to ensure that xj is in C² is to set up 3(K 1) constraints that match the values and both derivatives at the knots. Instead, we propose to fit the second derivative of the cubic spline (&)), a piecewise linear function, and then integrate it twice to get the desired cubic spline. As it is continuous, values at the knots (2), at knot tg). Together with the additional two parameters for the integration constants, not only reduce the number of parameters K + 3 (while still guaranteeing the trajectory is in C2) but, more importantly, we can control exactly the value of the second derivative. Importantly, as shown in Appenddix D.I] we can impose v2 eR to make sure that the first and second derivatives at the transition points are correct. Then, we optimize the additional knots and some of 2/ to minimize the error between the predicted and target r r-coordinates of the transition points (with an additional loss term for the signs of the first derivatives at the knots), We minimize this objective using L-BFGS-B (Liu & Nocedal, 1989) and Powell's method (Powell 1964) until we find solution where the error on the transition points is smaller than a user-defined threshold. it does not succeed, then we default to F 'inj"
describe it solely by its
Uk < O accordingly. We set some
5.2.2 UNBOUNDED PART OF THE TRAJECTORY
In this section, we define X' and describe how we can find the unbounded part of the trajectory given an unbounded motif and its properties, as well the coordinates of the last transition point and both derivatives at this point. We need them to ensure that x is twice differentiable at lend- First, we need to choose which properties we are interested in. They should describe the shape (and long-term behavior) of the unbounded motif in sufficient detail such that we do not need to see the underlying equation to visualize it. For instance, for motifs with horizontal asymptotes (s+-h and s-+)), we choose these to be h and f1/2 where = h is the horizontal asymptote and 11/2 is the time where the trajectory is in the middle between the last transition point and the h, i.e., atti/2) (a=(fem)) /2. In exponential decay, 41/2 fend would correspond to "half-life". Then, we need to come with a parametrization that would both guarantee the shape of the trajectory and allow to impose any possible properties (e.g.. h in s--h needs to satisfy h < xx(ena)). We parametrize 8-+h as st) 2(x((enn)) h)/(1 e + h where g is an appropriately defined cubic spline. We show how we can find such and why it has the desired properties. as well as properties and parameterizations for other motif', in Apppendix D.2] Importantly, it does not matter how complicated these parameterizations are, as they are not used by humans to understand the model. All information is directly available in F Fsem-
6 SEMANTIC ODE IN ACTION
In this section, we want to illustrate the usability of Semantic ODFs and highlight their advantages: semantic inductive biases (Section 6. 6. comprehensibility (Section 6.2), editing (Section 6.3) and flcxibility and robustness to noise (Section 6.4). To demonstrate the first three advantages, a case study of finding a pharmacokinetic model from a dataset of observed drug concentration curves. Such models are essential for drug development and later clinical practice. Details about the dataset can be found in Appenddix E.I. We will contrast our approach with one of the most popular ODE discovery al algorithms, SINDy (Brunton et al. 2016b). However, many of the observations will apply to other algorithms as well.
present
6.1 SEMANTIC INDUCTIVE BIASES
In Semantic ODEs, user can specify inductive biases about the semantic representation of the model. This is in contrast to the syntactic inductive biases, available in ODE discovery algorithms. Semantic inductive biases can be more meaningful and intuitive for users than syntactic ones as the relationship between the syntactic and semantic representation may be non-trivial. The role of syntactic biases is, predominantly, to ensure that the equations can be analyzed by humans. They not designed to accommodate prior knowledge about the system's 's behavior. Examples of inductive biases in SINDY and Semantic ODE are shown in Table I
The drug concentration curve describes the drug plasma concentration after administration as a function of time. Without any additional doses, we would expect the concentration decay 0
8

---

Published as a conference paper at ICLR 2025
Table 1: Examples of inductive biases in SINDy and Semantic ODE.
Syntactic inductive biases in SINDY
Semantic inductive biases in Semantic ODE
Autonomous system: whether the governing ODE system is time-invariant.
Library compasitions: The maximum number of motifs, the starting motif, the type of asymplotic behavior. The complexity of composition map: how many different compositions it should predict.
Library functions: whether include, e-g-. polynomials, irigonometric terms, expomentalllogaiimmiic terms, cross-termss Sparsity; the number of terms, strength of penulty terms such as 1.1 or 1.2 norms.
Complexity of the property maps: how many trend changes the property maps could have.
as t - oc. This is semantic inductive bias that can be easily put into Semantic ODE. We can enforce the last motif to always be + (decreasing, convex function with a horizontal asymptote) by removing all biologically impossible compositions from the library C'. Designing syntactic inductive biases based on the prior knowledge is more challenging. SINDy assumes á >i=1 aggSs, i.e., f is a linear combination of pre-specified functions. As choosing which terms should and should not appear in the equation is far from obvious, we choose a general library containing polynomial terms, exp, and trigonometric terms (see Appendix E.2 for details). We also consider different sparsity levels, from 1 to 15 non-zero terms.
6.2 COMPREHENSIBIL.IT
We have fitted Semantic ODE with the inductive bias described earlier and versions of SINDY with different sparsity constraints. The results can be seen in Table 2
Table 2: Results of fitting Semantic ODE and SINDy to the pharmacokinetic dataset.
Model
Syntactic biases
Semantic bisses
Symtactc representation
Semantic representation
In-demain
Out-domain
(e.
SINDY
Ougi, n S
NA
3.0 (t)
NA
0.222
0.024:o.05)
SINDY
5.56
(f)
NA
0.112
0.064
SINDY
2
OaBi, n S 5
NA
Equation (
NA
0.101 (0
16. 16.850
SINDY
2-- aegi.
NA
Equalton
NA
m.029
18 (.0)
SINDY
Zers n n 5
Equation
0.022l
(1.577
Semantic ODE
NA
leg IS 4, de ends with s..h
NA
Figure 7
0.016
o 0.033, 2003
Semantic ODE
NA
(A4-b, Sund, 9--b),
Pigure
0.018 (008) 0.015
We can see that Semantic ODE better fits the dataset than even the longest equations found by SINDy, Compact equations, e.g., .tt)= = 5.56 43. (tx( have poor performance. To improve it, we need to allow for much more complicated equations, such as Equation (9) in Appendix B, Such equations are very hard to analyze and, as a result, may be no more interpretable than black box model. They are also more prone to over-fitting, demonstrated by large out-domain error. More importantly, Semantic ODE can be directly understood by looking at its semantic representation in Figure 7
Composition map Fa S- -*4x
Property map for composition (5)
peints
Properiue
C35-
EXE
491
Figure 7: Semantic representation of Semantic ODE fitted to the pharmacological dataset.
The left side of Figure 7 tells us that for all initial conditions, the shape of the predicted trajectory is going to be the same, namely (84-6; 8_-,, As the composition map has only one branch, we have only one property map, and it is visualized on the right of the composition map. It consists of four subplots. Going from left to right, we have the t-coordinates of the transition points, x- coordinates of the transition points, values of the derivatives at the boundary transition points, and the properties of the unbounded motif. In our case, the unbounded motif is a convex, decreasing function with a horizontal asymptote, and it is described by the value of the horizontal asymptote and its "half-life", which is the í-coordinate of the point where the value of r is in the middle between the last transition point and the asymptote. By looking at this representation, we can readily see how the trajectory changes with respect to the initial condition. We see that the t-coordinates of the transition points remain fairly constant, whereas z-coordinates increase linearly. In particular, can see how the maximum of the trajectory (z(ti)) increases linearly from 0.6 up to 1.3. Arriving at similar observations about the discovered ODEs requires significant time and expertise. In particular,

---

Published as a conference paper at ICLR 2025
we do not know what the composition map of the discovered ODEs looks like. Thus, we cannot be sure whether the model behaves correctly.
6.3 EDITING
Often, the fitted model does not satisfy all our requirements. A very important requirement in pharmacology would be to make sure that the model is biologically possible. In particular, the predicted concentration should decay to as t - co. The horizontal asymptote of our model is close to O but not exactly. That is why the extrapolation error for ł > 1 is substantially higher than for t & 1. Fortunately, we can edit the property map directly. We can impose the value of the horizontal asymptote to be equal to 0 and retrain the model. The new property map can be seen in Figure 8 Importantly, as shown in the last row of Table 2 the extrapolation error dropped to levels comparable with the in-domain error. In two-step modeling, it may be challenging verify and impose such requirements for the predicted equations, especially the longer ones. As a result, we may end up with model that does not obey crucial domain-specific rules.
6.4 FLEXIBILITY AND ROBUSTNESS TO NOISE
As Semantic ODE does not assume that the trajectory is governed by closed-form ODE, it may fit systems that are not described by those. In particular, we show in Table 3 how, beyond standard ODEs (logistic growth model), it can fit systems governed by general differential equation =tt) f f(x(t), .t., where f does not have a compact closed-form representation, a multidimensional ODE where only one dimension is observed (pharmacokinetic model), delay differential equation (Mackey equation (integro-DE We compare our approach to SINDy Brunton et al. 2016b) and WSINDy (Reinbold et al.. 2020; Messenger & Bortz 2002aa) as implemented in Pys INDD library (de Silva et al.) [020] Kaptanoglu et al. (022). We consider variants constrained to terms (in a linear combination) ensure compactness (SINDy-5, WSINDy-
& Glass 1977), and an antegro-different
5) and where the number of terms is fine-tuned and may no longer be compact (SINDy, WSINDy).
We also include standard symbolic regression method, PySR (Cranmer 2020). adapted for ODE discovery and constrained to 20 symbols. We also compare with three black box approaches: Neural
ODE (Chen et al.] 2018), Neural Laplace (Holt et al.] 2022), and DeepONet (Lu et al.] 2020). Semantic ODE is more equally robust to noise and performs better than the methods constrained to compact equations in most cases. Moreover, its performance could possibly be further improved by incorporating semantic inductive biases and model editing discussed earlier. Additional experiments can be found in Apppndixx B) Details on experiments are available in Appendii E
Table 3: Comparison of Average RMSE obtained by different models. Average performance over 5 random seeds and data splits is shown with standard deviations in the brackets.
Logistic Growth
General
Pharmacokmmetic
Mackey Glass (DDE)
Iniegro-Dii
Method
noise
high noise
low nosse
high noise
low noise
hugh noise
noise
high noise
low notse
high noise
SINDs-5
0.012,
0.222
0.053, (15)
0.103:
0.089,co4)
0.230, mm
0.238
0.248,.05)
0.451
0,268:
WSINDy-5
0.030
0.222 are)
0.008;
[11102 mm8)
0.211:
01415)
1.2722 az3)
0.300l (e1)
0.150)
11.452: 315)
PySR 20
0.012; awm
0.224
0.078
0.115
0.053,
0.242,
0.261
0.288
0.027
0.383:
SINDy
0.218
0.068
0.115 m12)
0.020).co))
(.200)
0.2522
0:257,100)
0.315
0.248
WSINDY
0.020).com
0,217 ame)
0,062 (000
699)
0.688 ree)
(1223) (n))
0.20K) ass)
0.207.081)
0.132 (091
0.3000
Noural ODE
0.023, umss
0.197 (wss
0.029) 1065)
0.075
0.036,cce)
.um))
0.177
0.194|,000
0.073 .xT7
0.215,a09:
Newal Laplace
0.126 uoo
0,230),
0.108 year)
0.138 (124)
0.100).co2)
0.22%) un3
0.057 (03)
0.094 1e8)
0.075 msa)
0.2440.a14)
DeepONet
0.184
0.306
0.160 noz
0.105
0.058 c20)
0.212, (m))
0.107
0.132
0.100] m1s)
0.250/aa44
Semonic OnE
0.015
0.198
0.016.833:
0.068,
0.015
0.197
0.037
0.025,.0o 0.077
0.204,
7 DISCUSSION
Limitations and open challenges As Semantic ODE is the first model that allows for direct semantic modeling, we focused solely on 1-dimensional systems and we describe a possible roadmap to nigher-dimensions
settings in AAppendix G.2 The current definition of semantic representation assumes that the trajectory has finite composition, i.e., cannot have an oscillatory behavior (like sin). Of course, we could fit a periodic function on any bounded interval, but it would fail to predict the oscillatory behavior outside of it. We discuss more limitations in Appendix G.I
Direct semantic modeling as new way for modeling dynamical systems In this work we outlined the main principles of direct semantic modeling, discussed its advantages, and illustrated how it can be achieved in practice through Semantic ODE. We believe this approach can transform the way dynamical systems are modeled by shifting the focus from the equations to the system's 's behavior, making the models not only more understandable but also more flexible than other techniques.
10

---

Published as a conference paper at ICLR 2025
Ethics statement In this paper, we introduce novel approach for enhancing the comprehensibility of a dynamical system's 's behavior through direct semantic modeling, with a practical implementation called Semantic ODE. Improved tramsparency of machine learning models is crucial for tasks such as model debugging, ensuring compliance with domain-specific constraints, and addressing potential harmful biases. However, such techniques, if misused, can lead false sense of security in model decisions or be leveraged merely for superficial regulatory compliance. As our approach is applicable high-stakes domains like medicine and pharmacology, it is vital conduct thorough evaluation before deploying the model in such contexts. This evaluation must ensure that the model's behavior aligns with ethical considerations and does not support decisions that could negatively impact individuals' health and well-being.
Reproducibility statement All mathematical definitions provided in Secciios 4]and Jand Ap pendix A.2 The proofs are provided in Appendix D The implementation, including a block diagram and pseudocode, is discussed in Section 5 and Appendices C and [ The experimen- tal details are discussed in |Section and Appendix https://github. http:s//gggtubb...ooor//2/22222
All experimental code available at www.aaapryykssmattaaaaaa
This work was supported by Roche. We would like to thank Max Ruiz Luyten, Harry Amad, Julianna Piskorz, Andrew Rashbass, and anonymous reviewers for their useful comments and feedhack on earlier versions of this work.
Acknowledgments
REFERENCES
Takuya Akiba, Shotaro Sano, Toshihiko Yanase, Takeru Ohta, and Masanori Koyama. Optuna: next-generation hyperparameter optimization framework. In Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, 2019.
Henrik André-Jönsson and Dushan Z Badal. Using signature files for querying time-series data. In Principles of Data Mining and Knowledge Discovery: First European Symposium, PKDD '97 Trondheim, Norway, June 24-27, 1997 Proceedings 1, pp. 211-220. Springer, 1997.
Dimitris Bertsimas and Wes Gurnee. Learning sparse nonlinear dynamics via mixed- -ixeed-nteger opti- mization. Nonlinear Dynamics, 11 177))66855-660)) 2023.
L. Biggio, T. Bendinelli*, A. Neitz, A. Lucchi, and G. Parascandolo. Neural Symbolic Regression that Scales. In 38th International Conference on Machine Learning, July 2021.
Bongard and H. Lipson. Automated reverse engineering of nonlinear dynamical systems. Pro- ceedings of the National Academy of Sciences, 104(24):9943-9948, June 2007. ISSN 0027-8424, 1091-6490. doi: 10. (0.1073/ppaa..0099476614
Murray Bourne.
Solving
and Simultaneous Differential Equations.
2018
Steven L. Brunton and Nathan Kutz. Data-Driven Science and Engineering: Machine Learning, Dynamical Systems, and Control. Cambridge University Press, May 2022. ISBN 978-1-00-911563-
Steven L. Brunton, Joshua Proctor, and Nathan Kutz. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the National Academy of Sciences, 13(11))39322393) 3937, April 2016a. ISSN 0027-8424, 1091-6490. doi: 10. 1073/0nas.. 1517384113.
Steven L. Brunton, Joshua L. Proctor, and J. Nathan Kutz. Sparse Identification of Nonlinear Dynamics with Control (SINDYe). PFAC-PapersOnLine, 491871770-715, January 2016b. ISSN 2405-8963. doi: (0.10166..facoo.2216..22222
C. Butcher. Numerical Methods for Ordinary Differential Equations. John Wiley & Sons, August 2016. ISBN 19-1500-333
11

---

Published as a conference paper at ICLR 2025
Ricky T. Q. Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural Ordinary Differential Equations. In Advances in Neural Information Processing Systems, volume 31. Curran Associates, Inc., 2018.
J. T. Y. Cheung and G. Stephamopoulos. Representation of process trends Part I. A formal rep- resentation framework. Computers & Chemical Engineering, 14(4))495-510, May 1990. ISSN 0098-1354, doi: 110006000988-355490)888777
Miles Cranmer. PySR: Fast & parallelized symbolic regression in Python/Julia. Zenodo, September 2020.
Stéphane D' Ascoli, Sherre-Alexadrrt Kamienny, Guillaume Lample, and Francois Charton. Deep symbolic regression for recurrence prediction. In Proceedings the 39th International Conference on Machine Learning, pp.
PMLR, June 2022.
Brian de Silva, Kathleen Champion, Markus Quade, Lean-Christophe Loiseau, J. Kutz, and Steven Brunton. PySINDy: A Python package for the sparse identification of nonlinear dynamical systems from data. Journal of Open Source Software, 5149):2104, 2020. doi: 10.21 102//00/s/sss022 104.
N. Fritsch and R. E. Carlson. Monotone Piecewise Cubic Interpolation. SIAM Journal on Numerical Analysis, 17(2):23-2244 1980. ISSN 0036-1422
Yi Rang Han, Ping I. Lee, and K. Sandy Pang. Finding Tmax and Cmax in Multicompartmental Models. Drug Metabolism and Disposition, 46(11):1796-188)) November 2018. ISSN 0090-9556, 1521-009X. doi: 12.112//mm.111 18.082636.
Trevor Hastie and Robert Tibshirani. Generalized additive models. Statistical Science, 13))2997-318, 1986.
Samuel Holt, Zhaozhi Qian, and Mihaela van der Schaar. Neural Laplace: Learning diverse classes of differential equations in the Laplace domain. In International Conference on Machine Learning, pp. 8811-8833 PMLR, 2022.
Krzysztof Kacprzyk and Mihaela van der Schaar. Shape Arithmetic Expressions: Advancing Scientific Discovery Beyond Closed-form Equations. In Proceedings of The 27th International Conference on Artificial Intelligence and Statisties. PMLR, 2024.
Krzysztof Kacprzyk and Mihaela van der Schaar. Beyond Size-Based Metrics: Specific Complexity in Symbolic Regression. In The 28th International Conference Artificial Intelligence and Statistics, 2025.
Measuring Task-
Krzysztof Kacprzyk, Zhaozhi Qian, and Mihaela van der Schaar. D-CIPHER: Discovery of Closed- form Partial Differential Equations. In Advances in Neural Information Processing Systems, volume 36, pp. 27609-27644, December 2023.
Krzysztof Kacprzyk, Samuel Holt, Jeroen Berrevoets, Zhaozhi Qian, and Mihaela van der Schaar. ODE Discovery for Longitudinal Heterogeneous Treatment Effects Inference. In The Twelfth International Conference on Learning Representations, 2024a.
Krzysztof Kacprzyk, Tennison Liu, and Mihaela van der Schaar. Towards Transparent Time Series Forecasting, In The Twelfth International Conference on Learning Representations, 2024b.
Kadierdan Kaheman, Nathan Kutz, and Steven L. Brumton. SINDy-PI: A robust algorithm for
parallel implicit sparse identification of nonlinear dynamics. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 476(2242):0000229, October 2020. doi: 10. 00.1098/rsaa.220..7779
Alan A. Kaptanoglu, Brian M. de Silva, Urban Fasel, Kadierdan Kaheman, Andy J. Goldschmidt, Jared Callaham, Charles B. Delahunt, Zachary G. Nicolaou, Kathleen Champion, Cean-Chhittophe Loiseau, J. Nathan Kutz, and Steven L. Brunton. PySINDy: A comprehensive Python package for robust sparse system identification. Journal of Open Source Software, ((69):3994, 2022. doi: 10.21 1 (o5/joss.0399.
Diederik P. Kingma and Jimmy Ba. Adam: A Method for Stochastic Optimization, January 2017.
12

---

Published as a conference paper at ICLR 2025
John R. Koza, Genetic Programming: On the Programming Computers by Means of Natural Selection. Complex Adaptive Systems. MIT Press, Cambridge, Mass, 1992. ISBN 978-0-262- 11170-6.
Changhee Lee, Jinsung Yoon, and Mihaela van der Schaar. Dynamic-DeepHli: A Deep Learning Approach for Dynamic Survival Analysis With Competing Risks Based on Longitudinal Data. IEEE Transactions Biomedical Engineering, 67(1):122-133, January 2020. ISSN 1558-2531. doi: 10.1109/BME.22019
Dong Liu and Jorge Nocedal. On the limited memory BFGS method for large scale optimization. Mathematical programming. 45(1):503-528 1989.
JLEKS Lonardi and Pranav Patel. Finding motifs in time series. In Proc. of the 2nd Workshop Temporal Data Mining, pp. 53-68, 2002.
Lu Lu, Pengzhan Jin, and George Em Karniadakis. DeepONet: Learning nonlinear operators for identifying differential equations based on the universal approximation theorem of operators, April 2020.
Temry Lyons. Rough paths, Signatures and the modelling of functions on streams, May 2014.
Michael C Mackey and Leon Glass. Oscillation and chaos in physiological control systems. Science, 197(4300):277288, 1977.
Georg S Martius and Christoph Lampert. Extrapolation and learning equations. Sth International Conference Learning Representations, ICLR 2017- Workshop Track Proceedings, 2017.
Daniel A. Messenger and David M. Bortz. Weak SINDy: Galerkin- Selection. Multiscale Modeling & Simulation, 1933) January 2021a. ISSN 1540-3459, 1540-3467. doi: 10.1
Model
Daniel A. Messenger and David M. Bortz. Weak SINDy for partial differential equations. Journal of Computational Physics, 443:110525, October 2021b. ISSN 00219991. doi: 10.1... 110525.
Dr Mould and Rn Upton. Basic Concepts in Population Modeling, Simulation, and Model-Based Drug Development. CPT: Pharmacometrics & Systems Pharmacology, 1(9):6, 2012. ISSN 2163-8306. doi: 10.1008/psp.2112..4
Emre O. Neftci and Bruno B. Averbeck. Reinforcement learning in artificial and biological systems. Nature Machine Intelligence, 1(3):133-443, March 2019. ISSN 2522-5839. doi: 10.1038/
$42256-019-0025-4
Brenden K Petersen, Mikel Landajuela Larma, T Nathan Mundhenk, Claudio P Santiago, Soo K Kim, and Joanne Kim. Deep Symbolic Regression: Recovering Mathematical Expressions From Data via Risk-seeking Policy Gradients. In ICLR 2021, 2021.
Michael JD Powell. An efficient method for finding the minimum of a function of several variables without calculating derivatives. The computer journal, 7(2):155-162, 1964.
Steven Pruess. Shape preserving C cubic spline interpolation. IMA Journal of Numerical Analysis, 13(4):493-507, 1993.
Zhaozhi Qian, Krzysztof Kacprzyk, and Mihaela van der Schaar. D-CODE: Discovering Closed-form ODEs from Observed Trajectories. The Tenth International Conference Learning Representa- tions, 2022.
Markus Quade, Markus Abel, Kamran Shafi, Robert K. Niven, and Bernd R. Noack. Prediction of dynamical systems by symbolic regression. Physical Review F,, 94(1):012214, July 2016. ISSN 2470-0045, 2470-0053. doi: 10.1 103/Phyyeee
M. Raissi, Perdikaris, and G.E. Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378:686-707, February 2019. ISSN 00219991. doi: 10.1016//.jcp.201 10.045.
13

---

Published as a conference paper at ICLR 2025
Maziar Raissi and George Em Karniadakis. Hidden physics models: Machine learning of nonlinear partial differential equations. Journal of Computational Physics, 357:125-141, March 2018. ISSN 0021-9991. doi: 10.17..9,,
Patrick A. K. Reinbold, Daniel R. Gurevich, and Roman O. Grigoriev. Using noisy or incomplete data to discover models of spatiotemporal dynamics. Physical Review E, 101(1):010203, January 2020. ISSN 2470-0045, 2470-0053. doi: 10.1 S15S/PhysRevE.101.000222
Samuel H. Rudy, Steven Brumton, Joshua L. Proctor, and Nathan Kutz. Data-driven discovery of partial differential equations. Science Advances, 3(4):e April 2017. ISSN 2375-2548. doi: 10. 2e1112//caid...66 1602661......161.... 1602614.
Subham Sahoo, Christoph Lampert, and Georg Martius. Learning Equations for Extrapolation and Control. In Proceedings the 35th International Conference Machine Learning, pp. 4442 4450. PMLR, July 2018.
Michael Schmidt and Hod Lipson. Distilling Free-Form Natural Laws from Experimental Data. Science, 324(5923):81 April 2009. ISSN 0036-8075, 1095-9203. doi: 10. 166ssciece. 1165893.
Trevor Stephens. Gplearn: Genetic programming in python, with scikit-learn inspired and compati- ble api, 2022.
Sahar Torkamani and Volker Lohweg. Survey on time series motif discovery. WIREs Data Mining and Knowledge Discovery, 7(2):e1 112)=1119, 199, 2017. ISSN 1942-4795. doi: 10. 10.1002/widm.1 199.
Silviu-Marian Udrescu and Max Tegmark. AI Feynman: A physics-inspired method for symbolic regression. Science Advances, 6 6(16):euay6631, April 2020. ISSN 2375-2548. doi: 10.1 126ssciadv. aay2631.
Silviu-Marian Udrescu, Andrew Tan, Jiahai Feng, Orisvaldo Neto, Tailin Wu, and Max Tegmark. AI Feynman 2.0: Pareto-optimal symbolic regression exploiting graph modularity. 34th Conference on Neural Information Processing Systems (NeurIPS 2020), 2021.
Pierre François Verhulst. Recherches mathématiques sur la loi d' d'accroissement de la population. Hayez, 1845.
Julia Wilkerson, Kald Abdallah, Charles Hugh-Jones, Greg Curt, Mace Rothenberg, Ronit Simantov, Martin Murphy, Joseph Morrell, Joel Beetsch, Daniel J Sargent, Howard Scher, Peter Lebowitz, Richard Simon, Wilfred Stein, Susan Bates, and Tito Fojo. Estimation of tumour regression and growth rates during treatment in patients with advanced prostate cancer: A retrospective analysis. The Lancet Oncology, 18(1):143- 154, January 2017. ISSN 1470-2045. doi: 10.1016/ (1470-2045(16)30633
-ean-Baptiste Woillard, Brenda C. M. de Winter, Nassim Kamar, Pierre Marquet, Lionel Rostaing, and Annick Rousseau. Population pharmacokinetii model and Bayesian estimator for two tacrolimus formulations twice daily Prograf and once daily Advagraf. British Journal of Clinical Pharmacol- ogy, 71(3):391-40) March 2011. ISSN 1365-2125. doi:
Lexiang Ye and Eamonn Keogh. Time series shapelets: A new primitive for data mining. In Proceedings of the 15th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pp. 947-956, Paris France, June 2009. ACM. ISBN 978-1-60558-495-4 doi: 10. 0t.114//155701.15577222
Byoung-Kee Kee Yi and Christos Faloutsos. Fast time sequence indexing for arbitrary Lp norms. 2000.
14

---

Published as a conference paper at ICLR 2025
TABLE OF SUPPLEMENTARY MATERIALS
1.
Appendix A
notation and definitions
2.
Appendix
additional results
3.
Appendix
training the semantic predictor
4.
Appendix
trajectory predictor
5.
Appendix
experimental details
6.
Appendix
extended related works
7.
Appendix
additional discussion
A
NOTATION AND DEFINITIONS
A.1
NOTATION
Symbols used throughout this work can be found in Tables and
A.2 DEFINITIONS
In this section, we provide formal definition of some of the terms introduced in the main text.
From the work by Kacprzyk et al. (2024b).
Definition Let be set of intervals on I and be the set interval functions, i.e., real functions defined on intervals. A motif is binary relation between the set of interval functions F and the set of intervals I (i.c., s C F x Z). We denote (ø, i) E s as Φ ~ s and read as "ø on i has motif s". Each motif needs to be
well-defined, i.e., for any & E F, and any E I,
Φ ~ = i C do
(3)
translationnimvvriann, i.e., for any E I, and any F.
Φli siiss ~ do g) ( q) 4q
(4)
In the next definition, formally define the motifs Definition 4. Let C² and let 71.72 E R such that T1 12.
define semantic representation.
2|[r,,12] 8444 if E (f1,12) a(t) 0, iif) > 0
x|[r1,72] 84-6 Vt (21,1) b(t) 0, i(t) < 0
x|[ri,))] ~ 8-46 VW (21,12) a(t) 0, i(t) > 0
x|[r1,12] ~ __ vt E (21,2) a(t) < 0, j(t) <
z||r,, +00) ~ s $ vt E +o0) i(t) > 0, ø(t) > 0
z|/r + 20) ~ if Vt E (r1 + 00) +(f) > 0, ø(t) < 0 and lim x(t) too
o|[rı, +00)
if Vt E (01, i(f) < 0, ä(t) and fime-++0x(((
x o|[r,, +00) -. if 14 E +00) i(() < 0, =(t)
z||r 00) 8 E +-0) +(f) > 0, =(t) and lim xxt)
x|[ri, +001 ~ s_ if E +o0) i(f) < 0, i(t) > and lim x(t) E RR
Definition Let uj i be bounded trajectory parts on (fo, tema) and assume v = (c ,P)). where (qj:2 pi) the bounded part of the semantic representation (c. p). We say that 2 seemingly matches (cj= P.), and write it i (c), if for any transition point in (cj.P)
15

---

Published as a conference paper at ICLR 2025
Table 4: Notation used throughout the paper (Part 1).
Symbol
Meaning
M
The number of dimensions in dynamical system, M E N
m
Used for indexing dimensions, m. E (),...,M,
to
Start of the trajectory. to R Time domain, subset of IR, (fo. +00)
An M-dimensional trajectory, A I-dimensional trajectory, x Derivative of æ with respect to time t The number samples, C Used for indexing samples, d .J....,D. The number of measurements of sample d, Na N
RM
æ
IR
t :
d
Na
The time the th measurement of sample d, t!! E T
The n measurement of sample (taken at time ((d))
), yo e R
(d)
yo
vo
initial condition, the value A system of M ODEs, MRM++ IR M A single ODE, f IR?
fo, ato E RRM
f
f
(
The of continuous functions on
(T)
The set of twice continuously differentiable functions on T forecasting model predicting an M -ddimensionnl trajectory, F; C2
F
A forecasting model predicting I-dimensional trajectory, F: F C
Fsm
semantic predictor, predicts semantic representation of trajectory from the initial condition
Fij
A Irajectory predictor, predict trajectory from its semantic representation
Ca
A composition of trajectory ır A set of properties of trajectory The last transition point, fend E T The set of all possible compositions The set of all possible properties for composition e The set of all possible properties for all compositions Motifs, formally defined in Definition 4 Semantic representation of x Semantic representation of F
Ps
food
C
Pa
P
8=
(Ca, Pr)
(Cp,Pr) Pr)
O
Function composition
h
Value of the horizontal asymptote in motifs 8-+A and 4-h
f1/2
"half-life" property of motifs s-+A and +-h
Fpop
A property map, F prop R P
Foom
composition map, From C
Fle)
A property sub-map. Fprop + Pe
(c,p)
P, (c,P) (cr. x E C seemingly matches (c,p) E P. Definition 5
C conforms to (c,p)
σ (c,P)
u (1) (1),
. à (t) (t) if E {tu, tend or t is local extremum,
u (t) 8 (t) ift tend is an inflection point.
We also say seemingly matches i and write it uj ~ U).
B ADDITIONAL RESULTS
B.1 EQUATIONS DISCOVERED BY SINDY
Below are the equations discovered by different variants of SINDy in Section 6.2
*(1) -3.06x(1)
(5)
16

---

Published as a conference paper at ICLR 2025
Table Notation used throughout the paper (Part 2).
Symbol
Meaning
set of trajectories
21
Bounded part of 21 [to, fendl + R
Unbounded part of 4, |2 [fend +o0)
(q.P))
Bounded part of (c, p). i.e. (c, p) without the unbounded motif and its properties
(1S, (P)
Unbounded part of (c, p). i.e., the unbounded motif, its properties, the last transition point, and derivatives at
Fo
C trajectory predictor, F C. P Co
C trajectory predictor, Fug C P C²
C
Composition library, C' C
I
The maximum number of intervals/branches the composition map Feom: I E N
a(t) = 5.56 43)
(6)
*(t) 27.56t + 14.83 co)(2((t)) 12.36 13.79 ex78xx)(x(t))
(7)
a(t) = 39828.60 50679.071 - 1071000000001777704444443-333333 - 54591.66666 cos(t) + 1.58 (2n()
+ 16141.41 .11 (2) + 1.63cos(3rrtt) - sinn)))s - 1366.28 c3s
(8)
+(t) = 8347995.48 + 3 3350455.09t
3.75x(t))
36.16 cos(x(t))
4904904.06 st) - 8903838.66 cos + 12.19 co)((((t))
(9)
+ 889658.07 sin(t)) + 560285.98 (20(()) - 0.29 sinn(tt)))
2.81 cos(3(tt)) - 02220.222992 (3t) 4394.42 cos(t)) - 5.23 exp(((t))
B.2 SEMANTIC REPRESENTATION AFTER EDITING
Figure 8 shows semantic representation of Semantic ODE after editing performed in Section 6.3
Composition 6. t t * -.
Property for composition (5 pohnts
Transtion pelnts p-coechatte
Propartia unbbuudedd
13h-/6
20
Figure 8: Property map of the Semantic ODE after imposing the horizontal asymptote to be 0.
B.3 SEMANTIC REPRESENTATION OF THE LOGISTIC GROWTH MODEL
Figure 9 demonstrates semantic representation of the logistic growth model si(1) xx1((1 -2
x(0)
B.4 INSIDE OF FITTING A COMPOSITION MAP
Figure 10 shows the logarithm of the loss for each sample and for each considered composition while fitting composition map.
17

---

Published as a conference paper at ICLR 2025
Compesition map
Property map for composition (s S..a) Transtio pain lecandraten
Transtion pa riL centiatal
barvatives
Propatiin beburdde
bolivde suraboe
650-
/80
composkion Trans pents teceadiriten
Trenitonn pot.mopptt-o-oreiiiit
barvathes
Poppttio urbeuntel
Meeual 44 #2.00
3.14
0.0
20-
holpural
Property map composkion (5- Trans pesto leczadirston
Thendton pe centreel
beraative
Popetion unbeanter
3.00
11
18:
Figure 9: Semantic representation of the logistic growth model.
003
01
042
828
032
31
0.2
841
0.6
05
084
81
85
77
085
09
Br..t..F..d ml
Ge.bk mg
bo an ad
Le d... od
in add
Figure 10: Logarithm of loss for each sample for each considered composition while fitting composition map to the pharmacokinetic dataset in Section 6.2
B.5 EXTENSION TO MULTIPLE DIMENSION: PROOF OF CONCEPT
In this section, we show a proof of concept of how semantic modeling can be extended to multiple dimensions (as was described in Appendix G.2)
18

---

Published as a conference paper at ICLR 2025
33
SPAOA PAOTIOn
2.9
193
04
1.4
1.2
|predicted] (pradicted) R (aredicted) (data)
-2.4
0:3
6.1
511
0.8
3322
80
1b
0.6
(data)
R (data)
6380
0.0
0.0
0.2
0.4
0.5
0.8
1.0
0.30
3.32
2.1
0.305
3.20
Chank
331
10.30
3.30
Figure 11: Semantic ODE model extended to multiple trajectories and fitted to data governed by the SIR epidemiological model. We show three generalized additive models describing the horizontal asymptote of R (top), the maximum of / (middle), and the time when / is at its maximum (bottom).
We implemented the property maps described in Appendix G.2 (each property described generalized additive model) and fitted data following an SIR epidemiological model for different initial conditions. To simplify the problem specify the composition map (we only train property maps). We assume S' follows (s--o: 4n), 1 follows (8++0) 84-or
8-4),, and / follows
(844e) -.) The predicted trajectories and some of the property maps can be seen in Figure
The average RMSE on the test dataset is 0.019 for S trajectory, 0.011 for I and 0.014 for R. Note that the irreducible error on this dataset (caused by the added Gaussian noise) is 0.01. The shown property maps let us draw the following insights about the model:
The time when at its maximum mmx (I) is on average just below 0.2. 1o has a relatively large impact on tmax (I) by increasing it by 0.1 for very low Io decreasing it by 0.05 for very high Tg- The larger the Iu the faster the maximum is achieved.
So also has a negative impact on mmax (7) but it is much smaller (+00.02).
The maximum of / (denoted 7mxx) increases linearly with both So and To. This time So has slightly bigger impact (++0.1) compared To (+:0.44)
In both Emax ([) and Imax- the impact of Ro is insignificant.
The horizontal asymptote of R increases lincarly with all three initial conditions. In particular, the shape function associated with Rg has unit slope expected.
Interesting advantage of our approach is that even though the system is described by three variables, do not need to observe all of them to fit the trajectory (similarly the pharmacokinetie example in the paper). ODE discovery methods assume that all variables are observed which constrains their applicability in many settings.
19

---

Published as a conference paper at ICLR 2025
B.6 DUFFING OSCILLATOR
Chaotic systems usually have some kind of oscillatory behavior (i.e., cannot be described by finite composition). As discussed in Appendix G.T it means that chaotic systems currently beyond the capabilities of Semantic ODEs as they would not be able to correctly predict beyond the seen time domain. However, we could use it for a prediction on bounded time domain. We have compared different models on the Duffing oscillator. The results can be seen in Table 6
Table 6: Comparison of Average RMSE obtained by different models on the Duffing oscillator
datasets (with two noise settings). Average performance over 5 random seeds and data splits is shown with standard deviation in the brackets.
Duffing oscillator
Method
low noise
high noise
SINDy-5
0.278 (.032)
0.389
WSINDY-5
0.262 .33))
0.361(,072)
PySR-20
0.312(.049)
0.396
SINDY
0.284 0266)
0.386
WSINDy
0.263(.027)
0.339 (033))
Neural ODE
0.212
0.291 (.021)
Neural Laplace
0.176 (.032)
0.299
DeepONet
0.429 .084)
0.528 (666)
Semantic ODE
0.096,.02)
0. 0.236(.001)
B.7 REAL DATASETS
We compare the performance Semantic ODE against other baselines on two real datasets. The tumor growth dataset is based on the dataset collected by Wilkerson et al. (2017) based on eight clinical trials. We follow the preprocessing steps by Qian et al. (2022). The drug concentration dataset is based on data collected by (Woillard et al. 2011). The results are shown in Table 7
Table 7: Comparison of Average RMSE obtained by different models on two real datasets. Average performance over random seeds and data splits shown with standard deviation in the brackets. "Semantic ODE*" is variant of Semantic ODE where incorporated semantic inductive bias about the shape of the trajectory. We specified the composition to always be (s_+c: 8444)
for the tumor growth dataset and (84-c:8 sh+) for the drug concentration dataset. Note, the implementation of WSINDy we used cannot work with trajectories as sparse as in the drug concentration dataset.
Tumor growth (real)
Drug concentration (real)
SINDY-5
0.243
0.286
WSINDY-5
0.237 0.237(.015) (.015)
NaN
PySR-20
0.536(.346)
0.257
SINDy
0.249(.029)
0.286,.01)
WSINDY
0.236
NaN
Neural ODE
0.228
0.263
Neural Laplace
0.243
0.302(.022)
DeepONet
0.242(.016)
0. 0.265 (200) 020)
Semantic ODE
0.234(.019)
0.264(.o21)
Semantic ODE*
0.229
0.243(003)
B.8 BIFURCATIONS
We believe that our framework is uniquely positioned to perform quite well on systems exhibiting bifurcations (when small change to the parameter value causes sudden qualitative change in
20

---

Published as a conference paper at ICLR 2025
the system's 's behavior). In our framework, bifurcation occurs when the composition map predicts different composition. As discussed in Ampppddix G.2 the future, the semantic predictor may take input not only the initial conditions but also other auxiliary parameters. We can then represent the composition map as a decision tree that divides the input space into different compositions. This decision tree then informs where bifurcations occur.
We hope the following proof concept based on the current implementation demonstrates that viable approach. Instead of predicting trajectory from its initial condition, fix the initial condition to be always the same and predict a trajectory based on the parameter that we observe in our dataset. We generate the trajectories given the following differential equation
(10)
The initial condition z00) 1 and sampled uniformly from (-1,2). We choose the set of compositions to be (S+-R), (5-+h) and record the position of the bifurcation point found by our algorithm (as opposed to the ground truth). The mean absolute error for different noise settings can be seen in Figure 12 Note that the range of values of the trajectory is (0, 2), and even in high noise settings, the location of the bifurcation point can be identified.
0.08
Mean absolute error
0.06
0,02
0.00
Noise
Figure 12: Mean absolute error on the predicted bifurcation point for different noise settings.
B.9 GENERALIZATION
show in Section 6.3 how our model can extrapolate well to unseen time domains. this section, we show how we make Semantic ODE generalize to previously unseen data points (initial conditions) by extrapolating the property map. We observe that each property function in Figure looks upproximately like linear function. Thus, fit linear function each of these functions and then evaluate our model on initial conditions from range (1.0, 1.5). Note that our training set only contained initial conditions from (0, 1). The semantic representation of the resulting model be seen in Figure 13] We compare the performance of this model to ODEs from Table 2 The results can be seen in Table 8 We can that our model has suffered only a small drop in performance even though it has never seen single sample from that distribution. It also performs much better than any other ODE tested.
C TRAINING OF THE SEMANTIC PREDICTOR
Training Semantic ODE requires fitting its semantic predictor (trajectory predictor fixed). This is done in steps. First, train the composition map Fuom and then the property map Fprop- It
21

---

Published as a conference paper at ICLR 2025
Com posptetion map
Property map for composition (4+ -2
05-
15
Figure 13: Semantic representation of the semantic ODE made generalize to zu E (0,
Table 8: Results of fitting Semantic ODE and SINDy to the pharmacokinetic dataset.
Model
Syntactic blases
Semantic hinses
Syntactic representation
Semaniic
(0, de (1,1.5)
representation
SINDS
oagi,
NA
a(t)
3.0
NA
0.222ao.1)
0.240
SINDY
Orde $
NA
56
8.112
SINDY
1 oegi 1 S
NA
Equarion (2)
NA
0.301 (0023)
7.76444.00)
SINDy
1id0. 10
NA
NA
0.106 (.029%
SINDY
2011
Oedi. 15 NA
Pquarion
0.020)
0.20330.40)
Semaniic ODE
NA
(84-3)
NA
Pigure 13
0.018;
0.023, oos)
is possible to adjust the composition map before the property map is fitted or even provide your own composition map without fitting This constitutes one of the ways prior knowledge can be incorporated into the model. Then the dataset is divided into separate subsets, each for different composition-acoon according to the composition map. Then a separate property sub-map is trained on each of the subsets. A simple block diagram of training semantic predictor is shown in Figure 6
C.1 COMPOSITION MAP
To fit a composition map we start with a composition library C' C C which set of compositions we want to consider. One can use a default set of compositions up to a certain length or filter out impossible ones to steer and accelerate the search. Then for every sample in our dataset, we measure how well each of the compositions fits the trajectory by fitting the properties that minimize prediction error. This gives example of such matrix can be seen in Figure 10 We then use a dynamic programming algorithm to neighboring intervals,
matrix where each row is a sample and each column is composition. An
find the best split of zo into up / intervals, with different compositions that minimize the overall prediction error for the whole dataset. We also make sure that each interval not shorter than prespecified threshold and contains a minimum number of samples. In our implementation, we choose the intervals to be at least 10% of the length of the entire domain and contain at least two samples. The number / is chosen by the user (we use 3 in all our experiments). This procedure is described in Algorithm
C.2 PROPERTY MAP
Property map consists of few property sub-maps each trained on different subset of data. Property sub-map predicts the properties, i.e., the coordinates of the transition points, necessary derivatives, and properties of the unbounded motif. Every single property is predicted by different univariate function. To get these univariate functions we first choose a set basis functions B. By default, we choose B-Spline basis functions, identity and constant. Then using just [B] parameters and efficiently evaluate it using a single matrix-vecte multiplication. This gives us, so-called, raw properties. In practice, different transformations to ensure that the predicted properties make sense. For instance, transition points are in correct relation another or the derívatives have the correct sign. This procedure is summarized in Algorithm 2
parameterize univariate function
need to pass some of these functions through
!-coordinates Instead of predicting the (-coordinates directly, we predict the intervals between them and the t-coordinate of the last transition point. We use softmax to transform raw properties into positive values that add up to 1. We then multiply it by the t-coordinate of the last transition point (that was obtained by passing a raw property through a sigmoid function scaled to cover the interval of interest). We can then use cumulative sum over the intervals get the desired l-coordinates.
22

---

Published as a conference paper at ICLR 2025
Algorithm 1 Algorithm for learning the Composition Map Feom-
Input: Dataset {(=)0) t(d),
(] compositions library C", maximum number of branches Z,
trajectory predictor Fimj C. x P - C
Output: Composition map Feom IR C
1: procedure COMPUTELOSS(d, c, p)
2: loss < Na -.E
Po (cg(c.))(th))) y/ (
3: return loss
4: end procedure
5: Initialize loss table LossTable(d, c), DP table L(d, 6) and backtracking table Bid, c, b)
6: for d 1-10 to D do
7: for all c E C' do LossTable(d, c) minpeP. COMPUTELOSS(d.cc2)
8:
9: end for
end for
(2) .... S 5x) (P) (1) S 20
11: Sort samples so that 20
for d D down
13:
for all E C' do
for i to 7 do
15:
if d = D then
16:
L(d, + LosTTblle(d c)
17:
else
18:
i = then
19:
L(d,c,0) LossTable(d,c)) L(d 1,c,0 c) t
20:
else
22:
Lary + L(s
23:
Lywach + LossTableed,, min L(d c, -1)
e'sc
24:
if Laua < Lgwitch then
25:
c,d,c,i) + Lauy c,d, c
26:
27:
else
28:
I(d, c, Liwitch B(d, i) + d' corresponding to mine' L(d+1,c', i,
29:
30:
end if'
31:
end if
end if
33:
end for
34:
end for
35: end for
36: (c",i") = arg min o c)
37: Initialize composition sequence {c )
of + i
39: for d = to D do
11=100 Eym1 B(d, c
41:
if f
# c then
42:
iti-1
43:
end if 44: end for
45: Define composition map F (z) using switching points 2/(() where ca + c241
r-coordinates Instead of predicting the r-coordinates directly, predict the absolute value of the difference between z-coordinates of the consecutive transition points. We do that by passing the raw property through softplus function. Then based the monotonicity of given motif, we either
23

---

Published as a conference paper at ICLR 2025
Algorithm 2 Algorithm for Learning the Property Map prop.
Input: Dataset {at (d)
: [(() y trajectory predictor Finj C - CO, composition map
From I C, set of basis functions B
Output: Property map Prrop R -
for d= 1 D do
Evaluate basis functions: (dd) (x) (1)
"..."""")........ where by E B
2:
3: end for
4: Form matrix B E Dx|8| with rows ((d)
5: for 1 D do
6: Compute composition: c(d) = Feom (=0 (=)²²
7: end for
8: Partition dataset into subsets {De} where D. {d cld) c}
9: for each subset Dr do
10: Initialize parameter matrix W(e) RI where P is the number of raw properties
11: Extract corresponding basis evaluations B repeat
Compute properties: P. (c)
14:
Use Pa to predict the properties { (Pp (1/ JLED i.e.
(1)) )
15:
Coordinates of the transition points Derivatives
16:
Properties of unbounded motif
DD EdED, Na E (Fay(e, ppop(() (t(hi)) 5/: )
18:
L+
19:
Update parameters ((() to minimize T. (e.g.. L--FGS) until Convergence of optimization
21:
end for
22:
Define Fpro as Fproo (20) = F'prop
(rFFmm(su) (xo)
add or subtract this number from the previous r-coordinate to obtain all z-coordinates respecting the composition.
Derivative at the first transition point We pass the raw property through sigmoid function and then we scale and translate obtain value specific range as described in Table 9
Derivatives at the last transition point If there is at least one bounded motif in the composition, then the last transition point needs to be either an inflection point or a local extremum. As such, either the first the second derivative vanishes and does not need be trained. The trajectory during training is predicted through Fin that can only accept one derivative constraint at the last transition point. That why the other derivative (the one that does not vanish) implicitly determined by this trajectory predictor instead of being trained explicitly. We first find the trajectory and then look the non-vanishing derivative and use it as the value predicted by the property map.
Properties the unbounded motif All raw properties for the unbounded motifs are passed through the softmax make them positive. We then adjust them case-by-can basis make sure that the motif properties make sense. For instance, instead of predicting the value of the horizontal asymptote in s directly, predict the distance between the asymptote and z (ttena). As the predicted distance is positive, after subtracting it from s(/end), we get a valid value for the horizontal asymptote.
After we get all the properties, we pass it to FD to predict the trajectory. We calculate mean squared enror between the predicted trajectories and the observed ones. We also add two additional penalty terms. One discourages too much difference between derivatives at the transition points and the other discourages too large derivatives at the end. We calculate the total loss and update the parameters. In implementation, use the L-BEGS algorithm (Liu & Nocedal 1989).
24

---

Published as a conference paper at ICLR 2025
D TRAJECTORY PREDICTOR
D.1 BOUNDED PART OF THE TRAJECTORY
We decide to define A as a set of cubic splines. They are piecewise functions where each piece is defined as cubic polynomial. The places where two cubics are joined are called knots. Cubic splines require that the first and second derivatives at the knots be the same for neighboring cubics so that the cubic spline is guaranteed to be twice continuously differentiable. Cubic splines promising because they are flexible, and for a fixed set of knots, the equations for their values and derivatives linear in their parameters. Thus by just solving set of linear equations, it is straightforward to find a function that seemingly matches (c)= PI) (denoted as #1 ~ (c,P) and formally defined in Definition 5 in Appendix A.2). i.e., passes through the transition points Pl and has the correct first derivative values at local extrema and boundary points (to. tend) as well as correct second derivatives inflection points and the endpoint (fena). The challenge arises because xj ~ (q): Pt) may not imply an (c),)). This is illustrated in Figure 14 Appendix D. However, it possible make the implication hold by imposing additional conditions. We prove this in Theorem
Theorem Let u, on be bounded trajectory parts on (ta, tema) and assume ij = (q. P.). where (c.p) semantic representation. both of the following hold for every pair consecutive transition points (a, b) of oy then in ~ P() 7 "1 = (q,P)).
Vt (a, b) sign(ū((t) = sign(ö (t)), and
if neither of {a, () a local extremum, sign( (1)) (ign( {a,b)
Proof. See Appendix D.I.I
We come up with two different ways of imposing these conditions that lead us to develop two trajectory predictors: Fin C P C and Figj : C - C², We use Finj during training of Fuem as is fast and differentiable, but the found trajectory may not be in C2 (only continuous). At inference, we use F2 that is slower and not differentiable but ensures that the trajectory is in C².
D.1.1 PROOF OF THEOREM 1
We provide proof of Theorem below. But before we prove it, we need the following lemma.
Lemma 1. Let ui be bounded trajectory parts on tend) and assume da = (q)PP)) where (c,p) p) is semantic representation, " (c), and the composition of uj is different from q. Then there exists a pair of consecutive transition points (a, b) of oj such that
31 E (a, b) such that sign(ii (1)) * (ügn(ū (1)) or,
if neither of {a, b) a local extremum, 3t E {a, b} such that sign(ū((t) * sign ( (f))
Proof. First, we will show that there needs to be two consecutive transition points of i/ a, b, where "I ~ s|[a, b) and oe s| [[, bj. Suppose for every pair of consecutive transition points of i) denoted a, b, v ~ al[a, b] and uj ~ =I[a, b] for some bounded motif 8. Then necessarily u has the same composition as bj (as they are defined the same interval). Therefore if the composition of u is different from then there needs to be two consecutive transition points of i) a, b, where 01 s|[a, and / ~ sl[a, 6).
We now consider two cases. Either there is a motif s' * such that 2 ~ s'I[a, b] or there no such single motif.
Case 1.
Let us assume that uj s'I b), where s' + is different bounded motif. Then, u (q. we know that and s' have the same monotonicity. If vj is increasing on [a, b] then vi > u (a). which implies u (b) > uy (a). Thus ij is also increasing on (a, bj. Similarly, if oj is decreasing on [a, b]. Let us assume, without loss of generality, that vj is increasing on [a, b]. If is also convex then i needs to necessarily be concave (as s' # s). That means that Vt E (a, b) sign (i (t)) # (t). Thus one of the conditions in Lemma iis satisfied. Similarly, if oj is decreasing concave.
25

---

Published as a conference paper at ICLR 2025
Case 2.
Let us assume that there is no single motif s' * s such that u ~ s'1[a, b]. That means that has a non-empty set of transition points Q (a, b).
If or b is a local extremum of 2y then there exists E Q such that q is an inflection point of oI because we cannot have two local extrema consecutive transition points. Thus i changes curvature q and necessarily there exists point t (a, such that sign(i (ú (t)) J sign(U (t)).
If neither a nor b is local extremum of v/, then either there is an inflection point in Q and we arrive the same conclusion as before or there are no inflection points in Q. In that case, Q contains only local extrema. However, we cannot have two local extrema as consecutive transition points. Therefore Q contains only one local extremum. That means that sign(1 (a)) o sign(ú (b)). But have ( (a)) sign( èj (6)) or (a)
(6) In all cases, either sign(ū (a(((t(()) # sign(ú (a)) sign( (b)) (ign(((()) Which is what needed to show.
We can now prove Theorem
Proof. Let assume that u (c.P)). To show " = (c).P) is sufficient show that u has composition e- For contradiction, let us assume that the composition of uj is different than q- By Lemma 1 there exist two consccutive transition points of v. denoted a and such that 3t E (a, 6) such that sign (ū (t) * sign(ū (t) or, if neither of {a, b} sign(v((o) or sign(u' (b)) & (bb(). This contradicts the assumptions of our theorem. Namely, that for every pair of consecutive transition points (a, b) of 2) Vt E (a, b) sign(ū (t)) sign(ū (t)). and if neither of {a, b) is a local extremum, sign(ū (1) = sign(ū for t E {a,b). b). Thus the composition of u is 9 and, therefore, u = ()).P)).
local extremum, sign(u' (x((a) *
D.1.2 CO TRAJECTORY PREDICTOR
Range of values for the derivative at to To ensure that the trajectory found by F has the correct semantic representation, we need to constrain the values of the first derivative of x at lo. These values depend on the slope the line connecting the first transition point with the second one, i.e., we define slope as K transition point. They are presented in the table below.
=((t)=z(ta) h- The values also depend on the motif and the nature of the second
Motif
4
Range
8+++
inflection
(0, k)
8+-b
maximum
(1.5x, 3x)
84-b
inflection
(K, 3x)
8-
minimum
(3k, 1.5k)
inflection
(3n, k)
S S--b
inflection
(xs, 0)
Table 9: Allowed values for i+(o) enforced by FD
Trajectory conforms to the semantic representation In this section, we slightly relax the defini- tion of conformity by allowing the part of the trajectory between two inflection points to be modeled as a straight line with an appropriate slope. Even though it does not match any of the defined motifs, in this section, assume straight line matches both the convex and concave variants of motifs with the corresponding monotonicity. That means we can relax the first assumption in Theorem to Vl (a, b) sign(ū, (t)) sign(ū ()))o or (() 0 if both a, are inflection points. To prove that ( ) (q. ).). we first show that Flg (q. (q.Pt).
Lemma 2. Let (c, E C be semantic representation predicted by From. Then Figg (9, P\) (q.P)).
26

---

Published as a conference paper at ICLR 2025
Proof. As each motif is described by separate cubic polynomial, it is defined by four parameters. To predict the trajectory, we set four constraints. Two for the z-coordinates of the transition points, and two for the derivatives (one for each transition point). If the transition point is to then first derivative at ln to the value specified in p.. it is a local extremum then we set the first derivative at it O. If it is inflection point then we set the second derivative at it to O. The only condition left to satisfy is the second derivative at feod. if fend is a local extremum, or the first derivative at fend if it is an inflection point. However, as we use mij for training the semantic predictor Fsem we guaranteed that the "other" automatically matches the one specified in p,- Therefore, by Definition 5 Figg (9.P) (q.P).
set the
We will know prove that P (9, = (c)
Theorem 2. Let (c,p) E CC* x P be a semantic representation predicted by Frem- Then Fiag F P.)
(q)P)).
Proof. By Lemma 2. we get that Pc) ~ (q),P)). We will now use Theorem to show that F (9.P)) = (I) Let i4 = (q.P)) and let us denote F (mmj(q-22 P)) u- Let us take any pair of consecutive transition points of ej and denote them as a, 6.
Let us consider two cases. First, when j lo and second, where lg-
Case 1. + #lo lo
At least one of (a, b) is inflection point, points. Denote this point as That means that ü (4) (g) = 0. By definition of & t is cubic on la, bj. That means üj is a straight line on [a, b). passing through 0 at Now we are going to consider two subcases. Fither the other transition point a local extremum is another inflection point.
local extrema cannot be two consecutive transition
Case 1./ The other transition point is a local extremum.
Without loss of generality, let assume that b = q is an inflection point, is local extremum, and that un (a) < uj (b). That means that some point E (a, b). in (1) > O. We also know that ú (a) That means that at some point t (a,b), ü > As új is straight line passing through 0 we get that Vt E (a, b) ü (t) O. This means that Vt (a, b) sign(ū (f)) sign(ü (t)). Similarly for other cases (where is an inflection point, or 24 (a) 24 (b)).
Case 1.2 Both transition points are inflection points.
If and b are both transition points, then Vt (a, b) (t) passing through (a, u(a)) and (b, u(b)). As such (i = (b) = sign(y 24 (a)) sign(v)() (b) vj (a)) = ( (a)) = (ign( (b)).
0. This means in is straight line
Case 2. to
If = to then (a) = ú(to) is in the range described in Table 9 Without loss of generality, let us assume that to 0. Let describe i on [a, as Bzt + 32t Bit + So. Let denote u (0) = wo and 2/ (t,) un Let also denote the required first derivative at 0 as ug- As uj needs pass through both (0, uo) and (t1, u1) and needs to have ú (0) u/g, we get the following equations.
S u0
(11)
BBG + B2l + 814 + 90 = U1
(12)
31 = u/
(13)
As in Table we denote m(t))- u(to) as K. Then U1 = :lo + K. l. The second equation then reduces
4,
to
Bat + B2t1 + 2/ =
(14)
Now we have to go by all six cases in Table 9
Case 2.1 84 & is inflection point.
27

---

Published as a conference paper at ICLR 2025
As t is inflection point, we get that
6/311 + 2/3² = =0
(15)
That means
82 = 333tt
(16)
By substituting into the previous equation, we get
Bat 3,3 + u/g
(17)
which gives us
B3 = (u6 - R))
(18)
2t
82 = 3(24 - K2)
(19)
26
To satisfy the conditions of Theorem holds if
we need to ensure that ūi > 0 for all E (0,4). This
3(1/ k)
3(4/ k)
0
(20)
17
4
This is equivalent to
3(16 - k)(t (1) >
(21)
As < h1, this is equivalent to v/ - which equivalent to
m K
(22)
As our motif is 5++b, v/g 2 O. Thus 0 S uig specified in Table 9
Moreover, as neither Ly nor ( is a local extremum, we need to check the signs of the first derivatives. We do not need to check to. but we need ensure that (t1) > 0. But this follows from the fact that in (4o) > 0 and ün (1) >> 0 for all I E 41)
Case 2.2 84-6. bb. & is a maximum.
As tj is maximum, we get
3,st + 282t1 + 2/ =
(23)
That means
82 = - 1/ - -33tt
(24)
261
By substituting into Equation (14) we obtain
3st + 1/2 x (-ub- 3B3ttt)
+ ung
(25)
By rearranging, we get
8z 11 2k
(26)
t
82 3R 2nh
(27)
t
To satisfy the conditions of Theorem we need
we need to ensure that ū (t) < 0 for all f (1). That is
6 "1
2x
3x
2u/
t ++2
2
< 0
(28)
f
t1
which is equivalent to
(6ub - 12x) + (6n - 416)11 <
(29)
After rearranging
<311 61) < u/g (211 31)
(30)
28

---

Published as a conference paper at ICLR 2025
Let us denote t/t, as t'. Then the previous equation is equivalent to f'(31/ 6Å) < 20l 3x
(31)
As this is supposed to be true for all t' E (0, 1). two things need to be true.
20h 3x > 0
(32)
3u6 6R < 20/0 3k
(33)
This gives us
1.5k < u/ < 3n
(34)
as specified in Table 9
Case 2.3 8 b. 1 is an inflection.
We follow a similar first step in Case 2.1. We arrive at the conclusion that to satisfy the condition of Theorem need ensure that (t) < for all E (0,t1). This holds if
u/² > K
(35)
However, we also need to make sure that ú (t1) > 0. That is, we need to satisfy
21
3 1/g K)
1-
+ u
(36)
21
211
That is equivalent
r/g 3x
0
(37)
Thus un needs to satisfy u/ < 3k specified in Table 9
Case -- 4 is a minimum.
The steps are very similar those in Case 2.2. We arrive at equation Equation (31) but the inequality has a different direction. That is, we need
t((3ug - 6x) > 2u1 - 3x
(38)
to hold for all e' E (0, 1). This can only be true if
2uh 3x < 0
(39)
3ul - 6Λ > 20/ 3k
(40)
which gives us as specified Table
un < 1.5k
(41)
Case 2.5 +6. l is an inflection point.
The steps are analogous to those in 2.3. However, we need to ensure that (t()) > for all t E (0,t2)
and that ú (4) 0. This gives
3R < u/ < K
(42)
as in Table 9
Case 2.6 h. t is inflection point. Analogously to Case 2, we get
18 S 1/ <0
(43)
D.1.3 2 TRAJECTORY PREDICTOR
In this section, we describe F2 that predicts trajectory that is twice continuously differentiable.
29

---

Published as a conference paper at ICLR 2025
Why is it challenging? A natural idea to find cubic spline with the corresponding composition and properties would be to follow a similar approach as in Fuaj by choosing knots at the transition points. However, then the problem turns out to overdetermined, i.e., we are not guaranteed that the needed cubic spline exists. Indeed, given S motifs in the bounded part of the composition, we have conditions for each of the internal knots (the values for the two cubic, matching first and second derivative, and the value of one of the derivatives). This is because each internal knot is either local extremum (first derivative vanishes) or an inflection point (second derivative vanishes). We also specify the values of the derivatives at the endpoints, which gives overall 5(S 1) + 2 conditions, whereas we only have 4S parameters.
5$
The conclusion is simple: constraints is 3( [ - 1) + 2(S - 1) + 3K + 25 where K is the number of cubics, and is the number of motifs. By equating this number to 4K, we get the optimal number of cubies is: 2S. That means that we will need knots between transition points. This, however, poses a new challenge. Although the resulting cubic will have all the transition points and derivatives as required, i.e., will seemingly match the semantic representation, we are in no way guaranteed that it will have the correct composition! Although we would like the function between two transition points to have a fixed sign of both derivatives, nothing is stopping our function from adding additional trend changes in between. This is visualized in Figure 14.
need more knots. In fact, the general formula for the number of
[t: -0
((b) 010))=000)-0
25
soints of Addi onal infection points of
(tt)) pits) (1ta) =g
15-
fita) 9'th)
0.5
(())) ((((a)
(15) -000ol
A(ty) = ('(a))
0.0
62
0.4
0.6
3.8
3.0
Figure 14: f (c, (q.P)) but g # (c).P))-
To solve this problem, we propose a completely different approach to fitting cubic splines.
From cubic spline to a piecewise linear function and back this subsection, we show how we can describe a cubic spline in terms of its second derivative to reduce the number of parameters and control its second derivative. A cubic spline consisting of K cubics is uniquely described by its K + 1 knots and 4K coefficients (four per every cubic). However, the actual number of degrees of freedom is much smaller because the spline needs to be twice continuously differentiable at the knots. We observe that instead and then integrating it twice to get the required function (Figure 15). The second derivative of a cubic a linear function (which can be integrated analytically), so each piece is described by just two parameters. In fact, we can describe a piecewise linear function by just the values at the knots, ensuring continuity. The values at two consecutive knots then uniquely determine the linear functions. Together with the additional two parameters for the integration constants (c and d), we reduced the number of parameters to K + 3. Formally, we want to find ü E Co and then define u E C2 as
decrease the number of parameters by fitting the second derivative of the spline
u(t)= = ( de" + dt + d,
(44)
where ü is described by (to...., tK, Co...., Uk). This not only helps us decrease the number of parameters and constraints we need to impose but, crucially, allows us to exactly control sign(ü(t))
by just making sure each iq is either positive negative.
30

---

Published as a conference paper at ICLR 2025
2: RR
R:
K + kintts
X Iinos
500) Fital
50h 5053
(4 for each oubio fi)
panemetace
X' parameters
Upoo My and
Mbll Eii)) Gitl)-00s)) gih) gina
guarentze
Insagration guarantess
3OX conotaintt
construints
-1-
Figure 15: Traditionally, a cubic spline over + knots is described by 4K parameters and 3(K - constraints that guarantee that the function is twice continuously differentiable. Instead, we describe it using its second derivative, which can be parametrized by K + parameters and two integration constants. Integration guarantees the function is twice continuously differentiable without any additional constraints.
Fitting the second derivative We observed that the second derivative of a cubic spline is piecewise linear function that entirely defined by its values at the knots. We choose the knots as the transition points and add one knot between every two of them. The positions of these knots are additional parameters of our model. Given n transition points, we have 2n knots (denoted to.... 22n-2) and their associated values denoted to. ..... U2n-2. As each antiderivative is determined up constant, we need additional parameters c and as integration constants. Ultimately, is described by (to....: f2n-2, 2, Ub, ... 2, c, d). Overall, this gives 4n parameters. However, many of these parameters are directly determined by the semantic representation of æ. In particular, the following parameters are predetermined.
lz; for all E [0 n - 1] (I-coordinates of the transition points)
U2n-2 and all 122 where 12 is inflection point
c =tot d =fff
In addition, for every transition point l2 with specified first derivative (either a local extremum or the endpoint), the value U2--1 chosen to enforce the correct value of ú(t::). This can be done by observing that the value of i((2s) is equal to the sum of ú(t2z), where tzj is a previous transition point with a specified first derivative, and the integral of ü between those points. Namely,
12i ü(t) dt + ú(t2s)
il2)
(45)
This integral can be calculated analytically as ä is piecewise linear function. In the end, we are left with 2n 2 parameters obeying the following constraints.
12 < 121+1 < 121+2 for all E [0: n --2 -2]
sign(w:) sign(E(t.))
The goal is to find these parameters such that the resulting u(tzs) = 2(t2)) and sign(ú(tzi) (tt2)) and this can be formed as the following objective.
-((())))2(ttt) - =(tx))) + Amin(((te)) x (tt:),.)
(46)
where À is chosen to be very large. We then use L-BFGS-B and Powell to optimize this objective three times for different initial guesses. the maximum error on the transition points (max, (max;/u(lz) ((((2))=(((()))) ) (((2::)) is lower than some user-defined threshold (we
and the derivatives (max;|ú((2:) 0.001) then we have found a function with the correct semantic representation. If not then we default using Fo The threshold allows the user to choose trade-off between how smooth the trajectory needs to be and how faithful it needs to be to its semantic representation.
To make this concrete, we consider an example where the finite part of the compo- sition is (s+ (s+ b, tb, 8++6) (see Figure 16). It is described by 20 parameters
31

---

Published as a conference paper at ICLR 2025
C, (Co...., Eg. 0o. .o..... Ugg C, d), but most of them fixed (denoted by the red color in Figure 16). In par- ticular, to. t2, tq, to. by are the t-coordinates of the transition points. D2 = 0 as an inflection point and Ug as the last transition point. Moreover, 21 U5, 27 always chosen to ensure the first derivatives at 12, ls, is have the correct values and c (((n) = 1, d = s=lo) = 0. That leaves trainable param-
eters (Vo, 11, 12, 13, U3, ls, U6 17) constrained such that < b < 02 13 14 is ls l ls. vo, D2, 03 < O, VG > 0.
1.0
Values of [trainable]
Transition paints (fiedd)
Additional knots (trainable)
0.5
X
A2
0.5
-1.0
Figure 16: Example parametrization of u through its second derivative. has the composition (84-b, 8_ 8 4by 84+)). It described by 20 parameters (lo,..., Lg, 0o..... Ug C, d) where of them are fixed and 8 are trainable.
D.2 UNBOUNDED PART OF THE TRAJECTORY
In contrast to the bounded part of the trajectory, for the unbounded part have only predictor. It takes the unbounded motif and its properties, as well as the coordinates of the last transition point and both derivatives this point, and predicts trajectory jz. We need the derivatives guarantee that the trajectory is twice continuously differentiable at fend- The first challenge is to come up with useful properties that sufficiently describe the unbounded motif.
D.2.1 PROPERTIES OF UNBOUNDED MOTIFS
Below, we describe the properties of unbounded motifs one by one.
84+4 This motif represents an increasing, convex function. We decided to model it as a function exhibiting an exponential-like behavior. That is, we expect to double after fixed time interval. As it is true for a function such as u(t) = 3 x 23. to generalize the notion of this doubling time to settings with arbitrary initial conditions, including negative ones. In that case, the doubling time depends on t. However, we can make it approach a fixed value asymptotically. Let (t) be defined value such that a=(t (t)) 2,x(t). As j is increasing and convex, (t) is well-defined. Then we define our property as lim (t) and call it "asymptotic doubling time".
this case +) = x u(t). However, we need
8--u We describe s--. using an analogous property that call "negative asymptotic doubling time" defined in exactly the same way. The only difference is that, in this case, doubling makes the trajectory more negative and thus the value decreases.
84-0 This motif represents increasing, concave function. We decided to model it as a function that resembles a logarithm. Thus, we expect it to have an inverse behavior to the exponential function. As double the input, expect constant increase in the function's 's value. As with the previous examples, in general, this increase would depend on the current value of the trajectory. but we can make it converge / xc. Formally, define 7() number such that (x(21) (z(()) +1)) As
32

---

Published as a conference paper at ICLR 2025
a increases, this is well defined. Thus, we define the property lime-+20 ((t and call it "asymptotic incrementing factor".
--+u We describe s-++ using analogous property that we call "asymptotic decrementing factor". is defined as lim r(t), where y(t) is defined as number such that (a((t) = jr((t) - (t))
This motif represents a decreasing, convex function with horizontal asymptote. A natural property is the horizontal asymptote h. However, we decided to include one more property that describes how quickly the trajectory approaches that asymptote. We define 41/2 as the value of I where j is in the middle between the value of the last transition point and the asymptote. Formally, (2(1122) ++ We call "half-life" as would correspond to half-life if modeled as an exponential decay.
8-th
2
84-h We describe this motif using the same properties as s-+h- However, as the function increases, instead of "half-life" we call it "inverse half-life".
D.2.2 PARAMETRIZATION
Having defined the properties of interest in the previous section, we need to find a parametrization for each of these motifs that would allow to choose arbitrary property well as the last transition point and the two derivatives at this point.
We parametrize this motif 8+++
ja(t)= (2( + 0z(1 tend) + 01 lend) + Os
(47)
We can quickly calculate that
pplend = 01 + 0,5
(48)
(#(le = 0,02 + 04
(49)
[ = 0,65 + 20,
(50)
Now let us define (t) as earlier, i.e. (t) is the value such that jz(t + x+)) = (t). We are now
going to prove that 1om g()
log(2) which we denote as simply y. To do we first prove that 0g
lim; (t) - (a(t + 7.) = O.
lim (t)(()
2 [0 ((())(())))) + 0z(1 lema) + 0a(1 lend) - 0s]]
lim
1+000
po(t + 74)
O2 (t tww) lz(t + ? tend) + Oz(t 7x- -tend) + 0g
20 00t-
lim
2
(51)
627.
2
(2( loc2
From that it follows lime- 2px(() 1 7-) 0. As the function
2x((0)
is continuous in both
t and 7 and for every t >> Lend it is bounded for all y > O, it is absolutely continuous in From that it follows that limy-xxo g(() = g. 02
(2)
So we can set the parameters as follows.
log(2)
(52)
Yo
00 (pā(tend)
0,02)/2
(53)
04 (tend)
00000 0z
(54)
0, (tonnd) 01
(55)
33

---

Published as a conference paper at ICLR 2025
We have some freedom in choosing 81, but need to make sure that the function is indeed increasing and convex for all f fend- This is true 0g >> and 01 > 0. So we choose 81 to be
0 = min( (fena) ) /20,, (femm)) /22 1)
(56)
S--b
We parametrize this motif analogously to the previous one.
(=(t)= 0. Oy(t tend) 0s tend) + 0g
(57)
where
0 (2)
(58)
03 ((((ban)) 0203)/2
(59)
0 [ 0102
(60)
0g- (f(ood + 01
(61)
and 0, needs to be chosen to make sure that 0g and 04 positive. We choose it to be
01 min(-mi-a (atimm(tttena)/22000 i(tond)/205)
(62)
8+-u
We parametrize this motif
(x((() 0, log((2 lend)? + 0z(1 lend) + 1) + 04
(63)
where
01 7a/log(4)
(64)
03 = (#end) 101 = log (tena)))))))
(65)
02 02/2 2 (log( (4) (eae//mn//)/)))
(66)
04 = (=(tnn))
(67)
where is the property of 84-u that we call asymptotic incrementing factor.
is easy to verify that this function is indeed increasing and concave. there is bounded motif before it, it has to be $++e, and then fend is an inflection point. Thus
(tend) = 01 (20, 0% =0
(68)
as required. If there is no bounded motif before it, there is no way to set the second derivative at lend to any other value than O. This property fixed and not trained in the property map.
Let (t) be defined as earlier, i.e.,
7(1) = (=(2)) (z(()
(69)
We can compute it
(()) = .<2)) xx(()
01 log (f2(12221 lend) 2 + (2 lend) 1) O1 log(O2(( lend)² + Oz lend) + +1) (70)
(22 tend) + 02 tend) ) +
01 log
(1 lend)2 + 0z(l lend) +
Thus
[f02
lim
7(t) 81 log
2-000
02
(71)
01 log(4)
7
required.
34

---

Published as a conference paper at ICLR 2025
8-+u
This motif is parametrized analogously 4
(+()) = - log(O2(t - lema)² + 0z lend) + 1) + 04
(72)
where
0 7=/log(4)
(73)
83 -1
| (tend)) 001 - -log( loggll))))))))) (f)oa))))))))
(74)
82 = 0 - (log
12/2
(75)
04 = (a(aa
(76)
where 7o is the property of 8 that call asymptotic decrementing factor.
Let us define ry as earlier, i.e.,
((t) = po((t) (z((2))
(77)
Then
1
lim 7t) -0 log t+000
(78)
01 log(4)
7+
as required.
s 8-th
We parametrize this motif as
(x(t) h(g(t tend)
(79)
where
01
h(t) = 1 + et and g(t) is appropriately defined cubic spline.
+ 02
(80)
First, we need to make sure that p is actually decreasing and convex.
(a(()= = hhggt- lend) (ena)))(( - lend)
(81)
016¹
h(t1) (e' 1)2
(82)
So to satisfy (=(() < 0 for / 2 lend, need (1) > for 0. Let now look at the second
derivative.
a(t)= h(g(t tend)) tend) h(g(t tend))ā(t tend)
(83)
((t = (e' - l)et
(84)
(e' + 113
For > 0, this is always positive. Let us assume g(0) = O, then (t) > 0 for t > 0. To satisfy
(tt) > 0 for t tends need ğ(t) < 0 for t > 0. Let now look fend-
(=(eeed)) h((()) ko(0) 02 02
(85)
(fcnd)) h())((0) 0 -0)
(86)
[i( = (()))(()² + (()))(() = 0 6(0)
0
(87)
In addition, we also need to satisfy the properties, i.c.,
lim |z(()= = h 8-000
(88)
( = (Io(Fen)) + h2)/2
(89)
As is increasing,
lim (z(() O2 (-000
(90)
35

---

Published as a conference paper at ICLR 2025
Thus 02 h. From that follows that 81 ((x((fend) - h). From the "half-life" property, get
11) 1 + 6 fow) Fou)
+ h (p2((end) h2/2
(91)
From that, we get that g(ti/2 )nnd log(3). To summarize, we need to find such that
g(0) = 0
(92)
it (((end)
à(0) =
(93)
z x(tens) h
ī(0) = 0
(94)
g(1//2 fend) log(3)
(95)
We impose ī(0) = O, so that (f(andd) = which is always the case if there is bounded motif before. If it is the only motif, then one of the property maps is fixed and not trained.
To make sure that j(()) 2 0 for and straight line. Similar to approach in inAppeedixx D.I.3] describe the two cubics using second derivatives. Thus ğ is a piecewise linear function with knots at 0, f 22 such that
2 we define as piecewise function composed three cubic
ā(0) =
(96)
ğ(t1) 01
(97)
x((2)
(98)
(99)
and j(t) 0 for all t > 12. The goal is to now find t,, such that g is increasing, concave, and g(41/2 lend) log(3).
is always going to be concave if < 0. However, to make sure it stays increasing, we need to have (tz) > This requires g(0) + (t2 x v1)/2 > 0. tf - tead < 1.5 (o5loo3))/il then we set to 1122 - fend otherwise we set 12 1.5 log(3)/}(0) In both cases, we choose & = f2/ We used symbolic Python library sympy to arrive the following conclusions. In the first case, we set V1 as follows.
4 (11/2 lena) + log(3)
(100)
(t1/2 )
In the second case, we set it as follows.
(11/2 lend (0) + log(3)
1i Uj = 2 (11/2 lema) (1.5 (og(3) - (1.5 log (3))
(101)
From that, calculate the coefficients of the cubics and the slope and intercept of the straight line.
84-h This one analogous to 8-4h-
E EXPERIMENTAL DETAILS
All experimental code can found at (btteprs/ggtuuu..ccc///aarrrrrrp Semanticone
E. DATASETS
Pharmacokinetic model The pharmacokinetic dataset is based on the pharmacokinetic model
developed by WWillard et al. (201T) to model the plasma concentration of Tacrolimus. This model
36

---

Published as a conference paper at ICLR 2025
consists of system of ODEs described below.
dCcepot
kir Caepot
(102)
dt
dCmans
ku Caepot
hig Gransl
(103)
dt dCirmme dt dCuanss dt dCCent dt
knr Ceransl
kir Cirans2
(104)
kig crans2
Rig tramss
(105)
kn Crans3
((CL+ + Q) Ceeon + (Q C22
(106)
dCperi (Q Ceent/ /V1) (Q Cpeni V2)
(107)
dt
The values of the parameters and the initial conditions are presented in Table 10
Table 10: Parameters of the pharmacokinetic dataset.
Parameter
Value
CL
80.247
V1
486.0
79
V2
271
li
3.34
Caspor (ln)
10
Ceemt (4))
2o (Vi/1000)
peri (to)
V2/Vi Coea (to)
Cranst (to)
0.0
Cuaas (t)
0.0
Crans3 (f))
0.0
We create 100 samples of xo equally spaced between 0 and 20. We then solve the initial value problem to obtain a trajectory of Ceemt We then scale it back to appropriate units by multiplying by 1000 and dividing by V1. We observe each trajectory 20 equally spaced time points between 0 and 24. Then we scale the dataset by dividing the concentrations by 20 and the time points by 24. Finally, we add Gaussian noise with a standard deviation 0.01 for the low noise setting and o 0.2 for the high noise setting.
The out-domain dataset created similarly, but each trajectory is observed only to 0 and at 20 time points between 24 and 48 (1 and 2 after dividing by 24).
Logistic growth The logistic growth dataset is described by the following equation (Verhulst 1845).
a(t) = (1)(==(() =(t))
(108)
We create 200 samples xg equally spaced between 0.2 and We then solve the initial value problem to obtain trajectory for each. We observe each trajectory 20 equally spaced time points between 0 and Finally, we add Gaussian noise with a standard deviation o 0.01 for the low noise setting and 0.2 for the high noise setting.
General ODE
The general ODE dataset is descrihed by the equation
+(t) = (x)
(109)
where ff is not described by compact closed form expression. We describe f using probability density function of a 2-dimensional mixture of 3 Gaussians. The means, covariances, and weights for each Gaussian shown in Table 11
37

---

Published as a conference paper at ICLR 2025
Table 11: Parameters of the Gaussians used to define the general ODE dataset.
Weight
Mean
Covariance matrix
0.4
0.
[], 0]. [0. 1]]
-0.3
3. 3]
0.0 [0.0, 1J]
0.3
[-1,- 2]
[2, 02, (0, 2]]
We then define f(x, b) the value of the probability density function of this mixture multiplied by 50. We create 200 samples of Zg equally spaced between -3 and 3. We then solve the initial value problem obtain trajectory for each. We observe each trajectory at 20 equally spaced time points between 0 and 5. Then we scale the dataset by dividing x by 3 and the time points by 5. Finally, we add Gaussian noise with standard deviation 0.01 for the low noise setting and 0.2 for the high noise setting.
The Mackey-Glass dataset is described by the following Mackey-Glass equation (Mackey & Glass 1977
Mackey-Glass
a(1) = 000"x(( 7)
7xr(1
(110)
8r + x(t - 7)
We choose the following parameters: 0= 1, B - - 4.0, n 4. We create 200 samples of mo equally spaced between 1.0 and 3.0. We then solve the initial value problem to obtain trajectory for each. We observe each trajectory at 20 equally spaced time points between 0 and 30. Then we scale the dataset by dividing x by 3 and the time points by 30. Finally, we add a Gaussian noise with standard deviation = 0.01 for the low noise setting and a = 0.2 for the high noise setting.
Integro-DE
The integro-differential equation dataset is described by the following equation (Holt al. 2022 Bourne 2018).
#(1) = --2(() - x())ds
(111)
We create 100 samples of zo equally spaced between 1.0 and 1.0. We then solve the initial value problem to obtain trajectory for each. We observe each trajectory 20 equally spaced time points between 0 and 5. Finally, we add Gaussian noise with a standard deviation o = 0.01 for the low noise setting and 0.2 for the high noise setting.
F.2 METHODS
SINDY We use SINDy (Brunton al. 2016b) implemented in the PySINDy package (de Silva ct al,, 2020; Kaptanoglu et al. 2022). We pass the variable t as an additional dimension of the trajectory allow for time-dependent solution (not autonomous systems). We use the following library of functions:
1, xr, 1, xt, P, e", et, sin()), sin()))inn(()) cos(rr, cos(1), sin (2x),, (21(, c2s().r), (2)) sin( (31)
We use Mixed-Integer Optimized Sparse Regression (MIOSR) (Bertsimas & Gurnee 2023 for optimization as it allows us to choose a sparsity level-the number of terms in the equation. In our experiments, consider two variants of SINDy, that we denote SINDy and SINDy-5. SINDy-5 enforces the maximum number of terms to be 5 which would hopefully allow for analyzing the equation. In the SINDy variant, tuning (between 1 and 20). For both SINDy and SINDy-5, we tune the parameter a of MIOSR that describes the strength of L2 penalty (between 1e 3 and 1), and the derivative estimation algorithm. We choose between the following techniques: finite difference, spline, trend filtered, and smoothed finite difference as available PySINDy. The parameter ranges we consider for each of them are shown in Table 12)
choose the maximum number of terms during hyperparameter
We use WSINDY (Reinbold et al. 2020 Messenger Bortz 202 based on (Reinbold WSINDY
al. 2020) implemented in PySINDy package (de Silva al. 2020 Kaptanoglu al. 2022) We
38

---

Published as a conference paper at ICLR 2025
Table 12: Hyperparameter ranges for each of the derivative estimation methods.
Method
Hyperparameter ranges
finite difference
E [1....,5}
spline
E (le 3, order E 20, a E (le - 4,1) window Lealow fl,...,5)
trend filtered
smoothed finite difference
pass the variable I additional dimension of the trajectory allow for time-dependent solution (not just autonomous systems). We use the following library of functions:
1, 1,2² 2 xt, l e" e", sin(x), (in()))))sn)))))))))) sin((), cos(1)), cin()))oss(x),, cos
sin (2x), (2t), (2x) (2t), (in(3x), sin(
We use Mixed-Integel Optimized Sparse Regression (MIOSR) (Bertsimas & Gurnee 2023 for optimization as it allows us to choose a sparsity level-the number terms the equation. In experiments, consider two variants of WSINDy, that denote WSINDy and WSINDy- 5. WSINDy-5 enforces the maximum number of terms to be 5 which would hopefully allow for analyzing the equation. In the WSINDy variant, we choose the maximum number of terms during hyperparameter tuning (between and 20). For both WSINDy and WSINDy-5, we tune the parameter of MIOSR that describes the strength of the L2 penalty (between 1e and 1). We choose the parameter K (the number of domain centers) to be 200.
well-known symbolic regression method to ODE dis- covery by first estimating the derivative and then treating it as a label. We choose the derivative estimation technique and its parameters by hyperparameter tuning (as with SINDy) using the methods and parameter ranges in Table 12 We use the following operators and functions.
PySR We adapt PySR (Cranmer 2020
+, sin, exp, log(1 + [x]))
Following the advice the documentation, we also put a constraint to prevent nesting of sin functions. During hyperparameter tuming we allow PySR to search for 15 seconds and then we train the final model for minute. This bit larger but comparable time budget to Sematic ODE and much bigger time budget than SINDY and WSINDy require.
NeuralODE We implement NeuralODE (Chen et al. 2018 model using torchdi ffeq library. We parametrize the ODE as a fully connected neural network. The data is standardized before fitting. We set the batch size 32 and train for 200 epochs using Adam optimizer (Kingma & Ba 2017). We tune hyperparameters using Opt una (Akiba et al. 2019) for 20 trials. Ranges for the hyperparameters are shown in Table 13
Table 13: Hyperparameter ranges used for tuning Neural ODE,
Hyperparameter
Range
learning rate number of layers
(1e-5,le-1)
(1,3)
units in each layer (separately)
(16,128)
dropout rate weight decay activation function
(0.0,0.5)
(1e-6,1e-2
ELU, Sigmoid
DeepONet We implement DeepONet (Lu et al. 2020 using a fully connected neural network. The data is standardized before fitting. We train for 200 epochs using Adam optimizer (Kingma & Ba 2017). We tune hyperparameters using Optuna (Akiba et al. 2019) for 20 trials. Ranges for the hyperparameters shown in Table 14
39

---

Published as a conference paper at ICLR 2025
Table 14: Hyperparameter ranges used for tuning DeepONet.
Hyperparameter
Range
learning rate number of layers
(Ie-5,le-1)
(1,5)
number of hidden states
(10,100)
dropout rate weight decay batch size
(0.0,0.5)
(1e-6,1e-2)
{8, 16, 32)
DeepONet We implement Neural Laplace (Holt et al. (002)) using a fully connected neural network. The data is standardized before fitting. We train for 200 epochs using Adam optimizer (Kingma & Ba 2017). We tune hyperparameters using Optuna (Akiba al. 2019) for trials. Ranges for hyperparameters shown in Table 15
Table 15: Hyperparameter ranges used for tuning Neural Laplace.
Hyperparameter
Range
learning rate number of layers number of hidden states
(lc-5,le-1)
(2,5)
(10,100)
latent dimension dropout rate
(2,10)
(0.0,0.5)
weight decay
(1e-6,1e-2)
batch size
{8, 16, 32)
Semantic ODE
In all experiments in Section 6.4 we use full composition library containing all compositions up to 3 motifs (with the exception of logistic growth dataset where we only consider compositions up to 2 motifs). We choose the maximum number of branches for the composition map to be 3. Each univariate function the property maps is described as a linear combination of basis functions: constant, linear, and four B-Spline basis functions of degree 3. The property maps are trained using L-BFGS as implemented between derivatives to be 0.01 and we perform hyperparameter tuning of each property sub-mat map to find the optimal learning rate (between le-4 and 1.0) and the penalty term for the first derivative at the last transition point (between le-9 and le-1).
PyTorch. We fix the penalty term for the difference
E.3 BENCHIMRRKING PROCEDURE
The experimental results presented in Table 3 are result of the following procedure. For 5 different seeds, the dataset is randomly split into training, validation, and test datasets with ratios 0.7 0.15 0.15. For each method, performance on the test set. Each seed results in both a different split and a different random state of the algorithm (apart from SINDy, which is deterministic). As Semantic ODE splits the training dataset and trains separate property sub-maps, we pass the combined train and validation set to it. For each property sub-map, the subset of this set with the corresponding composition once again split into training and validation subsets. The validation set is used for hyperparameter tuning and for early stopping.
perform hyperparameter tuning for 20 trials and then report the
F EXTENDED RELATED WORKS
Symbolic regression and discovery of differential equations The discovery of differential equa- tions usually considered a part of broader area called symbolic regression. Symbolic regression is the area of machine learning whose task is to describe data using closed-form expression. Traditional symbolic regression has used genetic programming (Stephens) 2022] Cranmer (2020) for this take but recently neural network has also been utilized for that task. That includes representing the equation
40

---

Published as a conference paper at ICLR 2025
directly neural network by adapting the activation functions (Martius & Lampert 2017: Sahoo
iet al. 2018), using neural networks to prune the search space (Udrescu & Tegmark] [020 [ drrescu
et al. 2021. searching for equations using reinforcement learning (Petersen et al. 2021). or using
large pre-trained transformers (Biggio et al. 2021 D' Ascoli et al., 2022). Designing complexity met- rics and constraints that make sure the equations are simple enough analyze is itself a challenging research problem (Kacprzyk& & van der Schaar 2025). Standard symbolic regression can be adapted to ODE discovery by just estimating the derivative from data and treating it as a target (Quade et al. 2016). However, many dedicated ODE discovery techniques have been proposed. Among them, the most popular is SINDy (Brunton et al. 2b that describes the derivative as a linear combination of functions from a prespecified library. This was followed by numerous extensions, including implicit
equations (Kaheman et al., 2020), equations with control (Brunton et al.) 2016b), and longitudinal treatment effect estimation (Kacprzyk et al., (002aa). Approaches based on weak formulation of ODEs that allow to circumvent derivative estimation have also been proposed (Messenger & Bortz 202 Qian et al.] 2022). ODE discovery methods usually cannot be directly used to discover partial differential equations (PDEs) and thus many dedicated PDE discovery algorithms have been developed (Rudy et al., 2017; RRissi & Karniadakis] 2018 Messenger & Bortz, [2021b Kacprzyk let al. (0023). The challenge of finding compact and well-fitting closed-form expressions inspired Shape Arithmetic Expressions (Kacprzyk & van der Schaar 2024) that extend the prespecified set of well-known functions (e.g., exp or trigonometric functions) in symbolic regression by flexible and learnable univariate functions that do not have a compact symbolic representation but can be comprehended by looking at their graph.
Time series representation Semantic representation is closely related to the topic of time series representation. In particular, apart from the work by Kacprzyk et al. (2024b), there have been other works that try to symbolically describe the time series. In particular, Symbolic Aggregate approXimation (Lonardi & Patell 20022), Shape Description Alphabet JAndré-Jönsson & Badal 1997). or the triangular representation of process trends (Cheung & Stephanopoulos 1990) represent a time series sequence of symbols that resembles the definition of a composition. However, these methods are mostly used for data mining or classification rather than forecasting. Similarly, although shapelet- based methods (Ye & Keogh 2009) and motif discovery (Torkamani & Lohweg 2017) are concerned with the trajectory's 's shape, they are usually aimed at finding subsequences of a time series that represent the most important or repeating patterns. There are also numerous other nonsymbolic time series representation techniques ranging from Piecewise Aggregate Approximation (Yi & Faloutsos 2000) to signatures [Lyons. 2014).
Shape preserving splines As we describe the bounded part of our trajectory cubic spline, and we want it to have a specific shape, this may seem related to an area of machine learning called shape- preserving splines. Its main goal is to interpolate data while maintaining essential characteristics such as monotonicity, convexity, or non-negativity (Fritsch & Carlson 1980. Pruess 1993). Although their goal is to interpolate rather than reconstruct, and they rarely consider the exact positions of the local extrema and inflection points, ideas from this field may prove useful in designing better and more efficient trajectory predictors in the future.
Neural ODEs Neural ODEs (Chen et al. 2018) provide flexible way to model continuous
dynamical systems by parameterizing the derivative function of an ODE with a neural network. While this allows Neural ODEs to capture complex system behaviors, they often operate as black-box models, making the task of analyzing to obtain the semantic representation even more difficult than in the case of closed-form ODEs.
G ADDITIONAL DISCUSSION
G.1 LIMITATIONS
Finite compositions Throughout our work we implicitly assume that the compositions are finite. However, that is not always the case. In particular, any trajectory that has an oscillatory/periodic behavior has an infinite composition. For instance, sin. Although Semantic ODE can fit such trajectory on any bounded interval, it cannot yet predict it into the future. This could possibly be addressed by extending the definition of semantic representation by an additional layer that
41

---

Published as a conference paper at ICLR 2025
describes periodic trajectories. The idea is to take the infinite composition (infinite sequence of motifs) and represent it using some finite representation. For instance, a trajectory sin(t) would be described as meta-motif 1St-es8--od te: s+++e) repeating forever, where the meta-properties may include the "frequency" and "amplitude" that may also vary with time and be itself described using compositions. This extension still assumes there is some underlying pattern behind the infinite sequence of motif's that can be "compressed" into a shorter representation. If the infinite sequence of motif's algorithmically random, then motif-based representation may not be useful.
Long compositions
The current implementation may struggle if the ground-truth composition is finite but long. In that case, it may fall outside of the chosen set of compositions.
Chaotic systems Chaotic systems usually have some kind of oscillatory behavior (i.e., it cannot be described by finite composition). As discussed above, means that chaotic systems are currently beyond the capabilities of Semantic ODEs. They would not be able to correctly predict beyond the seen time domain. However, we could use for a prediction on a bounded time domain. We have included such experiment in Apppndii B.6 Although this performance would drop drastically if we tried to predict beyond the training time domain. In general, Semantic ODE can model big changes in the trajectory following small change in the initial condition (as is characteristic for chaotic systems). As the composition map is discontinuous at points where the composition changes, it can model a sudden change in the shape of the trajectory.
method achieves decent results,
Training of the composition map The composition map determines how the dataset is split and which property maps are trained. Thus, it is very important to do it well. Currently, we evaluate how well a particular composition fits a trajectory individually for every sample. However, then we fit property sub-maps where these properties vary smoothly between samples. Thu,s the neighboring samples should be taken into consideration when evaluating how well a particular composition fits a trajectory. If the same composition fits two neighboring samples well but for very different property values, this may not be a good match.
Ubounded motifs The type of unbounded motifs we choose and how we decide to parameterize them influences how well we can fit trajectories. For instance, deseribe 84++u behaving like an exponential function in the long term. This may not work well if the actual trajectory behaves like a quadratic function, for instance. To address it, we should have more motifs available, so that we can be more certain that we can find a well-fitting motif. However, this also makes training the composition map more difficult. This emphasizes why finding an effective and efficient training procedure for the composition map is essential.
One-dimensional systems The current implementation of Semantic ODE works only for one- dimensional systems. We describe a possible roadmap for realizing direct semantic modeling for multi-dimensional systems in Appendix G.2
G.2 ROADMAP FOR DIRECT SEMANTIC MODELING IN MULTIPLE DIMENSIONS
Although Semantic ODE realizes direct semantic modeling for one-dimensional trajectories, we believe direct semantic modeling can be successfully applied to trajectories with multiple dimensions. To realize that, we first need allow for multi-dimensional inputs to the semantic predictor. As the composition map is just classification algorithm and the property maps are just sets of static regression models, this should be possible. We could imagine composition map represented as a decision tree and the properties being predicted by, for instance, generalized additive models (GAMs) (Hastie & Tibshirani 1986). that way, we can still have an understandable model representation. The major challenge is to come up with efficient optimization procedure that would not only search this complex space but also make sure that every predicted semantic representation is indeed valid. An additional challenge is to keep the model in an appropriate form and design interface to edit the property maps as is done a Semantic ODE.
The easiest way model multiple trajectories (or multidimensional trajectories) is to model each trajectory dimension independently but conditioned on the multidimensional initial condition 2p e RM, i.e., for M trajectories, want to fit M forecasting models, F1,..., FMr such that for each
42

---

Published as a conference paper at ICLR 2025
m [M], Fm R- C². Each F.. can still be represented as Fraj Fllm) where Finj is the same described in the paper and scm R C P. sem can similarly be represented using a composition map and property maps for each predicted composition. A proof of concept of this approach is demonstrated in Appendix B.5 This approach should scale to multiple dimensions as is very modular. Each Frem can be analyzed independently. Each property map of F corresponding to a particular composition can be analyzed independently. And finally, each shape function for each of the predicted properties can be analyzed independently. We believe this modularity is one of the reasons why even systems with multiple dimensions can remain understandable and, more importantly, can be edited each change has a localized impact. In contrast, changing a single parameter in a system of ODEs may result in a change of all the trajectories in ways that may be difficult to predict without careful analysis.
Another extension would be to model some dimensions jointly. For instance, we could try to characterize the shape of awo-dimensional curve. This would require extending the current definition of semantic representation.
Computational cost The computational cost of extending semantic modeling to M dimensions depends on the exact approach taken. We will try to do our best guess assuming the approach remains as similar as possible to the one we describe in the paper.
With our current ideas, the added cost from having multiple dimensions comes mostly from the need to perform the same tasks M times. For M trajectories, we need to fit M composition maps and M property maps.
We do not expect the time to a single composition map to increase significantly. The main computational burden of the current implementation comes from fitting a composition to each sample to see how well it fits. For the same number of samples, the computation cost is going to be the same. Although there will be added cost find optimal decision boundaries in multiple dimensions (as opposed to one), we do not think this will be significant burden given current efficient classification algorithms.
Fitting of each property map will require M times the number of parameters. However, it is important to note that in most cases, the current property maps have just tens of parameters, this number is still likely to be a few hundred at most. It is possible that the number of property maps needed be fitted (for each of the compositions) is going to be larger. However, user would still like to narrow it down to reasonable number (likely less than 10) to make sure that the composition map is understandable.
Overall, we suspect the training time to increase at least linearly with the number of dimensions, but we do not expect any combinatorial explosion as often happens in symbolic regression, so for the number of trajectories we may be interested (probably less than 10), the computational costs should not be a major concern.
G.3 GRANULARITY OF INFORMATION
We believe the granularity of information (in terms of the shape of the trajectory, its maxima, minima, and asymptotic behavior) should be comparable between direct semantic modeling and other approaches, long as all approaches are framed as solutions to finding forecasting model. The key difference is that direct semantic modeling gives immediate access to the semantic representation of the model without any further analysis and allows us easily edit the model as well as incorporate semantic inductive biases during training.
There is, however, one way in which ODEs may offer more granular information by assuming a particular causal structure. Note: it may not be correct, and it makes them less flexible. For discretized ODE we can write the underlying causal graph as z(to) + x(to + At) + + (t). z(t). This generality is useful (as we show in Section 6.4) and it does not negatively impact the performance. However, this model does not allow us do interventions where change the value of x at particular t and predict the future. ODEs do allow for this kind of inference but, of course, there is no guarantee that this inference is correct. For instance, a system governed by delayed differential equation has a different causal structure.
Whereas model assumes more general r(66)
43

---

Published as a conference paper at ICLR 2025
We believe direct semantic modeling in the future be extended to accommodate some kinds of intervention, where they can be described as additional features beyond the initial conditions.
G.4 PARTIAL DIFFERENTIAL EQUATIONS
Partial differential equations (PDEs) longer describe trajectories function defined on the real line) but rather fields which are defined on multidimensional surfaces such as planes. Extending direct semantic modeling to PDEs requires developing new definition of the semantic representation of field. Moreover, as the initial condition is no longer single number, the semantic predictor needs to take whole functions as inputs (for both initial and boundary conditions).
G.5 FLEXIBILITY OF MOTIF-BASER SEMANTIC REPRESENTATION
In this section, we will prove the following representation theorem.
Theorem 3. Any function E C² (lo, +o0) whose second derivative vanishes at a finite number of points (i.e.. Kt : H(t) = 0H|< oc) can be described as a finite composition constructed from proposed motifs.
Proof. We call all points in {t : H(t) 0} the vanishing points. Let dnrst and bast be the smallest and largest such that #(t) I is positive or negative for all points in between consecutive vanishing points well as for all t < first and > Hlst-- That means
As ø vanishes only finite number of times, such numbers exist.
either strictly increasing or strictly decreasing in between consecutive vanishing points as well for all < lisa and f >> Hast That means there is, at most, one point between consecutive vanishing points where 1 = O. Similarly for t iirst and >> Hast That means there is finite number of points where è vanishes. we call such points extrema and denote the smallest and largest extremum as Giist and Clas- For all points between consecutive extrema is either positive negative (similarly for f < Cúiss and t Cias)). That means that is either strictly increasing or decreasing between consecutive extrema (similarly for < éûrs and > Elas). We can now take the union of all vanishing points and extrema and divide (4o, +00) into intervals such that on each interval, æ either strictly increasing or decreasing and is either convex or concave. Thus we can assign to each bounded interval one of the proposed bounded motifs. In the final (unbounded) interval, x can be strictly increasing and convex ($++u) or strictly decreasing and concave (s--w). is strictly increasing and concave, then either approaches real number (s+-h) or diverges to infinity -s+ If it is strictly decreasing and convex, then either it approaches a real number ($=+h) or diverges to infinity (8=+a). Thus we are able to describe the whole trajectory as a finite composition constructed from proposed motifs.
G.6 SHARING SAMPLES BETWEEN PROPERTY MAPS
Each Foro is trained separately for each composition. They do not share any information. However, believe some amount of sharing may improve performance where the two adjacent compositions are similar, and the exact boundary is uncertain. We consider this an interesting extension for future work.
44


Reviews:
[{"review_id": "LwbbDqVnN8", "reviewer": "Reviewer_rHqF", "summary": "This paper presents an approach to modeling the solutions (trajectories) of dynamical systems by bypassing the need to derive (discover) the closed-form of ordinary differential equations (ODEs). Instead of the conventional two-step approach (discovering an ODE and then analyzing it), the authors propose *direct semantic modeling*. \n\nThe promised contributions include:\n\n1. A shift from Syntactic to Semantic representations. The latter delivers a formal means to describing the behavior of trajectories and system dynamics. \n\n2. A direct Semantic Modeling framework, which essentially comprises a *trajectory predictor* that generates trajectories fitting the semantic description. This model alleviates the need for attaching a compact symbolic equation, and employs so-called \"motifs\" to describe trajectories, while remaining interpretable.\n\n3. Through a set of mostly 1D experiments, including a pharmacokinetic model, the authors demonstrate that the proposed semantic ODE offers improved performance, interpretability, and robustness to noise compared to traditional models like SINDy.", "questions": "It would be helpful if the authors reply to the following questions, which are linked to the previously identified weaknesses.\n\n1. How flexible is the proposed motif-based representation in accommodating unforeseen trajectory shapes or behaviors? What kind of motifs would the authors suggest adding if new trajectory patterns arise that fall outside the current motif set?\n\n2. Extending to higher-dimensional systems appears non-trivial within the current approach. Do the authors have any specific strategies or theoretical insights on how they might define and handle multi-dimensional semantic representations and motifs?\n\n3. Why did the authors choose not to include comparisons with neural ODEs or hybrid methods that combine neural networks with symbolic regression, which have become prominent in modeling complex systems? Could these methods offer competitive advantages in interpretability or performance?\n\n4. Do the authors believe that the proposed scheme offers the same granularity of information as i) equation discovery scheme and ii) solution discovery schemes making use of specific analytical functions? If so, in what sense?\n\n4. Since the model relies on motifs for capturing trajectory shapes, how would it perform in the presence of chaotic or unpredictable system dynamics? Have the authors tested the model’s performance on chaotic systems, and if so, what were the results?\nIf possible, testing on a few chaotic systems or adding a discussion on expected limitations in handling chaotic behaviors could clarify the approach’s robustness and generalizability. If chaotic dynamics prove problematic, the authors could suggest potential modifications for better handling non-smooth or irregular trajectories.\n\n5. What specific adjustments or extensions do the authors envision to broaden the adoption of Semantic ODE? Are there plans for enhancing the method to accommodate periodic or oscillatory behaviors that cannot be fully captured within finite compositions?", "limitations": "", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 4, "strengths": "The paper puts forth a scheme that moves away from the two-step modeling process of equation discovery and analysis, proposing instead that the semantic (behavioral) representation of a system can be derived directly from data, which in essence renders it closer to a solution discovery (as opposed to) equation discovery scheme. The concept of focusing on semantic rather than syntactic representation is original, as it allows modelers to interact with and interpret the behavior of systems without delving into complex symbolic mathematics. The originality here lies in both the reframing of the problem and the technical solution that bypasses the need for closed-form equations.\n\nThe paper demonstrates a clear and well presented theoretical development and practical implementation. Key symbols are defined in a notation table and important terms such as semantic representation, syntactic representation, and direct semantic modeling are well-defined, with illustrative examples and figures that clarify the difference between traditional ODE modeling and the proposed approach. The evaluation experiments are cover a range of dynamical systems (albeit low dimensional), including pharmacokinetic models, systems with delay differential equations, and integro-differential equations. The authors further validate the practical usability of Semantic ODE, with case studies and error analyses that underscore its flexibility, robustness, and potential for adaptation. \n\nThe potential significance of this paper is substantial. By removing the reliance on closed-form equations and focusing on the semantic representation of system behaviors, the proposed method has applications across multiple domains, where dynamical systems modeling is essential. Additionally, the demonstrated robustness to noise and the flexibility of Semantic ODE to adapt to various types of differential equations indicate that this method could be robust, at least for the type of low dimensional problems that are here treated.", "weaknesses": "Despite the interesting idea, which draws some ties to a (flexible) curve fitting approach, certain weaknesses are identified as follows:\n\n1. **Limitations Due to Motif-Based Semantic Representation**: While the paper argues that semantic modeling offers flexibility over traditional ODE discovery methods, the reliance on a predefined set of motifs to construct semantic representations introduces notable constraints. The specific motifs used limit the types of trajectories that can be effectively modeled, potentially narrowing the applicability of the method to systems with well-defined, predictable patterns. This approach risks excluding complex or novel trajectory shapes that fall outside the pre-selected motif set, thereby introducing a form of inductive bias that contradicts the method’s stated goal of flexibility.\n\n2. **Scalability and Extension to Higher Dimensions**: The current implementation focuses exclusively on one-dimensional systems, and while the authors mention potential future work on extending to higher-dimensional systems, they do not provide concrete strategies for achieving this. The direct semantic modeling approach, as presented, lacks a clear path for scaling motifs and trajectory predictions in multiple dimensions, which is crucial for broadening applicability in fields like biology, physics, and engineering where systems often involve multiple interacting variables.\n\n3. **Evaluation Against More Complex Baselines**: While the authors compare Semantic ODE to methods like SINDy, the experiments could be strengthened by benchmarking against more advanced and recent models in dynamical system solution discovery. Specifically, neural ODEs (e.g., Chen et al., 2018) and hybrid models that combine symbolic and neural approaches (like symbolic regression combined with neural networks) could serve as competitive baselines, as these methods have demonstrated flexibility and scalability, even in higher dimensions.  Do take note of the fact that such schemes as well, e.g. Neural ODEs, are linked to interpretability due to their ties to known dynamical system forms (e.g. state-state ODE formulation)\n\n4. **Conceptual Contribution**: At the end of the day, it is not so clear why using predefined motifs is really different to assumptions of classes of solutions (as SINDy does) that sever as a bases to define convexity or concavity (which is really what the motifs do). One should also keep in mind that equations discovery is not always a necessary step. It is, however, a very useful step for understanding the type of physics that is underlying the problem. How does the proposed scheme offer the same granularity of information?\n\n5. **Robustness to Unpredictable System Behaviors**: While the authors demonstrate that Semantic ODE is robust to noise, they do not address its performance in systems with unpredictable or highly chaotic behaviors. Since the semantic model relies on motif sequences, unpredictable behaviors could challenge the motif-based representation and disrupt accurate prediction or interpretation. Why not try this for example on nonlinear systems, such as the Duffing oscillator.", "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": "", "extra_scores": {}}, {"review_id": "3kU6i4wKwf", "reviewer": "Reviewer_DLBF", "summary": "This work compares syntactic and semantic representations of ordinary differential equations. In the second representation, the authors introduce direct semantic modeling via predefined bounded, unbounded, and horizontal asymptote dynamic motifs. This data-driven method is demonstrated through a one-dimensional ordinary differential equation. A semantic predictor is trained via a classification-base composition predictor and a properties prediction via a set of univariate functions, and fed into a trajectory predictor, which outperforms the traditional two-step process of equation discovery and analysis in terms of prior knowledge injection, model editing, feedback mechanism, and noise resistance. The syntactic inductive bias due to the analytic and human-understandable nature of equations can be mitigated and the semantic inductive bias can be cured by editing the model.", "questions": "How to technically treat singular points in a solution via the discrete representation since the determination of the positions of these points depend on the numerical grids? Does this process limit the model performance for solutions with a high density of singular points or short-wavelength features?", "limitations": "", "rating": 6, "confidence": 5, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "Compared to the two-step modeling approach exemplified by recently developed equation discovery models, the current work reports a data-driven model where the one-dimensional ordinary differential equation is segmented via transition points, and fitted by semantic motifs following the idea originally proposed in Kacprzyk 2023. The so-called direct semantic modeling approach allows more complexities beyond the analytical models and can integrate empirical knowledge via model editing. The definition of the problem and solution properties including the features of dynamic motifs are well formulated. Although the logic, workflow, and performance of the semantic predictor is demonstrated by using a one-dimensional model, the idea extends to higher dimensions.", "weaknesses": "Although in the discussion and appendix G.2, the authors proposed a roadmap for high-dimensional applications of the current model. However, there needs some elaborated discussion on the computational complexity and cost to achieve that. Even for the comparison made with the two-step model in the main text, the computational costs especially in problems with rising complexities (e.g. singularities) are not discussed.\nThe caption of Figure 1 is too brief, which makes the reading below very difficult (actually all figure captions are brief). The writing could be improved, for example, in lines 257-258, t_{i} and x(t_{i}) are not defined in the previous definition of p_{x}. Some content such as Fig. 10 may need to be moved to the main text since this is the key innovation in the current work. Fig. 7 may need to be merged into Fig. 6 for a better illustration of the improvement.", "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": "", "extra_scores": {}}, {"review_id": "y7Cs6USAR0", "reviewer": "Reviewer_DbZG", "summary": "This paper proposes a direct semantic modeling method, which is fundamentally different from the traditional two-step modeling, which first discovers the symbolic form of ODEs. This approach is evaluated through extensive experiments.", "questions": "1. How to ensure that the training subset of $F_{prop}$, which coresponds to one composition, have enough data? lf the dataset is divided into many subsets because of maybe many different intervals or different motifs, can $F fpropl$ be trained well with small datasets?\n2. Do different $F {propl$s of different compositions share weights? Or do you need to train as many $F fpropl$s as the compositions?\n3. I wonder if training with $F^0_{traj}$ but inference with $F^1_{traj}$ will introduce misalignment and error of the output.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 2, "presentation": 3, "contribution": 4, "strengths": "This algorithm is completely new and different from the previous works. The authors provide detailed definitions of many new concepts and provide many analyses of whether they are well-defined or not. With experiments on different tasks and equations, results show that this algorithm can outperform previous algorithms in many ways.", "weaknesses": "1. Note the capitalization of the title.\n2. Except for the traditional two-step modeling methods, there are also direct modeling methods for ODE trajectory simulation, like FNO, DeepONet, and many others. Considering the comparisons between them and Semantic ODE will be meaningful.\n3. This algorithm is quite novel in that it only provides semantics rather than closed-form expressions, but can you analyze more about what kinds of downstream tasks it has?\n4. l86, l475: Lack of a period at the end of the title.\n5.1210: (denoted $c x$)", "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": "", "extra_scores": {}}, {"review_id": "93xIeSG9Oo", "reviewer": "Reviewer_BR8f", "summary": "The authors propose direct semantic modeling of low-dimensional dynamical systems, which shifts away from the traditional two-step pipeline composed of first discovering a closed-form equation and then analyzing it.  The framework allows the incorporation of intuitive inductive biases into the optimization algorithm and the editing of the model’s behavior directly from data. The effectiveness of this framework is validated by extensive experiments including the logistic growth model, systems governed by a general differential equation, pharmacokinetic model, delay differential equation, and the integro-differential equation.", "questions": "1. Please elaborate on the advantages and disadvantages of direct semantic modeling compared to closed-form solutions in terms of model interpretability and generalizability.\n\n2. If possible, please demonstrate the model’s effectiveness on more complex data, such as the semantic representation learning and prediction of time series in experiments.\n\n3. How to efficiently identify bifurcation points, perform semantic modeling, and predict bifurcation behavior within dynamical systems using this framework when bifurcation behavior is present？\n\n4. If possible, the author could include a discussion on handling partial differential equations with this framework, noting that in addition to initial value problems, there are also boundary value problems.\n\n5. By the way, there are some minor typing mistakes in the article. On the second page, the differential operator symbol ‘d’ should be set in upright type.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "1. A novel approach, called direct semantic modeling, is proposed to learn the semantic representation of the dynamical system directly from data, eliminating the need for post-mathematical analysis.\n\n2. The method operates directly on the semantic representation, enables intuitive adjustments, and seamlessly incorporates constraints that mirror the system’s operational behavior, which enhances the modeling flexibility and boosts performance, as it bypasses the need for a concise closed-form equation.", "weaknesses": "1. Compared to methods that rely on closed-form solutions, the interpretability and generalizability of direct semantic modeling approaches appear to be less apparent. As the author mentioned in the main text: \" The primary objective of discovering a closed-form ODE, as opposed to using a black-box model, is to have a model representation that can be analyzed by humans to understand its behavior. Such understanding is necessary to ensure that the model behaves as expected.\".\n\n2. The author has showcased the method’s efficacy through a variety of examples of 1D ODE, such as the logistic growth model, systems governed by a general differential equation, pharmacokinetic model, delay differential equation, and the integro-differential equation, which is commendable. However, the lack of demonstration of more complex data such as experimental data with unknown dynamical equations may diminish its appeal to a broader audience and contribution.\n\n3. Identifying transition points (critical points) in dynamical systems that exhibit bifurcation behavior is inherently challenging, which may constrain the application of the framework to such systems, as their semantic representations tend to be more complex.", "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": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)