Figure 2 illustrates the general framework of the diabetic blood glucose prediction methods for multiple patients in clinical applications. First, heterogeneous CGM signals are aggregated into a raw data stream and delivered into the learning structure. Then, the learning structure employs data cleaning and fusion techniques to reorganize the input data. The processed data are then fed into learning methods for predictive modeling or other analytical purposes. Finally, the predictive outcomes and analytical findings can support clinical diagnosis, treatment strategies, and the implementation of precision medicine.

The heterogeneity and uncertainty are inherent in CGM data samples present a significant challenge to traditional personalized prediction models, impeding their generalization to multi-patient tasks Capobianco (2017). Consequently, the alignment of sample lengths becomes a necessity. Truncation (TRA) and padding (PAD) are two prevalent techniques employed to ensure uniformity in sample lengths Chen et al. (2020); Hammad et al. (2021). Truncation, which involves reducing data to a pre-determined length, inevitably leads to the loss of potentially crucial information. On the other hand, padding, which necessitates adding artificial data to attain a specified length, risks incorporating extraneous and meaningless data, potentially reducing the accuracy of predictions.

To address these issues, we proposed a subsequence repetition (SSR) method to align and reconstruct the heterogeneous CGM data, as illustrated in Figure 3. First, the input data are individually normalized to scalar values using the min-max normalization technique in order to minimize the impact of magnitude differences. The process of normalization and de-normalization is formulated as follows: D ̃ n = D n − min D n max D n − min D n , (2) D n = D ̃ n max D n − min D n + min D n , (3) where D ̃ ( n ) and D ⁽ⁿ⁾ denote the nth normalized data and input data, respectively. max (⋅) is the maximal values of D ⁽ⁿ⁾, and min (⋅) represents the minimal values.

