ABSTRACT

Worldwide, diabetes is affecting 370 million people, causing nearly five million deaths and absorbing more than 471 billion USD per year. Mathematical models have been developed to simulate, analyse and understand the dynamics of β-cells, insulin and glucose. In this paper, we consider the effect of genetic predisposition to diabetes on dynamics of β-cells, glucose and insulin. We assume that the β-cell dynamics is governed by the differential equation: . The model indicates different behaviours according to the presence or absence of genetic predisposition. In presence of predisposition (ε = 1), the model shows three equilibrium points: a stable physiological equilibrium point (G = 100, I = 20, β = 600), a stable trivial pathological equilibrium point (G = 600, I = 0, β = 0) and a saddle point (G = 250, I = 9.8, β = 129.36). In absence of predisposition (ε = 0), the model has only two equilibrium points: an unstable pathological equilibrium point (G = 600, I = 0, β = 0) and a stable physiological equilibrium point (G = 82.6, I = 23, β = 900). In order to see how physical activity, obesity and other factors affect insulin sensitivity, simulations are carried out with different values of insulin induced glucose uptake rate (c), β-cell maximum insulin secretory rate (d) and environmental capacity (K).

Keywords:Type 2 Diabetes, Mathematical Modelling, Physiological Equilibrium, Pathological Equilibrium

1. Introduction

Once associated with economic development and considered as a disease of the rich, diabetes is now affecting countries worldwide, threatening particularly low and middle-income countries. According to the last International Diabetes Federation (IDF) report, about 80% of the 371 million diabetics live in lowand middle-income countries; nearly five million deaths were due to diabetes; and more than 471 billion USD was spent on healthcare for diabetes in 2012 [1] . Due to its chronic nature with severe complications, diabetes needs costly prolonged treatment and care, affecting individuals and societies, and raising the equity problem between and within countries [2] . Although prevalence of diabetes is increasing all over the world, the Middle East and North Africa region has the highest comparative prevalence of diabetes (11%). Six of the top 10 countries with the highest prevalence of diabetes (in adults aged 20 to 79 years) are in this region: Kuwait (21.1%), Lebanon (20.2%), Qatar (20.2%), Saudi Arabia (20%), Bahrain (19.9%) and UAE (19.2%) [1] [3] . Consequently, the direct and indirect socio-economic burden of diabetes is exponentially increasing in this region [4] [5] .

During the last decades, a large number of mathematical models have been developed to simulate, analyse and understand the dynamics of glucose and insulin leading to diabetes. Bolie (1961) was among the pioneers in this field. He introduced a simple linear model, using ordinary differential equations in glucose and insulin [6] . The most frequently used model is the well known minimal model published in early eighties by Bergman et al. [7] . Variant and critical versions based on the minimal model were published by different authors, including Derouich and Boutayeb (adding physical effort) [8] , De Gaetano and Arino (using delayed differential equations) [9] , Li et al. (proposing a generalised model) [10] and Roy and Parker (extending the minimal model with free fatty acids) [11] . Other authors dealt with the dynamics of β-cells and the mechanisms leading to diabetes. Topp et al. published an interesting model on pathways to diabetes by considering β-cell mass, insulin, and glucose kinetics [12] . Hernandez et al. proposed an extension of the Topp model by adding the surface insulin receptor dynamics [13] . Using a mathematical model similar to the model of Topp et al., De Gaetano and collaborators introduced the concept of pancreatic reserve [14] . Gallenberger et al. proposed a model on the dynamics of glucose and insulin concentration connected to the β-cell cycle [15] . More generally and for modelling details, the reader is referred to recent reviews that dealt explicitly with different models using differential equations, delayed differential equations, integro-differential equations, stochastic differential equations, optimal control and other methods for glycaemic control, blood glucose monitoring and devices devoted to diabetic prevention [16] -[22] . In this paper, we propose a mathematical model stressing the effect of genetic predisposition to diabetes on the dynamics of β-cells, insulin and glucose.

2. Materials and Methods

2.1. Predisposition to Type 2 Diabetes

Obesity is generally associated with insulin resistance and hence represents a major risk factor for the development of type 2 diabetes [23] -[25] . However, not all obese people develop diabetes. In some of them, the β-cells mount a compensatory response to counter insulin resistance and the magnitude of compensation is probably genetically determined [25] . More generally, in lean and obese people not predisposed to develop diabetes, normoglycaemia is preserved since β-cells are able to sustain an adequate compensation to insulin resistance. In contrast, in people genetically susceptible to diabetes, the compensatory response is insufficient and β-cell decompensation ensues [23] . Following the discussion of Poitout and Robertson on glucolipotoxicity and their hypothesis on the progressive β-cell dysfunction in lean and obese phenotypes with genetic predisposition [23] , we propose a mathematical model of β-cell mass, glucose and insulin, extending the model of Topp et al. by adding the effect of genetic predisposition to diabetes.

2.2. The Mathematical Model

