Automatic and Structure-Aware Sparsification of Hybrid Neural ODEs with Application to Glucose Prediction
Authors: Bob Junyi Zou
Lu Tian
Venue: ICLR 2026 Poster (2026)
Keywords: Predictive Sparsity
Hybrid Neural ODE
Group LASSO
Glucose Prediction
Primary area: applications to physical sciences (physics, chemistry, biology, etc.)

Abstract: Hybrid neural ordinary differential equations (neural ODEs) integrate mechanistic models with neural ODEs, offering strong inductive bias and flexibility, and are particularly advantageous in data-scarce healthcare settings. However, excessive latent states and interactions from mechanistic models can lead to training inefficiency and over-fitting, limiting practical effectiveness of hybrid neural ODEs. In response, we propose a new hybrid pipeline for automatic state selection and structure optimization in mechanistic neural ODEs, combining domain-informed graph modifications with data-driven regularization to sparsify the model for improving predictive performance and stability while retaining mechanistic plausibility. Experiments on synthetic and real-world data show improved predictive performance and robustness with desired sparsity, establishing an effective solution for hybrid model reduction in healthcare applications.

Published as a conference paper at ICLR 2026
AUTOMATIC AND STRUCTURE-AWARE SPARSIFICATION OF HYBRID NEURAL ODES WITH APPLICATION TO GLUCOSE PREDICTION
Bob Junyi Zou
Lu Tian
Institute for Computational and Mathematical Engineering Stanford University Stanford, CA 94305
Department of Biomedical Data Science Stanford University Stanford, CA 94305
junyizouu
ABSTRACT
Hybrid neural ordinary differential equations (neural ODEs) integrate mechanistic models with neural ODEs, offering strong inductive bias and flexibilityy, and are particularly advantageous in data-scarce healthcare settings. However, excessive latent states and interactions from mechanistic models can lead to training ineffi- ciency and over-fitting, limiting practical effectiveness of hybrid neural ODBs. In response, we propose a new hybrid pipeline for automatic state selection and struc- optimization in mechanistic neural ODEs, combining domain-informed graph modifications with data-driven regularization to sparsify the model for improving predictive performance and stability while retaining mechanistic plausibility. Ex- periments on synthetic and real-world data show improved predictive performance and robustness with desired sparsity, establishing an effective solution for hybrid model reduction in healthcare applications.
INTRODUCTION
Hybrid modeling methods are receiving increased attention from the healthcare community because they combine inductive bias from mechanistic models with the flexibility of neural networks. These methods often prove especially valuable in small-data regimes--commonly found in healthcare, and medicine- by outperforming both fully black-box and purely white-box approaches in terms of predictive performance, robustness and interpretability (Ahmad et al. [2018] Du al. 2019: Mohan et al. 2019; Rackauckas al. 2020; Yazdani et al., 2020; Hussain et al. 202
Karniadakis et al.
2021
Qian et al. 2021 Sottlle et al., 2021 Zou et al. 2024).
In the field of dynamical system modeling. a significant category of hybrid approaches builds on neural ordinary differential equations (neural ODEs) (Haber & Ruthotto 2017] Chen et al.) 2018 Kidger 2021 which arose from the insight that deep residual neural networks can be formulated as continuous-time dynamical systems (Rico-Martinez et al.) 1992] Weinan (017). Neural ODEs well-suited for modeling dynamical systems as they offer continuous-tima representations with latent dynamics that integrate seamlessly into modern machine learning pipelines with automatic differentiation, making them both scalable and flexible. More recently, researchers have adapted crmmain-knowledgee inforgeeddnomeeddde relational inductive bias, often derived from
neural ODEs to incorporate mechanistic models causal graphs. This hybrid style- sometimes termed Graph Neural ODE (Poli et al. 2019), Neural Causal Model (Xia et al. 2021). Neural State Space Modeling (Hussain et al. mechanistic plausibility and interpretability while taking advantage of the flexibility of neural networks, thereby improving model performance and robustness in data-scarce settings.
2021 or Mechanistic Neural ODE (MNODE) (Zou et al.)
While hybrid neural ODEs show competitive performance in various healthcare and medical appli-
cations such as cardiovascular simulation (Grigorian et al.) 2024 Salvador et al., 2024). epidemic forecasting (Sottile et al. 2021 Huang et a 20224), disease progression and survival analysis (Dang
cet al.) 2023 Xiang et al., 2024), treatment effect estimation (Gwak et al. 2020 Zou et all 2024 and
pharmacology (Qian al. 2021 Hussain al. 2021), one persistent challenge in deploying them in

---

Published as a conference paper at ICLR 2026
practice is model reduction. Mechanistic models in physiology and medicine tend to become exces- sively large in attempts to capture wide-ranging and complex dynamics (c.g.. delays, heterogeneities, multt-compartment processes, etc.) and may contain dozens of latent states despite only handful of observable states. For instance, the state-of-the-e model for human carbohydrate-insuligg -llucose dynamies has more than 20 latent states, even though it only uses 2 input variables and less than observable state variables /Visentin et al. 2018). After hybridization, the added flexibility of neural components may render some latent states unnecessary or even detrimental when training data are scarce, as redundant states can significantly increase model variance, leading to over-fitting and undermining the benefits mechanistic models promise.
Traditional model reduction approaches in biochemistr- such as timescale separation (Michaelis Menten 1913, Johnson & Goody] ( and quasi-steady-state approximations (Schauer & Heinrich 1983 Bothe & Pierre 2010) ooten require deep domain expertise or extensive trial- and-error. On the other end of the spectrum, data-driven graph-based based model reduction offers a pathway solve this problem. In many healthcare application domains, mechanistic ODEs can be represented as reaction networks or directed graphs. where nodes denote system states and edges denote interactions (Hodgkin & Huxley 1952, Holz & Fahr 2001; [mmith et al., 2004 Canini & Perelson 2014; Man et al. 2014). In recent years, many solutions have emerged from the graph neural network (GNN) community for general graph pruning, using approaches such as topology- based node/edge selection (Spielman & Srivastaval 2008 Liu et al. 2023), learning-based sub-graph sampling (Wang et al., 2019 Zeng al. [2019] Zheng al.) 2020), or optimization-ben based graph sparsification (Li et al. 2020, Jiang et al. 2021 2202)). However, these reduction methods are typically data-driven and agnostic of any domain knowledge, and thus do not necessarily preserve key mechanistic structures or constraints. Furthermore, non- ron-gradieent-arr based reduction methods (e.g.. greedy search) can be prohibitively costly in computation for large, high-dimensional ODE systems. As result, a gap remains for computationally efficient solutions that reduce model complexity while preserving the mechanistic integrity and improving predictive performance for hybrid neural ODEs.
In this paper, we address this challenge by introducing hybrid, gradient- based algorithm for automatic state/edge selection and structure optimization in MNODEs. Our approach combines domain memain-kknwwler weeee--ooowwwownlleeeee ledge-informed graph modification with mix of L1 and L2 regularization that encour- ages graph sparsity to efficiently reduce model complexity. The graph modification step draws insights from classical reduction methods and graph theory to constrain the search space to mechanistically plausible sparse graphs that retain key topological structures. Meanwhile, the regularization step allows data-driven, gradient-based graph pruning during training, making the reduction process com- putationally efficient and adapted to observed data. By combining both mechanistic and data-driven elements, our reduction scheme integrates the best of both worlds and particularly well-suited for modeling complex dynamical systems in healthcare and medicine with limited data. Through extensive experiments on both synthetic and real-world data, we demonstrate that our algorithm outperforms other reduction strategies for MNODEs and also surpasses unreduced MNODEs and widely used black-box sequence models terms of predictive performance and robustness using less parameters. These findings highlight a promising path toward more efficient and effective hybrid modeling caattiioss-particllary in settings where high-quality data are scarce and model stability is erucial.
METHODOLOGY
2.1 Preliminary
Mechanistic controlled ODE system We define a mechanistic controlled ODE system M as a 4-tuple M = XS,
1. S = fsi,..., 8n is the set of state variables with cardinality n, and (t) : [0, +00) - R is a real-valued function of representing the value of state s; time t.
2. = { F' is the set of exogenous input variables with cardinality m, and xty(t) [0, +0o) - R is function of L representing the value of input xy at time l.
2

---

Published as a conference paper at ICLR 2026
3. F (fi...... fr} is the set of real-valued functions of S, X and that describe the system's temporal evolution:
(s;(t) fi (a(s)) pats) pats))(), (t),))
dt
where Spa(i) S S, ppt()) S X are subsets of state and input variables, respectively, on which the derivative of 8; with respect t depends, i.e., they are "parents" of s;.
4. So (0), ... sn (0)} is the set of initial conditions.
In addition, we can further split the state variable set into two disjoint subsets:
S observable states Sobs LJ latent states Stat-
Observable states are variables the system that can be directly measured through experiments or sensors. These are the quantities that can be collected and tracked over time. On the other hand, latent states are variables that are not directly accessible but still believed to play a role in system dynamics.
Directed graph representation of mechanistic ODE We define the directed graph representation of M = (S, X, F, So) as directed graph GM = (VM; EMs), whose node set and edge set are defined in the following way:
SMM {si...., I1, I1m..., (8). E EM + 8, Ept)) and (z), E Eιs + E
Specifically, (sy, 8i) E EMs means that the value of (t) influences the "direction" of (t))). Similarly, (I,, 84) E Eas means that the value of (t) influences the "direction" of (t).. Note that allow self loops-we can have (8;, 8;) E ΣM, if (1)/ dt depends on 8; (1) in the system. In the rest of the paper, we will the following definitions:
Relaxed directed acyclic graph (RDAG): We define a relaxed directed acyelic graph to be a directed graph with no directed cycles, except for self-loops. Note that the directed graph representations of mechanistic ODE systems are in general NOT RDAG.
2.2 Prediction task: time series forecasting with dynamic covariates
The main task we are interested in is to predict the future trajectory of observable state variables, given their observed history and both past and future exogenous input signals. More precisely, given:
(I) Past context: SP S (tx) J-mp) E [(p+1)x|Saal X (X(+)) E IR*x where L-p < ..... 1-1 < to = 0 are a set of discrete time stamps which observations of Sats and X are collected, and to is the beginning of the prediction window;
(2) Future inputs: X {X(+s) E Raxm, where 0 = t0 < < ... to are future prediction time stamps in the prediction window;
we want to predict ST. = (Soo (tw) E IS the future value of the observable states.
Data: The observed data consist of copies of {S obs: St Xt3 from multiple instances and the objective is use observed data to train an algorithm prospectively predicting observable state values in new instances based on history, {S X X²}.
2.3 Model architecture; mechanistic neural ODE (MNODE)
At a high level, MNODE follows the encoder-decoder sequence modeling paradigm, in which the encoder takes in historical context and output initial condition estimate of the latent states in the system and the decoder rolls out predictions based on the initial condition and future inputs:
Encoder( SP X²) = SL)()), (ecodee($(()).) x+) SF
where 5(()= (S... (0), SO) RD is initial condition estimate. In general, MNODE is compatible with any choice of encoder layer as long as the encoder can produce a reasonable estimate of the initial condition of the system. For the decoder layer, given the directed graph representation of the mechanistic ODE system Gm, future exogenous inputs X and an initial condition estimate S(0) E R", MNODE initializes node features in GM as Sto = $10), X = X and evolve state node features over time using a of ferw-forward neural networks [NN.). structured by GM
ds,(t) (Np (t), (t), (t), S(0) = Smm (0) = Encoder(Spmr , x.Pstt dt
(1)

---

We denote the resulting super-graph as GO = (Vo; F²). Each super-node in Va may collapse multiple observable state nodes into single "super-state" node, whose feature is defined as the concatenation of all observable node features within its MSCC. Let Sah C va be the set of super- nodes whose MSCCs contain least one observable state node, and X Va the set corresponding to X C V in G. For consistency, we similarly define S", S Ste pa(s)) and Xa pa(t)- This yields an Published as a conference paper at ICLR 2026
In practice, the solution of equation D can be approximated by a forward-Eule style discretization:
1h+1
AA (th+1 tn NN, pu(s)
X pa(s). tn),
(2)
where h = 0, 1,.... and we switched the notation from (at()) to Spnca) to emphasize the transition from continuous time-domain to discrete time grid. our implementation. we choose the encoder layer to be standard LSTM and the feed-forward neural networks [NN;) J to be standard MLPs.
2.4 Reduction Algorithm: Hybrid graph sparsification (HGS)
Step I: merging maximal strongly connected components Given directed graph representation of mechanistic ODE system G = (E, V) (since the dependency on the mechanistic ODE system is clear from context, we will omit the M subscript to simplify notations), we first collapse all maximal strongly connected components (MSCCs) in G into super-nodes to make it an RDAG (note that we allow self loops). This is implemented by first partitioning V into disjoint subsets of MSCCs C,:
V = =I Vi C,nCj = 0.
Next, we define a super-node set V by mapping each MSCC in V super-node in V
v = {cf CV, 1 isisk)
Then, to define edges between super-nodes, for each directed edge (u, E E, we add (ct, eg) to the super-edge set E", where C₁ and Cj are the two (not necessarily different) MSCCs contains u and v, respectively:
E = [(c), (u, E F, u u. U.C.,
RDAG representation of the mechanistic model, G
Rationale of step Transforming the original graph into an RDAG via collapsing the MSCCs reveals high-level causal structure of the system, simplifies the interpretation, and provides a topological ordering with acyclic structure that improves training stability, as feedback loops are known to cause blow-ups, exploding gradients and stiffness of the ODE system. With cycles, many complicated constraints on the system parameters are needed to simultaneously control stiffness, blow-ups and exploding gradients. Without cycles, the system Jacobian is upper triangular after proper rearrangement, and the corresponding eigenvalues are simply its' diagonal elements, allowing substantially fewer and simpler constraints on parameters ensure system stability. (see Appendix A8 for detailed discussions and examples). In addition, replacing each MSCC with a self-loop does not sacrifice much predictive power because neural networks are capable of approximating the effect of complex intra-component dynamics key motivation of hybrid modeling Raissi et al. (2019).
While we have chosen the default set-up of HGS step to collapse all MSCCs, users may, based on application needs and their own domain knowledge, choose not to collapse certain MSCCs, and Step 2 and of HGS will still be compatible in these cases. Causal interpretability can be preserved via temporal unfolding, where feedback loops are resolved into time-lagged dependencies between distinct temporal instances of the state variables.
Step Customization:
Step 2: augmenting graph with simpler shortcuts Next, we identify key mechanistic pathways and augment them with simpler shortcuts for potential model reduction. To this end, let Da.s be the set of nodes, whose removal disconnects and in Go:
Dr,s = {u E y u z, v + s, s no longer reachable from rr in Go after removing v):
and let G g,s be the sub-gran aub-graph induced by (I, 8} U Dr.e:
x.8 = (V2,. E Ve_ = (x, U Dr.sr E2, = ((u,v) E E u.u U E V
Define the partial transitive closure G of G to be:
Ea
{(1,2))
(2)8) Eg
G (Vr,, E2
E
Fa
((,,)), (I,
(z. s) E
E,S
4

---

Published as a conference paper at ICLR 2026
where Eg,, is the edge set of the transitive closure of Gh using the reachability relation of G Finally, we augment the original RDAG G" with the additional edges from partial transitive closures of pathway sub-graphs to form the augmented RDAG G
G (V*, F* F. E" U F..
and Guee will be the graph used for step 3. Appendix A2.1 shows an example of G vs G
Intuition and rationale of step 2 Intuitively, one may think of a physiological path as a student's high-school journey moving through grades 9 to 12. Normally, the student progresses step by step-9 + 10 + 11 + 12. A transitive closure adds all possible sskip-grade links, letting the student jump directly from grade
11 or 12, or from 10 to 12, as long as they always move higher grade (obey the reachability relations). A partial transitive closure is more cautious version: it allows some skipping but forbids overly aggressive jumps, like going straight from grade to 12. The idea is that, just as students progress different speeds, biological processes/pathway. in physiological systems also vary in how many intermediate states they pass through and can therefore often be better modeled with fewer latent states. For example, quasi-steady-stard approximations in chemical kinetics eliminate fast variables by assuming equilibrium. By adding shortcuts (transitive closure), the model gains flexibility capture these differences without discarding realistic reachability constraints.
Step 2 Customization Using a partial (rather than full) transitive closure a choice made by the authors in the context of glucose modeling because it prevents introducing direct input-output edges unsupported by the mechanistic model, and preserving some latent dynamics has been shown to be important (Dalla Man al., 2009). Similar to step rather than following this default set-up, users may choose to include full transitive closure or omit selected short-cut paths based on their needs.
Step 3: applying a mix of L1 and L2 regularization Given a processed RDAG G automatically remove redundant edges and nodes, a natural way is to associate weight with each edge and apply L1 penalty to shrink weights of redundant edges zero in the style of LASSO regularized regression. In the context of MNODE, straight- straight-forward formulation would be to modify Equation to:
ds; (t) dt
(N. @ (t)))) W Xpu(s) (t), O.),
where the ith neural network is parametrized by θ, W = (w( (u, E Ead is the set of edge-specific weights, and O stands for element-wise multiplications:
W O Sp (1) (u(s..))) s(t) s Spa))),, W X X (W(()))) xx1) x X (o))))
Given the mechanistic RDAG Ga defining the MNODE structure, NN parameter θ = (O....., On), and edge weights W. one may predict the state variable values (i.e. node features) ti...., ty with an initial condition estimate and future exogenous inputs Xa.F time. These predictions can be recursively calculated based on (2):
86 sp:th44 S..th + (ta+1 - t,. (W O
Sa,ta W @ xa pa(s))
pa(s)) txi O2), with sq. = =..
0.
We denote the resulting prediction of observable states Ša at ln by S (Sa.d, XaF. O W, Ga.s) to emphasize its dependence relevant model parameters, initial condition and exogenous input. We estimate encoder and decoder parameters simultancously by minimizing the weaa-squured-error loss function. To encourage graph sparsity while retaining identifiability, also place combination of La and L2 regularization on edge weights W and model weights θ to form the final loss function:
Σ h
S0.th obs
gaala (Sov (DaP; X Fw,,x,,F,WC.... Go,
+ +
Σ (u.
|wwool + Aa 10122. (3)
where D (S X²*P) are observed data available time 0, Su (:: is the encoder generating the initial condition of the system, and À is a penalty parameter. The I, penalty edge weight W is designed to encourage sparsity and the L2 penalty on decoder parameters is to boost identifiablility.
Equivalence to first-layer group LASSO: The above regularization is closely related to the idea of first-layer group LASSO mentioned (Cherepanova et al. 2023). It can be shown (see Appendix A2.2 that Equation iis equivalent to
S0.0h (5eato, X (T, (), Ga.)) 1, 1|- 2211 A
Σ
|
(v...)) 2/3
b
(u,e)E Es.e
5

---

Published as a conference paper at ICLR 2026
where λ3 = x 2- X2 Ô represent all non-firstllayer weights in the MLPs, and To.. "(a,e) ((a,,))(u.u)) is the (((a,)) scaled vector consisting of firttlayer-muliplicatinn weights associated with the edge (u, z). The E(,,))E (n.ag) 21//3 term is a variant standard group LASSO penalty encouraging group sparsity, and the vector I (a..e) = =0 + removing edge (u, v) E bas. Compared to the standard group LASSO penalty that raises |[] (w, ll to power of 1, a smaller exponent encourages stronger group sparsity with steeper gradient towards O. Operationally, the regularization parameters À are selected via K-fold cross-validation (CV).
Step Customization In step while the default set-up is penalize all edge weights W, if certain edges are deemed indispensable or more important by the user, they can be taken out of the regularization term or be given separately customized X.
2.5 Important note: implausibility of true support recovery should be pointed out that our method is not meant for true support recovery (i.e. recover the true underlying causal graph of the data generating process) because the expressivity of neural networks lead to equivalent MNODE models whose underlying graphs are different (Xia et al.) 2021), and the task of recovering the true graph structure is therefore theoretically implausible without making strong assumptions about ground truth. Our goal is efficiently induce model sparsity and facilitate the generation of data-driven hypotheses, which require further clinical validation, without making strong assumptions about the ground truth other than that it is more sparse than the original mechanistic prior.
3 RELATED WORK
Sparsity, LASSO, and Generalization LASSO-style penalties have long been used to induce sparsity. Applications range from residual networks (Lemhadri et al.
2021).
varying-coefficien
Wen et al. 2016), to feature selection in MLPs (Zhao et al. 2015 Sun al. 2016 Wang al. 2017 and local linear sparsity (Ross al. 2017). However, these methods often yield sparsity patterns that are hard to interpret. Our approach extends LASSO to MNODE while grounding the learned sparsity in mechanistic graph structures, providing interpretability and domain -amain-aligned insights. Prior work has shown that reduced neural networks can generalize as well as, even better than, their unreduced counterparts (Gale al. 2019 Bartoldson al. 2020 Hoefler al. 2021). For graph-based model reduction, You et al. (2020) found that mildly sparse relational graphs often outperform fully connected ones. Our results support this: MNODEs built from lightly pruned graphs perform better than those using dense graphs in low-data regimes.
models (Thompson et al. (023), and CNNs via group LASSO (Liu et al.
2015
Graph Sparsification via Optimization Our method is related to optimization-ban based graph sparsifi- cation approaches for GNNs. For example, Li et al. (2020) constrained the Lo norm of the adjacency matrix, Jiang al. (202l) used elastic net penalties, and Jiang al. (2023) applied exclusive group LASSO encourage neighborhood sparsity. Unlike these methods, we do not sparsify the adjacency matrix directly. Instead, we sparsify edge weights in message passing, reducing influence from less informative neighbors. Crucially, GNN sparsification methods are typically data-driven and ignore mechanistic structure. In contrast, our method is structure-awand aware:: it begins with domain- pruning step that restricts the search to physically plausible graphs. This makes our approach more suitable for hybrid models like the Graph Network Simulator (Sanchez-Gonzaler et al. 2020) that require mechanistic consistency. Finally, while Zou et al. (0024) used a greedy, stepwise reduction scheme, our method is more computationally efficient and yields better performanee-anllgouus to the gains of LASSO over stepwise selection in linear models (Hastie et al. 2020).
Landscape of Hybrid Modeling and Trade-off's There many variants of MNODE in the current landscape of hybrid modeling with different degrees of hybridization and strength of mechanistic

---

Published as a conference paper at ICLR 2026
prior. In Zou et al. 2022).. general form of mechanistic neural ODE is proposed as:
dS(t)
di = cm (S(t), ((), (()) + c:m (t)t) h000000000...0000001000000000010 ca + (1 - ca)1)
((t) = c18m + (1 - (S)), and dZ = fmnm (tt), (t))
dt
where fm. Pm and Am represent the functional form, model parameters and dependency structure (adjacency matrix) of the mechanistic prior, respectively, and fons. i = 1, 2 are neural networks. For C2 0, = 1, one recovers MNODE. The deep mechanistic simulator in Miller
al. 2020) uses C₁ 1, = 0. The neural closure learning model in Gupta & Lermusiaux
(2021) uses 6 0.5, C2 = C3 c4 0. The hybrid ODE model in Qian et al. (2021) uses
C1 = 0.5, C9 = 1, C3 = 0, while the standard black-box neural ODE uses C1 = C3 0, = 1. In general, the more mechanistic components, the more constraints on the function space of the resulting hybrid model. In data-limited regimes, such constraints help the model quickly learn reasonable representations but also undermine the model's ability capture complex real-world dynamics not represented by the mechanistic prior. On the other hand, the less mechanistic prior, the more flexible the hybrid model becomes. While the model gains more freedom to learn arbitrary patterns, it is also exposed to risks of overfitting and high variance. The challenge is to find the sweet spot: using just enough mechanistic guidance to keep the model grounded without too much constraining.
4 EXPERIMENTS
4.1 Experiments on synthetic data
General Setting We assume the data are generated by an unknown sparse dynamical system (ODE), but the mechanistic model is overly complex and contains redundant variables. We show that HGS produces models that are more sparse, predictive and robust than those obtained by existing methods.
Data generation We consider two sparsity regimes: true sparsity- feature have zero effect size, and quasi sparsity-redundant features have non-zero but small effect sizes. Our synthetic data are generated from the following two controlled ODE systems respectively:
ds1(t) (.5 1] (t),
True Sparsity:
dt
X
(s) 0.5(s; 1] + 4
21(1)
Quasi Sparsity:
dt
101-1"
==1
We generated time series samples of length == 60 based on the forward Fuler numerical integration scheme with time step Vt 0.05 over the domain t [0, 0.3] with zero initial conditions.
Training Sample Size To study the effect sample size on our method and validate its effectiveness limited data regime, for each sparsity regime, we generated 40 independent training sets of size 100 and 1000 respectively, as well as a held-out test set size 10,000.
Starting graph Assuming the true data generating process is unknown, we consider two settings in which the starting "mechanistic" graph contains redundant structures: (1) refined graph whose redundant part contains input nodes, latent node, and latent cycle; (2) comprehensive graph whose redundant part contains 6 input nodes, 3 latent nodes and 3 latent cycles (See Appendix A4. for illustrations). All redundant input variables are generated from independent N (0, 0.5).
Baseline models We selected the following baseline models: (1) Block-bor box sequence models
including: LSTM. Black-box neural ODE (BNODE), temporal convolution network (TCN) LLaa et al. 2016), Diagonal S4 (S4D) (Gu et al.) 2022) and vanilla transformer (Trans); (2) MNODE reduced by other reduction methods including: reduction (NR), NeuralSparse (NS) (Zheng et al.
7

---

Published as a conference paper at ICLR 2026
2020), exclusive group LASSO (EGL) (Jiang al. 2023), elastic net (EN) Jiang et al. 2021). random search (RD) and step-wise greedy search (GD). Implementation details are in Appendix A4
Evaluation and metries All evaluations are performed on the test set (see Appendix A5). We report RMSE with I-sigma standard error (SE) for predictive performance, Peak (worst-case) RMSE for robustness and Effective Number of Parameters (ENP, average number of parameters whose magnitude is > 10-3 under CV-selected hyperparameter setting) for model reduction. Additional metrics including MAPE, Peak MAPE and correlation are reported in A7.1
Results As shown in Figure 1 HGS outperforms black-box models small sample sizes, with the gap narrowing as sample size increases. At n = 1000, TCN surpasses HGS in RMSE, but HGS retains superior robustness. This reflects a known trade-off: with more data, the bias from regularization may outweigh its variance reduction benefits, though the latter still improves worst-case behavior. Compared to other reduction methods, HGS consistently achieves the best performance. However, when the input graph is already refined, the gains are modest- expected, since most reasonable methods perform well in this setting. On graphs with substantial redundancy, HGS's advantage becomes both large and statistically significant. In such cases, regularization alone struggles to recover signal, especially under quasi-sparsity. HGS's hybrid pruning, when paired with regularization, addresses this challenge effectively (See Appendix (7222 for an ablation study). It also yields the fewest effective nonzero parameters (ENPs), highlighting its strength in inducing sparsity,
4.2 Experiments on real-world data: blood glucose forecasting for T1D patients
Introduction and background For the real-world data experiment, we focus on modeling the dynamics in patients with Type diabetes (TID). TID patients have impaired insulin production and therefore rely on constant external insulin delivery to regulate their blood glucose level, making their glycemic regulation challenging dynamical system to model, especially during periods of physical activities. In this case, the latent states represent different physiological compartments within the human body and the observed state is blood glucose level. We choose our mechanistic model to the 2013 Version of the UVA-Padove model (Man al. 2014) (Appendix A6), which is apppopeo for modeling glycemic response in TID patients.
Data Our data are from the TID Exercise Initiative (TIDEXI) available to at
Riddell et al. 2023), which is
ott.ddo..oo/.....5411... 2tw.:0/0225332/0000222222 Jaeb Center for Health Research Foundation, 2022). After pre-processing (Appendix A3 the final data contain 342 time series from 105 patients. Each time series consists of 54 measurements, taken minutes apart, of blood glucose level and exogenous inputs including carbohydrate intake, insulin injection, heart rate and step count from 210 minutes before to 60 minutes after the onset of an exercise instance. We set the first 210 minutes as history window and the remaining 60 minutes as the prediction window.
Baseline models We also consider MNODE reduced by domain knowledge (DK) (Zou et al. (024).
Evaluation and metries We adopt repeated CV (see Appendix A5 to estimate various metrics of interest. Specifically, we split data randomly at the exercise instance level evaluate the average predictive performance on both intra- and inter-patient instances, reflecting the likely real-world deployment of the prediction algorithm to both existing and new patients. In addition, since intra- patient variability can be as large as inter-patient variability in TID modeling (Moscoso- (Moscoso-Vasquez et al.) 2016 Bell et al. 2021 Laguna et al., 2014), we anticipate similar results from patient-level cross- validation (see Appendix A4.2 for details and empirical evidence). In addition to all 6 metrics used in synthetic set-ups, we include model variance for measuring robustness and an clinical significance metric, Diagnostic Accuracy-accurag, of classifying patient as hyper (2 180 mg/dl), in-range (80-180 mg/dl) or hypo (S 80 mg/dl) using model predictions. All metrics, except Peak RMSE, Peak MAPE and ENP, are reported with 1-sigma SE.
8

---

Published as a conference paper at ICLR 2026
(a) Against black- box models (training size 100)
(b) Against black-box models (training size 1000)
PLAK SPMC
(c) Against other reduction methods (refined initial graph, training size = 100)
(d) Against other reduction methods (refined initial graph, training size- 1000)
524-
JAR CUnC
1960
2090
(e) Against other reduction methods (comprehensive initial graph, training size 100)
(1) Against other reduction methods (comprehensive initial graph, training size 1000)
Figure 1: Comparison across different evaluation settings. GD is omitted for comprehensive initial graph as training takes an unreasonable amount of time. Models are sorted from best to worst.
Results As shown in Figure 2 a), MNODE.HGS significantly outperforms traditional black-box models across all metrics, while using fewer parameters. This highlights the benefits of incorporating mechanistic structure and graph sparsification into model design. Figure b) compares HGS with alternative reduction methods applied to MNODE. HGS consistently yields better predictive performance, particularly in peak RMSE, suggesting greater robustness. To visualize the learned structures, we plot in Figure[(c) edge-weighted adjacency matrices of the reduced graphs, averaged over 10 runs. Unlike other methods, HGS not only promotes sparsity but also introduces new structural shortcuts that are otherwise inaccessible to regularization-basee approaches.
A Ablation study on model components To assess the contribution of each design step in Section 2 we conduct an ablation study using models trained with various subsets of the proposed pipeline. As shown in Figure d), removing any single step leads to marked drop in performance, underscoring the importance of all three components for the success of MNODE HGS.
Mechanistic interpretation of results HGS can yield interpretable and biologically plausible new insights. For example, HGS chooses to eliminate edges corresponding to glucagon feedback loops, which suggests that impaired glucagon response during hypoglycemia (Seaquist et al. 2013) may also persist during exercise-induced hypoglycemia-a novel hypothesis that could guide future investigations.

---

Published as a conference paper at ICLR 2026
1P01
Depeodi Aocasay
LOTE
114
B104
(a) Comparison against black-box models; the proposed MNODE. is in orange.
(b) Comparison against reduction methods; the proposed MNODE HGS is in orange.
NR
DK
GD
R.D
NS
(c) Hcatmap weighted adjacency matrices of reduced graphs produced by various reduction methods
1912
PAPS
(d) Ablation study model components. Red/yellow/green color indicates inclusion of step 1/2/3 of HGS.
Figure 2: Combined Results for Real-world Experiments. Models are sorted from best to worst.
BROADER IMPACT
We propose a three-step procedure to simplify the mechanistic graph underlying MNODEs, aiming to improve prediction, robustness and interpretability. This especially useful in biomedical domains, where models often involve complex biological processes and high-quality data are limited. By leveraging
graph refinement, structural pruning, and edge-weight sparsification, method produces compact and predictive models that align with mechanistic priors while reducing overfitting. This enables more transparent and data-efficient modeling, potentially accelerating discovery in systems biology, personalized medicine, and related fields. More broadly, our framework contributes to the effort to integrate domain knowledge into deep learning in a principled way. It provides insights into which components (e.g. cycles, delays) of mechanistic model are predictive, offering guidance for experimental focus and hypothesis generation.
10

---

Published as a conference paper at ICLR 2026
ACKNOWLEDGMENT
This publication based on research using data from Jaeb Center for Health Research Foundation that has been made available through Vivli, Inc. Vivli has not contributed to or approved, and is not in any way responsible for, the contents of this publication. We thank Alex Wang, Dessi Zeharieva, Emily Fox, Matthew Levine and Ramesh Johari for providing assistance with data acquisition and preliminary discussions on TID and model reduction.
REPRODUCIBILITY STATEMENT
Our data processing pipeline is described in Appendix A3 Model implementation and training details are provided in Appendix A4.2 Model evaluation details are provided in Appendix AS The mechanistic model used is described in Appendix A6 Our code is also submitted as supplementary materials.
LLM USAGE
LLM only used to aid and polish the writing the paper. We did not use LLM for any other purpose.
REFERENCES
Muhammad Aurangzeb Ahmad, Carly Eckert, and Ankur Teredesai. Interpretable machine learning in healthcare. In Proceedings of the 2018 ACM international conference on bioinformatics, computational biology, and health informatics, pp. 559-560, 2018.
Brian Bartoldson, Ari Morcos, Adrian Barbu, and Gordon Erlebacher. The generalization-staailiyy tradeoff in neural network pruning. Advances in Neural Information Processing Systems, 33: 20852-20864, 2020.
Emily Bell, Sabrina Binkowski, Elaine Sanderson, Barbara Keating, Grant Smith, Amelia Harray, and Flizabeth A Davis. Substantial intra-individual variability in post-prandial time to peak in controlled and free-living conditions in children with type diabetes. Nutrients, 13(1 (111))4154, 2021.
Dieter Bothe and Michel Pierre. Quasi-steady-staie approximation for reaction-diffusiin system with fast intermediate. Journal of Mathematical Analysis and Applications, 368(1):120-13)) 2010.
Lactitia Canini and Alan Perelson. Viral kinetic modeling: State of the art. Journal of pharmacoki- netics and pharnacodynamics, 11:33-443, 2014.
Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud. Neural ordinary differential equations. Advances neural information processing systems, 31, 2018.
Valeriia Cherepanova, Roman Levin, Gowthami Somepalli, Jonas Geiping, C Bayan Bruss, Andrew G Wilson, Tom Goldstein, and Micah Goldblum. A performance-dritn benchmark for feature selection in tabular deep learning. Advances in Neural Information Processing Systems, 36: 41956-41979, 2023.
Chiara Dalla Man, Mare D Breton, and Claudio Cobelli. Physical activity into the meal glu-
cose-insulin model of type diabetes: In silico studies, 2009.

