Benchmarking Continuous Time Models for Predicting Multiple Sclerosis Progression
Authors: Alexander Luke Ian Norcliffe
Lev Proleev
Diana Mincu
F Lee Hartsell
Katherine A Heller
Subhrajit Roy
Venue: Accepted by TMLR ()

Abstract: Multiple sclerosis is a disease that affects the brain and spinal cord, it can lead to severe disability and has no known cure. The majority of prior work in machine learning for multiple sclerosis has been centered around using Magnetic Resonance Imaging scans or laboratory tests; these modalities are both expensive to acquire and can be unreliable. In a recent paper it was shown that disease progression can be predicted effectively using performance outcome measures and demographic data. In our work we build on this to investigate the modeling side, using continuous time models to predict progression. We benchmark four continuous time models using a publicly available multiple sclerosis dataset. We find that the best continuous model is often able to outperform the best benchmarked discrete time model. We also carry out an extensive ablation to discover the sources of performance gains, we find that standardizing existing features leads to a larger performance increase than interpolating missing features.

Published in Transactions on Machine Learning Research (09/2023)
Benchmarking Continuous Time Models for Predicting Multiple Sclerosis Progression
Alexander Norcliffe*
olex. mhorckifee&&@maat com
University of Cambridge
Lev Proleev Google Research
lerp&ggoogle.c com
Diana Mincu Google Research
dmincu@@google.com com
Fletcher Lee Hartsell
edu
Duke University Health System
Katherine Heller
kheller@google.cco
Google Research
Subhrajit Roy Google Research for MSOACE
subherajttroy@gggonll..com
Reviewed on OpenReview: http:://ppaaeeriee.nn////ottt//p
Abstract
Multiple sclerosis is discase that affects the brain and spinal cord, it can lead to severe disability and has no known cure. The majority of prior work in machine learning for multiple sclerosis has been centered around using Magnetic Resonance Imaging scans or laboratory tests; these modalities are both expensive to acquire and can unreliable. In a recent paper was shown that disease progression can be predicted effectively using performance outcome measures and demographic data. In our work build this investigate the modeling side, using continuous time models to predict progression. We benchmark four continuous time models using publicly available multiple sclerosis dataset. We find that the best continuous model is often able to outperform the best benchmarked discrete time model. We also carry out an extensive ablation to discover the sources of performance gains, we find that standardizing existing features leads to a larger performance increase than interpolating missing features.
An associated video presentation available at: BqDLDK6P2H0
1 Introduction
Multiple Sclerosis (MS) is the most common demyelinating discase (Leray et al. 2016) affecting approxi- mately 2.8 million people worldwide (Lane et al. 2022 and being one of the leading causes of disability in
Work done as student rescarcher at Google Research from June-December 2022.
!ooresponning author.
!Data used in the preparation of this article were obtained from the Multiple Sclerosis Outcome Assessments Consortium (MSOAC). such, the investigators within MSOAC contributed the design and implementation the MSOAC Placebo database and/or oor provided placebo data, but did not participate in the analysis the data the writing of this report.

---

Published in Transactions on Machine Learning Research (09/2023)
young adults (Coetzee & Thompson 2020). It is lifelong disease of the nervous system that can lead to wide range of symptoms from latigue to problems with movement and vision. Currently there is no known cure but the symptoms can be treated with medication and physical therapy, where the aim of a clinician is to limit progression. Therefore, an effective application of machine learning methods in this space would be to predict progression of MS as well as informing interventions and treatment options to obtain the best prognosis for a given patient.
Prior machine learning works that investigate MS progression typically rely Magnetic Resonance Imaging (MRI) scans or laboratory tests (Pinto et a 2020 Zhao et al. 2017 Seccia et al., 2021). These modalities are both expensive to obtain (MRI scans cost on the order of thousands of US
2020
Tommasin et al.
dollars) and infrequent, added to this they are not always reliable (Whitaker et al.
1995
The progression
of MS has also been studied using clinical data (Lawton et al., 2015 Tozlu et al.
2019 De Brouwer et al.
2021) with models ranging from Latent Class Linear Mixed Models, to time-aware GRU.
To this end, in previous study Roy et al. (0022) demonstrated that it is possible to use both performance outcome measures and demographic data to predict disease progression. Performance outcome measures (POMs) are a time-series that contain either responses to questionnaires or results from physical tests to determine MS progression (Rudick et al. 2002). These tests are designed to evaluate the functions affected by MS, and include walking, balance, cognition and dexterity tests. Roy et al. (2022) showed that POMs and demographic data successfully predict progression for a range of models, from classical models such as logistic regression to deep learning models for sequences such as Temporal Convolutional Network (TCN) (Bai et al. 2018). However, these models are not suited to irregular time-series with missing features, and various techniques such imputing with zeros and binning are applied to make the data compatible with the models, which can add bias. In contrast, these realistic hurdles of irregular and missing measurements can be handled in principle by continuous time models. Continuous time models such as Neural ODEs (Chen et al. 2018) learn to model an instantaneous velocity of a state to calculate the time evolution. They have three key advantages: (1) they naturally handle irregular timo-sorics; (2) via the use of control signals missing features can be interpolated and therefore are imputed with more physically plausible values (Kidger et al. 2020); (3) continuous time models enforce the inductive bias that the data generated from dynamical process such as biological interactions (Alharbi & Rambely 2020). Neural Differential Equations have been applied in the
medical setting (Kidger et al.) 2020 Morrill et al.) 2021 Qian al. 2021 Seedat al.) 2022 Gwak et al,
2020), however to our best. knowledge we are the first to apply these to MS prediction using POMs. Our
primary contributions are: (1) We benchmark four continuous time models on public MS dataset, and see that the best contimuus model outperforms the previously best benchmarked discrete model. (2) We test the use of interpolation to fill missing values rather than fixed values, in this setting we find the difference is negligible. To carry this out we faced various challenges requiring technical solutions: (1) We adapt the framework Roy al. (2022) to handle continuous models which impute by interpolation. (2) We generalize existing method by Morrill et al. (2021) to train continuous models with irregular data in batches. (3) We introduce two new interpolation schemes using monotonically increasing cubics to interpolate missing values and be control signals for models. We view these as secondory contributions and emphasize that the main message of the paper is that continuous models can be applied successfully in the MS setting.
2 Background
2.1
Neural Differential Equations
Here we briefly introduce the concepts behind continuous time models. For more information we encourage the reader to see resources such as Chen et al.
(2018); Dupont et al. (201 Massaroli et al. (2020); Kidger
ct al. 2020)): Li et al. (2020): Kidger (2022). For hidden state h E RD. and time E R, Neural Ordinary 2018) model the instantaneous rale of change of h with x e + R", with parameters 0 E θ
Dilferential Equations (Neural ODEs) CChen et al.
respect to t using a neural network f
R" x
R
dh
= t(h,t,
(1)
dt
2

---

Published in Transactions on Machine Learning Research (09/2023)
For brevity we now use the convention that f(h,t, 0) is written fo(h, t). Given an initial condition h(la)= = ho, we can find the solution at general time h(()) using ODE solvers. The simplest is the Euler solver, where we start from the known (tt)) and repeatedly apply the update rule h(t+At) = h(t) + f(((t), t) this approaches the true solution as At approaches O. For more accurate solutions we can look to the existing numerics work and use advanced black box ODE solvers, for example Dormand-Prinee (Dormand & Prince 1980 to approximately solve the initial value problem and calculate the hidden state at later time
(())= = ODESolve(Jo ho, to, ().
2)
Crucially, we can differentiate the solutions of the ODE solve with respect to the parameters 0, the initial state ho: initial time to and the solution times t. This can be done by using the adjoint scnsitivity method (Chen) let al.] 2018 Pontryagin 1987) or by differentiating through the solver operations, since solver operations themselves are differentiable This allows us to backpropagate through the ODE solve and train any learnable parameters via gradient descent. Therefore, Neural ODEs can be used as layers in network, for example, as feedforward layer where we solve the ODE from to to a terminal time f1; or they can be used to output a series by calculating the hidden state at an ordered list of times {40, (,,,,,,,,)) [n).
One advantage of Neural ODEs is that they intrinsically deal with irregular time-series, without the need for bucketing. This is because the model is able to deal with continuous time by design by modeling the instantancous velocity, rather than discontinuous update. Moreover, they use the inductive bias that time-series is generated from an underlying dynamical process, that is, the instantancous rate of change of the current state depends on the stale itself. For example, in epidemiology, compartmental models (Tang el al. 2020) are used to determine how a disease spreads through a population; or in virology the concentration of virus, antibodics and infected cells can be modelled with difforential equations (Perelson et al., 1996). In our case, the progression of MS is driven by biological mechanisms so MS progression can also be modeled with differential equations. We do not know what these differential equations are, and therefore Neural ODFs are well suited to the problem by combining the inductive bias of an underlying dynamical system with the universal approximation power of neural networks (Cybenko 1989).
2.2 Interpolation and Control Signals
In many real world cases data is incomplete, to overcome this, missing values are imputed. Often this is done by filling with zero or the mean value of a fcature. In time-series applications we are able to use already observed values to plausibly fill in missing values based on the rest of the series. For example we may have a time-series ((ro, fo), (0,4)), (12,, f2)), we can linearly interpolate, giving ((zo, (01, (1), (12, t2)), where á = ø2-20 (4 - (o) + 20. More generally, we can calculate (t), which is an interpolation of the observed To
time-series (Davis 1975 Lunardi 2009). The values of X match the observations exactly, X ((t,). = Ina and
will also predict value for unobserved features. X() can be evaluated at all times not just observation
times, giving control signal based on the observations.
3 Continuous Modeling Framework
In this work we consider continuous models. The models are adapted for this setting, via our new framework building on that of Roy et al. (2022). We recommend secing Figure 1 and Section 3 of that work for information on this stage of the pipeline. Brielly, the preprocessing consists of converting raw data to a common representation specifically for discrete models, splits subjects into train/test sets, fills missing values with zeros and pads the time-series for batching; the postprocessing is able to calculate metrics for the whole population as well as subgroups. To adapt to the continuous setting, the models receive an interpolation of the dynamic features X (the POMs at different timestamps) known the control signal and a vector of the static features known the context (such as sex and ethmicity). They then make prediction at the given evaluation times, see more about the specific models and their architectures in Section 4
Passing an interpolation to a model. Since a control signal is contimuuus, we cannot pass the whole interpolation to model, as there are infinitely many points. This is instead achieved by making the control
3

---

Published in Transactions on Machine Learning Research (09/2023)
signals piece-wise cubic functions. So that between t, and 1,+1, (t) = ((t+ t.)² +- cc (2) d2. then the coefficients (a;, bi, Ca, d:) are passed to the model to compactly represent the interpolation.
A common requirement of Deep Learning training is being able to train using batches, this is non-trivial for irregular time-series. The original high level solution (Chen et al. 2018 for continuous models to merge all of the times in the different samples (concatenate and sort). The ODE is then solved for all of the merged times, then a binary mask is applied to extract the relevant predictions for each time-series. Since then, low-level implementations of parallel ODE solvers have been introduced [Lienen & Günnemann 2022 Kidger 2021 p to the Py (Paszke et al. 2019) and Jax libraries (Bradbury ct al., 2018). These use separate ODE solver for cach scries in the batch, all being run in parallel. This requires specialist parallel hardware- GPUs, which is lypically not the most restrictive requirement. Our models and experiments were implemented in TensorFlow and Keras (Abadi et al., 2015 Chollet et al, 2015), preventing us using these low level solutions. Instead of implementing this and diverging from the main aim of the paper, we opt for different high level solution. This is a generalization of the method given in Morrill et al. (2021) from Online Neural CDEs to all continuous methods. We use the integers And instead of interpolating the time-series with the true timestamps t, we interpolate with s giving X(s). We also interpolate the true time using 8, so that t == vss). Then, if the dynamics function we learn represents the rate of change of a state with respect to true time
Batching using Pscudo-Time
regularly spaced pseudo-time,
dh
dt = fo (h,t). t).
By adapting this to the s domain such that the rate of change with respect to pseudo-tima is
dh = fo (s))
de
ds
ds
we can show that ODESolve di ho, to = ODESolve 2h fo, 50, Note this is only true if ti is an in- creasing, surjective and piece-wise differentiable function (which can be achieved with certain interpolation schemes). This allows us to straightforwardly batch as we can solve all of the ODEs in the batch in the same regular& domain, and shift each to their respective true time domains. We use this strategy over merging the times since the interpolations are already calculated for the models to run, thus we are simply re-purposing them to achieve the same result. As well as this, we do not need to solve the ODE for all merged true times, only the shared regular pseudo-times, we do not store as many intermediate ODE solution values. This suggests that this method could be more memory efficient than the merging method and compute efficient for certain models, however, this was not explicitly tested. See Appendix [ for a full explanation of the theoretical gains as well as proof of equality & requirements of o.. As well as using pseudo-time, we must pad soquences to the same length, keeping track of which obscrvations are real and which are imputed in order to train in batches. To achieve this, an initial continuous model specific data preprocessing step is used.
3.1 Data Preprocessing
A subject from the data consists of vector of context features c, and timestamped time-series of sequence features st), label v(t) & sample weights w(t). The sample weights consist of is and Os and inform the loss function and evaluation metric if a prediction is used in the calculation, if w(t) = we care about how close ō(t) is to yt) (we can set w(t)= = 0 when we have a missing label). The sequences are irregular with varying lengths, both the sequences and contexts may have missing values. We wish to convert this so that all timo-serics are the same length, with missing features imputed and the irregular times dealt with so that we can train efficiently using minibatches. Therefore the preprocessing consists of these steps in order:
1. Standardize the sequences and the context feature-wise using only the observed values so they have mcan 0 and standard deviation 1.
2. Impute missing context features with a constant. We choose zero since that is now the mean.
3. Pad sequences to the same length. We fill: sequence features with missing values (nans); observation times with the last observed time; labels with zero; sample weights with zero.
4. Calculate cumulative sum of the number of observations in each sequence feature- wise, such that if there is missing value this count does not increase. This will be called the observation cound.
4

---

Published in Transactions on Machine Learning Research (09/2023)
denoted When impute later this provides information to the model that either have true observation or an imputed one, as required by (Kidger et al. 2020).
5. Concatenate (æ, o, f) along the feature dimension, so that these are the new sequence features. Convert this to a physically plausible datasct, where future observations do not affect the past (see Section 3.2 for explanation and two methods to do this).
Interpolate the new physically plausible time-scries using integers as the pseudo-time 8, accordance with our new batching strategy. Giving piece-wise cubic coefficients.
This produces series of coefficients that represent piece-wise cubics the interpolating function for each feature. It is computationally expensive to calculate the coefficients cach time they are needed, so they are calculated once and stored. They replace (æ, o, L) in the dataset as input to the model, since the sequence is easy to recover using the coefficients. After this preprocessing each subject in the data has a tuple of: ((coeffs, c). y, w²) Note that y and w are also longitudinal just like coeffs as we make predictions at every timestep. This preprocessing means every time-scrics the dataset is the same length (due to padding) and missing values filled (either with the the dataset along the first dimension giving us minibatches, using pseudo-times allows us to calculate the predictions of batch even for irregular times, and the sample weights account for padding and missing labels in the loss calculation. Therefore, we can now train a continuous model by repeatedly giving it a minibatch, calculating the loss and using an optimizer of choice. We can also calculate the performance metrics the same batched way.
interpolated values). After preprocessing we shuflle
3.2 Enforcing Physical Solutions
It is critical that a model's predictions are physically plausible, i.e. the model only uses the current and pre- vious measurements make a prediction. At inference lime we do not. have access to future measurements, so the model must be trained such that predictions at only use the measurements up to Control signals in general are not physical, they either the whole time-scrics to interpolate, or use the next observed value which is not physical when there are missing values (i.e. we skip missing value and use future observed value). To overcome this we enforce physical solutions in continuous models in two ways:
Reprocess the dataset, by making many copics of a scries (sec Reprocessing the Dataset below). Use continuously online control signals, called "recti" control signals (see Continuously Online In- terpolation below).
Reprocessing the Dataset. For each subject in series we make copics of the subject that for the i-th timestamp, we consider the time-series up until that point, and fill forward the i-th observation until the end. Then we set the sample weights to 0 everywhere except at point i where it is 1. This means that only previous observations are used to predict at timestamp i, and then by adjusting the sample weights we only consider the i-th label for that copy in the loss. We cannot consider earlier because then the model is using observations up until make predictions for points carlier than i, we don't predict later because it is relevant do so. For example, this time-series is transformed below:
I2
Ø1
Ø2
22
t
t
t
t1
y/
y1
y1
J/2
1/2
J73
1.0
0.0
0.0
0.0
1.0
0.0
t
to
t3
V1
Y/
Ø1
Ø3
1.0
1-0
1.0
t
l/2
21
3/2
y/3
0.0
0.0
1.0
Note that in the interest of space on the page we have wrapped the repeated subjects in the transformed time-series. We have highlighted in red where the copies are created, so the time series initially has shape 3 4 and the resulting series has shape x x 4, we are extending along the batch dimension. This is

---

Published in Transactions on Machine Learning Research (09/2023)
physical for all models because the actual data itself only contains previous observations, and is then the last measurement. repeated. The downside is that this can make the dataset increase size significantly without introducing any new information.
Continuously Online Interpolation (recti method). To avoid our dataset growing exponentially we can also use continuously online interpolation schemes also called the recti method (Morrill et al.] [021]) These extend the length of the timo-scries but not the number of samples. Let be the fill-forward of z, any missing observations then construct "recti" time-series of length 2n 1. Where X(2 1) (ăz, ye wy) for i e {I,..., n) and (2i) = (3,, 1+11 Vs. 0.0) for i E (1,.... n - 1}. For example, this time-scrics is transformed below:
replaced by the most recent observed value. For time-series of length we
X22
T
ij
22
T2
t1
t2
t3
t1
t2
t2
t3
t3
!/1
Y2
1/
y1
1/1
Y/
1/2
y/
1.0
1.0
1.0
1.0
0.0
1.0
0.0
1.0
As before we have highlighted in red where the repeats occur. The original shape of this series was 1 x 3 x 4 and becomes 1 x 5 x 4 so we are extending in the time dimension rather than mumber of subjects. Crucially with the data rectified, we must then use what is known as a discretely online interpolation scheme. These only use the next observation to interpolate and since the observations are repeated this remains physically plausible. Examples include Linear interpolation and Hermite Cubic interpolation (see Appendix F). With this, the model is similarly only able to use previous observations make a prediction. This method is only valid for models which already enforce physical solutions by design, such TCN. However, later obscrvations are still included in the time-series, models that are not physical by design (i.e. those that use whole time-series to predict at each point) require the method of reprocessing the data.
3.3 Interpolation Schemes
Here we give the interpolation schemes that are used in the contimuous models. We give short overviews of cach one here, for full information including equations, and advantages & disadvantages of cach one see Appendix The existing interpolation schemes used are described in full in Morrill et al. (2021). They are:
Linear: This linearly interpolates the points with a straight line between observations, the and coefficients are zero. can be usod as an increasing, surjective, pieco-wisc differentiable interpolation so can be used interpolate the time-channel for batching with pseudo-time.
Hermite Cubic: This takes the linear interpolation, and now the previously unused and b coefficients are used to smooth the interpolation. The first derivative is continuous at observations by enforcing the derivative to be equal to the gradient between that observation and the previous.
This ODE is faster to solve than with linear interpolation due to having a continuous derivative. the interpolation, the
Natural Cubic: This is a cubic interpolation that enforces contimuity derivative and the second derivative. This makes the interpolation the most "natural" for the observed points, however, the whole time-series is used to interpolate so is not physically plausible. Rectilinear: This interpolation uses the recti method of extending along the time dimension de- scribed previously, and then linearly interpolates those points. It is continuously online control signal. However, since it is linear it can still be slow to solve the ODE.
We also introduce and invostigate two novel interpolation schemos:
Monotonic Cubic:
This takes the Hermite Cubic described previously, but rather than enforcing a derivative that would be obtained by linearly interpolating backwards, we enforce zero derivative observation points. The purpose is that it is possible to interpolate incrcasing measurements with an increasing, surjective, piece-wise wise differentiable function. This allows to interpolate the time channel to guarantee batching with pseudo-time.
Recticubic: This uses the recti system of extending the length of the time-series. In contrast to the rectilinear interpolation, we then use Hermite cubics to interpolate the features, and monotonic cu-
RRepeating the label does introduce any bias since the sample weights are repeated labels.

---

Published in Transactions on Machine Learning Research (09/2023)
bics to interpolate the observation and time channels. This constructs smooth physically plausible interpolation that also respects the requirements for interpolation of the time channel.
Note, since we are using the pseudo-time time method to train in batches, we still require an increasing, surjective, pieco-wise differentiable interpolation of = (bs). So for the linear and rectilinear schemes we use linear interpolations of the features, observations and And for all other control signals we use the monotonic cubic interpolation for the observations and & and the named scheme for interpolating the features.
4
Models
In this work we consider four neural differential equation based models. They are adapted for this setting, recommend seeing their original papers for full description, here describe them implemented, these are also shown pictorially in Figure 1 which we also recommend following alongside the descriptions. All models receive (coeffis, c) as input, where e is the context and the coeffs can be used to create a control signal of the time-series, (s). For notational purposes, a c superscript on a function denotes the function takes the context as additional input, fā t) = f(h,t, t, c, 0). This is so we can keep parameters and context nominally separate from the hidden state and time, which dynamically change during the ODE solve.
4.1 Neural Controlled Differential Equations
The dynamics a Neural CDE (Kidger et al. 2020) uses the whole control signal (not just (b))) to constantly include information from the time-series in the dynamics.
dh ds = - ((s)
dX
ds
Where fi outputs a matrix since X and its derivatives are vectors. Note that since tà is the last dimension of X we implicitly include its derivative here for batching purposes. The model embeds the first observation to an initial condition in a latent space with eg, solves the above dynamics in the latent space and predicts the labels from the hidden state at a given time with d.
ho (()))
dh ho, 0,
h. ODESolve(
ds
(hi))
4.2 ODE-RNN
The ODE-RNN model (Rubanova et al. 2019 to differential equation between observations.
generalizes RNN models the hidden state evolves according
dh
de ds
s8, (b)))
ds
And then at observations the hidden state undergoes a discontimouu change. See the top right of Figure [ for a diagram of this, as well as the equations of the method below. The initial hidden state a vector of zeros. This evolves continuously between observations and at observations changes discontinuously based on its current value, the observation and context using L6. Finally predictions are made from the hidden state at given time with d
h-1 = 0
h; ODESolVe dh
ds h4-1, 1,
h = g8 X(i) de (H.))

