A computational model of the human glucose-insulin regulatory system☆

Keh-Dong Shiang, Fouad Kandeel

aDivision of Biostatistics, Department of Information Sciences, City of Hope National Medical Center, Duarte, CA 91010-3000, USA bDivision of Hematopoietic Stem Cell and Leukemia Research, City of Hope National Medical Center, Duarte, CA 91010-3000, USA cDivision of Diabetes, Endocrinology and Metabolism, City of Hope National Medical Center, Duarte, CA 91010-3000, USA Received 3 Aug 2010

INTRODUCTION

Diabetes mellitus is a chronic metabolic disorder characterized by abnormally high urine and blood glucose levels (i.e., hyperglycemia) due to insufficient insulin levels. Based on the statistics of American Diabetes Association (ADA)[1], approximately 23.6 million people, or 7.8% of the population, in the USA are afflicted with this disease. While an estimated 17.9 million people have been diagnosed with diabetes, regrettably, 5.7 million people (or nearly one quarter) are unaware that they have the disease. The total annual economic cost of diabetes in 2007 was estimated to be $174 billion. Diabetes can cause serious health complications including blindness, heart disease, kidney failure, stroke, nerve damage, and lower extremity amputations. According to the statistical report from Centers for Disease Control and Prevention (CDC), diabetes is the sixth leading cause of death in the USA[2]. Consequently, diagnosis, treatment, control and prevention of diabetes, are extremely critical in the current medical era.

In a normal subject, the beta-cells (β-cell) in the pancreas release insulin in response to rises in the level of glucose in the blood, which results in the storage of this source of energy as glycogen in the liver. Type Ⅰ diabetes, also called juvenile or insulindependent diabetes mellitus (IDDM), often manifests in childhood and may result from autoimmune destruction of insulin-producing β-cells of the pancreas. Thus, insulin hormone can no longer be produced. This type of diabetes is fatal without treatment with exogenous insulin to replace the missing hormone or providing patients with a functional replacement for the destroyed pancreatic β-beta cells, such as pancreas or islet-cell transplantation. Type Ⅱ diabetes (formerly called non-insulin-dependent diabetes mellitus, NIDDM, or adult-onset diabetes), a more widespread metabolic disorder, is primarily characterized by insulin resistance, relative insulin deficiency and hyperglycemia. Some cases of type Ⅱdiabetes also appear to be an autoimmune disease where the immune system attacks the β-cells, decreasing the function of producing insulin, while other type Ⅱdiabetes cases may simply result from excessive body weight that strains the ability of the β-cells to produce sufficient insulin. However, in both type Ⅰ and type Ⅱ cases, the human body loses its ability to regulate blood sugar, which causes a significantly negative effect on the patients′ quality of life or even be potentially fatal.

It is a common knowledge that blood glucose concentration in normal humans is maintained within a precise and stable range. Many external and internal factors affect the level of blood glucose such as food intake, rate of digestion, excretion, exercise, sleep, and psychological state. These individual or combinational influences constantly alter the physiological processes that regulate plasma glucose level. For instance, if blood glucose is elevated, after a regular meal (i.e., post-prandial), certain cells in the pancreatic islets of Langerhans named β-cells, will release the insulin hormone. The secreted insulin, then, leads to the uptake of glucose from the blood into the liver and other cells, such as muscle cells. Thus, blood glucose level will eventually go down to the normal range. On the other hand, blood glucose level may decrease imminently due to muscular activity, particularly when food intake is confined. This reduced level of blood glucose is immediately recognized by other sensitive pancreatic cells, the alpha-cells (α-cells). These cells then release glucagon that act on the cells of the liver to initiate the release of glucose. This results in blood glucose level elevating back to the normal range. Briefly, these islet-cell arguments establish the fact that the capacity to lower blood glucose depends on the responsiveness of the pancreatic beta-cells to glucose and the sensitivity of the glucose utilized by tissues to the released insulin. Thus, both pancreatic β-cell responsiveness and insulin sensitivity contribute to glucose tolerance. Low glucose tolerance in lean humans could be associated with diminished β-cell response to glucose, whereas low glucose tolerance in obese people could be associated with decreased insulin sensitivity. Furthermore, a shortage of plasma insulin and low glucose tolerance, resulting in a serious inability to lower blood glucose, will cause insulin resistance, which is the key symptom underlying the potential development of diabetes. However, to tackle diabetes disease and obesity problems, clinicians and researchers are now turning to mechanism-based mathematical models to reach quantitative diagnoses of glucose intolerance and insulin resistance, and also to predict the likely outcomes of therapeutic interventions. Their ultimate goal is to develop a mathematical model that can be used to accurately predict the outcomes and most successful treatment options for people who have diabetes.

Speaking of theoretical solutions to diabetic problems, we should mention the term, "mathematical model", which is a representative depiction of a real system via mathematical tools in these medicines′ golden years. The fundamental nature of a good mathematical model must be simple in design and exhibit the basic properties of the real system that we are attempting to simulate and understand. All welldeveloped models should be validated and tested against empirical data. In a practical sense, the quantitative comparisons of the model to the real system should lead to an improved mathematical model. The successful model can be applied to suggest the corresponding experiment to highlight a particular aspect of the weakness or problem, which may improve the method of data collection or the procedure of experimental processes. Thus, modeling itself is an evolutionary process, which is a evolving procedure in which something changes into a different but better form. Similarly, developing and using a successful mathematical model will guide us to learn more about certain simulating or existing processes rather than finding an entirely actual state of the system.

During the last few decades, a massive range of mathematical models, computer algorithms and statistical methods have been proposed in order to understand different aspects of diabetes, such as glucose metabolism, insulin kinetics, β-cell mass, and the glucose-insulin regulatory system. Several reviews have been devoted to mathematical models and diabetic disease[3-8]and are worthwhile to be referenced. Other than those reviewing journal articles, a pioneering work on modeling the glucose-insulin regulatory system and its ultradian insulin secretory oscillations can be traced back to Bolie[9]. In this pioneering study, a system of glucose-insulin regulation in terms of coupled differential equations of feedback was analyzed with the so-called critical damping criteria of a selfregulating feedback system (i.e., servomechanism theory). The secretion of insulin in the glucose-insulin endocrine metabolic system occurs in an oscillatory manner over a range of 50-150 min and is usually referred to as ultradian oscillations[10]. In 1965 and 1969, Ackerman et al.[11,12]pioneered the concept of offering their full exploratory work for the blood glucose regulatory response to the glucose tolerance test (GTT), which was governed by two coupled differential equations. In the following sections, we will introduce their conceptually illuminating model in greater detail, and also develop our computational model, which will be validated by using their model equations and other published experimental data and results.