---

Published as a conference paper at ICLR 2026
Ting Dang, Jing Han, Tong Xia, Erika Bondareva, Chloe Niegele-Brownn Jagmohan Chauhan, Andreas Grammenos, Dimitris Spathis, Pictro Cicuta, and Cecilia Mascolo. Conditional neural ode processes individual disease progression forecasting: a case study on covid-19, 19. Proceedings of the 29th ACM SIGKDD Conference on Knowledge Discovery and Data Mining, pp. 3914-3925, 2023.
Mengnan Du, Ninghao Liu, and Xia Hu. Techniques for interpretable machine learning. Communica- fions of the ACM, 63(1):68-77, 2019.
Trevor Gale, Erich Elsen, and Sara Hooker. The state of sparsity in deep neural networks. arXiv preprint arXiv:1902.09554. 2019.
Gevik Grigorian, Sandip V George, Sam Lishak, Rebecca J Shipley, and Simon Arridge. A hybrid neural ordinary differential equation model of the cardiovascular system. Journal of the Royal Society Interface, 21(212):20207700 2024.
Albert Gu, Karan Goel, Ankit Gupta, and Christopher Ré. On the parameterization and initialization of diagonal state space models. Advances in Neural Information Processing Systems, 35:35971 35983, 2022.
Abhinav Gupta and Pierre FJ Lermusiaux. Neural closure models for dynamical systems. Proceedings of the Royal Society A, 477(2252):20210044, 2021.
Daehoon Gwak, Gyubyeon Sim, Michael Poli, Stefano Massaroli, Jaegul Choo, and Edward Choi. Neural ordinary differential equations for intervention modeling. arXiv preprint arXiv:2010.08304, 2020.
Eldad Haber and Lars Ruthotto. Stable architectures for deep neural networks. Inverse problems, 34 (1):014004, 2017.
Trevor Hastie, Robert Tibshirani, and Ryan Tibshirani. Best subset, forward stepwise or lasso? analysis and recommendations based on extensive comparisons. Statistical Science, 35(4):579- 592, 2020.
Alan L Hodgkin and Andrew F Huxley. quantitative description of membrane current and its application to conduction and excitation in nerve. The Journal of physiology, 117(4))500, 1952.
Torsten Hoefler, Dan Alistarh, Tal Ben-Nun, Nikoli Dryden, and Alexandra Peste. Sparsity in deep learning: Pruning and growth for efficient inference and training in neural networks. Journal of Machine Learning Research, 22(241):1-124, 2021.
Martin Holz and Alfred Fahr. Compartment modeling. Advanced Drug Delivery Reviews, 48(2-3): 249-264, 2001.
Yuxi Huang, Huandong Wang, Guanghua Liu, Yong Li, and Tao Jiang. Neuralcode: Neural compart- mental ordinary differential equations model with automl for interpretable epidemic forecasting. ACM Transactions on Knowledge Discovery from Data, 2024.
Zeshan M Hussain, Rahul G Krishnan, and David Sontag. Neural pharmacodynamic state space modeling. In International Conference on Machine Learning, pp. 4510- PMLR, 2021.
Inc. Jaeb Center for Health Research Foundation. Type 1 diabetes exercise initiative: The effect of exercise glycemic control in type diabetes study. Vivli Center for Global Clinical Research Data, 2022. URL CRRL:htts:://// http://///...../0/.22///////////
12

