Calculating the Mean Amplitude of Glycemic Excursions from Continuous Glucose Data Using an Open-Code Programmable Algorithm Based on the Integer Nonlinear Method

The mean amplitude of glycemic excursions (MAGE) is an essential index for glycemic variability assessment, which is treated as a key reference for blood glucose controlling at clinic. However, the traditional “ruler and pencil” manual method for the calculation of MAGE is time-consuming and prone to error due to the huge data size, making the development of robust computer-aided program an urgent requirement. Although several software products are available instead of manual calculation, poor agreement among them is reported. Therefore, more studies are required in this field. In this paper, we developed a mathematical algorithm based on integer nonlinear programming. Following the proposed mathematical method, an open-code computer program named MAGECAA v1.0 was developed and validated. The results of the statistical analysis indicated that the developed program was robust compared to the manual method. The agreement among the developed program and currently available popular software is satisfied, indicating that the worry about the disagreement among different software products is not necessary. The open-code programmable algorithm is an extra resource for those peers who are interested in the related study on methodology in the future.


Introduction
Clinical researches have suggested that high glycemic variability may cause more serious damage to the body than high level stable blood glucose [1], which relates to the development of diabetic complication [2][3][4][5][6] and the increase of mortality in critically serious patients without diabetes [7,8]. In such circumstances, how to quantitatively evaluate the glycemic variability in diabetic blood glucose monitor is essential for the clinical diagnosis and treatment. Several indexes have been proposed for the quantitative evaluation of the glycemic variability, such as MBG, SDBG, IQR, LAGE, M-value, %CV, J-index, IGC, GRADE, MODD, LBGI, HBGI, ADRR, TI, LI, PGS, CONGA, and MAGE [9][10][11][12][13][14][15][16]. Presently, clinicians and researchers trend to choose MAGE as the preferred index [17,18], making it a "popular standard" in the quantitative evaluation of the within-day glycemic variability.
MAGE is an arithmetic average of either the upward or downward of all glycemic excursions exceeding the threshold (standard deviation of blood glucose (SDBG) obtained from all blood glucose concentrations within 24-hour period); the direction of the calculation is determined by the first countable excursion [19][20][21]. Following the original definition by Service et al. [19], clinicians and researchers applied the "ruler and pencil" graphical approach to calculate MAGE. However, this kind of manual approach is time-consuming and errorprone, when dealing with huge amount of data typically characterized by 288 observations in 5 minutes apart over a 24hour period, generated from continuous glucose monitoring (CGM). Thus, the development of computer-aided programs to calculate MAGE becomes an urgently needed task.

Computational and Mathematical Methods in Medicine
To address the urgent need above, several automatic programs have been developed [22][23][24][25], making great contribution to the automatic calculation of MAGE. However, it was pointed out that the agreement among them is poor [26]. Therefore, more studies are required to explore the automatic calculation of MAGE. Currently available programs can be roughly divided into two groups. One group has detailed descriptions about algorithm with graphical display, such as the automated algorithm described by Baghurst [22] and the computer software described by Fritzsche et al. (Fritzsche) [23]. The other group only provides executable software for automatic calculation, without detailed descriptions of the algorithms used in the software, such as the web-based application "GlyCulator" [24] and the Excel-based workbook "EasyGV" [25]. However, to the best of our knowledge, all of these methods did not provide programmable open codes, which are important resources required for peers to implement more studies on methodology in this field. Therefore, in this report, we developed a mathematical algorithm based on integer nonlinear programming method. Following the proposed mathematical method, a computeraided program named MAGECAA v1.0 was developed. The code of our program is open; if peers are interested in it, please contact xuefeiyu@smu.edu.cn for downloading. To validate the developed program, comparison study was implemented using blood glucose CGM datasets obtained from T1D, T2D, and gestational diabetes patients against the manual method (MAGEo) and other currently popular software products.