---

Published in Transactions on Machine Learning Research (09/2023)
X(s)
X(s)
Lover
AA ZZ
OA
(s)
10kw0))
h(s)
Lager
Neural CDE
ODE-RNN
X(s)
X(s)
ODE BNN Indies valu
denoni(k X00
h(s)
hss)
Prodetion Layer
Prelitex
Latent ODE
GRU-ODE
Figure Block diagram of the four continuous models. The interpolation of X(s) shows a lincar (dashod) and cubic (solid) interpolation. Diagrams best viewed zoomed in.
The interpolation rules mean that dw//ds zero outside the range of the time-series, the first update from i = - to 0 is only discontinuous one, since the ODE function is zero here. A technical detail is that for both the ODE-RNN and Latent ODE below, during preprocessing we create update mask for each subject at each timestep E {0, 1}, where 1 indicates we want to update the hidden state. The discontinnous update is instead h; = ( X(i)) + (1 u)h;-1- So for padded data the update zero.
4.3
Latent ODE
The Latent ODE extends the ODE-RNN model to encoder-decoder model that encodes entire sequences and then decodes them. We use ODE-RNN model to encode the whole sequence to single vector, taking the final value of the hidden state from the ODE-RNN. This is then fed into another ODE to decode, where the decoding ODE is given by
dh = [h
(b())) de
ds
The outputs of the decoding ODE are then fod through a final prediction layer dy.
ho ODERNNG (X) hidden value]
h; ODESolve vel dh
ds ho, 0,
(h,,)
When encoding the sequence with the ODE-RNN, we encode backwards such that the earliest observation has the largest effect on the initial condition of the decoding ODE, Using update masks in the backwards ODE-RNN particularly crucial here so that the padded values do not affect the encoding of the hidden state. Note that the Latent ODE encodes entire sequences, so this does not enforce physical solutions and cannot rectilinear recticubie rontrol signals.
8

---