Subsequently, the normalized time series data are segmented into fixed-length subsequences, each representing a continuous time period consisting of T steps (i.e., window size). The generations after segmentation can be denoted as D ̃ = { D ̃ # ( 1 ) ∈ R B ( 1 ) × T , … , D ̃ # ( N ) ∈ R B ( N ) × T } . D ̃ train and D ̃ test are used to denote the training set and test set, respectively, after data splitting. In the test set, we maintain a consistent number of segments across all samples. Nevertheless, due to the uncertainty of sample length, the number of split segments within the training set remains inconsistent across the samples. Hence, we employ SSR to achieve alignment of segments in the training set, which can be obtained by: (4) where s denotes the prediction steps, and T is the window size. The known subsequence segments are sequentially replicated until the total segment count equals β, where β represents the maximum number of segments within the training set. The reconstructed training data can be denoted as D ̃ # train = { X # train , Y # train } and the test data is D ̃ # test = { X # test , Y # test } .

The SSR method aligns the heterogeneous input data by utilizing the known data samples to fill the null subsequences. This method might increase redundant subsequence information, but it avoids losing the original temporal dependencies of the time series data. Both padding and truncation can easily affect the distribution of raw data.

The model architecture of the proposed HETER is illustrated in Figure 4. In HETER, we initially calculate the distance between each heterogeneous CGM record in the training set, then construct a sparse relation graph (SRGraph) based on the obtained distance matrix. SRGraph is capable of capturing the intricate relationships among patients by arranging the input data in a graph structure. This is particularly crucial when considering the heterogeneity of individuals’ blood glucose levels influenced by varying lifestyle choices, treatment regimens, and genetic factors. Additionally, SRGraph excels in uncovering latent correlations between samples or periods across both spatial and temporal dimensions. This makes SRGraph particularly suitable for multiple patients’ BGP, where correlations may exist but are not explicitly evident in the data.

Due to the inconsistent sample length among patients, the distance cannot be calculated on a point-to-point basis. Consequently, we employ the dynamic time-warping (DTW) approach to align the sample lengths of the raw input data Li et al. (2020). This method is a widely used approach for measuring similarity between time series Li (2021). It aims to minimize the distance between two time series while aligning their sample lengths. The process can be expressed as follows: γ i , j = d i s t a n c e D i , D j = min Λ ∑ s 1 , s 2 ∈ Λ d i s t a n c e d i , s 1 , d j , s 2 , (5) where γ _i,j represents the calculated distance between samples i and j. Λ = [λ ¹, … , λ ^κ ] denotes the optimal alignment path between two samples. The alignment path λ = (s ₁, s ₂) must satisfy 1 ≤ s ₁ ≤ L _i and 1 ≤ s ₂ ≤ L _j , where s ₁ and s ₂ represent the time points of the samples, and L _i and L _j denote the sample lengths of D ⁽ⁱ⁾ and D ^(j), respectively. The path begins at λ ¹ = (0, 0) and ends at λ ^κ = (L _i−1, L _j−1). Both d _i ∈ D ⁽ⁱ⁾ and d _j ∈ D ^(j) represent observations from the CGM record. The distance matrix for the input time series can be expressed as Γ = [γ _1,2, … , γ _N,N−1].

Then, the formal definition of the proposed SRGraph can be described below.

Definition 1 (SRGraph). The SRGraph is formulated as G = ( V , E ) , where V ∈ R N is the set of nodes (i.e., CGM samples), and E ∈ R 2 × N × K is the set of edges (i.e., the relationships between CGM samples). N is the sample size and is equal to the number of nodes in the graph. K is a hyperparameter that denotes the number of items we aim to retrieve in the top positions of the distance ranking.

Definition 2 (Node Neighborhood). Let v _i , v _j ∈ V denote nodes in the graph, and E i , j = ( v i , v j ) to denote an directed edge pointing from v _i to v _j . The neighborhood of a node v _i is represented as N ( v i ) = { v j ∈ V | ( v i , v j ) ∈ E } .

Definition 3 (Adjacency Matrix). The adjacency matrix is a mathematical representation of a graph. This study sorts the distances and selects the top-K minimal distances for each sample to construct the adjacency matrix. This matrix is denoted as A ∈ R N × K . If ( v i , v j ) ∈ E , then A _i,j > 0, and if ( v i , v j ) ∉ E , then A _i,j = 0.

In general, CGM samples are treated as nodes in the graphs. The edges of the graph are constructed based on the distance relationships between these samples. With the information about nodes and edges, we generate the graph adjacency matrix, which is subsequently learned by the graph structure learning module.

The proposed graph learning module intends to integrate the information of a node with that of its neighboring nodes to capture the temporal dependencies in a graph. Two stacked graph convolutional networks (GCNs) Zhang et al. (2019) are employed in the proposed HETER to embed the information of the SAGraph. This can be expressed as follows: R g 1 = W g 1 ⋆ X # = M ̃ − 1 2 A ̃ M ̃ − 1 2 X # W g 1 + b g , (6) R g 2 = W g 2 ⋆ D r o p o u t φ R g 1 , p g , (7) where symbol ⋆ indicates the graph convolution operation. X # is the input for proposed model. M ̃ ∈ R T × T denotes the diagonal matrix. The element M ̃ i , i represents the degree of node i, corresponding to the number of edges connected to node i. A ̃ = A + I u ∈ R N × ( K + 1 ) denotes an adjacency matrix with self-loops. A is the adjacency matrix of SAGraph and I _u represents a unit matrix. R g ( 1 ) is the embedding graph representation. R g ( 2 ) is the output graph representation. p _g is dropout probability. φ(⋅) is rectified linear unit (ReLU) activation. W _g represents the learnable weighting matrix. b _g is the learnable bias.

The potential topological characteristics between patients are extracted by the graph learning module Wu et al. (2020). The representation enhancement module is designed to exploit further the key information in both the temporal and spatial dimensions.

First, a convolution component is employed to filter out the useless information in spatial dimensions, the Dropout technique also be integrated to prevent overfitting: R c o = D r o p o u t φ W c o * R g 2 , p c o , (8) where W _co is the learnable convolution kernel, R c o represents the output from the convolution component, and * denotes the traditional convolution operation. p _c is the dropout probability for the convolution component.

A limitation of graph convolutional networks is their inability to capture the continuous dynamical patterns from input CGM data. This is because the information propagation of graph convolutional networks occurs merely through nodes, ignoring the continuity of temporal features within nodes. To address this issue, temporal attention (TA) and a recurrent learning component are incorporated into the proposed HETER. Temporal attention is applied to highlight the key time steps in a period: S t a = s o f t m a x W t a 1 φ W t a 2 R c + b t a 2 + b t a 1 , (9) where S t a denotes the output attention score of TA. W _ta and b _ta are the learnable weighting matrix and bias, respectively.

HETER utilizes the attention score to reweight the output of the convolution component, which is then delivered to the gated recurrent unit (GRU)-based recurrent learning component to further extract the temporal dynamics. The information propagation processing in the recurrent learning component can be formulated as: r t = σ W r h t − 1 ; R c o ⊙ S t a + b r , (10) u t = σ W u h t − 1 ; R c o ⊙ S t a + b u , (11) c t = σ W c r t ⊙ h t − 1 ; R c o ⊙ S t a + b c , (12) h t = 1 − u t ⊙ h t − 1 + u t ⊙ c t , (13) where h _t is the hidden state at t steps of GRU unit. u _t , r _t , c _t denotes the update gate, reset gate, and candidate hidden state of recurrent learning component, respectively. σ(⋅) indicates the sigmoid activation function. [;] represents the concatenation operation, and ⊙ is the element-wise multiplication operation.

The output layer of the proposed HETER method consists of two multilayer perceptron networks, which are designed to encode the result into the desired output shape. To avoid the issue of vanishing or exploding gradients induced by excessive nonlinearity, the model’s linear input representation is added to the product of the output layer to yield the final prediction results: R o = S M L P T M L P D r o p o u t h t , p o , (14) O ̂ # = R o ⊕ W o X # + b o , (15) where R o is the output from the generation layer, and O ̂ # represents the predictive values. p _o is the dropout probability, and ⊕ denotes the element-wise addition operation. W _o and b _o are the weight matrix and bias, respectively, used to generate the linear representation of the model input. The final results Y ̂ # are obtained after de-normalization by Eq. (3). Here, SMLP (⋅) and TMLP (⋅) correspond to Spatial-LP and Temporal-LP, respectively.

To verify the effectiveness of the proposed HETER, we select the following five methods as baselines for comparison. These methods have already achieved promising performance in personal diabetic blood glucose prediction or universal time series tasks.

1) Vanilla LSTM Rabby et al. (2021) is a fundamental long short-term memory (LSTM) model that utilizes its internal gating mechanisms to capture temporal dynamics.

