Integrating Expert ODEs into Neural ODEs: Pharmacology and Disease Progression
Authors: Zhaozhi Qian
William R. Zame
Lucas M Fleuren
Paul Elbers
Mihaela van der Schaar
Venue: NeurIPS 2021 Poster (2021)
Keywords: COVID-19
Neural ODE
Hybrid Model
Expert-designed ODE

Abstract: Modeling a system's temporal behaviour in reaction to external stimuli is a fundamental problem in many areas. Pure Machine Learning (ML) approaches often fail in the small sample regime and cannot provide actionable insights beyond predictions. A promising modification has been to incorporate expert domain knowledge into ML models. The application we consider is predicting the patient health status and disease progression over time, where a wealth of domain knowledge is available from pharmacology. Pharmacological models describe the dynamics of carefully-chosen medically meaningful variables in terms of systems of Ordinary Differential Equations (ODEs). However, these models only describe a limited collection of variables, and these variables are often not observable in clinical environments. To close this gap, we propose the latent hybridisation model (LHM) that integrates a system of expert-designed ODEs with machine-learned Neural ODEs to fully describe the dynamics of the system and to link the expert and latent variables to observable quantities. We evaluated LHM on synthetic data as well as real-world intensive care data of COVID-19 patients. LHM consistently outperforms previous works, especially when few training samples are available such as at the beginning of the pandemic.

Integrating Expert ODEs into Neural ODEs:
Pharmacology and Disease Progression
Zhaozhi Qian University of Cambridge zq2240cam.ae ac.
William R. Zame UCLA zame0econ. ucla. edu
Lucas M. Fleuren
Amsterdam UMC
w.Flluren@anstteddamucc
Paul Elbers Amsterdam UMC
Mihaela van der Schaar University of Cambridge
p elberr@@amsterdammnc nl
UCLA
The Alan Turing Institute mv4720cam ac. uk
Abstract
Modeling system's temporal behaviour in reaction to external stimuli fun- damental problem in many areas. Pure Machine Learning (ML) approaches often fail in the small sample regime and cannot provide actionable insights beyond predictions. A promising modification has been to incorporate expert domain knowledge into ML models. The application we consider predicting the pa- tient health status and disease progression over time, where wealth of domain knowledge is available from pharmacology. Pharmacological models describe the dynamics of carefully-chosen medically meaningful variables in terms of systems Ordinary Differential Equations (ODEs). However, these models only describe limited collection of variables, and these variables often not observable in clinical environments. To close this gap, we propose the latent hybridisation model (LHM) that integrates system of expert-designed ODEs with machine-learner Neural ODEs to fully describe the dynamics of the system and to link the expert and latent variables to observable quantities. We evaluated LHM synthetic data as well as real-world intensive care data of COVID-19 patients. LHM consistently outperforms previous works, especially when few training samples available such as at the beginning of the pandemic.
1 Introduction
Understanding the temporal evolution of a dynamical system is the central problem in many areas. The Machine Learning (ML) approach to this problem has been to learn collection of latent variables and construct a dynamical model of the system directly from observational data. While ML has achieved strong predictive performance in some applications, it has central weaknesses. The first is that it requires large datasets. The second is that the latent variables that the ML approach identifies often have no physical interpretation and do not correspond to any previously-ident -identiiiiddd quantities.
One approach to dealing with these weaknesses has been to incorporate expert domain knowledge into ML models. Most of the work using this approach has focused incorporating high-level knowledge about the underlying physical system, such as conservation of energy [9, 32, 99], independence of mechanism [63], monotonicity [59], or linearity [34]. In addition, there have been attempts to integrate domain-specific. "expert models" into ML models to create "hybrid" models. Most of this
work has employed expert models that directly issue predictions [54, 86, 91, 93] extract useful
features from the raw measurements [41].
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

(A) Pharmacology
(B) Machine learning model
(C) Hybrid model
a
Z
a
a
?
?
X
X
laboratory setting
clinical setting
Figure I: Dependency structure of the three models designed for the laboratory clinical settings. Dashed nodes represent unobservable variables. The expert variables are observable the laboratory setting but not in the clinical setting. The pharmacological model dless not contain the links to the clinical variables x.
The approach taken in this paper begins with expert model in the form of system of Ordinary Differential Equations (ODEs) and integrates that expert model into a system of Neural ODEs [14]. The specific problem address is that of predicting disease progression and health status over time; the specific expert model(s) come from Pharmacology [42] - but we believe our approach may be much more widely applicable. For a number of diseases, available pharmacological models, built the basis of specialized knowledge and laboratory experiments, provide a description of the dynamics of carefully-chosen medically meaningful variables in terms of ODEs that govern the evolution of these states [20, 2, 30]. However, these models are typically not directly applicable in clinical environments, because they involve too few variables fully describe a patient's health state [80, 38], because the expert variables which the models employ may be observable in the laboratory setting but clinical environments [29, 6]. and because the relationships between the expert variables and clinically observable quantities is not known [25]. We will give a example later in Section 3.5.
This paper proposes a novel hybrid modeling framework, the Latent Hybridisation Model (LHM), that imbeds a given pharmacological model (a collection of expert variables and the ODEs that describe the evolution of these variables) into larger latent variable MI. model (a system of Neural ODEs). In the larger model, we use observational data to learn both the evolution of the unobservable latent variables and the relationship between measurements and all the latent variables the expert variables from the tharmacologica model and the latent variables in the larger model. The machine learning component provides links between the expert variables and the clinical measurements, the underlying pharmacological model improves sample efficiency, and the expert variables provide additional insights to the clinicians, A variety of experiments (using synthetic and real data) demonstrate the effectiveness of our hybrid approach.
2 Problem setting
We consider a set of hospitalized patients [N] = (1...,,NT N} over a time horizon [0. 7); t 0 represents the time of admission and t T represents the maximal length of stay. The health status of each patient i is characterized by collection of observable physiological variables x;()) E RR D E N+; because the physiological variables may include vital signs, bloodwork values, biomarkers, etc., (I) is typically high-dimensional vector. Although the physiological variables are observable, they are typically measured only at discrete times, and with error. To avoid confusion, we distinguish the measurements of these variables yt) from the true values; i.c.
(1) xt) + 6ic
where the independent noise term EØ accommodates the measurement error (modeling e, autocorrelated stochastic process is left as a future work). For illustrative purposes, we also assume that e, follows Normal distribution N (0, of), but any parametric distribution could be easily accommodated. We denote the measurement times for each patient as T, = {t.1, f........ We write (t) R A N+ for the treatments the patient receives. Some treatments (e.g. intravenous medications) are continuous; others (e.g. surgical interventions) are discrete, so some components of (1) may be continuous functions but others are (discontinuous) step functions.
convenient to write A1 12] (t)(1) L 12) and [1 (2] (yx(1) S15 12, T.) } for the treatments and measurements (respectively) during the the time window [&1, t2]. Note that

---