---

Published as a conference paper at ICLR 2026
Bo Jiang, Beibei Wang, Jin Tang, and Bin Luo. Geens: Graph elastic convolutional networks for data representation. IEEE transactions on pattern analysis and machine intelligence, 44(9):4935-4947, 2021.
Bo Jiang, Beibei Wang. Si Chen, Jin Tang, and Bin Luo. Graph neural network meets sparse representation: Graph sparse neural networks via exclusive group lasso. IEEE Transactions on Pattem Analysis and Machine Intelligence, (5(10):12692-120))) 2023.
Kenneth Johnson and Roger Goody. The original michaelis constant: translation of the 1913 michaelis-menter paper. Biochemistry, 50(39):8264-8265 2011.
George Em Karniadakis, Ioannis G Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang. Physics-informed machine learning. Nature Reviews Physics, 3(6):422-440, 2021.
Patrick Kidger. On Neural Differential Equations. PhD thesis, University of Oxford, 2021.
Suyong Kim, Weiqi Ji, Sili Deng, Yingho Ma, and Christopher Rackauckas. Stiff neural ordinary differential equations. Chaos: An Interdisciplinary Journal of Nonlinear Science, 31(9), 2021.
Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. International Conference Learning Representations (ICLR), 2015.
Alejandro Laguna, Paolo Rossetti, F Javier Sapudia-Blasco, Josep Vehi, and Jorge Bondia. Identification of intra-patient variability in the postprandial response of patients with type 1 diabetes. Biomedical Signal Processing and Control, 12:39-46, 2014.
Colin Lea, René Vidal, Austin Reiter, and Gregory D. Hager. Temporal convolutional networks: A unified approach action segmentation. In Gang Hua and Hervé Jégou (eds.), Computer Vision - ECCV 2016 Workshops, pp. 47-54, Cham, 2016. Springer International Publishing. ISBN 978-3-319-49409...
Ismael Lemhadri, Feng Ruan, Louis Abraham, and Robert Tibshirani. Lassonet: A neural network with feature sparsity. Journal of Machine Learning Research, 22(127):1-29, 2021.
Jiayu Li, Tianyun Zhang, Hao Tian, Shengmin Jin, Makan Fardad, and Reza Zafarani. Sgen: A graph sparsifier based on graph convolutional networks. In Advances Knowledge Discovery and Data Mining: 24th Pacific-Asia Conference, PAKDD 2020, Singapore, May 11-14, 2020, Proceedings, Part 24, pp. 275-287. Springer, 2020.
Baoyuan Liu, Min Wang, Hassan Foroosh, Marshall Tappen, and Marianna Pensky. Sparse convolu- tional neural networks. In Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 806-814, 2015.
Zirui Liu, Kaixiong Zhou, Zhimeng Jiang, Li Li, Rui Chen, Soo-Hyun Choi, and Xia Hu. Dspar: An embarrassingly simple strategy for efficient gnn training and inference via degree-based sparsification. Transactions Machine Learning Research, 2023.
Chiara Dalla Man, Francesco Micheletto, Dayu Lv, Mare Breton, Boris Kovatchev, and Claudio Cobelli. The UVA/PADOA Type diabetes simulator: New features. Journal of diabetes science and technology, 8(1):26-34, 2014.
Leonor Michaelis and Maud L Menten. Die kinetik der invertinwirkung biochem 49: 333-369 Find this article online, 1913.
Andrew C Miller, Nicholas J Foti, and Emily Fox. Learning insulin-glucose dynamics in the wild. In Machine learning for healthcare conference, pp. 172-197. PMLR, 2020.
13