The Proposed Mathematical
Algorithm. Let ( ) ( = 1 , 2 , . . . , ) represent the discrete blood glucose values obtained in CGM; then is a discrete function defined in time set { 1 , 2 , . . . , }. Let SDBG be the standard deviation of ( ) ( = 1 , 2 , . . . , ). A graph depicting the glycemic variability can be formed by connecting all the discrete values of function . When the difference of a peak and an adjacent nadir exceeds SDBG, the corresponding peak is labeled as a valid peak. The key point for the calculation of MAGE is to correctly count the valid peak or nadir.
Amplitude is the difference of functional values in a peak and a nadir of the function graph. A valid amplitude is labeled and counted when it is bigger than SDBG. To compute MAGE, these valid amplitudes of function should be firstly searched. From the mathematical perspective, a peak or nadir should be an extreme point of the function. These extreme points related to valid amplitudes are called valid extreme points. The MAGE computation problem can be solved by calculating all valid extreme points.
Thus, the above MAGE computation problem can be transformed to an integer nonlinear programming (INLP) problem.
Suppose that * is the maximal with respect to the optimal solution of function (2), and then the optimal solution { * 1 , * 2 , . . . , * * } represents the valid extreme points to compute MAGE.
If the extreme point * 1 is the local minimum point, then If the extreme point * 1 is the local maximum point, then Computational and Mathematical Methods in Medicine 3 Now, our object is to get the optimal solution { * 1 , * 2 , . . . , * * } of function (2). If the INLP problem is solved by enumeration algorithm, the amounts of different subsequences { 1 , 2 , . . . , } should be ( ). For 3 ≤ ≤ , the amounts of all subsequences are It is difficult to directly solve function (2) when is large. A faster optimization algorithm is usually used to solve the above INLP problem. For example, a penalty function algorithm [27,28] was used here to transform the above INLP problem (2) to an unconstrained optimization problem as follows. Set The new unconstrained integer optimization problem becomes arg min 1 , 2 ,..., where 3 ≤ ≤ and and ( = 1, . . . , −1) are penalty coefficients and tend to be +∞.
When the optimum value * > 0, problem (9) has no optimal solution. Differential evolution (DE) algorithm, a faster optimization algorithm proposed by Storn and Price [29], is a simple but powerful population-based stochastic search technique to solve global optimization problems over continuous domains. Many researchers modified DE algorithm to improve its performance when it was applied to a specific problem [30][31][32][33]. The idea of a modified DE algorithm proposed by Lin [34] to solve mixed-integer nonlinear programming problem is used here.
DE searches for a global optimal point in andimensional hyperspace. Let be the -dimensional search space of the INLP problem under consideration. The DE evolves a population of NP -dimensional individual vectors, that is, solution candidates, N = ( 1 , . . . , ) ∈ , = 1, . . . , NP, from one generation to the next. The evolution begins with a randomly initialized population of -dimensional integer parameter vectors in space {1, 2, . . . , } . In each vector, integer parameters are sorted in an ascending order. Each vector forms a candidate solution to the unconstrained optimization problem. At each generation , DE employs the mutation and crossover operations to produce a trial vector U for each individual vector N , also called target vector, in the current population.
The details of the employment of the DE algorithm to solve the above INLP problem are as follows from (a) to (d).
(a) Initialization. A randomly initialized population is created to cover the entire search space uniformly as in the following form: where is a random number in the range [0, 1] and NINT [B] is expressed as the nearest integer vector to real vector B.
(b) Mutation Operation. For each target vector N at generation , randomly sample three other individuals N 1 , N 2 , and N 3 from the same generation, where 1 , 2 , and 3 are random and mutually different integers generated over the range [1, ], which should be different from the current trial vector's index . Then an associated mutant vector V = (V 1 , . . . , V ) can be generated by using strategy: where is a factor in [0, 1] for scaling differential vectors.
(c) Crossover Operation. The crossover operation is applied to each pair of the generated mutant vector V and its corresponding target vector N to generate a trial vector U = ( 1 , 2 , . . . , ).
where CR ∈ [0, 1] is a crossover constant that is determined by users. rand is a randomly chosen index in [1, ] which ensures that gets at least one parameter from V . The integer parameters of target vector N are also sorted from small to big.
(d) Selection Operation. The trial vector U is compared to its corresponding target vector N using the greedy criterion to decide whether a member of generation + 1 existed or not. If vector U yields a smaller cost function value (U ) than (N ), then N +1 is set to U ; otherwise, the old vector is retained. The operation is expressed as follows:  The above (b) to (d) steps are repeated till the evolution times arrived to certain number (general = 200); from the last evolutionary generation, the individual vector N with the smallest value of the objective function is the optimal solution of the problem of the current -dimensional extreme point combination. The algorithm will search from 3 to till the objective function obtains the minimum value so as to obtain the optimal extreme point combination { * 1 , * 2 , . . . , * * } and then use formulas (5)∼(6) to calculate the MAGE value. Figure 1 shows the entire algorithm. was developed. The MAGE calculation program can be described as a process that selects valid extreme points from a time-ordered set of glucose concentrations whose adjacent differences are all greater than the threshold (typically 1 SDBG obtained from 24-hour period blood glucose concentrations). It can be summarized as finding the optimum vector combination solution of the valid extreme points and using INLP to establish the mathematic method, which can be solved by differential evolution (DE) algorithm. Once all valid extreme points of countable excursions have been identified, the MAGE is determined by MAGE + or by MAGE − , depending on the direction of the first countable excursion. For more in-depth understanding, it depends on the first valid extreme point of the vector combination, because that point indicates the direction. In addition, the average of both MAGE + and MAGE − , designated as MAGEa, is also calculated. MAGE-CAA v1.0 is based on INLP and has several different outputs: SDBG, MAGE + , MAGE − , MAGE, MAGEa, and so forth. Besides it also needs plot to show all valid extreme points joined by straight lines; MATLAB (MathWorks5, USA) is chosen as the programming environment accordingly. Generally, some data points could not be extreme point according to the mathematical definition, like points shown in Figure 2, which shows turning points and + 1 which are equal in values (not being extreme points) but there is an amplitude from peak point − 1 to nadir point , so we combine points and + 1 as one point + 1 or remove point and then let point + 1 be a local maximum point. If the points appear to be opposite, then point +1 should be a local minimum point.
MAGECAA v1.0 consists of the following major modules: (1) import CGM data and calculate the SDBG as the threshold; (2) identify all extreme points; (3) find the optimum vector combination solution of valid extreme points; and (4) display the calculated parameters and plots.