Do f] and A. represent histories at time t while A, T) and Di(t: T) represent treatment plans and predictions, respectively. Our objective to predict the future measurements given the history and treatment plan A: A:llo T): [P 7) v. A.). tol, A. :T)
Historical observations
Treatment plan
Understanding this distribution will allow us to compute both point estimates and credible intervals (reflecting uncertainty). Note that is important for the clinician to understand uncertainty order to balance risk and reward. When the context is clear, we will omit the subscript i and the time index L.
3 Method
3.1 The pharmacological model
We begin with a pharmacological model which describes the dynamics of a collection of "expert" variables z" (4) ER Each expert variable captures a distinct and medically-meaningg aspect of the human body, e.g. the activation of immune system. The pharmacological model describes the dynamics as a system of Ordinary Differential Equations (ODEs):
こ (t) (z² (t), O"));
where we have written (t) for the time derivative zf. The functional form of f IR XRA is specified but the unknown parameters 8 (e.g., coefficients) need to be estimated from data." 1.2
R
is important to note that the system of ODEs (3) describes dynamies that are self-contained, in the sense that the time derivatives 2 (L) depend only on the current values of the expert variables z (1) and the current treatments at)) and not on histories or on other variables. To ensure that this obtains, it may be necessary to limit the scope of the model and limit attention to single system of the body (or perhaps to several closely related systems) [19]. As consequence of these limitations, the expert variables will usually not give full picture of the health status of the patient and will usually not account for the full array of observable physiological variables x(f) [29].
3.2 The latent hybridisation model: linking expert variables with measurements
As we have already noted, pharmacological models are typically developed and calibrated in the laboratory, where the expert variables can be directly measured - in patients, in laboratory animals, or even in vitro (Figure A). In clinical environments, the expert variables are frequently not observed (Figure B and C). To use the pharmacological models in clinical environments, we must establish links between the expert variables z (t) and the clinical measurements yt)). To do this we introduce additional latent variables (t) E M and posit the following relationship between the latent variables z", z and the observable physiological variables x():
x (1) = g(z" (1), 2" (1), a()); 7)
The function g REXMXA + R is neural network with (unknown) weights y, and maps the latent space to the "physiological space". Here we also allow the treatment a(t) to affect the mapping between the latent zf), z" (t) and the observable x(t). This is a standard design in the state space modeling literature, which we conform to (e.g. Equation in [8] and Equation and 2 in [85]). We posit that the dynamics of the latent variables z" (t) follow system of ODEs governed by its current values zn (t), the treatments a(t) and the current values of the expert variables (t).
2mhh (1) = 1T (z"
(2" (z""()) (t),
xf (1), a(t;; 0m)),
The function fm :R - R is a neural network with (unknown) weights 6th Equations (3)-(5) specify the dynamics of LHM. It is convenient to write z(t) = (t) z z²'(tt (t)) for the vector of
all latent variables and = (0", 8m 7, o) for the set of all (unknown) coefficients.
'The system in Equation (3) is quite general; appropriate choices of expert variables allow it to capture both high-order ODEs and dime-dependent ODEs [67].
"Some care must be taken because systems such as (3) do not always admit unique global solutions. In practice, the pharmacological models are sufficiently well behaved that global solutions exist and are unique. Although closed-form solutions may not be available, there various efficient numerical methods for solution.
3

---

2
Numerically solve ODE f", f
Sample
8
Encoder
Latent space
Measurements y/o : to]
3
Apply mapping g
0/0 T
History
"Reconstruction" to
Prediction T
Figure 2: Illustration of the training and prediction procedure.
The coefficients θ will be learned from data. However, even after these coefficients learned, the initial state of the patient z: (0) is still unknown. In fact, the variation in initial states reflects the heterogeneity of the patient population. the coefficients and the initial state were known, the entire trajectory of z; given the treatments could be computed (numerically). Because we have assumed that the noise/errors e, are independent, we have
Mto T] TT D. : 2 (,) T]. O, Vto <T
(6)
However, because the initial state is unknown, it must be learned from the measurements (t).
LHM would reduce to pure latent neural ODE model [14, 72] if we omitted the expert variables (Figure B). However, that would amount to discarding prior (expert) information and so is evidently undesirable. Indeed, as have noted in the introduction, incorporating this prior (expert) information into our hybrid model.
approach is driven by the idea of
In the current work, assume that the pharmacological model in Equation 3 is correct. In practice, the model might wrong two ways. The obvious way is that the functional form of f might be misspecified (e.g. a linear model might be specified when the truth is actually nonlinear). Many existing techniques can address such misspecification and could integrated into LHM [35, 64, 98]; see the discussion in Appendix A.7. Alternatively, it might be that the system of expert variables is not self-contained, and that their evolution actually depends on additional latent variables, we leave this more challenging problem for future work.
Practical extensions to LHM such as including static covariates and modeling informative sampling are discussed in Appendix A.7.
3.3 Independent and informative priors
It may be challenging pinpoint the exact value of the latent variables z based observations (e.g. due to measurement noise or sampling irregularity). For this reason, we quantify the uncertainty around z using Bayesian inference. In what follows, assume the initial states 2;( of patients are independently sampled from prior distribution (, ~ Po. Two points are worth noting.
Independent Priors. We use independent prior distributions on the expert variables x and the latent variables zzz, i.c. (()( = (0)) x PPz" (0)). This guarantees that information in zz (0) does not duplicate (any of the) information in z" (0), which captures our belief that the latent variables are incremental to the expert variables. addition, independent priors are also commonly used in Bayesian latent variable models such as variational autoencoders (VAEs) [47, 37J.
Informative Priors The prior distribution on the expert variables (z" (z)) should reflect domain
knowledge. Such knowledge is usually available from previous studies in Pharmacology [38]. Using an informative prior tends to improve the estimation settings [52]. Moreover, the expert variables usually take values in specific ranges (e.g. [0, 10] [42])
latent variables, especially in small-sample
and going beyond the valid range may lead to divergence. The informative prior can encode such prior knowledge to stabilize training.

---

3,4 Model training and prediction via amortized variational inference
Given the training dataset D {():[o A; TEN we use amortized variational inference
(AVT) to estimate the global parameters @ and the unknown initial condition 2,(0) [97]. Figure 2 presents a diagram of the training procedure. We start by learning a variational distribution to approximate the posterior (0) D.[0= 77, A, T]). As is standard in AVI [47, 97], we use
Normal distribution with diagonal covariance matrix as approximation:
(0, T), 7) = ((4, Ej); 14, Σ, =
A, T) ②)
Here the parameters μ4, Σ; are produced by an inference network (also known an encoder) e(-) with trainable weights ø. When the context clear, we will denote the variational distribution defined by Equations (7) as Qa- The evidence lower bound (ELBO) for the global parameter Q and the inference network parameters Ø is defined as
ELBO(O, Φ) = E2 [logP()][0: ( 77,44(0:),2000 (),) - P
(8)
To compute the ELBO for given Q and Φ. we sample z(0 Qa and numerically solve the ODEs to obtain (t). Vt E [0, T] (Figure 2: Steps 1, 2). Then, we compute the inner log-likelihood function using the mapping g and the noise distribution (Equations (1) and (4)). Finally, we use Monte Carlo sampling to evaluate the KL divergence term: Ex [log (0) [ logIPo (z))) This is because the informative prior Po may not have an analytical KL divergence (unlike the standard Normal prior used in previous works [47, 72). We optimize ELBO by stochastic gradient ascent and update all parameters jointly in an end-to-end manner (detailed in Appendix A.3).
The prediction procedure follows the same steps as illustrated in Figure 2. For a new patient with history D.[0= to] A, tol. we first estimate the variational posterior Qo using the trained encoder.
From Equation (6), we can estimate the target distribution in Equation (2) as:
E2 (I)))(()
to), A, T])], Vt to.
where Attty T| is a future treatment plan. The outer expectation can be approximated by Monte Carlo sampling from Qo and the inner probability is given by the likelihood function.
Choice of variational distribution and encoder. The training procedure above is agmostic to the exact choice of variational distribution and encoder architecture. We choose the Normal distribution to make fair comparisons with the previous works [14, 72]. . For the same reason, we use the reversed time-aware LSTM encoder proposed in [14]. In the Appendix A.4, we show additional experiments with more complex variational distributions, i.c. Normalizing Flows [70].
3.5
Using LHM to provide clinical decision support
In order for clinicians to optimally treat patients, they need to predict the progression of disease given the treatments. Although machine learning models may demonstrate feature importance [16, 4], they do not uncover the relationships between those features and the underlying pathophysiology. LHM can provide the missing link between clinical observations and disease mechanisms. In combination with clinical reasoning, it can provide treating clinicians with decision support in several complementary ways.
First, LHM can inform the clinicians about the values of the expert variables (t) that cannot be observed in the clinical environment but important for prognosis, choice of treatment, and anticipation of complications. For example, understanding and predicting immune response is pivotal when deciding on immunosuppresive therapy in the treatment of COVID-19: an extreme immune response may lead to potentially fatal cytokine storm [24], but a suppressed immune response may be equally dangerous in case of (secondary) infection [48. 84J. However, because immune response is not directly observable in the clinical environment, clinicians must rely on proxies such as C-reactive protein (CRP) for inflammation [58]; by their very nature, such proxies are noisy and highly imperfect measures of the desired values.
Secondly, LHM can provide the clinician with predictions of the disease progression given the treatments, enabling the clinicians to design the best treatment plan for the patient at hand.
Finally, LHM can bridge the gap between the laboratory and clinical environments, helping to align model output with clinical reasoning, and thus to bring models to the patient bedside and also to foster translational research [78, 28].
5

---

4 Related works
Hybrid models. Hybrid models combine a given expert model with ML. [88]. Depending on the type and functionality of the expert model, various approaches have been proposed. Residual Models and Ensembley use expert models that can issue predictions directly [54, 86, 89, 91, 93]. A Residual Model fits a ML model to the residuals of the expert model while an Ensemble averages the ML and expert predictions. Feature Extraction makes use of an expert model that extracts useful features from the measurements [41] J; an ML model then uses these features to make predictions. These methods are not suitable for our setting because our expert model is an ODE that governs the latent variables (5); it does not issue predictions of measurements nor does it extract features that the MI. model can use. Appendix A.5 Table 2 summarizes these approaches.
ML inspired by physics uses physical laws to guide the design of architectures [77, 96, 82], loss functions [94, 26], and weight initialization [68]. Examples include Hamiltonian neural networks [32, 9, 99], which reflect the conservation of energy. These models utilize general physical laws rather than a specific expert model, and are rather different than the hybrid models discussed above.
Works involving expert ODEs. Some recent works also propose to integrate expert ODEs into the learning system. It is worth comparing with their problem settings and proposed methods. [53] assumes that the observations are directly generated by the expert variables x(t) = g(z" (t)) (compare with Equation 4 of this work). It does not involve any machine-learned latent variable or neural ODE, Hence, [53] tends to fall short when the expert variables alone are inadequate to predict the observations. The method in [56] differs from LHM in the type of expert ODE used and the way to incorporate machine-learned latent variables. In [56]. the expert ODE depends on an unknown time-varying parameter O(t). i.e. (1) = [z (4), (()) (compare with Equation 3). It treats (l) as a latent variable and uses ML techniques to infer its value. In contrast, LHM assumes that the expert ODE is self-contained, i.c. with no dependency on any unknown time- varryin parameters or variables. The method in [51] addresses different problem setting, where the expert model is mis-specified. It introduces ML component fm to additively correct for the discrepancy between the expert ODE and the true dynamics f, i.e. f f f It thus relates to the residual models discussed above. Finally, the method in [39] deals with discrete time series and incorporates expert difference equations. It is nontrivial generalize the model and the inference algorithm in [39] to the continuous time setting, which this paper focuses on.
Neural ODEs. Neural ODEs approximate unknown ODEs by neural network [14], frequently using standard feed-forward networks. ODE'VAE uses an architecture with identity blocks to approximate second-order ODEs [95] and GRU-ODE uses an architecture inspired by the Gated Recurrent Unit (GRU) [21, 15]. Neural ODEs and extensions have achieved state-of tfhthc-rt performance a variety of problems involving irregularly-sappe
time series data [72, 21, 45]. We discuss other approaches to learning unknown ODEs from data in Appendix A.5.
Mechanistic models. Mechanistic models are widely applied in sciences such as Physics [81], Epidemiology [92, 31], and Pharmacology [30,42]. These models use ODEs describe a system's continuous-time evolution, possibly under external interventions. The dynamics of the system is
specified deterministically through the governing ODEs; e.g., Equation (3). When the mechanistic models are developed based on experimental data, they may have causal interpretation [76] (Appendix A.5). LHM focuses on the prediction problem and it leverages the scientific or causal knowledge encoded in the expert ODE f to issue principled predictions.
Latent variable models. Latent variable models are widely used in discase progression modeling [87, 4]. These models attempt to infer a set of latent variables to predict complex high-dimensional measurements. The latent variables sometimes have high-level interpretations (e.g. cluster mem- bership), but do not usually correspond to any well-defined and clinically meaningful physiological variable. Moreover, without informative priors, the latent variables can usually be identified only up to certain transformations (e.g. permutation of cluster labels [62]). By contrast, LHM involves medically meaningful expert variables driven by known governing equations and following informative priors.
5 Experiment and evaluation
Here we present the results of two experimental studies, one with simulated data and one with real data. In both experiments, we study the effect of dexamethasone treatment for

---

COVID-19 patients. Both studies are modeled on the real-life treatment of COVID-19 pa- tients in the ICU. The implementation of LHM and the experiment code are available at https://gi /tgt:///////.///..//// wrgss//aaaersshaaalllrrprstttooo
or https: : //gi thub. com/
5.1 Simulation study
In this simulation, use LHM to predict the results of single dexamethasone treatment. Each patient & will receive a one-time treatment; with dosage dN ~ uniform(0, 10] mg and time s; ~ uniform (0, 14]. Our objective is to predict future measurements.
Datasets. We generated variety of datasets to evaluate the model performance under different scenarios. To evaluate how the number of clinical measurements affects performance, we generated datasets with D = 20, 40 or 80 observable physiological variables x. For each dataset, we set the number of un-modeled states z" according to the number of variables in x to be A D/100 2, 4 or 8 (respectively). (We made this choice to reflect the fact that a larger number of physiologica variables often necessitates larger number of un-modeled states.) We consider time horizon of T T= 14 days; this is the median length of stay in hospital for Covid-19 patients [69]. After setting M and D, we generate the data points within dataset independently.
We use a pharmacological model adapted from [18] that describes five expert variables (F 5) under dexamethasone treatment for COVID-19 patients. We use the same model in the real-data experiment. We specify the model and the expert variables Appendix A. The un- modeled states z" are governed by the nonlinear ODE fm shown in equation (5); the true physiological variables x are generated by the function g in Equation (4). The specifications of fr and g are provided in Appendix A.4. For each patient i, each of the components of its initial condition (;( were independently drawn from an exponential distribution with rate À = 100. Measurement noises are drawn independently from Eyt (N,0 o²), for 0.2, 0.4 or 0.8; Equation (1). We first simulate all daily measurements at t = 1, 1,2,... T, and then randomly remove measurements with probability 0.5 to proxy the fact that measurements are made irregularly.
Prediction task For a given patient we use the measurements | : 4o] up to some time to and the treatment plan for that particular patient to predict the future measurements Vf 77. (Note that treatment may have occurred prior to time to or may be planned following time to.) To evaluate the performance under different lengths of observed history, we use to 5, 10 or 12 days. We evaluate the overall prediction accuracy over & E [4o= : 7'] with Root Mean Squared Error (RMSE). We evaluate the uncertainty calibration in the same time window using the Continuous Ranked Probability Score (CRPS). The evaluation metrics for each future time step is shown in Appendix A.4.
Training and Evaluation. We partition each dataset into a training set, validation set, and testing set. We consider training sets consisting of No- = 10, 100, or 500 data points; each validation set has 100 data points and each testing set has 1000 data points.
Benchmarks. We compare the performance of our method (LHM) with the performance of four other methods: latent Neural ODE (NODE), the original Pharmacology model (Expert), the residual model (Residual), and the ensemble model (Ensemble) of Expert and NODE, described below. The details of the optimization and hyper-parameter settings are reported in Appendix A.4.
NODE involves Z latent variables z(t) E RZ whose evolutions f (z(t), a(t); 9Tm), where f is a neural network with trainable weights 6th. The number of latent variables Z is typpe--aramettr that is set to be greater than M - E, which is the number of true latent variables. 3 NODE predict the physiological variables as ý (t) = INv (z(1) a(t); 7N), where ON is a neural network with trainable weights 7N. Expert is given the true governing equation (3), which describes the expert variables z (t), but no un-modeled latent variables. We use a neural network
described by 2(t)
gE with trainable weights YE to predict the physiological variables: (((() = (z (t), a(t); TE).
Residual Given a trained Expert model, we calculate its residuals r(t) y yI) (t((). Then a
NODE is trained to predict the residuals. The final prediction is F + F(l). Ensemble makes
prediction as (t) + (t) where y and Ez(() are the predictions issued by NODE and Expert respectively. The ensemble weights Wie and Wis are learned on the validation set to minimize the prediction error.
'e found that the performance of NODE is not sensitive to the exact choice of Z so long as it is sufficiently larger than M E. (This is consistent with findings reported in the literature [22].)
7

---

Prediction accuracy: RMSE
1.0
1.0
1.0
0.8
0.8
0.8
0.6
0.6
0.6
10
100
500
0.2
0.4
0.8
Sample size
Latent dimension M
Noise level 0
Uncertainty calibration: CRPS
0.8
0.8
0.8
NODE
Expert
0.6
0.6
0.6
Recidual
Ensemble
0.4
0.4
0.4
LHM
10
100 Sample sizw
500
0.2
0.4
0.8
Latent dimension M
Noise level
Figure 3: Simulation Results. Prediction performance on future measurements 3/55 14) given the observed
history 3/0 : 4] measured by RMSE (first row) and CRPS (sccond row). The three columns show the results under (1) different training sample sizes No. (2) different numbers of mo-dodeeed variables M, and (3) different noise levels a. The shaded areas represent 95% confidence intervals.
Results. Figure shows the predictive performance under various samples sizes N, number of latent variables M, and noise levels (results for other settings are reported in Appendix A.4). To contextualize the scale of the performance metric, we report that a naive baseline (predicting a constant number of all individuals and all time steps) achieves an RMSE around 1.0. All methods outperform this naive baseline and LHM achieves the best overall performance. Expert does not perform well because it leaves out the latent variables 2 ((). As result, its performance does not consistently improve when more training samples become available. In contrast, NODE is flexible and fully data-driven. Although its performance improves with increased sample size, NODE is less sample efficient and it achieves worse performance for all the sample sizes considered. NODE is more robust to the increase in .M because it treats the number of latent variables a hyper-parameter We observe that the performance gap between NODE and LHM decreases when more latent variables z" zt"(t) are added. This is because increasing M reduces the proportion of the variables E/M that can be explained by the expert model. Both Residual and Ensemble achieve performance gains over NODE and Expert alone, but they under-perfore LHM because they perform averaging directly in the output space rather than trying to infer z in the latent space.
5.2 Real-data experiments
In this experiment, we use real data to evaluate the predictive performance of LHM and to illustrate its utility for decision support in realistic clinical setting that closely tracks the actual treatment of COVID-19 patients in ICU.
Dataset. We used data from the Dutch Data Warehouse (DDW), multicenter and full-admission anonymized electronic health records database of critically ill COVID-19 patients [27]. Up until March 2021, DDW has collected the health trajectories for 3464 patients admitted to intensive care units (ICU) across the Netherlands. However, even if we use the entire DDW for training, the sample size is still relatively small compared to what is typically used by ML (tens or hundreds of thousands of samples [40]). Furthermore, patients are even scarcer at the early stage of pandemic, arguably when a decision support tool is most needed: only 607 patients were admitted at the first peak (by April 2020). After applying the cligibility criterion detailed in Appendix A.6, we obtained a dataset of 2097 patients whose discase progression is characterized by an sampled time series of 27 physiological variables (Appendix A.6). These variables capture the vital signals, respiratory mechanies, and biomarkers that are crucial for clinical decisions. In addition, we also included 11 static variables that are known to affect the progression of COVID-19, e.g. BMI (Appendix A.6).
Prediction task. Denote to as the time when the patient received the first dose of dexamethasone (we set (o = 24 for untreated patients). We use the history up to 24 hours before lo. V1(4 24 ln). to predict the future 1 4o + 24H] over a time horizon H = 1,3 or 7 days. We use No = 100, 250, 500 or 1000 patients for training, 97 for validation, and 1000 for testing. The pharmacological model and the prior distribution of expert variables are detailed in Appendix A.I.
8

---

Patient A
C-Reactive Protein
Immune Response
25:
Dex Concentration
400
8
300
200
4-
100
2
0
0
0
50
100
150
200
0
50
100
150
200
50
100
150
200
Patient B
X1: C-Reactive Protein
zf:
Immune Response
25:
Dex Concentration
400
8
300
4-
200
4-
100
0
0
50
100
150
200
0
50
100
150
200
50
100
150
200
Patient
X1: C-Reactive Protein
z Immune Response
z5: Dex Concentration
400
8
300
4-
200
100
2
0
50
100
150
200
50
100
150
200
50
100
150
200
Hours after admission
Figure 4: The observed measurements and the inferred expert variables for three illustrative patients. Left: The observed physiological variable (t): C-Reactive Protein. Middle: The inferred expert variable =f (t): the immune response viral infection. Right: The inferred expert variable = (t): dexamethasone concentration. Vertical dotted lines mark the times dexamethasone injections.
Table 1: Prediction accuracy (RMSE) on COVID-19 intensive care data under different training sample sizes N. Prediction horizon H 24 hours. The standard deviations shown in the brackets.
Method \No
100
250
500
1000
Expert
0.718 (0.71)
0.704 (0.02)
0.702 (0.02)
0.713 (0.01)
Residual
0.958 (0.63)
1.003 (0.03)
0.717 (0.05)
0.635 (0.04)
Ensemble
0.707 (0.60)
0.657 (0.05)
0.628 (0.05)
0.599 (0.05)
NODE
0.662 (0.65)
0.659 (0.02)
0.644 (0.05)
(0.04)
ODE2VAE
0.674 (0.62)
0.666 (0.02)
0.643 (0.02)
0.619 (0.02)
GRU-ODE
0.722 (0.60)
0.673 (0.05)
0,623 (0.05)
0,601 (0.05)
Time LSTM
0.706 (0.63)
0.649 (0.03)
0.600 (0.03)
0.631 (0.02)
LHM
0.633 (0.51)
0.605 (0.02)
0.529 (0.02)
0.511 (0.02)
Benchmarks. In addition to all the benchmarks introduced in Section 5.1, we compared the results with two extensions of NODE, GRU-ODE and VAE which achieved strong performance in medical time series prediction [95,21]. We also used the Time L.TM strong baseline [7].
Results. The main results are shown Table I (additional results are shown in Appendix A.6). LHM consistently outperformed the benchmarks. Its performance with No- 100 samples is close to the pure ML approaches' performance with No = 500 samples. As the sample size increases from 100 to 1000, the predictive accuracy of LHM improves by 19% while Time LSTM improves by 11% and NODE by less than 5%. A larger improvement rate suggests LHM adapts to the newly available data faster, which is important when the samples achieved poor performance because it is unable to capture the full array of clinical measurements.
scarce. As expected, the standalone expert model
LHM in action. Here we show how LHM can support clinical decisions beyond predicting future clinical measurements. Managing the level of immune response pivotal when deciding on im- munosuppresive therapy for COVID-19 patients [24, 84, 48]. This is challenging task because

---

the "right" level of immune response varies across patients [43]: for most patients, we would like reduce the immune response avoid cytokine storm and consequent organ failure [24], but for patients with other infections (e.g., secondary bacterial or fungal infection), we would like to keep their immune systems activated [48, 17]. Because immune response not directly observable the clinical settings, clinicians resort to unspecific inflammatory markers such as C-Reactive Protein (CRP) [58]. However, better markers of the immune response such as the cytokine Type IFNs can be measured in the laboratory setting and have been included in the pharmacological model as an expert variable z5. Moreover, the immune response affected by dexamethasone concentrations in the lung tissue (the expert variable zg). These concentrations are not easily routinely measured in a clinical setting, and are therefore not available to elinicians.
Figure shows the measurements of CRP X1, and the inferred immune response z and dexametha- sone concentration 25. The two expert variables are inferred by a trained LHHM using the first five days of observations. We selected three representative patients based on the treatment regimen they re- ceived. Patient A is representative of the 59.8 99.%% of the population who did not receive dexamethasone; Patient C is representative of the 12.2% who received dexamethasone according to the guidelines [60], and Patient is representative of the 28.0% of the population who received dexamethasone but whose treatment was not according to the guidelines.
For patient A, the initial level of CRP was moderately high, but then rose and peaked at about 100 hours after admission. In the absence of contraindications, clinician might begin dexamethasone treatment at this point, but LHM predicts that immune activity of would decrease afterwards even without treatment. (The right panel is blank because dexamethasone was never administered.)
Patient B admitted the ICU with a very high level of the inflammatory marker CRP. Two doses of dexamethasone were given in rapid succession, preceding a decline in both CRP I and immune activity of. However, after the dexamethasone depleted in the patient's body, the expert model predicts that the immune response will pick up again. This is reflected in the re-occurrence of the high CRP level.
Patient was admitted to the ICU with a moderately high level CRP. LHM also inferred a high level of immune activity 2f the initial stage. Inflammation treatments, which has an immunosuppressive effect; so 21 and zf display the same downward trend. Because dexamethasone concentrations falls rapidly within 24 hours after treatment, it repeated at 24 hour intervals, as seen in the right panel. This pattern is clinically expected under the treatment regimen.
greatly reduced after dexamethasone
6 Discussion and future work
This paper has focused on single disease (COVID-19), single treatment (dexamethasone), and single expert model. The ultimate goal is build model that encompasses a variety of diseases, variety of treatments and multiple expert models. This is a challenge for future work.
10

---

Acknowledgments and Disclosure of Funding
This work was supported by the US Office of Naval Research (ONR) and the National Science Foundation (NSF, Grant number: 1722516).
We would like to thank anonymous reviewers as well as members of the vanderschaar-lab for many insightful comments and suggestions.
We would like thank the following individuals for providing domain expertise in Pharmacology and designing the expert ODEs.
Bernhard Steiert, PhD, Pharmaceutical Sciences, Roche Pharma Research and Early Development (pRED), Roche Innovation Center Basel, Basel, Switzerland
Richard Peck, MA, MB, BChir, FFPM, FRCP, Pharmaceutical Sciences, Roche Pharma Research and Early Development (pRED), Roche Innovation Center Basel, Basel, Switzerland
We would like to thank the following individuals for curating and sharing the Dutch Data Warehouse (DDW) for critically ill COVID-19 patients.
From collaborating hospitals having shared data:
Diederik Gommers, MD, PhD, Department of Intensive Care, Erasmus Medical Center, Rotterdam, The Netherlands,
Olaf L. Cremer, MD, PhD, Intensive Care, UMC Utrecht, Utrecht, The Netherlands,
Rob J. Bosman, MD, ICU, OLVG, Amsterdam, The Netherlands,
Sander Rigter, MD, Department of Anesthesiology and Intensive Care, St. Antonius Hospital, Nieuwegein, The Netherlands,
Evert-Jan Wils, MD, PhD, Department of Intensive Care, Franciscus Gasthuis & Vlietland, Rotterdam, The Netherlands,
Tim Frenzel, MD, PhD, Department of Intensive Care Medicine, Radboud University Medical Center, Nijmegen, The Netherlands,
Dave A. Dongelmans, MD, PhD, Department of Intensive Care Medicine, Amsterdam UMC, Amster- dam, The Netherlands,
Remko de Jong. MD, Intensive Care, Bovenij Ziekenhuis, Amsterdam, The Netherlands,
Marco Peters, MD, Intensive Care, Canisius Wilhelmina Ziekenhuis, Nijmegen, The Netherlands,
Marlijn J.A Kamps, MD, Intensive Care, Catharina Ziekenhuis Findhoven, Findhoven, The Nether- lands,
Dharmanand Ramnarain, MD, Department of Intensive Care, ETZ Tilburg, Tilburg, The Netherlands,
Ralph Nowitzky, MD, Intensive Care, HagaZiekenhuis, Den Haag, The Netherlands,
Fleur G.C.A. Nooteboom, MD, Intensive Care, Laurentius Ziekenhuis, Roermond, The Netherlands,
Wouter de Ruijter, MD, PhD, Department Intensive Care Medicine, Northwest Clinics, Alkmaar, The Netherlands,
Louise C. Urlings-Strop. MD, PhD, Intensive Care, Reinier de Graaf Gasthuis, Delft, The Netherlands,
Ellen G.M. Smit, MD, Intensive Care, Spaarne Gasthuis, Haarlem en Hoofddorp. The Netherlands,
D. Jannet Mehagnoul- MD, PhD, Intensive Care, VieCuri Medisch Centrum, Venlo, The Netherlands,
Julia Koeter, MD, Intensive Care, Canisius Wilhelmina Ziekenhuis, Nijmegen, The Netherlands,
Tom Dormans, MD, PhD, Intensive care, Zuyderland MC, Heerlen, The Netherlands,
Cornelis P.C. de Jager, MD, PhD, Department Intensive Care, Jeroen Bosch Ziekenhuis, Den Bosch, The Netherlands,

---

Stefaan H.A. Hendriks, MD, Intensive Care, Albert Schweitzerziekenhui,, Dordrecht, The Nether- lands,
Sefanja Achterberg, MD, PhD, ICU, Haaglanden Medisch Centrum, Den Haag, The Netherlands,
Evelien Oostdijk, MD, PhD, ICU, Maasstad Ziekenhuis Rotterdam, Rotterdam, The Netherlands,
Auke Reidinga, MD, ICU, SEH, BWC, Martiiizickenhuis Groningen, The Netherlands,
Barbara Festen-Spanjee MD, Intensive Care, Ziekenhuis Gelderse Vallei, Ede, The Netherlands,
Gert B. Brumnekreef, MD, Department of Intensive Care, Ziekenhuisgroep Twente, Almelo, The Netherlands,
Alexander D. Cornet, MD, PhD, FRCP, Department of Intensive Care, Medisch Spectrum Twente, Enschede, The Netherlands,
Walter van den Tempel, MD, Department of Intensive Care, Ikazia Ziekenhuis Rotterdam, Rotterdam, The Netherlands,
Age D. Boelens, MD, Anesthesiology. Antonius Ziekenhuis Sneek, Sneek, The Netherlands,
Peter Koetsier, MD, Intensive Care, Medisch Centrum Leeuwarden, Leeuwarden, The Netherlands,
Judith Lens, MD, ICU, UJselllnnd Ziekenhuis, Capelle aan den IJssel, The Netherlands, Roger Rietschote, Business Intelligence, Haaglanden MC, Den The Netherlands,
Harald J. Faber, MD, ICU, WZA, Assen, The Netherlands,
A. Karakus, MD, Department of Intensive Care, Diakonessenhuis Hospital, Utrecht, The Netherlands,
Robert Entjes, MD, Department of Intensive Care, Admiraal De Ruyter Ziekenhuis, Goes, The Netherlands,
Paul de Jong, MD, Department of Anesthesia and Intensive Care, Slingeland Ziekenhuis, Doetinchem, The Netherlands,
Thijs C.D. Rettig, MD, PhD), Department of Intensive Care, Amphia Ziekenhuis, Breda, The Nether- lands,
M.C. Reuland, MD, Department of Intensive Care Medicine, Amsterdam UMC, Universiteit van Amsterdam, Amsterdam, The Netherlands,
Laura van Manen, MD, Department of Intensive Care, BovenIJ Ziekenhuis, Amsterdam, The Nether- lands,
Leon Montenij, MD, PhD, Department of Anesthesiology, Pain Management and Intensive Care, Catharina Ziekenhuis Findhoven, Findhoven, The Netherlands,
Jasper Bommel, MD, PhD, Department of Intensive Care, Erasmus Medical Center, Rotterdam, The Netherlands,
Roy van den Berg, Department of Intensive Care, ETZ Tilburg, Tilburg, The Netherlands,
Ellen Geest, Department of ICMT, Haga Ziekenhuis, Den Haag, The Netherlands, Anisa Hana, MD, PhD, Intensive Care, Laurentius Ziekenhuis, Roermond, The Netherlands,
B. van den Bogaard, MD, PhD, ICU, OLVG, Amsterdam, The Netherlands,
Prof. Peter Pickkers, Department of Intensive Care Medicine, Radboud University Medical Centre, Nijmegen, The Netherlands,
Pim van der Heiden, MD, PhD, Intensive Care, Reinier de Graaf Gasthuis, Delft, The Netherlands,
Claudia (C.W.) van Gemeren, MD, Intensive Care, Spaarne Gasthuis, Haarlem en Hoofddorp. The Netherlands,
Arend Jan Meinders, MD, Department of Internal Medicine and Intensive Care, St Antonius Hospital, Nieuwegein, The Netherlands,
12

---

Martha de Bruin, MD, Department of Intensive Care, Franciscus Gasthuis & Vlietland, Rotterdam, The Netherlands,
Emma Rademaker, MD, MSc, Department of Intensive Care, UMC Utrecht, Utrecht, The Netherlands,
Frits H.M. van Osch, PhD, Department of Clinical Epidemiology, VieCuri Medisch Centrum, Venlo, The Netherlands,
Martijn de Kruif, MD, PhD, Department of Pulmonology, Zuyderland MC, Heerlen, The Netherlands,
Nicolas Schroten, MD, Intensive Care, Albert Schweitzerziikenhuiis Dordrecht, The Netherlands, Klaas Sierk Arnold, MD, Anesthesiology. Antonius Ziekenhuis Sneek, Sneek, The Netherlands,
J.W. Fijen, MD, PhD, Department of Intensive Care, Diakonessenhuis Hospital, Utrecht, The Nether- land,
Jacomar J.M. van Koesveld, MCU, IJsselland Zickenhuis, Capelle aan den [Jssel, The Netherlands,
Koen S. Simons, MD, PhD, Department of Intensive Care, Jeroen Bosch Ziekenhuis, Den Bosch, The Netherlands,
Joost Labout, MD, PhD, ICU, Maasstad Ziekenhuis Rotterdam, The Netherlands,
Bart van de Gaauw, MD, Martiniziekenhuis, Groningen, The Netherlands,
Michael Kuiper, Intensive Care, Medisch Centrum Leeuwarden, Leeuwarden, The Netherlands,
Albertus Beishuizen, MD, PhD, Department of Intensive Care, Medisch Spectrum Twente, Enschede, The Netherlands,
Dennis Geutjes, Department of Information Technology, Slingeland Ziekenhuis, Doetinchem, The Netherlands,
Johan Lutisan, MD, ICU, WZA, Assen, The Netherlands,
Bart Grady, MD, PhD, Department Intensive Care, Ziekenhuisgroep Twente, Almelo, The Netherlands,
Remko van den Akker, Intensive Care, Adrz, Goes, The Netherlands,
Sesmu Arbous, MD, PhD, Intensivist, LUMC, Leiden, The Netherlands,
Tom A. Rijpstra, MD, Department of Intensive Care, Amphia Ziekenhuis, Breda, The Netherlands,
Roos Renckens, MD, PhD, Department of Internl Medicine, Northwest Clinics, Alkmaar, the Netherlands,
From collaborating hospitals having signed the data sharing agreement:
Daniël Pretorius, MD, Department of Intensive Care Medicine, Hospital Jansdal, Harderwijk, The Netherlands,
Menno Beukema, MD, Department of Intensive Care, Streekziekenhuis Koningin Beatrix, Winter- swijk, The Netherlands,
Bram Simons, MD, Intensive Care, Bravis Ziekenhuis, Bergen op Zoom en Roosendaal, The Nether- lands,
A.A. Rijkeboer, MD, ICU, Flevoziekenhuis, Almere, The Netherlands,
Marcel Aries, MD, PhD, MUMC+, University Maastricht, Maastricht, The Netherlands,
Niels Gritters van den Oever, MD, Intensive Care, Treant Zorggroep. Emmen, The Netherlands,
Martijn van Tellingen, MD, EDIC, Department of Intensive Care Medicine, afdeling Intensive Care, ziekenhuis Tjongerschans, Heerenveen, The Netherlands,
Annemieke Dijkstra, MD, Department of Intensive Care Medicine, Het Van Weel-Bethesda Zieken- huis, Dirksland, The Netherlands,
Rutger Raalte, Department of Intensive Care, Tergooi hospital, Hilversum, The Netherlands,
From the Laboratory for Critical Care Computational Intelligence:
13

---

Tariq A. Dam, MD, Department of Intensive Care Medicine, Laboratory for Critical Care Com- putational Intelligence, Amsterdam Medical Data Science, Amsterdam UMC, Vrije Universiteit, Amsterdam, The Netherlands,
Martin E. Haan, MD,, Department of Intensive Care Medicine, Laboratory for Critical Care Com- putational Intelligence, Amsterdam Medical Data Science, Amsterdam UMC, Vrije Universiteit, Amsterdam, The Netherlands,
Mark Hoogendoorn, PhD, Quantitative Data Analytics Group, Department Computer Science, Faculty of Science, VU University, Amsterdam, The Netherlands,
Armand R.J. Girbes, MD, PhD. EDIC, Department of Intensive Care Medicine, Laboratory for Critical Care Computational Intelligence, Amsterdam Medical Data Science, Amsterdam UMC, Vrije Universiteit, Amsterdam, The Netherlands,
Paul W.G. Elbers, MD, PhD, EDIC, Department of Intensive Care Medicine, Laboratory for Criti- cal Care Computational Intelligence, Amsterdam Medical Data Science, Amsterdam UMC, Vrije Universiteit, Amsterdam, The Netherlands,
Patrick J. Thoral, MD, EDIC, Department of Intensive Care Medicine, Laboratory for Critical Care Computational Intelligence, Amsterdam Medical Data Science, Amsterdam UMC, Vrije Universiteit, Amsterdam, The Netherlands,
Dagmar M. Ouweneel, PhD, Department of Intensive Care Medicine, Laboratory for Critical Care Computational Intelligence, Amsterdam Medical Data Science, Amsterdam UMC, Vrije Universiteit, Amsterdam, The Netherlands,
Ronald Driessen, Department of Intensive Care Medicine, Laboratory for Critical Care Computational Intelligence, Amsterdam Medical Data Science, Amsterdam UMC, Vrije Universiteit, Amsterdam, The Netherlands,
Jan Peppink, Department of Intensive Care Medicine, Laboratory for Critical Care Computational Intelligence, Amsterdam Medical Data Science, Amsterdam UMC, Vrije Universiteit, Amsterdam, The Netherlands,
H.J. de Grooth, MD, PhD, Department of Intensive Care Medicine, Laboratory for Critical Care Computational Intelligence, Amsterdam Medical Data Science, Amsterdam UMC, Vrije Universiteit, Amsterdam, The Netherlands,
From Pacmed:
Robbert C.A. Lalisang, MD, Pacmed, Amsterdam, The Netherlands,
Michele Tonutti, MRes, Pacmed, Amsterdam, The Netherlands,
Daan P. de Bruin, MSc, Pacmed, Amsterdam, The Netherlands, Sebastiaan J.J. Vonk, MSc, Pacmed, Amsterdam, The Netherlands,
Mattia Fornasa, PhD, Pacmed, Amsterdam, The Netherlands,
Tomas Machado, Pacmed, Amsterdam, The Netherlands,
Michael de Neree tot Babberich, Pacmed, Amsterdam, The Netherlands,
Olivier Thijssens, MSc, Pacmed, Amsterdam, The Netherlands,
Lot Wagemakers, Pacmed, Amsterdam, The Netherlands,
Hilde G.A. van der Pol, Pacmed, Amsterdam, The Netherlands,
Tom Hendriks, Pacmed, Amsterdam, The Netherlands, Julie Berend, Pacmed, Amsterdam, The Netherlands,
Virginia Ceni Silva, Pacmed, Amsterdam, The Netherlands,
Robert F.J. Kullberg, MD, Pacmed, Amsterdam, The Netherlands,
Taco Houwert, MSe, Pacmed, Amsterdam, The Netherlands,
Hidde Hovenkamp, MSc, Pacmed, Amsterdam, The Netherlands,
14

---

Roberto Noorduijn Londono, MSc, Pacmed, Amsterdam, The Netherlands,
Davide Quintarelli, MSc, Pacmed, Amsterdam, The Netherlands, Martijn G. Scholtemeijer, MD, Pacmed, Amsterdam, The Netherlands,
Aletta A. de Beer, MSc, Pacmed, Amsterdam, The Netherlands, Giovanni Cina, PhD, Pacmed, Amsterdam, The Netherlands, Willem E. Herter, BSc, Pacmed, Amsterdam, The Netherlands, Adam Izdebski, Pacmed, Amsterdam, The Netherlands,
From RCCnet:
Leo Heunks, MD, PhD, Department of Intensive Care Medicine, Amsterdam Medical Data Science, Amsterdam UMC, Vrije Universiteit, Amsterdam, The Netherlands,
Nicole Juffermans, MD, PhD, ICU, OLVG, Amsterdam, The Netherlands,
Arjen J.C. Slooter, MD, PhD, Department of Intensive Care Medicine, UMC Utrecht, Utrecht University, Utrecht, the Netherlands,
From other collaborating partners:
Martijn Beudel, MD, PhD, Department of Neurology, Amsterdam UMC, Universiteit van Amsterdam, Amsterdam, The Netherlands,
References
[1] Abul K Abbas, Andrew H Lichtman, and Shiv Pillai. Cellular and molecular immunology E-book. Elsevier Health Sciences, 2014.
[2] Balaji M Agoram, Steven W Martin, and Piet H van der Graaf. The role of mechanism- basedd pharmackknnetcc-phamacaayaaaadnd (pk-pd) modelling in translational research of biologics. Drug discovery today, 12 2223-24):1088-0224 2007.
[3] Ahmed M Alaa, Scott Hu, and Mihaela Schaar. Learning from clinical judgments: Semi- markov-modulated marked hawkes processes for risk prognosis. In International Conference on Machine Learning, pages 60-69. PMLR, 2017.
[4] Ahmed M Alaa and Mihaela van der Schaar. Attentive state-space modeling of disease pro- gression. In 33rd Conference on Neural Information Processing Systems (NeurIPS 2019), 2019.
[5] Cedric Archambeau, Dan Cornfrrd Manfred Opper, and John Shawe- Taylor. Gaussian process approximations of stochastic differential equations. In Gaussian Processes Practice, pages 1-16. PMLR, 2007.
[6] Jeffrey Aronson and Robin E Ferner. Biomarkers-a a general review. Current protocols in pharmacology, 76(1):9-23, 2017.
[7] Inci M Baytas, Cao Xiao, Xi Zhang, Fei Wang, Anil K Jain, and Jiayu Zhou. Patient subtyp- ing via time-aware Istm networks. Proceedings of the 23rd ACM SIGKDD international conference on knowledge discovery and data mining, pages 65-74, 2017.
[8] Yoshua Bengio and Paolo Frasconi. An input output hmm architecture. Advances in neural information processing systems, pages 427-434, 1995.
[9] Tom Bertalan, Felix Dietrich, Igor Mezić, and Ioannis G Kevrekidis. On learning hamiltonian systems from data. Chaos: An Interdisciplinary Journal of Nonlinear Science, 19112) 2019.
[10] Josh Bongard and Hod Lipson. Automated reverse engineering of nonlinear dynamical systems. Proceedings of the National Academy of Sciences, 104(24):9943-9944 2007.
[11] Steven I. Brunton, Joshua L. Proctor, and Nathan Kutz. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the national academy of sciences, 113(15):3932-393) 2016.
15

---

[12] Giuseppe Carleo and Matthias Troyer. Solving the quantum many-body problem with artificial neural networks. Science, 355(6325):602-6666 606. 2017.
[13] Gang Chen, Yingtao Zuo, Jian Sun, and Yueming Li. cerrort-eecco-mmcciee----- reduced- order model for limit cycle oscillation prediction of nonlinear aeroelastic system. Mathematical problems in engineering, 2012, 2012.
[14] Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud. Neural ordinary differential equations. arXiv preprint arXiv: 80666736, 2018.
[15] Kyunghyun Cho, Bart Van Merriënboer, Caglar Gulcehre. Dzmitry Bahdanau, Fethi Bougares, Holger Schwenk, and Yoshua Bengio. Learning phrase representations using rmn encoder- decoder for statistical machine translation. CoRR, abs/1406. 1078, 2014.
[16) Edward Choi, Mohammad Taha Bahadori, Joshua A Kulas, Andy Schuetz, Walter F Stewart, and Jimeng Sun. Retain: An interpretable predictive model for healthcare using reverse time attention mechanism. Advances in Neural Information Processing Systems, pages 3512-3520, 2016.
[17] Michael J Cox, Nicholas Loman, Debby Bogaert, and Justin O' O'Grady. Co-infections: potentially lethal and unexplored in covid-19. The Lancet Microbe, 1(1)xe11, 2020.
(18) Wei Dai, Rohit Rao, Anna Sher, Nessy Tania, Cynthia J Musante, and Richard Allen A prototype qsp model of the immune response to sars-cov-2 for community development. CPT: pharmacometrics & systems pharmacology, 10(1):18-22) 29, 2021.
[19] Meindert Danhof. Systems pharmacology-toward. the modeling of network interactions.
European Journal of Pharmaceutical Sciences, 94:4-14, 2016.
[20] Meindert Danhof, Joost de Jongh, Elizabeth CM De Lange, Oscar Della Pasqua, Bart A Ploeger, and Rob A Voskuyl. Mcchanism- baaeddd---dddd---e pharmacokinetic- paaemaccokiitic-hhammoylmmmmmm modeling: biophase distribution, receptor theory, and dynamical systems analysis. Annu. Rev. Pharmacol. Toxicol., 47:357 2007.
[21] Edward De Brouwer, Jaak Simm, Adam Arany, and Yves Moreau. Gru-ode-bayes: Continuous modeling of pooaadcally--obervee time series. In 33rd Conference on Neural Information Processing Systems (NeurIPS 2019), 2019.
[22] Emilien Dupont, Arnaud Doucet, and Yee Whye Teh. Augmented neural odes. In 33rd Conference on Neural Information Processing Systems (NeurIPS 2019), 2019.
[23] Justin C Earp, Nancy A Pyszczynski, Diana S Molano, and William J Jusko. Pharmacokinetics of dexamethasone in rat model of rheumatoid arthritis. Biopharmaceuties & drug disposition, 29(6):366-37 3000:3666-7,,00000 2008
[24] David C Fajgenbaum and Carl H June. Cytokine storm. New England Journal of Medicine, 383((3):2255-2D)) 2273,, 2020.
[25] James D Falvey, Teagan Hoskin, Berrie Meijer, Anna Ashcroft, Russell Walmsley, Andrew S Day, and Richard B Gearry. Disease activity assessment in ibd: clinical indices and biomarkers fail to predict endoscopic remission. Inflammatory bowel diseases, 21 2444:8244-811, 2015.
[26] Ferdinando
Fioretto, Terrence WK Mak, and Pascal Van Hentenryck. Predicting ac optimal power flows: Combining deep learning and lagrangian dual methods. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 34 01, pages 630-637, 2020.
[27] Lucas M Fleuren, Daan P de Bruin, Michele Tonutti, Robbert CA Lalisang, and Paul WG Elbers. Large-scale icu data sharing for global collaboration: the first 1633 critically ill covid-19 patients the dutch data warehouse. Intensive care medicine, 47(4):478-481, 2021.
[28] Lucas M Fleuren, Patrick Thoral, Duncan Shillan, Ari Ercole, Paul WG Elbers, and Right Data Right Now Collaborators Mark Hoogendoorn Ben Gibbison Thomas LT Klausch Tingjie Guo Luca F. Roggeveen Elconora L. Swart Armand RJ Girbes. Machine learning intensive care medicine: ready for take-off? Intensive care medicine, 46: 16611416-488 2020.
[29] Richard Frank and Richard Hargreaves. Clinical biomarkers in drug discovery and development. Nature reviews Drug discovery, 2(7):566-580, 2003.
[30] Rudolf Gesztelyi, Judit Zsuga, Adam Kemeny-Beke Balazs Varga, Bela Juhasz, and Arpad Tosaki. The hill equation and the origin of quantitative pharmacology. Archive for history of exact sciences, 66(4) 2012.
16

---

[31]
Michael A Gilchrist and Akira Sasaki. Modeling host host-hooaastt parasite coevolution: a nested approach based on mechanistic models. Journal of Theoretical Biology, 218(3)228 2002.
[32] Sam Greydanus, Misko Dzamba, and Jason Yosinski. Hamiltonian neural networks. In 33rd Conference on Neural Information Processing Systems (NeurIPS 20/9), 2019.
[33] RECOVERY Collaborative Group. Dexamethasone hospitalized patients with covid-19. New England Journal of Medicine, 384(8):693-700) 2021.
[34] Vincent Le Guen and Nicolas Thome. Disentangling physical dynamics from unknown factors for unsupervised video prediction. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 11474-1 1484, 2020.
[35] Franz Hamilton, Alun L Lloyd, and Kevin B Flores. Hybrid modeling and prediction of dynamical systems. PLoS compulational biology, 13(7)=e1005655, 2017.
[36] Jiequn Han, Arnulf Jentzen, and E Weinan. Solving high-dimensional partial differential equations using deep learning. Proceedings of the National Academy of Sciences, 15((4))88505- 8510, 2018
[37] Irina Higgins, Loic Matthey, Arka Pal, Christopher Burgess, Xavier Glorot, Matthew Botvinick, Shakir Mohamed, and Alexander Lerchner. beta-vae: Learning basic visual concepts with constrained variational framework. ICLR, 2017.
[38] Nick Holford, Young-A Heo, and Brian Anderson. A pharmacokinetic standard for babies and adults. Journal of pharmaceutical sciences, (22()):24112955)) 2013.
[39] Zeshan Hussain, Rahul G Krishnan, and David Sontag. Neural pharmacodynamic state space modeling. In Proceedings of the Thirty-Eight Eighth International Conference on Machine Learning (ICML), 2021.
[40] Alistair EW Johnson, Tom Pollard, Lu Shen, H Lehman Li- Li-Wei, Mengling Feng, Mohammad Ghassemi, Benjamin Moody, Peter Szolovits, Leo Anthony Celi, and Roger G Mark. Mimic-iii, a freely accessible critical care database. Scientific data, 3(1):1-9, 2016.
[41] Anuj Karpatne, William Watkins, Jordan Read, and Vipin Kumar. works (pgnn): An application in lake temperature modeling. arXiv preprint arXiv: 11311, 2017.
neural net-
[42] Bertram G Katzung. Basic and clinical pharmacology. Mc Graw Hill, 2012.
[43] Marla J Keller, Elizabeth Kitsis. Shitij Arora, Jen-Ting Chen, Shivani Agarwal, Michael J Ross, Yaron Tomer, and William Southern. Effect of systemic glucocorticoids on mortality or mechanical ventilation in patients with covid- 19. Journal of hospital medicine, 15(8):489 493, 2020.
[44] Shabaana A Khader, Maziar Divangahi, Willem Hanekom, Philip C Hill, Markus Maeurer, Karen W Makar, Katrin D Mayer-Barber, Musa M Mhlanga, Elisa Nemes, Larry S Schlesinger, et al. Targeting innate immunity for tuberculosis vaccination. The Journal of clinical investiga- tion, 129(9):3482-3411 2020.
[45] Patrick Kidger, James Morrill, James Foster, and Terry Lyons. Neural controlled differential equations for irregular time series. In 34th Conference on Neural Information Processing Systems (NeurIPS 2020), 2020.
[46] Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint 141.::1 2014.
[47] Diederik P Kingma and Max Welling, Auto-encoding variational bayes. arXiv preprint 13Xiv:1312.61 14, 2013.
[48] Philipp Koehler, Matteo Bassetti, Arunaloke Chakrabarti, Sharon CA Chen, Arnado Lopes Colombo, Martin Hoenigl. Nikolay Klimko, Cornelia Lass-Flöri, Rita O Oladele, Donald C Vinh, et al. Defining and managing covid-19 1--as-iiiss pulmonary aspergillosis: the 2020 ccmm/isham consensus criteria for research and clinical guidance. The Lancet Infectious Diseases, 2020.
[49] William G Kramer, Richard Lewis, Tyson C Cobb, Wilbur Forester, James A Visconti, Lee A Wanke, Harold G Boxenbaum, and Richard H Reuning. Pharmacokinetics of digoxin: comparison of a two-and three-compartment model in man. Journal of pharmacokinetics and biopharmaceutics, 2(4):299-33 312, 1974.
17