In order to determine whether or not a patient has pre-diabetes or diabetes, health care providers usually conduct a fasting plasma glucose (FPG) test or a GTT. The ADA recommends the FPG test because it is easier, faster, and less expensive. Therefore, in the following, several quantitative assessment methods are briefly introduced. Unquestionably, due to the advantage of skipping GTT procedures, fasting plasma glucose level is simpler and quicker to measure, and its measurement is more acceptable to patients than any glucose tolerance test.

The homeostasis model assessment (HOMA) is an index of insulin resistance (IR). It was developed by Matthews et al.[13,14]and derived from the product of FPG and fasting plasma insulin (FPI) divided by a constant, 22.5.

where glucose is given in mg/dL and insulin is given in μU/mL. In this equation, the constant 405 should be replaced by 22.5 if the glucose is reported in mmol/L.

Since hepatic glucose production (HGP) is the main determinant of FPG concentration, and FPI concentration is the main regulator of HGP, HOMA-IR index is practically a measure of hepatic IR. For easy interpretation, lower HOMA-IR values indicate greater insulin sensitivity, whereas higher HOMA-IR values indicate lower insulin sensitivity (i.e., IR).

Another way to see this HOMA index function, another index, insulin sensitivity (IS), is defined as

By applying the same fasting values, pancreatic β-cell function (HOMA β-cell) can be estimated by the evaluation form:

where the unit of FPI is μU/mL, and the unit of FPG is mmol/L.

Moreover, a worth mentioning index, Quantitative Insulin Sensitivity Check Index (QUICKI)[15], is derived by calculating the inverse of the sum of logarithmically expressed values of fasting glucose and fasting insulin:

in which G0is the fasting glucose level and I0is the fasting insulin value. Many investigators and researchers believe that QUICKI is superior to HOMA in determining insulin sensitivity, although those two values correlate well.

Other than the above fasting-value methods, the simplest and widely used test for detection of diabetes is GTT. For this test, a subject fasts for 12 h, and is then given a large amount of glucose. During the next few hours, blood samples are drawn and glucose levels are measured and recorded. By fitting the GTT data to a mathematical model proposed by Ackerman et al.[11,12], the diagnosis information from the model can be applied to indicate which subject has diabetes. This model is described by a differential equation system where the variables are deviations of glucose levels from the subject′s baseline value in blood (in the morning after fasting overnight) and the similar deviation of insulin concentration. The differential equation system governing this GTT model is expressed as follows:where pi(i=1, 2, 3, 4), are positive constants, J is the rate of glucose infusion from the intestines (or intravenously), g(t) is the difference between blood glucose concentration G(t) and its baseline value G0, and i(t) is the difference between plasma insulin concentration I(t) and its baseline value I0, as shown in the following equations.

The reason for this variable transformation is because we are usually more interested in the difference values (i.e., fluctuations or excursions) of glucose and insulin (i.e., the relative values, not the absolute values). The diagram of this two-compartment model is shown in Fig. 1.

Fig. 1 Sketchy diagram of the mathematical model of the glucose-insulin regulatory system.

The construction of the model equations (1.5) and (1.6) is based on the following hypotheses:

1) Each variable, g and i, has various influences upon the appropriate change speed with a negative feedback (i.e., utilization or clearance) process, which is shown as the parameters -p1and -p4in both equations.

2) An increase in blood gluose levels provokes an increase of insulin secretion, which is expressed as the positive feedback (i.e., stimulation) parameter +p3in the second equation.

3) An increase of hormone insulin secretion leads to a reduction in blood glucose levels, which is formulated as the negative feedback (i.e., utilization) parameter-p2in the first equation.

The differential term dg/dt is defined as the change in blood glucose difference with respect to the change in time. Similarly, the expression di/dt is defined as the change in plasma insulin difference with respect to the change in time. To illustrate this method, we have simulated the glucose-insulin dynamics system for both normal and diabetic subjects in the Results section. For the purpose of algorithm integration, the section "Computational methods and theory" will describe the mathematical formulation in greater detail. It has to be noted here that the analytical demonstration of an appropriate qualitative and quantitative behavior for this mathematical model is the great departure point for the subsequent experimentally clinical and theoretically numerical determination of an optimal key parameter for diagnosing the disease.

It has been reported that glycemic variability (or brittleness) scores can be assessed by both the Lability Index (LI)[16,17]and the Mean Amplitude of Glycemic Excursions (MAGE)[18-20]methods.

This LI score provides a measure of blood glucose variability in diabetes and is based on the square of change in levels of blood glucose from one reading to the next, divided by the time interval and summed for a week. A value of LI for each of 4 w is derived based on computing the following sum for each week in the period:

where glucose (mmol/L) is the ith reading of the week taken at time Houri. The upper bound, N, is the total number of readings in 1 w. The minimum and maximum time intervals used are 1 h and 12 h, respectively. Consequently, the median LI in type Ⅰdiabetes control subjects (n = 100) was 223 m(mol/L)2/(h·week) (25thto 75thinterquartile range: 130-329). Patients who received islet transplant (n = 51) had a median LI value of 497 m(mol/L)2/(h·week) (25thto 75thinterquartile range: 330-692) before transplantation. After transplantation, their median LI values became 40 m(mol/L)2h·week (25thto 75thinterquartile range: 14-83). It is clearly seen that there is a large decrease in LI scores after islet transplantation. This indicates that islet transplantation is effective in curing typeⅠdiabetes and also results in better glucose control (i.e., smaller blood glucose swings) for these patients, at least, for the first few years after their transplantations. It is also shown that the LI scoring system does provide a metric to make a quantitative comparison between patient groups and also complement the clinical assessment of glucose variability in diabetic patients.

The MAGE is another measurement and a reflection of how much blood glucose increases or decreases throughout a day. It measures the amplitude of the daily "large" blood glucose excursions. The question is "how large is the amplitude?" The answer is that the intraday glycemic excursions with their amplitudes are greater than one standard deviation (SD). Theoretically, the MAGE requires at least 14 blood glucose measurements over consecutive 48 h before and 2-h after breakfast, lunch and dinner, and at bedtime with an optional measurement at 3 AM. A glycemic excursion is then computed as the absolute difference in peak and subsequent nadir (or vice versa) glucose values, with the direction (decrease - peak to nadir vs increase - nadir to peak) determined by the first quantifiable excursion in the 48 h. All excursions greater than one SD of the 7+ glucose readings for the day were summed and divided by the number of qualified excursions to give the MAGE score in mg/dL (or mmol/L) glucose. To emphasize major glucose swings and eliminate minor ones, however, excursions less than one SD are ignored. It is recognized that the MAGE scores are lower in healthy research participants than in diabetic research participants. In other words, the lower the MAGE score, the less severe the swings of the blood glucose levels. Roughly speaking, MAGE score average is approximately about ＜ 90 mg/dL in healthy subjects and about ＞ 150 mg/dL in diabetic patients or subjects with poor glycemic control. As shown in the study by Ryan et al.[16,17], blood glucose excursions expressed as MAGE score were significantly lower after islet transplantation.