Statistical Analyses.
Spearman's correlation analysis was applied in evaluating the relationship between the MAGE values obtained by different methods with respect to the same patients. Bland-Altman plots were used to represent the agreement of the methods [35,36]. < 0.01 was considered significant.  Figure 4, Spearman's correlation analysis identified that there is a highly significant linear correlation between MAGEc and MAGEo ( = 0.998, resp.; < 0.01), and the mean difference found in Bland-Altman plot was −0.03 ± 0.21, which was statistically significantly small. This result indicated that the developed program was robust compared with the manual method.

Agreement of MAGECAA v1.0 with Fritzsche and EasyGV.
To evaluate the agreement of MAGECAA v1.0 with two currently available popular software products, that is, Fritzsche and EasyGV, we did pairwise comparison of MAGE among them. Table 1 shows the correlation coefficients for MAGE calculation between the calculators, which ranged from 0.926 to 0.987 ( < 0.01 for all). Figure 5 shows the Bland-Altman plot among the three software products; the dashed lines represent the 95% confidence limit of the differences between the two methods, and the solid line demonstrates that the mean difference between the methods is close to 0, indicating that little difference existed among them. The mean differences are 0.14 ± 0.55 (MAGEc versus Fritzsche), 0.52 ± 1.17 (MAGEc versus EasyGV), and 0.38 ± 1.12 (Fritzsche versus EasyGV), showing the good agreement between these three software products. To explore whether MAGEa may be a more useful index than MAGEc which depends on MAGE + or MAGE − , the relationship between MAGEc and MAGEa was tested by using all the 248 CGM measurements via Spearman's correlation analysis and Bland-Altman analysis. As shown in Figure 6, the correlation coefficient was = 0.993 ( < 0.01) and the mean difference was −0.05 ±0.48. The results showed that little difference was observed between using MAGEa and MAGEc, indicating that MAGEa may be a useful index for evaluating glycemic variability.

Discussions
We developed a computer-aided open-code program named MAGECAA v1.0 based on INLP algorithm for automatic calculation of MAGE. Compared with the existing methods, the proposed novel method turns to search the optimal solution of the combination of extreme points from overall CGM measurements instead of searching adjacent extreme points from the beginning to the end step by step. As for the computational time, if used for one person, the proposed method is comparable with currently available methods; if used for batch calculation, the proposed method is more powerful. The programmable open codes are useful for the study of methodology for automatic calculation of MAGE. The comparison study using MAGECAA v1.0 against manual method indicated that the agreement is satisfied. The pairwise comparison study between MAGECAA v1.0 and two other available software products, Fritzsche and EasyGV, based on Spearman's correlation analysis and Bland-Altman plots, demonstrated that the agreements between them met the requirement. Our study showed that the worry about the disagreement among the currently available popular software products is not necessary, which is different from the proposal by Sechterberger et al. [26]. The reason may be the increased amount of CGM data used in our research.
The MAGE value depends on MAGE + or MAGE − , following the direction of the first accountable glucose excursion. Considering the fact that only one direction of glucose excursion is adopted in current popular software, unavoidably resulting in the omission of the other directions of valid  Computational and Mathematical Methods in Medicine excursions, we implemented an extra experiment in which data from both directions of the valid glycemic excursions are utilized. MAGEa was used to represent the mean of MAGE + and MAGE − for the calculation of MAGE. As shown by our data, a close linear correlation between MAGEa and MAGEc was observed, indicating that the difference between MAGE + and MAGE − is significantly small. Therefore, we proposed that MAGEa might be another suitable parameter to quantify glycemic variability.
To conclude, an open-code software program named MAGECAA v1.0 for automatic calculation of MAGE based on a mathematical algorithm has been proposed and evaluated. The programmable open codes are useful for further methodology study in the future. The comparison study indicated that the agreement among the proposed software and existing software is satisfied, and the worry about the disagreement among currently available different popular software products is not necessary.