---

Published as a conference paper at ICLR 2026
Senthilkumar Mohan, Chandrasegar Thirumalai, and Gautam Srivastava. Effective heart disease prediction using hybrid machine learning techniques. IEEE access, 7:81542-81554, 2019.
Marcela Moscoso- Moscoso-Vasquez, Patricio Colmegna, and Ricardo S Sanchez-Pena. Intra-patient dynamic variations in type diabetes: A review. In 2016 IEEE Conference on Control Applications (CCA), pp. 416 IEEE, 2016.
Michael Poli, Stefano Massaroli, Junyoung Park, Atsushi Yamashita, Hajime Asama, and Jinkyoo Park. Graph neural ordinary differential equations. arXiv preprint arXiv: 1910:353, 2019.
Zhaozhi Qian, William Zame, Lucas Fleuren, Paul Elbers, and Mihaela van der Schaar. Integrating expert ODEs into neural ODEs: Pharmacology and disease progression. Advances in Neural Information Processing Systems, 14:11364-1130 1383, 2021.
Christopher Rackauckas, Yingbo Ma, Julius Martensen, Collin Warner, Kirill Zubov, Rohit Supekar, Dominic Skinner, Ali Ramadhan, and Alan Edelman. Universal differential equations for scientific machine learning. arXiv preprint arXiv:2001.0438. 2020.
Maziar Raissi, Paris Perdikaris, and George E Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational physics, 378:686-707. 2019.
Ramiro Rico-Martinez, K Krischer, IG Kevrekidis, MC Kube, and JL Hudson. Discrete-vs.
continuous-tima nonlinear signal processing of Cu electrodissolution data. Chemical Engineering Communications, 118(1):25-48, 1992.
Michael C Riddell, Zoey Li, Robin L Gal, Peter Calhoun, Peter G Jacobs, Mark A Clements, Corby K Martin, Francis Doyle III. Susana Patton, Jessica Castle, et Examining the acute glycemic effects of different types of structured exercise sessions Type diabetes in a real-world setting: The Type diabetes and exercise initiative (TIDEXI). Diabetes care, 46(4):704-713, 2023.
Andrew Ross, Isaac Lage, and Finale Velez. The neural lasso: Local linear sparsity for interpretable explanations. In Workshop on Transparent and Interpretable Machine Learning in Safety Critical Environments, 31st Conference on Neural Information Processing Systems, volume 4, 2017.
Matteo Salvador, Marina Strocchi, Francesco Regazzoni, Christoph M Augustin, Luca Dede', Steven A Niederer, and Alfio Quarteroni. Whole-heart electromechanical simulations using latent neural ordinary differential equations. NPJ Digital Medicine, 7(1):90, 2024.
Alvaro Sanchez-Gonzalez, Jonathan Godwin, Tobias Pfaff, Rex Ying, Jure Leskovec, and Peter Battaglia. Learning simulate complex physics with graph networks. In International conference on machine learning, pp. 8459-8468. PMLR, 2020.
M Schauer and R Heinrich. Quasi-steady-state approximation in the mathematical modeling of biochemical reaction networks. Mathematical biosciences, 65(2):155-170, 1983.
Elizabeth Seaquist, John Anderson, Belinda Childs, Philip Cryer, Samuel Dagogo-Jack, Lisa Fish, Simon R Heller, Henry Rodriguez, James Rosenzweig, and Robert Vigersky. Hypoglycemia and diabetes: a report of a workgroup of the american diabetes association and the endocrine society. The Journal of Clinical Endocrinology & Melabolism, 1845))1885-1889 2013.
Bram W Smith, J Geoffrey Chase, Roger Nokes, Geoffrey M Shaw, and Graeme Wake. Minimal haemodynamic system model including ventricular interaction and valve dynamics. Medical engineering & physics, 16(2):131-139, 2004.
14

---

Published as a conference paper at ICLR 2026
Peter D Sottile, David Albers, Peter E DeWitt, Seth Russell, JN Stroh, David Kao, Bonnie Adrian, Matthew Levine, Ryan Mooney, Lenny Larchick, al. Real-time electronic health record mortality prediction during the COVID-19 pandemic: A prospective cohort study. Journal of the American Medical Informatics Association, 28( (8(11):334--236 2021.
Daniel A Spielman and Nikhil Srivastava. Graph sparsification by effective resistances. In Proceedings of the fortieth annual ACM symposium on Theory of computing, pp. 563-568, 2008.
Kai Sun, Shao-Hsuan Huang, David Shan-Hill Wong, and Shi-Shang Jang. Design and application of variable selection method for multilayer perceptron neural network with lasso. IEEE transactions neural networks and learning systems, 1366)1138-1336, 2016.
Ryan Thompson, Amir Dezfouli, and Robert Kohn. The contextual lasso: Sparse linear models via deep neural networks. Advances in Neural Information Processing Systems, 36:19940-19961, 2023.
Roberto Visentin, Enrique Campos-Náñez, Michele Schiavon, Dayu Lv, Martina Vettoretti, Mare Breton, Boris Kovatchev, Chiara Dalla Man, and Claudio Cobelli. The uva/padova type diabetes simulator goes from single meal to single day. Journal of diabetes science and technology, 12(2): 273-281, 2018.
Jian Wang, Chen Xu, Xifeng Yang, and Jacek M Zurada. A novel pruning algorithm for smoothing feedforward neural networks based group lasso method. IEEE transactions neural networks and learning systems, 29(5)22012-200)) 2024, 2017.
Lu Wang, Wenchao Yu, Wei Wang, Wei Cheng, Wei Zhang, Hongyuan Zha, Xiaofeng He, and Haifeng Chen. Learning robust representations with graph denoising policy network. 2019 IEEE International Conference on Data Mining (ICDM), pp. 1378-1383. IEEE, 2019.
Ee Weinan. A proposal on machine learning via dynamical systems. Communications in Mathematics and Statistics, 1(5):1-11, 2017.
Wei Wen, Chunpeng Wu, Yandan Wang, Yiran Chen, and Hai Li. Learning structured sparsity in deep neural networks. Advances in neural information processing systems, 29, 2016.
Joseph M Worsham and Jugal K Kalita. A guide to neural ordinary differential equations: Machine learning for data-driven digital engineering. Digital Engineering, pp. 100060, 2025.
Kevin Xia, Kai-Zhan Lee, Yoshua Bengio, and Elias Bareinboim. The neusal-neural connection: Expressiveness, learnability, and inference. Advances in Neural Information Processing Systems, 34:10823-10836, 2021.
Jinlin Xiang, Bozhao Qi, Mare Cerou, Wei Zhao, and Qi Tang. Dn-ode: Data-driven neural-ode modeling for breast cancer tumor dynamics and progression-free survivals. Computers Biology and Medicine, 180: 108776, 2024.
Alireza Yazdani, Lu Lu, Maziar Raissi, and George Em Karniadakis. Systems biology informed deep learning for inferring parameters and hidden dynamics. PLOS computational biology, 16(11): e1007575, 2020.
Jiaxuan You, Jure Leskovec, Kaiming He, and Saining Xie. Graph structure of neural networks, In International Conference Machine Learning, pp. 10881-10891. PMLR, 2020.
Hanqing Zeng, Hongkuan Zhou, Ajitesh Srivastava, Raigopal Kannan, and Viktor Prasanna. Graph- saint: Graph sampling based inductive learning method. arXiv preprint
2019.
15

---

Published as a conference paper at ICLR 2026
Lei Zhao, Qinghua Hu, and Wenwu Wang. Heterogeneous feature selection with multi-modal deep neural networks and sparse group lasso. JEEE Transactions on Multimedia, 17(11):1936-1948 2015.
Cheng Zheng, Bo Zong, Wei Cheng, Dongjin Song, Jingchao Ni, Wenchao Yu, Haifeng Chen, and Wei Wang. Robust graph representation learning via neural sparsification. In International Conference on Machine Learning, pp. 1158-11468.. PMLR, 2020.
Bob Junyi Zou, Matthew E Levine, Dessi P Zaharieva, Ramesh Johari, and Emily B Fox. Hybrid square neural ode causal modeling. arXiv preprint arXiv:2402.17255 2024.
A
APPENDIX
INDEX OF APPENDIX CONTENTS
A.2 2 Additional Details about HGS
A.3 Real-world Data Pre-pres proceessssing
A.4 Experimental Details
A.5 Evaluation
A.6 Mechanistic Model for Glycemic Regulation
A.7 7 Additional Results
A.8 Stability Analysis for Cyclic Systems
A.9 Tabulated Results
A2
ADDITIONAL DETAILS ABOUT HGS
A2.1
ILLUSTRATION OF HGS
Figure shows illustrative example of step 1 and 2 of HGS for simple graph.
A2.2 PROOF FOR GROUP LASSO EQUIVALENCE
To make the connection, note that in the context of MIP model computation, applying scaling parameter to input feature is equivalent to re-scaling all the first-layer weights linked to the corresponding feature by a factor w. Therefore, suppose we let ((,,)) be the vector consisting of mouplayee-ulllplliaaiinn weights associated with the edge (u, v) Eu.s, and é be the vec- tor consisting of all non-firt-layer-mulliiiiee
model parameters all MIPss, we can rewrite ga.th (,, O W. G).) and the regularization term in as obs
Sa (1:3 (1 (m,)) (1000) E E 6) 1, Ga.e)
and Σ (Aalura + A2 (1.)) (1 + x211@112
Eo.e
respectively, where 1 is the set of unit edge weights corresponding to the unweighted specification.
Thus, letting T... W(u 0(n,v) and I [I (u.,))) Eacc), the loss function can be reparametrized as:
Σ h
S Soata (DaP. ;), x.F. :T, θ), 1, 1, 1000
16

---