A well-known blood glucose meter, continuous glucose monitoring system (CGMS), is an FDA-approved device that records subjects′ blood glucose levels throughout the day and night. In other words, CGMS is used to provide continuous "real-time" readings about trends in blood glucose levels[21]. This may allow users to know the levels of their glucose and whether they are rising or falling and to intervene by eating food or taking insulin to prevent them from going too high or too low. Clinically, blood glucose regulation can be evaluated based on the CGMS device. DirecNet study group[22]has pointed out that no simple test currently exists for assessment of glycemic variability in patients with diabetes. They have done the first report to compare eight-point testing with CGMS as a means to evaluate glycemic control. Despite the much larger number of measurements with CGMS than with eightpoint testing, the overall mean glucose levels were nearly identical. However, these few-sample-point device and CGMS are the essential tools in measuring the blood glucose levels, and the recorded data are used for physiologic analysis to control the glucose variance.

Recently, the Minimal Model (MM), proposed by the team of Bergman and Cobelli[23-35], is one of the most informative computational methods for studying glucose and insulin kinetics in metabolism. The MM for glucose kinetics is illustrated in Fig. 2.

Fig. 2 The compartmental flow chart of the minimal model for glucose kinetics.

In this figure, I(t) is the plasma insulin level, and Ibrepresents its basal level; G(t) is the plasma glucose level, and its basal level is denoted as Gb. The coupled differential equations corresponding to the glucose minimal model are expressed as

with G(0)=G0and X(0)=0. In these equations, X(t) is the interstitial insulin at time t. There are a total of four unknown parameters in this model: G0, p1, p2, and p3, which are also defined with the units and brief descriptions are provided below:

p1[min-1], glucose effectiveness, p1=SG, the rate of net glucose utilization without dynamic insulin response (i.e., insulin-independent).

p2[min-1], rate constant expressing the spontaneous decrease of tissue glucose uptake ability.

p3[min-2(μU/mL)-1], insulin-dependent increase in tissue glucose uptake ability.

G0[mg/dL], theoretical glycemia at time 0 after the instantaneous glucose bolus.

The MM describes the time-course of glucose plasma concentrations, depending upon insulin concentrations and the new intermediate variable X, representing the "insulin activity in a remote compartment". This synthetically contrived and physiologically inaccessible variable X plays a transitive role between blood glucose and plasma insulin. It is clearly shown that this virtual variable in equation (1.11) has replaced the position of the insulin variable in equation (1.6).

It is commonly recognized that MM is a well-known and successful mathematical model in simulating glucose metabolism and insulin kinetics. This is also why we cannot skip the introduction of the model. In our current study, however, we focused mainly on presenting a new approach of parameters estimation on glucose-insulin regulatory system. Since MM has been extensively studied by many clinicians and researchers, we have restricted ourselves to only the coupled ordinary differential equations based on the Ackerman′s model. Unquestionably, the application of MM using our proposed approach will be carried out in future trials.

COMPUTATIONAL METHODS AND THEORY

Glucose-insulin dynamics

The majority of mathematical models proposed in the literature were dedicated to the dynamics of glucose-insulin, including oral glucose tolerance test (OGTT), intravenous glucose tolerance test (IVGTT), and frequently sampled intravenous glucose tolerance test (FSIVGTT). In 1970, various mathematical models have been applied to evaluate glucose disappearance and insulin-glucose dynamics[36-38]. Particularly, Bolie[9], among the pioneers in this field, used two coupled ordinary differential equations to constitute the following simple model:

where G = G(t) represents blood glucose concentration, I = I(t) represents plasma insulin and p, a1, a2, a3, and a4are parameters.

In 2005, Wu[39]proposed the following model for self-management of type Ⅱ diabetes:

where x represents blood glucose level over the baseline at time t and ω0is the system natural frequency. The postprandial blood glucose excursion can be considered as a hormone regulated springy system. The food intake is treated as a bolus injection of glucose and expressed as the impulse force f(t). Effects of exercises and hypoglycemic medication are lumped as the damping factor, β. The differential equation of such an oscillatory system, which is used to describe postprandial blood glucose excursions, can be found in many physics and mathematics textbooks.

Generally speaking, preprandial blood glucose levels are generally fluctuating with relatively insignificant magnitudes, and, thus, can be approximated as a flat level. If the impulse force f(t) takes the form of the Dirac delta function, Fδ(t-0) with F being a food intake dependent parameter, the solution of equation (2.3) is given by

where

is the frequency of the system. As you can see from equation (2.4), it is described as a three-parameter model, F, ω and β. Implications of these three parameters not only could reveal distinctive characteristics of diabetic and normal individuals but also provide guidelines to adjust diabetic patient′s treatment and lifestyle.

It has been demonstrated in the study of Ackerman et al.[40]that the response of blood-glucose concentration (G) as a function of time (t) is represented adequately by an equation involving four parameters in the equation:

The values of these four parameters are defined by the four experimental measurements (see Table 1) usually made in an ordinary GTT. On the basis of measurements in over 750 persons in their study, it was suggested that the value of ω0could be used to distinguish normal from diabetic subjects more closely than any other parameters. It should also be noted that the solution of the oscillatory pattern in equation (2.6) has the same form as the expression in equation (2.4).

As indicated earlier, the mathematical models proposed by Ackerman[11,12,40], Bolie[9]and Wu[39]are the simplest models and have been indisputably useful in physiological diabetes research and served as a starting point for many other models. Therefore, our present work extends the scope of the public-recognized model by introducing a novel approach using the perturbation method for parameters estimation and the multiple-shooting method for solving two coupled ordinary differential equations on glucose-insulin regulatory system.

The derived model equations in this study are written as the forms in equations (1.5)-(1.8), which describes the blood glucose regulatory system during a GTT. Actually, the starting equations are expressed as these simple forms:

where t is time, G is blood glucose concentration, and I is plasma insulin concentration. The impulse function, J(t), is basically a delta function representing the large dose of glucose given initially in the GTT after the subject has fasted. In this section, we only briefly list some key formulas and equations for the purpose of attracting readers′ focus.

After a lengthy derivation, the deviation g(t) of a subject′s blood glucose concentration satisfies the second-order differential equation:

which can be rearranged as

In physics terminology, if a frictional (i.e., damping -α term) force proportional to the velocity does exist, the harmonic oscillator will be described as a "damped" harmonic oscillator. When the damping factor, α, is greater than zero, the dynamic system will decay, but may or may not oscillate, depending on the relation between the damping factor α and the natural frequency ω0. As a consequence, the solution g(t) will have three possible outcomes:

However, all those curve tendencies may not be reflected as the same patterns on other damped oscillatory systems, as they really depend upon their amplitudes, angular frequencies and damping factors. Consequently, the period of oscillation is finally formulated as