Published in Transactions on Machine Learning Research (09/2023)
4.4
GRU-ODE
The GRU-ODE (De Brouwer et al. 2014h. The dynamics are given by
2019 is continuous analogue of the Gated Recurrent Unit CCho el al.
o (w (s) + Wrah Wroc + b,) (() (s) + W
bz) g = tanh (W g g. (s) + gh h) + Woc by) dh
(W
du
ds (1 - 2) (g h)
ds
for given weight matrices and bias vectors. These dynamics also depend on time since (bss) is the last dimension of (s). We then follow the same steps as the Neural CDE, embed with eg, solve the dynamics in the latent space and predict with dz-
ho (x)())
h. ODESolve dh ds
ho, 0,
0 (h.)
All embedding and dynamics functions are given by multilayer perceptrons with ReLU hidden activations Emnbedding layers use one hidden layer and dynamics functions use two. The embedding layers use no final non-linearity, the dynamies functions use tanh to improve stability of the hidden dynamics. All predic- tion layers are single linear layers which end in non-linearities specific to the problem: sigmoid for binary classification, softmax for multiclass classification and no non-linearily for regression.
5 Experiments
Our experiments are designed to investigate how well the continuous baseline methods can perform in the same MS setting non-continuous methods. Therefore, test the models on the same dataset from Roy ct al. 2022) and compare against the best performing model from that paper for a given task.
Typerparameters. We found hyperparameters using grid search and 10 fold cross-validation. The final hy hyprpaarmeeter configurations and values lested are given in Appendix B Using cross-validattion gives both different parameter initializations and train data splits to obtain uncertainty estimates on performance. Models are trained for epochs with the Adam optimizer (Kingma Ba 2015), learning rates and batchsizes are hyperparameters given in Appendix B
Dataset Description. The Multiple Sclerosis Outcome Assessments Consortium (MSOAC) (Rudick al. 2014h is partnership whose goal is to collect, standardize, and analyze data about MS. One such dataset which we use is the Placebo dataset, which consists of 2465 separate patient records from 9 clinical trials. Since this is public dataset, Institutional Review Board (IRB) approval is not required. The features contain demographics, medical history, MS specific data and POMs. The labels are constructed from Expanded Disability Status Scale (EDSS) scores (Kurtzke 1983). These are clinically annotated scores to determine the severity of the MS, the scores range from to 10 in steps of 0.5. EDSS scores are partitioned into four distinet subsets: 0 1 for no disability, 1.5 2.5 for mild disability, 3 4.5 for moderate disability, and 5 10 for severe disability. As in Roy et al (2022) the aim is to predict the EDSS score for patient for four given prediction horizons. Due to the irregularity of clinical visits these are split into prediction windows: O 6 months; 6 12 months; 12 18 months and 18 24 months, over the window the EDSS scores taken in clinic are averaged to EDSSSmean, if labels are stiill missing we set the sample weight w to 0 so the prediction is not used to calculate metrics or the loss. We look at four subtasks, for each one we are trying to answer a question about the EDSS score in the future prodiction window for cach subject: predicting the EDSS score
Except for GRU-ODE where the dynamics are restricted as described.

---

Published in Transactions on Machine Learning Research (09/2023)
directly a regression problem; predicting if EDSS > 3 as a binary classification problem; predicting if EDSS > 5 as a binary classification problem and predicting EDSS severity category as multiclass classification problem. These final three are clinically insightful tasks, since they indicate the category of disability as well as if there has been a significant change in the MS in a given window. When predicting EDSS directly we report root mean squared error (RMSE) and train with mean squared error. For the binary classification tasks, due to class imbalance, we report the area under the precision recall curve (AUPRC) and train with binary cross entropy. When predicting the EDSS as a severity category we calculate the AUPRC for each class in a one vs rest approach and report the mean across classes and train with the cross entropy.
Rather than using all features we use a reduced set of PASAT, SDMT, NHPT, T25FW, EDSS, Age and Sex. These correspond to Paced Auditory Serial Addition Test, Symbol Digit Modalities Test, Nine Hole Peg Test, Timed 25-Foot Walk and the current EDSS score. These are the majority of the functional tests, removing most demographics and all questionnaires from the full feature set.. We use a reduced set because interpolation and training are computationally expensive, the cost of the interpolation is (n where "f is the number of features, nt is the number of subjects and n, is the length of the series. Additionally, for every call of the Neural CDE dynamics function there is matrix-vector multiplication of cost Onny). We also initially found continuous models were unstable to many features, we carry out feature ablation in Appendix C.2 where we confirm this reduced set performs best. We specifically use these features because they have the lowest attrition compared to other features, they are approximately 10 times more abundant at the start of the study increasing to roughly 1000 times more abundant at 9 months (see Figure 2 of Roy et al. 2022) for a chart tracking the feature sparsity over time).
5.1 Results and Analysis
Overall Performance. We compare against the previously best performing model TCN, the performance has been attributed to the model's awareness of temporal evolution, in particular for predicting long term time-series. The results are given in Table [I We see that in ten out of sixteen cases the Neural CDE outperforms TCN, often quite substantially. In some cases the Neural CDE performs worse than TCN, in particular for predicting the EDSS severity category it is consistently worse. We hypothesize this is due to the features used. Whilst Neural CDE performs worse with all features, it is likely that a different subset is optimal for this task. We also see that the other continuous models are occasionally able to beat TCN, but this isn't consistent. This shows that being continuous does not guarantee performance gains, but it also does not lead to performance decrease (Neural CDE was able to beat TCN, and the other three methods typically beat all other methods from Roy et al. (2022)) ). The insight gained is that the correct model with the right hyperparameters must still be used for the problem, we have also further supported the use of Neural CDE over other continuous models. Neural CDE outperforming ODE-RNN, Latent ODE and GRU-ODE is consistent with the findings in (Kidger et al. 2020), this is likely due to Neural CDEs having provably more representation power acting on sequences than the other continuous models (see Kidger et al. (2020) for these theorems). To further support the use of Neural CDE we also newly benchmark RNN Appendix A where we see that RNN performs worse than TCN and Neural CDE.
Standardization and Interpolation Ablation. Here we investigate how the performance changes based on whether we standardize features or interpolate missing values. When we don't 't impute with interpolated values we do so with zero instead. Note that interpolation must still be done to construct control signals X(s), this ablation is about how we impute the missing values. Results for the Neural CDE are given in Table g We see that for Neural CDEs standardizing the values has the largest positive effect on performance. This is likely due to standardization transforming X into a manageable range, and by extension d /ds, that the dynamics of the Neural CDE are well regularized. We run the same ablations on the other models in Appendix C Tables [| 13 they support this hypothesis since they have far less significance between standardizing and not standardizing features and their dynamics do not end in the matrix vector multiplication dX/ds, only the vector scalar multiplication de//ds. We see that interpolating or filling with zeros does not have a significant effect; when standardization is fixed the difference between interpolating and filling with zero is almost negligible. We carry out the same ablation of the other continuous models in Appendix C and see similar results, therefore standardization is crucial for performance but interpolating missing values is not.
10

---

Published in Transactions on Machine Learning Research (09/2023)
Table 1: MSOAC Values against the best performing model from (Roy et al, 2022 All metrics are AUPRC, except for the mDSS task, where we use RMSE. Best values are bold and values better than TCN but not best underlined.
Prodiction Task
Prediction
Neural CDE
ODE RNN
Latent ODE
GRU ODE
TCN
Window
0- mo 12 mo
1.264 +: 0.055
1.408 + 0.059
1.768 + 0.139
1.617 0.072
2.312 0.271
1.650 + 0.067
1.627 = 0.055
1.886 + 0.248
1.776 0.071
2.236 0.207
EDSS
18
1.725 = 0.074
1.652 = 0. 066
1.826 + 0.123
1.810 0.090
2.066 0.086
18 24 mo
1.666 + 0.128
1.587 + 0.078
1.741 + 0.100
1.707 + 0.097
1.834 0.107
0- mo
0.909 + 0.014
0.908 + 0.010
0.869 + 0.018
0.871 0.017
0.865 0.021
EDSS
12 mo
0.820 + 0.027
0.852 + 0.016
0.823 0.024
0.819 0.030
0.788 0.025
>
18 mo
0.768 + 0.031
0.797 + 0.027
0.766 + 0.032
0.776 0.027
0.730 0.041
18 24
0.703 + 0.038
0.742 + 0. 0.035
0.697 + 0.044
0.703 + 0.044
0.666 0.062
0 mo
0.848 + 0.035
0.872 + 0.021
0.807 + 0.032
0.788 0.026
0.794 : 0.036
EDSS
12 mo
0.722 + 0.039
0.793 = 0.039
0.751 + 0.026
0.740 0.047
0.710 0.045
>
18 mo
0.669 + 0.037
0.730 + 0.031
0.678 + 0.056
0.663 0.052
0.643 0.070
18 24
0.632 + 0.037
0.658 = 0.056
0.594 + 0.047
0.602 + 0.068
0.552 0.077
0- mo
0.782 + 0.028
0.688 + 0.009
0.566 + 0.028
0.580 + 0.029
0.569 0.036
EDDSSmaa Severity
6 12 mo
0.709 + 0.044
0.672 + 0.016
0.601 + 0.024
0.596 0.021
0.555 0.038
18 mo
0.674 + 0.037
0.648 + 0.022
0.574 : 0.027
0.589 0.028
0.538 0.032
Category
18 24
0.632 + 0.037
0.616 = 0.025
0.540 + 0.018
0.557 0.031
0.513 0.055
Table 2: Ablation of the effects of standardization and interpolation of features for Neural CDE MSOAC. All metrics are AUPRC, except for the EDSSmman Lask, where we use RMSE.
Prediction Task
Prediction Window
W Standardize W Interpolate
W Standardize
WO Standardize W Interpolate
WO Standardize WO Interpolate
WO Interpolate
0 6 mo 6 12 mo
1.408 + 0.059
1.406 + 0.072
4.508 + 2.124
3.463 = 1.142
1.627 + 0.055
1.633 + 0.075
3.561 +:: 1.289
3.677 + 1.718
EDSS mean
12 18 mo
1.652 + 0.066
1.638 = 0.071
3.633 + 1.099
3.797 + 1.246
18 24 mo
1.587 + 0.078
1.578 + 0.073
3.976 + 2.341
4.047 + 1.886
0 6 mo
0.908 + 0.010
0.908 + 0.013
0.729 + 0.071
0.732 +: 0.078
EDSS mean
6 12 mo
0.852 + 0.016
0.853 = 0.022
0.662 = 0.060
0.697 = 0.033
12 18 mo
0.797 + 0.027
0.798 + 0.025
0.595 + 0.087
0.619 + 0.107
18 24 mo
0.742 + 0.035
0.742 + 0.037
0.571 + 0.117
0.564 = 0.071
0 6 mo
0.872 + 0.021
0.868 + 0.022
0.587 + 0.097
0.612 + 0.090
EDSS mean
6 12
0.793 + 0.039
0.799 + 0.036
0.560 + 0.112
0.481 + 0.114
12 18 mo
0.730 + 0.031
0.731 + 0.040
0.441 + 0.117
0.438 + 0.147
18 24 mo
0.658 + 0.056
0.658 = 0.052
0.380 = 0.126
0.397 = 0.141
0 mo
0.688 + 0.009
0.688 = 0.018
0.437 + 0.082
0.413 0.082
EDSS mean
6 12 mo
0.672 + 0.016
0.670 + 0.019
0.448 + 0.069
0.438 + 0.080
As Severity
12 18 mo
0.648 + 0.022
0.648 = 0.022
0.438 = 0.050
0.441 = 0.074
Calegory
18 24 mo
0.616 + 0.025
0.616 + 0.020
0.464 + 0.049
0.453 + 0.058
Interpolation Scheme Ablation. Here we investigate how the performance changes when using different interpolation schemes. For this we standardize the fcatures and take the best performing models for each scheme. This also tests the two methods for enforcing physical solutions the method of reprocessing the dataset is contained in the Linear, Hermite Cubic, Mono Cubic and Natural Cubic interpolations, the recti
11

---

Published in Transactions on Machine Learning Research (09/2023)
Table 3: Ablation of which interpolations work best for Neural CDE on MSOAC. All metrics are AUPRO, except for the mDSS task, where we use RMSE.
Prediction
Prediction
Natural
Hermite Cubic
Monotonic Cubic
Recti-
Recti-
Linear
Task
Window
Cubic
Linear
Cubic
0 mo
1.409 0.077
1.409 = 0.067
1.409 + 0.059
1.411 0.072
1.418 = 0.062
1.408 0.059
6 12 mo
1.627
0.055
1.633 0.091
1.644 + 0.080
1.644 0.077
1.639 + 0.062
1.651 0.100
EDSSmean
18
1.672 1 0.079
1.671 1 0.063
1.673 1 0.079
1.652 0.066
1.663 1 0.087
1.676 0.054
18 24 mo
1.602 + 0.089
1.597 = 0.067
1.587 = 0.078
1.592 0.060
1.608 = 0.090
1.689 + 0.080
0- mo
0.902 0.012
0.903 = 0.014
0.903 + 0.011
0.903 + 0.010
0.908 = 0.010
0.904 0.012
EDSSmean
6 12 mo
0.847 0.016
0.848 = 0.018
0.850 + 0.017
0.848 0.020
0.846 = 0.021
0.852 0.016
18 mo
0.794 0.027
0.795 0.027
0.797 = 0.033
0.797 0.027
0.792 = 0.024
0.795 0.024
18 24
0.733 + 0.042
0.740 + 0.042
0.737 + 0.036
0.738 0.036
0.742 + 0.035
0.738 + 0.042
0- mo
0.866 0.018
0.865 = 0.032
0.869 + 0.018
0.864 + 0.022
0.872 + 0.021
0.865 0.018
EDSSmean
12 mo
0.788 0.030
0.789 1 0.033
0.790 1 0.026
0.788 0.025
0.791 1 0.032
0.793 0.089
18 mo
0.724 0.034
0.730 0.031
0.728 + 0.027
0.727 0.044
0.726 = 0.039
0.722 0.030
18 24 mo
0.658 0.043
0.654 + 0.055
0.647 + 0.052
0.632 0.050
0.658 + 0.056
0.658 0.042
0 mo
0.686 0.014
0.683 1 0.016
0.686 1 0.012
0.686 1 0.013
0.688 1 0.009
0.686 1 0.017
EDSSmean As Severity
6 12 mo
0.663 0.018
0.669 = 0.017
0.672 = 0.016
0.664 + 0.015
0.668 = 0.016
0.666 0.025
18 mo
0.638 0.020
0.648
0.022
0.644 0.024
0.643 + 0.026
0.645 + 0.020
0.643 0.018
Category
18 24 mo
0.616 0.021
0.616 1 0.025
0.612 1 0.023
0.613 1 0.018
0.614 + 0.028
0.614 0.027
method is contained in the Rectilinear and Recticubic interpolations. The values for the Neural CDE are given in Table 3 We see that there is essentially no difference between the interpolation choices. Each one performs best on a task, and all values are within standard deviation of each other. We see the same results for ODE-RNN, Latent ODE and GRU-ODE whose results are given in Appendix [ This is consistent with the findings in Moorill et al. 2021) where the different interpolation schemes perform as well as each other since the learnable part of the Neural CDE makes up for the different interpolations.
6 Conclusion
We investigated the use of continuous time models to predict the progression of multiple sclerosis from per- formance outcome measures and demographic data. We compared four continuous time models against the previous best discrete time model temporal convolutional network (TCN) - on a public multiple sclerosis dataset (MSOAC). We saw that Neural Controlled Differential Equations are able to beat TCN on the ma- jority of experiments (ten out of sixteen), showing a well tuned continuous time model is able to outperform compete with the current best discrete time model. This thorough investigation shows adjusting the modeling choices leads to changes in experimental results and lays the foundation for further research into the application of continuous time models for MS progression. After extensive ablation studies we found that standardization is vital to achieve competitive results using the Neural CDE, whereas the exact interpolation scheme is not. To adapt the framework of Roy et al. (2022) we also introduced two new control signals and a way to efficiently batch irregular time-series via these control signals as secoudary contributions.
Limitations. We saw that three of the four continuous time models (ODE-RNN, Latent ODE, GRU- ODE) were not able to consistently outperform TCN, demonstrating modeling in continuous time does not guarantee performance gains, and that choosing the correct model still crucial step the prediction pipeline. We also saw that continuous time models required feature selection step in order to reach their peak performance, adding further resource demands to the training.
Future Work.
The main aim of this work was to demonstrate that continuous time models can be a competitive option for modeling MS progression. We achieved this on Neural CDEs, and our next steps
12

---

Published in Transactions on Machine Learning Research (09/2023)
are to contimue to test Neural CDEs more MS datasets such as Floodlight (Baker et al. 2021 | Mosaic (Roy et al. 2023). We also plan to investigate how the Neural CDE scales to more features, and making it more stable in this scenario. From an engineering standpoint, we plan to improve the interpolation speed by parallelizing and optimizing the interpolation process. Finally, in light the objective of clinicians limit progression of MS, we plan to include intervention modeling our setup.
MS
13

---

Published in Transactions on Machine Learning Research (09/2023)
Broader Impact Statement
In our work we invostigate using continuous time models to model the progression Multiple Sclerosis. Our work has shown that in some scenarios Neural rontrolled Differential Equations can outperform the current best performing discrete time model Temporal Convolutional Network for longitudinal prediction.
Applications. We view positive applications of our work, since this could be used to help practitioners treat Multiple Sclerosis and also demonstrates for future rescarch that continuous time models be effective for this problem. Since this work is in early stages it must not be deployed yet, and more research and evaluation is required first: tool of this nature needs to be tested to the same level of scrutiny as drug trials. If deployed, it must be used alongside the expert opinion of trained practitioners, and if in disagreement a further expert assessment is sought out. Crucially, in its current state a trained medical practitioner's judgement is more reliable than the Neural CDE model.
Datasets. We use a publicly available dataset in this work: Multiple Sclerosis Outcome Assessments Consortium (MSOAC) Rudick al., 2014) (https://c-path. wwwptsh..rgg//oooooeee/o wrggpproggam///ncaeeee Since the dalaset is publicly available Institutional Review Board approval is not required.
Code Release
Our code does not exist in isolation but as part of a larger code base containing proprietary code. As such our code not publicly available at this time, we plan to open-source in the future. We have included implementation details in the Appendix, including library specific caveats in Appendix and their fixes aid reproduction of our results.
Author Contributions
Lead author. Implemented the continuous models in TensorFlow. Implemented the interpolation schemes. Conceptualised and implemented the continuous model specific preprocessing. Implemented the irregular batching and wrote the proofs in the Appendix. Lead the writing of the paper and author responses.
Alexander Norcliffe:
Lev Proleev:
Ran the experiments. Provided significant. conceptual input. Provided significant assistance implementing continuous models and continuous preprocessing. Implemented the general preprocessing pipeline part of RRoy et al.] 2022). Reviewed code. Contributed to the writing of the paper.
Diana Mincu: Provided assistance implementing continuous models and contimuous preprocessing. Pro- vided significant conceptual input. Implemented the general preprocessing pipeline as part of (Roy et al. 2022). Reviewed code. Contributed to the writing of the paper and the author responses.
Fletcher Lee Hartsell:
As an expert in Multiple Sclerosis, provided very specialised input during project
meetings.
Katherine Heller:
Provided input during project meetings.
Subhrajit Roy: Senior Author. Conceptualisation of the use of continuous methods in the Multiple Sclerosis setting. Lead project to completion. Provided significant conceptual input. Provided assistance implementing continuous models and continuous preprocessing. Implemented the general preprocessing pipeline part of (Roy et al. 2022). Reviewed code. Contributed to the writing of the paper and the author responses.
Acknowledgments
We would like to thank the anonymous reviewers and action editor Qibin Zhao for their time and effort review the paper. We'd like to thank Patrick Kidger for introducing Alexander Noreliffe to Subhrajit Roy
14

---

Published in Transactions on Machine Learning Research (09/2023)
which helped to start this project. We'd like to thank James Morrill and Patrick Kidger for offering guidance the use of different control signals in Neural CDFs. We'd like to thank our colleague Mercy Asiedu for providing a detailed review of the paper before submission. We'd also like to thank our colleague Wenhui Hu for help reviewing code. We'd finally like to thank our colleagues Emily Salkey, Eltayeb Ahmed and Chintan Ghate for interesting discussion.
15

---

Published in Transactions on Machine Learning Research (09/2023)
References
Martín Abadi, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng Chen, Craig Citro, Greg S. Cor- rado, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, lan Goodfellow, Andrew Harp, GeofTrey Irving, Michael Isard, Yangging Jia, Rafal Jozelowicz, Lukasz Kaiser, Manjunath Kudlur, Josh Levenberg, Dandelion Mané, Rajat Monga, Sherry Moore, Derek Murray, Chris Olah, Mike Schuster, Jonathon Shlens, Benoit Steiner, Πlya Sutskever, Kunal Talwar, Paul Tucker, Vincent Vanhoucke, Vijay Vasudevan, Fernanda Viégas, Oriol Vinyals, Pete Warden, Martin Wattenberg, Martin Wicke, Yuan Yu, and Xiaoqiang Zheng. TensorFlow: Large-scale machine learning heterogeneous systems, 2015. URL https://www.ter www.wrwwwwtttrro....oooooooooo.. org/ Software available from tensorflow.org.
Sana Abdulkream Alharbi and Azmin Sham Rambely. A new ode-based model for tumor cells and immune system competition. Mathematics, 8(8):1285, 2020.
Shaojie Bai, Zico Kolter, and Vladlen Koltun. An empirical evaluation generic convolutional and recurrent networks for sequence modeling. arXiv preprint arXiv:1803. 01271, 2018.
Mike Baker, Johan van Beok, and Christian Gossens. Digital health: Smartphone-based monitoring of multiple sclerosis using floodlight. Nature Portfolio. URL https://www. nature. 00412-0, 2021.
James Bradbury, Roy Frostig, Peter Hawkins, Malthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman- and Qiao Zhang. JAX: compos-
able transformations Python-NumPy programs, 2018. URL http::/githee http://gggtuuuucoo./.e.e wwt.comm/ooooo comm/googe//ax
Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. Advances in neural information processing systems, 31, 2018.
Kyunghyun Cho, Bart Van Merriënboer, Dzmitry Bahdanau, and Yoshua Bengio. On the properties of neural machine translation: Encoder-decoder approaches. arXiv preprint 1rXiv:1409.1255, 2014.
François Chollet et al. Keras. https://keras: 10. 2015.
Timothy Coetzee and Alan J Thompson. Atlas of ms 2020: Informing global policy change, 2020.
George Cybenko. Approsimation by superpositions of sigmoidal function. Mathematics of control, signals and systems, 2(4):303-331 314, 1989.
Philip J Davis. Interpolation and approcimation. Courier Corporation, 1975.
Edward De Brouwer, Jaak Simm, Adam Arany, and Yves Moreau. Gru-ode-bayes Continuous modeling of sporadically-ooserveed time series. Advances in neural information processing systems, 32, 2019.
Edward De Brouwer, Thijs Becker, Yves Moreau, Eva Kubala Havrdova, Maria Trojano, Sara Eichau, Serkan Ozakbas, Marco Onofrj, Pierre Grammond, Jens Kuhle, et al. Longitudinal machine learning modeling of ms patient trajectories improves predictions of disability progression. Computer methods and programs in biomedicine, 208:106180, 2021.
John R Dormand and Peter Prince. A family of embedded rungo-kutta formulae. Journal of computational and applied mathematics, 6(1):19 26, 1980.
Emilion Dupont, Arnaud Doucet, and Yee Whye Tch. Augmented neural odes. Advances in neural informa- tion processing systems, 32, 2019.
Dachoon Gwak, Gyuhyeon Sim, Michael Poli, Stefano Massaroli, Jaegul Choo, and Edward Choi. Neural ordinary differential equations for intervention modeling. arXie preprint x=Xi:22010008944, 2020.
John J Hopfield. Neurons with graded response have collective computational properties like those of two- state neurons. Proceedings of the national academy of sciences, 81(10))3088-3000 1984.
Patrick Kidger. On Neural Differential Equations. PhD thesis, University of Oxford, 2021.
16

---

Published in Transactions on Machine Learning Research (09/2023)
Patrick Kidger. On neural differential equations. arXiv preprint arXiv:2202.02433 2022.
Patrick Kidger, James Morrill, James Foster, and Terry Lyons. Neural controlled differential equations for irregular time series. Advances in Neural Information Processing Systems, 33:6696 33:6696-6707, 2020.
Diederik P Kingma and Jimmy Ba. Adam: method for stochastic optimization. International Conference on Learning Representations, 2015.
John F Kurtzke. Rating neurologic impairment in multiple sclerosis: an expanded disability status scale (edss). Neurology, 33( 33(11):1444 1444, 1983.
Jo Lane, Huah Shin Ng, Carmel Poyser, Robyn M Lucas, and Helen Tremlett. Multiple sclerosis incidence: A systematic review of change over time by geographical region. Multiple Sclerosis and Related Disorders, 63:103932, 2022.
Michael Lawton, Kate Tilling, Neil Robertson, Helen Tremlett, Feng Zhu, Katharine Harding, Joel Oger, and Yoav Ben-Shlomo. A longitudinal model for disease progression was developed and applied to multiple sclerosis. Journal of clinical epidemiology, 68(11):1335 1365, 2015.
Emmanuelle Leray, Thibault Morcau, Agnès Fromont, and Gilles Edan. Epidemiology of multiple sclerosis. Revue neurologique, 172(1):3-15 13, 2016.
Xuechen Li, Ting-Kam Kam Leonard Wong, Ricky TQ Chen, and David K Duvenaud. Scalable gradients and variational inference for stochastic differential equations. In Symposium on Advances in Approzimate Bayesian Inference, pp. 1-28. PMLR, 2020.
Marten Lienen and Stephan Günnemann. arXiv:2210.113 12377, 2022.
torchode: A parallel ode solver for pytorch. ar preprint
Ernest Lindelöf. Sur l'application de la méthode des approximations successives aux équations différentielles ordinaires du premier ordre. Comptes rendus hebdomadaires des séances de l'Académie des sciences, 116 (3):454-457, 1894.
Alessandra Lunardi. Interpolation theory, volume 9. Edizioni della normale Pisa, 2009.
Stefano Massaroli, Michael Poli, Jinkyoo Park, Atsushi Yamashita, and Hajime Asama. Dissecting neural odes. Advances in Neural Information Processing Systems, 33:3952 33:3952-3963, 2020.
James Morrill, Patrick Kidger, Lingyi Yang, and Terry Lyons. Neural controlled differential equations for online prodiction tasks. arXiv preprint arXiv:2106. 110.11000.. 2021.
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjic Bai, and Soumith Chintala. Pytorch: An imperative style, high- performance deep learning library. In Advances Neurul Information Processing Systems 32, pp. 8024 8035. Curran Associates, Inc., 2019. URL http://papers http://paperr..nuuiss.cooo wwwr.is..c/paappaeea 2015-pytorchs aatl-pyyorrca----mmpaatttttteeee hht.appppaaaieeetty.eeeaggaeaeee
deep- www.ddeeeepppneeeee learning learning-library. pdf
Alan S Perelson, Avidan U Neumann, Martin Markowitz, John M Leonard, and David D Ho. Hiv-l dynamics in vivo: virion clearance rate, infected cell life-span, and viral generation time. Science, 271 (5255):1582 1586, 1996.
Mauro F Pinto, Hugo Oliveira, Sónia Batista, Luis Cruz, Mafalda Pinto, Inês Corrcia, Pedro Martins, and César Teixeira. Prediction of disease progression and outcomes multiple sclerosis with machine learning. Scientific reports, 10(1):1-13, 2020.
Lev Semenovich Pontryagin. Mathematical theory of optimal processes. CRC press, 1987.
17

---

Published in Transactions on Machine Learning Research (09/2023)
Zhaozhi Qian, William Zame, Lucas Fleuren, Paul Elbers, and Mihacla van der Schaar. Integrating expert odes into neural odes: pharmacology and disease progression. Advances in Neural Information. Processing Systems, 34:11364 11383, 2021.
Subhrajit Roy; Diana Mincu, Lev Proleev, Negar Rostamzadeh, Chintan Ghate, Natalie Harris, Christina Chen, Jessica Schrouff, Nenad Tomašov, Fletcher Lee Hartsell, et al. Disability prediction in multiple sclerosis using performance oulcome measures and demographic data. In Conference on Health, Inference, and Learning, pp. 375-390 PMLR, 2022.
Subhrajit Roy, Diana Mineu, Lev Proleev, Chintan Ghate, Jennifer Graves, David Steiner, Fletcher Hartsell, and Katherine Heller. Performance of machine learning models for predicting high-severity symptoms in multiple sclerosis. 2023.
Yulia Rubanova, Ricky TQ Chen, and David K Duvenaud. Latent ordinary differential equations for irregularly-sampled time series. Advances in nearal informalion processing systems, 32, 2019.
RA Rudick, G Cutter, and S Reingold. The multiple sclerosis functional composite: a new clinical outcome measure for multiple sclerosis trials. Multiple Sclerosis Journal, 8(5):359-33) 2002.
Richard A Rudick, Nicholas LaRocca, Lynn D Hudson, and Msoac. Multiple sclerosis outcome assessments consortium: genesis and initial project plan. Multiple Sclerosis Journal, 20(1):12-14 17, 2014.
Ruggiero Seccia, Daniele Gammelli, Fabio Dominici, Silvia Romano, Anna Chiara Landi, Marco Salvetti, Andrea Tacchella, Andrea Zaccaria, Andrea Crisanti, Francesca Grassi, et al. Considering patient clinical history impacts performance of machine learning models in predicting course of multiple sclerosis. PloS one, 15(3)xe0230219 2020.
Nabeel Seedat, Fergus Imrie, Alexis Bellot, Zhaozhi Qian, and Mihacla van der Schaar. Continuous-ting modeling of counterfactual outcomes using neural controlled differential equations. International Confer- ence on Machine Learning, 2022.
Lu Tang, Yiwang Zhou, Lili Wang, Soumik Purkayastha, Leyao Zhang, Jie He, Fei Wang, and Peter X-K Song. revicw of culti-compartment infectious discase models. International Statistical Review, 88(2): 462 513, 2020.
Silvia Tommasin, Sirio Cocozza, Alessandro Taloni, Costanza Gianni, Nikolaos Potsas, Giuseppe Pontillo, Maria Petracca, Serena Ruggieri, Laura De Giglio, Carlo Poszilli, al. Machine learning classifier to identify clinical and radiological fcatures relevant to disability progression in multiple sclerosis. Journal of Neurology, 268(12):4834 4845, 2021.
Ceren Tozlu, Dominique Sappey-Marinier, Gabriel Kocevar, François Cotton, Sandra Vukusic, Françoise Durand-Dubief, and Delphine Maucort-Boulch. Modeling individual disability evolution in multiple scle- rosis patients based on longitudinal multimodal imaging and clinical data. bio pp. 733295, 2019.
JN Whitaker, HF McFarland, P Rudge, and SC Reingold. Outcomes assessment in multiple sclerosis clinical trials: a eritical analysis. Multiple Sclerosis Journal, 1(1):37 1995.
Yijun Zhao, Brian C Healy, Dalia Rotstein, Charles RG Guttmann, Rohit Bakshi, Howard L Weiner, Carla F Brodley, and Tanuja Chitnis. Exploration of machine learning techniques in predicting multiple sclerosis discase course. PloS one, 12(4)0c01748 2017.
18

---

Published in Transactions on Machine Learning Research (09/2023)
A Additional RNN Baseline
Here we present results for recurrent neural network (RNN) (Hopficld 1984 on the MSOAC data, as an additional deep learning baseline deisgned for sequences not tested in Roy et al. (2022). The RNN consists of a ReLU multilayer perceptron that acts on a hidden state and observation, so that the hidden state at observation n n+1 given by hn+1 = (he, In+1), with h-1 = O. Predictions at time by are made by taking the hidden state, using a learnable prediction linear layer on it, and applying a final activation relevant the task, ((,,)) = de (h,,). We train until convergence using the Adam optimizer with batchsize of 1 (as is done for TCN in Roy et al. (0022) and learning rate 0.001. We vary the number of hidden layers in the MLP, and the width of the MLP hidden layers, which is also the width of the hidden state. Hyperparameters are found using grid search with 10 fold cross-validation (listed in Appendix B), the 10 models trained on these folds of training data are evaluated to obtain test metric mean and standard deviations. We give the results for RNN, Neural CDE and TCN in Table
Table 4: RNN and NCDE MSOAC values against the best performing model from (Roy et al. 2022 (TCN). All metrics are AUPRC, except for the EDSSmean task, where we use RMSE. Best values are in bold and values better than TCN but not best are underlined.
Prediction
Prediction
Neural CDE
TCN
RNN
Task
Window
0- mo
1.264 + 0.055
1.408 0.059
1.318 + 0.049
6 12 mo
1.650 + 0.067
1.627 0.055
1.681 + 0.068
EDSS mean
12 18 mo
1.725 = 0.074
1.652 0.066
1.801 + 0.044
18 24 mo
1.666 = 0.128
1.587 + 0.078
1.691 + 0.083
0- mo
0.909 + 0.014
0.908 + 0.010
0.898 + 0.010
EDSS mean >3 3
6 - 12 mo
0.820 + 0.027
0.852 + 0.016
0.818 + 0.023
12 18 mo
0.768 + 0.031
0.797 0.027
0.744 + 0.036
18 24 mo
0.703 + 0.038
0.742 + 0.035
0.704 + 0.030
0- mo
0.848 + 0.035
0.872 + 0.021
0.856 + 0.022
EDSS menan >5 5
6 - 12 mo
0.722 = 0.039
0.793 + 0.039
0.738 + 0.031
12 18 mo
0.669 + 0.037
0.730 0.031
0.560 + 0.023
18 24 mo
0.632 + 0.037
0.658 + 0.056
0.472 + 0.076
0- 6 mo
0.782 + 0.028
0.688 + 0.009
0.704 + 0.013
EDSS mean
6 12 mo
0.709 + 0.044
0.672 + 0.016
0.638 + 0.017
As Severity Category
12 - 18 mo
0.674 + 0.037
0.648 0.022
0.600 + 0.026
18 24 mo
0.632 + 0.037
0.616 0.025
0.577 + 0.018
We see that RNN is rarely able to beat TCN, but it does happen. However, RNN never the best performing model. To quantify this, if we take the rating across the 16 tasks (1 being best performing, being worst), the average ratings for TCN, Neural CDE and RNN are 1.83 + 0.73, 1.50 + 0.71 and 2.69 + 0.46 respectively. This shows that Neural CDE is the best performing, followed by TCN, and then RNN. Crucially it further supports the use of Neural CDEs in this context.
B Implementation and Training Details
All models are implemented using TensorFlow (Abadi et al. 2015 and Keras (Chollet et al. 2015 using the TensorFlow ODE solver https: https://www.tensorrloo
tfp/ math/ mathhode ode, Our code is not publicly available currently, however, we plan to open-source in the fulure and include full implementation details here for reproducibility purposes. All models are trained for 50 epochs with the Adam optimizer (Kingma & Bal (2015). Hyperparameters are found using grid search with 10 fold cross-validation. We found two library specific caveats with the FlosorFlow ODE solver, we describe them and provide the fixes in Appendix G.
tensorflow www.orrgggooeeee
19

---

Published in Transactions on Machine Learning Research (09/2023)
B.1 Hyperparameters
Here give the different hyperparameters of each model, the values used in the grid search and the final value found. These are given for Neural CDE in Table [ ODE-RNN in Table 6 Latent ODE in Table 7 GRU-ODE in Table [ and RNN in Table 9 In some cases the selected hyperparameter varied between the
tasks, in particular for latent widths and MLP hidden widths. We record these specific values in Table 10
Table 5: Table of hyperparameters for Neural CDE.
Hyperparameter
Candidate Values
Best Value
Embedding Layer Learning Rate
0.001 [10, 20, 40) 1
0.001
Embedding Layer Hidden Width
Varies 1
Embedding Layer No. Hidden Layers
NCDE Learning Rate
0.001
0.001
NCDE Hidden Width
[10, 20, 40] 2 0.01
Varies
NCDE No. Hidden Layers
2
Prediction Layer Learning Rate
0.01
Latent Size
[5, 10, 20]
Varies 100
Batchsize
50, 100, 200]
Table 6: Table of hyperparameters for ODE-RNN.
Hyperparameter
Candidate Values
Best Value
Discontinuous Update Learning Rate
0.001
0.001
Discontinuous Update Hidden Width
[10, 20, 40] 1
Varies 1
Discontimuous Update No. Hidden Layers
ODE Learning Rate
0.001
0.001
ODE Hidden Width
[10, 20, 40] 2 0.01
Varies
ODE No. Hidden Layers
2
Prediction Layer Learning Rate
0.01
Latent Size
[5, 10, 20]
Varies
Batchsize
[50, 100, 200]
100
Table 7: Table of hyperparameters for Latent ODE.
Hyperparameter
Candidate Values
Best Value
ODE-RNN Encoder Discontimuous Update Learning Rate
0.001
0,001
ODE-RNN Encoder Discontimuous Update Hidden Width
[10, 20, 40) 1
Varies 1
ODE-RNN Encoder Discontinuous Update No. Hidden Layers
ODE-RNN Encoder ODE Learning Rate
0.001
0.001
ODE-RNN Encoder ODE Hidden Width
[10, 20, 40) 2
Varies
ODE-RNN Encoder ODE No. Hidden Layers
2
ODE Decoder Learning Rate
0.001
0.001
ODE Decoder Hidden Width
[10, 20, 40] 2
Varies
ODE Decoder No. Hidden Layers
2
Prediction Layer Learning Rate
0.01
0.01
Latent Size
[5, 10, 20]
Varies
Batchsize
[50, 100, 200]
100
20

---

Published in Transactions on Machine Learning Research (09/2023)
Table 8: Table of hyperparameters for GRU-ODE.
Hyperparameter
Candidate Values
Best Value
Embedding Layer Learning Rate
0.001
0.001
Embedding Layer Hidden Width
[10, 20, 40)
Varies
Embedding Layer No. Hidden Layers
2 0.001
2
GRU-ODE Learning Rate
0.001
Prediction Layer Learning Rate
0.01
0.01
Latent Size
[5, 10, 20]
Varies
Batchsize
[50, 100, 200]
100
Table 9: Table of hyperparameters for the RNN.
Hyperparameter
Candidate Values
Best Value
Hidden Width
[32, 64]
Varies
No. Hidden Layers
[2, 3]
Varies
Table 10: MSOAC Hyperparameters that varied between tasks, For the continuous models, the first value is the selected Latent Width and the second value is the selected Hidden Width of the models. For the RNN, the first value is the hidden width, the second value is the number of hidden layers.
Prediction Task
Prediction Window
Neural
ODE
Latent ODE
GRU
RNN
CDE
RNN
ODE
0- 6 mo
40,
20
40,
10
40,
5
40,
64,
2
6 12 mo
40,
20
40,
5
40,
5
40,
20
64,
2
EDSS mean
12 18 mo
40, 20
40,
5
40,
5
40,
20
64,
2
18 24 mo
40, 10
40,
5
40,
10
40,
10
64,
2
0- 6 mo
40, 20
20,
10
10,
20
40,
20
32,
2
EDSS mean
6 12 mo
20, 20
20,
10
40,
5
40,
20
64,
2
>
12 18 mo
40, 20
20,
5
10.
20
40,
20
64,
2
18 24 mo
20, 20
40,
5
10,
20
40,
20
32,
2
0- 6 mo
40, 20
10,
10
10,
10
40,
20
32,
3
E.SS mean
12 mo
40,
20
40,
5
20,
20
40,
20
64,
2
> 5
12 18 mo
40,
20
40,
5
10,
20
40,
20
64,
2
18 24 mo
40,
20
10,
20
10,
20
40,
20
64,
2
6 mo
40, 20
40,
40,
10
40,
20
64,
2
EDSS
6 12 mo
40,
20
40,
5
40.
5
40,
20
64,
2
As Severity
12 18 mo
40,
20
40,
5
10.
20
40,
20
64,
2
Category
18 24 mo
40, 20
20,
10
10.
20
40,
20
64,
2
21

---

Published in Transactions on Machine Learning Research (09/2023)
C Ablation Studies
In Scction we carried out ablations investigating how standardization, interpolating missing values and interpolation scheme affect the performance of Neural CDEs. In this section we do the same for ODE-RNN, Latent ODE and GRU-ODE.
C.1
Standardizing and Interpolating
We look at the effects of standardizing features and interpolating missing features for ODE-RNN in Table 11 Latent ODE in Table T2 and GRU-ODE in Table [3] We see that on the whole ODE-RNN does best when standardizing features and filling missing values with zeros, Latent ODE does not have a clear superior scheme and GRU-ODE does best when standardizing features and filling missing values with interpolated values. However for these three methods the differences are small, whereas with Neural CDEs there was clear distinction between standardizing and not standardizing features. This is likely due to the models' dynamics functions. The dynamics function for the Neural CDE ends with dX/dss, whereas all others end with dø/ds. In MSOAC the true times are fairly small and therefore so is dy//ds. However, X can be quite large, and therefore so is d'//ds. And so standardizing X aids Neural CDEs more because d'/dd in the dynamics is explicitly shifted to a more manageable range, whereas d /ds was already in that range.
Table 11: Ablation of the effects of standardization and interpolation of features for ODE-RNN on MSOAC. All metrics are AUPRC, except for the EDSSmmen task, where we use RMSE.
Prediction
Prediction Window
Standardize W Interpolate
W Standardize WO Interpolate
WO Standardize W Interpolate
WO Standardize WO Interpolate
Task
0 6 mo
1.768 + 0.139
1.743 + 0.145
1.725 + 0.167
1.756 + 0.138
12 mo
1.886 0.248
1.870 = 0,103
1.842 + 0.167
1.821 = 0.163
EDSS
18 mo
1.826 0.123
1.802 + 0.075
1.825 + 0.117
1.827 0.086
24 mo
1.741 0.100
1.709 = 0.101
1.717 = 0.066
1.727 = 0.100
0 6 mo
0.869 + 0.018
0.873 = 0.009
0.866 + 0.026
0.864 = 0.025
EDSS >
6 12 mo
0.823 + 0.024
0.825 + 0.029
0.815 + 0.030
0.816 = 0.027
12 18 mo
0.766 0.032
0.762 + 0.034
0.758 + 0.051
0.768 = 0.023
24 mo
0.697 0.044
0.695 = 0.042
0.701 0.038
0.698 0.032
0 6 mo
0.807 + 0.032
0.818 + 0.027
0.792 + 0.054
0.806 + 0.027
EDSS mean > 5
6 12 mo
0.751 + 0.026
0.757 1 0.040
0.726 + 0.028
0.728 + 0.049
12 18 mo
0.678 = 0.056
0.682 = 0.035
0.655 = 0.029
0.662 = 0.050
24 mo
0.594 0.047
0.600 + 0.044
0.588 + 0.053
0.593 0.055
0- 6 mo
0.566 + 0.028
0.578 + 0.024
0.564 + 0.025
0.564 + 0.038
EDSSmean As Severity Category
6 12 mo
0.601 + 0.024
0.593 1 0.036
0.588 + 0.018
0.583 + 0.033
12 18 mo
0.574 0.027
0.576 = 0.026
0.577 = 0.027
0.577 = 0.026
24 mo
0.540 0.018
0.545 + 0.034
0.542 + 0.028
0.542 = 0.022
22

---

Published in Transactions on Machine Learning Research (09/2023)
Table 12: Ablation of the effects of standardization and interpolation of features for Latent ODE on MSOAC. All metrics are AUPRC, except for the EDSS mean task, where we use RMSE.
Prediction Task
Prediction
W Standardize W Interpolate
W Standardize WO Interpolate
WO Standardize W Interpolate
WO Standardize WO Interpolate
Window
0 mo
1.617 0.072
1.615 + 0.088
1.585 + 0.079
1.596 0.094
12 mo
1.776 0.071
1.784 = 0.081
1.754 0.099
1.775 = 0.073
EDSS mean
18 mo
1.810 0.090
1.763 + 0.094
1.735 + 0.100
1.758 0.084
18 24 mo
1.707 + 0.097
1.683 = 0.094
1.644 + 0.090
1.653 0.074
0 6 mo
0.871 + 0.017
0.879 + 0.017
0.866 + 0.025
0.867 + 0.018
EDSS
12 mo
0.819 0.030
0.824 = 0.026
0.818 0.034
0.819 = 0.020
12 18 mo
0.776 0.027
0.770 = 0.040
0.777 0.037
0.777 = 0.024
24 mo
0.703 0.044
0.713 = 0.044
0.711 0.032
0.712 = 0.035
0 6 mo
0.788 + 0.026
0.786 = 0.031
0.786 = 0.025
0.775 = 0.037
EDSS >
6 12 mo
0.740 + 0.047
0.737 + 0.035
0.732 + 0.039
0.741 + 0.037
12 18 mo
0.663 0.052
0.665 = 0.049
0.676 = 0.040
0.679 = 0.035
24 mo
0.602 0.068
0.613 = 0.064
0.602 = 0.058
0.604 = 0.063
0 6 mo
0.580 + 0.029
0.590 + 0.020
0.585 + 0.018
0.591 + 0.034
EDSSmeat As Severity
6 12 mo
0.596 + 0.021
0.599 = 0.018
0.600 = 0.032
0.609 = 0.018
12 18 mo
0.589 0.028
0.593 = 0.026
0.599 = 0.026
0.610 = 0.024
Category
24 mo
0.557 0.031
0.563 + 0.034
0.573 + 0.025
0.571 0.024
Table 13: Ablation of the effects of slandardization and interpolation of features for GRU-ODE on MSOAC. All metrics are AUPRC, except for the EDSS task, where we use RMSE.
Prediction Task
Prediction
W Standardize W Interpolate
W Standardize WO Interpolate
WO Standardize
WO Standardize WO Interpolate
Window
W Interpolate
0 mo
2.312 0.271
2.295 + 0.362
2.115 0.374
2.043 0.188
12 mo
2.236 0.207
2.252 = 0,183
2.044 0.135
1.977 = 0.119
EDSS mean
18 mo
2.066 0.086
2.066 + 0.138
2.009 + 0.112
2.001 0.096
18 24 mo
1.834 + 0.107
1.842 + 0.116
1.880 + 0.069
1.908 + 0.097
0 6 mo
0.865 +: 0.021
0.867 + 0.019
0.848 + 0.024
0.851 + 0.029
EDSS menn
12 mo
0.788 0.025
0.791 = 0.023
0.780 0.033
0.753 = 0.026
12 18 mo
0.730 0.041
0.727 + 0.047
0.657 + 0.040
0.659 0.062
24 mo
0.666 0.062
0.650 = 0.039
0.577 = 0.089
0.566 0.074
0 6 mo
0.794 + 0.036
0.805 + 0.031
0.756 + 0.035
0.761 + 0.031
EDSS >
6 12 mo
0.710 + 0.045
0.698 + 0.042
0.671 + 0.061
0.640 = 0.037
12 18 mo
0.643 0.070
0.620 + 0.060
0.531 + 0.132
0.559 = 0.081
24 mo
0.552 0.077
0.524 + 0.093
0.425 + 0.164
0.393 0.091
0- 6 mo
0.569 + 0.036
0.558 + 0.032
0.510 + 0.021
0.513 + 0.036
EDSSmean As Severity
6 12 mo
0.555 + 0.038
0.545 1 0.029
0.492 + 0.033
0.500 + 0.039
12 18 mo
0.538 0.032
0.548 = 0.038
0.496 - 0.043
0.496 = 0.028
Category
24 mo
0.513 0.055
0.511 + 0.038
0.442 + 0.033
0.440 0.040
23

---

Published in Transactions on Machine Learning Research (09/2023)
C.2 Interpolation Schemes
We investigate which interpolation schemes are best for ODE-RNN in Table 14 Latent ODE in Table 15 and GRU-ODE in Table I6 As was the case for Neural CDEs, we see that there is minimal difference between the choices interpolation scheme.
Table 14: Ablation of which interpolations work best for ODE-RNN on MSOAC. All metrics are AUPRC, except for the EDSS task, where we use RMSE.
Prediction
Prediction
Natural
Hermite
Monotomic Cubic
Recti- Linear
Recti-
Linear
Task
Window
Cubic
Cubic
Cubic
0 mo
1.805 + 0.186
1.820 = 0.232
1.939 0.290
1.975 0.230
1.889 + 0.294
1.768 0.139
6 mo
1.898 1 0.204
1.807 1 0.099
1.917 0.196
2.012 0.215
1.915 1 0.186
1.886 0.248
EDSSmean
12
1.826 + 0.123
1.896 0.174
1.875 0.147
1.855 0.128
1.852 + 0.136
1.835 0.143
18 24
1.741 + 0.100
1.754 + 0.108
1.766 + 0.103
1.801 + 0.136
1.833 + 0.175
1.757 + 0.136
0 mo
0.863 + 0.020
0.868 + 0.024
0.864 0.024
0.865 = 0.016
0.863 + 0.017
0.869 + 0.018
EDSSmean
mo
0.025
0.810 1 0.034
0.823 0.024
0.820 1 0.027
0.029
0.814 0.026
>
mo
0.028
0.760 +: 0.042
0.756 0.027
0.763
0.766
0.032
0.732 0.026
18
0.695 + 0.032
0.692 +: 0.036
0.673 0.038
0.694 0.039
0.692 + 0.034
0.697 0.044
0- mo
0.803 1 0.030
0.804 1 0.039
0.807 0.032
0.807 1 0.027
0.802 TT 0.028
0.805 1 0.020
EDSSmean
6 12 mo
0.742 = 0.047
0.732 = 0.056
0.751 0.026
0.736 = 0.031
0.730 + 0.060
0.744 + 0.043
>
mo
0.678 0.056
0.653 0.041
0.672 0.053
0.662
0.675 + 0.040
0.663 0.042
18
0.386 + 0.046
0.592 0.045
0.594 0.047
0.590 0.050
0.581 + 0.049
0.572
0- 6 mo
0.566 + 0.028
0.562 # 0.036
0.556 + 0.040
0.559 = 0.030
0.558 + 0.024
0.537 + 0.028
EDSSmean As Severity
6 12 mo
0.591 + 0.041
0.579 + 0.034
0.580 0.024
0.583 + 0.046
0.601 t 0.024
0.586 + 0.021
12 mo
0.359 + 0.019
0.560 + 0.042
0.559
0.025
0.574 0.027
0.570 1 0.020
0.573 1 0.025
Category
18
0.339 + 0.026
0.538 0.029
0.530 0.029
0.538 0.030
0.534 + 0.032
Table 15: Ablation of which interpolations work best for Latent ODE on MSOAC. All metrics are AUPRC, except for the EDSSmean task, where we use RMSE. Recall that since Latent ODE encodes a whole sequence and then decodes to predict, the rectilincar and recticubic control signals cannot be used.
Prediction Task
Prediction Window
Natural Cubic
Hermite Cubic
Monotonic Cubic
Linear
0 - 6 mo
1.619 + 0.078
1.620 + 0.095
1.643 + 0.080
1.617 0.072
6 - 12 mo
1.806 + 0.077
1.776 + 0.071
1.819 + 0.132
1.785 + 0.058
EDSS mean
12 mo
1.816 0.116
1.810 + 0.090
1.810 0.089
1.838 0.151
18 24 mo
1.726 + 0.081
1.717 + 0.094
1.747 0.115
1.707 + 0.097
0 - 6 mo
0.863 + 0.024
0.866 + 0.017
0.868 + 0.023
0.871 0.017
EDSS > 3
6 - 12 mo
0.818 + 0.029
0.817 4: 0.019
0.819 + 0.026
0.819 +: 0.030
12 mo
0.776 = 0.027
0.769 + 0.022
0.765 = 0.038
0.763 - 0.033
18 24 mo
0.699 0.054
0.700 + 0.047
0.692 + 0.045
0.703 + 0.044
0 6 mo
0.779 + 0.044
0.788 + 0.026
0.780 + 0.057
0.779 + 0.040
EDSS > 5
6 - 12 mo
0.733 + 0.039
0.719 + 0.037
0.718 + 0.046
0.740 +- 0.047
12 - 18 mo
0.654 + 0.038
0.651 + 0.048
0.656 + 0.044
0.663 + 0.052
18 mo
0.581 + 0.050
0.574 4 0.061
0.590 + 0.047
0.602 + 0.068
0 6 mo
0.580 + 0.029
0.568 + 0.035
0.573 + 0.030
0.559 + 0.024
EDSS As Severity Category
6 - 12 mo
0.595 + 0.026
0.587 + 0.023
0.596 + 0.021
0.589 + 0.031
12 18 mo
0.587 + 0.030
0.589 + 0.028
0.586 + 0.033
0.585 - 0.031
18 mo
0.557 + 0.021
0.555 + 0.033
0.557 + 0.031
0.556 + 0.033
24

---

Published in Transactions on Machine Learning Research (09/2023)
Table 16: Ablation of which interpolations work best for GRU-ODE on MSOAC. All metrics are AUPRO, except for the mDSS task, where we use RMSE.
Prediction
Prediction
Natural
Hermite
Monotonic Cubic
Recti-
Recti-
Linear
Task
Window
Cubic
Cubic
Lincar
Cubic
0 mo 6 12 mo
2.312 0.271
2,421 = 0.199
2.570 + 0.393
2.381
0.445
2,434 = 0.430
2.512 0.630
2.335 0.342
2.275 0.167
2.236 0.207
2.290 0.247
2.324 + 0.379
2.468 0.313
EDSSmean
18
2.077 0.183
2.066
0.086
2.114 1 0.130
2.148
0.229
2.126 1 0.161
2:241 0.287
18 24 mo
1.863 + 0.105
1.834 = 0.107
1.870 = 0.124
1.863 0.094
1.884 = 0.075
1.914 + 0.127
0- mo
0.854 0.081
0,860 = 0.012
0.855 + 0.038
0.855 0.030
0.865 = 0.021
0.847 0.027
EDSSmean
6 12 mo
0.785 0.029
0.773 0.038
0.776 + 0.047
0.787 0.027
0.774 = 0.036
0.788 0.025
18
0.694
0.062
0.730 0.041
0.704 = 0.096
0.723 0.049
0.723 = 0.042
0.663 0.060
18 24
0.666 + 0.062
0.632 + 0.112
0.631 + 0.071
0.559 163
0.661 + 0.062
0.623 0.100
0- mo
0.773 0.052
0.794 + 0.036
0.782 + 0.029
0.761 + 0.041
0.790 = 0.043
0.778 0.052
EDSSmean
12 mo
0.673 0.061
0.702 0.054
0.710 1 0.045
0.695 0.034
0.649 0.063
0.038
18
0.630 0.063
0.624 0.093
0.623 + 0.082
0.643 = 0.070
0.070
0.639 = 0.046
0.571 0.137
24 mo
0.511 0.060
0.463 + 0.184
0.523 + 0.069
0.552 0.077
0.476 + 0.117
0.514 0.117
0 mo
0.552 0.037
0.552 1 0.021
0.367 1 0.036
0.560 T. 0.026
0.569 1 0. 0.036
0.559 0.031
EDSSmean As Severity
6 12 mo
0.541 0.037
0.544 = 0.028
0.531 = 0.050
0.541 0.033
0.540 = 0.037
0.555 0.038
18 mo
0.538
0.498 0.036
0.504 0.062
0.534 0.032
0.513 + 0.059
0.509 0.058
Category
18 24
0.513 0.055
0.464 1 0.057
0.484 1 0.054
0.496 0.069
0.509 1 0.040
0.479 0.049
C.3 Feature Ablation
In Section 5] it was stated that used reduced feature set for the continuous models. We use Paced Audi- tory Serial Addition Test, Symbol Digit Modalities Test, Nine Hole Peg Test, Timed 2o-Foot Walk, current EDSS score, Age and Sex. We focus on most of the functional tests, removing most of the demographics and all questionnaires. In order to evaluate these features against the full feature set, we test Neural CDE on the full feature set and the reduced set in Table we see that Neural CDE is consistently better on the roduced feature set, and significantly so.
Table 17: Feature Ablation Neural CDE. *; Not evaluated.
Prediction Task
Prediction Window
Neural CDE
Neural CDE
Reduced Features
All Features
0 mo
0.908 : 0.010
0.756 + 0.113
EDSSmean
6 12 mo
0.852 + 0.016
0.711 + 0.089
12 18
0.797 + 0.027
0.695 + 0.058
18 24 mo
0.742 + 0.035
0.641 + 0.068
0 mo
0.872 + 0.021
0.609 + 0.085
EDSS menn
6 12 mo
0.793 + 0.039
0.517 + 0.112
12 18 mo
0.730 + 0.031
0.455 + 0.116
18 24 mo
0.658 + 0.056
0.445 + 0.122
0 mo
0.688 + 0.009
0.487 + 0.067
EDSS As Severity
6 12 mo
0.672 : 0.016
0.517 + 0.065
12 18 mo
0.648 + 0.022
0.482 + 0.071
Category
18 24
0.616 + 0.025
25

---

Published in Transactions on Machine Learning Research (09/2023)
We further hypothesize that one of the reasons NCDE becomes unstable to large amount of features is because the dynamics function is given by matrix-vector multiplication. The dynamics of the NCDF is ds which is a matrix-vector multiplication, whereas the other three models have de (h, (())) d which is vector-scalar multiplication. This would mean, since all functions fÿ have final tanh activation, the NCDE dynamics are bounded by = =Σ; =E,d, dXx/ds, whereas the other models dynamics are bounded by =db//ds. This explains why NCDE is unstable to a large amount of features, because the magnitude of the velocity of the hidden state grows with the number of features. One possible solution would be to divide the NCDE dynamics by the dimensionality of the vector X, which would make the dynamiis agnostic to the number of features. It further explains the importance of standardizing X, keeping dX/ds ds in a manageable range to avoid numerical instability in the NCDE trajectories.
given by dh zs = [h, (s)) dx
dynamics function given by dh
D Batching Irregular Time-Series
In Section we introduced the pseudo-time method for training continuous models using batches. This relies the equality
ODESolve(f, ho, to t) = ODESolve(f. ho, s),
where f(h, fh, (s)) de and that we interpolate t with an increasing, surjective, pieco- diffcrentiable de
function so that t (s). We prove this below.
We carry out an ODE solve over an infinitesimal pseudo-time step from Su to so + As. We have to = ()())), deh
to + Al = (s0 + As). For small enough As we have A = z As. We solve the ODE for the initial condition h(s))= ho, and dh ds
f(h, s). Carrying out the ODE solve over the infinitesimal step we have:
h(so + As) h(su) + (oou))+(((((s)))) f(h(so), so)As
-h(so) + ((h(80), (50)) du ds
As
ho f(oo, tt)
We also carry out a second ODE solve from to to to + At, where the initial condition is h(to) hg and dh t), Carrying out the ODE solve over this infinitesimal step have:
hito + At) = ho + f(ho: At
Giving us h(to + At) = h(so + As). ODEs are fully determined by their dynamics and initial condition so we can repeat the above step making h(so + As) the new initial condition. This can be repeated for the whole s domain, completing the proof. We demonstrate the batching procedure pictorially in Figure 2
3
From the proof we require in the limit that At = de As to carry out the infinitesimal solve, since At 2 O and As 2 0, we require ds to both exist (differentiable) and be 2 0, so we need an increasing function. As well as this we need to be able to carry out this step anywhere on the interval so &i must be surjective. Finally, we do not need to be globally differentiable, since for whole ODE solve can split the interval into non-overlapping sub-intervals and vo the ODE sequentially on those. Therefore & only needs to be piece-wise differentiable. Examples of interpolations that fulfill this criteria are the linear and mono cubic interpolations.
why we need vi to be increasing, surjective and piece-wis differentiable. We de
del
Computational and Memory Requirements. We use meeudo-times to achieve batching as a simple solution that software agnostic and reuses already calculated interpolations from the Neural CDE. However, there is at least some theoretical evidence to suggest it would be more memory efficient than the method of taking unions of time-series and more compute efficient specifically for Neural CDEs. This is also supported by literature.
On page 60 of (Kidger 2021), Section 3.2.1.3, Kidger remarks while discussing the scheme of Morrill et al. h021h for Neural CDEs: "Remark 3.12. This was actually a mistake we made in 1Kid+20a) the above
The Picard-Lindeed theorem requires that the dynamics function be Lipschitz continuous Lindeläf 1894.. This achieved with ReLU activations.
26

---

Published in Transactions on Machine Learning Research (09/2023)
5
Figure 2: Visualization of the reparameterization from s to t. In the top left is a regular time-series where represents pseudo-time (this panel is rotated 90 anti-clockwise). In the top right is how we map from to = (s)), in the bottom right have the irregular time-series where & is true time. When solve include de d; in the dynamics function, which is effectively using the top right figure as way to transform from s to t. Each time-series has a different version of this since the time-series are irregular, but these are all calculated before training so can be included easily for efficient batching.
procedure not done, in favour of alternate (storage-inefficient) scheme that involved taking the union over the times 1; needed to handle each batch of data." That is, instead of the method that propose (which is generalization of (Morrill et al. 2021)), ) the times for each series in the batch are merged & sorted, then the combined system is solved for all times and only certain outputs are used in the loss calculation using a binary mask. So the existing literature supports this method over the union method of batching.
We can also show this by considering the complexities of the two methods. In the following, assume that all time-series in the dataset are length Re, the analysis can be straightforwardly adapted when this is not the case. Assume that all initial times are the same, which is requirement of merging the time-series to train in batches. We also assume that no two time-series share any evaluation times, if they do then these are combined when merging the sorted lists.
We start by considering the method of taking unions. First the union must be taken over sets of times in the batch and sorted. Then a binary mask created to indicate whether an output used the loss or metric calculation. This has memory cost ((frn)) where n, is the batchsize, i.c. for each series the batch (size n4) need a binary mask of length n&Mt, the total number of merged times. Further to this, there is further memory cost during the ODE solve, here we calculate the hidden state at nyn, times, which must be stored (for loss calculation) for each hidden state in the batch giving further memory requirement of On)rma)) where ni, is the size of the hidden state. Combining these two costs we have On + 1)) memory cost.
Now we consider the computational cost. Assuming that during training there shuffling of the series create new batches each epoch, then we must merge and sort the time-series at every iteration. Merging the sorted lists is cost (nnn))) there are nn/ng iterations per epoch where n, is the size of the train set,,
27

