Individualized Dosing Dynamics via Neural Eigen Decomposition
Authors: Stav Belogolovsky
Ido Greenberg
Danny Eytan
Shie Mannor
Venue: NeurIPS 2023 poster (2023)
Keywords: personalized medicine
dosing dynamics
sequential prediction
stochastic differential equations
Kalman filter
recurrent neural networks
medical drug control

Abstract: Dosing models often use differential equations to model biological dynamics. Neural differential equations in particular can learn to predict the derivative of a process, which permits predictions at irregular points of time. However, this temporal flexibility often comes with a high sensitivity to noise, whereas medical problems often present high noise and limited data. Moreover, medical dosing models must generalize reliably over individual patients and changing treatment policies. To address these challenges, we introduce the Neural Eigen Stochastic Differential Equation algorithm (NESDE). NESDE provides individualized modeling (using a hypernetwork over patient-level parameters); generalization to new treatment policies (using decoupled control); tunable expressiveness according to the noise level (using piecewise linearity); and fast, continuous, closed-form prediction (using spectral representation). We demonstrate the robustness of NESDE in both synthetic and real medical problems, and use the learned dynamics to publish simulated medical gym environments.

Individualized Dosing Dynamics via Neural Eigen Decomposition
Stav Belogolovsky" stav. bww.@gmail......mmmmoom......i.m com
Ido Greenberg" gido@campus technion ac i1
Danny Eytan aaliary.colicgaaan com
Shie Mannor shie@ee. technion ac.il ii
Abstract
Dosing models often use differential equations to model biological dynamics. Neu- ral differential equations in particular can learn to predict the derivative of process, which permits predictions at irregular points of time. However, this temporal flex- ibility often comes with high sensitivity to noise, whereas medical problems often present high noise and limited data. Moreover, medical dosing models must generalize reliably over individual patients and changing treatment policies. To address these challenges, we introduce the Neural Eigen Stochastic Differential Equation algorithm (NESDE). NESDE provides individualized modeling (using a hypernetwork over patient-level parameters); generalization to new treatment policies (using decoupled control); tunable expressiveness according to the noise level (using piecewise linearity); and fast, continuous, closed-form prediction (us- ing spectral representation). We demonstrate the robustness of NESDE in both synthetic and real medical problems, and use the learned dynamics to publish simulated medical gym environments.
Introduction
Sequential forecasting in irregular points of time is required many real-world problems, such as modeling dosing dynamics of various medicines (pharmacodynamics)) Consider patient whose physiological or biochemical state requires continuous monitoring, while blood tests are only available with limited frequency. Pharmacodynamics models often rely on an ordinary differential equation models (ODE) for forecasting, Additional expressiveness can be obtained via customized learned models, such as neural-ODE, which learns to predict the derivative of the process [Chen et al. 2018 Liu et al. 2019]. By predicting the derivative, neural-OD can make irregular predictions at flexible time-steps. unlike regular models that operate in constant time-steps (e.g., Kalman filter, Kalman 1960) and recurrent neural networks, Rumelhart et al. 1986]).
However, real-world forecasting remains challenge for several reasons. First, the variation between patients often requires personalized modeling. Second, neural-ODE methods are often data-hungry: they aggregate numerous derivatives provided by non-linear neural network, which often sensitive to noise. Training over large dataset may stabilize the predictions, but data is often limited. Third, most neural-ODY methods only provide a point-estimate, while uncertainty estimation often critical in medical settings. Fourth, for every single prediction. the neural-ODD runs a numeric ODE solver,
Department of Electrical and Computer Engineering, Technion-Israd Institute of Technology, Haifa, Israel Department of Physiology and Biophysics, Faculty Medicine, Technion Israel Institute of Technology, Haifa, Israel Nvidia Research
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