---

[50] Toni Lassila, Andrea Manzoni, Alfio Quarteroni, and Gianluigi Rozza. Model order reduction in fluid dynamies: challenges and perspectives. Reduced Order Methods for modeling and computational reduction, pages 235-273, 2014.
[51] Vincent Le Guen, Yuan Yin, Jérémie Dona, Ibrahim Ayed, Emmanuel de Bézenac, Nicolas Thome, and Patrick Gallinari. Augmenting physical models with deep networks for complex dynamics forecasting, In Ninth International Conference on Learning Representations ICLR 2021, 2021.
[52] Michael D Lee and Wolf Vanpaemel. Determining informative priors for cognitive models. Psychonomic Bulletin & Review, 25(1):114-127, 2018.
[53] Ori Linial, Neta Ravid, Danny Eytan, and Uri Shalit. Generative ode modeling with known unknowns. In Proceedings of the Conference on Health, Inference, and Learning, pages 79-94, 2021.
[54] Dehao Liu and Yan Wang. Multi-fidelity plysics-c neural network and its application in materials modeling. Journal of Mechanical Design, 141(12), 2019.
[55] Carl M Metzler. Usefulness of the two-compartment open model in pharmacokinetics. Journal of the American Statistical Association, 66(333):49-53, 1971.
[56] Andrew C Miller, Nicholas Foti, and Emily Fox. Learning insulin-glucose dynamics in the wild. In Machine Learning for Healthcare Conference, pages 172-197. PMLR, 2020.
[57] John B Moore and Carl H June. Cytokine release syndrome in severe covid-19, 19. Science, 368(6490):473-474 2020.
[58] Richard Mortensen. C-reactive protein, inflammation, and innate immunity. Immunologic research, 24(2):163-177) 2001.
[59] Nikhil Muralidhar, Mohammad Raihanul Islam, Manish Marwah, Anuj Karpatne, and Naren Ramakrishnan. Incorporating prior domain knowledge into deep neural networks. 2018
IEEE International Conference on Big Data (Big Data), pages 36-45. IEEE, 2018.
[60] NIH. Therapeutic management of adults with covid-19.
attps://wwt Accessed:
@wvid19treatmeen iegrreattennguudellenee nih, http/thherppeticc-aaaeeeee 2021-05-26.
[61] Alexander Norcliffe, Cristian Bodnar, Ben Day, Jacob Moss, and Pietro Lio. Neural ode
processes. ICLR, 2021.
[62] Ari Pakman, Yueqi Wang, Catalin Mitelut, JinHyung Lee, and Liam Paninski. Neural clustering processes. In International Conference on Machine Learning, pages 7455-7465. PMLR, 2020.
[63] Giambattista Parascandolo, Niki Kilbertus, Mateo Rojas-Carulla, and Bernhard Schölkopf. Learning independent causal mechanisms. In International Conference on Machine Learning, pages 4036-4044. PMLR, 2018.
[64] Eric Parish and Karthik Duraisamy. A paradigm for data-driven predictive modeling using field inversion and machine learning. Journal of Computational Physics, 305:758-774, 2016.
[65] Judea Pearl. Causal diagrams for empirical research. Biometrika, 82(4))669 1995.
[66] Judea Pearl. Causality. Cambridge university press, 2009.
[67]
Lawrence Perko. Differential equations and dynamical systems, volume 7. Springer Science & Business Media, 2013.
[68] Jordan S Read, Xiaowei Jia, Jared Willard, Alison P Appling, Jacob Zwart, Samantha Oliver, Anuj Karpatne, Gretchen JA Hansen, Paul C Hanson, William Watkins, et al. Process-guided deep learning predictions of lake water temperature. Water Resources Research, 55(11):9173- 9190, 2019.
[69] Eleanor M Rees, Emily S Nightingale, Yalda Jafari, Naomi R Waterlow, Samuel Clifford, Carl AB Pearson, Thibaut Jombart, Simon R Procter, Gwenan M Knight, CMMID Working Group, et al. Covid-19 length of hospital stay: a systematic review and data synthesis. BMC medicine, 1-(1):1-22, 2020.
[70] Danilo Rezende and Shakir Mohamed. Variational inference with normalizing flows. In International Conference Machine Learning, pages 1530-1538. PMLR, 2015.
18