According to the general diabetologist′s theory (detailed results and references in section Results), healthy individuals should satisfy the condition of T=2π/ω0＜ 4 h, conversely, while the reversed inequality of the period T indicated diabetic cases. As a matter of fact, this conditional situation actually occurs frequently in practice - the blood glucose concentration of a non-diabetic who has just absorbed a large amount of glucose will be at or below fasting level in 2 h or less.

Nonlinear Least Squares Fitting and Gauss Newton Method

The linear regression model is a system of linear equations that can express the model using data matrix X, response vector Y and parameter vector P. The ithrow of X and Y will contain the x and y value for the ithdata element. Thus, the model can be written as a detailed matrix form:

where ε is the normally distributed random error term with the expected value of a column vector of zeros and variance σ2In, where Inis the n×n identity matrix.

The pure matrix notation becomes

where

Y is an n-by-1 vector of responses X is the n-by-m design matrix for the model, where m = 2

P is a m-by-1 vector of parameters, α and β, where m = 2

ε is an n-by-1 vector of random errors

For nonlinear least squares fitting to a number of unknown equations′ parameters, linear least squares fitting could be applied iteratively to a linear form of the function until convergence is achieved. Now we consider the nonlinear models of the form

The least squares method still applies, with the sum of squared residuals (SSR) function to be minimized with respect to the parameter vector P. Please note that SSR is sometimes called RSS as it stands for residualsquared-sum. The least square estimate of P,, is the set of parameters that minimizes the sum of squared residuals:

To obtain the normal equations, we consider the partial derivatives of SSR() with respect to each parameter, and set them equal to zero. This gives us a system of m normal equations. Each normal equation is given by differentiating SSR(),Unfortunately, the partial derivatives of nonlinear functions are also functions of the parameters; thus, an explicit solution (i.e., analytical solution) forcannot be obtained. We will have to get the solution by numerical optimization.

The basic concept of the Gauss-Newton optimization method can be summarized as:

1) Use a Taylor series expansion to linearize a nonlinear function;

2) Apply least-square theory to obtain new estimate of the parameters that move in the direction of minimizing the SSR.

1) Start with initial estimated values for each adjustable parameter in the equations.

2) Generate the data trajectory for the entire curve based on the initial input values. Independent random number generator and numerical integration subroutine are required for implementation.

3) Calculate the sum-of-squares of the vertical distances of the data points from the curve. Or compute the weighted sum-of-squares if weighting factors are included.

4) Adjust the parameters to make the curve come closer to the data points by evaluating the corresponding cost (or objective) function.

5) Stop the calculations and go to step 6 when the adjustments make virtually no difference in the sumof-squares or any available stopping criteria are satisfied; if not, repeat simulation steps 2-4 with newly selected initial values.

6) Report the best-fitted parameter values.

Please note that steps 2, 3 and 4 are the only ones whose details differ from method to method.

The Perturbation Search Method

We should remark here that gradient-based methods frequently encounter divergent difficulties, and convergent case only occurs under certain restricted conditions. In addition, gradient-based methods usually pose the risk of being trapped in a local minimum, rather than in the global minimum. If inappropriate values of the initial parameter are selected and run for the entire simulation because they yield the best quality of data fit for the predefined working range, the incorrect conclusion will be drawn quickly without letting developers having any further tests using other initial parameter sets. In common practice, local methods optimize the cost function directly with respect to initial values and parameter vector. This optimization scheme is often called initial-value approach or alternatively named "single-shooting method". Due to the presence of multiple minima in some cases, convergence of local optimization methods to the global minimum is usually limited to a rather small domain in search space. In contrast, global methods have an essentially larger convergence domain but the computational cost increases rapidly. At this point, it has to be mentioned that one of the stable and effective global searching methods is the multiple-shooting algorithm developed by Bock[41-44]; a brief description will be given in the next subsection. The basic idea of the algorithm is to consider the task of the single initial-value approach as a multi-point boundary value problem. Mostly importantly, this multiple-shooting algorithm is used to incorporate the proposed perturbation-search method in our current study. It has been concluded that single-shooting method steps quickly in local minima or runs into a divergent situation, whereas multiple-shooting method performs better than single-shooting approach at the expense of higher computational costs.

As a matter of fact, our concept of the perturbation search method was inspired by the prospect of basic perturbation theory. The principle of perturbation theory is to study dynamical systems that are small perturbations of integrable systems. Perturbation theory is a technique in quantum mechanics for finding approximate solutions to real-world problems. It has a wide variety of applications throughout chemistry, physics and mathematics. In quantum mechanics, we often need to solve the Schrodinger equation that corresponds to a given physical potential or Hamiltonian. For a real-world potential or Hamiltonian, finding a solution analytically is almost impossible. In perturbation theory, a complicated potential is broken down into a solvable potential (without perturbed terms), plus a small-change remainder. Basically, we solve the first part (without perturbed terms) of the equation as the reference solution, and then consider how the second part (with perturbed terms) influences our answer. In general, perturbation theory solves such problem in two steps. First, we obtain the eigenfunctions and eigenvalues of the unperturbed Hamiltonian. Second, these eigenfunctions and eigenvalues are corrected to account for the perturbation′s influence. As a matter of fact, perturbation theory gives these corrections as an infinite series of terms, which become increasingly smaller for well-behaved dynamic systems. Quite frequently, the corrections are only taken through first or second order. As long as the perturbation is small compared to the unperturbed potential or Hamiltonian, perturbation theory will tell us how to correct the solutions to the unperturbed problem and to approximately account for the influence of the perturbation. It has to be indicated that this theory has no direct relationship with our mathematical model, but it does serve as a stimulant to creative thought of our ODE′s parameters estimation approach. Recently, we have presented a perturbation-based estimate algorithm for parameters of coupled ordinary differential equations and its application from chemical reactions to metabolic dynamics[45].

Perturbation search is a systematic direct search method for solving optimization problems that does not require any information about the gradient of the cost function. However, it does require the presence of pre-defined perturbed values for each parameter pi. Briefly speaking, the algorithm searches a set of symmetric points around the current parameter values, looking for the points where the cost function is lower than the cost function at the previous parameter set. This systematic direct search method can be used to solve a variety of parameters′ estimation problems that are not well suited for standard optimization algorithms, including problems in which the cost function is discontinuous, nondifferentiable, stochastic or highly nonlinear.

While the gradient of cost functions is evaluated, the derivatives of the cost function can be intuitively known and simply formulated as

Thus, it is computed numerically through the finite difference method. The choice of ?can be rather influential, and generallyshould be taken as some fractional part of pj, such as 0.01 or 0.1.

