Skip the Equations: Learning Behavior of Personalized Dynamical Systems Directly From Data
Authors: Krzysztof Kacprzyk
Julianna Piskorz
Mihaela van der Schaar
Venue: ICML 2025 poster (2025)
Keywords: dynamical systems
differential equations
ODE discovery
interpretability
transparency
verification
Primary area: applications->chemistry_physics_and_earth_sciences
TL;DR: An interpretable framework for trajectory prediction from static features: a decision tree identifies the curve’s global shape, while generalized additive models assign the key properties such as start, peak, or inflection.

Abstract: While black-box approaches are commonly used for data-driven modeling of dynamical systems, they often obscure a system's underlying behavior and properties, limiting adoption in areas such as medicine and pharmacology. A two-step process of discovering ordinary differential equations (ODEs) and their subsequent mathematical analysis can yield insights into the system's dynamics. However, this analysis may be infeasible for complex equations, and refining the ODE to meet certain behavioral requirements can be challenging. Direct semantic modeling has recently been proposed to address these issues by predicting the system's behavior, such as the trajectory's shape, directly from data, bypassing post-hoc mathematical analysis. In this work, we extend the original instantiation, limited to one-dimensional trajectories and inputs, to accommodate multi-dimensional trajectories with additional personalization, allowing evolution to depend on auxiliary static features (e.g., patient covariates). In a series of experiments, we show how our approach enables practitioners to integrate prior knowledge, understand the dynamics, ensure desired behaviors, and revise the model when necessary.

Skip the Equations: Learning Behavior of Personalized Dynamical Systems Directly From Data
Krzysztof Kacprzyk'
Julianna Piskorz Mihaela van der Schaar'
Abstract
ied, for instance, in bifurcation theory (Blanchard al., 2012)), but we emphasize the fact that we usually have many such parameters, we can observe them, and they are often linked to personal characteristics of individual peo- ple, A primary example of a personalized dynamical system (PDS) the human body, which exhibits lot of variability when it comes to drug metabolisms, disease progressions, or treatment effects. Beyond physiological systems, other examples of PDS may involve the atmosphere (e.g., tem- perature or pollution) based on particular location or the animal population in a particular environment.
While black-box approaches are commonly used for data-driven modeling of dynamical systems, they often obscure system's underlying behav- ior and properties, limiting adoption in areas such as medicine and pharmacology. A two-step pro- cess of discovering ordinary differential equations (ODEs) and their subsequent mathematical anal- ysis can yield insights into the system's dynam- ics. However, this analysis may be infeasible for complex equations, and refining the ODE meet certain behavioral requirements can be challeng- ing. Direct semantic modeling has recently been proposed to address these issues by predicting the system's is behavior, such as the trajectory's 's shape, directly from data, bypassing post-hoc mathe- matical analysis. In this work, we extend the original instantiation, limited to one-dimensional trajectories and inputs, to accommodate multi- dimensional trajectories with additional personal- ization, allowing evolution to depend on auxiliary static features (e.g., patient covariates). In se- rics of experiments, we show how our approach enables practitioners to integrate prior knowledge, understand the dynamics, ensure desired behav- iors, and revise the model when necessary.
Discovery of Closed-Form ODEs In practice, it is of- ten essential to ensure the model's behavior aligns with domain requirements. For instance, when develop- ing drugs, researchers must verify that the pharmacokinetic model (Mould & Upton, 2012) is biologically plausible (e.g., ensuring non-negativity and eventual decay of drug concentrations). Moreover, dosing recommendations may depend on the drug's maximum concentration and its timing (Han et al., 2018). 8). An effective way to achieve this level of understanding is through the discovery of closed-form ordinary differential equations (ODEs) (Bongard & Lipson, 2007), where an algorithm identifies a concise mathemat- ical representation that can then be analyzed by experts. More formally, closed-form function is a function that can be represented by mathematical expression consisting of finite number of variables, constants, arithmetic opera- tions (+, +), and some well-known functions such as trigonometric functions. Initial approaches to ODE discov- ery (Bongard & Lipson, 2007; Schmidt & Lipson, 2009) employed genetic programming (Koza, 1992) to search the space of possible equations. More recently, ODEs are often represented linear combinations of terms from a pre- specified library (Brunton et al., 2016b) which makes the search much more efficient. This approach (generally called SINDy) has been extended to various settings, including implicit equations (Kaheman et al., 2020), equations with control (Brunton et al., 2016a), and partial diffcrential equa- tions (Rudy al., 2017). The extended discussion about related works can be found in Appendix E.
1. Introduction
Personalized Dynamical Systems In machine learning (MI.), modeling dynamical systems is core activity with extensive applications in fields such as physics (Rudy et al., 2017), biology (Chen et al., 2019), engineering (Brunton & Kutz, 2022), and medicine (Lee et al., 2020). We call a dynamical system personalized when its evolution fune- tion depends on certain observable features. It is similar to dynamical systems that depend on a parameter (as stud-
Universsity of Cambridge. Cambridge, UK. Correspondence to: Krzysztof Kacprzyk <kk751 @cam.ac.uk>.
Proceedings of the 42nd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).
Limitations of Two-Step Modeling Traditionally, to un- derstand dynamical system's 's behavior, an ODE needs to

---

Skip the Equations: Learning Behavior of Personalized Dynamical Systems Directly From Data
Data
EPISODE
Equation discovery
Xmax
@Cdepot dt. dCwansi dc dGransz de dCranna dt dGrat dt
-ki, Gdupps
X
Kisfidepot
Ftranst
x(t) = t++0000
lim
Analysis
KepCeranst Kepherancz
maa
KepCaranez KepCerance
(CL. + Crentt
QCeert V22
Ky Crrnnsa
V/
Challenging edit
pari QGent
Cpert
dt V2
V22
Syntactic representation
Semantic representation
(equations)
(behavior)
Figure 1. In two-step modeling, to learn the behavior of dynamical system, need first discover closed form goveming equations and then perform a time-consuming analysis to arrive at the semantic representation of the model. EPISODE allows us to skip equation discovery and prodict the semantic representation direetly from data. The figure shows part of EPISODE trained on the Tacrolimus datasct (see Section 5). where we know that the trajectory will always have the shown shape, and its properties, such as the maximum of the trajectory, described generalized additive models. Each of the individual shape functions can be scrutinized and edited if necessary. allowing for verification and enforcement of certain conditions, such as lim; xrt) = O.
be first discovered and then analyzed by humans (two-step modeling). This process has many limitations. It is time- intensive, demands substantial mathematical expertise, and can be intractable for complex equations. Moreover, when the discovered ODE fails to meet certain behavioral require- ments, revising is difficult because the link between the form of the equation (its syntactic representation and its behavior (its semantic representation) is rarely simple.
tion 3, we describe the main ingredients of model: the composition map implemented as decision tree, and the property map constructed from a set of generalized additive models. This complex architecture requires a novel training algorithm, which we describe Section 4. In Section we present case study of how EPISODE can discover novel population pharmacokinetic model with better perfor- mance than expert-designed model from the literature. We demonstrate how our approach enables practitioners to integrate prior knowledge, understand the dynamics, verify the model's properties, and revise if necessary. Finally, we demonstrate the flexibility of our approach series of experiments Section
Direct Semantic Modeling To overcome these limita- tions, direct semantic modeling (DSM), has recently been proposed (Kacprzyk& & van der Schaar, 2025). It predicts system's semantic representation (its behavior) directly from data, bypassing post-hoc analysis and enabling greater flexibility. This allows for more intuitive model editing and incorporation of prior knowledge. However, the current instantiation of direct semantic modeling, Semantic ODE, limited one-dimensional scenarios without personaliza- tion (apart from the initial condition).
2. Problem Formalism and Background
In this section, we formulate the problem of time series forecasting from static features which provides formal- ism for modeling PDS that generalizes ODEs. Then, we deseribe the principles of DSM and the theory behind Se- mantic ODEs on which we build in Section 3.
Contributions and Outline In this paper, we develop EPISODE (Editing and Personalization In Semantic Ordinary Differential Equations) that applies the principles of DSM framework to modeling multi-diimnssioal PDS (Figure 1). In Section 2, we formalize modeling PDS as time series forecasting from static features (extending ODE discovery), describe the principles of DSM, and present the theoretical background necessary for our approach. In Sec-
Notation and Terminology Let tend E 14 U {+oo} and T (to 0. fend) be the time domain. Let M N. We say æ T + RM is an M-dimensional trujectory if for each m E [M] (where [M] = {l..... ........) M),, Xm E CT)(T), i.c. Fm is twice continuously differentiable function on T. We also denote it as æ E C²). We call any ML model whose output is trajectory a forecasting model.
2

---