Table I: A summary of the features of NESDE. In Sections
NESDE is tested on synthetic and real medical data. In addition, each component of NESDE studied experimentally specified. Challenge
6
Solution
Experimental evaluation
Hyper-network with high-level features input
Hyper-network ablation (Appendix E.1
Individualized modeling
Regularized dynamics: piecewise-linear with complex eigenvalues
Varying train size (Section Appendix E.2 varying sparsity (Appendix
Sample efficiency
E.3
Probabilistic Kalman filtering
Uncertainty estimation
NLL evaluation (Sections
Spectral representation with closed-form solution
Time flexibility (Appendix interpretability (Appendix
E.
Fast continuous inference
Decoupling control from other inputs
Out of distribution control (Section Appendix (E6
Control generalization
along with multiple neural network calculations of the derivative. This computational overhead in inference may limit latency-sensitlv applications.
fifth challenge comes from control. the framework of retrospective forecasting control signal (drug dosage) is often considered part of the observation [De Brouwer et al., 2019). However, this approach raises difficulties the control is observed at different times or more frequently than other observations. If the control is part of the model output, it may also bias the train loss away from the true objective. Finally, by treating control and observations together, the patterns learned by the model may overfit the control policy used in the data and generalize poorly new policies.
Generalization to out-of-distribution control policies is essential when the predictive model supports decision-maker as the control policy may be affected by the model. Such decision-maken is an important use-case of sequential prediction: model- based reinforcement learning and control problems require reliable model [Moerland et al. 2020 Angermueller et al. 2019)) in particular in
risk-sensitive control Yuu et al. 2021 Greenber et al. 2022 Greenberg and Mannor 202T]
Section 4 introduces the Neural Eigen-SDE algorithm (NESDE) for continuous forecasting, which is designed address the challenges listed above. NESDE uses a hypernetwork to provide an individu- alized stochastic differential equation (SDE), regularized to be piecewise-linear and represented in spectral form. The SDE derives probabilistic model similar to Kalman filtering, which provides uncertainty estimation. Finally, the SDE decouples the control signal from other observations, to discourage the model from learning control patterns that may be violated under -fi-trriiittiitoo control policies. Table I summarizze all these features.
Section tests NESDE against both neural-ODE methods and recurrent neural networks. NESDE demonstrates robustness to both noise (by learning from little data) and outtr--distiibution control policies. In Appendix [F.5] the SDE model of NESDE is shown enable potential domain knowledge and provide interpretability via the predicted SDE cigenvalues. Appendix E.4 demonstrates the disadvantage of discrete methods in continuous forecasting.
In Section 66] NESDE demonstrates high prediction accuracy in two medical forecasting problems with noisy and irregular real-world data: (1) blood coagulation prediction given Heparin dosage, and (2) prediction of the Vancomycin (antibiotics) levels for patients who received it.
Contribution:
We characterize the main challenges in continuous forecasting for medication dosing control.
We design the novel Neural Eigen-SDE algorithm (NESDE), which addresses the challenges as summarized in Table I and demonstrated empirically in Sections and
We use NESDE to improve modeling accuracy in two medication dosing processes. Based on the learned models, simulate gym environments for future research of healthcare control.

---

1.1 Related Work
Classic filtering: Classic models for sequential prediction in time-series include ARIMA mod- els [Moran and Whittle [95T] and the Kalman filter (KF) [Kalman| T960]. The KF provides probabilistic distributions and in particular uncertainty estimation. While the classic KF is limited to linear dynamics, many non-linear extensions have been suggested [Krishnan et al. 2015 Coskun et al., 2017 Revach et al.) 2021 Greenberg et al. [021]. However, such models are lypically limited constant prediction horizon (lime-step. Longer-horizon predictions are often made by applying the model recursively (Herrera et al. 2007] Bontempi et al. 2013). This poses a significant challenge to many optimization methods [Kolen and Kremer, 2001 as also demonstrated in Appendix [E.4]
Limited types of irregularity can also be handled by KF with intermittent observations
Park and
Sahai 2011 Sinopoli et al. 2004 periodical time-steps (Li al. 2008
Recurrent neural networks: Sequential prediction is often addressed via neural network models, relying on architectures such as RNN [Rumelhart et al., 1986). LSTM (Hochreiter and Schmidhuber 1997] and transformers [Vaswani et al. 2017]. LSTM, for example, is a key component in many SOTA algorithms for non-linear sequential prediction (Neu et al. 202T]. LSTM can be extended to a filtering framework to alternately making predictions and processing observations, and even to provide uncertainty estimation [Gao et al. [019]. However, these models are typically limited to constant time-tteps, and thus suffer from the limitations discussed above.
Neural-ODE models: Parameterized ODE models can be optimized by propagating the gradients of a loss function through an ODE solver (Chen et al., 2018 Liu et al.] 2019 Rubanova et al., 2019]. By predicting the process derivative and using an ODE solver in real-time, these methods can choose the effective time-steps flcxibly. Uncertainty estimation can be added via process variance prediction [De Brouwer et al., 2019]. However, since neural-ODE methods learn a non-linear dynamics model, the ODE solver operates numerically and recursively on top of multiple neural network calculations. This affects running time, training difficulty and data efficiency as discussed above. While neural-ODE models have been studied for medical applications with irregular data [Lu let al. [222], simpler models are commonly preferred in practice. For example, the effects of Heparin on blood coagulation is usually modeled by either using discrete models (Nemati et al. 2016 or manually based on domain knowledge [Delavenne al. 2017].
Our method uses SDE with piecewise linear dynamics (note this is different from piecewise linear process). The linear dynamics per time interval permit efficient and continuous closed-fore forecasting of both mean and covariance, |Schirmer et al. [2022) also rely on a linear ODE model, but only support operators with real-valued eigenvalues (which limits the modeling of periodic processes), and do not separate control signal from observations (which limits generalization to out- of of-distribution control). Our piecewise linear architecture, tested below against alternative methods including De Brouwer et al. [2019) and Schirmer et al. [2022]. is demonstrated to be more robust to noisy, sparse or small datasets, even under out-of out-of-distribttion control policies.
2 Preliminaries: Linear SDE
We consider particular case of the general linear Stochastic Differential Equation (SDE):
dx = xt) + ū(t)] + (t)))
(1)
where : + R" is time-dependent state; R=x* is fixed dynamics operator; ū :R + R
is the control signal: and dW: :R + R is Brownian motion vector with covariance Q E R
General SDEs can be solved numerically using the first-order approximation (x At (t) or using more delicate approximations (Wang and Lin. 1998). The linear SDE, however, and in particular Eq. ( can be solved analytically [Herzog] 2013]:
[ [')(()- S ('
X(t) (t) (p(to)--x(6)
(2)
where (to) is an initial condition, and (+L the eigenfunction of the system. More specifically, if V is the matrix whose columns {2 are the cigenvectors of A, and A is the diagonal matrix
3

---

whose diagonal contains the corresponding eigenvalues
= {Ai then
-(01 (1
@(1) - Ve
Ui ed
Un edta
(3)
If the initial condition given X ((6) Nuu, Eo), Eq. (() N (n((), E()))
becomes
(4)
(tt) (t) (b(tt) uo 2(T) ((()d-
21t) (())
to
where ('(t) ((to) - E (tt)) J ((T) (())))((((()-))
Note that if'v X < 0 and 0, have (1)) Eq. (4) may produce a complex solution; Appendix exxplain how to use a careful parameterization
1000 0 (stable system). In addition, if X complex,
only calculate the real solutions.
3 Problem Sctup: Sparsely-Observabbe SDE
We focus on online sequential prediction of process (t) IR... To predict (1a) a cer- tain lu. we can use noisy observations (t) (at given times t < tu), well as control signal u(t) E R (vt < to); offline data of Y and u from other sequences; and one sample of con- textual information C' E RRde per sequence (cap- turing properties of the whole sequence). The dynamics of Y are unknown and may vary between sequences. For example, sequences may represent different patients, each with its own dynamics; C' may represent patient infor- mation; and the objective is "zero-shot" learning upon arrival of a sequence of any new patient. both irregular and sparse: they are received at arbitrary points of time, and are sparse in comparison the required prediction frequency (i.e., continuous forecasting, as illustrated in Fig. (D].
SDE
SDE (control pack at (-
observettio
amo
time
(a)
(b)
Figure 1: Samples of sparsely observed SDEs: the Brownian noise and the sparse observations pose jor challenge learning the underlying SDE dynamics. Efficient learning from external trajectories data re- quired, as the current trajectory often does not contain sufficient observations.
addition, the observations within sequence are
To model the problem, we assume the observations (t) to originate from an unobservable latent process (t) R (where n > is hyperparameter). More specifically:
dx(()= = Fe (1((), (()). Y(1) X())mm 4(1) Y(t) + vo)
(5)
where Fc is stochastic dynamics operator (which may depend on the context C); Y is simply the first coordinates of X; Y is the corresponding observation; and (t) is its i..d Gaussian noise with cero-mean and (unknown) covariance Rç E Rm (which may also depend on C¹). Our goal is predict Y, where the dynamics Fe are unknown and the latent subspace of X is unobservable. In cases where data of Y is not available, we measure our prediction accuracy against Ý The control u(t) is modeled separately from Y, is not part of the prediction objective, and does depend on X.
4 Neural Eigen-SDE Algorithm
Model: In this section, we introduce the Neural Eigen-SDE algorithm (NESDE, shown in Algorithm ! and Fig. 2. NESDE predicts the signal (1) of Eq. (5) continuously any required point of time T.
It relies piecewise linear approximation which reduces Eq. (5) into Eq. (I:
(6) where I, (t,, (44) is time interval, dW is Brownian noise with covariance matrix Qi, and A, E guxn B Rx Qi E ax E 1 form the lincar dynamics model corresponding to the
Vt E. dx((t) = [A (X(t) - a) + B- u(t)] dW(()

---

Output
Inputs
NESDE
Mila + 1206)
Ohsercation
lflor
Hypermet Imputs
Hypernet
Control
Inital coodittion
u(t)
Elgan 5DD
pilC 89)
V.A.O Run.
Cutput
V.A Q.B
Context
Heparet
Updane dynamics Ie..y
prediction
Time
Time
Figure 2: NESDE algorithm: Hypernet determines the SDE param- eters; Eigen-SDE solver uses them to predict the next time-interval;; the filter updates the state upon arrival initiates a new interval. For more frequent updates of the dynamics, the initial condition becomes the last prediction (dotted arrow).
Figure 3: Zoom into the individualiza- tion module of NESDE: Hypernet uses the context and the estimated state determine the SDE parameters for the solver.
new observation, which
interval I,. terms of Eq. ! we substitute A = A; and Bu A.a. Note that if A, is stable system and u = 0, the asymptotic state is ju(t) 1-000, a. To solve Eq. (6 within every I,, NESDE has to learn the parameters (A, Qi,
The end of - 1 typically represents one of two events: either an update of the dynamics A (allowing observation. A new observation time 4,
the piecewise linear dynamics), the arrival of triggers update of X (t:) according the conditional distribution X (4+)) (t:) (this is particular case of Kalman filtering, as shown in Appendix A). Then, the prediction continues for Z; according Eq. 66. Note that once (to) is initialized have a Normal distribution, it remains Normally- distributed throughout both the process dynamics (Eq. 4 and the observations filtering (Appendix(A) This allows NESDE to efficiently capture the distribution of X (t). where the estimated covariance represents the uncertainty.
Eigen-SDE solver (ESDE) spectral dynamics representation: A key feature of NESDE is that A, only represented implicitly through the parameters V, X defining its eigenfunction ( of Eq. (3) (we drop the interval index i with a slight abuse of notation). The spectral representation allows Eq. (4) to solve X (tt analytically for any & E Z; at once: the predictions are not limited to predefined times, and do not require recursive iterations with constant time-steps. This particularly useful in the caarsely-observable setup of Section 3 as it lets numerous predictions be made at once without being "interrupted" by a new measurement.
The calculation of Eq. requires efficient integration. Many SDE solvers apply recursive numeric integration [Chen et al. 2018 De Brouwer al., 2019]. In NESDE, however, thanks the spectral decomposition, the integration only depends on known functions of instead of ((t) (Eq. (4h), hence recursion is not needed, and the computation constant over an interval I; (or has any other maalyically-integraal form), Appendix B shows how calculate the integration analytically. Piecewise constant is common, for example, when the control is updated along with the observations.
be paralleled. Furthermore, if the control is
In addition to simplifying the calculation, the spectrum of A, carries significant meaning about the dynamics. For example, negative cigenvalues correspond to a stable solution, whereas imaginary ones indicate periodicity, demonstrated in Fig, 4 Note that the (possibly-comples cigenvalues must be constrained to represent a real matrix A;- The constraints and the calculations in the complex space detailed in Appendix
Note that if the process X (t) actually follows the piecewise linear model of Eq. (6), and the model parameters are known correctly, then the Eigen-SDE solver trivially returns the optimal predictions. The complete proposition and proof are provided Appendix
Propasition (Eigen-SDie SDE solver optimality). If (t) follows Eq. (6) with the same parameters used by the Eigen-SDE solver, then under certain conditions, the solver prediction at any point of time optimizes both the expected error and the expected log- ligllikcl
5

---

NESDE
prediction
prediction
prediction
-5
-10
NESDE
-15
10
time
time
time
(a) Complex A
(b) Real A
(c) Imaginary A
Figure 4: Sample trajectories with different types of dynamics (and random control signal that is not displayed).
As demonstrated in Appendix E.S NESDE directly estimates the dynamics À and reveals their type.
Proof Sketch. If the process X(t) follows Eq. 6 the Eigen- SDE solver output corresponds the true distribution (() ~ (u(()), 2()) for any 1. Thus, both the expected error and the expected log-likelihood are optimal.
Updating solver and filter parame- ters: NESDE provides the parameters V,A, Q, B, a to the Eigen-SDE solver, well the noise R to the observation filter. As NESDE assumes a piecewise linear model, it separates the time into intervals I = (t,, 1++1) (the interval length is hyperparameter), and uses dedicated model to predict new parame- ters at the beginning ( of every interval.
Algorithm 1 NESDE
Input: context C: control signal u(t); update times E R², prediction times (Pz, HI,E Initialize: μu, 2, a, < Prior(C)
for I; in I: do V. Q B, a, R +- Hypernet(C,,4.E)
for t in Pz, do
Pt: E, < ESDE( E, u, L,,,,;V,,,Q B) // Eg.. (4)
predict: Y, N a, E, R)
The model receives the current state (4, and the contextual information C, and returns the parameters for I, Specifically, use Hypernet [Ha al.) 2016]: one neural network 91 with net- work parameters O, returns the weights of a second network G2. i.e., W = (C: Then, 92 maps the current state estimation to a set of dynamics parameters: (V, X, Q, B. a, R) (2(:;; W) 92(X; 91 (C; θ)) (Fig. 3). The dynamics (2(X; W) may be updated more often than the individualized weights W In (C; O). as the context C' may change with lower frequency. For the initial state, where X is unavailable,
if given observation 8, then Σ < Filter(u,, E,, R. V.) 1/ Appendix A end if
end for end for
learn state prior from C/ by a dedicated network; this prior helps NESDE to function as a "zero-shot" model.
The Hypernet module implementation gives control over the non-linearity and non-stationarity of the model. The importance of the Hypernet module for individualized modeling is demonstrated in ablation test in Appendix E.1 In our current implementation, only V, X, Q are renewed every time interval. a (asymptotic signal) and R (obseryation noise) are only predicted once per sequence, assume they independent of the state. The control mapping B is assumed to be global parameter.
Training:
The learnable parameters of NESDE are the control mapping B and Hypernet's parameters (which in turn determine the rest the parameters). To optimize them, the training relies on dataset of sequences of control signals {ueeg (t;) {)] eequj and (sparser and possibly irregular) states and observations (YYeeg())) Yseg (en (if is not available, we use Y instead the training target).
The latent space dimension and the model-update frequency At are determined as hyperparameters. Then, we use the standard Adam optimizer [Diederik P. Kingma 2013) to optimize the parameters with respect to the loss NLL(j)
log P(Y (t,) (T((t, (I,,)) (where ,l, Σ predicted by NESDE sequentially from u, Y). Each training iteration corresponds to a batch sequences of data, where the NLL is aggregated over all the samples of the sequences. Note that our supervision for the training is limited to the times the observations, even if we wish to make more frequent predictions in inference.

---

Table 2:
Test errors the irregular synthetic benchmarks, estimated over 5 seeds and 1000 test trajectories per seed, with standard deviation calculated across seeds.
Complex dynamics cigenvalues
Real dynamics eigenvalues
Model
MSE
OOD MSE
MSE
OOD MSE
LSTM
0.23 = 0.001
0.589 + 0.02
0.381 + 0.002
2.354 = 0.84
GRU-ODE-Bayee
0.182 = 0.0004
0.361 = 0.044
0.219 = 0.0004
0.355 0.005
CRU
0.233 + 0.0054
0.584 + 0.009
0.231 + 0.001
0.541 + 0.026
NESDE (ours)
0.176 0.0001
0.178 + 0.001
0.222 + 0.0005
0.332 0.005
As demonstrated below, the unique architecture of NESDE provides effective regularization and data efficiency (due to piecewise linearity). along with tunable expressiveness (neural updates with controlled frequency). Yet, it is important to note that the piecewise linear SDE operator does limit the expressiveness of the model (e.g.. in comparison to other neural-OD -DDE models). Further, NESDE is only optimal under restrictive set of assumptions, specified in Proposition
5 Synthetic Data Experiments
In this section, test three main aspects of NESDE: (1) prediction from partial and irregular observations, (2) robustness to out-of-distribution control (OOD), and (3) sample efficiency. We experiment with data of simulated stochastic process, designed to mimic partially observable medical processes with indirect control.
The simulated data includes trajectories of 1 i-dimensional signal Y. with noiseless mea- surements at random irregular times. The goal is to predict the future values of Y given its past observations. However, Y is mixed with a latent (unobservable) variable, and they follow linear dynamics with both decay and periodicity (i.c., complex dynam- ics eigenvalues). In addition, we observe control signal that affects the latent variable (hence affects Y, but only indirectly through the dynamics). The control negatively corre- lated with the observations: u, b.-0.5-V. Y,. by ~ O0, 0.5) is a piecewise constant addi- tive noise (changing 10 times per trajectory).
0.5
2.8
0.4-
NESDE
1.6
NESDE
GRU-ODE
GRU-ODE
0.3
LSTM
MSE
0.8
LSTM
CRU
CRU
0.2
0.2
6
Observation within trajectory
Observation within trajectory
Same control distribution (b) Out of distribution control
Figure 5: MSE vs. number of observations so far in the trajectory, in the complex dynamics setting, for: (a) standard test set, and (b) test set with out-of-distribution control policy. 95% confidence intervals are calculated over seeds.
As baselines for comparison, we choose recent SDE-based methods that provide Bayesian uncertainty estimation: GRU-ODE-Bayer [De Brouwer et al. 2019) and CRU (Schirmer et al. 2022 (despite terminology, both implement an SDE rather than ODE). In these methods, concatenating the control signal to the observation results in poor learning, the control becomes part of the model output and dominates the loss function. To enable effective learning for the baselines, we mask-on out the control from the loss. Additionally, design dedicated LSTM model that supports irregular predictions, as described in Appendix F22
Out-of-distributin control (OOD): We simulate two benchmarks one with complex eigenvalues and another with real cigenvalues (no periodicity). We train all models on a dataset of 1000 random trajectories, and test a separate dataset with different trajectories that follow the same distribution. In addition, we use an OOD test dataset, where the control is positively correlated with the observations: u, b; + 0.5 Y, This can simulate, for example, forecasting of the same biochemical process after changing the medicine dosage policy.
2 and Fig. 5a summarize the prediction errors. Before changing the control policy, NESDE achieves the best accuracy in the complex dynamics, and is on par with GRU-ODE-Bayes in the real dynamics. Notice that CRU, which relies on a real-valued linear model in latent space, is indeed sub-optimal under the complex dynamics, compared to NESDE and GRU--DE-Bayes The LSTM presents high errors in both benchmarks.
7

---

Once the control changes, all models naturally deteriorate. Yet, NESDE presents the smallest deterioration and best accuracy in the OOD test datasets for both complex and real dynamics. In particular, NESDE provides a high prediction accuracy after mere 2 observations (Fig. [b),, making a useful zero-shot model. The robustness to the modified control policy can be attributed to the model of NESDE in Eq. ( which decouples the control from the observations.
In similar setting in Appendix E.6 the control used in the training data has continuous knowledge of Y. Since the model only observes Y in a limited frequency, u carries additional information about Y. This results in extreme overfitting and poor generalization to different control policies for all methods except for NESDE, which maintains robust OOD predictions in this challenging setting,
Sample efficiency: We train each method over datasets with different number of trajectories. Each model is trained on each dataset separately until convergence. As shown in Fig. 6 NESDE achieves the best test accuracy for every training dataset, and learns reliably even from few as 100 trajectories. The other methods deteriorate significantly in the smaller datasets. Note that in the real dynamics, LSTM fails regardless of the amount of data, as also reflected in Table 2
GRU-ODE-Bayes achieves the best sam- ple efficiency among the baselines. In Appendix E.2 we use a benchmark from the study of GRU-ODE-Bayer itself [De Brouwer et al., 2019), and demonstrate the superior sample efficiency of NESDE in that benchmark as well. Appendix extends the notion of sample efficiency to sparse trajectories: for a constant number of training trajectories, it reduces the num- ber of observations per trajectory. NESDE demonstrates high robustness to the amount of data in that setting as well.
1.6
1.6
NESDE
NESDE
GRUODE
GRUODE
MSE
0.8
CRU
LSTM
MSE
0.8
CRU
LSTT
0.4
F.3
200
600
800
200
400
500
800
Train size
(a) Complex dynamics
(b) Real dynamics
Figure 6: Test MSE train data size. 95% confidence
intervals calculated over 1000 test trajectories.
Regular LSTM: Appendix E.4 extends the experiments for regular data with constant time-steps. In the regular setting, LSTM provides competitive accuracy when observations are dense. However, if the signal is not observed every time-step, LSTM deteriorates significantly in unobserved time-steps (Fig. 12c,, possibly because gradients have to be propagated over many steps. Hence, even in regular settings, LSTM struggles to provide predictions more frequent than the measurements.
Medication Dosing Regimes
As discussed in Section | many medical applications could potentially benefit from based methods. Specifically, we address medication dosing problems, where observations are often sparse, the dosing is a control signal, and uncertainty estimation is crucial. We test NESDE on two such domains, against the same baselines as in Section 5 (GRU-ODE-Bayee CRU and an irregular LSTM). We also add naive model with "no-dynamies" (predicts the last observed value).
The benchmarks in this section were derived from the MIMIC-VV dataset [Johnson et al., 2020]. Typically to electronic health records, the dataset contains a vast amount of side-information (c.g., weight and heart rate). We use some of this information as an additional input for each model according its structure ccontext-featuree for the hyper-networt of NESDE, covariates for GRU- state variables for CRU, and embedding units for the LSTM). Some context features correspond to online measurements and occasionally updated. In both domains, we constraint the process eigenvalues À to be negative, to reflect the stability of the biophysical processes. Indeed, the spectral representation of NESDE provides with natural way incorporate such domain knowledge, which often cannot be used otherwise, For all models, in both domains. we use 60-10-30 train-validttiontees data partition. See more implementation details in Appendix F
6.1 Unfractionated Heparin Dosing
Unfractionated Heparin (UH) is widely used anticoagulant drug. It may be given in continuous infusion to patients with life-threatening clots and works by interfering with the normal coagulation cascade. As the effect is not casily predicted, the drug's 's actual activity on coagulation is traditionally
8

---

Table Test mean square errors (MSE) negative log-likelihood (NLL, for models that provide probabilistic prediction) in the medication-dosing benchmarks.
UH Dosing
Vancomycin Dosing
Model
MSE
NLL
MSE
NLL
Naive
613.3 = 13.48
112.2 + 16.4
LSTM
482.1 6.52
92.89 + 11.3
GRU-ODE-Bayes
491 6.88
4.52 + 0.008
80.54 + 11.8
6.38 0.12
CRU
450.4 8.27
4.49 0.012
76.4 + 12.8
3.87 + 0.2
NESDE (ours)
411.2 + 7.39
4.43 0.01
70.71 = 12.3
3.69 = 0.13
monitored using a lab test performed on a blood sample from the patient: activated Partial Thrombo- plastin Time (aPTT) test. The clinical objective is to keep the aPTT value in certain range. The problem poses several challenges: different patients are affected differently; the aPTT test results are delayed; monitoring and control are required in higher frequency than measurements; and deviations of the aPTT from the objective range may fatal. particular, underdosed UH may cause clot formation and overdosed UH may cause internal bleeding [Landefeld et al. 1987]. Dosage rates are manually decided by a physician, using simple protocols and trial-and-erro over significant amounts of time. Here we focus on continuous prediction key component for aPTT control.
Following the preprocessing described in Appendix F.1 the MIMIC-IV dataset de- rives 5866 trajectories of a continuous UH control signal, an irreuularly-obsevve aPTT signal (whose prediction is the goal), and 42 context features. It is known that UH does not affect the coagulation time (aPTT) directly (but only through other unobserved processes, Delavenne et al. 2017);; thus, we mask the control mapping to have no direct effect on the aPTT met- ric, but only the latent variable (which can be interpreted as the body UH level). The control (UH) and observations (aPTT) are one-dimensional (m = 1), and we set the whole state dimension to 4.
70
preditioon coservaters
Dose (g) vit Dlood (ugimil
SPTT [se]
predictions obsonationn
Ull [U/min]
50
1000
Time [min]
Time [min]
(a) Heparin
(b) Vancomycin
Figure 7: A sample of patients from (a) the UH dosing dataset, and (b) the VM dosing dataset. The lower plots correspond to medication dosage (UH (a) and VM in (b)). The upper plots correspond to the continuous prediction of NESDE (aPTT levels in (a) and VM concentration in (b)), with 95% confidence intervals. In both settings, the prediction every point relies on all the observations up that point.
6.2 Vancomycin Dosing
Vancomycin (VM) is an antibiotic that has been in use for several decades. However, the methodology of dosing VM remains subject of debate in the medical community [Rybak et al. 2009]] and there is
a significant degree of variability in VM dynamics among patients [Marsot et al. 2012] The dosage of VM is critical; it could become toxic if overdosed [Filippone al. [017], J, and ineffective as
antibiotic if underdosed. The concentration of VM in the blood can be measured through lab test, but these tests often infrequent and irregular, which fits into our problem setting.
Here, the goal is to predict the VM concentration the blood at any given time, where the dosage and other patient measurements are known. Following the preprocessing described in Appendix F. dataset derives 3564 trajectories of VM dosages at discrete times, blood concentration of VM (m 1) at irregular times, and similarly to Section 6.1 42 context features. This problem is less noisy than the UH dosing problem, as the task is to learn the direct dynamics the VM concentration, and not the effects of the antibiotics. The whole state dimension is set to 2, and also mask the control mapping B to have no direct effect on the VM concentration, where the latent variable that directly affected could be viewed as the amount of drug within the whole body, which in turn affects the actual VM concentration in the blood.
6.3 Results
Fig. displays sample trajectories predicted by NESDE in both domains. As summarized in Table 3 NESDE outperforms the other baselines in both UH and VM dosing tasks, terms of both square errors (MSE) and likelihood (NLL). For the UH dosing problem, Fig. g also presents the errors

---

vs. prediction horizon (the time passed since the last observation). Evidently, NESDE provides the best accuracy in all the horizons. Its smoothness allows NESDE obtain high accuracy in short horizons, and its robustness allows generalization to aut-of-distributinn horizons: while most of the data corresponds horizons of 5-7 hours (see Fig. 14 the appendix), NESDE provides reliable prediction at other horizons as well. By contrast, LSTM and GRU-ODE-Bayes have difficulty with short horizons; they only become competitive with the naive model after hours. CRU provides more robust predictions, but is still outperformed by NESDE.
Despite the large range of aPTT levels in the data (e.g., the top 5% are above 100s and the bottom 5% are below 25s), 50% of all the predictions have errors lower than 12.4s an accuracy level that is considered clinically safe. Fig. g shows that indeed, up to hours after the last lab test, the average error is smaller than 10s.
Error [s]
NESDE
GRU-ODE
Running time: All Heparin and Vancomyein experi- ments were run on single Ubuntu machine with eight 19-10900X CPU cores and Nvidia's RTX A5000 GPU. A single training process of NESDE required 8-24 hours (with variance mostly attributed to early stopping cri- terion). LSTM required similar running times (6-24 hours), whereas GRU-ODE-Bayee and CRU required 36-72 hours to train.
LSTM
CRU
1000
1200
1400
Time since last observation [m]
Figure aPTT prediction errors in the UH problem, vs. the time passed since the last aPTT test.
7 Conclusion
Motivated by medical forecasting and control problems, we characterized a set challenges in modeling continuous dosing dynamics: sample efficiency, uncertainty estimation, personalized modeling, continuous inference and generalization to different control policies. To address these challenges, introduced the novel NESDE algorithm, based stochastic differential equation with spectral representation. We demonstrated the reliability of NESDE in a variety of synthetic and real data experiments, including high noise, little training data, contextual side-ing oile-ifformtiion and out-of-distributie control signals. In addition, NESDE demonstrated high prediction accuracy after few observations, making a useful zero-shot model.
applied NESDE to two real-life high-noise medical problems with sparse and irregular measure- ments: (I) blood coagulation forecasting in Heparin-treated treated patients, and (2) Vancomycin levels prediction in patients treated by antibiotics. In both problems, NESDE significantly improved the prediction accuracy compared alternative methods. This demonstrates NESDE's s advantage in continuous forecasting of controlled SDE systems, which include medication dosing and controlled chemical and biological processes.
As demonstrated in the experiments, NESDE provides robust, reliable and uncertainty- continu- ous forecasting. This paves the way to development of decision making in continuous high-noise decision processes, including medical treatment, finance and operations management. Future research may address medical optimization via both control policies (e.g., to control medication dosing) and sampling policies (to control measurements timing, e.g., of blood tests).
Acknowledgements:: This research was supported by VATAT Fund to the Technion Artificial Intelligence Hub (Tech.AI).
10

---

References
Christof Angermueller, David Dohan, David Belanger, Ramya Deshpande, Kevin Murphy, and Lucy Colwell. Model-based reinforcement learning for biological sequence design. In International conference on learning representations, 2019.
Gianluca Bontempi, Souhaib Ben Taieb, and Yann-Adl Le Borgne. Machine learning strategies for time series forecasting. Lecture Notes in Business Information Processing, 138, 01 2013. doi:
10.10077997--34442-3338844888
Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud. Neural ordinary differential equations. arXiv preprint 18X.vv1806.07366 2018.
Huseyin Coskun, Felix Achilles, Robert DiPietro, Nassir Navab, and Federico Tombari. Long short-term memory kalman filters.recurrent neural estimators for pose regularization. ICCV, 2017.
URL wttttt::ttgtthui com//elleucialsstnke ICCV2017
Edward De Brouwer, Jaak Simm, Adam Arany, and Yves Moreau. Gru-ode- bayes: Continuous modeling of cawuudicall--bbeeeeee time series. H. Wallach, H. Larochelle. A. Beygelzimer, F. d' d'Alché-Buc, E. Fox, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019. URL httpltttts:/ppcccdiinnee
pdf
X. Delavenne, E. Ollier, S. Chollet, F. Sandri, J. Lanoiselée. S. Hodin, A. Montmartin, J.-F. Fuzellier, P. Mismetti, and L. Gergelé. Plammackkineiichhamacyyyauuaaua dosing during cardiopulmonary bypass. BJA: British Journal of Anaesthesia, 118(5):705-712 May 2017. ISSN 0007-0912. doi: 10. 10.1099/bja/aax0044, URL WtRphtpps: //doi. org/10. 1099/bba.09///aa//// aex044
model for unfractionated heparin
Jimmy Ba Diederik P. Kingma. Adam: A method for stochastic optimization. ICLR, 2015. URL
a/12/ 6980
Morri L. Eaton. Multivariate Statisties: Vector Spuce Approach. John Wiley and Sons, 1983.
Edward J Filippone, Walter K Kraft, and John Farber. The nephrotoxicity of vancomycin. Clinical Pharmacology Therapeutics, 102(3):459-485 469, 2017.
Chang Gao, Junkun Yan, Shenghua Zhou, Bo Chen, and Hongwei Liu. Long short-term memory- based recurrent neural networks for nonlinear target tracking. Signal Processing, 164, 05 2019. doi: 10. hitp:ijsippro.201..055
Ido Greenberg and Shie Mannor. Detecting rewards deterioration in episodic reinforcement learn- ing. In Proceedings of the 38th International Conference on Machine Learning, volume 139, pages 13842-38533 PMLR, 18-24 Jul 2021. URL http:tttss//prroeediii mlr A....peess/v139
Freenberg21a.htmm
Ido Greenberg, Netanel Yannay, and Shie Mannor. Noise estimation is not optimal: How to use kalman filter the right way. arXiv preprint luXiv:2104.02372, 2021.
Ido Greenberg, Yinlam Chow, Mohammad Ghavamzadeh, and Shie Mannor. Efficient risk-averse reinforcement learning. Advances in Neural Information Processing Systems, 2022.
David Ha, Andrew Dai, and Quoc V Le. Hypernetworks. arXiv preprint arXiv:1609.09106, 2016.
Luis Herrera, Hector Pomares, I. Rojas, Alberto Guillén, Alberto Pricto, and Olga Valenzuela. Recursive prediction for long term time series forecasting using advanced models. Neurocomputing. 70:2870-2880 10 2007. doi: 10. 00.10166..eeccom.2006.444444 20..0.....
Florian Herzog, Stochastic differential equations. 2013. URL htt./tppee///eeaaccc//aaaaaa///e http://mssa....oooooooooooaaaaaa
Etocchastic-Systems/SS5
Sepp Hochreithreeea and Jurgen Schmidhuber. Long short-term memory. Neural Computation, 1997. URL
https://direct. http://iirrect.....eeeue http:///ww....l.ceoocccoocoocata Short- Toart-Terma Merr-y

---

A Johnson, L Bulgarelli, T Pollard, Horn. LA Celi, and R Mark. Mimic-iv, 2020.
R. E. Kalman, A New Approach to Linear Filtering and Prediction Problems. Journal of Basic Engineering, 82(1):35-445 03 1960. ISSN 0021-9223. doi: 10.1115/1.366255.. URL https
//doi. org/10 /0/.//..// 1115/1 3662552
John F. Kolen and Stefan C. Kremer. Gradient Flow in Recurrent Nets: The Difficulty of Learning LongTerm Dependencies, pages 237-243, Wiley-IEEE Press, 2001. doi: 00.1109/97847055440037 ch14.
Rahul G. Krishnan, Uri Shalit, and David Sontag, Deep kalman filters, 2015.
Seth Landefeld, E Francis Cook, Margaret Flatley, Monica Weisberg, and Lee Goldman. Identifica- tion and preliminary validation of predictors of major bleeding in hospitalized patients starting anticoagulant therapy. The American journal of medicine, 82(4):703-7133 1987.
Weihua Li, Sirish L.. Shah, and Deyun Xiao. Kalman filters in non-uniformly sampled multirate systems: For fdi and beyond. Automatica, 44(1):199-208, jan 2008. ISSN 0005- 1098. doi: 10. aut.://aaaammmmca000....... 05.000 009
URL http::///ooa..rrr.ooooooooaaaoa. wtt.://...utttttiii.ooooo 2007
Xuanqing Liu, Si Si, Qin Cao, Sanjiv Kumar, and Cho-Jui Hsieh. Neural sde: Stabilizing neural ode networks with stochastic noise. arXiv preprint arXiv: 1900002355,3 2019.
James Lu, Kaiwen Deng, Xinyuan Zhang, Gengbo Liu, and Yuanfang Guan. Neural-ode for pharmacokinetics modeling and its advantage alternative machine learning models in predicting new dosing regimens. Iscience, 24(7): 102804, 2021.
Amélie Marsot, Audrey Boulamery, Bernard Bruguerolle, and Nicolas Simon. Vancomycin. Clinical pharmacokineties, 51(1):1-13, 2012.
Thomas M. Moerland, Joost Broekens, and Catholijn M. Jonker Medde-baa learning: A survey. CoRR, fbss20006166122 2020. URL htttpss//arxxvv.or///a/////////t 16712
A. Moran and Peter Whittle. Hypothesis testing in time series analysis, 1951.
Shamim Nemati, Mohammad Ghassemi, and Gari Clifford. Optimal medication dosing from suboptimal clinical examples: A deep reinforcement learning approach. In 2016 38th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC), pages 2978 2981. IEEE, 2016.
Dominic A. Neu, Johannes Lahann, and Peter Fettke. A systematic literature review on state- of-the-art deep learning methods for process prediction. CoRR, abs/2101.09320,, 2021. URL
https://arxiv. www.:xiwww..oroaeooooooo htg//bb//201..99339
Se Yong Park and Anant Sahai. Intermittent kalman filtering: Eigenvalue cycles and nonuniform sampling. Proceedings of the 2011 American Control Conference, pages 3692-3697, 201 doi: 00.1009AACC.2011559999
Guy Revach, Nir Shlezinger, Xiaoyong Ni, Adria Lopez Escoriza, Ruud G. van Sloun, and Yonina C. Eldar. Kalmannet: Neural network aided kalman filtering for partially known dynamics, 2021.
Yulia Rubanova, Ricky Q. Chen, and David Duvenaud. Latent ordinary differential equations for irreuularly-sappled time series. In H. Wallach, H. Larochelle, A. Beygelzimer, F. d' d'Alché-Buu E. Fox, and R. Garnett, editors, Advances Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019. URL /ttettps://peee http://proeeeen...e/ee/// neur ips. wttpppppp...pppeeeeeeeep cc/pppprr/219///lee #26684555777770066c0004444144343
pdf
David Rumelhart et Learning representations by back-propagating errors. Nature, 1986. URL
https://wwit nature. com.aarticles/3225552333
12

---

Michael Rybak, Ben Lomaestro, John C. Rotschafer, Jr. Moellering, Robert, William Craig, Mariannn Billeter, Joseph R. Dalovisio, and Donald P. Levine. Therapeutic monitoring of vancomyen in adult patients: A consensus review of the American Society of Health-System Pharmacists, the Infectious Diseases Society of America, and the Society of Infectious Diseases Pharmacists. American Journal of Health-Systemm Pharmacy, 66(1):82-4) 01 2009. ISSN 1079-2082, doi: 10.21 4444..044444
Mona Schirmer, Mazin Eltaybb, Stefan Lessmann, and Maja Rudolph. Modeling inregular time series with continuous recurrntt units. In International Conference on Machine Learning, pages 19388-19405. PMLR, 2022.
B. Sinopoli, L. Schenato, M. Franceschetti, K. Poolla, M.I. Jordan, and S.S. Sastry. Kalman filtering with intermittent observations. JEEE Transactions Automatic Control, doi: 10.1 10.1109/TAAC.200488411..
2004.
Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N. Gomez, Lukasz Kaiser, and Illia Polosukhin. Attention all you need, 2017.
Yi-Jen Wang and Chin-Teng Lin. Runge-kutta neural network for identification dynamical systems in high accuracy. IEEE Transactions Neural Networks, 922)294-307, 1998.
Chao Yu, Jiming Liu, Shamim Nemati, and Guosheng Yin. Reinforcement learning in healthcare: A survey. ACM Computing Surveys (CSUR), 55(1):1-36, 2021.
13

---

Appendices
Table of Contents
A Observation Filtering: The Conditional Distribution and the Relation to Kalman Filtering
15
B Integrator Implementation
15
C The Dynamics Spectrum and Complex Figenfunction Implementation
16
Solver Analysis
Extended Experiments
17
E.1
Ablation Study for Patient Individualization
17
E.2
Comparison to ODE-based Methods
18
F,3
Sparse Observations
19
E.4
Synthetic Data Experiments with Regular Observations
20
E.5
Interpretability: Inspecting the Spectrum
21
Model Expressiveness and Overfitting
21
F
Medication Dosing Prediction: Implementation Details
22
F1
Data preprocessing
22
F22
LSTM Baseline Implementation
22
F3
Extended Results
23
14

---

A Observation Filtering: The Conditional Distribution and the Relation to Kalman Filtering
As described in Section 4 the NESDE algorithm keeps an estimated Normal distribution of the system state (t) any point time. The distribution develops continuously through time according the dynamics specified by Eq. ( except for the discrete times where an observation Y (1) received: in every such point of time, the (1) estimate is updated to be the conditional distribution (L)
The conditional distribution can be derived as follows. Recall that X (p, 2) (we remove the time index t as we focus now on filtering at single point of time). Denote X (Y, Z) where Y Rm (similarly to Eq. (5)) and Z R"-m;
Calculating the conditional Normal distribution:
and similarly, jr ((Y,#z) and
Eyy Eyz Ezy Ezz
Σ
First consider a noiseless observation (R = 0): then according to Faton 1983], the conditional
distribution X|Y = Y is given by X = -)
and Z ~ N(z M'zz). where
112 == pz Iz 2 -uy)
Ezzz Ezz
Ezy E Eyz
Exx
In the general case of R * O, we can redefine the state to include the observation explicitly: X = () = Y, Y, T)T, where is, E of X are adjusted by fy = My. Eyy = Eyy + R. Eer = R and Σ Σyz. Then, the conditional distribution can be derived as in the noiseless case above, by simply considering the new observation as a noiseless observation of X =
The relation to the Kalman filtering:
The derivation of the conditional distribution is equivalent to the filtering step of the Kalman filter [Kalman 1960 where the (discrete) model
X1+1 = A X2 + Ww
(wr N(0, Q))
Y H + 14
(M,,N(0, R)).
Our setup can be recovered by substituting the following observation model H E R=xn, which observes the first m coordinates of X and ignores the rest:
1
H
0
and the Kalman filtering step is then
K= EH (HEH + R1 n'
Hμ)
E KHE
Note that while the standard Kalman filter framework indeed supports the filtering of distributions upon arrival of new observation, its progress through time is limited to discrete and constant time-steps (see the model above), whereas our SDE-based model can directly make predictions to any arbitrary future time t.
B Integrator Implementation
Below, we describe the implementation of the integrator of the Eigen-SDE solver mentioned in Section 4
Numerical integration given u(t): In the presence of an arbitrary (continuous) control signal u(!).
it is impossible to compute the integral that corresponds with u(t) (Eq. (2) analytically. On the
15

---

other hand, u(t) is given in advance, and the eigenfunction, (t)), is a known function that can be calculated efficiently any given time, By discretizing the time to any fixed At, one could simply replace the integral by sum term
[ Φ(T) u(()dd
t + At)ulto + At)
while this sum represent I-to calculations, can be computed efficiently, as does not require any At recursive computation, both (() and u(() are pre-determined, known functions. Each element of the sum is independent of the other elements, and thus the computation could be parallelized.
Analytic integration: The control u is often constant over any single time-interval I (e.g., when the control is piecewise constant). In such cases, for given interval [to, t] in which u(t) = uy, the integral could be solved analytically:
[ (
J. ---Ar-- uIdT=
[" dTV dTV du (e- -M- - - to
one might notice that for large time intervals this form is numerically unstable, to address this issue, note that this integral is multiplied (Eq. ② by ((1) = VeA hence we stabilize the solution with the latter exponent:
(t)) (c-Aco - ~u = (CAl--to)-sse c - =V (cA((-to)- -
achieve a numerically stable computation.
In addition to the integral over ul(), we also need to calculate the integral over Q (Eq. (4 this case, Q is constant, and the following holds;
[" - - +
( ((()))(((7)-))))
e e
10
where denotes the Hadamard product, and
2X1
d1 +
An
A
An
A1
2An
In this form, possible solve the integral analytically, similarly to the integral of the control signal, and again, we use the exponent term from (tt) to obtain a numerically stable computation.
C The Dynamics Spectrum and Complex Eigenfunction Implementation
The form the eigenfunction matrix as presented in Section is valid for real eigenvalues. Complex eigenvalues induce slightly different form; firstly, they come in pairs, i.e., if
+ bi is an eigenvalue A (Eq. ( then こ = a bi (the complex conjugate 2) is an eigenvalue of A. The corresponding eigenvector of is complex well, denote it by Ureal + Vim i, then (the complex conjugate of v) is the eigenvector that correspond to z. Secondly, the eigenfunction matrix takes the form:
(1) = ear
|'reat cos(bt) Vim
sin(btt Vim cos(bt) + Ureat sin(bt)
For brevity, we consider only the elements that correspond with 2, Σ. To parametrize this form, we use the same number parameters (each complex number need two parameters represent, but since they come in pairs with their conjugates we get the same overall number) which are organized differently. Mixed eigenvalues (e.g., both real and complex) induce a mixed eigenfunction that is concatenation of the two forms. Since the complex case requires different computation, we leave the number of complex eigenvalues to be hyperparameter. Same as for the real eigenvalues setting, it is possible to derive an analytical computation for the integrals. Here, it takes different form, as
16

---

the complex cigenvalues introduce trigonometric functions to the eigenfunction matrix. To describe the analytical computation, first notice that:
cos( sin(bt)
sin cos(bt))
T (t) eat Ureal
Vim
1
and thus:
/cos(bt) sin(bl)
sin(bt))
$(1)
cos(bl)
Ureat
Uim
Note that here we consider a two-dimensional SDE, for the general case the trigonometric matrix is a block-diagonal matrix, and the exponent becomes diagonal matrix in which each element repeats twice. It is clear that similarly to the real eigenvalues case, the integral term that includes u (as shown above) can be decomposed, and it is possible to derive analytical solution for exponent multiplied by sine or cosine. One major difference is that here we use matrix product instead of Hadamard product. The integral over Q becomes more tricky, but can be separated and computed as well, with the assistance of basic linear algebra (both are implemented in our code).
D Solver Analysis
Below more complete version of Proposition Land its proof.
Proposition& 2 (Figen-SIDF solver optimality: complete formulation). Let (() be signal that follows Eq. (6) for any time interval I, = [&,, t.+1]. and u(t) a control signal that is constant over Ij for any i. For any consider the Eigen-SDE solver with the parameters corresponding to Eq. (6) (for the same Z,). Assume that the first solver (i = 0) is initialized with the true initial distribution (0) N(u, Eu), and for i 2 1, the i'th solver is initialized with the i 1'th output, along with an observation filter if an observation was received. For any interval and any time E I, consider the prediction (t) ~ N (p(t), E(t) of the solver. Then, ju() minimizes the expected square error of the signal (t), and (t) maximizes the expected log- ligeeiilldoodd of X(t).
Proof. We prove by induction over that for any i and any t E I,, (() corresponds to the true distribution of the signal X (t).
For i = 0, (t,) = (()) corresponds to the true initial distribution, and since there are no "interrupting" observations within To, then the solution Eqs. (2) and (4) of Eq. (6) corresponds to the
true distribution of X (t) for any [t,, 1++1)- Since constant over To, then the prediction (t) of the Eigen-SDE solver follows Eq. 4 accurately using the analytic integration (see Appendix [ note that if u were not constant, the solver would still follow the solution up a numeric integration error). Regarding f1, according to Appendix A (t,) corresponds to the true distribution of X (t1) after conditioning on the observation (t,) (if there was an observation at f otherwise, no filtering is needed). This completes the induction basis. Using the same arguments, if we assume for an arbitrary 2 0 that (t,)) corresponds to the true distribution, then x() corresponds to the true distribution for any E I, = [t,, tr+1), completing the induction.
Now, for any t, since (t) N(u(t), E{t)) is in fact the true distribution of Xt),, the expected square error E [()((t)] E [(u X())² is minimized by choosing μ1 p(t); and the expected log-likelihood E [e(()] E (log
Σ)] is maximized by = ((),, E = E()).
E
Extended Experiments
E.1
Ablation Study for Patient Individualization
To provide an insight over the importance of the dynamics-indde individualization, we perform an ablation study for the hypernetwork module. We use the same medical benchmarks as in Section 3 and fit a version of NESDE with neutralized hypernetwork module. In particular, we fix the context inputs of the module to be a vector of 1s, and thus prevent any propagation from the context features to the model's output. The results are presented in Table 4 and show great degradation in model
17

---

Table 4: Test mean square enrors (MSE) and negative bog-likelhhod (NLL) in medication--osing benchmarks.
This is an extension of Table 3 with the additional results of NESDE without the hypernetwork.
UH Dosing
Vancomycin Dosing
Model
MSE
NLL
MSE
NLL
Naive
613.3 + 13.48
112.2 + 16.4
LSTM
482.1 + 6.52
92.89 + 11.3
GRU-ODE-Bayes
491 + 6.88
4.52 + 0.008
80.54 + 11.8
CRU
4.49 0.012
76.4 12.8
8.87
NESDE no hypernet
529.7 + 13.34
5.42
87.32 11.57
11.57
3.73 0.13
NESDE (ours)
411.2 = 7.39
4.43 = 0.01
1
3.69 + 0.13
performance in the UH-dosing benchmark, approving that the hypernetwork indeed utilize the information within the context features. In the Vancomycin dosing benchmark, while still observe a degradation comparing to NESDE, the version of NESDE without hypernetwork still outperforms LSTM in terms of MSE and the rest of the baselines (except NESDE) in terms of NLL.
E.2
Comparison to ODE-based Methods
2019 a recent ODE- based method that can provide uncertainty estimation (which is typical requirement in medical applications). Similarly to other recent ODE-based methods [Chen et al. 2018), GRU-ODE-Bayes relies non-linear neural network model for the differential equation. GRU-ODE-Bay relatively poor prediction accuracy in Section [] which may be partially attributed to the benchmark settings. First, the benchmark required GRU-ODE-Bayer to handle control signal. As proposed in De Brouwer et al. [2019), we incorporated the control as part of the observation space. However, such control--bberrvatio mix raises time synchrony issues (c.g., most training input samples include only control signal without observation) and even affect the training supervision (since the new control dimension in the state space affects the loss). Second, as discussed above, the piecewise linear dynamics of NESDE provide higher sample efficiency in face of the 1000 training trajectories in Section
Section compares NESDE to GRU-ODE-Bayee De Brouwer et al.
presents
03
Figure 9: sample test trajectory of the sparsely-observable process. The observations and NESDE predictions (based on training over 400 trajectories) are presented separately for each of the two dimensions of the process.
In this section, explicitly study the sample efficiency of NESDE vs. GRU-ODE-Bayen in problem with no control signal. Specifically, we generate data from the GitHub repository De Brouwer et al.) [2019) The data consists of irregular samples of the two-dinensional Ornstein- Uhlenbeck process, which follows the SDE
dx; 0(H de)d + adWI,
where the noise follows Wiener process, which is set in this experiment to have the covariance matrix
1 0.5
Coc
0.5
The process is sparsely-observed: we use a sample rate of 0.6 (approximately 6 observations for 10 time units). Each sampled trajectory has time support of 10 time units. The process has two dimensions, and each observation can include either of the dimensions or both of them. The dynamics
18

---

of the process are linear and remain constant for all the trajectories: however, the stable "center" of the dynamics of each trajectory (similarly in Eq. (6)) is sampled from a uniform distribution, increasing the difficulty of the task and requiring to infer in an online manner.
Fig.P 9 presents sample of trajectory observations along with the corresponding predictions of the NESDE model (trained over 400 trajectories). Similarly to De Brouwer et al. [2019), the models are tested over each trajectory by observing all the measurements from times S. and then predicting the process at the times of the remaining observations until the end of the trajectory.
GRU_ODE
GRU_ODE
-02
NESDE
NESDE
-0.3
0.040
-0.4
0.35
-0.5
-0.5
0050
100
200
400
800
1600
100
200
400
800
1600
train _ize
train size
GRU_ODE
-0.4
-0.5
-0.5
GRU_ODE
0.025
-0.7
0.000
tmee
time
Figure 10: Top: losses of NESDE and GRU-ODE-Bayes over the OU benchmark. along with confidence intervals of 95% over the test trajectories. NESDE demonstrates higher data efficiency, its deterioration in small training datasets is moderate comparison GRU-ODE-Baees Bottom: errors vs. time, given 400 training trajectories, where all the test predictions rely on observations from times I 4. The advantage of NESDE becomes larger as the prediction horizon is longer.
To test for data efficiency, we train both models over training datasets with different numbers of trajectories. As shown in Fig. TO, the sparsely-observal setting with limited training data causes GRU-ODE-Bayes to falter, whereas NESDÈ learns robustly in this scenario. The advantage of NESDE over GRU-ODE-Bii -ayes increases when learning from smaller datasets (Fig. [0] top), when predicting for longer horizons (Fig. 10 bottom). This demonstrates the stability and data efficiency of the piccewise linear dynamics model of NESDE in comparison to non-linear ODE models.
E.3 Sparse Observations
This experiment addresses the sparsity of each trajectory. We use the same benchmark as Section and generate 4 train datasets, each one contains 400 trajectories, and test set of 1000 trajectories. In each train-set, the trajectories have the same number of data samples, which varies between datasets (4,6,8,10). The test-set contains trajectories of varying number of observations, over the same support. For each train-set, we train all the models until convergence, and test them. Fig. presents the MSE over the test set, for both the complex and the real eigenvalues settings. It is noticeable that even with very sparse observations, NESDE achieves good performance. Here, GRU-ODE -DDE-Bays appears to be more sample-efficiens than CRU and LSTM, but it is less sample efficient than NESDE.
19

---

0.6
NESDE
NESDE
GRUODE
GRUODE
MSE
LSTM
MSE
LSTM
CRU
CRU
0.3
0.2
10
Train data density [observations per trajectory]
Train data densily
[observations per trajectory]
(a) Complex dynamics
(b) Real dynamics
Figure 11: Test MSE vs. train observations eeservation--pe--tajeeoooooo per-trajectory. 95% confidence intervals are calculated over 1000
test trajectories.
E.4 Synthetic Data Experiments with Regular Observations
While NESDE (and ODE-based models) can provide predictions at any point of time, a vanilla LSTM limited the predefined prediction horizon. Shorter horizons provide higher temporal resolution, but this comes with a cost: more recursive computations are needed per time interval, increasing both learning complexity and running time. For example, if medical measurements are available once per hour while predictions are required every 10 seconds, the model would have to run recursively 360 times between consecutive measurements, and would have to be trained accordingly in advance. We use the synthetic data environment from Section 5 in the complex dynamics setting, and test both regularly and -istr--distribution control (see Section 5}. Here, we use L.STM models trained with resolutions of 1, 8 and 50 predictions per observation. All the LSTM models receive the control and the current observation Y as an input, along with boolean bo specifying observability: in absence of observation, LSTM layer, with 32 neurons between the two. To compare LSTMs with various resolutions, we work with regular samples, 10 samples, one at each second. The control changes in 10 seconds' resolution, and contains information about the true state.
set Y 0 and b, O. The models consist of a linear layer top of an
In Fig. 12c we present sample trajectory (without the control signal) with the predictions of the various LSTMs and NESDE. It can be observed that while NESDE provides continuous. smooth predictions, the resolution of the LSTMs must be adapted for good performance. As shown in Fig. 12a all the methods perform well from time 3 and on, still, NESDE and the low-resolution variants of LSTM attain the best results. The poor accuracy of the high-resolution LSTM demonstrates tradeoff in recursive models, moreover, GRUODE shows similar behavior in this analysis, which may hint on the recursive components within GRUODE.
the accuracy-vs-resoluee
NESDE: NSE-0.19
2.0
NESDE NSE-0.19
MSE=0043
signal
CRU: MSC-0,20
CRU: MSE-
LSTM 1:50
LSTH (1:1h KSE 18
1.0
LSTN [1:1]: NSE-O.21
LSTN
[ael:: NaEHO.22
prediction
0.4
LSTH 11:501: 35
LSTN (1:50): MSE-0.42
NESDE
0.4
observation
0.2
2
8
10
2
6
time
time
time
(a) One-step prediction
(b) Out distribution control
(c) Sample trajectory
Figure 12: MSE for predictions, relying the whole history of the trajectory for (a) the test set, and (b) cuttof-disttibution test set. The uncertainty corresponds to ..95-confidence-inteeee
over 1000 trajectories. (c) Sample trajectory and predictions. The LSTM predictions are limited to predefined times (c.g.. LSTM 1:1 only predicts at ohservation times). but their predictions are connected by lines for visibility. The shading corresponds to NESDE uncertainty (note that the LST does not provide uncertainty estimation).
The out-of-distribution test results (Fig. 12b) show that a change in the control policy could result in major errors; while NESDE achieves errors which are close to Fig. 12a the other methods deteriorate their performance. Notice the scale difference between the figures. The high-resolution LSTM and the bDE-based methods suffer the most, and the low-resolution variants of the LSTM, demonstrate
20

---

robustness to the control change. This result is similar to the results we present in Section 5 although here we see similarities between the variants of the LSTM and the ODE-based methods.
E.5 Interpretabillty; Inspecting the Spectrum
In addition explicit predictions at flexible times, NESDE provides direct estimation of the process dynamics, carrying significant information about the essence of the process.
For example, consider the following 3 processes, each with one observable variable and one latent vari-
able: A = -0.5 =i) with the corresponding eigenvalues λ ~ -0.75 + 1.98;; A₂ = ( (50.5 0.5
with A ~ (-1.3, -..1 and A3 = (2 =3) with A R #1.71.. As demonstrated in Fig. the three processes have substantially different dynamies: roughly speaking, real negative eigenvalues correspond to decay, whereas imaginary eigenvalues correspond to periodicity.
For each process, we train NESDE over dataset of 200 trajectories with 5-20 observations each. We set NESDE to assume an underlying dimension of n = 2 (i.e., one latent dimension in addition to the 1 observable variable); train it once in real mode (real eigenvalues) and once in complex mode (conjugate pairs of complex eigenvalues); and choose the model with the better NLL over the validation data. Note that instead of training twice, the required expressiveness could be obtained using n = 4 in complex mode (see Appendix C however, in this section we keep n = 2 for the sake of spectrum interpretability.
As the processes have linear dynamics, for each of them NESDE learned to predict a consistent dynamics model: all estimated eigenvalues are similar over different trajectories, with standard deviations smaller than 0.1. The learned eigenvalues for the three processes are \₃ -0.77 = 1.98;;
d (-0.7, -0.19)T= ;; and À -0.03 + 0.83i. That is, NESDE recovers the eigenvalues class (complex, real, or imaginary), which captures the essence of the dynamics - even though it only observes one of the two dimensions of the process. The eigenvalues not always recovered with high accuracy, possibly due to the latent dimensions making the dynamics formulation ambiguous.
E.6 Model Expressiveness and Overfitting
is well known that more complex models are capable to find complex connections within the data, but are also more likely to overfit the data. It is quite common that data that involves control is biased or affected by confounding factors: a pilot may change his course of flight because he saw storm that was off-the-radar; physician could adapt his treatment according to some measure that is off-charts. Usually, using enough validation data could solve the overfitting issue, although sometimes the same confounding effects show in the validation data, which results in model that is overfitted to the dataset. When targeting a model for control adjustment, is important that it would be robust to changes in the control; model that performs poorly when facing different control is unusable for control tuning. To exemplify an extreme case of confounding factors in the context of control, we add a correlation between the control (observed all times) the predictable measure (observed sparsely), in particular at times that the predictable is unobserved. We harness the same synthetic data benchmark in Section and use regular time samples, and the same LSTM baselines as in Appendix E.4 but here we generate different two types of control signals:
1. Same Distribution (SD): each time the control u(t) b 0.8 Y..
2. Out of Distribution (OOD): at each time I, the control u(() = b. + 0.8 Y..
by is random piecewise constant and Y, is the exact value of the measure we wish to predict. The first type used to generate the train and the test sets, additionally we generate an -iitt--distributtin test-set using the second type. We observe in Fig. 13 that GRU-ODE-Bun Bayes and the high-resolution L.TT achieve very low MSE over the SD as seen during training. CRU also achieves very low MSE, although not as much. The results over the OOD data show that the high performance over SD came with a cost - the better a model is over SD the worse it is over OOD. The results of LSTM 1:1 are not surprising, it sees the control signal only at observation-times so it cannot exploit the hidden information within the control signal. NESDE does not ignore such information, while maintaining the robustness w.r.t. control.
21

---

10
0.6
NESCE: MSE+0.17
1MPT:: MS--0 P8
CRLL-ODE: MISEOO
GRU ODE: MS5-5.56
0.5
CRUE MSE
0.4
USTM (3:33 NSS
LSTM (1:33: MSC-0
LSTH (1.8). MSE-2.67
2.3
1.5T (2 500: MS2-0.01
(1:50) NSE-2.75
0.2
0.1
10
(a) Same control distribution (b) Out of distribution control Figure 13: MSE for predictions under regular time samples, where the control signal is correlated to the Shows the results for set that has the
measure wish to predict, even times when it unobserved. same correlation between the control and the predictable measure the train set. (b) present the MSE for a different test set, with different correlation. Notice the different scales of the graphs.
F
Medication Dosing Prediction: Implementation Details
Below, we elaborate on the implementation details of Section 6
F1 Data preprocessing
We derive our data from the MIMIC-IV dataset [Johnson et al., 2020 available under the PhysioNet Credentialed Health Data License For the UH dosing dataset, extract the patients that were given UH during their intensive care unit (ICU) stay. We exclude patients that were treated with discrete (not continuous) doses of UH, or with other anticoagulants; or that were tested for aPTT less than two times. The control signal (UH dosing rate) is normalized by the patient weight. Each trajectory of measurements set to begin one hour before the first UH dose, and is split in the case of 48 hours without UH admission. This process resulted with 5866 trajectories, containing a continuous UH signal, an guregularly-observe aPTT signal, and discretized context features. Note that we do not normalize the aPTT values.
Heparin:
Vancomycin:
The VM dosing dataset derived similarly, from patients who received VM during their ICU stay, where we consider only patients with at least 2 VM concentration measurements. Fach trajectory begins the patient's 's admission time, and also split in the case of 48 hours without VM dosage. Additionally, we add an artificial observation of 0 at time t = 0, as the VM concentration is 0 before any dose was given (we do not use these observations when computing the error).
General implementation details: For each train trajectory, we only sample some the observations, to enforce longer and different prediction horizons, which was found to aid the training robustness. Hyperparameters (c.g., learning rate) were chosen by tria--and-eroor with respect to the validation-sed (separately for each model).
Context variables C' are used in both domains. We extract 42 features, some measured continuously (e.g., heart rate, blood pressure), some discrete (e.g., lab tests, weight) and some static (e.g., age, background diagnoses). Each feature is averaged (after removing outliers) over a fixed time-interval of four hours, and then normalized.
F.2
LSTM Baseline Implementation
The LSTM module we use as a baseline has been tailored specifically to the setting:
1. It includes an embedding unit for the context, which is updated whenever a context is observed, and an embedded context is stored for future use.
2. The inputs for the module include the embedded context, the previous observations, the control signal and the time difference between the current time and the next prediction time.
3. Where the control signal is piecewise constant: any time it changes we produce predictions (even though no sample observed) that are then used as an input for the model, to model the effect of the UH more accurately.
22

---

We train with the same methodology we use for NESDE where the training hyperparameters chosen by the best performance over the validation data.
The model contains two fully connected elements: one for the context, with hidden layers of size 32 and 16-dimensional output which is fed into a I anh activation; the second one uses the LSTM output to produce a one-dimensional output, which is fed into ReLU activation to produce positive outputs, its size determined by the LSTM dimensions. The LSTM itself has an input of 19 dimensions; 16 + 1 + 1 + 1 for the context, control, previous observations and the time interval predict. It has hidden size 64 and two recurrent layers, with dropout of 0.2. All the interconnections between the linear layers include ReLU activations.
Architecture for the medication dosing benchmarks:
Architecture for the synthetic data benchmarks: Here, there is no context, then the model contains one fully connected element that receives the LSTM output and has linear layers of sizes 32 and with a Tanh activation between them. The LSTM has an input of dimensions; for the state, control signal, and the time interval predict. has hidden size of 32 and two recurrent layers, with dropout of 0.2.
F.3
Extended Results
The figures below present more detailed information for the experiments discussed in Section 6
Samples in the data
600
400
1250 1500
Time since last observation [m]
Figure 14: Histogram prediction horizons the UH dosing data (Section 6). Notice that the peak of the
histogram around 6 hours (360 minutes) corresponds to the accuracy peak of the LSTM and GRU-ODE-Bayes
in Fig. g
23


Reviews:
[{"review_id": "pVQcZtkEoW", "reviewer": "Reviewer_uGSy", "summary": "This paper presents a method for predicting time series at irregular time points. Its approach systematically integrates control signals, a convenient attribute for modelling dosing regimes in medicine. The latent process follows a piecewise linear stochastic differential equation (SDE) with control. The SDE is parameterized using its spectral representation, allowing for effective and analytical propagation. The parameters of the SDE are learned by a Hypernet. Experimental validation of the method is conducted on artificial data and two medication dosing problems.", "questions": "--", "limitations": "The piece-wise linear assumption on the observation sequence $Y(t)$ seems a restrictive assumption and limits the type of dynamics that can be modeled with the method. (From my understanding, the model learns a piece-wise linear SDE for $X(t)$ where the observations $Y(t)$ is simply a subset of the features of $X(t)$ , i.e. $Y(t)=X(t)_{1:m}$  . Therefore, the sole non-linearity allowed in the observation sequence come from the piece-wise linear formulation.) A discussion on when the piece-wise linear assumption fails compared to other (more expressive) models is encouraged. ", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "- The paper introduces a method that effectively models control, with experimental results indicating an improvement over baseline models lacking explicit control.\n- The well-structured experimental study underscores the applicability and value of the method for two specific medical dosing scenarios.\n- The authors propose an efficient technique to compute dynamics by parameterizing the system's eigenfunction, thereby enabling fast and analytical computation.\n- This parameterization strategy not only facilitates the integration of prior knowledge on the dynamics, as evidenced in the dosing experiments, but also enhances interpretability—which is a useful attribute in practise", "weaknesses": "- The related work sections seems comprehensive, though neural SDE approaches could also be discussed as the model is based on a SDE latent process.\n- Explicitly stating which equations represent ESDE() and Filter() in Algorithm 1 would enhance clarity.\n- Although the experimental study appears to be well-executed and thoughtfully designed, the paper could be enhanced by including more baseline comparisons. Specifically, a comparison with NeuralSDE models would be intriguing, as they also yield uncertainty estimates and aren't confined to the piece-wise linearity constraint.\n- Minor: My understanding is that the context does not change over time, though Figure 2 suggest this on the first sight.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their helpful feedback. We will modify the manuscript as discussed below to address the reviewer’s comments. In particular, we will clarify the selection of baselines, which seems to be the main concern in the review. If our responses are satisfactory, we will appreciate it if the reviewer can re-evaluate our contribution accordingly.\n\n**Related Work and Baselines**\n\nAs the reviewer suggests, our main baselines are indeed from the family of Neural-ODE. Specifically, we looked for Neural-ODE baselines that both handle noise and provide uncertainty estimation (so that NLL can be measured). We selected GRU-ODE-Bayes and CRU, which satisfy both properties, and have available code on GitHub. Note that despite the different ODE/SDE terminology, both methods practically implement a Neural-SDE, which must not be confused with [1] that uses noise as regularization for Neural-ODEs (without uncertainty estimates). We will clarify this in the manuscript. Thanks for pointing to it.\n\n**Explicitly stating which equations represent ESDE() and Filter() in Algorithm 1 would enhance clarity.**\n\n We will state the equation for ESDE()  (Eq. 4). For the Filter(), we will refer to Appendix A.\n\n**Minor: My understanding is that the context does not change over time, though Figure 2 suggest this on the first sight.**\n\nThe context changes over time, although in fixed time-intervals. We will add a statement that would address this unclarity.\n\n__________\n\nRegarding the limitations of the piecewise linear assumption:\nThe piecewise linear assumption highly depends on the interval length (Algorithm 1, and lines 207-216), which is a hyperparameter that controls how often the (linear) dynamics change to another (linear) dynamics. As it is a hyperparameter, it could be adjusted to the specific domain -- small intervals for “highly nonlinear” domains and large intervals for domains which behave close to linear. The additional expressiveness comes at the cost of more computations, and recursive updates, which could harm the overall performance, as we see in the nonlinear models we compare to, hence the importance of adjusting it.\n\n\n__________\n[1] Liu, Xuanqing, et al. \"Neural sde: Stabilizing neural ode networks with stochastic noise.\" arXiv preprint arXiv:1906.02355 (2019)."}, {"review_id": "5qB1nMPnCZ", "reviewer": "Reviewer_vmHv", "summary": "This paper proposes a novel network approach, named Neural Eigen Stochastic Differential Equations (NESDE), to solve sequential prediction problems mainly in medical dosing control field of applications. ", "questions": "$\\mathbf{1}.$ See $\\mathbf{1}.$ in $\\mathbf{Weakness}$.\n\n$\\mathbf{2}.$ The experiment with synthetic data seems to support, e.g., the sample efficiency of the proposed method. However, it lacks training efficiency/computational complexity in the experiments. Is it possible to quantify the computational cost when comparing with other methods?   \n\n$\\mathbf{Trivial}.$ What is the point of comparison in Figure 1?  ", "limitations": "This paper does not provide broader impact or limitation statements, but i do not have any concerns on potential negative societal impact.", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "$\\mathbf{1}.$ The proposal of NESDE is novel, along with the hypernet that determines the subsequent model parameters.\n\n$\\mathbf{2}.$ This paper is very well written. The motivation, i.e., to address challenges in sequential prediction of medical dosing control, and the limitation of current works, is very well explained. The related works and background fundamental theories are also very well stated with intuitive illustrations. \n\n$\\mathbf{3}.$ The theoretical analysis is comprehensive and sound (also seen in $\\mathbf{Appendix}$.) \n\n$\\mathbf{4}.$ The empirical analysis is detailed and thorough. Author designed the experiments to validity the algorithmic merits and also real life application in the medical dosing field. ", "weaknesses": "$\\mathbf{Note: }.$ It is a unfamiliar application field to me, and thus it is hard to determine whether or not the proposed method is indeed a significantly novel approach in the field. \n\n$\\mathbf{1}.$ The major concern is whether or not the proposed method has a broader application. It seems that the direct counterpart is neural ordinary differential equations, and is it possible to discuss and compare the general applicabilities of NODEs versus NESDEs?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their positive and detailed feedback.\n\n**General applicability of NESDE and Neural-ODE**\n\nNESDE constitutes a framework for forecasting individualized dosing dynamics. As described in Section 3, we consider the general problem of forecasting the dynamics of a controlled SDE. This problem is in fact a wide family of problems, including medication dosing regimes, and general controlled chemical or biological processes. Many such processes require control: fermentation in the wine-making process, irrigating and fertilizing a crop (agriculture), and monitoring chemical reactions under uncertainty. Neural-ODE on the other hand could be used for various applications (e.g., classification, regression, forecasting), when combined in different architectures. \nIn the current version of our paper, in Sections 1, 1.1, we discuss the applicability of NODE-based methods to real-life problems, and in particular to the problem of forecasting dosing dynamics. In Section 5 we compare NESDE to GRU-ODE-Bayes, which is a NODE-based method for forecasting the dynamics of an SDE.\n\nFollowing your comment, we will add a short discussion, comparing the applicability of both methods, in Section 1.1. Thanks for bringing this to our attention.\n\n**The experiment with synthetic data seems to support, e.g., the sample efficiency of the proposed method. However, it lacks training efficiency/computational complexity in the experiments. Is it possible to quantify the computational cost when comparing with other methods?**\n\nThe computational complexity could be measured in run-time, and we will add it to the appendix. Note that compared to the neural ODE methods, NESDE enjoys better run-time.\n\n**What is the point of comparison in Figure 1?**\n\nThe purpose of Figure 1 is to illustrate the difficulty of the combination of sparse observations with noisy dynamics. We present two types of SDEs, not comparing them. We will clear this out in the figure's caption."}, {"review_id": "fMmj0a4rem", "reviewer": "Reviewer_hUXX", "summary": "The paper presents a piece-wise linear neural SDE with control for biomedical forecasting. The key contribution is combining rigorous linear SDE solutions with hypernet parameterisation.", "questions": "The imaginary/cyclicity argument could have been elaborated more in experiments. Is medical data cyclical? I'm not sure. \n\nHow come in fig3 we get same performance with 2 or 8 datapoints? ", "limitations": "No issues", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "The paper is fantastically written and very clear. \n\nThe idea of linearising SDEs is relatively novel, and the overall model is sensible and principled.\n\nThe results are good, and show that a relatively simple model construction is perhaps ideal for medical data. ", "weaknesses": "The method is somewhat incremental: it combines linear SDEs with hypernets.\n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their positive feedback.\n\n**Regarding the value of the method**\n\nNESDE aims to solve a hard problem – forecasting and filtering a sparsely-observed SDE with an independent control signal. To address the accompanying challenges, NESDE combines several components. First, in this context, models that can capture complex behaviors are prone to overfitting (see the OOD experiments, Section 5). NESDE enables to inject domain knowledge and control the model’s complexity/timescale by tuning its hyperparameters; size of the latent process, type (imaginary/real/complex) of the eigenvalues, stability of the system, partition to context/state/control, and the interval length. By tuning these hyperparameters, NESDE can capture complex behavior and provide forecasts which are more robust to OOD control. Second, the ability to update future forecasts, given the current observation, is important. The filtering mechanism of NESDE adapts to new observations online, which translates to better predictions. Finally, by using the analytical solution, NESDE does not require the usage of an integrator, which is computationally expensive, and mandatory for Neural-ODE based methods.\n\n\n**The imaginary/cyclicity argument could have been elaborated more in experiments. Is medical data cyclical? I'm not sure.**\n\nThe imaginary argument is a general capability of NESDE, but we did not find it useful in the dosing problems we addressed. We do include experiments with complex eigenvalues (which are cyclic and decaying) in Section 5 and in the appendix. In appendix section E.1 we show how NESDE handles pure imaginary eigenvalues.\n\n**How come in fig3 we get same performance with 2 or 8 datapoints?**\n\nThe number of datapoints correspond to the filtration process: the model learns the dynamics from “offline” data, then infers “online” data in real-time. Here, there is no observation noise, then at some point (after approx’ 2 observations) the uncertainty of NESDE becomes solely the noise of the dynamics (with an additional small error) and thus cannot be greatly improved."}, {"review_id": "irBkslqyu3", "reviewer": "Reviewer_a6Bv", "summary": "This develops a robust SDE algorithm for medical forecasting. This is demonstrated on both synthetic and real data.", "questions": "SDEs are not my area of expertise. I expect the other reviewers will have more enlightening questions and suggestions.", "limitations": "No societal ramifications, but no limitations mentioned as far as shortcomings of the model are concerned.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The problem solved is clearly compelling and the performance is evaluated on a wide range of datasets, where it exhibits superior performance to reasonable alternatives. It would appear to be an overall robust method.", "weaknesses": "It would be nice if some of the components of the model such as hypernet were defined here rather than briefly described, even though I understand the desire to save space. It makes it so that this paper is more difficult to understand/review.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their comments. If our response is satisfactory, we kindly ask the reviewer to update the score. \n\n**Components Description**\n\nIn the current version, we briefly describe the components in Section 4. The observation filter is further elaborated in Appendix A, the ESDE is elaborated in Appendices B, C.\nFor the Hypernet, we will add to the appendices a complete presentation of the architecture, in addition to a block diagram (attached as a pdf file in the author rebuttal) that specifically explains how the Hypernet module functions within NESDE. Thanks for pointing it out."}]

Author rebuttal:
We are excited that all the reviewers recommended accepting the paper, and thank them for their helpful comments. The reviewers found that the “problem solved is clearly compelling” (a6Bv), and the paper is “fantastically written and very clear” (hUXX). They describe our proposed method by “the proposal of NESDE is novel” (vmHv), “effectively models control” (uGSy), and “also enhances interpretability—which is a useful attribute in practise” (uGSy). The empirical evaluation “is detailed and thorough” (vmHv) and “underscores the applicability and value of the method” (uGSy). The results “show that a relatively simple model construction is perhaps ideal for medical data” (hUXX).

Meta-review: Four reviewers evaluated the submission and scored it 7666. The agreed strengths are

* The approach (NESDE and the hypernet) is novel and interpretable
* Experimental study is carefully designed and executed, and shows good performance in two specific medical dosing scenarios

Based on the positive feedback, I recommend acceptance.
Decision: Accept (poster)