Speaking of the gradient function or the "derivative" term in mathematics, we should also mention that, although the term of the gradient of cost function, such as the expression in the above equation, is not meaningful, it has already conveyed some sense of the "directional derivative". In light of equation (31), it is also worthwhile to note that instinctive understanding of the principles of metabolic dynamics that govern how the trajectories or cost functions move in various situations is vital to success in parameters searching. As everyone knows, in calculus, the derivative is a measurement of how a function changes when the values of its inputs change. Loosely speaking, a derivative can be thought of as how much a quantity is changing at some given point. To interpret this term, a directional derivative, a partial derivative of f that measures its variation in the direction of the coordinate axes, should be introduced here, which enables us to find the rate of change of a function of two or more variables in any direction. For example, if f is a function of x and y, then its partial derivatives measure the variation in f in the x direction and the y direction. If we use three variables, we can define the directional derivatives in a similar manner. For instance, the directional derivative of f at (x0, y0, z0) in the direction of a unit vectoris

Now, suppose we have a function f of three variables, x, y, and z, and we consider all possible directional derivatives of f at a given point. These directional derivatives will provide us the rates of changes of f in all possible directions. The question is in which of these directions the function f changes the fastest and what the maximum rate of change is. According to equation (4.3), we have the answer: the maximum value of the directional derivativeisand it only occurs whenis pointing to the same direction as the gradient vector. A comparison between equation (4.1) and its similar equation (4.2) does enlighten us to pursue our perturbation search algorithm, even though the variables used in these two equations are different - one describes the variation of ODE parameter, and the other one represents the displacement on axis.

In numerical analysis, finite difference method does play a significant role. They are one of the simplest ways of approximating a differential operator and are extensively used in solving differential equations. The difference quotient is defined as the increment of the value of a function divided by the increment of the independent variable, which includes three common forms: forward, backward and central differences. A forward difference of function is an expression of the form f(x+h)-f(x). A backward difference of function uses the function values at x and x-h, instead of the values at x+h and x, which is f(x)-f(x-h). The central difference of function is given by f(x+h/2)-f(x-h/2). All those three function′s increments divided by the independent variable′s increment, h, will give three different approximated derivative values. The centered difference approach is more accurate than both forward and backward difference formulas. Conceptually, the central difference approach uses two points on both sides of x, so it can reach the 2ndorder accuracy if function f is expended in a Taylor series around the point x. It has to be noted here that this concept does give rise to motive ideas and becomes the source of inspiration for us to extend the perturbation theory idea to the estimation of ODE′s parameters.

In fact, numerical optimization with the perturbation search method in the current study is a "trial and error" process. We begin with a single parameter as the starting example. Starting from an initial set of parameter values for which the initial cost function is known, a new set of parameter values is calculated by perturbing each of the initial parameter values, and the corresponding cost functions will be evaluated. More specifically, three trials, two perturbation trials plus one trial without perturbation, for each iteration run are simulated by choosing the algebraic sign (-), reversing of sign (+) and no sign of each base parameter values. The two computed cost function values SSR(p -△p) and SSR(p+△p), corresponding to negative and positive perturbations, plus the unperturbed cost function SSR(p) value are compared to the cost function value at the initial point or previous point. If a cost function value is less than the previous minimum cost function value, SSRmin, the current set of parameter values for which the cost function has its least value becomes the new starting point for the next trial. If neither value is less than the previous SSRminvalue, the previous parameter values and its previous SSRminwill remain the same for the next trial.

The multiplication rule can be used to determine the total number of outcomes in a sequence of events. In a sequence of n events, the first one has k1possibilities and the second event has k2and the third has k3, and so forth. Thus, the total number of possibilities of the sequence will be the product of all possibility values,

k1·k2·k3…….

It can be concluded that the 2-possibility n events produce a total of 2npossibilities. Similarly, if there is a sequence of n = 5 events which have three possibilities for each, it will yield the total possibility number 3n= 35= 243. The programming logic of this iterative parameter combination scheme is illustrated as the pseudo-code of 5-parameter-loop in Fig. 3.

Fig. 3 The algorithm for looping all five parameters.

It is worthwhile to mention that the parameter optimization can also be performed by a systematic reducing parameter search range. For instance, the relatively-wide initial search region is selected and used for running a number of iterations until the cost function reaches a constant-like value. Then, a reduction in the size of the search region at subsequent iteration will be carried out. The next iterations will follow the same procedures except that the search ranges of parameters are smaller. Basically, we hope the accuracy of the parameters will be enhanced through this range-reducing approach. From the practical point of view, the range reduction factor can be proposed as any appropriate expression based on the results of a trial run and output of each iteration. These steps are repeated until almost no variation of the estimated parameters or any predefined convergent criteria of the parameters and cost function are met. However, this adjusted approach of parameter range reduction can be implemented differently to suit the properties each individual ODE system.

Multiple-Shooting Method

In the initial-value approach of the ODE system, initial guesses for equation parameters are chosen and the dynamical equations are solved numerically. The parameters can be identified as those minimizing the cost function. For many dynamical simulations, this approach is numerically unstable by yielding a diverging trajectory or stopping in a local minimum. This problem can be circumvented by a multipleshooting algorithm developed by Bock[41-44]. The basic idea of the algorithm is to consider the task as a multipoint boundary value problem. The fitting time interval [ti, tf] is partitioned into k subintervals,

For each subinterval, local initial values are introduced as additional parameters. The dynamical equation is integrated piecewise and cost function is evaluated and minimized as exact in the initial value approach. While the dynamical parameter vector P is constant over the entire interval, the local initial values are optimized separately in each subinterval. Starting guesses for them are appropriately chosen to match the observations. This approach leads to an initially discontinuous trajectory, but is close to the measurements. The final trajectory will eventually be continuous as the computed solution at the end of one subinterval matches the local initial values of the next one. This condition is taken into account as equality constraints in the optimization procedure.

To achieve a good estimating of the domains of ODE parameters, we propose an optimal search with subinterval constraints and developed an algorithm to enforce it. The searched time interval is divided into a number of subintervals, which is the number of data points minus one. We have extended Bock′s concept by replacing the one single-shooting SSR value in the initial-value approach by a summed value of "multiple" subintervals′ SSR values. To illustrate the application of this method, a conceptual multiple-shooting method diagram is shown in Fig. 4. As you can see from this figure, we sum over all the SSR components, which are represented as green lines with double arrows, of each subinterval.

Fig. 4 Sketchy diagram illustrates the application of multiple shooting method.

Coupled Ordinary Differential Equations Solver

To illustrate the solving method of ODE system, a system of first-order ODE is expressed as

where n = 1, 2, 3, ……, N

In the dynamics of time series, the independent variable x is usually replaced with time t. Therefore, a simple but general vector form of a system ODEs with the initial condition is formulated as

where y is the state vector of the system, t is time and θ = [p, y0] is a parameter vector consisting of all p parameters appearing in the right hand side of the ODE system and all initial-value vector y0.