---

[71] Maria Isabel Ribeiro. Kalman and extended kalman filters: Concept, derivation and properties. Institute for Systems and Robotics, 43:46, 2004.
[72] Yulia Rubanova, Ricky TQ Chen, and David Duvenaud. Latent ordinary differential equations for iregularly-sampled time series. In 33th Conference on Neural Information Processing Systems (NeurIPS 2019), 2019.
[73] Donald B Rubin. Causal inference using potential outcomes: Design, modeling, decisions.
Journal of the American Statistical Association, 100(469):32-3331 2005.
[74] Samuel H Rudy, Steven L Brunton, Joshua L Proctor, and Nathan Kutz. Data-driven discovery of partial differential equations. Science Advances. 3(4)xe1602614 2017.
[75] Andreas Ruttor, Philipp Batz, and Manfred Opper. Approximate gaussian process inference for the drift function in stochastic differential equations. In Advances in Neural Information Processing Systems, pages 2040-2048. Citeseer, 2013.
[76] Bernhard Schölkopf, Francesco Locatello, Stefan Bauer, Nan Rosemary Ke, Nal Kalchbrenner, Anirudh Goyal, and Yoshua Bengio. Toward causal representation learning. Proceedings of the IEEE, 10965)x12-6334 2021.
[77] Kristof T Schütt, Pieter-Jan Kindermans, Huziel E Sauceda, Stefan Chmicla, Alexandre Tkatchenko, and Klaus-Robert Müller. Schnet: A continuous-filter convolutional neural network for modeling quantum interactions. arXiv preprint arXiv: (ttv::/706.856000 08566, 2017.
[78] Duncan Shillan, Jonathan AC Sterne, Alan Champneys, and Ben Gibbison. Use of machine learning to analyse routinely collected intensive care unit data: a systematic review. Critical Care, 23(1):1-11, 2019.
[79] Justin Sirignano and Konstantinos Spiliopoulos. Dgm: A deep learning algorithm for solving partial differential equations. Journal of computational physics, 375:1339-1364, 2018.
[80] Simone MC Spoorenberg, Vera HM Deneer, Jan C Grutters, Astrid E Pulles, GP Voorn, Ger T Rijkers, Willem Jan W Bos, and Ewoudt MW van de Garde. Pharmacokinetics of oral vs. intravenous dexamethasone in patients hospitalized with community-acquiree pneumonia. British journal of clinical pharmacology, 78(1):7 (171))88777)) 83, 2014.
[81] William James Stronge. Impact mechanics. Cambridge university press, 2018.
[82]
Nathanicl Thomas, Tess Smidt, Steven Kearnes, Lusann Yang, Li Li, Kai Kohlhoff, and Patrick Riley. Tensor field networks: Rotation-ant translation-equivariant neural networks for 3d point clouds. arXiv preprint vrXiv:1802.082/13 2018.
[83] Bruno M Tomazini, Israel Maia, Alexandre B Cavalcanti, Otavio Berwanger, Regis Rosa, Viviane C Veiga, Alvaro Avezum, Renato D Lopes, Flavia R Bueno, Maria Vitoria AO Silva, al. Effect of dexamethasone on days alive and ventilator- free in patients with moderate or severe acute respiratory distress syndrome and covid-19; the codex randomized clinical trial. 2020.
Jama, 324(11):
[84] Judith van Paassen, Jeroen S Vos, Eva M Hoekstra, Katinka MI Neumann, Pauline C Boot, and Sesmu M Arbous. Corticosteroid use in covid-19 patients: systematic review and meta-analysis on clinical outcomes. Critical Care, 24(1):1-22, 2020.
[85] Michel Verhaegen. Identification of the deterministic part of mimo state space models given in innovations form from onput-output data. Automatica, 30(1):61-74 1994.
[86] Jian-Xun Wang, Jin-Long Wu, and Heng Xiao. Physics-informed machine learning approach for reconstructing reynolds stress modeling discrepancies based on dns data. Physical Review Fluids, 2(3):034 2017.
[87] Xiang Wang, David Sontag, and Fei Wang. Unsupervised learning of disease progression models. In Proceedings of the 20th ACM SIGKDD international conference on Knowledge discovery and data mining, pages 85-94, 2014.
[88] Jared Willard, Xiaowei Jia, Shaoming Xu, Michael Steinbach, and Vipin Kumar. Integrating physics-based modeling with machine learning: A survey. arXiv preprint arXiv:2003.04919. 2020.
[89] Jin-Long Wu, Heng Xiao, and Eric Paterson. Physics-infrrmed machine learning approach for augmenting turbulence models: A comprehensive framework. Physical Review Fluids, 3(7):074 2018
19