We employ the mean absolute error (MAE), mean absolute percentage error (MAPE), and root mean square error (RMSE) to provide a comprehensive evaluation the performance of the predictive model. The MAE provides a straightforward metric of average error magnitude. The MAPE offers a scale-independent error measure in percentage terms, which is critical given the varying blood glucose levels across different contexts. The RMSE, emphasizing larger errors, assumes importance in our study, particularly considering the potential health consequences of substantial deviations in blood glucose levels. Moreover, these criteria are frequently adopted metrics for evaluating the accuracy of blood glucose predictions Xie and Wang (2020); Yang et al. (2023), which can be expressed as follows.

1) Mean absolute error

The experiments use the most recent 36 h (144 data points) of data for testing, while the remaining data were allocated for training. All experiments were conducted five times to improve the reliability of the results. The optimal hyperparameters were determined using the grid search method. The models were trained using the Adam optimization algorithm Kingma and Ba (2015), with mean squared error (MSE) employed as the loss function. The training epochs and learning rate are adjusted to approach their optimal states for each method. Supplementary Table S2 provides the detailed settings for the hyperparameters of each method. We implement the models using Porch (v.1.12.1) Paszke et al. (2019). The experiments for baseline methods and the proposed method were carried out on a server equipped with Intel(R) Xeon(R) Gold 6226R CPU (2.90 GHz) with 128G memory and were accelerated by two NVIDIA RTX A6000 GPUs.