The model

Glucose, variable x, enters the blood circulation in two ways: from digested food containing carbohydrates and from the liver where it was stored as glycogen (‘glycogenolysis’). Glucose may also be newly synthesized from precursors like certain amino acids and glycerol (‘gluconeogenesis’). During an overnight fast low insulin levels favour hepatic glycogenolysis and hepatic and renal gluconeogenesis [17]. We neglect renal glucose production and consider only the glucose inflow from digesting a meal and glucose that is released from liver, i.e. hepatic glucose production which combines both glycogenolysis and gluconeogenesis. This last flow is assumed to be proportional with the concentration of glucagon in the blood, with proportionality constant ε. This condition may hold if the glycogen level in liver is sufficiently high. Glucose leaves the blood circulation to be used for energy expenditure in cells and by being stored in liver and muscle tissue (‘glycogen synthesis’) when the blood glucose level is high. This last process can be modelled as a second order chemical reaction leading to a decrease of blood glucose at a rate δxy, in which the coefficient δ is related to insulin resistance R = 1/ δ [18]. The dynamics of glucose thus reads (1) where I(t) denotes the inflow of glucose in the blood circulation through the wall of the gastrointestinal tract and E(t, x) the transformation of blood glucose into energy that is expended by the body.

Insulin, variable y, is produced by the β-cells in the pancreas at a rate that is assumed to depend linearly upon the blood glucose level. Insulin is cleared from the blood circulation by liver and kidneys at a constant rate. Furthermore, insulin decays linearly in blood with a half-life in the order of 5 minutes. Matthews et al. [18] modelled the homeostatic relation between blood glucose and insulin concentrations at low glucose concentrations with the insulin concentration y in mU/L. Using their result we formulate a differential equation for the dynamics of insulin: (2ab) where the decay parameter takes the value λ = 8.3 h^-1. The parameter β represents the β-cell function. In [18] a value of 3.5 mmol/L in plasma is found for x₀. Taking into account that our data is in whole blood we set x₀ = 3 mmol/L. Furthermore, from Fig 1 of [18] it is found that at a 100% β-cell function β/λ = 4.2 giving β = 34.9 mU h^-1mmol^-1.

Glucagon, variable z, is produced by the α-cells in the pancreas. A low blood glucose level corresponds with an increased glucagon production. In [19] it is stated that insulin plays an intermediating role in this respect: a high insulin production inhibits glucagon production by α-cells. For the expression of the production of glucagon as it depends on the insulin production we take an exponential function. It requires only one parameter (m) and exhibits a close to linear decay. Moreover, for larger values of the insulin production it will tend to zero while remaining positive. This leads to a differential equation for glucagon dynamics: (3ab)

The parameter α denotes the maximal production of glucagon by the α-cells. Glucagon has similar to insulin a half-life on the order of minutes; the parameter γ represents this process of decay. The expression (3b) indicates that glucagon production takes larger values when the insulin production is decreasing at lower glucose concentrations.

For the function I(t) the profile of a Gamma-distribution is taken with a multiplicative factor w denoting the total amount of glucose in the meal: (4)

It is assumed that during digestion the expenditure of calories is constant. Moreover, the fraction of expended glucose based calories depends on the blood glucose concentration: the higher the concentration the larger the fraction. The following functional relation is taken (5) so that at the upper boundary of feasible blood glucose concentrations, x = 10 mmol/L, the fraction of glucose based expended calories is about 99%. The dynamics of (1) is steered by the forcing term I(t; w, a, k) so it acts at a time scale of hours. The dynamics of (2) and (3) act at a time scale of minutes with half-lives on the order of 5 minutes [20–21]. This large difference in times scales between (1) and at the other hand (2) and (3) allows a quasi-steady state approximation [13] meaning that for (2) and (3) the equilibrium can be taken as an estimate for the solution, so (6ab)

Fitting parameters from data

The variables y and z can be removed from (1) by using (6ab); this changes (1) into (7) with (8abc)

It is noted from (8a) that from blood glucose data the product of the parameters δ and β can be estimated using the estimate of ω. However, the individual parameters are not identifiable. If we leave out the intake of food and the effect of physical exertion, the dynamics of blood glucose is governed by the last two terms of the right hand side of (7). These make the system tend to an equilibrium x = x_b, called the baseline. It is remarked that other hormones such as cortisol and adrenalin still may play a role and affect the estimation of the parameters, in particular σ, so that in models with more detail the values of these parameters will change.

Next we construct with (7) a blood glucose profile brought about by meal 1 given in Table 1 and Fig 5 of [7], where during digestion of a fast absorption meal venous blood glucose values were measured from 23 nondiabetic subjects. The meal had a total glucose content of w = 50.5 g, or 280.3 mmol. In the figure average values are given.

We assume that q = 100 mmol/h being equivalent with a glucose based energy expenditure of about 75 kcal/h. For the parameter n we choose the value 1.75, see [22]. At the starting and end point of the measurements this system is in the equilibrium state meaning that the effect of the forcing I(t; w, a, k) is not felt. It corresponds with a baseline value x_b = 4.8 mmol/L. For given parameter values the solution of (7) is found by numerical integration. The parameters a, k, ω and σ are varied in order to find the values of these parameters for which the least square method yields the best solution. Using the Matlab based software package Grind (www.sparcs-center.org/grind) the best fit is found for (9) see Fig 1. Using (8a) it is concluded that at a 100% β-cell function (β = 34.9) the parameter δ = 2. We introduce the notion of insulin resistance by R = 1/δ = 0.5 h mU / L.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Fitting the solution of differential Eq (7).