---

Published in Transactions on Machine Learning Research (09/2023)
merging and sorting has computational cost me) over the whole training, where ne is the number of epochs. For each solve we calculate the output of the ODE at 1674 times, which means during the backward solve, we have to stop the solve explicitly this many times to adjust the adjoint state (which undergoes discrete jumps at measurement points, see (Chen et al.) 2018) for an explanation of this). There is a cost of restarting an ODE solve (calculating the initial step size etc.) which if we treat as constant gives us (ngm)) additional computational cost. This is done for every iteration of every epoch which gives a cost of On..igne) So combining with the cost of sorting we have additional (2 to train.
Now we can compare to the pseudo-time method starting with the memory complexity. We need to now store cubic coefficients of the interpolation. This gives four coefficients for cach point in each scries giving memory complexity (4n,m))) If we consider the memory requirement of the solve, we need to store the hidden state at all times during the solve which is reduced to (Igneenn) Adding these two gives the total memory cost ((n]n,(An,/ng) + ma/m))). Comparing to the union method we directly compare n2 + 1 for the union method to 4n./n] + na/nb: where we have disrogarded the common factor of nīn,, Whilst this is dependent on specific values of nt; ny and Rs; we can see that as the batchsize increases the pseudo-time method will become more memory efficient than the union method (as supported by (Kidger 2021)).
Next we look at the computational complexity. As given in Section 5 to calculate the interpolation vb we have to interpolate each time-series individually which gives cost ((nyn,m)) where ny is the number of features. This is a single preprocessing step so does not need to be done at each training iteration, so we do not consider it here. If we consider the computational cost during training, we only need to solve the ODE for n, times so the computational cost of restarting the ODE using the adjoint method for every iteration of
every epoch is (n..m, Finally we have to calculate (s(s) and de at every function evaluation during
dh = fh, ((s) ds As described in Appendix E this involves taking a floor operation, the solve to calculate de and calculating a cubic. Which is constant, but must be done at each function evaluation. This is done for every sample in the batch (size n³) every iteration of training giving ((n.n.m,,) where nn is the number of dynamies function evaluations during the solve. Combining these we have cost Omnni,me/n no + n.n,mn)) Comparing this to the union method which had cost O(2nn we depend on nn vs n4. Using advanced solvers n, is independent of n4, so the computational cost is inconclusive and depends on nh, i.e. how complex the hidden dynamics are. A significant exception is for Neural CDEs, the best performing model we benchmarked. Here the computational cost for the pseudo-time method is unchanged. However, for the union method the dynamics include 4x interpolate the features). Now we include the computational cost. associated with calculating ds or de at every step. We also include the computational cost of calculating the interpolation coefficients. So now the union method involves all the computational costs that the pseudo-time method does, it also includes a further O ny,)) cost associated with merging sets of sorted times and restarting the backward solve at more times. So theoretically we can see that for Neural CDEs the pseudo-time method is computationally more advantageous than the union method.
see that the computational costs
dr (where the true time is used to dX dx
So we have theoretical evidence to suggest the pseudo-time method is more memory efficient than the union method, but the computational complexity is problem dependent (how many function evaluations are needed during the ODE solve), except for Neural CDEs where the pscudo-time mothod is theoretically faster. Again note that this has not been experimentally verified but is supported by Kidger (2021).
E Calculating Interpolations
In this section describe how a general control signal is calculated without using any specific interpolation scheme and in Appendix F we discuss the interpolation schemes themselves. We use piece-wise cubics to interpolate, we don't use polynomials of arbitrary order since they can oscillate significantly.
Assuming we use integers as pseudo-time to interpolate, given a point s to evaluate at, a control signal will first find the largest integer which is smaller than 8 (floor function), and then X(s) can be calculated using
((s) as -i)3 + be(s 1)2 + ca(s -i) + d;.
(3)
28

---

Published in Transactions on Machine Learning Research (09/2023)
And its derivative can be calculated using
dX 3a,(s - 1)2 2b,(s 1) 6s
(4)
It is possible to expand these brackets, and then we calculate new coefficients that
X()) ā,s + b,s² + 2s + dj.
However, we found that this can be numerically unstable thus we use the method given by Equation 3 With this, the interpolation step finds the coefficients 0,, bi, G and d,, which change depending on the interpolation scheme.
To deal with missingness, consider one feature of partially observed time-series with n measurements: ((ro, 0), (0, (), 2), (0, Here are missing measurements of this feature at and 3. 3. So remove these:
3), (Ii 2:: 2, 2), (2 1))
(( 0): (x2, 2), ( 2), 1,, 11)
We then run the interpolation on what is left getting b, coefficients for these observations. Finally to fill in the missing observations, we must continue the previous cubic piece. For example; between s = 3 and = 4, we would like to continue the cubic that is between polynomials to match at all s:
and 3. To do this we require the two
a,(s - 1)3 (e(s 2)2 cs i) + d; as41 (i + 1)) b.+ (i 1))² c (i+ d++1 Solving this, to find the coefficients where have missing value take linear combination of the previous coefficients
a:+1
Q4
bi+1
b;
C++1
0
C₁
de+1
[
d3
For -1 n, ie. outside the range of the timo-scries, set 5 and d-1 do d, dn-1- That is, outside the range of the time-series the control signal predicts either the first or the last observed values as constant, with zero derivative.
Note, since we have no knowledge about the missingness of features, or the irregularity of the time-series we must carry out this procedure separately for each feature each time-series. This slow process and therefore part of data preprocessing and the coofficients are stored.
F Interpolation Schemes
Here we give information about the specific interpolation schemes used, including the positives and negatives of using each. The majority of this is summarised from (Morrill et al. 2021). Importantly, for the t channel we must use a control signal that is increasing, surjective and piece-wise wise differentiable. We would also like to use an increasing interpolation for the observation counts since these can only increase and we can enforce this inductive bias. In the following we give the steps for finding the coefficients to interpolate, this assumes general set of times, rather than the integer pseudo-times we use. Therefore these can be used generally, and the specific case of the integers is trivial change.
F.1 Linear
This is the simplest interpolation scheme, we carry out linear interpolation between two observed points,
a
0
b;
0
=
(I1+1 - x4 2x)/ tx)
dg
29

---

Published in Transactions on Machine Learning Research (09/2023)
Pros and Cons. This is the simplest interpolation scheme, and intuitively the casiest to understand. If the observations are increasing then the control signal is an increasing function, therefore this can be used to interpolate time and observation counts. However, it does not have a continuous first derivative so is slow to use in an ODE solver, also misreprcsents reality, it unlikely features change in such a way that the derivative discontinuous at points.
F.2 Hermite Cubic
Hermite interpolation is a smooth version of linear interpolation, where rather than discontinuities in the derivative at observations these are smoothed by making the derivative continuous there. Let the derivative at each point be equal to the derivative that would be calculated by linearly interpolating with the previous observation, m, = (x: - 11 - 14-1) for i 2 1. Then set mg = m.. With this we can calculate the coofficients
a
(-2(x;;1 - x;)) + (ti+1 - (:)(mm - m m241))((t-111) (3(0141 - xii - (li+1 - 44)mmm - m)+ (14+1))((444+ 1,)²
C;
my
d
24
Pros and Cons. This is smoothed out version of the linear interpolation, is faster to use in solve. Additionally, it is smooth it is a better inductive bias for interpolating features. However, it cannot be used to interpolate time or observation count since it is not increasing, Also whilst having a continuous first derivative, does not have a continuous second derivative.
F.3 Monotonic Cubic
We desire scheme with continuous first derivative to interpolate time and observation count. However, Hermite cubics are not increasing functions and picce-wise lincar does not have continuous derivatives. To overcome this we propose a new interpolation scheme, the monotonie cubic (mono cubic). This is a Hermite cubic but we enforce zero derivative at every observation, therefore it is an increasing function if the observations are always increasing. has continuous first derivative so we can use it place of linear and it is increasing and surjective so we can use in place of the Hermite interpolation.
The formula. for the coefficients is found by taking the formula for the Hermite coefficients and setting m; 0 for all i.
aj -((xi-1 :x4)/(1++ 8 (t.//(t+1 tz)² 0
b,
3(1,+1
de
Pros and Cons. The mono cubie is both smoother than the linear interpolation and is increasing unlike the Hermite cubic. Therefore it is good to use for interpolating time or obsorvations counts. However, is not as good for interpolating features, since zero derivative is enforced at the observation points which not realistic inductive bias
F.4 Natural Cubic
The natural cubic is a piece-wise cubic that is twice-contimnously differentiable everywhere. So the function, its first derivative and its second derivative are all continuous at and between observation points.
We start by constructing a system of linear equations for a vector
k t2)+ f441) k4t1(t, - 1+11 Fi-1 3q Fi Fi+1
t, - t, t, - t,+1
This for 0 n 1. We also have boundary conditions ko 0 and kn 0. This system of linear equations can be solved using the tridiagonal algorithm which is O(n), rather than O(n²) for Gaussian
30

---

Published in Transactions on Machine Learning Research (09/2023)
climination. After finding the vector k its elements can be used to find the coefficients.
di
(ki+1 ki)/6(li+1 - 1;)²
b;
k,/2
Ci (r++1 - (1:41-
xx) - 14) (ks+1 + 2k:) - t,)
d;
24
From the above formula we see that the second derivative of the interpolation at point i is h; Additionally, from the formula for & we see that the derivative point is what we expect from linear control signals (#e41 - [t,/ f.) with corrective term to make the first and second derivatives continuous.
Pros and Cons. This is the smoothest of the interpolants with continuous second derivative. This makes the ODE solve faster. However, it is not physical since it uses the whole series to find the coefficients part of the tridiagonal solve. This means natural cubic can only be used when we enforce physical solutions using the method of copying subjects, this cannot be used with recti control signals.
G Library Specific Caveats
Our implementation uses TensorFlow and Keras, for which a differentiable ODE Solver exists. This can be found at httpps//www. tensorflow. htttp:://ww.teessorfoo... wtt.://rgg./.prbbbbb////aappaaap Loopi_docs/a hhocs//ythooo//t//mmat///// /te//////od/ We found two main caveats using this that do not. necessarily exist other libraries such as TorchDifiFq (PyTorch) We provide the straightforward fixes here to aid
(Chen et al.) 2018 Diffrax (JAX) (Kidger 2022).
implementation.
G.1 Running an ODE Backwards Time
For models such as the Latent ODE, an ODE solved backwards time, this when we solve from t to where t 2 to. The TensorFlow solver does not handle this by design, it assumes S to (time is ordered). We able to overcome this by solving related ODE, with ODE function. If previously
dh
di = t(h,t)
we use a new ODE
dh j(h,t) = -hh, -t)
dt
And now solve from -t, to tto with the new ODE function
h(to) ODESolve(f, h(t,), -t1, -to)-
And since -1 5 lb, the Tensorflow ODE solver can solve this. This requires that
ODESolve(f, h(tz) t1, -to) (itt)
which we prove below.
Let 7 f= -t, and solve a new ODE for h(i), where h(tt) h(ti). We carry out solve over an infinitesimal step from & to & Af giving
h(6 + Aî) -h(T)) F(h(H),
=h(E)) - -(h(E), -h(11 (()()))
Which is equal by definition to h(t, + At), (where t is negative). As with the proof of minibatching Appendix D repeatedly applying this gives h((t)) h(()), completing the proof. So we can solve an ODE backwards by solving related ODE forwards with initial and final times multiplied by -1.
31

---

Published in Transactions on Machine Learning Research (09/2023)
G.2 Repeated Calls to ODE Solver
The TensorFlow ODE solver is able to take sorted list of times to solve up to, and it functions fully in this scenario. Therefore for models like Neural CDE and GRU-ODE we do not need to repeatedly call the solver, we can just ask it to solve for times of interest. For the ODE-RNN and by extension Latent ODE (as it uses an ODE-RNN encoder), we must repeatedly call the solver. This is because we discontimuously change the state at observations, and so we must stop the ODE solve, update the hidden state and call the solve again.
Repeatedly calling the solver prevents the model from backpropagating, preventing the models from training. In fact it doesn't it simply prevent the models from training, but the program stops running. The fix to this problem to explicitly pass parameters into the ODE function, and then also pass these as "constants" to the ODE solver, so that knows to calculate their gradients. We include simple example code below:
This code firstly defines the ODE function that explicitly has parameters written out and named. We demonstrate two hidden layer ReLU multilayer perceptron with tanh final non-linearity (as is the case with dynamics functions). Following that we show how it is used with TensorFlow solver.
import tensorflow as tf import consorflou_rrbabbl11y as tfp
class Repeatable0DEFunc (tf. layers.
def ceinit__ssell x_dim, b_dim): super () init__()
self w1 tf. Variable (tw.addtw..oomalaaaaaaaaaaaa___/ h_dim)), trainable=True)
self.bi tf oaa.apaeeeeaaaaaa.o.ooooaaaaaaaa
htw.aaaee.....aaaaa.oooooaaaaaaa (ha. trainable=True)
self.w2 tf vaccaballec....raaaooo...ooo.loa
h h_dim)), trainable=True)
self.b2 tf Wariable (erarpolloff.f.aaooooo..ooooooad
trainable=True)
self.wt w3 = tf rrrfaaileff..rnddmm.oommmallll (shpprrehhshahppeeehehhe__a x_dim)), trainable=True)
self.b3 = tf. Ww..aaraalett.....cceaa normal (shape= trainable=True)
def call(self, t, x, w1, b1, w2, b2, w3, b3):
x tf. nn rewu(tf relml(ff..aammmll (x, w1) b1)
x if..n.reeuu (etluttt..mtmm( matmul (x, w2) b2)
tf.math. SanahhtanhCeee (tf ..mmutllll (x, w3) b3)
return x
ode fn Repeatable0DEFuncGG (x_dim, h_dim)
x next tfp.mathh ode. DormandPrince() (
ode_fn.
t_now,
x_now,
solution_timos=/ [t_next],
constants
"w1": ade_fn.wi.
"b1": ode bd,_fn.bl,
"w2": ode ode_fn.we
"b2": b2,
"w3": ode fn,33,
"b3": ode b3,
},
). states [-1]
32


Reviews:
[{"review_id": "vxzmXvvQPz", "reviewer": "Reviewer_7QTB", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "See cell below", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "Contribution proposes a kind of reparametrization in the context of NeuralODE to have efficient batching during learning with irregular time series\nBy comparing different NeuralODE solvers on the task of MS disease progression prediction, authors show that Neural Controlled Differential Equations can outperform causal discrete models like temporal convolutional networks (TCN). The 3 other NeuralODE solvers do not perform as well.\nAlso to deal with continuously defined data given discrete inputs in a prediction context, smooth and causal interpolations are investigated. The recticubic method is promoted.\nThe work builds heavily on Roy et al. (2022)\n\n--\nSorry for the late review, I had not submitted it in the right place...", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "The paper requires a good practical experience with NeuralODE solvers to be fully grasped. Certain technical paragraph should expanded and rewritten with a bit more pedagogy:\n\n“And then at observations undergoes a discrete change.” Can you rephrase?\n“This evolves continuously between observations and at observations changes discretely based on its current value, the observation and context using gφc .” this needs rephrasing too. The full paragraph would benefit from editing to be honest.\nThe paper claims some computational gains but no details or comparison in running time is presented. I would suggest to be more explicit on computational costs.\n\nAlso I would suggest to add a pseudo-code of the best method to be able to follow what are input parameters, how data are batched, how gradient updates are done etc.\n\nI would appreciate some code in supplementary material to be able to experiment myself with the approach.\n\nTypos:\n\n“biological mechanisms so “ -> biological mechanisms “\nIt is critical that a model is casual -> It is critical that a model is causal\nFor example,this -> For example, this\nNote the Latent ODE … -> Note that the Latent ODE\nan further expert -> an expert", "broader_impact_concerns": "None", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "1GcPiehh8R", "reviewer": "Reviewer_QEEx", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "# Strengths\n\n1.  _The problem addressed by the work is of high importance_\n2.  The improvement over the results shown in the paper published in CHIL 2022 are notable (however see weakness #1 )\n3.  A detailed exposition of the techniques used in the experiments. I would conclude that the paper is well structured if not for weakness #4.\n\n\n# Weaknesses\n\n1.  *Inconclusive demonstration of performance of the proposed approach.* It is unjustifiably assumed that TCN is the state of the art model and comparison with it demonstrates the point. However, checking the previous paper that this work extends, one can see that the competing models are all weak (neither is specifically a time series model). Thus, beating TCN does not seem to lead to any conclusions about performance of the NeuralCDE model investigated in this paper. *The claims of outperforming the current state of the art are unjustified.* How about looking into classical forecasting time series models and their ensembles (see a couple of links below for the starting point)? Maybe look into Cyclic Boosting as well? I am not sure why this is not done.\n    -   [Example classical forecasting models](https://github.com/Nixtla/statsforecast/tree/main/experiments/m3)\n    -   [a package for forecasting models](https://github.com/Nixtla/statsforecast/tree/main)\n2.  *Highly specific to the proposed application.* Although the application to MS problem is extremely important, the paper, despite being technical, feels more appropriate in a specialized venue than at TMLR.\n3.  *Limited innovation.* Most of the ideas used in the paper, including the interpolation to create continuous signal, have already been published in the preceding paper. The main novelty is in application of NODE approach to the POM data for MS prediction, while previously it was a TCN. The elaborate interpolation schemes did not show a large effect in the ablation study, leaving the performance gain to the use of NODE.\n4.  Sufficiently convoluted writing. This is minor, but definitely makes the paper harder to read. Some of the terms appear before their definitions.\n    1.  Mathematical notation appears in the text without sufficient explanation of the meaning, e.g. \"c\" on page 4, \"recti\" on page 4 in item 2, what are the weights w, etc.?\n    2.  The experiments are not sufficiently described. How is the forecasting problem turned into classification? I guess that a prediction produced by the model at each time point is compared with the actual value or EDSS and each such comparison is considered as a classification event, but I am unable to confirm this in the paper besides asking the authors. Do the authors have EDSS available even where measurements are missing? Many unanswered questions\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper presents an approach to predicting the progression of multiple sclerosis (MS) modeling changes through continuous time via the NODE. The framework adapts existing continuous models to handle irregular and missing data through a preprocessing step and two new interpolation schemes. The preprocessing step standardizes sequences, imputes missing context features, pads sequences to the same length, and interpolates the new causal time-series using integers as pseudo-time. The interpolation schemes, including a monotonic cubic scheme for the time channel, are designed to interpolate missing values and act as continuous control signals for models. The authors demonstrate that this approach, which intrinsically handles irregular time-series and imputes missing values with physically plausible values (smooth), often outperforms a discrete model in benchmark tests on a public MS dataset. The paper underscores the potential of continuous models and machine learning methods in predicting MS progression and informing interventions and treatment options for patients.\n", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "1.  Justify why is this an ML paper rather than an application paper suitable to a specific venue (this is the most important concern)\n2.  Perform a thorough comparison with proper classical and modern time series forecasting models in the same settings with and without the interpolation schemes.\n3.  Better explain the differences with the CHIL paper in the text\n4.  Address weakness #4 above.\n", "broader_impact_concerns": "No concerns", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "fELntyuert", "reviewer": "Reviewer_v21S", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strength: \n\n- There are sufficient ablation studies which demonstrates the contribution of each technical details on the performance. \n\n- The literature review and comparison among different methods are clear and straight forward. \n\nWeakness:\n\n- The definition of label is not quite clear. Do the positive patients in 6-12 month also include patients who has already been positive in 0-6 month, or they are purely new on-set patients? The progression may be better defined in the framework of survival analysis instead of classification.\n\n- Lack of novelty. The paper only practices several existing methods on a new task. Also, the analysis on what leads to such discrepancy among different models is insufficient. ", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper benchmarks the multiple sclerosis prediction tasks with several continuous time models and extensively experiments different interpolation methods. The paper shows that the continuous time model outperforms discrete time models, and standardizing features outperforms interpolating features.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "The authors should provide more details on the definition of the label and inclusion criterion for each tasks. Also, more insights on the models will be helpful, since the paper does not have contribution on methodology.", "broader_impact_concerns": "The paper proposes a model for healthcare prediction. The performance of the model should be carefully reviewed before being deployed in practice.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "8NiVc9IkPM", "reviewer": "Reviewer_eSNb", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The presented work appear to be overall focused on application of CDEs, whereas the technical innovations and contributions are incremental and limited. For an application-focused paper, the empirical results are also relatively limited to support the stated contributions.\n\n1. The two main technical points, interpolate and mini-batching strategies, were motivated to have memory and computational gains. These however were not supported or ablated in experiments. Additionally, the computational requirements needed to pre-process the splines over the datasets should be discussed (which typically can be costly).\n\n2. The ablation results do not show significant difference between the interpolation schemes, which make it difficult to assess the contribution of this technical point.\n\n3. Parallel-time solvers that decouple the batch already exist. (torchode https://openreview.net/pdf?id=uiKVKTiUYB0). Relation between the presented pseudo-timing training and these related works should be discussed. \n\n4. The baseline models included, both in terms of discrete time models and continuous time models, are limited. \n\n5. Other than CDE, the rest of the NODE models do not show any improvements over the TCNs. This contradicts the motivation of the work about the advantage of continuous time models.\n\n6. Two schemes which enforce the models to only look at past data are introduced, but there is no data to show how they differ in terms of performance.\n\n7. The reduced feature set outperforming all features warrants further investigations and reasoning.\n\n8. Standardization works better appears to be one of the takeaways. The author posit it is because standardization takes X to a manageable range, however ddid not extend on it further. The effect of interpolation upon standardization also seems nominal, which the authors have pointed out. More insights on why could be helpful\n\n9. \tMore details on ψ can be helpful. How is ψ determined? Its nature? Requirements? Will any arbitrary increasing function work? To me, it does not feel like it will? So, which functions work? Does the choice of function affect modelling in any way?\n\n\n\n\nMinor comments.\n\t•\tUsing c to represent both context vector and spline coefficients is confusing at times.\n\t•\tThe word “causal” is used in a sense that it might not necessarily imply causation or the effect, as is common in literature surrounding do-calculus and causal statistics. ", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper presents the use of continuous time models, particularly neural C/ODEs, for longitudinal modeling of multiple sclerosis. This is an application-oriented studies with some methodological innovations in interpolation schemes and mini-batching strategies.Experiments were conducted to compare CDE and a variety of neural ODE variants to temporal convolutional networks (TCNs).", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Please see the detailed list above for a) necessary addition of experiments in ablation and baselines (1/4/6), 2) addition of clarification with respect to contribution of works (1/2/3/5), and 3) further investigations on key aspects of the experimental works (7/8/9)", "broader_impact_concerns": "No concern.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: The presented work appear to be overall focused on application of CDEs. For an application-focused paper, the empirical results are relatively limited to support the stated contributions.

After author's rebuttal, the reviewers still have major concerns on lacking comparison with more traditional methods for forecasting other than deep learning, and the limited novelty.  The comparing with a single method claiming it as a state of the art and refering to a previous comparison that did not use strong baselines either is not convincing. 

As noted that the authors have reduced their claim, and this paper's main claim is that benchmarking continuous models on this MS dataset and comparing to the previous best benchmarked model (TCN), any SOTA claims are removed.  Accordingly, I would like to recommend acceptance.  The authors are strongly encouraged to further revise the paper to address the comments related to empirical results.
Decision: Accept (poster)