Understanding how biomarker distributions evolve over time is a central challenge in digital health and chronic disease monitoring. In diabetes, changes in the distribution of glucose measurements can reveal patterns of disease progression and treatment response that conventional summary measures miss. Motivated by a 26-week clinical trial comparing the closed-loop insulin delivery system t:slim X2 with standard therapy in children with type 1 diabetes, we propose a probabilistic framework to model the continuous-time evolution of time-indexed distributions using continuous glucose monitoring data (CGM) collected every five minutes. We represent the glucose distribution as a Gaussian mixture, with time-varying mixture weights governed by a neural ODE. We estimate the model parameter using a distribution-matching criterion based on the maximum mean discrepancy. The resulting framework is interpretable, computationally efficient, and sensitive to subtle temporal distributional changes. Applied to CGM trial data, the method detects treatment-related improvements in glucose dynamics that are difficult to capture with traditional analytical approaches.
Figures & tables
Figure 1: Individual participant analysis for the bivariate model (glucose and its first derivative) using K=5 Gaussian components. Each row corresponds to one participant (IDs 20, 58, and 82). Left: raw CGM time series showing glucose concentration (mg/dL) over the observation period. Middle left: estimated weight trajectories α(t)=(αk(t))k=15 learned by the neural ODE, representing the evolution of the mixture proportions over normalized time t∈[0,1] . Middle right and right: contours of the fitted bivariate Gaussian mixture density at the initial and final times, summarizing the joint distribution of glucose level and its rate of change.
Figure 2: Comparison of weight trajectory dynamics between Treatment ( green ) and Control ( red ) groups for the d=2 bivariate model with K=5 mixture components. Each panel shows the evolution of component weights αs(t) between weeks 20–26 over normalized time t∈[0,1] . Group means are shown as thick dashed lines. The shaded bands represent a statistical envelope around the mean (e.g., between the 5th and 95th percentiles).
Figure 3: Predicted glucose density distributions between weeks 20 and 26, comparing Treatment ( green ) and Control ( red ) groups for the d=2 model with K=5 mixture components. The marginal density over glucose concentration and speed is computed by drawing samples from each participant’s GMM, estimating the sliced-Wasserstein barycenter in 2D, and converting the barycenter samples into a smooth density via a Gaussian KDE on the grid. The top row corresponds to the Control group and the bottom row to the Treatment group. In both rows, the first column shows the initial distribution at week 20, the second column displays the final distribution at week 26, and the third column presents the difference between the two distributions.
Figure 4: Wild Bootstrap MMD test p -values comparing Treatment vs. Control groups over time for the d=2 model with K=5 mixture components. The dashed black line indicates the significance threshold α=0.05 . The colors correspond to the different components: blue (component 1), orange (component 2), green (component 3), red (component 4), and violet (component 5).
Figure 5: Quantile curves (median and 25%–75% bands) of the change in GMM mixture weights for each of the K=5 components over time, relative to their initial value, for Treatment ( green ) and Control ( red ) groups. Each panel shows the temporal evolution of a component’s weight deviation from baseline.
Appendix figures & tables8 assets
Supplementary material from the paper’s appendix.
Appendix
Figure 6: Pointwise L2 -error over time for d=1 . We compare our model , KDE , MAF , and GAMLSS . Curves represent averages over 100 independent runs. Errors corresponding to the discrete-time MMD stage are shown as ■ . L2 -errors are approximated by Monte Carlo integration.
Figure 7: Pointwise L2 -error over time for d=10 . We compare our model , KDE and MAF . Curves represent averages over 100 independent runs. Errors corresponding to the discrete-time MMD stage are shown as ■ .
Figure 8: Comparison of weight trajectory dynamics between Treatment ( green ) and Control ( red ) groups for the univariate model ( d=1 ) with K=3 mixture components. Each panel shows the evolution of component weights αs(t) between weeks 20–26 over normalized time t∈[0,1] . Group means are shown as thick dashed lines. The shaded bands represent a statistical envelope around the mean, where darker shading indicates a higher density of trajectories.
Figure 9: Quantile curves (median and 25%–75% bands) of the change in GMM mixture weights for each of the K=3 components over time, relative to their initial value, for Treatment ( green ) and Control ( red ) groups in the univariate model ( d=1 ). Each panel shows the temporal evolution of a component’s weight deviation from baseline.
Figure 10: Predicted glucose density distributions for the univariate model ( d=1 ) with K=3 mixture components. The curves represent the marginal density over glucose concentration, computed as the Fréchet mean (1D Wasserstein barycenter) across individuals within each group. The top row compares the Treatment ( green ) and Control ( red ) groups at the initial observation (week 20, left) and final observation (week 26, right). The bottom row displays the difference between the final and initial distributions for the Control group (left) and the Treatment group (right).
Figure 11: Wild Bootstrap MMD test p -values comparing Treatment vs. Control groups over time for the univariate model ( d=1 ) with K=3 mixture components. The dashed black line indicates the significance threshold α=0.05 . The colors correspond to the different components: blue (component 1), orange (component 2), and green (component 3).
Univariate ( d=1 )
Param.
s=1
s=2
s=3
μs
109.38
178.38
275.16
σs2
421.86
613.47
1692.82
Bivariate ( d=2 )
Param.
s=1
s=2
s=3
s=4
s=5
μs
[90.74−0.0081]
[129.24−0.0152]
[173.03−0.0165]
[229.530.0902]
[305.20−0.1327]
Appendix
Table 1: Calibrated global statistics for univariate and bivariate data.
Proposed Two-Stage Estimator
Step 1: per-time MMD mixture fit
Step 2: time-series neural ODE fit
Mixture components ( K )
5
Vector field MLP
2 layers, width 64
Kernel bandwidth ( σ )
median heuristic
Integration horizon ( T )
1.0
Ridge regularization ( λ )
10−2
Solver step size
0.01
Adam learning rate
10−3
Adam learning rate
10−3
Iterations / inner grad. steps
50 / 20
Max epochs
2000
Appendix
Table 2: Hyperparameters for the proposed estimator and baselines.
Gaussian Mixture Models represent data distributions as weighted combinations of Gaussian components for clustering, density estimation, and multimodal prediction.
University of Virginia, School of Data Science, Charlottesville, VA 22903, USA · University of Virginia, Department of Pediatrics, Charlottesville, VA 22903, USA