---

[90] D Xiao, CE Heaney, L Mottet, F Fang. W Lin, IM Navon, Y Guo, OK Matar, AG Robins, and CC Pain. reduced order model for turbulent flows in the urban environment using machine learning. Building and Environment, 148:323-337, 2019.
[91] Tianfang Xu and Albert Valocchi. Data-driven methods to improve baseflow prediction of a regional groundwater model. Computers & Geosciences, 85:124-136, 2015.
[92] Zifeng Yang, Zhiqi Zeng, Ke Wang. Sook-San Wong, Wenhua Liang, Mark Zanin, Peng Liu, Xudong Cao, Zhongqiang Gao, Zhitong Mai, et al. Modified seir and ai prediction of the epidemics trend of covid-19 in china under public health interventions. Journal of thoracic disease, 12(3):165, 2020.
[93] Kun Yao, John E Herr, David Toth, Ryker Mckintyre, and John Parkhill. The tensormol-0.1 model chemistry: neural network augmented with long-range physics. Chemical science, 98):222612269, 2018.
[94] Alireza Yazdani, Lu Lu. Maziar Raissi, and George Em Karniadakis. Systems biology informed deep learning for inferring parameters and hidden dynamics. PLoS computational biology, 16115):e 2020.
[95] Cagatay Yildiz, Markus Heinonen, and Harri Lähdesmäki. Ode2vae: Deep generative second order odes with bayesian neural networks. In Proceedings of the 33rd Conference on Neural Information Processing Systems (NeurIPS 2019), 2019.
[96] Leonardo Zepeda-Núñez Yixiao Chen, Jiefu Zhang, Weile Jia, Linfeng Zhang, and Lin Lin. Deep density: circumventing the kohn-sham equations via symmetry preserving neural networks. arXiv preprint 1912.00775 2019.
[97] Cheng Zhang. Judith Bütepage, Hedvig Kjellström, and Stephan Mandt. Advances in variational inference. IEEE transactions on pattern analysis and machine intelligence, 41(8):2008 2000220666 2018.
[98] Liang Zhang, Gang Wang, and Georgios B Giannakis. Real-time power system state estimation via deep unrolled neural networks. In 2018 IEEE Globul Conference Signal and Information Processing (GlobalSIP), pages 907-911. 1. IEEE, 2018.
[99] Yaofeng Desmond Zhong, Biswadip Dey, and Amit Chakraborty. Symplectic ode-net: Learning hamiltonian dynamics with control. In International Conference Learning Representations, 2019.
20