Among numerical integration methods, the Euler method is the simplest integration method for solving ordinary differential equations, which is also called first-order Runge-Kutta method or Taylor series expansion of order one. The next common integration method is the second-order Runge-Kutta method, and its truncation error is higher than the one in the Euler method. However, detailed formulations and their applications of these methods are easily found in the literature and many numerical method books, which can also be clearly referenced in References[45,46].

Program flowchart

The flowchart of our numerical optimization algorithm for ODE′s parameter estimations is presented in Fig. 5. This flowchart outlines the order in which these functional elements in our computational approach were sequentially performed.

Program requirements and specifications

The implementation of the computerized perturbation search methodology is straightforward. However, two functions or subroutines are required to accomplish the entire simulation. The first program is a coupled differential equation solver using any numerical integration method, such as the Runge-Kutta method. The second program is a random number generator. Fortunately, these two programs are available and easily found in most programming languages and scientific software packages. It should be mentioned that R is a free software environment for statistical computing and graphics, which provides a wide variety of statistical linear and nonlinear modeling. Therefore, our code of simulation model was written in R, but other languages, such as Fortran, C/C++, Java and Visual Basic, or other commercial mathematical and statistical software packages, such as Matlab®and SAS®, could also be used. In addition, the hardware platform is IBM-compatible personal computer (PC) with Microsoft®Windows 2000 and XP operating systems.

RESULTS

Here we give an example with four sets of data points to depict the glucose-insulin regulatory dynamic trajectories[47]. Table 1 shows an example for two subjects whose blood glucose and plasma insulin deviation values were sampled at four time points, initial-0 h, 1 h, 2 h, 3 h.

Based on those data, we estimated the four parameters, p1, p2, p3, and p4, by using our perturbation search algorithm and multiple-shooting method. In addition, those computed parameter values are also applied to make fitted curves compared to the original sample data points (see the circles in Fig.6). As you can see from this figure, the theoretically simulated curves perfectly match the experimental data points.

As you can see from the graphs and the computed results, the four parameter values, which vary from person to person, determine the evolution of the trajectory over time. Both subjects started with the positive glucose variation 80 mg/dL (i.e., 80 mg/dL above the fasting baseline), which were immediately measured right after the ingestion of a large dose of glucose. As you can tell from those glucose curves, there is a rapid drop at 1 h time for case 1; in contrast, there is a relatively slower decrease for case 2. As a whole, the glucose trajectory of case 1 oscillates faster than that of the case 2, and it also quickly reaches the equilibrium state (i.e., fasting level). For the same reason, a similar faster variation of the plasma insulin concentration for case 1 produces quicker fluctuations (i.e., higher frequency). This reflects that case 1 has a shorter period of time on recovering from an initially large glucose excursion to an equilibrium fasting level. This result thus indicates that case 1 is relatively normal (i.e., non-diabetic) compared to case 2.

The R-based code fragments and outputs are listed below for your reference.

Fig. 5 The flowchart of the current simulation program.

Our computed oscillation periods are Tsubject1=179.88 min and Tsubject2=289.51 min. In Jorge Cruz′s paper[47], their calculated oscillation periods are Tsubject1=179.9 min and Tsubject2=289.9 min. These extremely close (i.e., consistent) results imply that the first subject is normal and the second one is diabetic. This also indicates that our computation method is correctly validated.

As you can see from these R-based programming calculated results, the values of alpha (α) are less than the values of omega zero (ω0), so the oscillations types we have demonstrated are categorized as the underdamping type, α2＜ ω0

Fig. 6 Simulations of the evolution of blood glucose and insulin concentration.

Table 1 Four measurements (deviations from baselines) of blood glucose and plasma insulin concentration in two subjects.

2. However, the detailed damping "type" classification is referenced in our study note.

Based on the final solutions of the three-type glucose variation, those curves are plotted and compared. The R-based code fragments and graphical outputs are shown in the following two diagrams, Fig.7 and Fig.8.

Fig. 7 Examples of damped oscillatory motion with arbitrary constants ω=0.001 and α=0.005.

The choices of two coefficients of the solutions, c1=10 and c2=0.1, are arbitrarily selected, but there is only one condition that needs to be met in this simulation - all the simulated data should be started at the almost same initial point, such as 10 in these plots. Thus, based on the same starting point, we can visualize and see the difference on how those metabolic trajectories move.

As you can see from Fig. 8, the over-damping oscillation represented by a green curve is dramatically elevated to an unrealistic high level. In practice, this situation in plasma glucose response to GTT actually can never occur. On the other hand, we should suggest that the curve patterns in both the under-damping and critical-damping types could play a reasonable role on the platform of glucose metabolic simulations.

Fig. 8 Examples of damped oscillatory motion with arbitrary constants ω=0.01 and α=0.005.

DISCUSSION

Monitoring blood glucose concentration is crucial to the successful diagnosis and treatment of people with diabetes. The earlier the diagnosis of the disease the better the chances of controlling it with insulin and helping the subject live longer. One simple test for diagnosing diabetes is GTT, where the subject ingests a large amount of glucose, and then has his or her blood monitored for about 4-6 h following the glucose administration. The computational model in the present work has investigated this diabetic GTT metabolic test one step further, with respect to the already published Ackerman[11,12,40], Bolie[9], Wu[39], and Cruz[47]models, although they have provided a widely-accepted tool for clinicians and researchers to evaluate the excursions of blood glucose levels and to diagnose the existence of diabetes.

In the section of Results, we present a motivating example from the online published articles[47]regarding the theoretical diagnosis of diabetes. The classification of the dynamical system′s behaviors, such as oscillating with period of time, T, usually depends on some conditions or relations between these system′s parameters, such as damping factor and angular frequency.

The most notable motion from our simulations and graphs is the under-damping oscillation, which is described by the product of the exponential term and sinusoidal cosine function. This product solution gradually decreases the amplitude of the oscillation until it reaches the equilibrium state (i.e., the term of "fasting state" or "baseline" in diabetes field). The exponential term acts as a boundary constraint for the amplitude of the sinusoidal function, resulting in a gradual decrease of oscillation. Another important concept from the graph is that the period of the oscillation does not change, even though the amplitude is constantly decreasing. Therefore, the period, T=2π/ω0, is the good index to be computed in this investigation and other studies. Other than that reason, based on other researchers′ previous experience and evidence, experimental testing of the model showed that the parameter,, was quite robust and proved to be a good indicator of diabetes. In particular, healthy individuals satisfied the criteria:

To simply interpret this mathematical expression, the blood glucose concentration of a non-diabetic who has just absorbed a large amount of glucose will be at or below fasting level in 4 h or less. On the contrary, the reversed inequality indicates diabetic case.