Skip the Equations: Learning Behavior of Personalized Dynamical Systems Directly From Data
2.1. Time Series Forecasting From Static Features
Forecasting Using Semantic and Trajectory Predictors We assume that the main reason for discovering a closed- form ODE is to obtain a model with a concise syntactic representation (symbolic form) that can be subsequently analyzed to obtain the semantic representation of the under- lying forecasting model deser--dseirittion of how trajectory's be- havior, such as its shape and properties, changes as we vary the input. Recently, DSM has been proposed as an alterna- tive to this two-step process for modeling low-dimensional dynamical systems (Kacprzyk & van der Schaar, 2025), Instead of discovering an equation from data and then ana- lyzing it to obtain its semantic representation, this approach generates the semantic representation directly from data. In general the forecasting model F defined as composition F = Fry Frem- Fsem directly predicts the semantic rep- resentation of the trajectory (description of its behaviors), which is passed to the trajectory predictor Fus that matches this description to trajectories in given hypothesis space. No post-hoc mathematical analysis is required as the se- mantic representation of F can be directly accessed through Fsom The model can be easily edited to enforce a specific change in the semantic representation. Incorporating prior knowledge is also streamlined and more intuitive.
Definition We choose to formulate modeling PDS as time series forecasting from static features. In this setting, we have a dataset D = {v(d) (the) yn )))) I where each (d)
yn E R represents a noisy measurement of some ground
truth trajectory ædd) E (T) at time point (.4)) E T. and E V C TR. are the static features. The goal is to find a forecasting function F: V - (T) that predicts an M -ddimensional trajectory x T - R M from the static features E V. These features may include things like the patient's s covariates, the type and dose of the treatment, the city's location, or initial values of any of the M trajectories.
Relationship with ODEs This setting accommodates and extends the dominant approach to modeling dynamical sys- tems: ODEs. Each system of ODEs f defines forecasting model FODE through the initial value problem (IVP), i.e., for each initial condition =0) Io E RM FODE maps 2o to M trajectories governed by f satisfying this initial condi- tion. Thus, ODE discovery can be treated as a special case of fitting a forecasting model when the static features contain only the initial measurements, i.e., v = zo. Importantly, time series forecasting from static features also allows for trajectories that not governed by fully-observed ODEs or which are only defined for bounded domain (fend F +oc). It is also useful when we want to predict the whole tra- jectory without access to the first observation (v does not contain æo), and where we have more variability than can be captured solely by the initial condition.' 1
2.2.1. SEMANTIC REPRESENTATION OF TRAJECTORY
To realize DSM, we need to formalize the semantic represen- tation of a trajectory that allows us to uniquely match it to an actual trajectory through trajectory predictor Fraj. Seman- tic ODEs employ such a definition for a one-dimensional trajectory. We will introduce it below as it will form a starting point for our own definition. The semantic represen- tation of trajectory x: : (0, femi) + R (note that the original formulation only allows for fend = +00, so we extend it here) is described as a pair (cz, or is the composition of x. i.c., its overall shape and Pa is the set of properties describing this shape.
2.2. Direct Semantic Modeling
The discovery of closed-form ODEs is a common way of modeling dynamical systems in a way that potentially al- lows us to understand the underlying behavior of the system through a careful analysis of the found mathematical equa- tion. However, not only do ODEs not cover all forecasting settings, but the two-step process of first discovering equation and then analyzing it has many limitations. A framework of DSM has recently been proposed to address these limitations when fitting forecasting models. In our work, we want to apply this framework to the particular problem of time series forecasting from static features. First, outline the difference between syntactic and semantic representation. Then, we explain the main principle of pred-cccittion using a semantic predictor. Finally, we introduce the definition of semantic representation used in Semantic ODEs, which we use as a starting point in our method.
Composition More formally, a composition e, is a se- quence of motifs where each motif deseribes the shape of the trajectory on particular interval, for instance "increas- ing and strictly convex". We have four bounded motifs and six unbounded motifs. Each of them is of the form $i i.c., is described by two symbols (cach + or -) and one letter from (b, u, h). The symbols describe the signs of, respectively, the first and second derivatives. The letter b denotes a motif for a bounded time domains, whereas both letter and u are used for unbounded time domains. The letter h is used for motifs that have horizontal asymptote. For instance, 8- is an unbounded motif describing a func- tion that is increasing (+), strictly concave (-), and with a horizontal asymptote (h). An unbounded motif is always the last motif in the composition (if fend = +oo). We also call such composition unbounded. Compositions without an
'ee note, however, that this last point could be accommodated by discovering ODEs with control and treating static features as constant control trajectories.
3

---

Skip the Equations: Learning Behavior of Personalized Dynamical Systems Directly From Data
unbounded motif (i.e., when lend E R) are called bounded. We call set of all possible compositions C.
Each Dimension Separately We propose to treat each dimension of the trajectory separately. As all information needed to predict the whole M-dimensional trajectory is encoded in the static features v, in particular v is sufficient to predict any of the M dimensions. Our model does not need to take the current values of the remaining dimen- sions as input, which significantly simplifies the modeling (and understanding) process. That means that our fore- casting model F Funj Fiem consists of M submodels { ('m m== each taking E V and predicting a univariate trajectory : ' R. Then each Fm. can be defined as Fm F mij F(m) where Funj is the original trajectory pre- dictor used in Semantic ODEs and F(m) is the new semantic predictor that takes v E V as input and outputs a semantic representation of the m. dimension of the trajectory, i.e., (Ce.. Pr.... ). where (cz. P2) is the semantic representation of dieediinniioall trajectory as defined in Section 2.2.1. As for each trajectory dimension m E [M], we train separate semantic predictor Fdmm) Throughout the rest of this section, we fix m and drop the superscript m to improve readability. All that follows describes only one of M semantic predic- tors. In contrast the semantic predictor in Semantic ODEs, each Fsem takes a multidimensional vector v E V input instead of a single number. Thus, the comprchensibility of Fsem is in no way guaranteed and requires careful choice of the underlying architecture. As in Semantic ODEs, we describe F.em as a pair (Feom: Frou),, where the former is a composition map and the latter a property map. We describe them in the next sections. A block diagram showing how a single semantic predictor works is shown in Appendix C.
Properties Ps is the set of properties describing the com- position Cz. This includes the coordinates of the transition points (points between two consecutive motifs, the point at t = 0 and at fend if fend E IR). In addition, it also con- tains information on the first derivative at the first transi- tion point and the first and second derivative at the last transition point. Finally, if fend +00, the set of prop- erties also contains the properties of the unbounded mo- tif. For instance, one of the properties of motif 4 is the value of the horizontal asymptote. Given a trajectory x with a composition (s--b. +i)), its set of properties is Ps = (ln, (lo); (1 (())) (t1)), ((11), J, (1/2),, where each x(t)) is a transition point, and (h, t1/2) are the properties of sh+ We denote the set of all proper- ties for composition as Pe and define P = Ucec Pe. An example of trajectory and its semantic representation can be seen in Figure 2.
Semantic representation of x
Composition: (B S.+.)
Trensition points:
x
(0.0.3).(0.3.0.55..
(0.5.09
Derivatives:
((n)) - 0.2,x(t.) = 55,
x(r))
(8--h. S4-b. 5_ 2
composition
of
motif
h = 0,0 0.9
3.1. Composition Map
Figure2 2. An example of trajectory and its semantic representa- tion. The semantic representation describes the composition of the trajectory, (8++b, 84 8--,, 8-++)) the coordinates of the transition points, some of the derivatives, and properties of the unbounded motif s
The goal of the composition map From to predict one part of the semantic representation of the trajectory, namely the composition C', where C' c C is set of all admissible compositions chosen by the user. We can interpret F From as a classification algorithm that takes the static features E V and outputs one of the classes in C. To make the composi- tion map as transparent as possible, we implement as a classification decision tree of small depth (Breiman, 2017). This allows us understand how the shape (composition) of the trajectory depends on the static features. An example of such a composition map can be seen in Figure 3, where we show a trained composition map on the Tumor dataset (details in Appendix 1. Each internal node is a simple condition checking if one of the features is above or below a certain threshold. The leaves (terminal nodes without children) predict one of the admissible compositions.
3. EPISODE
The goal of this work is to realize the DSM framework for the problem of time series forecasting from static fea- tures and thus extending the original instantiation, Seman- tic ODEs, that only works for -nediimnsioonll dynamical systems where the only input is the initial condition (i.e.,
M = 1, v = zo). In a setting where the input is one- dimensional, every property is described by a univariate function that can usually be plotted and understood directly through its graph. For multivariate inputs, the task is much harder as we need to employ and adapt transparent ML models for multivariate regression and classification. We first discuss the high-level structure of our model and then describe its two main ingredients: multivariate composition and property maps.
3.2. Property Map
Property Sub-Maps The goal of the property map Fpop : V - P describe the properties of the particular compo-
4

---

Skip the Equations: Learning Behavior of Personalized Dynamical Systems Directly From Data
Below each plot, we show the marginal distribution of the corresponding feature in the training dataset. This gives an indication of the confidence in each segment of the function.
Dose
s 0.19
> 0.19
Dose
0.050
0.050
Shape funzton
Shops fundttn
D.125
0.025
(Seau)
S 0.48
> 0.48
0.000
0.000
0.025
0.025
Triquenty
20
60
80
60
Figure 3. Example of composition map found by our algorithm. Not only did our method find a dosc of the drug to be the most im- portant feature determining the trajectory behavior, hut also deter- mined the thresholds where the shape of the trajectory changes. For low doses, the tumor grows exponentially (s+4... for higher doses, it decreases at the beginning but then relapses (s (s-4b, 84-4), and finally, for very high doses the tumor is predicted to decay (*-+))-
0.1
0.025
Shaoe
-0.025
-01
Trequenccy
Tresueeryy
0.2 int lumor volume
0.3 0.4
sition of the trajectory. As our trajectory may have a differ- ent composition depending on the static features v (as pre- dicted by the composition map), the property map consists of few property sub- maps, each for different composition predicted by the composition map. For instance, a compo- sition map shown in Figure 3 requires three property sub- maps. We denote them as Fpiopt" F'poopt and prop" prop
dosage
Figure 4. Shape functions of a GAM describing att1) of the com- position (8-+b, 8++u] (i.c., the minimum of the trajectory) trained on the Tumor dataset. Rach shape function is such that its expected value over the marginal distribution of the feature is O. The bias (expected value) equal to 0.212. We see that age and weight do not impact the minimum tumor volume. The other two func- tions let us observe that the minimum increases linearly with the increase in the initial tumor volume and decreases as the dose of the drug is increased.
(v) = F
where each "roop" V - Pa. Formally (v
Property Functions as GAMs Each property sub- map Pipp : V + Pe is set of property func- tions, each predicting different property from Pa- For instance if
4. Training
need 8 prop- erty functions for each of the property of Pa ((o),, x(41, ((()), ((1),(((1) d). where d is the "asymptotic doubling time". To keep these functions trans- parent, we characterize as many of them as possible as gener- alized additive models (GAMs) (Hastie & Tibshirani, 1986) or simple transformations of GAMs. GAM is fully trans- parent model for a static regression setting where the output is sum of univariate functions (called shape functions), each taking a different feature. It can be easily compre- hended by just plotting all shape functions. As there are no interactions between features, each of the shape functions can be analyzed independently. They proved to be very flexi- ble and transparent models in many real-world settings (Lou et al., 2012; Caruana et al., 2015). Using GAMs, we can, for instance, describe the x-coordinate of the second transition point of (s_ +++ the minimum of the trajectory as a((1)) -K and the individual g 's can be plotted as seen in Figure 4. Note that when we visualize a GAM, we center each shape function such that its expected value over the marginal dis- tribution is O. These contributions get absorbed into the bias term, which equals the expected value (go- = E[)(())).
(s-+b, s++u],
In this section, we describe the training process of a single Fm = Funj a Fim)) i.e., a forecasting model that takes v E V and outputs a one-dimensional trajectory æ. We train each Fm separately for all m [M]. For the remainder of this section, we fix and drop the superscripts (m) to improve readability. As in Semantic ODEs, the trajectory predictor Funj is fixed and not trainable. Its job is to match semantic representation (c, p) E x P to trajectory such that (C2. P2) = (c,p) p). A block diagram depicting the training procedure of a single semantic predictor is shown in Appendix C.
#((1),
Sequential and Parallel Training As described in Sec- tion 3, Frem = (Feom, Fprop). To train Frem we first train From and then train Ppprop. i.c., we train each of the property submaps Prou for all the compositions that are predicted by Feom- After the composition map is fitted, all property sub-maps can be trained in parallel.
k= gk (vx) + 9u. where gu E IR is a bias term
4.1. Training the Composition Map
From is a classification decision tree that takes a vector of static features and predicts the composition c E C².
5

---

Skip the Equations: Learning Behavior of Personalized Dynamical Systems Directly From Data
Matching Compositions to Observed Trajectories We start, similarly to Semantic ODE, by checking how each of the compositions in C" matches each of our D samples. For each sample (v)(4) (th5), yh 1)==)) (where yi is the th
functions (e.g., B-Spline basis functions (de Boor, 1968)).
We define the raw shape function of the hthh feature for the 1th property as 9g = -- Wi,k, where Φb, is 8th basis function for the h.th feature.
coordinate of yn"") and composition E C' we calculate eld)(c defined as
From Raw Shape Functions to Properties We can now get the properties from raw shape functions by appropriate transformations. For instance, if we need to have a positive property function (which is equivalent to having all shape functions and bias positive as we show in Appendix B), it is enough to apply a softplus function, denoted soft, and defined as soft, (z) log(1 [1 e*), each of the raw shape functions. The resulting GAM for the 1th property function would be Σx= soft, (gk (vv)) + soft, (330).. We discuss other kinds of transformations in Appendix C.
Na Σ
(d min pEPe Nu
( - vn vi=t)²
We search the space Pe using the L-BFGS algorithm. To improve the fit, before minimizing the mean squared error (MSE) between the F (c. p) and yd( we first minimize the Dynamic Time Warping distance (Müller, 2007). We ob- served that this procedure performs better, especially when the composition is long.
Fitting tree Fitting classification decision tree requires a dataset of static features and classification labels. In our setting, instead of a single label for each sample, we have whole vector of [C'| real numbers eld), Instead of min- imizing the classification error, we want to minimize the following objective: b. Ed_n e [Faum (FFam((((()) (v( To find decision tree that minimizes such an objective, we perform greedy search where, at each node, we evaluate all possible splits for all possible features (similarly to CART (Breiman, 2017)) and choose the one that gives the smallest error. In practice, we put few additional constraints to prevent a composition map that is unnecessarily more complex with- out a significant improvement in performance. We provide more details in Appendix C.3.
Optimization Before optimization, we precompute (H,)(1)) for all & E [K] and b E [B]. This allows us to efficiently calculate the shape functions as simple matrix products. Then the model is trained to find W and 3 that minimize
D
N₂ (((h))) -
(2)
n=1
where
② == Fing (Fom((u(t)) Fr pmop (vl4); W W.
(3)
We optimize it using the L-BFGS algorithm (Liu & Nocedal, 1989).
4.2. Training Property Sub-Map
Given composition map Feom. we partition our training set depending on the composition predicted by From. Then for each E Im(Paom), we fit separate property sub-map Pinp (d) N (v))((((()) = c). The
5. Case Study: Novel Pharmacokinetic Model for Tacrolimus
on D. = {(v(()) (t(e) Mn property map is defined by weights W E RLLKK xB and biases B E RR,, where B is the number of basis functions describing each of the shape functions of a GAM, and L is the number of properties. To simplify the presentation, we assume that all variables are numerical.
In this study, we demonstrate the main advantages of our proposed approach while developing a novel population pharmacokinetic (PopPK) model for Tacrolimus. That in- cludes incorporating semantic inductive biases (Section 5.2), understanding the dynamics (Section 5.3), verifying the properties (Section 5.4), and revising the model to guaran- tee no unexpected behavior (Section 5.5). For that task, we use a real Tacrolimus dataset (Woillard et al., 201 1) con- taining patient covariates, the dose of the drug, and drug plasma concentration measured over time (for more details, see Appendix D.1). We not only obtain a model perform- ing better than the discovered ODEs and the expert-found PopPK model from the literature, but we can certify that model is biologically plausible.
Property functions need to satisfy many constraints to make sure that the properties predicted by them actually describe a possible trajectory. For instance, if we have a composition (s=+h), we need to make sure that the derivative at is negative. Or =(tu) =(t2) ) if the first motif is increasing. Our goal is to satisfy as many of these constraints as possible by construction. Because of that, we introduce a concept of raw functions. Given L properties, we have K X L raw shape functions. Each of them is represented as a linear combination of basis
Raw Shape Functions
Imppeemnntation of EPISODE is available https:// i thub. https//ccommrrsstoooe-oeo wwwpttoooo-kkkkkkk//kkk ISDDD..
6

---

Skip the Equations: Learning Behavior of Personalized Dynamical Systems Directly From Data
5.1. Tacrolimus and PopPK Models
Tacrolimus is an immunosuppressive drug primarily used to prevent organ rejection after kidney transplantation. PopPK models predict the drug concentration over time for an indi- vidual with particular characteristics. They are crucial for understanding the variability in drug concentrations within a patient population. This information is essential for op- timizing drug dosing to improve drug safety and efficacy. Although black-box models have been used for this task, concerns remain about their lack of transparency and guar- antees of appropriate behavior under every condition. ODE discovery can potentially deliver a transparent model, but getting insights. verifying the properties, and refining the model is challenging, as we illustrate below.
(S++br S4- S.
Figure 5. Composition map found by our algorithm trained on the Tacrolimus dataset.
The equation found by SINDY constrained to 5 terms (SINDy-5) is the following.
+(1) 00.21 - 21.543/ - 6.42x(1))
(4)
- 0.01 1z + 19.065t2
Even though the equation is compact, it is not straightfor- ward to understand the system dynamics without a careful mathematical analysis. Even answering a simple question of whether the concentration will increase at first cannot be answered directly. Comprehending what the rest of the tra- jectory looks like even more challenging and the problem compounds for longer expressions.
5.2. Prior Knowledge
EPISODE allows us to incorporate semantic inductive bi- ases, i.e., prior knowledge we have about the behavior of the dynamical system. This is in contrast to syntactic inductive biases (i.e., information about the structural form of the equations) available in ODE discovery algorithms.
5.4. Verification
PopPK models describe the drug plasma concentration af- ter the drug is administered as a function of time. Thus we expect that the concentration should decay to O as t - 00 (assuming no additional doses). We incorpo- rate this information directly in our model by choosing the composition library C' to only contain compositions where the last motif is s (decreasing, convex, and ap- proaching horizontal asymptote). As also do not expect many trend changes, we constrain our model to the following four compositions: (8=+h), (8--b, (8+-6, 6-, sh),), (8++6: 6,, b-
Expert PopPK models have a structural form inspired by biological processes (e.g., movement of the drug between different organs, natural decay, excretion). This usually guarantees that the models are biologically plausible. This, however, is not guaranteed in data-driven models. For the purpose of this section, we are just going to focus on two conditions that a biologically plausible PopPK model needs to adhere to.
The concentration should initially increase after the drug is administered and then fall as the drug decays.
Designing syntactic inductive biases is more challenging because, in many complex settings, we often do not have much knowledge about the structural form of the equation describing the phenomenon. In PDS like this one, it is far from obvious whether the term (t) or 2 should appear in the final equation.
The concentration should decay to O.
With our approach, verifying these two conditions is straightforward. The composition map of the found model is very simple and assigns the composition (5 84-6. s_th) to all samples. This is precisely the shape we want. The function first increases, then de- creases, and approaches horizontal asymptote. We can also verify whether it decays to O. We can look at the GAM that predicts the horizontal asymptote (all shape functions can be seen in Figure 6 in Appendix A. We can see that our model does not yet satisfy this condition. We will rectify this the next subsection.
5.3. Understanding Dynamics
To understand the dynamics of our dynamical system we need to access the semantic representation of our model. Fortunately, this is readily available in the semantic predictor Fsem- No additional mathematical analysis is necessary. We can see that the composition map very simple as just predicts the composition (8444, 84- 8- +h) for all samples (see Figure 5). We can now be sure that, matter for which patient we use our model, and what drug dose we choose, the trajectory is always going to have correct shape. This, however, is not the only thing we need to verify about our model. We discuss it further in Section 5.4.
Verifying these two conditions for a closed-form ODE is potentially possible, but it may be time-consuming. Let us consider Equation (4) found by SINDY-5. Assume that the two conditions above hold. Then for some small æ there exists t' large enough such that we have i(t) < 0 and 0 S x(t) < T2 for all > t'. Then for > t', i(t) > 0.321 - 21. 21.543 6.422, 0..00...4000 + 19.06512 +
7

---

Skip the Equations: Learning Behavior of Personalized Dynamical Systems Directly From Data
+00 as 1 - +00 So, there has to be > 1 such that
Table I. Results of fitting different models on the real Tacrolimus dataset. Performance on a test dataset. Standard deviation of the performance over the samples in brackets. Error measured using Equation (5). Lower is better.
a(t) > 0, which contradicts our conditions.
5.5. Editing
Editing the model is useful for both improving performance and satisfying some domain-speciff criteria. In the previous section, we established that the horizontal asymptote of our model is not what it should be. We can just take one of the property maps and set it to O. However, this may deteriorate performance, as it was not enforced during training. Instead, we can first encourage the asymptote to be close to and only then set it to 0 exactly. As the semantic predictor directly predicts the semantic representation of the trajectory, we can put soft constraints any aspect of the semantic representations, such as the coordinates of the transition points or the properties of the unbounded motif. In this case, we refit our model after adding penalty h?, where is the value of the horizontal asymptote. After training, we can see that the bias term (the expected value) is -0.0019, which is much closer what wanted (all shape functions can safely
Black box models
Closed-form expressions
NeuralODE ANODE
0.263
SINDy-5 SINDy-20
0.337 (.144)
0.268
0.331
LatentODE
(.281
Direct semantic modeling
WSINDy-5
0.348 (.182)
EPISODE
0.255
WSINDD--20
0.526(,631)
EPISODE-verified
0.257
Expert PopPK
0.351
ODEs and the expert model from the literature. It is also slightly better than unverifiable black-box approaches.
6. Flexibility
Although the most crucial aspects of our method are its transparency, verifiability, and editability, it is instructive whether it can also deliver a good performance in practice.
be seen in Figure in Appendix A.1). Now, set the whole property map to O without a significant drop in performance (see Figure 8). Finally, we implemented a method to prune shape functions whose contribution (after being centered) is smaller than le 5. Part of the model visualized in Figure 1, and the property funetions for the maximum and the time it is achieved can be seen in Figures 9 and 10 in Appendix A.2.
Baselines In this section, we compare our method against two ODE discovery approaches, SINDy (Brunton al., 2016b) and WSINDy (Reinbold al., 2020; Messenger Bortz, 2021a) (as implemented in PySINDY (de Silva et al., 2020), in two different variants: up to 5 terms (compact) and up to 20 terms (not compact). Smaller equations are easier to interpret but may not be as flexible as long expres- sions. We also compare our approach to three black box approaches: Neural ODE (Chen et al., 2018), ANODE (Dupont et al., 2019), and Latent ODE (Rubanova et al., 2019) probabiibaii variant).
5.6. Performance
The performance of different methods, with different con- straints, can be seen in Table 1. We use the following error metric (lower is better).
N₂
EPISODE For EPISODE, we choose a generous set of compositions without much prior knowledge. We consider all compositions up to length 4 (or length 8 for the Bike sharing dataset) except for compositions where two con- secutive transition points are inflection points. This ends up being, respectively, 26 and 34 compositions. To demon- strate the power of how easily we can incorporate prior knowledge into our model, we also compare against our method with inductive bias about the possible compositions of the trajectory. We call this variant EPISODE* The ex- act composition libraries used for each dataset are shown in Appendix D.3. Although the training process consists of two steps (training the composition map and then the property maps), the whole pipeline is end-to-end, and there is no manual intervention in experiments performed in this section.
1 D
Σ NA
1
1 M
F(u(d) (1))0
(d) 2) yr (2]
(5)
Here, F is the predictive model (e.g., EPISODE), D is the number of samples, M is the dimensionality of the system, No is the number of measurements of the th sample (o(d), (152) yn
(d) N=n). We choose this metric because for
M- = 1 it reduces to the standard mean RMSE, i.c..
D
No Σ
(F(old)) (1ttt)
D
(used in Semantic ODE paper), and we normalize by M so that it is easier to compare results between systems of different dimensionalities. The error is calculated on a held- out dataset of samples, whose trajectories are observed at the same time points as the ones in the training dataset. We observe that our method performs better than the found
We test the methods on a variety of datasets. To show that we can model traditional multi-dimensional
Datasets
8

---

Skip the Equations: Learning Behavior of Personalized Dynamical Systems Directly From Data
Table Comparison between different methods tested on several synthetic and real datasets. means that the method has not converged to a solution in at least one trial. The table shows means from 5 random dataset splits and seeds. Standard deviations in brackets, Error measured using Equation (5). Lower is better.
SIR
III
PK
Tumor
Tacrolimus real
Bike sharing
a = 0.01
0.1
o 0.01
= 0.1
0.01
0.1
0.01
0.1
real
Black box models
NeuraIODE
0.021,000
0. ,011
0.188; c09)
0,214,,00 (000)
0.156 oo6)
0.182 005)
0.300, .nna)
(.334
0.267, oss)
0.178,
ANODE
0.017,00m (00))
0.001
0.188; cos)
0,214 0o5)
0.147 (000
0.184 0099
(,40
0.335, (000)
(.250,o01))
0.154
LatentDDE
(.02)
0.101 (300))
0.183
0.145 0 (.203
0.172 .000)
..310
0.338, (92)
0. 0.2600(o1z) 013)
(.149
Closed-forde expressions
SINDy-5
0.014
0.135, a00r
0.176
0.284.00m)
0.211
0.249 .o10)
0.096, oor;
0.151,,
(.310
0.267,.10)
SINDy-20
0.015
(.180
(.062
0.138,
0.3
(.250
WSINDY-5
0.026,oom)
0. (o00))
0.246
0.279
0.309
0.094 cor;
0.145
0 1909
WSINDY-20
0,016
0. 144,000
0.04
0.129(.008)
0.816
0.310
Direct semantic modeling
EPISODE
0.026 (00))
0.088 away
0,143 (007)
0.088
0. 0.132
0.144, s002:
0.250...4)
0.192 (.25)
EPISODE
0.012
0. (001)
0.089.coo .00)
(007)
0.095,,030
0.146
0.098(.014
0.165
0.256g.077) 0007)
0. (.15
systems, we test the models on the SIR epidemiological model and HIV model (Hill et al., 2018). We also test on synthetic tumor growth (Tumor) and pharmacokinetic system (PK) as well as on real Tacrolimus and Bicycle sharing datasets. For synthetic datasets, we consider two different noise settings: low (a = 0.01) and high (a = 0.1). Details about all datasets can be found in Appendix D.1.
making DSM applicable in settings with multi-dimensional inputs, some of the original issues remain. That includes restriction to finite compositions which prevents modeling systems with oscillatory or periodic behavior. Fitting the composition map remains the most cone-connuming part of the training due to the necessity of matching every composi- tion to each sample. Future research should focus on finding more efficient methods of fitting the composition map.
Results The results can be seen in Table 2. The error met- ric used is the same as in Section 5.6, i.e., Equation (5). We show sample trajectories predicted by EPISODE on those datasets in Figure 11 in Appendix A.3. We see that EPISODE or EPISODE outperforms other approaches most datasets. Performance on the SIR datasets shows that it can compete even on problems where ODE-based meth- ods have an advantage (a system of ODEs describing the dynamics the search space of SINDy). is promising that EPISODE can outperform even black-box approaches in certain settings. This can be explained by some of the induc- tive biases playing the role of regularizers (e.g., enforcing smoothness or small number of inflection points).
A Acknovledgements
This work was supported by Roche and AstraZeneca. We would like to thank Harry Amad, Benjamin Lapostolle, Vic- tor Baillet, and anonymous reviewers for their useful com- ments and feedback on earlier versions of this work.
Impact Statement
This paper presents a new method to improve understanding of a dynamical system's 's behavior through direct semantic modeling in personalized dynamical systems. Enhancing transparency in machine learning models is essential for tasks such as debugging, meeting domain-specific regula- tions, and mitigating potential biases. However, if not used correctly, these techniques may create a misleading sense of security regarding model outputs or be used merely for nom- inal compliance with regulations. Given that our method is relevant to critical fields like medicine and pharmacology, it is crucial to perform a comprehensive evaluation before implementing the model in these sensitive areas. Such as- sessments should confirm that the model's 's performance aligns with ethical standards and does not result in decisions that could harm individuals' health and well-being.
7. Discussion
Applications EPISODE revolutionizes the modeling of personalized dynamical systems by making them transpar- ent, editable, and intuitive. Unlike traditional ODE discov- ery. EPISODE skips complex analysis, directly predicting system behavior for rapid validation and real-time adjust- ments. This empowers experts to refine models dynamically, ensuring accuracy, reliability, and trust-crucial for high- stakes applications. By bridging machine learning with human intuition, EPISODE makes scientific modeling not just powerful, but truly usable.
Limitations and Open Challenges Although EPISODE resolves one of the biggest shortcomings of Semantic ODEs,
9

---

Skip the Equations: Learning Behavior of Personalized Dynamical Systems Directly From Data
References
Chen, R. T. Q., Rubanova, Y., Bettencourt, J., and Duve- naud, D. K. Neural Ordinary Differential Equations. In Advances in Neural Information Processing Systems, vol- ume 31. Curran Associates, Inc., 2018.
Akiba, T., Sano, S., Yanase, T., Ohta, T., and Koyama, M. Optuna: A next- eext-generation hyperparameter optimization framework. In Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, 2019.
Chen, Y., Angulo, M. T., and Liu, Y.Y.. Revealing Complex Ecological Dynamics via Symbolic Regression. BioEs- says, 41(12):190006) 2019. ISSN 1521-1878. doi: 00.1002/bies.201100066
Bertsimas, D. and Gurnee, W. Learning sparse nonlinear dynamies via maxed-integed optimization. Nonlinear Dy- namics, 11 111(7):6585-660 6604, 2023.
Cranmer, M. PySR: Fast & parallelized symbolic regression in Python/Julia. Zenodo, September 2020.
Biggio, L.,, Bendinelli T., Neitz, A., Lucchi, A., and Paras- candolo, G. Neural Symbolic Regression that Scales. In 38th International Conference on Machine Learning, July 2021.
D' Ascoli, S., Kamienny, P.-A., Lample, G., and Charton, F. Deep symbolic regression for recurrence prediction. Proceedings of the 39th International Conference on Machine Learning, pp. 4520-4536. PMLR, June 2022.
Blanchard, P., Devaney, R. L., and Hall, G. R. Differential Equations. Cengage Learning, July 2012. ISBN 978-1- 133-38808-1.
de Boor, C. On uniform approximation by splines. Journal of Approximation Theory, 1(2):219-235, September 1968. ISSN 00219045. doi: 10. N0.1006/0021-99444 990000221--9456689900----9
Bongard, J. and Lipson, H. Automated reverse engi- neering of nonlinear dynamical systems. Proceedings of the National Academy of Sciences, 104(24):9943- 9948, June 2007. ISSN 0027-8424, 1091-6490. doi: 10. (0073//mas.00994760000
de Silva, B., Champion, K., Quade, M., Loiseau, J.-C., Kutz, J., and Brunton, S. PySINDy: A Python package for the sparse identification of nonlinear dynamical systems from data. Journal of Open Source Software, 5(49):2104, 2020. doi: 10./os0/ 021044.
Breiman, L. Classification and Regression Trees. Routledge, 2017.
Dupont, E., Doucet, A., and Teh, Y. W. Augmented Neural ODEs. In Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
Brunton, S. L. and Kutz, N. Data-Driven Science and Engineering: Machine Learning, Dynamical Systems, and Control. Cambridge University Press, May 2022. ISBN 978- 178-1-00-911563-6
Fanaee-T, T, H. Bike sharing. UCI Machine Learning Reposi- tory, 2013.
Brunton, S. L., Proctor, J. L., and Kutz, J. N. Discov- ering governing equations from data by sparse identi- fication of nonlinear dynamical systems. Proceedings of the National Academy of Sciences, 113(15):3932- 3937, April 2016a. ISSN 0027-8424, 1091-6490. doi: 10.1073/pnas. 1517384411.
Han, Y. R., Lee, P. I., and Pang, K. S. Finding Tmax and Cmax in Multicompartmental Models. Drug Metabolism and Disposition,
November 2018. ISSN 0090-9556, 1521-009X. doi: 10.1124/dmd.1 18. 082636.
Hastie, T. and Tibshirani, R. Generalized additive models. Statistical Science, 1(3):297-318, 1986.
Brunton, S. L., Proctor, J. L., and Kutz, J. N. Sparse Identi- fication of Nonlinear Dynamics with Control (SINDYc). IFAC-PapersOnLLine, 49(18):71--715, January 2016b. ISSN 2405-8963, doi: 10.
Hill, A. L., Rosenbloom, D. I. S., Nowak, M. A., and Sili- ciano, R. F. Insight into treatment of HIV infection from viral dynamics models. Immunological Reviews, 285(1): 9-25, 2018. ISSN 1600-065X. doi: 10.1111/imr 12698.
Caruana, R., Lou, Y., Gehrke, J., Koch, P., Sturm, M., and Elhadad, N. Intelligible Models for HealthCare: Pre- dicting Pneumonia Risk and Hospital 30-day Readmis- sion. In Proceedings of the 21th ACM SIGKDD In- ternational Conference Knowledge Discovery and Data Mining, pp. 1721-1730, Sydney NSW Australia, August 2015. ACM. ISBN 978-1-4503-3664 2-64-.. doi: 10. (45/2783258.2788611
Holt, S., Qian, Z., and van der Schaar, M. Deep Genera- tive Symbolic Regression. The Eleventh International Conference on Learning Representations, 2023.
Holt, S. L., Qian, Z., and van der Schaar, M. Neural Laplace: Learning diverse classes of differential equations in the Laplace domain. In International Conference on Machine Learning, pp. 8811 8832. PMLR, 2022.
Chen, R. T. Q. Torchdiffeq. 2018.
10

---

Skip the Equations: Learning Behavior of Personalized Dynamical Systems Directly From Data
Kacprzyk, K. and der Schaar, M. Shape Arithmetic Expressions: Advancing Scientific Discovery Beyond Closed-form Equations. In Proceedings of The 27th Inter- national Conference on Artificial Intelligence and Statis- tics. PMLR, 2024.
Lee, C., Yoon, J., and van der Schaar, M. Dynamic- DeepHit: A Deep Learning Approach for Dynamic Sur- vival Analysis With Competing Risks Based on Longitu- dinal Data. IEEE Transactions on Biomedical Engineer- ing, 67(1): 67(1):122-133, January 2020. ISSN 1558-2531. doi: 10.1 1..109/TBME.2211.92990
Kacprzyk, K. and van der Schaar, M. No Equations Needed: Learning System Dynamics Without Relying on Closed- Form ODEs. In The Thirteenth International Conference on Learning Representations. 2025.
Liu, D. C. and Nocedal, J. On the limited memory BFGS method for large scale optimization. Mathematical pro- granming, 45(1):503 -28,, 1989.
Lou, Y., Caruana, R., and Gehrke, J. Intelligible models for classification and regression. In Proceedings of the 18th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD 12. pp. 150-158, New York, NY, USA, August 2012. Association for Computing doi: 10.1145/
Kacprzyk, K., Qian, Z., and van der Schaar, M. D-CIPHER: Discovery of Closed-form Partial Differential Equations. In Advances in Neural Information Processing Systems, volume 36, pp. 27609-276670 2764,, December 2023.
Kacprzyk, K., Holt, S., Berrevoets, J., Qian, Z., and van der Schaar, M. ODE Discovery for Longitudinal Het- erogeneous Treatment Effects Inference. In The Twelfth International Conference on Learning Representations, 2024a.
Machinery. ISBN 978-1-4503-1462-6 2339530.2339556.
Martius, G. S. and Lampert, C. Extrapolation and learning equations. In 5th International Conference on Learning Representations, ICLR 2017-Workstie Track Proceed- ings, 2017.
Kacprzyk, K., Liu, T., and van der Schaar, M. Towards Transparent Time Series Forecasting. In The Twelfth International Conference Learning Representations, 2024b.
Messenger. D. A. and Bortz, D. M. Weak SINDy: Galerkin- Based Data-Driven Model Selection. Multi- scale Modeling & Simulation, 19(3):1474-144) January 2021a. ISSN 1540-3459, 1540-3467. 20M1343166.
doi: 10.1137/
Kaheman, K., Kutz, J. N., and Brunton, S. L. SINDy-PI: A robust algorithm for parallel implicit sparse identification of nonlinear dynamics. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 47602242):20200279, October 2020. doi: 10.1098/rspa. 2020.0279.
Messenger, D. A. and Bortz, D. M. Weak SINDy for partial differential equations. Journal of Computational Physics, 443:1 10525, October 2021b. ISSN 00219991. doi: 10. (016/...p.22211.1053
Kaptanoglu, A. A., de Silva, B. M., Fasel, U., Kaheman, K., Goldschmidt, A. J., Callaham, J., Delahunt, C. B., Nicolaou, Z. G., Champion, K., Loiseau, J.-C., Kutz, J. N., and Brunton, S. L. PySINDy: A comprehensive Python package for robust sparse system identification. Journal of Open Source Software, 7(69)33994, 2022. doi: 10.21105/joss.0330 03994..
Mould, D. and Upton, R. Basic Concepts in Population Modeling, Simulation, and Model-Based Drug Develop- ment. CPT: Pharmacometrics & Systems Pharmacology, 1(9):6, 2012. ISSN 2163-8306. doi: 10.1038/psp.20122
Müller, M. Dynamic time warping. Information retrieval for music and motion, pp. 69-84, 2007.
Petersen, B. K., Larma, M. L., Mundhenk, T. N., Santiago, C. P., Kim, S. K., and Kim, J. T. Deep Symbolic Regres- sion: Recovering Mathematical Expressions From Data via Risk-seeking Policy Gradients. In ICLR 2021, 2021.
Kermack, W. O., McKendrick, A. G., and Walker, G. T. A contribution to the mathematical theory of epidemics. Proceedings of the Royal Society of London. Series A, Containing Papers of Mathematical and Physi- cal Character, 115(772):700-71, January 1997. doi: 00.1098/ssaa.192..011
Qian, Z.,, Kacprzyk, K., and van der Schaar, M. D-CODE: Discovering Closed form ODEs from Observed Trajec- tories. The Tenth International Conference on Learning Representations, 2022.
Kingma, D. P. and Ba, J. Adam: A Method for Stochastic Optimization, January 2017.
Quade, M., Abel, M., Shafi, K., Niven, R. K., and Noack, B. R. Prediction of dynamical systems by symbolic re- gression. Physical Review E, 94(1):012214, July 2016. ISSN 2470-0045, 2470-0053, doi: 10.1I../PhysRevE.. 94.012214.
Koza, J. Genetic Programming: On the Programmning of Computers by Means of Natural Selection. Complex Adaptive Systems. MIT Press, Cambridge, Mass, 1992. ISBN 978-0-262-11170-5

---

Skip the Equations: Learning Behavior of Personalized Dynamical Systems Directly From Data
Raissi, M. and Karniadakis, G. E. Hidden physies models: Machine learning of nonlinear partial differential equa- tions. Journal of Computational Physics, 357:125-144 March 2018. ISSN 0021-9991 doi: 10. jo.. 11.039.
Woillard, J.-B., de Winter, B. C. M., Kamar, N., Mar- quet, P., Rostaing, I.., and Rousseau, A. Population pharmacokinetic model and Bayesian estimator for two tacrolimus formulations-twerr daily Prograf and once daily Advagraf. British Journal of Clinical Pharmacol- ogy. 71(3):391-42,, March 201 I. ISSN 1365-2125. doi: 10.1 1365-2125. 201................
Reinbold, P. A. K., Gurevich, D. R., and Grigoriev, R. O. Using noisy incomplete data to discover models of spatiotemporal dynamics. Physical Review E, (1) 010203, January 2020. ISSN 2470-0045, 2470-0053. doi: 10.1105/PhysRevE.10.. 101.000000.
Rubanova, Y., Chen, R. T. Q., and Duvenaud, D. K. Latent Ordinary Differential Equations for Iregularly-Sampler Time Series. In Advances in Neural Information Process- ing Systems, volume 32. Curran Associates, Inc., 2019.
Rudy, S. H., Brunton, S. L., Proctor, J. L., and Kutz, J. N. Data-driven discovery of partial differential equations. Science Advances, 364):166 April 2017. ISSN 2375-2548, doi: 211116//ciad..1622666
Sahoo, S., Lampert, C., and Martius, G. Learning Equations for Extrapolation and Control. In Proceedings of the 35th International Conference on Machine Learning, pp. 4442-4450. PMLR, July 2018.
Schmidt, M. and Lipson, H. Distilling Free-Form Natural Laws from Experimental Data. Science, 324((923))81 1- 85, April 2009. ISSN 0036-8075, 1095-9203. doi: 10. 1126/science. 165893.
Stephens, T. Gplearn: Genetic programming in python, with scikit-learn inspired and compatible api, 2022.
Udrescu, S.-M. and Tegmark, M. AI Feynman: A physics- inspired method for symbolic regression. Science Ad- vances, 6( 6(16))ca April 2020. ISSN 2375-2548.
doi: 10.1 10.11166/caaa..ay226...1
Udrescu, S.-M., Tan, A., Feng, J., Neto, O., Wu, T., and Tegmark, M. Al Feynman 2.0: Pareto-optimal symbolic regression exploiting graph modularity. 34th Confer- ence Neural Information Processing Systems (NeurIPS 2020), 2021.
Wilkerson, J., Abdallah, K., Hugh-Jones, C., Curt, G., Rothenberg. M., Simantov, R., Murphy, M., Morrell, J., Beetsch, J., Sargent, D. Scher, H. Lebowitz, P. Simon, R., Stein, W. D., Bates, S. E., and Fojo, T. Estimation of tumour regression and growth rates dur- ing treatment in patients with advanced prostate can- cer: A retrospective analysis. The Lancet Oncology, 18(1):143-15) January 2017. ISSN 1470-2045, doi: 10. 477-2004511630633.4
12

---

Skip the Equations: Learning Behavior of Personalized Dynamical Systems Directly From Data
Table of supplementary materials
1. Appendix A: additional empirical results
2. Appendix B: theoretical results
3. Appendix C: implementation details
4. Appendix D: experimental details
5. Appendix E: extended related works
6. Appendix F; additional discussion
A. Additional Results
A.1. Horizontal Asymptote
Figure 6 shows GAM describing the horizontal asymptote of the composition (8++6, - 8-+h) of the original EPISODE trained on the Tacrolimus dataset. Figure 7 shows how this property function changes as we add a penalty term encouraging the asymptote to be close to O. Finally, Figure 8 shows the property function after it is set to to make the model biologically plausible.
Shize fancionn
Shape Lancton
Shape Lancion
Shape fanctin
-9,00025
0.025
-3,00350
-0.1
Reasency Shape
0.00075
-02
-0.025
-03
-0,05
Requency
10
10-
2.5
5.0
7.5
10.0
0.75
100
DOSE
DV_0
SEX
WEIGHT
HT
0.02
0.235
Share tunddod
0.01
0.230
12.00005
8.00
0.225
Shape
00.0
-0.01
-0,00010
0.02
Freesency
Frequency
Fresend
25-
25
7.5
10.0
12.5
15.0
200
4:20
420
802
CREAT
CTP
TORM
Figure 6. The original GAM describing the horizontal asymptote of the composition (s++8, 84-bb by 8 trained on the Tacrolimus datasct. make biologically plausible, should have an asymptote equal We describe how fix it Section 5.5
A.2. Maximum of the Trajectory
Figure shows the property function describing the maximum of the trajectory for the model trained the Tacrolimus dataset. Figure 10 describes the time which it is achieved. Both GAMs have been pruned remove shape functions that contribute less than le - 5.
A.3. Examples Trajectories
Figure 11 shows sample trajectories predicted on three datasets used in Section 6.
B. Theory
To realize direct semantic modeling, we need the following three ingredients.
Definition of the semantic representation of M-dimensional trajectory & E C²
13

---

Skip the Equations: Learning Behavior of Personalized Dynamical Systems Directly From Data
I=-12
0.02
0,000
2.5
Shape fundion
Shope furdiee
0.00
00.005
2.0
0.0
Sh
-0.020
-0.1
00.15
Froouino
Froouing
Frouuccy
10
10
20
SEX
1e-19
0.10
-0.00 00.00000
Shape lunciann
00075
Shape furetonn
0.05
00100
00000
00325
0.05
0.00185
00000
-0.0
-0.00 -0.00200
Frcquenes
Ficeuena
7.5
10.0
125
15.0
200
400
600
800
4B
CREAT
OVP
FORM
Figure 7. GAM describing the horizontal asymptote of the composition ($4-4b; 84- 8--b, 8-44) trained on the Tacrolimus dataset after adding penalty encouraging the horizontal asymptote be close 0. We now safely 0 without worrying much about lost performance.
2. Trajectory predictor Firjj that output a trajectory given semantic representation of this trajectory.
3. Semantic predictor Fsom that outputs semantic representation of æ given vector of static features V.
More formally, the semantic representation of an M -ddimensional trajectory æ is defined as follows.
Definition B.1 (Semantic representation of æ). Let æ be an A -ddimensional trajectory. The semantic representation of pair (Cz, P=), where Ca and pz vectors of entries defined (c=) Cz.. and (p=)... Prm
We can similarly define Funj-
Definition B.2 (Trajectory predictor). The th coordinate of the trajectory predicted by the trajectory predictor Fina defined as [m)] (c. Finj(CC::Pn. Pmm)
Of course, such trajectory predictor is consistent, i.e., the semantic representation of the trajectory predicted by the trajectory predictor is equal to the original semantic representation.
Proposition B.3 (Consistency of trajectory predictor). Let (c,p) be semantic representation of some trajectory. Let (c (c. p). Then (Cz, Pa) (c.
Proof. From Definition B.2, Im = F((m)
uaj (c. p) (Cm, Pm). From Fm.. = Furaj (Cms Pm) and by the original definition of Finj, we get that (Crm, PPr..) = )C From Definition B.I, Cr... -Czz))n and prmm. (p )m. Thus (Cx) = Can and (px) m Pm for each m. Therefore Pc)) (c.
Positive GAMs
Proposition B.4. Let h(() (v())) + be GAM. Let assume that h is positive and continuous for all possible imputs from a hyperbox impul domain V = (01, bi) x x... x [ax, bx). Then can be represented using positive shape functions and a positive bias term.
Proof. Let 7 = minrelaa Ih (v) for each &s E [K] If h is positive and continuous for all possible inputs from a hyperbox then there exists such that h(v) > for all v E V. Let us define new shape functions
IN (vx) gN (vx) 7k </(K +
Then the new bias term needs to be equal to
K Σ
Tu 9u
(ūx(vx) (ux))
(8)
Am
14

---

Skip the Equations: Learning Behavior of Personalized Dynamical Systems Directly From Data
0.050
0.050
0.490
0.450
0.030
Shape lundion.
Shape fundlee
Shope Tundlen
Shop: funclan
Shapp fandlee
3.000
0.000
0.203
0.000
0.000
13 030
10,00
-0.00
-0.50
-0050
2.5
5.0
7.5
10.0
0.25
0.50
0.75
80
100
30
40
DOSE
DV.
SEX
WeiGH
HT
0.050
0.050
0.050
0.050
0.050
Shaps Tunddon
Shaps Tunddin
Shaps funduun
Shaps fundlun
3.000
0.000
0.000
0.000
0.000
-3.025
-0.025
-0.025
-0.023
030
10.000
9..50
-9.50
-0.550
Foauency
Feqqueayy
800
118
CREAI
1DRM
Figure 8. GAM describing the horizontal asymptote of the composition (8++0: 84 +h) trained the Tacrolimus dataset after setting it to O.
Thus 9u-
90
K E (-76 + +e/(K + 1)). This means
K
K
Do 90
Th
(9)
K+1
k=
As V is hyperbox then it must hold that Σ 7/ + 90 > c.
That means that
K = e/(K +
go xx
(1-K K+ K+1
+ 00
(10)
Thus all shape functions and the bias term are positive.
C. Implementation
C.1. Block diagrams
Figure 12 shows a block diagram of how a single semantic predictor works. Figure 13 shows the overall training procedure of single semantic predictor.
C.2. Code
The code of the implementation and the is EPISODE and https::////...c htt///gttttt con./vanderrschaaraaa//EE...
C.3. Composition map
To find a decision tree that minimizes such an objective, we perform greedy search where, at each node, we evaluate all possible splits for all possible features and choose the one that gives the smallest error. Consider a leaf of a partially fitted tree and subset of [D] denoted I such that all (dd) for d E I satisfy all conditions from root to this leaf. We define the current error at the leaf to be EII) = min.ec [2 Exere [c] [c]. The goal is to find feature split that divides into Z, and Lg
In such that EZ that is unnecessarily more complex without a significant improvement in performance. For instance, we set maximum depth of the tree or a minimum relative improvement for split to occur. We also penalize longer compositions. We set minimum number of samples a leaf so that we have enough data to train a property sub-map for the composition predicted by that leaf. Finally, we train trees for different maximum depths and choose the one with the lowest validation loss to prevent overfitting.
[Z [Z E(Z)+ + [2T E(I2) < (Z)). In practice, we put a few additional constraints to prevent a composition map
15

---

Skip the Equations: Learning Behavior of Personalized Dynamical Systems Directly From Data
1e
le
2.0
Shape tunction
Shape function
0.50
function
1.5
0.05
1.0
0.00
0.00
0.5
0.05
-0.25
0.0
Frequncyy
Prequency
Frequency.
10
10
25
2.5
7.5
10.0
0.25
0.50
0.75
7.5
12.5
15.0
DV_0
1.050
0.05
Shape function
0.05
1.025
0.00
0.00
Blas
1.000
0.975
-0.10
Frequency
Frequency
25
25
1
FORM
Figure 9. GAM describing the maximum value the composition ($4-4b, 84-b, 8_ 8-bb, 8-h) trained on the Tacrolimus dataset. This property function has been pruned to remove shape functions that contribute less than 1e
C.4. Property map
Property functions need satisfy many constraints to make that the properties predicted by them actually describe possible trajectory. For instance, property functions describing æ=(o) and sx(1) need to satisfy ==(o) < xx(1) )) for all possible if the first motif is increasing. To achieve that, we design the property maps so that they satisfy most of these constraints by construction. For instance, instead of fitting the property function for x(4) directly, we fit a GAM for the difference between zto) and z(t1). As long as we make this difference positive (or negative the first motif is decreasing), we can add this GAM the one describing x((a) and get a valid property function. How can we get a GAM that is always positive? We observe that if GAM positive for all possible inputs (from a hyperbox input domain), then it is always possible represent it using only positive shape functions and positive bias (proof in Appendix B). Thus, can limit ourselves to only such shape functions. Ultimately the property function for xtf1) evaluated at is given by
K
( (ww)) + soft, (8i) xx(o)
soft
(1I)
This not only guarantees that xx41) > xlo) but, crucially, as x((0) is also described by (unconstrained) GAM, x(41) is still a GAM where the th shape function is given by
soft, (of + (wx)
(12)
and the bias term is equal to (S1) + Su-
Works Well With Categorical Features Because GAMs work very naturally with categorical features (gk. just assigns different value to each category), this makes our model much more appropriate than ODEs settings where categorical features are present. ODEs, by default, can only work with numerical variables, and one-ho hot encoding categorical features may result in complicated equations.
16

---

Skip the Equations: Learning Behavior of Personalized Dynamical Systems Directly From Data
0,010
Staped tancon
0015
Shapp cunchan
0.002
Shape funciin
0.000
0.005
0.01
Shopp
-0.001
-0010
Reavency
Requency
25
25
25-
25
5.0
7.5
10.0
0.30
0.75
DOSE
DV_0
SEX
FORM
0.0000
0,0875
Figure GAM describing at which the maximum is obtained of the composition (#+46, 84-b, 8--h, 8. 4h) trained on the Tacrolimus dataset. This property function has been pruned to remove shape functions that contribute less than 1e 5
D. Experimental Details
D.1. Datasets
SIR The SIR dataset is a synthetic dataset based on an SIR epidemiological model (Kermack et al., 1997). a
three-dimensional dynamical system governed by the following system of ODEs,
dS
3ST
(13)
dt
d/
381 7 dR
(14)
dt
77,
(15)
dl
and initial conditions
S(0) = So
(16)
I(0) = Iu
(17)
R(0) Ro-
(18)
We set 3 0.3 and 7 = 0.1. We sample So uniformly from (0.6, 1.0), 1o uniformly from (0.01, 0.1, Ro uniformly from (0.0, 1.0). We generate 500 trajectories, each measured 20 equally spaced time points between and 100. Finally, we add Gaussian noise to each measurement with 0.01 and divide the time points by 100 scale them to interval (0, 1). The set of static features contains only the initial conditions, i.e., (So, Iu. Ro).
PK The PK dataset is synthetic dataset based on a pharmacokinetic model published by Woillard et al. (2011). It predicts the drug concentration over time based on patient covariates and the dose of the drug. It has form of a system of ODEs and covariate models describing the parameters of the ODEs in terms patients' covariates. It consists of the following system
17

---

Skip the Equations: Learning Behavior of Personalized Dynamical Systems Directly From Data
SIR dataset
16
0.
(sSeereci
64-
tomevel
Kppnckod
077
PK dataset
1.7-
Can Corcuttztinn gpodeeds
14
Canderrad (pncieed)
32-
1.00
13-
300
0.56
0.
0.0
0.00
0.3
02
3.4
3.4
3.3
13
0.3
3.2
0.1
3.5
38
HIV dataset
2.00
Tralenent
Tepsdited
1:2-
3.5
1.00
0.30
0.0
Figure 11. Sample trajectories predicted by EPISODE* top row: SIR dataset (a = 0.01), middle row: PK dataset (a 0.01), bottom row: HIV dataset = 0.01).
of ODEs.
dCdepot
hig Catepoo
(19)
dt
dCiranst
kig Cdepor
Ain Guanstl
(20)
dt
dCirans2
hin Cirunsl
Crranss
(21)
dCurans3
kig Cirunzz
Niramss
(22)
d.ccatt
ki Cranss
((CT+ Ocent V1) (Q Openi V2)
(23)
dt
dCperi
(Q Ceem VV) (Q Cperi/ / V22)
(24)
18

---

Skip the Equations: Learning Behavior of Personalized Dynamical Systems Directly From Data
(=m) senr
Fcom
(
gg
Pprop
+( prop
( prop kun
m
m
JAK
WMR
0.0,04,,.88114
0.3,0.7,0.0
2.8,2.9
Composition:
Transition points:
(18.4).(11.114))
Derivatives:
x(t)) 0.3, (G,) 07, =
Properties of Sa-a!
2, 2.9
Semantic representation of Xm
Figure Block diagram depicting single semantic predictor
The initial conditions are:
depoot (0) UDOSE
(25)
Otranst (0)
(26)
Ciramsz (0) 0
(27)
Ctransa (0) 0
(28)
V/1
Coent(O) (0) = UDVO X
(29)
1000
Operi (0) 1/2
x Czent (0).
(30)
Vi
The parameters are described as follows.
UHT
IICCL
CL TVCL X
CYPCL PCyP
(31)
35
Vi = TVVI x STV1 FTORM
(32)
Q = TVQ
(33)
1/2 TVV2
(34)
bi TVKTR STKTR STORM
(35)
19

---

Skip the Equations: Learning Behavior of Personalized Dynamical Systems Directly From Data
For each (p((), ((0), y(at) D:
decision Foom
D
For each E C': [cd min perc Na
that minimizes
Fcom
1
|l (()m)().))(rdd)) -
[a] 2 Pram (v(()
(d)
Fsem
Divide D according to Feom
De
Fit property submap
Fprop
For each e; E Im(Feoom::
Pprop
De (1o(t),(0)) ybu)) ((0))
(D(1) ci
Dez
Fit property submap
prop
argmin w.B
Σ Na
|iraj (From (Dl()) Firrpp ( ww.)) y(u))
Figure 13. Block diagram depicting the general training procedure
The constants are defined as follows.
TVCL 21.2
(36)
TVV1 486
(37)
TVQ = 79
(38)
TVV2 271
(39)
TVKTR 3.34
(40)
HTCL = 1.14
(41)
CYPCL 2.00
(42)
STKTR 1.53
(43)
STVI 0.29
(44)
(45)
We generate 500 patients by uniformly sampling the static features from the ranges described in Table 3. Then we solve the initial value problem to get all trajectories but return (observe) only Ceent We then convert to proper units by multiplying by 1000 Each trajectory is measured 20 equally spaced time points between and 24. Then we divide the time points by 1.
24 and the trajectories by 20 get them in range (0, 1). We also divide upvp by 20 corresponds to the initial value. Finally, we add Gaussian noise with 0.01.
Tumor The Tumor dataset is synthetic dataset used by Kacprzyk et al. (2024b). It is based on the tumor growth model proposed by Wilkerson et al. (2017). is described by the following equation.
x(t) o exp(-dt) (1 +) exp(gt)
(46)
Where Φ, d,g g are defined as follows.
9 go (Vage/ / /2..0)0.5
d do * Udosage/ Wweight
(47)
Φ = 1/(1 + exp( Udosage Φu))
20

---

Skip the Equations: Learning Behavior of Personalized Dynamical Systems Directly From Data
Table Feature ranges in PK dataset
Feature
Range
"DOSE
(1,10)
VDVO
(0,20)
USEX
0,1
OWEIGITT
(45,110)
UHT
(20,47)
VHN
(6,16)
UCREAT
(60,830)
UCYP
0,1
UFORM
0,1
And where 9u, do, do set to
9o 2.0
dg 180
(48)
6 = 10.
We generate 500 samples by uniformly sampling the static features from the feature ranges described in Table 4. Each trajectory is measured at 20 equally spaced points on (0, 1). Finally, we add Gaussian noise with = 0.01.
Table 4, Feature ranges in Tumor dataset
Feature
Range
Vage
(20,80)
Uweight
(40,100)
Uintitial
(0.1,0.5)
Gdose
(0.0,1.0)
Tacrolimus We take the drug concentration curves of two tacrolimus formulations obtained by (Woillard et al., 2011). The static features include sex, weight, hematocrit, hemoglobin, creatinine, dose, CYP3A5 genotype, and formulation. We also supplement it with the initial concentration of the drug in the blood. After preprocessing, we end up with 90 samples, each with 10 measurements on average. As not all trajectories are measured at exactly the same time points, we perform linear interpolation for each trajectory and sample it at the following time points: 0, 0.33, 0.67, 1, 1.5, 2, 4, 6, 9, 12. Finally, we divide the time points by 12 and the trajectories by 20 bring both them into (0, 1) range. We also divide the initial condition by 20.
Bike sharing We use real bike-sharing dataset from the UCI repository (Fanace-T, 2013). We perform preprocessing so that each sample corresponds to one day measured 24 times (every hour). As features Uweathersits Utemp= Vhum and Uwindspeed change throughout the day, we take their average make them into static features. We also round Wweathersi to keep it categorical variable. We also add the value at = 0 as additional feature. At the end we have the following features: Uno- We divide the trajectories and Ux, by 500 to scale it and the time
Useason Umonth, Uworkingday, tweathersit. Ttemp Chum, Uwindspood-
points by 23 to bring them to (0, 1) range.
D.2. Baselines
SINDy We use SINDy (Brunton et al., 2016b) as implemented in PySINDy package (de Silva et al., 2020; Kaptanoglu et al., 2022). We pass the variable t and the static features v as control variables. We use a library with all polynomial terms up to the second degree, including interaction between different variables. We use Mixed- Integer Optimized Sparse Regression (MIOSR) (Bertsimas & Gurnee, 2023) for optimization allows us to choose sparsity level- the number of
21

---

Skip the Equations: Learning Behavior of Personalized Dynamical Systems Directly From Data
terms in the equation. In our experiments, we consider two variants of SINDy, that we denote SINDY-5 and SINDy-20, enforcing the maximum number of terms be, respectively, and 20. We tune the parameters of MIOSR that describes the strength of L2 penalty (between le 3 and 1), and the derivative estimation algorithm using Optuna (Akiba et al., 2019) for 30 trials. We choose between the following techniques: finite difference, spline, trend filtered, and smoothed finite difference as available in PySINDy. The parameter ranges we consider for each of them are shown in Table 5.
Table 5. Hyperparamcter ranges for each of the derivative estimation methods.
Method
Hyperparameter ranges
finite difference
k E (1.....5) 5)
spline
8 S (1e-3,) 3,-3,,1. 3, order E {0, 1, 2), a E (le-4,1) 4,1) window length E {l....,5}
trend filtered
smoothed finite difference
WSINDy We use WSINDy (Reinbold et al., 2020; Messenger & Bortz, 202 based on (Reinbold et al., 2020) as implemented in PySINDy package (de Silva et al., 2020; Kaptanoglu al., 2022). We pass the variable t and the static features U as control variables. We use library with all polynomial terms up to the second degree, including interaction between different variables. We use Mixed- Integer Optimized Sparse Regression (MIOSR) (Bertsimas & Gurnee, 2023) for optimization as allows us to choose a sparsity level-the number of terms in the equation. In our experiments, we consider two variants of WSINDy, that we denote WSINDy-5 and WSINDy-20, enforcing the maximum number of terms to be, respectively, 5 and 20. We tune the parameter a of MIOSR that describes the strength of the L2 penalty (between 1e and 1) using Optuna (Akiba et al., 2019) for 30 trials. We choose the parameter K (the number of domain centers) to be 200.
NeuralODE We implement NeuralODE (Chen et al., 2018) model using torchdiffer library (Chen, 2018). We parametrize the ODE neural network that takes æ, t, and inputs and outputs an M-dimensional vector. We then append vector of zeros of size K to account for the static features being constant (having derivative 0). The network consists of fully connected layers. The categorical variables are one-hot-encoded before training and the data is standardized. We set the batch size to 32 and 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 are shown in Table 6.
Table 6. Hyperparameter ranges used for tuning Neural ODE.
Hyperparameter
Range
learning rate
(le-5,le-1)
number of layers
(1,3)
units in each layer (separately)
(16,128)
dropout rate
(0.0,0.5)
weight decay activation function
(1e-6,1e-2)
ELU, Sigmoid
ANODE We implemented ANODE (Dupont et al., 2019) similarly to NeuralODE with the difference that we augment the x by a few dimensions and evolve this extended dynamical system. We set the initial conditions for these additional dimensions to O. We train and tune in the same way as NeuralODE with the only difference being that we also tune the number of additional dimensions. We consider a range between I and 10.
LatentODE We implement a non- non- pon-brbbabilistic variant of LatentODE (Rubanova et al., 2019). It is similar to NeuralODE with the difference being that the ODE is solved in the latent space, so we also learn an encoder and decoder. We tune hyperparameters using Optuna (Akiba et al., 2019) for 20 trials. Ranges for the hyperparameters shown in Table 7.
D.3. EPISODE
The composition libraries used in EPISODE* depend on the particular dataset and are as follows.
22

---

Skip the Equations: Learning Behavior of Personalized Dynamical Systems Directly From Data
Table 7. Hyperparameter ranges used for tuning Latent ODE.
Hyperparameter
Range
learning rate
(Ie-5,1e-1)
number of layers in the encoder
(1,3)
units in each layer of the encoder (separately)
(16,128)
number of layers in decoder units in each layer of decoder (separately)
(1,3)
(16,128)
number of layers in the derivative
(1,3)
units each layer of derivative (separately)
(16,128)
latent dimension dropout rate weight decay activation function
(1,10)
(0.0,0.5)
(1c-6,1e-2)
ELU, Sigmoid
SIR. m 1: {(8- 8- (-))), 2: {(8+40: (4 m 3: {(++ 8+-h)-
PK. m 1: (844b. b- 8-th)
Tumor m = { ((+++0) (8-46) (s (s=+h);
Tacrolimus, (844b)
St
Bike sharing. 1: ((5_.,,
bb)s (8 (844b: 84 4bs 8-4)+
D.4. Computational Cost and Complexity
As mentioned in Section 7, training the composition map is the most time-consuming process, as it requires preprocessing step of fitting every admissible composition to each sample (this can be parallelized). The actual time to fit the decision tree is negligible in our experiments. This preprocessing step has time complexity O(DMIC') where D is the number of samples, M is the dimensionality of the trajectory, and [C'| is the number of admissible compositions. Crucially, this time does not depend on K (the dimensionality of static features V/).
Table 8 shows the time needed train the composition map for each dataset. We currently computationally intensive procedure of minimizing a Dynamic Time Warping (DTW) distance and then the standard MSE error. There is trade-off between the accuracy of fit and the computation time. Practitioners may opt for more approximate fits (e.g., by not using DTW).
Table Time needed to train the composition map for each of the datasets. All experiments were performed on 18-core Intel Core 19- 10980XE with 60GB of RAM.
Dataset
Noise
Time (min)
compositions
I samples
Time per sample (s)
Time per composition (s)
Time per single fit (s)
SIR
0.01
233.67
26
500
28.04
539.24
0.36
STR
0.1
178.66
500
21.44
412.28
0.27
PK
0.01
68.73
500
8.25
158.
0.32
PK
0.1
66.65
26
500
8.00
153,80
0.31
Tumor
0.01
77.33
26
500
9.28
178.45
0.36
Tumor
0.1
66.37
26
500
7.96
153.16
0.31
Tacrolimus
real
8.65
26
90
5.77
19,96
0.22
Bike sharing
real
462.27
42.35
815.77
1.25
HIV
0.01
259.33
26
500
31.12
598
0.40
HIV
0.1
194.33
26
500
23.32
448.45
0.30
The advantage of training the composition map separately from the property map is that the errors table e(4) [c] (Equation (1))
can be computed once, saved, and then reused for different training settings of the composition decision tree and the property maps. This makes introducing changes the model much quicker.
23

---

Skip the Equations: Learning Behavior of Personalized Dynamical Systems Directly From Data
Table 9 shows the overall computational cost for all methods all datasets (including five runs and hyperparameter tuning). Note that we calculate the errors table c( [c] (Equation (1)) once for all five seeds and then subsample the results based on the split.
Table 9. Computation time (in minutes) for all experiments (five runs with hyperparameter tuning). All experiments were performed on 18-core Intel Core 10-10800E with 60GB of RAM.
SIR
PK
Tumor
Tacrolimus
Bike sharing
HIV
Noise
0.01
a 0.1
0.01
0 = 0.1
0.01
= 0.1
real
real
0 = 0.01
or = 0.1
NeuralODE
86.07
215.03
94.09
521.42
49.90
56.70
9.57
146.17
265.32
152.75
ANODE
94.79
185.77
94,36
291.19
56.09
62.32
10.95
187.07
139.76
140.12
LateniODE
67.63
149.08
56.64
64.77
46.66
51.05
6.92
94.78
79.56
57.19
SINDY-5
14.87
26.11
40.91
44.43
8.96
23.39
14.03
25.22
81.38
91.08
SINDy-20
15.99
23.21
57.63
82.76
8.19
19.51
32.68
60.86
105.65
93.98
WSINDy-5
12.74
13.58
26.41
13,00
10.83
11.16
32,06
15.62
42.56
29,51
WSINDy-20
11.60
13.36
115.38
10.39
38.83
39.06
27.96
EPISODE (property mups)
100.20
75.64
118.07
44.96
31.62
55.84
35.12
104.84
166.23
169.67
EPISODE (composition maps)
233.67
178.66
68.73
66.65
77.33
66.37
8.65
462.27
259.33
194.33
EPISODE (whole)
333.87
254.30
186.80
111.61
108.95
122.21
43.77
567.11
425.55
364.00
EPISODE (property maps)
98.04
82.87
55.15
53.35
23.02
24.42
34.34
85.57
152.04
156.83
EPISODE* (composition mups)
8.99
6.87
2.64
2.56
11,90
10.21
0.33
27.19
26.60
19,93
EPISODE* (Whole)
107.03
89.74
57.79
55.91
34.92
34.64
34.67
112.76
178.63
176.76
E. Extended Related Works
Symbolic Regression and Discovery of Differential Equations The identification of differential equations is often viewed part of the larger field known as symbolic regression. This area of machine learning aims to express data through a closed-form mathematical expression. Traditionally, genetic programming has been employed for symbolic regression (Stephens, 2022; Cranmer, 2020), but recent advancements have seen the integration of neural networks. This involves methods such as directly representing equations within neural networks by modifying activation functions (Martius Lampert, 2017; Sahoo et al., 2018), leveraging neural networks to narrow down the search space (Udrescu & Tegmark, 2020; Udrescu et al., 2021), exploring equations through reinforcement learning (Petersen al., 2021), and employing large pre-trained transformers (Biggio et al., 2021; D' Ascoli al., 2022; Holt et al., 2023). Standard symbolic regression can be adapted for ordinary differential equation (ODE) discovery by estimating derivatives from data to treat them as targets (Quade et al., 2016). Despite this, several specialized ODE discovery techniques have emerged. Among these, the Sparse Identification of Nonlinear Dynamical Systems (SINDy) (Brunton et al.. 2016b) stands out, describing derivatives as linear combinations of functions from a predetermined library. This has led to various extensions, including those for implicit equations (Kaheman et al., 2020), controlled equations (Brunton et al., 201 and treatment effect estimation (Kacprzyk et al., 2024a). Alternative approaches utilize weak formulations of ODEs, circumventing the need for derivative estimation (Messenger& & Bortz, 2021a; Qian et al., 2022). Additionally, methods specifically designed for discovering partial differential equations (PDEs) have been developed, as standard ODE discovery methods cannot be directly applied to PDFs (Rudy et al., 2017; Raissi & Karniadakis, 2018; Messenger & Bortz, 2021 Kacprzyk et al., 2023). The challenge of finding comprchensible and well-fitting fosed-form expressions inspired Shape Arithmetic Expressions (Kacprzyk & van der Schaar, 2024) that extend the prespecified set of well-known functions (e.g.. trigonometric or exponential functions) in symbolic regression by flexible and learnable univariate functions. Those functions do not necessarily have a compact symbolic representation but are supposed to be understood by looking at their line plot.
Black Box Approaches for Modeling Dynamical Systems Black-box approaches to modeling dynamical systems have gained significant attention due to their ability learn complex dynamics without requiring explicit governing equations. Neural Ordinary Differential Equations (Neural ODEs) (Chen al., 2018) introduced a model that parameterizes differential equations using neural networks, enabling flexible trajectory modeling. Building upon this, Augmented Neural ODEs (Dupont et al., 2019) addresses the limitations of Neural ODEs in handling complex dynamics by expanding the state space, allowing for more expressive representations. Latent ODEs (Rubanova et al., 2019) further extend this framework by inferring latent representations of time-series data. Neural Laplace (Holt et al., 2022) generalizes Neural ODEs by learning solutions in the Laplace domain, offering improved performance for systems governed by diverse types of differential equations.
24

---

Skip the Equations: Learning Behavior of Personalized Dynamical Systems Directly From Data
F. Additional Discussion
F.I Sensitivity the Choice C'
The results in Table 2 (EPISODE) are for a generous set of compositions. We consider all compositions up to length (or length 8 for the Bike sharing dataset) except for compositions where two consecutive transition points are inflection points. This results in 26 and 34 compositions, respectively. Thus, EPISODE achieves competitive results for relatively large choice of compositions (and minimum prior knowledge). EPISODE* compositions, constraining the set to between 1 and compositions depending the dataset (details in Appendix D.3). We can see that it significantly improves the performance in some settings (Tumor dataset), but often has no significant impact (Bike sharing, Tacrolimus, HIV). In these settings, EPISODE is able to autonomously determine which compositions should be used in the composition map.
prior knowledge about the shape of the
F.2. Extending to Periodic or Oscillating Trajectories
EPISODE could, in the future, be extended to periodic or oscillating trajectories by extending the definition of semantic representation to accommodate periodic compositions. Let us assume that we know that the trajectories are described by an infinite composition consisting of repeating ($+-e: 8- 8- 84+c). Then we can describe this composition segment using properties like "amplitude" or "frequency". For each trajectory in our training dataset, we can extract additional trajectories describing how these properties change time. We
model these auxiliary trajectories using framework similar to EPISODE. However, we leave a full exposition of this idea for future paper.
25


Reviews:
[{"review_id": "v9lOCPegIK", "reviewer": "Reviewer_evCm", "summary": "The paper proposes a method called EPISODE for learning personalized dynamical systems (PDS) without explicitly discovering ordinary differential equations (ODEs). As mentioned in this paper, traditional approaches to modeling dynamical systems involve first identifying closed-form equations and then analyzing their properties. In contrast, EPISODE bypasses this two-step process by directly predicting a semantic representation of the system’s behavior (including the trajectory’s shape and key quantitative properties) from static inputs or covariates. The method utilizes decision trees to predict the shape of trajectories and generalized additive models (GAMs) to model their quantitative properties.\n\n## update after rebuttal\nI am quite satisfied with the rebuttal.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Dear Reviewer evCm,\n\nThank you very much for your positive review! We are glad you found our method novel and our claims well-evidenced. We reply to your questions first, and then we discuss other weaknesses and suggestions.\n\n### Q1. Runtime complexity of the composition map training\n\nAs mentioned in the limitations, training the composition map is the most time-consuming process as it requires a preprocessing step of fitting every admissible composition to each sample (this can be parallelized). The actual time to fit the decision tree is negligible in our experiments. This preprocessing step has time complexity $O(D M |\\mathcal{C}'|)$ where $D$ is the number of samples, $M$ is the dimensionality of the trajectory, and $|\\mathcal{C}'|$ is the number of admissible compositions. Crucially, this time does not depend on $K$ (the dimensionality of static features $\\mathcal{V}$).\n\nWe have now included a table showing the time needed to train the composition map for each dataset (see [here](https://imgur.com/a/90rECtR)). We currently use a computationally intensive procedure of minimizing a Dynamic Time Warping (DTW) distance and then the standard MSE error. There is a trade-off between the accuracy of fit and the computation time. Practitioners may opt for more approximate fits (e.g., by not using DTW).\n\nThe advantage of training the composition map separately from the property map is that the errors table $e^{(d)}[c]$ (Eq. 1) can be computed once, saved, and then reused for different training settings of the composition decision tree and the property maps. This makes introducing changes to the model much quicker. The training time of the property maps (for five seeds and with tuning) can be found [here](https://imgur.com/a/4jhWn1H). \n\n**Actions taken**: Added a discussion on the computational complexity in a newly created Appendix F (Additional Discussion) and tables with computation times in Appendix D.\n\n\n### Q2. Sensitivity to the initial choice of compositions\n\nThe results in Table 2 (EPISODE) are for a generous set of compositions. We consider all compositions up to length 4 (or length 8 for the Bike sharing dataset) except for compositions where two consecutive transition points are inflection points. This ends up being, respectively, 26 and 34 compositions. Thus, EPISODE achieves competitive results for a relatively large choice of compositions (and minimum prior knowledge). EPISODE* uses prior knowledge about the shape of the compositions, constraining the set to between 1 and 4 compositions depending on the dataset (details in Appendix D.3). We can see that it significantly improves the performance in some settings (Tumor dataset) but often has no significant impact (Bike sharing dataset or Tacrolimus). In these settings, EPISODE is able to autonomously determine which compositions should be used in the composition map. This finding generalizes to the additional experiments that we have run (see [here](https://imgur.com/a/5nnGzqi)).\n\n**Actions taken**: Added a discussion on the sensitivity of the method to the choice of the set of admissible compositions in a newly created Appendix F (Additional Discussion).\n\n### Q3. Automatic strategy for composition selection\n\nAs described above, the algorithm can usually autonomously determine which compositions should be used in the composition map. Thus, a viable approach is to start with a large number of compositions and let the algorithm narrow it down. In the future, we would like to design a fast algorithm that provides an approximate fit between a composition and a trajectory that narrows down the search space before the actual fitting.\n\n### Q4. Periodic or oscillatory behavior\n\nEPISODE could, in the future, be extended to periodic or oscillating trajectories by extending the definition of semantic representation to accommodate periodic compositions. Let us assume that we know that the trajectories are described by an infinite composition consisting of repeating $(s_{+-c},s_{--c},s_{-+c},s_{++c})$. Then we can describe this composition segment using properties like \"amplitude\" or \"frequency\". For each trajectory in our training dataset we can extract additional trajectories describing how these properties change over time. We can now model these auxiliary trajectories using a framework similar to EPISODE. However, we leave a full exposition of this idea for a future paper.\n\n**Actions taken**: Added a discussion on extending EPISODE to infinite compositions in a newly created Appendix F (Additional Discussion).\n\n### W1. Finite compositions\n\nPlease see our response to Q4.\n\n### W2. Complexity of composition map training\n\nPlease see our response to Q1.\n\n### Minor\n**Actions taken**: We have now added explicit guidelines on interpreting GAM plots in Appendix E.\n\n\n---\n\nWe hope we have addressed all your concerns. Please let us know if you have any additional questions, we are more than happy to address them.\n\nKind regards,\n\nAuthors", "claims_and_evidence": "The paper generally supports its claims with clear evidence. The key claim is that EPISODE offers a transparent, editable, and interpretable alternative to traditional equation discovery or black-box methods. This claim is supported by three aspects:\n- Demonstrations on synthetic and real datasets (e.g., SIR model, Tacrolimus).\n- Comparisons against established baselines.\n- The detailed case study shows the practical benefit of transparency, interpretability, and editability.\nHowever, some claims regarding scalability to more complex dynamical systems or handling oscillatory behaviors are acknowledged as limitations.", "methods_and_evaluation_criteria": "In my opinion, the methods and evaluation criteria are very reasonable and relevant. EPISODE is rigorously compared with both closed-form ODE discovery methods (SINDy variants) and black-box approaches (Neural ODEs variants), using metrics such as mean squared error (MSE). The authors also carefully incorporate practical interpretability considerations into their evaluation, demonstrating the method’s practical benefits through real-world medical datasets.", "theoretical_claims": "Yes. But the paper does not primarily rely on complex theoretical claims.", "experimental_designs_or_analyses": "Yes. The experimental designs and analyses are sound and clearly structured. Experiments effectively illustrate the method’s strengths in interpretability and performance. The authors also conduct clear ablation studies and robustness checks.", "relation_to_broader_scientific_literature": "The paper positions itself clearly within the broader literature on dynamical system modeling and equation discovery (e.g., references to SINDy, Neural ODEs, ANODE, Semantic ODEs).", "essential_references_not_discussed": "The paper adequately references and discusses existing works central to understanding its contribution. No immediate essential references seem missing.", "other_strengths_and_weaknesses": "Strengths:\n- Novel combination of direct semantic modeling and GAMs, resulting in transparent and interpretable modeling. Simple but effective.\n- Practical relevance demonstrated through pharmacokinetic applications.\n\n\nWeaknesses:\n- Method restricted to finite compositions, limiting application to oscillatory or periodic systems.\n- The complexity and runtime involved in the composition map training could limit scalability.", "other_comments_or_suggestions": "The writing and presentation are generally clear. Minor typos (e.g., spacing, punctuation) could be corrected upon careful proofreading. More explicit guidelines on interpreting GAM plots might further benefit readers unfamiliar with GAMs.", "questions_for_authors": "1. What is the runtime complexity of the composition map training as the dataset size and dimensionality grow? \n2. How sensitive is the method to the initial choice of compositions? \n3. Could an automatic or semi-automatic strategy for composition selection apply?\n4. How to possibly adapt or extend EPISODE to systems with periodic or oscillatory behavior?", "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": {"ethical_review_concerns": "N.A."}}, {"review_id": "hpQSMGvtMk", "reviewer": "Reviewer_g6tH", "summary": "The paper proposes EPISODE, a framework for learning the behavior of personalized dynamical systems without requiring explicit equation discovery. Instead of the traditional two-step approach of identifying ODEs and then analyzing them, EPISODE directly predicts the semantic representation from data. The paper extends prior direct semantic modeling approaches to accommodate multi-dimensional trajectories and integrate personalization via auxiliary static features.", "questions": "", "limitations": "", "rating": 2, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Dear Reviewer g6tH,\n\nThank you very much for your review. We appreciate your time and effort spent reviewing our paper. We are glad you found our conceptual framework solid, the case study well-motivated, and the comparisons comprehensive. We address your comments below.\n\n### Reliance on inductive biases\n\n> The method seems reasonable but highly relies on the incorporation of semantic inductive biases.\n\nThank you for this comment, as it allows us to clarify the experimental setup in our paper. Although our method *can* incorporate semantic inductive biases, it *does not require* this information. The results in Table 2 (EPISODE) are for a generous set of compositions without much prior knowledge. We consider all compositions up to length 4 (or length 8 for the Bike sharing dataset) except for compositions where two consecutive transition points are inflection points. This ends up being, respectively, 26 and 34 compositions. Thus, EPISODE achieves competitive results for a relatively large choice of compositions and minimum prior knowledge. Often, the results are not significantly different for EPISODE*, where we incorporate such knowledge (see Tacrolimus or Bike sharing dataset). Prior knowledge is, however, useful in reducing the training time, satisfying requirements, and improving extrapolation.\n\n**Actions taken**: We have now clarified the amount of prior knowledge used in experiments in Section 6.\n\n### End-to-end\n\n> Additionally, the realization is not end-to-end, requiring manual verification and editing ...\n\nWe would like to clarify that, although manual verification and editing are big advantages of our approach, they are *not necessary*. Although the training process consists of two steps (training the composition map and then the property maps), the whole pipeline *is* end-to-end. There is *no manual intervention* in experiments performed in Section 6.\n\n**Actions taken**: We have now clarified the absence of manual intervention in experiments in Section 6.\n\n### Evaluation metric\n\nThank you for pointing this out! We have now described the used metric in Sections 5 and 6. The error is given by\n\n$$ \\mathcal{L} = \\frac{1}{D} \\sum_{d=1}^D \\sqrt{\\frac{1}{N_d} \\sum_{n=1}^{N_d} \\left(\\frac{1}{M} ||\\mathbf{F}(\\mathbf{v}^{(d)})(t_n^{(d)})-\\mathbf{y}_n^{(d)}||_2^2\\right)} $$\n\nwhere $\\mathbf{F}$ is the predictive model (EPISODE), $D$ is the number of samples, $M$ is the dimensionality of the system, $N_d$ is the number of measurements of the $d^{\\text{th}}$ sample $(\\mathbf{v}^{(d)}, (t\\_n^{(d)}, \\mathbf{y}\\_n^{(d)})\\_{n=1}^{N\\_d})$. We choose this metric because for $M=1$ it reduces to the standard mean RMSE, i.e., $\\frac{1}{D}\\sum_{d=1}^D \\sqrt{\\frac{1}{N_d} \\sum_{n=1}^{N_d} (F(\\mathbf{v}^{(d)})(t_n^{(d)})-y_n^{(d)})^2}$  (used in Semantic ODE paper) and we normalize by $M$ so that it is easier to compare results between systems of different dimensionality.\n\nThe error is calculated on a held-out dataset of samples, whose trajectories are observed at the same time points as the ones in the training dataset.\n\n**Actions taken**: Described the used metric in Sections 5 and 6.\n\n### Computational cost\n\nWe have now added a [table](https://imgur.com/a/4jhWn1H) to Appendix D showing the overall computational cost of all approaches (including five runs and hyperparameter tuning). We also show the individual times to train the composition map and the property maps. All experiments were performed on an 18-core Intel Core i9-10980XE with 60GB of RAM. \n\nAs mentioned in the limitations, training of the composition map is the most time-consuming process as it requires fitting every admissible composition to each sample (this can be parallelized). Note that we perform this preprocessing *once* for all five seeds and then subsample the results based on the split. A more detailed breakdown of the time needed to train the composition map is provided in the [following table](https://imgur.com/a/90rECtR). We discuss it further in our response to Reviewer evCm (Q1). \n\n**Actions taken:** Added computation times in Appendix D.\n\n### Performance\n\nWe have run additional experiments, including a new dataset with the viral dynamics model of HIV and settings with higher noise. All of the results can be seen [here](https://imgur.com/a/5nnGzqi).\n\nEPISODE* achieves nearly perfect performance on the SIR dataset, demonstrating superior noise robustness to SINDy whose constrained search space should have given it an advantage on this simple problem. Both variants of EPISODE show superior performance on the PK and HIV datasets (in both low and high noise settings) compared to both ODE discovery and black-box models. Both variants also achieve superior performance to ODE discovery methods on the two real datasets (Tacrolimus and Bike sharing).\n\n---\n\nWe hope we have addressed all your concerns. Please let us know if you have any additional questions, we are more than happy to address them.\n\nKind regards,\n\nAuthors", "claims_and_evidence": "The paper claims to provide an alternative to equation-based modeling by learning behavior directly from data, which is supported by a detailed introduction to the algorithm and experimental results on both synthetic and real-world datasets.", "methods_and_evaluation_criteria": "The method seems reasonable but highly relies on the incorporation of semantic inductive biases. Such incorporation benefits the learning process but, at the same time, depends on prior knowledge about the behavior of the dynamical system, raising questions about its scalability to complex systems.\n\nAdditionally, the realization is not end-to-end, requiring manual verification and editing, which raises concerns about efficiency.\n\nThe performance comparison is first mentioned in line 420, but I did not find an introduction to the evaluation metric. Moreover, based on the results shown in Tables 1 and 2, the method does not significantly outperform other baselines.", "theoretical_claims": "The paper provides a solid conceptual framework, and no obvious issues are found with the theorem.", "experimental_designs_or_analyses": "The case study on Tacrolimus is well-motivated. The comparison with different model classes is comprehensive. However, the computation efficiency are not clearly discussed.", "relation_to_broader_scientific_literature": "The work is closely related to recent efforts in physics-informed machine learning and symbolic regression.", "essential_references_not_discussed": "NAN", "other_strengths_and_weaknesses": "See above.", "other_comments_or_suggestions": "See above.", "questions_for_authors": "See above.", "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": "BCQXSMp6js", "reviewer": "Reviewer_htLB", "summary": "The paper tackles the modeling of personalized dynamical systems — that is, dynamical systems whose trajectory evolves conditioned on a set of static parameters, such as the initial condition in ordinary differential equation (ODE) systems, together with other \"personal\" covariates, such as one patient's weight or age, which affect the dynamics of e.g. drug metabolism in the patient. To do this, the paper extends the very recent Semantic ODE framework of Kacprzyk and van der Schaar (2025), which break down the shape of ODE trajectories into a *composition* (a sequence of *motifs*  such as \"increasing and convex\" or \"decreasing and concave\") and their *properties* (the specifics of these motifs, such as the locations and heights of local maxima). \n\nSemantic ODE is trained to predict the easily interpretable semantic composition and its properties that best describe one-dimensional noisy trajectories, given the initial condition of the system. The authors extend Semantic ODEs in two directions:\n\n1. First, they allow for the modeling of systems of dimension larger than one, by making use of a  \"channel independent strategy\" — that is, by treating each component of the target multidimensional system as independent. \n2. Second, they allow the composition map and properties sub-maps to process multidimensional covariates. \n\nThe authors test their methodology in different synthetic and real world datasets that correspond to very different, personalized dynamical systems, and compare their performance against that of symbolic and neural ODE models.\n\n*References*:\n- Kacprzyk and van der Schaar (2025): No Equations Needed: Learning System Dynamics Without Relying on Closed-Form ODEs, ICLR\n\n## Update After Review\n\nThe authors' replies to my questions and concerns were satisfactory. Hence, I increased my score.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Dear Reviewer htLB,\n\nThank you so much for such a comprehensive review. We deeply appreciate your time and attention spent on our paper. We are glad that you found our methodology interesting, comparisons comprehensive, and claims justified. We answer the six weaknesses in your summary, then the two questions. Finally, we address your comment about the incremental nature of our work.\n\n### W1. Datasets with more than one dimension\nWe have now added experiments on the viral dynamics model of HIV [1] ($M=3$) in two noise settings, where the features include not only the initial conditions but also four additional parameters. We have also tested our method on the SIR model with more noise. Overall, we have added three additional multidimensional datasets (all results can be seen [here](https://imgur.com/a/5nnGzqi).\n\nThe main challenge of extending Semantic ODE to multidimensional systems lies in the fact that these trajectories no longer can depend on a single initial condition $x_0\\in\\mathbb{R}$ but on a multidimensional $\\mathbf{x}_0 \\in \\mathbb{R}^M$ (where $M>1$). Thus, the main contribution is the extension to multidimensional inputs, which we demonstrate in all our experiments.\n\nWe hope that this explanation and the additional experiments provide sufficient evidence for our claim.\n\n**Actions taken**: Added additional experiments on multidimensional datasets in Appendix A.\n\n### W2. Details of $F_{\\text{traj}}$\n**Actions taken**: Added description of the trajectory predictor in Appendix E as well as details on compositions, motifs, and properties to make the paper self-contained. \n\n### W3. Similar examples in the introduction\n**Actions taken**: Replaced some of the examples and references in the introduction to limit overlap with the paper on Semantic ODEs.\n\n### W4. Predicted trajectories\n**Actions taken**: Included plots of sample trajectories predicted by EPISODE in Appendix A. These can be seen [here](https://imgur.com/a/IZHZ6l2).\n\n### W5. Metrics used\nPlease see our response to Reviewer g6tH (Evaluation metric).\n\n**Actions taken**: Described the used metric in Section 5.\n\n### W6. Computational cost\n\nWe have now added a [table](https://imgur.com/a/4jhWn1H) to Appendix D showing the overall computational cost of all approaches (including five runs and hyperparameter tuning). We also show the individual times to train the composition map and the property maps. All experiments were performed on an 18-core Intel Core i9-10980XE with 60GB of RAM. \n\nAs mentioned in the limitations, training of the composition map is the most time-consuming process as it requires fitting every admissible composition to each sample (this can be parallelized). Note that we perform this preprocessing *once* for all five seeds and then subsample the results based on the split. The actual time to fit the decision tree is negligible. A more detailed breakdown of the time needed to train the composition map is provided in the [following table](https://imgur.com/a/90rECtR). We discuss it further in our response to Reviewer evCm (Q1). \n\n**Actions taken:** Added computation times in Appendix D.\n\n### Q1 Time horizon\nEach dataset has a specific time horizon that is subsequently scaled to $(0,1)$. The details are in Appendix D.1. However, for unbounded motifs, our model can predict for any time $t\\geq0$. We evaluate each model on a on a held-out dataset of samples, whose trajectories are observed at the same time points as the ones in the training dataset.\n\n### Q2 Robustness\n**Actions taken**: Added additional five experiments on datasets with higher noise ($\\sigma=0.1$) in Appendix A. EPISODE shows robustness to noise in these settings. The results can be seen [here](https://imgur.com/a/3pjzwQS)\n\n### Minor\n**Actions taken**: Fixed the typo on page 4 (it should have been $\\mathbf{v} \\in \\mathcal{V}$).\n\n### Incremental contributions\n\nSemantic ODE only works for one-dimensional trajectories predicted from one-dimensional inputs. We have extended it to more realistic settings. In particular, we\n1. Extended the composition map from interval partition to a decision tree and designed a novel optimization algorithm for fitting it.\n2. Extended the property functions from univariate functions to GAMs and designed a novel optimization algorithm for fitting GAMs such that the predicted properties are consistent with the composition.\n3. Improved composition map training by using Dynamic Time Warping distance.\n4. Accommodated bounded compositions ($t_{\\text{end}}\\in\\mathbb{R}$).\n5. Accommodated categorical variables as features.\n6. Implemented a visualization tool to visualize the fitted GAMs.\n\n**References**\n\n1. Hill, A. L., Rosenbloom, D. I., Nowak, M. A., & Siliciano, R. F. (2018). Insight into treatment of HIV infection from viral dynamics models. Immunological reviews, 285(1), 9-25.\n\n---\n\nWe hope we have addressed all your concerns. Please let us know if you have any additional questions, we are eager to address them.\n\nKind regards,\n\nAuthors", "claims_and_evidence": "The authors claim that their methodology enables practitioners to:\n\n1. readily integrate prior knowledge, which is done via the pre-selection of the composition library (as in Section 5.2);\n2. easily understand the main aspect of the dynamics, through the interpretability of their motif (as in Section 5.3);\n3. ensure desired behavior and revise the model when necessary (as in Section 5.4 and 5.5).\n\nThe authors also claim that their model is flexible enough to be competitive in forecasting tasks. \n\nOverall, their results in Sections 5 and 6 align well with these claims. However, **most of these claims also apply to the original Semantic ODE framework on which this work is based**.  \n\nThe proposal extends Semantic ODE to multidimensional systems with multidimensional covariates. **Although their target datasets do feature different covariates, they mostly consists of one-dimensional processes**. Furthermore, the authors never display their target datasets nor their predictions, so one can only rely on the scores in Tables 1 and 2, **which are also not described**.", "methods_and_evaluation_criteria": "The methodology proposed by the authors is interesting. Their results demonstrate that it represents an alternative to ODE modeling which, as the author point out, does not necessarily lead to easily interpretable representations of the dynamics. \n\n*Regarding evaluation*, I do not manage to find what metric is reported in Tables 1 and 2.  From the introduction one gathers that it must correspond to some divergence between target and predictions up to some horizon time (i.e. a forecasting task). **None of these details are provided** (or at least, I could not find them).", "theoretical_claims": "There are no major theoretical claims in this work.", "experimental_designs_or_analyses": "The organization of Section 5 nicely illustrates how the proposal addresses all the claims made in the introduction. Section 6 compares the proposal against well-established ODE-based methods, and Appendix D contains the details about the target datasets and the training details of all baselines. \n\nUnfortunately, and as I wrote above, there is little to no information regarding the metrics reported in Table 1 and 2.", "relation_to_broader_scientific_literature": "The authors heavily rely on the work by Kacprzyk and van der Schaar (2025), which introduced Semantic ODE. Indeed, the proposal extends Semantic ODE to multidimensional systems with several static covariates, *which required novel implementation of the composition map and properties sub-maps*. These extensions are important and have great applicability. \n\nHowever, only one target dataset has dimension larger than one. The paper would benefit from experiments in other higher dimensional systems — even if such system will often display periodic or periodic-like behavior, which is hard to model with the proposed methodology, as acknowledged by the authors.\n\nFurthermore, there are some details, like the specifics of the trajectory predictor $F_{traj}$, which the authors could include into the Appendix. Otherwise the reader is forced to read Kacprzyk and van der Schaar (2025) to understand the proposal. In my opinion, an ICML paper should be self-contained. Specially given that the work by Kacprzyk and van der Schaar (2025) just came out.\n\nFinally, the introduction section borrows many lines from Kacprzyk and van der Schaar (2025). This is to be avoided.", "essential_references_not_discussed": "**I do not know about any fundamental paper missed by the authors**. On the contrary, the authors  compared their methodology against different ODE-based models, which are widely used by the dynamical system and machine learning community. What is more, since the proposal's most crucial aspects are its transparency, verifiability, and editability (as written by the authors), I do not think it is necessary to compare their method with the most recent, transformer-based forecasting models that do not rely on continuous-time representations.", "other_strengths_and_weaknesses": "I think the main weakness I find is that the paper could be seen as incremental, for it relies too much on the Semantic ODE framework of Kacprzyk and van der Schaar (2025).\n\nLet me summarize the main weaknesses I highlight above:\n1. The authors only study one target dataset of dimension larger than one. This is a weakness because one of the main claims of the authors is that they extend Semantic ODE to multidimensional systems.\n2. The specifics of the trajectory predictor $F_{traj}$ are missing. The paper shoul be self-contained.\n3. The introduction borrows many lines from Kacprzyk and van der Schaar (2025). \n4. The authors do not display the target trajectories nor their predictions. The time series forecasting community often includes such plots, because they help readers better understand the performance of the models under study. It would also illustrate the connection between the motifs predicted by the model and the actual data.\n5. There is no description of the metrics used to compute the scores in Tables 1 and 2 (or at least, I did not find them).\n6. The authors give no details (or I did not find them) regarding the computational cost of their proposal. For example, how long does it take to optimize their models on the different datasets as compared to the baselines?", "other_comments_or_suggestions": "On the second column, first paragraph of page 4 the authors use both $v \\in \\mathcal{V}$ and $\\boldsymbol{v} \\in \\mathcal{V}$ to refer to the covariates. Unless I am misunderstanding something, they always correspond to vectors, right?", "questions_for_authors": "1. How far into the future you predict the dynamics of the systems you study? \n2. The strength of the noise corruption is set to 0.01. How does the model respond to different noise corruptions? Is it as robust as Semantic ODE?", "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:

Meta-review: This manuscript presents a method to identify personalized dynamics through a compositional set of motifs that remove the need for formulating parametrized closed-form equations that are then adapted to each individual. The work is based on decision trees for the compositional maps and generalized additive models. Key to the motivation is the ability to work on multi-dimensional systems (as opposed to prior work in this area that has focused on 1D time series). This is accomplished through separating and treating each dimension separately. 

While the reviewers made note of some weaknesses, such as the potential to scale and limitations in some cases, e.g., periodic systems, these were viewed as minor weaknesses and the authors were able to respond to most of them. I therefore believe that this work would be a good addition to the ICML conference.
Decision: Accept (poster)