The data (o) is from a meal which has rice as the main component (Fig 5 of [7], ∆ Meal 1) with w = 280.3 mmol. The parameter estimates are given in (9).

https://doi.org/10.1371/journal.pone.0169135.g001

In S1 Appendix, a sensitivity analysis for the peak value and the base line has been carried out. Two methods were utilized: Pearson Rank Correlation Coefficients (PRCC values) and extended Fourier Amplitude Sensitivity Testing (eFAST).

Dependence upon β-cell function and glucose content of the meal

The dependence of the peak value x_max upon the size w of the carbohydrate component of the meal is worth to be consided over a larger range of these parameters. In Fig 2A the results are given. There is a linear relation that holds for w within the interval [0, 500]. For a meal with given value w^(m) this line is determined by two points: (w^(m), ) and (0, x_b), where x_b is the baseline value. Such a relation only holds for meals with a comparable composition of carbohydrates.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Peak value x_max as a function of the size w of the carbohydrate component of the meal and the parameter ω.

(a) Parameter w is varied while the other parameters are fixed (solid line). The peak value that corresponds with the meal (w = 280.3) is given by (●). A linear regression (dashed line) is made after deleting the values w ≥ 500: x_max = 4.7 + 0.01 w. (b) Dependence of x_max upon ω. The value of x_max for (9) is given by (●). For ω = 2 the peak value is 12.8 (■); this can be brought down to 9.7 by doubling ω (□), see the section on closed-loop blood glucose control.

https://doi.org/10.1371/journal.pone.0169135.g002

The dependence of x_max upon the parameter ω plays a main role in the control of the blood glucose level, see the sensitivity analysis in S1 Appendix. It is noted, see Fig 1, that a value of ω = 8.30 yields for the given meal a blood glucose concentration that stays within the safe interval [3, 10]; in Fig 2B the corresponding peak value is indicated by (●). In case the β-cell function falls short and a dosage of exogenous insulin is mandatory, it is important to know the decrease of the peak x_max as a function of the amount of exogenous insulin. In case of diabetes (ω = 2) a peak value of 12.8 will occur for the meal that was presented (■). This can be brought down to 9.7 by doubling ω (□). In the next section it is worked out how this technically can be achieved.

Closed-loop blood glucose control

In the case of diabetes exogenous insulin can be brought into the blood circulation in order to maintain a safe blood glucose level. The external control variable u(t) introduced in (2a) changes (6a) into (10) so that (7) becomes (11) with (12abc)

There is a vast literature on the control problem of finding the function u(t) that keeps the blood glucose concentration at a safe level [2, 23–24]. One group of investigations MPC (Model Predictive Control) takes a mathematical model as a starting point [4–6]. The nonlinear model (11) could be used for that purpose in a control algorithm such as the Kalman filter [25]. However, being aware that a high blood glucose level is due to a low value of the parameter ω, a natural choice is to take (13) with ρ ≥ 0 and b(x, x₀) given by (2b). It is noted that from the total inflow of insulin in the blood circulation at all times a fixed fraction (14) comes from the infusion device. This type of external feedback is similar to the power control mechanism of the pedelec (electric bicycle). The approach differs from control algorithms such as PID (Proportional-Integral-Derivative control) that use a set point, e.g. the baseline [5–6]. If our model were perfect, then ρ could be kept at a constant value being sufficiently large to let the blood glucose concentration stay at all times within the safe range. Moreover, the body will experience coherent fluctuations of glucose and insulin in the blood circulation comparable with those of a nondiabetic person. This may bring about favourable effects such as a reduced risk of (further) insulin resistance [26]. If the infusion consists of long-acting insulin with a slower decay in the blood circulation (λ₁ < λ), then (13) changes into (15)

Additional functions can be incorporated in the insulin infusion device, such as range control by cutting off dosage below a certain blood glucose level and a temporal increase of ρ at extremely high blood glucose values. For β > 0 fine tuning can also be done by letting the amplification factor μ depend upon x and its time derivative : (16)

By increasing the amplification factor for > 0 the system will anticipate a larger need of insulin and compensate for delays in the control system.

In order to further work out the closed-loop control for a diabetic case, we take the case (9) with now ω = 2, assuming that there is no insulin resistance. Thus, we have δ = 2, so that κ = 0.24 and β = 8.3. The blood glucose peak value of 12.8 in Fig 2B is brought back to 9.7 (ω = 4) for ρ = 8.33. Thus, in that case the insulin infusion has to be at about 25% of the production by the pancreas of a healthy person. For meals with a different value of w the value of ρ, that brings the peak value x_max below 10 mmol/L, can be found by using the linear dependence of the peak value upon the size w of the meal (Fig 2A).

The parameter κ contains the unidentifiable parameter δ. This problem is solved by repeating the intake of a meal, e.g. the one with parameter values (9), with an external control (13) for which ρ takes a given fixed small positive value. Then the new estimate of the parameter ω will take a larger value ω_c, so κ = (ω_c—ω)/ρ. This yields with (12c) a value for the insulin resistance R = 1/δ.