In practice, for a particular subject, the accuracy of the diagnosis can be benefited from multiple (i.e., repeated) measurements as the reliable ranges of the parameters, the values of ′mean±standard deviation′ or even confidence intervals, can be accurately estimated. Thus, clinicians can safely evaluate the diagnostic conditions based on the reliably estimated parameters and the above formulas′ criteria. However, due to the constraints of limited resources and the uneasy-fitted scheduled procedure time, GTT is usually performed only once within several months or a year. Nevertheless, our non-repeated GTT calculated result suggests that the dynamic response of diabetic subjects to typical hourly glycemic challenges extends over a much broader range in time domain than that of non-diabetic individuals. In other words, the oscillation period is higher in diabetic subject and lower in normal subject.

Ideal studies for characterizing blood glucose dynamics should include blood glucose measurements at a sufficiently regular sampling frequency to detect all variations over a period of several days (i.e., in MAGE evaluations), or even several weeks long (i.e., in LI measurements), while the subject is undergoing a wide range of activities, such as meals, exercise and sleep. However, our simulations on the GTT studies were based on an ideal case with only four sample points at the hourly time ranges, (0-3 h), which are not sufficient enough to fit them to an appropriate curve that describes the entire blood glucose dynamics as well as to examine the special features of glucose excursions. Theoretically speaking, there are at least eight well-recorded sample points to monitor patient′s glucose regulatory phenomena.

Lastly, based on the above two equations as well as our calculated results in section Results clearly indicate that the Ackerman model and also the extended approach in the current study may play an important role in diagnosing diabetes, or even provide clinicians and scientists in endocrinology and metabolism fields insight into the transition nature of human′s metabolic mechanism from normal to impaired glucose tolerance.

In conclusion, from a mathematical point of view, the conventional parameters estimation methods, such as the Gauss-Newton method, usually have the drawback of singular inverted matrix. The main purpose of our approach is to develop a simple, stable and formula-derivation-free computational program and, most importantly, to avoid divergence problem. Furthermore, there are no bounded constraints applied, and thus the movements of all equations′ parameters during the entire computer simulation are absolutely free. Clinically speaking, in order to determine whether or not a subject has pre-diabetes or diabetes, health care providers usually conduct a GTT and apply standard HOMA and QUICKI methods. The aim of this study is to explore how to interpret those laboratory data as well as enhance the Ackerman mathematical model and how numerical analysis and computational iteration program are developed to search and identify coupled ODEs′ parameter values.As presented in four recently published articles[45,48-50]the proposed mathematical models and computation methods were successfully applied to have effectively direct measurements of insulin hormones secretion and dynamics in glucose metabolism. The present work addressed the application of multiple-shooting method and perturbation-search algorithm to biomedical problems, particularly to endocrinology and metabolism fields, which involves two coupled differential equations with four parameters describing the glucose-insulin regulatory system. Our demonstration showed that our approach could be used to practically evaluate the subject′s glucose deviations and diagnose the diabetic disease as well as to theoretically serve as a promising starting point on ODE′s parameters estimation. Hopefully, the current computer simulation can be studied further to provide investigators an effective tool and key information for the design of clinical studies that involve blood glucose dynamics and also for the possible development of new glucose monitoring systems. On the other hand, the extended and enhanced approach in the current study could play an important role in diagnosing diabetes, or even provide clinicians and scientists in endocrinology and metabolism fields insight into the transition nature of human metabolic mechanism from normal to impaired glucose tolerance. Furthermore, the mathematical models can be used to suggest better treatments, such as to improve scheduling of insulin injections, to regulate patient′s diet, and have a solid pancreas or novel islet-cell transplant treatment. Due to both new experimental clinical studies and improved mathematical models, future treatments will certainly evolve to better regulate this metabolic disease′s problem in diabetic patients. Lastly, it should be emphasized that our integrated approach of perturbation search and multiple-shooting methods provides an attractive feature to numerical estimation of the coupled differential equations, although parameters estimation has been recognized as a well-known challenge and accomplishment in mathematics and statistics.

ACKNOWLEDGEMENTS

The authors wish to thank Professors Jeffrey Longmate and Ravi Bhatia at City of Hope for many fruitful and encouraging discussions on the topic and generous financial supports. The authors are also greatly indebted to the anonymous referees for their valuable comments of an earlier draft, which have resulted in an improved presentation of this paper.

[1] American Diabetes Association (ADA), All about diabetes, Online ～available from the World Wide Web: http://www.diabetes.org/home.jsp.

[2] Centers for Disease Control and Prevention (CDC), Diabetes Public Health Resource, Basics about diabetes, Online ～available from the World Wide Web: http://www.cdc.gov/diabetes/faq/basics.htm.

[3] Boutayeb A, Chetouani A. A critical review of mathematical models and data used in diabetology. Biomed Eng Online 2006;5:43-51.

[4] Mari A. Mathematical modelling in glucose metabolism and insulin secretion, Curr Opin Clin Nutr Metab Care 2002;5:495-501.

[5] Bergman RN, Finegood DT, Ader M. Assessment of insulin sensitivity in vivo. Endiocrine Reviews 1985;6:45-86.

[6] Bergman RN. The minimal model: yesterday, today and tomorrow. In: Minimal model approach and determinants of glucose tolerance. Lovejoy, Boston, Louisiana State University Press 1997;7:3-50.

[7] Bergman RN. The minimal model of glucose regulation: a biography. In: Mathematical modeling in nutrition and health. Kluwer Academic/Plenum, 2001.

[8] Pattaranit R, van den Berg HA. Mathematical models of energy homeostasis. J Royal Society Interface 2008;5:1119-35.

[9] Bolie VW. Coefficients of normal blood glucose regulation. J Appl Physiol 1961;16:783-8.

[10] Simon C, Brandenberger G. Ultradian oscillations of insulin secretion in humans. Diabetes 2002;51(Suppl 1) 258-61.

[11] Ackerman E, Gatewood LC, Rosevear JW, Molnar GD. Model studies of blood glucose regulation. Bull Math Biophys 1965; 27:21-37.

[12] Ackerman E, Gatewood LC, Rosevear JW, Molnar GD. Blood glucose regulation and diabetes. In: Concepts and Models of Biomathematics. Marcel Dekker 1969;131-56.[13] Matthews DR, Hosker JP, Rudenski AS, Naylor BA, Treacher DF, Turner RC. Homeostasis model assessment: insulin resistance and beta-cell function from fasting glucose and insulin concentrations in man. Diabetologia 1985;28:412-9.

[14] Wallace TM, Levy JC, Matthews DR. Use and abuse of HOMA modeling. Diabetes Care 2004;27:1487-95.

[15] Katz A, Nambi SS, Mather K, Baron AD, Follmann DA, Sullivan G, Quon MJ. Quantitative insulin sensitivity check index: a simple, accurate method for assessing insulin sensitivity in humans. J Clin Endocrinol Metab 2000;85:2402-10.