The mathematical model proposed is the following:

(1)

(2)

(3)

Most of mathematical models for glucose dynamics assume that the concentration of glucose in the blood is determined by a differential equation of the form [7] -[9] [11] -[13] . In other words, the concentration of glucose increases by a rate (glucose production by liver and kidneys) and decreases by a rate (independent of insulin) and a rate representing the glucose uptake due to insulin sensitivity.

However, as indicated by Li et al., the relation assumes a mass-action law while this term may be replaced by a more general expression [10] . Consequently, we use the Michaelis function, assuming that insulin-dependent net glucose tissue uptake takes the form: where is the half saturation constant, and is the maximum capacity. We assume thus that the variation of glucose is governed by the differential equation:.

For insulin dynamics, we keep the same expression used in the model of Topp, assuming that the net insulin secretion rate can be modelled as a sigmoidal function of glucose level while insulin clearance is of the form

. Using a Hill function of the form as a sigmoid function reaching half its maximum at, the insulin dynamics takes the form: where is the mass of pancreatic β-cells (mg), and secretion of insulin is supposed uniform for all β-cells, at the same maximal rate.

Finally, as stressed earlier, we assume that the β-cell dynamics depends on genetic predisposition to diabetes and is, consequently governed by the differential equation:

This equation indicates that when ε = 1 (predisposition to diabetes) we simply get the expression used in the model of Topp, when ε = 0 (no predisposition to diabetes) we get a logistic equation in and for ε between 0 and 1 we get different combinations according to the level of predisposition. In the first right-hand term, is the growth rate of β-cells and K (mg) represents the carrying capacity while the second right-hand term is the second-degree polynomial in used by Topp, meaning that the β-cell mass decreases outside the polynomial zeros and increases between the zeros (when the polynomial has two real zeros).

3. Results and Discussion

3.1. The Model Behaviour

Taking ε = 1 and using the values of parameters given in Table 1 yields results similar to those given by Topp, namely a physiological equilibrium point (G = 100, I = 20, β = 600), a trivial pathological equilibrium point (G = 600, I = 0, β = 0); and a saddle point (G = 250, I = 9.8, β = 129.36).

When ε = 0, our model has only two equilibrium points. A physiological equilibrium point (G = 82.6, I = 23, β = 900) and a trivial pathological equilibrium point (G = 600, I = 0, β = 0). For ε = 0.5, the model has three equilibrium points comparable to those obtained with ε = 1.

Stability analysis (see Appendix) shows that the physiological equilibrium point is stable for all values of ε. The pathological equilibrium point is stable when ε = 1 or 0.5 and more generally for values of ε different of 0 but it is unstable when ε = 0. The third equilibrium point is unstable when it exists (for values of ε different of 0) (Table 2).

3.2. Discussion

Assuming predisposition to diabetes (ε = 1), the qualitative behaviour of our model is similar to that of the models developed by Topp et al. [12] , Hernandez et al. [13] , and De Gaetano et al. [14] represented by three equilibrium points: a stable physiological steady state (basal values of glycaemia and insulin); a trivial stable pathological severe diabetes steady state corresponding to a hyperglycaemic state with zero levels of β-cell mass and insulin, and an unstable saddle point with intermediate values of glycaemia, insulin and β-cell mass. Topp et al. and Hernandez et al. defined the bifurcation pathway to diabetes as a combination of parameters leading to the elimination of the physiological and saddle fixed points and making the pathological equilibrium point as a global attractor [12] [13] . De Gaetano et al. indicated that, when the saddle point exceeds a threshold, glucose toxicity prevails and β-cell mass eventually declines to zero, leading to the pathological steady state with zero insulin production [14] . The main difference of our model with respect to the previous models is the introduction of the effect of genetic predisposition to diabetes. Indeed, when ε = 0 (no predisposition to develop diabetes), our model gives only two equilibrium points: the trivial pathological point and the physiological steady state. Moreover, the pathological equilibrium point is unstable, leaving the physiological equilibrium point (with more realistic values for insulin (23) and β (900)) as the global attractor indicating that diabetes will not be developed.

When (ε = 1), the use of the Michaelis function in the first differential equation indicates an improvement of the results (compared to those of Topp) in the case of the physiological equilibrium point (Figure 1). Indeed, the

value obtained for β (600) seems more realistic than the value obtained by Topp et al. (300). In fact, according to Hernandez et al., the average mass of β-cells in a normal individual is around 850 mg [13] . Another study devoted the β-cell mass in European subjects found that β-cell mass was in average about 35% lower in type 2 diabetics (573 mg) than in non diabetic subjects (888 mg) [26] [27] .

3.3. Simulation with Different Values of Parameters

In order to see how the parameter in Equation (1) affects insulin sensitivity, a simulation was carried out with different values of c, d and K.