Reviews:
[{"review_id": "xeSlvl60UOW", "reviewer": "Reviewer_62Te", "summary": "This paper proposes an architecture to mix Neural ODE models and more structured ODE coming from expert knowledge. They evaluate their approach on a dataset from COVID-19 patients using an expert ODE of the dexamethasone.", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "I think this is a good paper with an interesting idea to mix expert knowledge and ODE models. However, the experiments section focuses on a single ODE /dataset, on which the results are compelling. Importantly, to understand if the performance gain comes from the expert ODE or from a better behaved ODE, extra experiment is needed. \n\n**Model formulation**\n- Why is it beneficial to add the current treatment value as an input to the observation function. Indeed, the treatment is already encoded in the model through the ODE of the latent process (equations 3. and 5.)\n- $a(t)$ is considered continuous, so it requires some assumption to cast any type of treatment in this form. This would need to be discussed.\n\n**Related works**\n- Related work on mechanistic/pharmacology models embedded in ML architectures : I think you don't cite reference [1]. Maybe too recent but would be good to discuss the difference with your work.\n\n- I think your work is very much linked to this paper as well : [2]\n\n**Experiments :**\n\n- From the synthetic experiment, I would not be surprised that the Expert performs worse than LHM even if M=0, simply because the number of parameters is higher for LHM. It would be good to have expand the plots of Figure 3 with M=0, as this would give a better view of the MSE obtained for larger M.\n\n- Similarly as in the paper of Hussain et al. [1], it would be nice to assess the effect of a mispecification of the expert ODE. Indeed, it might be that part of the improvement over an unconstrained Neural-ODE is the fact that part of the ODE is well-behaved [3] (low Lipschitz constant). In practice this would just mean to encode in the model a \\emph{wrong} parametric form for $f^e$. This would help narrowing down the source of the performance gain.\n\n[1] Neural Pharmacodynamic State Space Modeling Zeshan Hussain, Rahul G. Krishnan, David Sontag\n[2] Linial, O., Ravid, N., Eytan, D., & Shalit, U. (2021, April). Generative ODE modeling with known unknowns. In Proceedings of the Conference on Health, Inference, and Learning (pp. 79-94).\n[3] How to Train Your Neural ODE: the World of Jacobian and Kinetic Regularization, Finlay et al.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "They have. I think one limitation is the source of improvement of the model, which would require an extra experiment.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "2.5"}, "extra_text": {}}, {"review_id": "mgdbCjYVP9U", "reviewer": "Reviewer_yPRn", "summary": "A VAE model in which the decoder includes an ODE describing the evolution of some non-directly observable quantities, and the encoder learns mappings from the observed measurements to the initial states of the non-observed quantities, as well as to the ODE parameters. The model is used for forecasting in interesting application that seeks to forecast response to treatment of COVID-19 patients. \n\nI have some concerns with respect to the novelty of the paper, since it seems to be a special case of recent work, as well with respect to how the model might choose to treat the ODE part of it. In the general case the latter might simply be ignored since there is no protection against something like that baked in the model. ", "questions": "", "limitations": "", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper presents an interesting application in forecasting the results of different treatments on patients with COVID-19.\nIt proposes a learning architecture based on a VAE model where the decoder contains an expert model, a pharmacological ODE, \nthat describes the evolution of different quantities of interest, collectively denoted by $\\mathbf z^e(t)$, such as immune \nresponce to infection, concentration of the drug in the body, viral load. Typically these variables are not directly observable \nin the routine clinical setting, but only in laboratory.  Only the structure of the ODE is known but not its parameters $\\theta$.\n\nThe clinician has direct access only to a set of observable variables, $\\mathbf x(t)$, or more precisely to their noisy measurements, \n$\\mathbf y(t)$, measuring quantities such as temperature, heart rate, and more. In addition the clinician also provides the treatment\nwhich is also modelled as a temporal sequence $\\mathbf a(t)$.  Since the non-observable variables $\\mathbf z^e(t)$ \nof the pharmacological ODE explain only one part of the process that produces the observational data these are completed with an \nadditional set of variables $\\mathbf z^m(t)$. \n\nThe VAE model is structured as follows: \n* the encoder receives the observation data $\\mathbf y(t)$ and treatment $a(t)$ and encodes them to the parameters of the ODE, $\\theta$, and the \ninitial states $\\mathbf z^e(0)$  $\\mathbf z^m(0)$; these are then respectivelly passed through the pharmacological ODE and a learned ODE to produce \nthe complete sequences $\\mathbf z^e(t), \\mathbf z^m(t)$\n* the decoder receives the latent variables and reconstructs the input signal. \n\nIn forecasting time the model receives the $\\mathbf y(t)$ and $\\mathbf a(t)$ up to the present time, \nprovides the posterior, and the decoder will use the sampled initial states and parameters of the\npharmacological ODE, to conditionally generate the future observation trajectory given also a future \ntreatment plan.\n\n\nThe paper provides two sets of experiments one on syntetic data generated using the known pharmacological model and one on real world COViD-19 patient \ndata. The model is compared against a number of baselines that involve a variant in which the pharmacological ODE is removed (NODE), one which is using \nonly the pharmacological ODE (Expert), two that rely in combinations of them, and a additional number including autoregressive models. \n\n\nThe paper is quite well written and easy to follow. I have though a number of concerns. \n\nThe main one has to do with whether indeed the model is really using the pharmacological ODE and the resulting $\\mathbf z^e(t)$ latent variables in its generative\nprocess. Since the latent space has two components $\\mathbf z^e(t)$ and  $\\mathbf z^m(t)$ what is there that guarantees that all the information flow will not happen\ntrhough $\\mathbf z^m(t)$ and the decoder $g$ will learn to ignore it, reducing the model to the NODE baseline. The empirical results show that probably this is not \nhappening, since the proposed model outperforms systematically NODE. However there is nothing in the model architecture and objective function that would prevent that\nfrom happening in the general case. The end result would be that we will have a model that while we believe it is based in domain knowledge it has learned to circumvent\nthat knowledge. \n\nIn fact such a setting has been extensively studied in Takeishi, Kalousis, Physics-Integrated Variational Autoencoders for Robust and Interpretable \nGenerative Modeling, 2021, which seems to subsume LHM as a special case. In that work a number of regularisers have been proposed in an effort to guarantee that the \nODE will not be ignored and the additional non ODE latent variables will be used in as minimal way as possible. Moreover there is no assumption about the completeness \nof the provided ODE, i.e. the model can work with incomplete knowledge. \n\nVery relevant work, though not in a generative setting and thus no interpretable latent variables, is also Yin et al, Augmenting Physical Models with Deep Networks for \nComplex Dynamics Forecasting, ICLR 2021. That work learns a dynamics model to be used for forecasting by learning augment {\\em incomplete} ODEs. Its model has two dynamics \ncomponents; one directly relies on the incomplete ODE and one is data driven and additively completes the ODE based dynamics. Here too the model makes sure, through appropriate \nregularisation, that the data-driven component will only capture information that is not already contained in the physics component, making sure in that manner that it does \nnot overwrite the ODE based compelent.  \n\n\nExperiments: \nIn the synthetic experiments, since we have complete access to all data, including the non-observed latents, it would have been very informative:\n*  to show how well the encoder is able to predict the $\\theta^e$ parameters of the ODE, \n*  and compare the forecasted/completed curves of $\\mathbf y(t)$, $\\mathbf z^{e}$, produced by the model with the real ones. Such a figure, fig 4, is \nonly partially given for $\\mathbf z^{e}$ for the real data, which is expected since there we do not have access to the real $\\mathbf z^{e}$. However \nthere too it would be quite interesting to see how the forecasted $\\mathbf y(t), t>t_0$ values compare with the ones that actually occur in the real data. \n* the paper also mentions a forecasting horizon of 1, 3, and 7 days, the results of which I found in the appendix,  but it is not clear to which horizon \nthe table 1 in the main paper corresponds to. Again some visualisations of the actual dynamics would have also been quite useful here, to see how the\nperformance scales as we move further to the future. \n\n\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Given the application domain, that of forecasting for medical purposes, it is rather important to make sure that the model operates as it is intended. This goes back to my comment above on whether indeed the model makes use of its ODE or simply chooses to ignore it, since as the model is there is no explicit control to avoid such a scenario. The end result would be a model that is not driven by the pharmacology model, the latent variables $\\mathbf z^e(t)$ would be meaningless. That does not necessarily mean that the forecasting performance would be poor, it would rather be that of a classical data driven model.\nAs a consequence one would need to include mechanisms to avoid such case from happening, see related work above, as well as mechanisms that are able to detect the extend to which such a case happens, e.g. some sensitivity analysis of $g$ with respect to $\\mathbf z^e(t)$ and $\\mathbf z^m(t)$. \n \nGiven the sensitive application domain such a discussion is more than relevant.  ", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "6"}, "extra_text": {}}, {"review_id": "VihnYAeVj1_", "reviewer": "Reviewer_Cu63", "summary": "This paper introduces an approach to integrate expert information into learned ODE models of observed dynamic processes. Motivated by applications in healthcare (particularly in pharmacology and disease progression) the expert information is provided as a verified dynamics model of lab-measured features that may not be readily measured in the normal clinical practice. However, the inductive biases afforded on the physiological dynamics, conditioned on possible treatment interventions, are considered worthwhile in the task of predicting future patient measurements with the learned ODE model. The two models are brought together by inferring the latent variables of the unobserved \"expert\" measurements and how they may affect the latent state generating the observed measurements. The proposed hybrid approach is evaluated on synthetic as well as real ICU data, comparing with several contemporary nODE appraoches.", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "\n\n#### **Originality**\nWhile the use of expert models or inductive biases to accelerate or improve model performance is not new, the use of a factored latent space between joint ODEs is unique and interesting. The authors do a great job listing and comparing to contemporary methods and approaches in the literature. One omission that should be included and considered is that of Miller, et al (2020 MLHC; reference below) which learns an improved dynamics model of the Insulin-Glucose interaction dynamics from a pre-defined ODE and observed patient data. I would be interested to hear from the authors how they feel their proposed approach differs from this work.\n\n\n#### **Quality and Clarity**\nI found the paper to be very well written with an appropriate level of technical detail and formalization such that the construction of the proposed LHM was easy to follow. I particularly enjoyed the measured discussion that was included to discuss both potential benefits of LHM in practice (in particular, Section 3.5 was well formed) as well as the current limitations of the approach (e.g independence of noise terms vs. autocorrelation, possible expert model misspecification, etc). I did have a few questions however and look forward to the author’s clarifications that they are able to provide.\n\nWas the claim on line 51 that the inferred expert variables provide additional insights to clinicians ever evaluated? If not, the hypothetical nature of this claim should be clarified. Or is the evidence of this statement provided in the analysis presented around Figure 4? \n\nIn the modeling of the expected measurements over the future time interval, does the length of that interval affect the confidence or uncertainty associated with the predictions? While these quantities are referred to toward the end of Section 2, the experimental results do not fully demonstrate the effect of prediction horizon on model confidence. It’s unclear what the choice of $t_0$ has on the downstream performance of the predictions. Does the model appropriately extrapolate over long horizons? How do those predictions improve with longer histories ahead of the prediction interval? Is there a maximum future time window for which the predictions are reliable?\n\nThis line of questioning led me to wonder what the specified prediction time was for the ICU dexamethasone experiments included in Section 5.2? When were the predictions made? At the first treatment? At distinct time points?\n\nThe use of the term “expert variables” is quite vague throughout the paper. More concrete examples of what these variables may be in the context of the conditions investigated in the paper would be really helpful to get a sense of what the proposed LHM is expected to learn. This would be most helpful in the early stages of the paper. One particular area could be at the outset of Section 3.2. \n\nIn the final paragraph of the “Datasets” subsection of Section 5.1, a range for $\\sigma$ is provided. However, in Section 3.2 this variable is included in the learnable parameters $\\Theta$? Are the experiments in Section 5.1 investigating the effects of a fixed $\\sigma$ parameter?\n\n\n#### **Significance**\nWhile I don’t feel especially qualified to comment on the significance of this paper. I do feel that the contributions put forward by the work are potentially of high impact as far as the type of expert modeling needed for LHM is available and appropriately specified for the desired use case. \n\n\n##### **Additional references**\nMiller, Andrew C., Nicholas J. Foti, and Emily Fox. \"Learning Insulin-Glucose Dynamics in the Wild.\" Machine Learning for Healthcare Conference. PMLR, 2020.\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "As mentioned above, the authors sufficiently address the current technical limitations of the proposed LHM approach. There are no expressed considerations or limitations in view of potential applications of this model in live situations. However, at this time this research appears to be largely establishing a proof of concept for circumstances when a fully validated (and trusted) expert model exists and LHM serves to only adapt and improve that model based on observed data (without considerations toward equitable performance across possible patient subgroups).", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "10"}, "extra_text": {"ethical_concerns": "I found there to be no ethical concerns with the proposed approaches and methods used in this paper. One potential area could simply be the use of unverified or poorly calibrated models in clinical settings. This paper makes no claim for immediate adoption nor is it anywhere close to being applied. "}}, {"review_id": "85UdvI_snxG", "reviewer": "Reviewer_59iy", "summary": "In this work the authors proposed the latent hybridisation model to predict disease progression from physiological measurements. The proposed model combines ordinary differential equations designed by experts and the neural ODE approach to achieve prediction accuracy and the interpretability of the expert variables. The performance was evaluated on a synthetic data as well as a real-world dataset. Improved performance was demonstrated comparing with several competing models.", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "I think the authors demonstrated a nice idea to integrate Pharmacological models as prior knowledge for a generic machine learning model. The work is well-presented, and I enjoyed reading the article. \nThe implementation seems follow standard approach of variational inference and limited technical contribution was presented. The evaluation demonstrates improvement in performance comparing to several strong baselines on a real-word ICU dataset. It would be nice to have another dataset to validate the improvement. The authors also demonstrated the benefit in interpretability with qualitative analysis.\n\nDetailed comments:\n\nHow are the parameters in the expert ODEs determined? Can the authors provide a table for the values of ODE parameters?\n\nCan the authors provide interpretation on the values on RMSE, i.e, how good/off a predicted track is for RMSE ~ 0.6. It would be very helpful to plot a few time series comparing the predicted ones v.s. ovserved ones. Also, it would be helpful to provide the definition or reference for RMSE / CRPS in line 254-255.\n\nCan the authors provide details on the results presented in Fig 4? Specifically, are the selected trajectories from the training or test split? Are all observed time points used for reconstructing the time courses? I am confused that the blue curve in x1 looks like a perfect fitting to the datapoints with no prediction errors. \n\nAre the time courses of Z2 in figure 4 just the convoluted responses of drug treatments? If so, it seems not as informative as claimed.\n\nIs there any evidence supporting that z1 represents the real response of cytokine Type I IFN?", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors adequately addressed the limitations", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}]

Author rebuttal:

Meta-review: First, thanks to the authors for this engaging submission on an important topic.  While the reviewers expressed some concern about understanding the method, the back and forth between the reviewers and authors appears to have addressed most of those concerns. 

Reviewer 62Te was concerned about the source of the improvement, given all of the moving pieces in the proposed approach.  The additional simulation study has shed more light on the behavior of the proposed algorithm. 

Reviewer yPRn voiced concerns about novelty, and how the proposed approach is related to two recent physics + ML algorithms.  The authors addressed these concerns in their response, and I encourage them to more explicitly compare these methods to their own in the manuscript. 

Decision: Accept (Poster)