[16] Ryan EA, Lakey JRT, Rajotte RV, Korbutt GS, Kin T, Imes S, et al. Clinical outcomes and insulin secretion after islet transplantation with the Edmonton protocol. Diabetes 2001;50:710-9.

[17] Ryan EA, Shandro T, Green K, Paty BW, Senior PA, Bigam D, et al. Assessment of the severity of hypoglycemia and glycemic lability in type 1 diabetic subjects undergoing islet transplantation. Diabetes 2004;53:955-62.

[18] Service FJ, Molnar GD, Rosevear JW, Ackerman E, Gatewood LC, Taylor WF. Mean amplitude of glycemic excursions, a measure of diabetic instability. Diabetes 1970;19:644-55.

[19] Service FJ, O′Brien PC, Rizza RA. Measurements of glucose control, Diabetes Care 1987;10:225-37.

[20] Selam JL. How to measure glycaemic instability. Diabete Metab 2000;26:148-51.

[21] Klonoff DC. Continuous glucose monitoring. Diabetes Care 2005;28:1231-9.

[22] Diabetes Research in Children Network (DirecNet) Study Group. Eight-point glucose testing versus the continuous glucose monitoring system in evaluation of glycemic control in type 1 diabetes. J Clin Endocrinol Metab 2005;90:3387-91.

[23] Bergman RN, Ider YZ, Bowden CR, Cobelli C. Quantitative estimation of insulin sensitivity. Am J Physiol 1979;236: E667-77.

[24] Toffolo G, Bergman RN, Finegood DT, Bowden CR, Cobelli C. Quantitative estimation of beta cell sensitivity to glucose in the intact organism: a minimal model of insulin kinetics in the dog. Diabetes 1980;29:979-90.

[25] Pacini G, Bergman RN. MINMOD: a computer program to calculate insulin sensitivity and pancreatic responsivity from the frequently sampled intravenous glucose tolerance test. Comput Meth Prog Biomed 1986;23:113-22.

[26] Steil GM, Volund A, Kahn SE, Bergman RN. Reduced sample number for calculation of insulin sensitivity and glucose effectiveness from the Minimal Model, Diabetes 1993;42:250-6.

[27] Caumo A, Cobelli C. Hepatic glucose production during the labeled IVGTT: estimation by deconvolution with a new minimal model. Am J Physiol 1993;264:E829-41.

[28] Caumo A, Cobelli C, Finegood DT. Minimal model estimate of glucose effectiveness: role of the minimal model volume and of the second hidden compartment. Am J Physiol 1998;274:E573-6.

[29] Caumo A, Giacca A, Morgese M, Pozza G, Micossi P, Cobelli C. Minimal model of glucose disappearance: lessons from the labelled IVGTT. Diab Med 1991;8:822-32.

[30] Caumo A, Vicini P, Cobelli C. Is the minimal model too minimal? Diabetologia 1996;39:97-1000.

[31] Cobelli C, Pacini G, Toffolo G, Sacca L. Estimation of insulin sensitivity and glucose clearance from minimal model: new insights from labeled IVGTT. Am J Physiol 1986;250:E591-8.

[32] Cobelli C, Thomaseth K. The minimal model of glucose disappearance: optimal input studies. Math Biosci 1987;83:125-55.

[33] Cobelli C, Toffolo G, Ferrannini E. A model of glucose kinetics and their control by insulin, compartmental and noncompartmental approaches. Math Biosci 1984;72:291-315.

[34] Cobelli C, Pacini G. Insulin secretion and hepatic extraction in humans by Minimal Modeling of c-peptide and insulin kinetics. Diabetes 1988;37:223-31.

[35] Nucci G, Cobelli C. Models of subcatuneous insulin kinetics: a critical review. Comput Meth Prog Biomed. 2000; 62:249-57.

[36] Della C, Romano MR, Voehhelin MR, Seriam E. On a mathematical model for the analysis of the glucose tolerance curve. Diabetes 1970;19:145-8.

[39] Wu HI. A case study of type 2 diabetes self-management. Biomed Eng Online 2005;4:4.

[37] Serge G, Turcogl M, Varcellone G. Modelling blood glucose and insulin kinetics in normal diabetic and obese subjects. Diabetes 1973;22:94-7.

[38] Srinivasan R, Kadish AH, Sridhar R. A mathematical model for the control mechanism of free-fatty acid and glucose metabolism in normal humans. Comp Biomed Res 1970;3:146-9.

[40] Ackerman E, Rosevear JW, McGuckin WF.A mathematical model of the glucose tolerance test. Phys Med Biol 1964;9:202-13.

[41] Bock HG. Numerical treatment of inverse problems in chemical reaction kinetics, in: K. Ebert, P. Deuhard, W. Jager (Eds.), Modelling of Chemical Reaction Systems, Springer, Heidelberg, NY, 1981, pp. 102-125.

[42] Bock HG. Recent advances in parameter identification techniques for ordinary differential equations, in: P. Deuflhard, E. Hairer (Eds.), Numerical Treatment of Inverse Problems in Differential and Integral Equations, Birkhauser, 1983, pp. 95-121.

[43] Bock HG. Randwertproblemmethoden zur parameteridentifizierung, in: systemen nichtlinearer differentialgleichungen, PhD Thesis, Universita¨t Bonn, 1987.

[44] Bock H, Plitt K. A multiple shooting algorithm for direct solution of optimal control problems, in: Proceedings of the 9thIFAC World Congress Budapest, Pergamon, Oxford, 1984.

[46] Press WH, Teukolsky SA, Vetterling WT, Flannery BP. Numerical recipes in C, The Art of Scientific Computing, Second Edition, Cambridge University Press, 1992.

[47] Cruz J, Barahona P. An interval constraint approach to handle parametric ordinary differential equations for decision support, Online ～ available from the World Wide Web: http://centria.di.fct.unl.pt/～pb/papers/cruz_ekdb99.pdf and http://centria.di.fct.unl.pt/～pb/papers/cruz_cp99.pdf

[45] Shiang KD. A perturbation-based estimate algorithm for parameters of coupled ordinary differential equations, applications from chemical reactions to metabolic dynamics. Comput Meth Prog Biomed 2009;94:118-42.

[48] Trujillo-Arriaga HM, Román-Ramos R. Fitting and evaluating the glucose curve during a quasi continuous sampled oral glucose tolerance test. Comput Biol Med 2008;38:185-95.

[49] Fernandez M, Villasana M, Streja D. Glucose dynamics in Type I diabetes: Insights from the classic and linear minimal models. Comput Biol Med 2007;37:611-27.

[50] Morbiducci U, Tura A, Grigioni M. Genetic algorithms for parameter estimation in mathematical modeling of glucose metabolism. Comput Biol Med 2005;35;862-74.