Taking into account the results of studies showing that physical exercise can increase insulin sensitivity by 36% [13] , we used a simulation with the corresponding increased value of c (0.98) and obtained physiological equilibrium points (G = 100, I = 14.7, β = 441) with lower values of insulin and β-cell mass when ε = 1 and (G = 72, I = 18.5, β = 900) with a lower value of glucose when ε = 0 (Figure 2) (Table 3).

On the other hand, insulin resistance due to obesity and/or other factors has been shown to decrease insulin induced glucose uptake by 50% - 100% in both diabetic and non-diabetic subjects [12] [13] . In this direction, a simulation was used with a value of c decreased by 60% (c = 0.288), giving the physiological equilibrium point (G = 100, I = 50, β = 1500) with normoglycaemia associated with higher values of insulin and β-cell mass when ε = 1 and the equilibrium point (G = 130, I = 41.4, β = 900) when ε = 0 with values of insulin and β-cell mass lower than that of the case ε = 1 but with a value of glucose greater than normoglycaemia (Figure 3) (Table 4). At a first glance, the value of G = 130 when ε = 0 seems contradictory with the fact of no predisposition to diabetes. However, the adaptation of β-cell suggests that, in this case, an increase of the β-cell mass (K) or the β-cell maximum insulin secretory rate (d) will assure near normoglycaemia as indicated by the following simulations:

1) Following Rahier et al. who found that the β-cell mass increased with BMI in both diabetic (25%) and nondiabetic individuals (20%) [26] , a third simulation was used with high level of insulin resistance and high carry-

ing capacity (K) giving the physiological equilibrium point (G = 100, I = 50, β = 1500) when ε = 1 and (G = 116, I = 45, β = 1125) when ε = 0 (Table 5).

2) Assuming a higher β-cell insulin secretory rate as a compensatory response to insulin resistance, a fourth simulation was used with high level of insulin resistance and a higher level of the β-cell maximum insulin secretory rate (d = 60), giving the physiological equilibrium point (G = 100, I = 50, β = 1080) when ε = 1 and (G = 110, I = 47, β = 900) when ε = 0 (Table 6).

3) Finally, a fifth simulation with high level of insulin resistance compensated by both a higher carrying capacity (K) and a higher level of the β-cell maximum insulin secretory rate was used, giving the physiological equilibrium point (G = 100, I = 50, β = 1080) when ε = 1 and (G = 98, I = 51, β = 1125) when ε = 0 (Figure 4) (Table 7).

4. Conclusions

In this paper, we considered the effect of genetic predisposition to type 2 diabetes on the dynamics of β-cells, insulin and glucose. Our model shows that in absence of predisposition to diabetes (ε = 0) the pathological equilibrium point is unstable and evolution will be towards the stable physiological equilibrium state (G = 82.6, I = 23, β = 900).

Assuming that physical exercise can increase insulin sensitivity by 36% whereas overweight/obesity and other factors may decrease insulin induced glucose uptake in both diabetic and non diabetic subjects, simulations were carried out with different values of the parameters c, d and K. The model shows that increased insulin sensitivity (higher value of c) due to physical exercise leads to a physiological equilibrium point with normoglycaemia and a low value of β-cell mass when ε = 1 and a value of glucose approaching normoglycaemia (G = 72) and a normal value of β-cell mass when ε = 0. On the other hand, insulin resistance (lower value of c) can be compensated by an increased β-cell insulin secretory rate (d) or by an increased level of β-cell carrying capacity (K). In particular, when both the values of K and d are increased, the model gives a physiological equilibrium point with nearly the same values of G, I and β for all values of ε. Consequently, our model stresses, as shown by experimental studies, the fact that not all people with insulin resistance develop diabetes. More interestingly, our model confirms the idea that non diabetic and especially pre-diabetic people may avoid the evolution towards type 2 diabetes by acting on the risk factors, especially physical activity and overweight/obesity.

Author’s Contributions

The three authors participated equally to this paper. They read and approved the final manuscript.

Competing Interests

The authors declare that they have no competing interests.

References

Appendix

The equilibrium points satisfy the equations:

(4)

(5)

(6)

When ε = 1, we can easily compute the three equilibrium points of the system:; and.

Where

When, the trivial equilibrium point is easily obtained, whereas the non trivial equilibrium points are computed using numerical approximation by Maple:

When,.

When, and.

Stability Analysis

The Jacobian matrix of the system is given by:

Consequently, stability analysis of the trivial equilibrium point is straightforward.

When ε = 0, P₀: (G = 600, I = 0, β = 0) is unstable since has two negative eigenvalues: λ₁ = −b and λ₂ = f but the third eigenvalue (λ₃ = r) is positive.

When ε = 1, P₀: (G = 600, I = 0, β = 0) is stable since has three negative eigenvalues: λ₁ = b, λ₂ = −f and.

When ε = 1, P₀: (G = 600, I = 0, β = 0) is stable since has three negative eigenvalues: λ₁ = b, λ₂ = −f  and.

For the non trivial equilibrium points, eigenvalues are numerically computed using Maple (Table 2).

NOTES

^*Corresponding author.