The mathematical model
Suppose that C = C(t) and D = D(t) represent the numbers of diabetics with and without complications, respectively, and let N = N(t) = C(t) + D(t) denote the size of the population of diabetics at time t (see Nomenclature). Then, as was noted earlier, N(t) ~ 3% of the world population. Let I = I(t) denote the incidence of diabetes mellitus. The model parameters to be incorporated are μ (the natural mortality rate), λ (the probability of a diabetic person developing a complication), γ (the rate at which complications are cured), ν (the rate at which diabetic patients with complications become severely disabled) and δ (the mortality rate due to complications).
A schematic representation of the model is shown in Figure 1.
Figure 1: The mathematical model
The diagram shows that I = I(t) cases are diagnosed in a time interval of length t and are assumed to have no complications upon diagnosis. In that same time interval, the number of sufferers without complications, D = D(t), is seen to decrease by the amounts μD (natural mortality) and λD (sufferers who develop complications), and to increase by the amount γC (sufferers whose complications are cured). During this time interval, the number of diabetics with complications is increased by the afore-mentioned amount γC and by the amount μC (natural mortality), νC (patients who become severely disabled and whose disabilities cannot be cured) and δC (those who die from their complications).
These rates of change are formalized by the ordinary differential equations (ODEs)
which, since N(t) = D(t) + C(t), give rise to the initial-value problem (IVP)
C'(t) = -(λ + θ)C(t) + λN(t), t > 0; C(0) = C0 (1)
N'(t) = I(t) - (ν + δ)C(t) - μN(t), t > 0; N(0) = N0 (2)
where θ = γ + μ + ν + δ, 
and C0, N0 are the initial values of C(t) and N(t), respectively.
In the case when the probability of a diabetic person developing a complication, λ, is constant, the model equations (1), (2) are linear in C(t) and N(t): this linear model will be discussed in the following paragraph. The non-linear model corresponding to a variable λ will be considered by the authors in another paper more devoted to numerical analysis.
The linear case
The critical point and its stability property
The probability of developing a complication, λ, will be estimated to have the constant value
The initial-value problem (1), (2) may consequently be written in matrix-vector form as
x'(t) = Ax(t) + b(t), t > 0;x(0) = X0 (4)
in which
Suppose that I is the steady-state value of the incidence, then the model reaches its critical point when dC/dt and dN/dt given in (1) and (2) vanish simultaneously, that is when
λN - (λ + θ)C = 0, (6)
I - μN - (ν + δ)C = 0. (7)
Solving (6) and (7)gives
and
The eigenvalues of the matrix A, χ1, χ2, are the roots of the quadratic equation (the characteristic equation)
χ2 + (λ + θ + μ)χ + μ(λ + θ) + λ(ν + δ) = 0. (9)
The discriminant, Δ, of this equation is given by
Δ = (λ + θ + μ)2 - 4 [μ(λ + θ) + λ(ν + δ)]
and, recalling that θ = ν + μ + δ + γ it follows that
Δ λ + μ + δ + γ + 2ν.
Solving (9) gives
and it is then easy to check that, when the parameters of the model are such that
(a) Δ > 0, χ1, χ2 are both real and negative;
(b) Δ = 0, χ1 = χ2 are real and negative;
(c) Δ χ1 and χ2 are complex conjugate with negative real parts.
it may be concluded, therefore, that the critical point (C*, N*) of (1), (2), given by (8), is stable.
Numerical solution and stability
It may be shown that the solution x(t) of the IVP (4) satisfies the recurrence relation
where l > 0 is an increment in t (the time step). This recurrence relation may be used to generate x(tn + 1) in terms of x(tn), thus monitoring C(t) and N(t) at the discrete points t = tn = nl(n = 0, 1, 2, ...).
One very simple way of estimating x(t + l) is to approximate to second order the integral in (11) by the trapezoidal rule, viz.
and then to replace, also to second order, exp(lA) in (11) and (12) by its (1,1) Padé approximant
exp(lA) = (E - 1/2lA)-1 (E + 1/2lA) (13)
where E is the identity matrix of order two.
Denoting by Xn the numerical approximation to x(tn) calculated using (11), (13), it may be shown, by substituting (12) with (13) in (11) and then by pre-multiplying by (E - 1/2lA), that
(E - 1/2lA)Xn + 1 = (E + 1/2lA)Xn (14)
+
1/2l[(E - 1/2lA)bn + 1 + (E + 1/2lA)bn];n = 0, 1, 2, ...
where Xn = (Cn, Nn)T, T denoting transpose, and bn = (0, In)T with In = I(tn). It may then be shown that Cn + 1 and Nn + 1 (n = 0, 1, 2,...) may be determined by solving the algebraic equations given by
(Method 1) (1 + 1/2l(λ + θ))Cn + 1 - 1/2lλNn + 1 =
[1 - 1/2l(λ + θ)]Cn + 1/2lλNn - 1/4l2λ(In + 1 - In) (15)
and
1/2l(ν + δ)Cn + 1 + (1 + 1/2lμ)Nn + 1 =
-1/2l(ν + δ)Cn + (1 - 1/2lμ)Nn
+1/2l(1 + lμ)In + 1 + 1/2l(1 - lμ)In (16)
assuming convergence, Cn + 1 = Cn = C, Nn + 1 = Nn = N and In + 1 = In = I, say, then equations (15) and (16) become
(λ + θ)C - λN = 0, (17)
(ν + δ)C + μN = I, (18)
respectively. Obviously (17) and (18) are the same as (6) and (7) and so the fixed point (C+, N+) of the numerical solution sequence (Cn, Nn), n = 0, 1, 2,... is the same as the critical point (C*, N*) of the linear initial value-problem.
For comparison purpose, the IVP (1), (2) was also solved using the well-known Euler method (a first-order method) given by
(Method 2) Cn + 1 = [1 - l(λ + θ)]Cn + lλNn (19)
Nn + 1 = -l(ν + δ)Cn + (1 - lμ) + lINn (20)
The method1 is unconditionally stable whereas the Euler method is conditionally stable provided ([33])
Numerical experiments
Taking I(t) = I to be constant equations (15) and (16) simplify to
(1 + 1/2l(λ + θ))Cn + 1 - 1/2lλNn + 1 =
[1 - 1/2l(λ + θ)]Cn + 1/2lλNn (22)
1/2l(ν + δ)Cn + 1 + (1 + 1/2lμ)Nn + 1 =
-1/2l(ν + δ)Cn + (1 - 1/2lμ)Nn + lI, (23)
respectively, for n = 0, 1, 2, .... In the numerical experiments, I was given the value 6.106yr-1 and the parameters ν, δ, μ, γ, were given the value shown in Table 2. the critical values of C and N were then calculated from (8) and were found to be
C* = 47000000 and N* = 61100000. (24)
Using Matlab, four numerical experiments were carried out taking as initial conditions
C0 = C* ± 500 and N0 = N* ± 500. (25)
A time step of l = 0.01yr-1 was used and the solution to the IVP (1), (2) was computed by solving (22), (23) for n = 0, 1, 2, .... Using all combinations of initial conditions, the computed solution converged to the values of C* and N* given in (24). By way of example, the fixed points, C+ and N+, to which the numerical solution converged are shown in Table 3.
The initial conditions in (25) are close to the steady-state solutions C* and N*. Other initial conditions, further from C* and N*, will converge to the same values of C+ and N+ for the same value of l, though these values will be reached at different times.
Retaining the parameters values shown in Table 2 the effect of the choice of time step was monitored in a series of 11 further experiments. The fixed-point values, C+ and N+ to which convergence occurred are shown in Table 3, where it may be seen that, for the larger values of l (l ≥ 0.5yr), there is very close agreement with the critical-point values given in (24).
For l ≤ 2yr (approximately), the two methods give similar results but for l > 2.5yr (approximately (21)) the Euler method diverged. The values of C+ and N+ using Euler method are given in Table 3.
It may be concluded from these results that the Euler method may be used with confidence if the diabetic population is to be monitored at time intervals up to approximately two years using the linear model (1), (2). However, to monitor the population less frequently the numerical method (Method 1) should be used.
Answers to all your Biology Questions