Published as a conference paper at ICLR 2026
X
Partial Transitive Closure
Sut
X
Sat
X
Swt
Sut
Su
Sux
Collapse MSCC
Suse
Share
Sur
S000
Weight
Swt
Figure 3: Step and 2 of the Hybrid Graph Sparsification Algorithm
(1.,))11
Az
101 1 +
Σ
(elW((uoo) + A2
(4)
w ((,,0)
(u.
Note that (w,e) (u, v) Ea.e only appears in the last term of (4), One may minimize this loss with respect to W, while fixing other parameters and obtain the minimizer
(2
(2)g (1.,)
W (m..))
X
Substituting ((,,o) back into 4 gives the desired loss function with respect to only I, and
B
obs Sa.tn (ga.fo, X. (T, (r,)).1 1| de 221l0|1 A
Σ
[[ (v.)) 12/3
2/3
cases, fo
(u,e)EEu.
we replace (w) by w (m,oo= then the same derivation can show that the equivalent loss function becomes
|tc obs
5th (B)
G
patients,h h
+ 121 +2 A1de
Σ
| (ewll
(,
with a standard group LASSO penalty on I.
A3
APPENDIX: REAL-WORLD DATA PRE-PROCESSING
A3.1
SELECTION
We select patients on open-loop pumps with age under 40 and body mass index (BMI) less than 30. From exercise instances recorded by these patients, we focus on the time window between 210 minutes before exercise onset and 60 minutes after exercise onset, and select those that satisfy the following conditions:
17

---

Published as a conference paper at ICLR 2026
1. The exercise lasted for at least 30 minutes.
2. There are no missing blood glucose readings (recorded every 5 minutes) or heart rate reading (recorded every 10 seconds) in the time window of interest.
3. There is at least one carbohydrate intake in the time window of interest.
With the above selection criteria, we end up with 324 exercise instances from 105 patients.
A3.2 FEATURES, UNITS AND INTERPOLATION
For each selected exercise instance, we use the following features derived from the TIDEXI raw data:
Feature
File FA FACM
Data Field
Unit
UVASim Unit
CGM Reading (Blood Glucose Concentration)
FATEST
mg/dL
mg/dL
Basal Insulin Rate
FATEST
U/oou
U/min
Bolus Insulin
FACM FAMLPM VS
FATEST
U
U/min
Dietary Total Carbohydrate
FATEST
mg/min
Verily Heart Rate
VSCAT
bmp
NA
Verily Step Count
FA
FATEST
count
NA
Table 1: Summary of features and their units in TI DEXI and UVA/Padova
Feature
Notation
Ascending time
i-th entry of CGM reading and its time-stamp i-th entry of basal insulin rate and its time-stamp
(gi, t
Yes
(a;, 12)
Yes
i-th entry of bolus insulin and its time-stamp
(bi, t
Yes
i-th entry of carb and its time-stamp
(ma, mm
No in original Data, need sorting
i-th entry of Verily HR and its time-stamp i-th entry of Verily Step Count and its stamst
(hiy
Yes
Yes
Table 2: Notations for various features
A3.3 INTERPOLATING BASAL FLOW RATE
Since the basal flow rate is already RATE, interpolating it is relatively straightforward We approxi- mate basal flow rate with step function where the magnitude and time of jumps are determined by a, and 1f respectively:
0g/60
to <t<t
15i<N
fa(t)
ON
/60
/6
IN <
0
otherwise
where the unit for fa is U/min, the same as the units used for insulin delivery rate UVA/Padova; and is the total number of basal insulin rate data entries.
A3.4 INTERPOLATING BOLUS INSULIN RATE
The recorded data are dosage and the corresponding time. To convert dosage to rate, we refer to real-world experiences where most open loop pumps inject a dose of bolus insulin at a rate of 1.5U/min. We ignore basal insulin when computing the effective bolus rates because basal insulin rates contribute negligibly to the maximum rate
18

---

Published as a conference paper at ICLR 2026
First, to account for cases of overlapping bolus doses (a new dose is applied while the previous dose has not been completely delivered), we pre-process the bolus insulin data with the following trick: two bolus doses (by, (2), (bj+1, 4)+1) overlap (i.e. (5+1 < 19 + b//1.mmin,, it is as if we only introduced one dose of bolus insulin of b + by+1 time t.. We can recursively apply this trick to combine all overlapping doses (see Algorithm for details), after which we end up with new list bolus insulin data (b), 12) in which we can assume there are no overlapping doses.
Algorithm 1 Pre-process Bolus:
Input: Bolus Insulin Data B {(b). 1.)) where by is U and og min.
i=1
Initialize B an empty list
while <- length (3) do
if 0 < +1 + /// then
B(i) = (b, + (+11,9)
delete (b,+1, 11 from in place
else
i
Append B(i) B
end
end while Return
Next, we define the continuous-time time bolus delivery function ft( as
1.5
(b5/1.5)min otherwise
st<
fo(t)
0
where the unit for ft) U/min as well.
A3.5 INTERPOLATING CARBOHYDRATI INTAKE RATE
Suppose the patient consumed carbohydrates at a given time in the TIDEXI dataset. We assume constant meal consumption rate of 45000 milligrams per minute:
(1) ||4000,, M {(m,, (m)) DatalU,, S15 Cn + +ms/45)
A3.6 INTERPOLATING HEART RATE AND STEP COUNT
We interpolate heart rate and step count using rolling window average with a window size of minutes:
(a) = [H(t))
Σ
h
H = [(h, th) E Datalt $ t S t + 5min}
(h.A)EEH())
fo(1) [V(O)
Σ
v, V [1v ()) Datalt & t 15min)
(M)))))())
A3.7 CHOICE OF TIME GRID AND DISCRETIZATION
Since our algorithm uses forward-Fuler style numerical integration scheme to solve the continuous neural ODE, we need to discretize the continuous, interpolated input features before we can use them. We first choose a suitable discrete time grid for each exercise instance. For each selected exereise instance, we focus the time window from 210 minutes prior to the start of exercise, to 60 minutes after the start of exercise. We choose discrete time grid the CGM measurement time stamps
19

---

Published as a conference paper at ICLR 2026
within the time window of interest, and since CGM measurements consistently taken in 5 minute increments, this splits the time window into minute intervals (with each interval containing one measurement) and we obtain 54 discrete time steps:
t, = = iAt, i= 0=.....,53 53,
5min.
Given this time grid, the processed data used by all the models in the main paper are described in T Note that we do not interpolate CGM readings and use them they are recorded, which is made possible by choosing the time grid to match the CGM recording time stamps t; lg t In
Feature
Definition T. ty
Time Stamp
CGM Reading
= 9:
Average Insulin Injection Rate
TTR:
at S fa(t) + fo(t) de
1
Average Carb Intake Rate
M, At pl fm(() dt
Average Heart Rate
H: (t,tt)
Average Step Count
V ft(t))
Table 3: Notations for various features
conclusion, each exercise instance, after above pre-processing, is turned into 5-dimensional time series with 54 time steps of the form:
(G. IIR,, M, H,, Vi
In Figure we provide graphical illustration of the data pre-processing pipeline.
-000
in mgdi
bacal
interpolated flow Uthour
KaWr solus in Intergilated Bolus Rote uniin
Intrepolated
Meall cart in
Inteepolatee meet hake righmin Seoothhed meal rate migitnan
Pan Step Count Processed Step Count
-240
220
130
150
220
-00
30
Time frim acnerise moeetiii
Figure 4: An illustration of raw and pre- process data of one of the exercise instances.
A4
APPENDIX: EXPERIMENTAL DETAILS
A4.1
SYNTHETIC MECHANISTIC GRAPH
In 5 provide illustration of the mechanistic graph used in the synthetic experiments.
20

---

Published as a conference paper at ICLR 2026
X
X
X
X
Sat
X2
X,
Sut
Sut
Sue
X
Sut
S
S,
Sur
Sux
Star
X
SyMt
SaLt
Sut
Sut
Refined "Mechanistic" Graph
True Sparsity Graph
Comprehensive "Mechanistic" Graph
Figure 5: An illustration of the mechanistic true graphs used in the synthetic experiments
A4.2 EXPERIMENTAL SET-UPS
Synthetic data set-up
Synthetic data are generated with the following script:
def gen_syn_data
rng=np oogemp.ranoom.dffalllt..mmm (seed=seed)
n=60
(1, n, reshape (n,
data=[]
for k in range (train_size):
[=((++))0000tpp.xpggg1///////// ( (
+rng. nerganormall (0, 0.5, (n,1)) for i in range
for in range (3): // change to 6 for comprehensive graph append (rng. normanormllllm (0, 0.5, 5, (n,1))) (x, axis=1)
dt=5e-2
s1=[0]
s2= [0]
s3= [0]
v= [0]
for in range(n)
if exp==1:
else:
v. v.append (|100110|1||11111000000000000010
-0..00++4(4.3330055-11111111)1
use apapped (appppddf(((ndcdxxxxx1100001a110
4ee44x
21

---

Published as a conference paper at ICLR 2026
-0.5* (v[-1]-1)))
for comprehensive graph
sample==n..concaaee ([np. reshape (n,1)), axis=-1) apta.a case==np.array (data)
noi wolbse-rng.staaddddddomm.. (aperrmaalsiiiiieeeeeeeeee (size=cases shape, dtype= float64') cass=nn.concatennan [cases, noise i, ;,]]], axis=-1)
return cases
Real-world data set-up For each model mentioned in the experiment section, here offer detailed description of the corresponding computational method and the hyper-parameters used. Throughout this section we are given an exercise time series of 54 time steps (corresponding to 54 minute intervals that made up the time window starting from 210 minutes prior to exereise termination, and ending at 30 minutes after exercise termination) and features (corresponding to CGM reading, insulin, carbohydrate, heart rate, and step count, in that order). We denote the first feature (CGM reading) as $1, and the other 4 features time steps and subscript indicate feature indices. For example, = is the 1st input feature of (carbohydrate) discrete time step t = 12-
We use superscript to indicate discrete
Our goal is to predict the CGM trace during the first 60 minutes following exercise onset corresponding to the output s1:12 E R and therefore we set the number of prediction steps to be 12 for all models. We further split the given time series into historical context D (S ubs XPit--ait-i E R starting glucose s/o E R. inputs during exercise f IR 12x4 (twelve inputs that recorded 1 time step ahead of the expected outputs). We use s 12 to indicate the CGM trace predicted by models, for model states and z for black-box latent states, h, for the final hidden state and cell state of the LSTM initial condition learner. For ease of computation and without loss of generality, we set the At term in forward-Eun Eulerr style discretization to be 1 for all relevant models, and thus we omit it in the equations.
Data Splits Since synthetic data are generated independently, each training set is split simply by the default index into 4 subsets for cross voos--vlli,tion, and the first split is used for re- re-training and obtaining the final model after the optimal hyper- parameter setting has been obtained. For real- world data, since there might be correlation between adjacent samples in the original pre-processed data, for each repetition, we generated random permutations with:
rng=np.. random. default_rng (seed=2024)
perms=np. zeros (repeats, cases casess shape [0] dtype=" int32
for in range (repeats):
perms [i] rrgg... wwwpppeetatttt.. ion (cases. [0])
print (cases. shape)
and apply the permutation before the standard index-based 3-fold train-validation split. Again, the first split used for obtaining the final model with optimal hyper-parameee setting.
Why Instance Level Random Permutation In real-world deployment, the model is expected to be applied to both existing patients (with new instances) and new patients, making cross-validation based on random instance-level splitting more appropriate, as naturally includes predictions for both seen and unseen patients. In contrast, splitting strictly at the patient level would ignore the model's ss ability to generalize within the same patient across different instances nontrivial and practically important challenge in diabetes modeling, Furthermore, it is well recognized that intra-patient variability-the variation in glucose dynamics across different instances of the same
22

---

Published as a conference paper at ICLR 2026
patient-can can be high inter-patient variability in the diabetes research community. In other words, predicting future glucose trajectories for an existing patient can be as challenging predicting glucose level for a new, unseen patient. To quantify this phenomenon, we computed the root- mean-squared difference (RMSD) of mean and standard deviation of glucose levels both within and across patients in our data. Let e; and 8; denote the mean and standard deviation of glucose in instance and define iintrr {(i,j)
< J. instances and j are from different patients) and J, instances and j are from the same patient). The intra- aitra-patient RMSD of mean glucose values, i.e., / (a (e, e,, was 54.24, compared an inter- patient RMSD of 50, i.e.. (r.j) (e; e,)². Similiarly, the intra- and inter-patient RMSD of glucose standard deviations (v (iJ) (s, - s,)² and (a.)) (s, - 8y)²) were 25.75 and 22.67, respectively. Those summaries empirically confirming that intra- intra-patient variability is comparable to inter-patient variability and result from patient-level cross- validaaiin would likely be similar to that from instance-level cross-validation.
Nintr {(i.j)
Initialization and Optimizer For all experiments, we use the Adam optimizer (Kingma & Ba 2015) to perform stochastic gradient descent. We initialize all model weights with the Py PyTorch default setting with seed 2024 + r, where increases by 1 starting from for each experiment repetition (40 repeats for synthetic and 10 repeats for real-world experiments) unless state otherwise the detailed model computation algorithm below.
Training Epochs In synthetic experiments, we train for 600 epochs and pick the epoch with best validation loss. For woaldworld experiment, We train for 200 epochs and pick the epoch with best validation loss.
Hyper-parameter Search We use grid search to tune hyper-parameters including learning rate and dropout rate. When choosing the grid, we restrict the search space to areas where the models have less than 20,000 parameters and we also try to limit the number of grid points to be less than 50. This is to make sure the computational cost of the experiments is capped at reasonable level for small data sets and individual users. The grid used for each model in each experiment will be provided below together with model descriptions.
We also enforced constraint to simplify the optimization: if a latent node has only one incoming edge and one out-going edge, i.e., #1 - U - u2, then those edges share common weight W Wa,ug"
Weight-sharing
Computing Resources and Time Usage All experiments are performed machines with Ubuntu 22.04 operating system, Xeon Gold 6148 CPU and single Nvidia 2080 Ti RTX 11GB GPU. Wall- clock run-time for both real-world and synthetic experiments range from hour per repetition to 12 hours per repetition, depending on the size of the starting graph and the type of algorithm used, with MNODE GL and S4D being the fastest ones and MNODE GD and transformers being the slowest ones.
A4.3 MNODE WITHOUT REDUCTION (MNODE_NR)
We take the directed graph representation of Muva- denoted as Guva construct MNODE_NR.
Hyper-parameters for real-world Experiments The LSTM has 2 layers and |Viaal hidden dimen- sion (i.e. this corresponds to the number of nodes in the directed graph) and for the same reason we set the 1st features of the initial condition state vector to be observed initial glucose level sf". All ML.Ps have 2 hidden layers and 16 hidden units with dropout O and activation ReLu. We place Ta
23

---

Published as a conference paper at ICLR 2026
Algorithm 2 MNODE_NR
Input: historical context DP starting glucose fo exogenous inputs X the directed graph representation (as defined in Section 2) of the UVA PAA-aadova model Guva- At h, LSTM(DE)
Sto h/0]
for i = Q do
9,+1 8; At-
Ui+1 81+1
end for
for i 11:99 do
St. 54 AA MLPs(S, x Guva)
end for Output: bu
regularization on all MLP parameters with penalty learning rate Ir. We tune these hyper-pameten with grid search on the following grid:
A2 and train the model with
A = {10 i 3.... x Ir = {10- {100,,10' 10 3).
A4.4 MNODE REDUCED BY DOMAIN KNOWLEDGE (MNODE_DK)
We take the reduced UVA-Padova model from Appendix C of Zou et al. (2024), and use its directed graph representation denoted Gruva construct MNODE DKR. Other than the graph used by MNODE, model computation exactly the same the MNODE NR.
Hyper-parameters for real-world Experiments The LSTM has 2 layers and Viwval hidden dimen- sion (i.e. this corresponds to the number of nodes in the reduced graph) and for the same reason we set the 1st features of the initial condition state vector to be observed initial glucose level s/r. All MLPs have 2 hidden layers and 16 hidden units with dropout O and activation ReLu. We place L2 regularization on all MLP parameters with penalty hyper-panee prprrpparreeter X2, and train the model with learning rate Ir. We tune these hyper-pameters with grid search on the following grid:
A = (10-² ,,....83 = (10-2 10-
A4.5 MNODE WITH HYBRID GRAPH SPARSIFICATION (MNODE HGS)
We omit the implementation detail of MNODE_HGS here as it is already discussed in great detail in section 2
Hyper-parameters for real-world Experiments The LSTM has 2 layers and | hidden dimen- sion (i.e. this corresponds the number of nodes in G% and for the same reason set the 1st features of the initial condition state vector to be observed initial glucose level 8f0. All MLPs have hidden layers and 16 hidden units with dropout and activation ReLu. We place L, regularization on all the edge weights with penalty hyper-paramend parameter \1, Lg regularization on all MLP parameters with penalty syper-parameter X2, and train the model with learning rate (r. We tune these hyper- pameees with grid search on the following grid:
A1 {10 10 10 7) A2 (10- 6. 8) x Ir {10 10-3).
Synthetic Experiments We also set 8/" s and uses the modified graph obtained from the given mechanistic graph. We keep everything else the same and the hyper-parameter search grid is:
A (10-5 10-7) A2 {10- = 6....,8) x br (10-² (0-2,10-3) 10-33.
24

---

Published as a conference paper at ICLR 2026
A4.6 MNODE REDUCED BY EXCLUSIVE GROUP LASSO (MNODE_EGL)
MNODE EGL uses the same model architecture as MNODE HGS except (1) MNODE EGL uses Guva (Viva; Euva) instead of Gava (2) the regularization term is defined as:
Σ
Σ
DE Vau
(u,
We tune hyper-pameters with grid search on the following grid:
= {10- 3,.....,, 8} x lr (10-2, (
A4.7 MNODE REDUCED BY ELASTIC NET (MNODEJEN)
MNODE_EN uses the same model architecture as MNODE_HGS except (1) MNODE_EGL. uses Guva (Vuva; Euva) instead of Gava, (2) the regularization term is defined as:
Σ
)1(u)(,0)) A2m((.))
(u,0)E Fur
Hyper-parameterr for real-world Experiments on the following grid:
We tune these hyper-pameters with grid search
A = {10-3, 10 10-00 x A2 = {10-* 1000 i= (), x f10-2,11000000000000010111
2,1100000 20.0...
A4.8 MNODE REDUCED BY NEURAL SPARSE (MNODE
The neural sparse algorithm tries to learn distribution (parameterized by neural networks) form which good sub-graphs are samples. Its computation is given below:
Algorithm MNODE.NS
Input:
historical context D ,starting glucose s exogenous inputs fou the directed graph
representation (as defined in Section 2 of the UVA-Padova model Guva, Al = 1, K size of the sub-graphs to sampled
Index edges in Euva by some order
Initialize model parameter Initialize model parameter e ~ (O0, 1),
NO0 1),
E R Evol E BK Fonal
и softmax(())
exp ((log(T) log( - log(s))) 10)
W = Σ, m:;
Round w to 2 decimal places
Construct Ey {e, E E w, > 0}. and construct G/vva (Vivay E Run MNODE NR with DP. s eto Givv and return output
for real-world Experiments The paramttasssssers and model configuration for the MNODE_NE part of the model is identical to the implementation of MNODE NR above. In addition to >2 and Ir, tune K, the maximal number of edges to be included in the sampled sub- -rbbgraph. These hyper-parameters are tuned with grid search on the following grid:
A {10-) | ...... 8) lr 2 10-1003 $3 K {12, 16, 20,
Hyper-parameters for Synthetic Experiments parameter search grid is:
We keep everything else the same and the hyper-
A (10 Ir ( K ( .... 10).
25

---

Published as a conference paper at ICLR 2026
A4.9 MNODE REDUCED BY GREEDY SEARCH (MNODE.GD)
The greedy search is implemented the standard step-wise backward way; at each iteration, the algo- rithm considers all existing edges and choose the one whose removal leads to the most improvement in validation loss, and stop when no edge's s removal improves validation loss. Specifically: Since
Algorithm 4 MNODE.GD
Input: historical context DP starting glucose s exogenous inputs (fo--1, the directed graph representation of the UVA-Padova model Guva. At = Index edges in Luva by some order
Stop=FALSE
MinLoss= 107
while Stop = FALSE do for e, Euva do
Construct G² = (Vwa; Euva {ei))
Run MNODE NR with DP, s X G", record validation loss as 4 end for if min( min({1)))1) (1))(E))) < MinLoss then
MinLoss= min(f4))er (Emm)
Evva = Euv (Eargmin(1.))
else Stop=TRUE
end if
end while
Construct G = (Vuva. Euva)
Run MNODE.NR with DP, stp Xlao-1, ,," and return its output
G
the greedy algorithm can be extremely slow, we do not tune hyper-parameters and instead fix the Le penalty hyper-paramettr to be 10-6 and learning rate be 10-3. All other settings the MNODE.NR.
A4.10 MNODE REDUCED BY RANDOM SEARCH (MNODE.RD)
The random search randomly picks 5 sub-graphs that has contains - percent edges and select the best one. Specifically:
Algorithm 5 MNODE RD
Input: historical context DF starting glucose 8f0, exogenous inputs X fo the directed graph representation of the UVA-Padova model Guva, number random sub- -rbpgraphs R = 5, sub graph edge ratio P = (0.1, 0.2, 0.4). At =
for E P do Uniformly sample R sub-gran aub-graphs of Giva that have the number of nodes and p) Fuval number of edges, denote them as Gp.1...., Gp.R for 1: R do
Run MNODE NR with DP, t0, fog--1 Gpory record validation loss as 4p.r
end for end for
Return the output of MNODE NR with Gargmin(tp.r?
for real-world Experiments As in greedy search, random search is also slow on real-world data. Therefore, do not tune hyper-parameters and instead fix the L2 penalty hyper- parameter be 10-3 and learning rate be 0.02. All other settings the same as MNODE NR.
26

---

Published as a conference paper at ICLR 2026
Hyper-parameters for Synthetic Experiments In the synthetic case, the graph is smaller and can afford tune A2 over {10-6). The learning rate is fixed 10-3 All other settings the same as real-data.
A4.11 BNODE
For BNODE and the subsequent black-bon models, we point the reader to implementations referenced in the associated citations in the main paper. Here we describe our implementation.
Algorithm 6 Black-box Neural ODE Model
Input: historical context DP, ,ttrringg glucose s/r, exogenous inputs Xt 1, At
h, c LS))
Sta 1s[0]
S =
for i = -1:quu
St: Stt-1 + At MLP
end for Output: g/
Hlyperpparameter for real-world Experiments The LSTM has 2 layers and d hidden dimension. Note that here the hidden dimension of LSTM also determines the state dimension of the neural ODE, which is tunable hyperparameter. All MLPs have hidden layers and 16 hidden units with dropout and activation ReL.U, trained with learning rate of Ir we tune these hyper-parameters with grid search on the following grid:
{6. 12, 18) x (0, (0,00.10.2]
[10 10-3 100-1.
Hyper-parameterss for Synthetic Experiments We keep all the settings the same the real world experiments
A4.12 TCN
Algorithm 7 Temporal Convolutional Network Model
Input: historical context DF, ,starting glucose 5 exogenous inputs X Xte--1,
fl 0 E IR9
1 X' concatenate( X, X, dim features)
t
segin concatenate( DP. X', dim time)
sefout TCN(seqa))
81 Linear( (SeQut)
Output: fii
pyrer--ara for real-world Experiments The TCN model is taken directly from the code repository posted on https////ggghubb..oooooooaaaaaaa
py] with input size set to 5, number of channels set to list of n copies of m, kemel size set to and dropout set to trained with learning rate Ir. We tune these hyper-parameters with grid search the following grid:
(2,3) (16,32} l - (2,3,4) x (0.0.1,0.) x lr = 210-2,10 10-3 10 10-1).
A4.13 LSTM
27

---

Published as a conference paper at ICLR 2026
Algorithm 8 Long Short Term Memory Model
Input:
historical context Dº ,starting glucose s exogenous inputs X hi e Encoder LSTM(D")
Set initial hidden state and cell state of Decoder LSTM to h, respectively seqoot; bg. 4q = Decoder LSTM(X
$1 Linearr(seean
Output: sf
Hyper-parameete for real-world Experiments Both Encoder and Decoder LSTM have n layers and hidden states with dropout set to a, trained with learning rate ir. We tune these with grid search on the following grid:
n = (2.3) x (6, 12, 18) x a =
x lr (10 2010 2,10
Hyper-parameters for Synthetic Experiments ments.
The settings are the same as real-world experi-
A4.14 TRANSFORMER
Algorithm 9 Transformer Model
Input: historical context DP, starting glucose s exogenous inputs X true output $1 41 (needed during training)
X
R2 t x/
X" concatenate X, X, dim features)
encoder in concatenate DP, X', dim time) (concatenating all inputs to form a masked context)
if Model in Training Mode then
decoder_in in co 81 - (expected output shifted to the right)
decode Transformer( wwn.oreeeeeeeeeeerreeeeeee.e in, decoder in,
mask)
end if
Model in Evaluation Mode then
decoder_in concatenate( (s 0 R²-1)
for = 1 9 do decode Transformer(encoder in, decoder in)
decoder_in,:: decode_out;
end for decode_out Transformer((nccoderiin decoder.in)
end
shig $
Linear(decode..uu
Output: sfa:
Sapeee-prraeetrs for real-world Experiments We use the transformer model provided by the pytorch nn class, and its syper-parameters are set as follows: d_model set to d, number of attention heads set to number of encoder layers set to 2, number of decoder layers set to 2, the dim fiidfoedforwrrd is set and dropout is set to a, trained learning rate of 10-3, We tune the paremetas with the following grid:
= (8,16) x n [4,8) x {16, 32) x {0,0.1)
28

---

Published as a conference paper at ICLR 2026
A4.15 S4D
Algorithm 10 S4 Diagonal Model
Input: historical context DP, starting glucose s exogenous inputs X
X 0 R29
t = $1
X' concatenate X, X, dim = features)
sefinn concatenate DP X', dim time)
setfnn Linear(Segia
seqin = Transpose Seqin: 1,2)
seqout S4D( seq.1)
sefout Transpose Segout, 1, 2)
$1 dhi Linear( Linear(seqout)
Output: s!
TIDEXI Experiments We take the S4D model directly from the following github repos-
itory
com http:///ww................ttt.t/
py and its hyper-parameters are set as: d_model set to d, d_state set dropout set to trained at learning rate of Ir. We tune the hyper-parameters with the following grid:
(4, 6,8 {m = (32 a = (0, 0.1,0.)) x = {10- 2,,00 3,0
A5 APPENDIX: EVALUATION
In this section we describe how evaluate the performance of various learning algorithms. Suppose training data D {=1, [ Z are i.i.d. samples from some distribution P, and learning algorithm M - F maps of arbitrary size a function fD. We also have an evaluation loss function (which may not the training loss function) L F Z + R U {0} that maps a learned function and a test sample to non-negative value. We are primarily interested in evaluating the expected prediction error (also known as generalization error) of the algorithms on unseen training and test data:
EPE(M, L, P) = E2..P.D< 2)].
In our setting, can be written as (S Xa.P ...FF ça.F, and the learned function maps the input (Sar Xa.P. X².F) to the output ( S xa.PP .x.FF) to approximate s
A5.1 SYNTIEETIC EXPERIMENTS
In the synthetic experiment, first generate test set D ..,,... _M} with sufficiently large sample size M 10000, where $,, 2.M are i.i.d samples from distribution P. We then generate K = 40 copies of training data, Dr, k = 1, ..... K, of size N = 100/ 1001000 each, perform the experiment K rounds with the K training sets and K different random seeds, and obtain the corresponding prediction function Jn.. We choose I. to be the standard mean squared error loss function and RMSE estimator for the k-th experiment round is computed as:
RMSE
1
Σ
lIfo. (Sap
X...P Xa.F) S.a.F
M
The reported RMSE is the average RMSE over the K rounds:
K
RMSE
RMSE
K
29

---

Published as a conference paper at ICLR 2026
and the la Monte carlo standard error is computed as
s.e.
RMSE RMSE,)
K(K - 1) R=1
A5.2 REAL WORLD EXPERIMENT
In this section we describe how we evaluate the performance of a learning algorithm using observed data only.
Note that when the loss is the mean squared error, we have the following variance-bais decomposition of EPE:
EPE Ez,pllifn(Sth, X xx.P) S 112) variance + bias"
noise
variance = Ez,n [[ xar, xaa)) En EnL obs
Xap X
bias
Ez {|EEoloossan X xx.)) xxp, - -
Soa X X
obs
noise = Ez
E[S
Sa. X xx.P) s ubs
In a synthetic experiment, one may sample from P as much as possible and estimate all the above quantities as described section A In reality, however, unknown and we only have limited data. Therefore, we have to resort to techniques such as cross validation or bootstrap. While the standard K-fold cross validation can give unbiased estimator for EPE, it cannot estimate variance or bias because each sample only enters the test set once. Due to the implausibility of observing repeated samples with the same input in this case, make the ideal assumption that there is noise in predicting state variables, i.e.,
E[S obs I X Xx.) S
so that we can estimate bias and variance separately. Under this assumption, noise zero and EPE = bias" + variance. Specifically, in this paper we use a modified K-fold cross-validation to construct unbiased estimators for variance and bias. First, split D into disjoint subsets of equal size D = _ km Dr, denote (- = D Dr. Then, each ( consists of K 1 disjoint subsets D,
D( U D., (EI(-R)
where 11-10 {i fillsis K,i k). Our estimator for variance is:
K
1
Σ
Σ
Σ
Fe, (S X xx.P)) - SDC- (Sa. Xa.P XaPP Xa.))
variance
D
K-1
A=
(S
SEI
ID (SaP Xap Xa.P
K-1 Σ
So. (SaP Xa.P Xa.P
VET
bias 2
fDi- (SaP P XaP, ...P)
D
k=
(S
MSE variance + bias
Σ
Σ
Jo. (So. X.P XaPP xa.) obs
ga.P
N(K
k=
(S.X)ED& ie
30

---

Published as a conference paper at ICLR 2026
MEAL MODEL
GARTRO AESTEBAA TRACT
Pasno GRUESSS
Rate Appeorance
Eggrstion
LIVER
GLACORR
AND TISSUE
Procectin
INSUIN
BUSULES
Degradation
Raae AppaRInnCA
OLUCARCN
Degsasstio
Dec aiton
Plaane Okoupon
GLUCAGON DELIVERY
Apcenrince
Insulin Pump Model
Sensor Model
Controller
Figure 6: UVA/Padova Simulator S2013, taken from Figure of Man et al. (2014)
RMSE
MSE
We perform the experiment R 10 rounds using the above procedure, each round with different permutation of training data and a new random seed, to obtain R estimates of variance and RMSE, and report their average. Note that this method effectively estimates the variance and bias of prediction algorithm trained with N/K rather than samples. As a consequence, it tends to underestimate the prediction performance. Standard errors are computed in the same way by deviding standard deviation by the square root of R.
A6 UVA- SIMULATOR S2013
Here provide the exact full UVA-Padova S2013 model equations. Variables that given meaningful interpretations are model parameters.
A6.1 SUMMARY DIAGRAM
At a high level, UVA-Padova can be summarized by the diagram in Figure 6 which is taken from Figure in Man et al. (2014) It divides the complex physiological system into 10 subsystems, which are linked by key causal states such as Rate of Appearance, Endogenous Glucose Production and Utilization. Next, we will introduce each subsystem one by one and also explain the physiological meanings behind state variables.
31

---

Published as a conference paper at ICLR 2026
A6.2 GLUCOSE SUBSYSTEM
G EGP + Ra Uu - E' - k1 G + k2Gr
(5)
G. -- G1 K2G.
(6)
Gp/VG
(7)
Gp: Plasma Glucose G. Tissue Glucose, EGP: Endogenous Glucose Production Rate, Ra Rate of Glucose Appearance, Un: Insulin-independent Utilization Rate, Usa Insulin-dependent Utilization Rate, E Excretion Rate, V& Volume Parameter, Plasma Glucose Concentration
A6.3 INSULIN SUBSYSTRM
I = - (mz + ma) + m14 + Rai
(8)
I, -mm1 m3) m24p
(9)
I
(10)
T Plasma Insulin, T Liver Insulin, Rai Rate of Insulin Appearance, Vi Volume Parameter, Plasma Insulin Concentration
A6.4 GLUCOSE RATE OF APPEARANCE
Qst Qatol + Qsto2
(11)
Qotol = -k SpriQetol + 5
(12)
Qat22 kempt (Qeto) Qoto2 Qot1l
(13)
Quat Kaba Qgue + Kemmpe (Qrto) Q
(14)
Ra = fkabs Qgue(BBW) (BW)
(15)
kempt (Qsto) kmin (kmaz - kmin (aQ.to abD) tanh (BQ,to BcBD + 2)/2
(16)
Qstor: First Stomach Compartment, Qsto2: Second Stomach Compartment, Qout: Gut Compartment, Carbohydrate Ingestion Rate
A6.5 ENDOGENOUS GLUCOSE PRODUCTION
EGP = kpp Kp2Gp - k + E
(17)
X k2 (X- I,]
(18)
i, kk; - -T)
(19)
H -Ki X H
+ kH max(H H Hb)
(20)
L:; Remote Insulin Action EGP, ""; Glucagon Action EGP, I, Remote Insulin Concentration.
H Plasma Glucagon Concentration, H₁: Basal Glucagon Concentration Parameter
32

---

Published as a conference paper at ICLR 2026
A6.6 GLUCOSE UTILIZATION
Ua Fens
(21)
(Vmo + Vm + risk)]Gr
Uid
Kmo + Gt.
(22)
Km Gt.
X + P2u(1 - Ib)
(23)
0
G. G
risk
10(log(G) log(Gr)) 10(log(Gaa) log GT))2r²
Gth G Gb
(24)
G < Gth.
Fens: Glucose Independent Utilization Constant, X: Insulin Action on Glucose Utilization, I, Basal
Insulin Concentration Constant, risk Hypoglycemia Risk Factor, G2 Basal Glucose Concentration Parameter, Gt Hypoglycemia Glucose Concentration Threshold.
A6.7
RENAL EXCRETION
É kei max(G, kez, 0)
(25)
A6.8
SUBCUTANEOUS INSULIN KINETICS
Rai kal + Kaz
(26)
[sc1 = - -Kks + Kal /scl + IIR
(27)
isc2 kaJscl ha2tsc2
(28)
Icc1 First Subcutaneous Insulin Compartment, Isce: Second Subeutaneous Insulin Compartment,
IIR Exogenous Insulin Delivery Rate
A6.9
SURCUTANEOUS GLUCOSE KINETICS
C.=-IGG+TT.G
(29)
Ge Subcutancous Glucose Concentration
A6.10
GLUCAGON SECRETION AND KINETICS
it - SRH Raz
(30)
SRH = SR + SR
(31)
SR p [SR max (oz(Gth GG) SR HtHH000 0)]
G 2 Go
SR (o(C,n -G) SR 0)
(32)
G G
1+1
SR = 17 max(-
(33)
SR
First Glucagon Secretion Compartment, SRd, Second Glucagon Secretion Compartment, Basal Glucagon Secretion Parameter, Raw: Rate of Glucagon Appearance
SR;;
A6.11
SUBCUTANEOUS GLUCAGON KINETICS
Hoct (km1 H42) Hins
(34)
Hoc2 = khillset kna
(35)
Ray = kha H
(36)
33

---

Published as a conference paper at ICLR 2026
He First Subcutancous Glucagon Compartment, Hazz Second Subcutaneous Glucagon Compart-
ment, Hin Subcutaneous Glucagon Infusion Rate.
A7
EXTRA RESULTS
A7.1 ADDITIONAL METRICS FOR SYNTHETIC EXPERIMENTS
we observed similar trends and patterns in model performance for metrics including MAPE, Peak MAPE and Pearson Correlation as those observed in main text for RMSE.
Comments As shown Figure 7 to Figure 12
Figure 7: Comparing against black-box models, training size 100
FIAE MAPE
Figure 8: Comparing against black-box models, training size = 1000
MA
3.
Figure 9: Comparing against other reduction methods, refined graph, training size 100
A7.2 ABLATION STUDY ON HGS UNDER LIMITED DATA, TRUE SPARSITY REGIME WITH COMPREHENSIVE STARTING GRAPH
Comments To better understand why regularization based methods perform poorly on comprehensive graph, also performed ablation study on HGS model components and the results are shown in Figure 13 We that graph modification (step step 2) alone regularization (step 3) alone
34

---

Published as a conference paper at ICLR 2026
Figure 10: Comparing against other reduction methods, refined graph, training size 1000
MAV
PIALMAN
Figure 11: Comparing against other reduction methods, comprehensive graph, training size 100
Figure 12: Comparing against other reduction methods, comprehensive graph, training size 1000
is not effective (and even worse than no reduction) at improving predictive performance or robustness. It requires both to achieve the desired outcome.
Figure 13: Ablation study on model components, training size = 100, true sparsity, comprehensive graph
A8 INSTABILITY OF DYNAMICAL SYSTEMS WITH CYCLES: A TOY EXAMPLE
Here we discuss how cycles/loops in dynamical systems can lead to instability. There are three sources of numerical stability: blowing-up, exploding gradient and stiffness.
35

---

Published as a conference paper at ICLR 2026
A8.1 BLOWING-UP
Consider simple state system with circular dependence:
tsp (1) = as1 (t) + bs2 dt ds2 (I) (t) + d82 (t)
dt
The system will blow up when the Jacobian of the system,
=
has eigenvalues with positive real parts. That means that if one solves the system using the forward Euler method with a step size h, then the system with the update rule
$1 + kh) 52 + kh)
$1(t)
(hJ + T)*
s2
would blow up as k increases. Let us then consider the two eigenvalues of J:
1 2 (a + d + V(a d)2 + 4be)
(J)
be max(0, ad), then A+ (J) positive number making the system unbounded, even if and are negative.
On the other hand, suppose the system is acyclic except for self-loops (cxactly what our step is doing), say ds (t) dt depends on s2(t) but ds2(1)/dl (t) does not depend on ()) anymore, then J becomes
a b 0
d
with a and d as its eigenvalues. Then, the system is allowed to freely model how s2(t) affects dsj (t)/dt without concerns of explosion (no constraint on b), as long as a and are negative.
A8.2 EXPLODING GRADIENT
Now suppose we discretize the system and solve it with RNNs, then the gradient will also behave like (J + 1)* which will cause the RNNs to have exploding gradients when & is large.
A8.3 STIFFNESS
It is also well-known in both the physics and the neural ODE community that stiffness of the system can also cause numerical instability (Kim et al. 2021 Worsham & Kalita 2025).
Consider the dissipative version of the 2-state system with a circular dependence:
ds (t)
(t) + bb2((t)
dt ds2(1) dt
cs (t) -s²(t)
The eigenvalues of the system are
(J) - vbc.
< be 1, then A and the system blowing-ulleppooding gradient are under control. However, the system can still suffer from stiffness, which often defined as the ratio between the magnitude of the fastest to slowest stable rates:
A_
1 + vbc
A+
Vbc
36

---

Published as a conference paper at ICLR 2026
Table 4: Black-box comparison (mean + standard error 40 trials) (Sample Size N 100)
Model
RMSE
Peak RMSE
ENP
True Sparsity
MNODE
0.1039 + 0.0013
0.4124
757.55
BNODE
0.1237 + 0.0021
1.2451
3686.78
S4D
0.1130 + 0.0012
0.8247
452.78
TCN
0.1214 0.0049
1.8787
13671.23
Transformer
0.1699 + 0.0039
1.2015
10892.65
LSTM
0.1191 0.0021
0.9048
8510.78
Quasi Sparsity
MNODE
0.1056 + 0.0014
0.4971
760.73
BNODE
0.1249 + 0.0021
0.7979
3441.73
S4D
0.1137 = 0.0012
0.6873
456.38
TCN
0.1274 + 0.0056
2.4566
13154.93
Transformer
0.1698 + 0.0039
1.1801
10929.11 17
LSTM
0.1194 + 0.0021
0.8659
8338.50
Table 5: Black-box comparison (mean + standard error over 40 trials) (Sample Size N = 1000)
Model
Peak RMSE RMSE
ENP
True Sparsity
MNODE
0.0981 + 0.0010
0.2369
730.25
BNODE
0.1009 + 0.0011
0.5123
2219.60
S4D
0.0986 0.0010
0,5624
602.63
TCN
0.0875 + 0.0010
0.7938
14360.29
Transformer
0.2437 0.0066
1.2028
10696.88
LSTM
0.1003 + 0.0011
0.7060
8982.25
Quasi Sparsity
MNODE
0.0987 + 0.0010
0.2377
727.25
BNODE
0.1014 + 0.0011
0.5887
2330.20
S4D
0.0991 + 0.0010
0.4444
647.00
TCN
0.0896 + 0.0011
0.6814
14173.83
Transformer
0.1420 + 0.0027
0,8703
10576.15
LSTM
0.1008 0.0011
0.6457
8981.73
Therefore, the system stiffness can still blow up, if be t 1. On the other hand, without circular dependence, both eigenvalues are -1 and the stiffness is always regardless of the value of b.
In summary, to ensure the stability of the ODE system, complex constraints model parameters are needed for system with cycles than for system without cycles.
A9 TABULATED RESULTS
A9.1 SYNTHETIC DATA EXPERIMENTS
We attach the tabulated version of experiment results on synthetic data in Table to Table 9
A9.2
REAL-WORLD DATA EXPERIMENTS
We attach the tabulated version of experiment results on real-world data in Table 10 and Table I1
37

---

Published as a conference paper at ICLR 2026
Table Reduction method comparison (mean + standard error over 40 trials)
(Sample Size N 100) Refined Initial Graph.
Model
RMSE
Peak RMSE
ENP
GL
0.1039 + 0.0013
0.3124
757.5 .55
EGL
0.1225 = 0.0021
0.4826
1771.40
NS
0.1070 + 0.0014
0.3833
1046.95
EN
0.1225 0.0021
0.4820
1770.38
RD
0.1098 = 0.0016
0.6919
1208.05
GD
0.1114 0.0016
0.7097
1425.48
NR
0.1233 + 0.0020
0.6476
1643.83
Table Reduction method comparison (mean + standard error over 40 trials) (Sample Size N 1000) Refined Initial Graph.
Model
RMSE
Peak RMSE
ENP
GL
0.0987 + 0.0010
0.2377
727.25
EGL
0.1058 + 0.0013
0.3055
1772.03
NS
0.1045 = 0.0012
0.3492
1034.40
0.1066 + 0.0013
0.3101
1771.90
RD
0.1038 + 0.0012
0.3423
1154.30
GD
0.1040 = 0.0011
0.5046
1420.48
NR
0.1078 = 0.0013
0.3805
1641.48
Table Reduction method comparison (mean + standard error over 40 trials)
(Sample Size N 100) Comprehensive graph.
Model
RMSE
Peak RMSE
ENP
GL
0. + 0.0013
0.3892
1468 25
EGL
0.1415 0.0031
0.8537
5098.98
NS
0.1845 + 0.0039
1.1209
3565.75
EN
0.1794 0.0050
0.9047
5094.78
RD
0.1236 + 0.0022
0.5097
3905.25
GD
0.3129 0.0063
1.8013
5053.78
NR
0.1235 0.0020
0.4509
4933.95
Table Reduction method comparison (mean + standard error over 400 trials)
(Sample Size N 1000) Comprehensive graph.
Model
RMSE
Peak RMSE
ENP
GL
0.0992 + 0.0010
0.2339
1415.03
EGL
0.1287 0.0022
0.5938
5097.38
NS
0.1758 + 0.0035
0,6184
3881.53
EN
0.1287 0.0022
0.5397
5095.23
RD
0.1278 + 0.0021
0.4183
4410.73
GD
0.3415 0.0058
1,5982
5048.73
NR
0.1165 + 0.0017
0.3507
4946.93
38

---

Published as a conference paper at ICLR 2026
Table 10: Predictive performance (mean + standard error over 10 trials).
Model
RMSE
MAPE
Corr
Acc.
MNODE NR
36.19 1 0.33
0.230 0.002
0.649 0.006
0.760 1 0.004
MNODE _DK
36.58 1 0.60
0.229 0.003
0.657 0.006
0.765 0.004
MNODE HGS12
35.92 0.31
0.227 + 0.002
0.657 0.004
0.768 0.002
MNODE_ HGSI
HGSDE
3 0.40
0.229 + 0.003
0.648 + 0.005
0.764
MNODE.HOG
35.96 0.30
0.229 0.002
0.650 + 0.005
0.766
MNODE JGG
0.22
0.223 + 0.002
0.682 0.003
0.786
MNODE H
36.13 0.41
0.229 0.002
0.651 0.005
0.768
MNODE HGS13
3 0.37
0.227 0.002
0.654 + 0.007
+ 0.004
MNODE HGS23
0.228 + 0.002
0.669 + 0.003
0.769 +
MNODE EGL
35.86 0.25
0.227 0.001
0.650 + 0.005
0.760
MNODE_EN
35.93 0.29
0.227 0.002
0.651
0.760
MNODE
3 0.28
0.227 + 0.002
0.655 1 0.005
MNODE_ _T
+
0.229 + 0.002
0.655 + 0.007
0 0.003
MNODE RD
35.78
0.227 + 0.003
0.662 + 0.002
0.769
BNODE
37.08 0.25
0.260 0.002
0.666 0.003
0.759 0.003
S4D
42.91 0.39
0.283 1 0.003
0.629 0.005
4.000.
LSTM
40.69 + 0.39
0.266 1 0.003
0.666 0.003
0.733 + 0.003
TCN
11.09 0.44
0.277 1 0.003
0.672 + 0.003
0.725 4 0.002
Transformer
46.29 0.51
0.283 + 0.003
0.592 - 0.004
0.664 0.009
Table 11: Complexity and peak-error metrics. Variance shows mean + standard error over 10 trials; ENP and peak metrics are means only.
Model
Variance
ENP
Peak RMSE
Peak MAPE
MNODE NR
125.8 + 17.9
10684
189.6
1.401
MNODE DK
149.2 + 34.3
6956
183.3
1.464
MNODE_HGS12
119.9
9033
202.7
1.469
MNODE_ 1GS1
136.1 21.5
8848
177.9
1.398
ANODE_HGS?
115.4
10643
193.7
1.426
MNODE HGS
7
7551
123.4
1.222
MNODE HGS3
1
7966
167.8
1.313
MNODE HGSI3
96.6 17.2
7735
169.3
1.354
MNODE HGS23
96.4 19.9
8054
166.6
1.354
MNODE EGL
92.0 12.4
8326
169.4
1.321
MNODE EN
10..7
8548
169.0
1.239
MNODELNS
96.5
8730
161.3
1,330
MNODE GD
124.7 23.8
8861
260.6
1.349
MNODE
107.0 22.4
8955
184.8
1.322
BNODE
332.5 17.4
8596
190.9
1.541
S4D
291.4
8099
194.3
1.667
LSTM
317.8 27.9
8102
178.6
1.437
TCN
384.5 1 23.4
8261
161.6
1.803
Transformer
509.0 1 27.8
8122
210.2
1.953
39


Reviews:
[{"review_id": "hOwemCnGCH", "reviewer": "Reviewer_Sy1g", "summary": "This paper addresses a key bottleneck in hybrid neural ODEs (i.e., over-parameterization and over-fitting arising from combination of mechanistic model complexity, addition of data-driven components and sparse data) by proposing an automatic structure-learning pipeline. The method combines domain-informed graph modifications with data-driven regularization to identify essential states and interactions in mechanistic neural ODEs, thereby improving training efficiency, prediction accuracy, and robustness. The authors evaluate performance on synthetic and a real healthcare dataset, showing clear improvements over existing baselines while yielding sparser and more interpretable models. Overall, this is a strong paper. The problem is relevant and timely, the methodology is technically sound, and the experiments convincingly support the claims.", "questions": "Model class expressiveness: Does the proposed sparsification framework assume linearity in the mechanistic model parameters? If so, how does this limit applicability to nonlinear systems common in healtchare and other fields?\n\nGeneralization to larger mechanistic graphs: How does performance scale with increasing state dimensionality or interaction complexity? Any expected limitations?\n\nPractical deployment: What guidance can you provide for practitioners on choosing sparsity regularization strengths and graph-modification thresholds?\n\nTraining efficiency: Can you quantify computational savings (e.g., reduced training time, memory use) achieved through the learned sparsity?\n\nInterpretability: Are there examples where the discovered reduced model aligns with clinical or mechanistic knowledge? If so, including brief commentary would strengthen the healthcare relevance.", "limitations": "", "rating": 8, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 4, "strengths": "Originality: The paper addresses an important challenge in hybrid neural ODEs: automatic reduction of mechanistic state space and structure. The authors propose a creative combination of domain-guided graph pruning and regularization-driven sparsification, contributing to the emerging literature on structure-learning for scientific ML. Although I must admit I am not an expert in graph-based techniques so I am not able to assess the novelty of these claims, but the authors acknowledge a large body of literature.\n\nQuality: The authors present a solid methodological formulation with thoughtful integration of mechanistic priors and data-driven sparsity and strong empirical validation across both synthetic and real healthcare datasets, confirming generality. Moreover, they perform a comprehensive comparison against state-of-the-art baselines.\n\nClarity: The problem is well-motivated, and writing is generally clear and well-structured. The authors make clear statements about what this method is and what it is not (e.g., does not guarantee true model discovery). Overall, the model architecture and training procedures are communicated effectively, multiple repetitions are performed for all results to confirm robustness and claims are aligned with results.\n\nSignificance: This paper positions hybrid modeling for more deployment-ready use in data-scarce healthcare contexts, which is a domain where interpretability and robustness are essential. In addition, automated model reduction has broad relevance beyond the presented domain and opens a promising direction for robust structure learning in hybrid models.", "weaknesses": "Figures readability: Many figures contain very small titles/legends/axis labels, making interpretation difficult. Increasing font sizes and improving layout will improve accessibility and impact.\n\nInterpretation depth: While results are strong, the text describing figures and drawing insights is brief. Slightly reducing methodological narrative to provide richer interpretation would improve readability and scientific impact.\n\nLimited discussion of learned functional forms: The paper does not clearly articulate whether the recoverable mechanistic relationships must be linear w.r.t. parameters (as in LASSO-type sparsity) or whether the framework can support nonlinear mechanistic hypotheses. Clarifying this assumption and its implications for general healthcare systems would strengthen the positioning.\n\nComputational considerations not fully reported: Training time, scalability to larger mechanistic graphs, and computational savings from sparsification could be discussed more concretely in Appendix. Would this be a concern in the practical implementation and use in the types of healtcare applications of interest?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "q36HJ32Ces", "reviewer": "Reviewer_mjPg", "summary": "The paper targets an application of Hybrid ODE models in healthcare. Hybrid neural ODEs combine neural ODEs with mechanistic models to incorporate domain relational knowledge. However, mechanistic models in physiology and medicin tend to become excessively large when capturing wide ranging dynamics.  These models may contain  many latent states for a handful of observable states, and interactions among states can lead to inefficiency and overfitting. Model reduction techniques often require domain knowledge. Data-driven reduction techniques, on the other hand, are completely devoid of any domain knowledge and do not preserve mechanistic structure or constraints.\n\nThe paper proposes a method for automatic state selection and structure optimization in mechanistic neural ODEs by combining graph modification based on domain knowledge with data-driven regularization that improve predictive performance and stability while remaining mechanistically plausible. The approach combines domain-knowledge informed graph modification with a mix of L1 and L2 regularization . The graph modification bases itself on classical reduction methods to retain key topological structures. While the regularization step allows gradient based pruning during training, making the processes efficient.\n\nThe architecture of the model is an encoder-decoder on where the encoder takes in history and produces an initial state. The decoder takes the initial state and graph representation of the ODE system and future inputs evolves the state features as a Neural ODE.\n\nThe actual reduction algorithm has three steps. Step 1 merges maximally strongly connected components into supernodes with self loops. Step 2 augments the graph with short-cuts using partial transitive closure, and step 3 applies a data-driven mixture of L1 and L2 regularization.\n\nExperiments are performed on healthcare data showing improved prediction and robustness with desired sparsity.", "questions": "In the motivation of the method the paper states that the method allows domain-knowledge informed graph modification. I don’t follow how the method allows this in the steps of the method. How is this domain-knowledge used in the experiments?", "limitations": "", "rating": 4, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "The paper is very detailed regarding the method and well organized.\n\nThe method combines classical and data-driven approaches for the state reduction and sparsification problem.\n\nExperiments are done both on synthetic and real data. Ablation experiments are done for all the steps of the method.\n\nThe paper targets biomedical domains where methods must work with complex models and limited data and interpretability is significant.", "weaknesses": "Rationale of step 1 says that transforming the original graph into and RDAG reveals essential causal structure. What is the evidence for this? It also says that this improves training stability. Has this been shown in the experiments? What is the reason for this improvement in stability?\n\nHow do self-loops allow more flexible modeling of intra-component dynamics? It is said to be a key principle of hybrid modeling (line 190), but I am unsure how this is established.\n\nThe formulation of step 2 is too formal which leads to a decrease in clarity of the paper. It would be useful to start from intuition and build from there. For instance, it would be useful to first describe what is the intuitive description of a partial closure.\n\nThe introduction to the paper claims that the method allows incorporating domain knowledge reduction of the graph. As far as I can tell in steps 1 and 2, the graph reduction methods are more structural (not data-driven) rather than incorporating domain knowledge.\n\nTypos etc.\n\nline 167. “maximally strongly” should be maximal strongly …\nThere appears to be a formatting problem where line spacing has changed starting from page 5.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "oXFcL6SYGu", "reviewer": "Reviewer_egCh", "summary": "The authors propose a sparsification method for Neural ODEs parameterized by DAGs / relaxed DAGs, and benchmark on predictive accuracy for held-out samples from time-series data.", "questions": "To my understanding, only the MNODE based methods are given the prior graph determining the underlying dynamics?  In that case it seems the comparisons to the black box NODE models are somewhat uninformative, without encoding that prior into the model (even as a form of regularization).", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "The authors get good results on their synthetic and empirical benchmarks, and their graph reduction seems like a reasonable solution to the problem.  The comparison to other graph reduction methods makes the argument in favor of the current setup more convincing.", "weaknesses": "I’m a little confused by the motivation of the setup.  In most instances I’m familiar with, the incentive for parameterizing a causal model or a dynamical system with a graph is interpretability, but the authors insist that the setup isn’t identifiable and the purpose is only for predicting future states.  In that case, what is the rationale for insisting on DAG / RDAG structure at all?  One could still encode prior knowledge by constraining dependence of one state variable on others with masking.\n\nThe third part of the model seems to just be the inclusion of L1/L2 regularization, which I think is straightforward enough to not be a novel portion of the proposed method.\n\nThe highly varying y-axes across graphs makes it much much harder to gauge which performance improvements are substantial and which ones are almost negligible.  I would suggest the authors use tables in some cases so a reader can more easily tell.  For example, many bar plots in Figure 2b are somewhat misleading because the scale is so small.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "Oj5gt8R9Q4", "reviewer": "Reviewer_5nbh", "summary": "This paper presents HGS, a method for simplifying the structure of mechanistic neural ODEs. The authors point out that the mechanistic models often used may be too complex for the small datasets we typically have in fields like healthcare, which can lead to overfitting. Their proposed solution includes three steps: first to modify the model's graph structure based on its topology by collapsing cycles and adding shortcuts, and then to use regularization to learn which connections to prune. The method is evaluated on synthetic data and a real world glucose forecasting task.", "questions": "Please address the weaknesses above.", "limitations": "", "rating": 4, "confidence": 4, "soundness": 1, "presentation": 3, "contribution": 2, "strengths": "1. The work is well motivated. Trying to bridge the gap between complex mechanistic models built by experts and data-driven methods is a critical research area, and the paper does a good job of framing the problem.\n2. The idea to first perform a structural simplification of the graph before applying a data-driven pruning is quite clever. It's a nice way to inject some domain agnostic heuristics to constrain the learning problem.\n3. The set of experiments seems comprehensive. The authors have benchmarked their method against a good range of strong baselines, and the ablation study clearly shows that each part of their pipeline contributes to the final result.", "weaknesses": "My main concerns are with the real-world application, and I would need the authors to address these points before I could reconsider my score.\n\n1. My primary concern is how the model handles patient variability. The T1DEXI dataset includes 105 different people, and it's a physiological fact that glucose dynamics differ significantly between individuals. From my reading of Equation 1, the MNODE learns a single set of dynamics for everyone, i.e., an \"average patient\" model. This seems to sidestep an important challenge in this domain. Could you clarify if this is indeed the case, and if so, explain the rationale behind this choice? To be clear, this has been discussed extensively in recent years ([1-4] for example).\n2. The paper mentions that the cross validation splits were created from a \"random permutation\" of the 342 time series. To avoid data leakage, it's standard practice in clinical ML to split data at the patient level (i.e., all data from a single person stays in one fold). Could you confirm whether your evaluation followed this practice? If not, the reported performance might not accurately reflect generalization to new, unseen patients.\n3. The paper uses the term MNODE, but it seems that the functional forms of the UVA-Padova model are discarded, with only the graph structure being retained. The actual dynamics are then learned by MLPs. This is a very weak form of mechanistic prior. It would be helpful to discuss the trade-offs here.\n4. It's worth noting that in the synthetic experiments, the TCN model starts to outperform HGS on the RMSE metric as the sample size grows to N=1000. While HGS remains more robust in terms of Peak RMSE, this suggests that its strong inductive bias might become a disadvantage when more data is available. Could you comment on this trade-off?\n\n[1] Generative ODE Modeling with Known Unknowns\n\n[2] Physics-Integrated Variational Autoencoders for Robust and Interpretable Generative Modeling\n\n[3] Learning Physics Constrained Dynamics Using Autoencoders\n\n[4] CONFIDE: Contextual Finite Difference Modelling of PDEs", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)