Geometric Singular Perturbation Analysis of a Multiple Time-scale

Model for Diabetes and COVID-19 Comorbidity

CHONTITA RATTANAKUL1,2, YONGWIMON LENBURY1,2,

NATHNARONG KHAJOHNSAKSUMETH1,2, CHARIN MODCHANG3,4

1Department of Mathematics, Faculty of Science, Mahidol University, Rama 6 Rd., THAILAND

2Centre of Excellence in Mathematics, MHESI, Bangkok 10400, THAILAND

3Department of Physics, Faculty of Science, Mahidol University, Rama 6 Rd., THAILAND

4Centre of Excellence in Physics, MHESI, Bangkok 10400, THAILAND

Abstract: More and more information on mortality and morbidity indicates that in order to fight the COVID-19

pandemic, it is important to focus our attention on comorbidities. Several reports evidence of how many elderly

patients who become severely ill exhibit underlying illness such as cardiovascular disease, kidney disease, tumor,

and more to our special attention here, type 2 diabetes. Better understanding of the mechanism underlying the

comorbidity between different diseases requires merging models of systems across different time-scales. The

model homogenization across multiple spatial and time scales poses an important challenge to researchers in the

field of medical science. An approach that has been found relatively efficient in the analysis of such models is the

use of singular perturbation technique. Here, we study a differential equation model system with multiple time

scales which describes the diabetes and COVID-19 comorbidity. It tracks the changes in levels of plasma glucose,

insulin, and functional-cells, incorporating insulin resistance and inflammation responses. The model is analyzed

with the geometric singular perturbation technique, by which conditions on the system parameters may be derived

to identify regions in which the system exhibits different dynamic behavior, whether the system would be stable,

or eventually oscillate in a sustained fashion. Discussion of these conditions allows us to better understand how

comorbidity mediates the development of life-threatening symptoms in a diabetic patient in order that proper care

and treatment may be prescribed.

Key-Words: Diabetes model, COVID-19, comorbidity, geometric singular perturbation analysis

Received: March 31, 2022. Revised: July 4, 2022. Accepted: August 22, 2022. Published: September 5, 2022.

1 Introduction

The spread of coronavirus across the globe has devel-

oped into a major pandemic which impacts all aspects

of the human society. A great deal of research effort

has been devoted to modeling and prediction of the

progress of coronavirus infection and its impacts of

social distancing and other measures in various coun-

tries hit hard by the pandemic, for example in [1], [2]

and [3].

According to Smith [4], medical doctors have been

investigating an unexpected coronavirus complica-

tion that can appear a couple of weeks or sometimes

months, after the initial COVID-19 infection. In-

creasing number of patients have been reported to

develop life-threatening symptoms that necessitate

prompt medical attention. This new complication of

corona virus is dangerous and persists throughout a

patient's life since, although both type 1 and type 2 di-

abetes can be managed with treatment, diabetes can't

be cured, which means some patients will need ex-

pensive insulin therapy for the rest of their lives [4].

Although most patients will survive COVID infec-

tion, many of them experience long-term symptoms

that frequently require medical attention. Moreover,

quite a few even progress to multisystem inflamma-

tory syndrome, according to Smith[4].

According to Campbell[5], patients with type 1 di-

abetes are more likely to die of Covid-19 than those

with type 2. Also, NHS research confirmed that di-

abetes in coronavirus patients significantly increases

the risk of their dying. "Almost one in three of all

deaths from coronavirus among people in hospitals

in England during the pandemic have been associated

with diabetes," according to the study[5].

In an attempt to gain a better understanding of such

complication, we consider here a model system of dif-

ferential equations that describes the glucose-insulin

dynamics, based on that proposed by De Gaetano et

al. [6] and those mentioned in the work on the compu-

tational patient having diabetes and COVID by Barbi-

ero and Lió[7], incorporating the effects of inflamma-

tory response and insulin resistance, which we have

modified to take into account the drastically

different time scales observed in the time variations

WSEAS TRANSACTIONS on BIOLOGY and BIOMEDICINE

DOI: 10.37394/23208.2022.19.20

Chontita Rattanakul, Yongwimon Lenbury,

Nathnarong Khajohnsaksumeth, Charin Modchang

E-ISSN: 2224-2902

176

Volume 19, 2022

of the concentrations of glucose, insulin, the number

of functional β-cells and the inflammatory response.

The presence of glucose (glycaemia) and insulin in

one's blood is regulated through a negative feedback

loop in which plasma glucose stimulates β-cells to re-

lease insulin giving rise to insulin-mediated increase

in tissue glucose uptake while liver gluconeogenesis

and glycogenolysis are decreased [6].

In 2001, a nonlinear mathematical model of the

glucose-insulin control mechanism was proposed by

Lenbury et al.[8], taking into account the role of β-

cells in maintaining and regulating insulin level in hu-

man's bloodstream. They utilized a gastrointestinal

term to model the absorption of glucose by the intes-

tine, and the entry of glucose into the human's plasma,

under the assumption that this process declines expo-

nentially with time. The model is analyzed using the

geometric singular perturbation method by which var-

ious conditions on the system parameters are obtained

to identify different dynamic behavior exhibited by

the system, including the possibility of limit cycle be-

havior in the system model which closely simulates

sustained oscillatory time series often observed in ex-

perimental data. To take into account the temporal ab-

sorption of glucose, they used a sinusoidal term with

the means to study how the patients respond when

subject to ambulatory-fed conditions. The resulting

model equations are as follows.

dt =r1zy −r2zx +c1z

dt =R3N

z−R4x+C2+w(t)

dt =R5(y−ˆy)(T−z) + R6z(T−z)−R7z

where xand yare the differences in the plasma insulin

and glucose concentrations from their respective fast-

ing (basal) concentrations, and zrepresents the densi-

ty of the pancreatic β-cells in the proliferative phase.

In 2008, Giang et al.[9] analyzed a delay model

of the glucose-insulin system proposed by Palumbo

et al. in [10]. They proved persistence, as well as

existence and stability of a unique positive equilibri-

um point of the model system. Moreover, their work

[9] illustrates uniform persistence and global stabil-

ity of equilibrium solutions. Using omega limit set

of a persistent solution and the full time solution that

they derived, they were able to investigate the effect

of delays in terms of oscillating solutions. In addition,

global stability and slowly oscillating behavior were

shown to be possible provided suitable conditions on

the system parameters are satisfied. The model equa-

tions under consideration are as follows.

G(t) = −KxgG(t)−KxgiG(t)I(t−τi) + Tgh

I(t) = −KxiI(t) + TiGmax

f(G(t−τg))

where f(G) = Gγ

Gγ+Gγ

∗

with sup

G≥0

f(G) = 1,G

and Iare, respectively, glucose and insulin concen-

trations, with corresponding delays of τgand τi.

In the same year, Kardar et al.[11] made use of an

advisory/control algorithm in aide of out-patients suf-

fering from insulin dependent diabetes mellitus. Uti-

lizing Mamdani type fuzzy logic controllers, their ad-

visory/control algorithm combines expert knowledge

about diabetes treatments in the attempt to regulate

the glucose level in bloodstream, while the patien-

t is experiencing disturbances in glucose concentra-

tion due to food intakes or there are fluctuations in

the measured plasma glucose level due to measure-

ment errors.

Also in 2008, De Gaetano et al.[6] studied a mod-

el of the pancreatic islet compensation, and identified

certain basic qualitative characteristics of the model's

solutions. They simulated the model's performance

over one's lifetime under different conditions. They

then compared their model with two previous mod-

els of diabetes progression, and concluded that the

proposed model provided a realistic and robust de-

scription of the glucose-insulin process in healthy and

diabetic person. More importantly, by varying the

size of one of the system parameters, their model can

be reduced to either "fast", providing short-term pre-

diction, or "slow" subsystems, which describes long-

term progression.

In 2015, Cao et al.[13] studied a model which con-

sists of four states of type 2 diabetes, assuming no in-

put or output. Using a model consisting of partial

differential equations, a general well-posedness result

was obtained and the exponential stability of

dynamic solution was shown. The steady-state

solution predicts the diabetics' stable distribution

probability[13].

The models mentioned above are useful in answer-

ing a number of crucial questions concerning the pro-

gression of diabetes mellitus. However, in the model

analyses, they did not take into account the fact that

the state variables under consideration vary with time

at extremely different time scales.

Moreover, we observe that the model of De Gae-

tano et al.[6] has not properly accounted for the re-

ports that the ability of the cells to absorb or use up

blood sugar appears to decline as the glucose concen-

tration in the bloodstream increases to a high level.

It is mentioned in [14], that in a patient with predi-

WSEAS TRANSACTIONS on BIOLOGY and BIOMEDICINE

DOI: 10.37394/23208.2022.19.20

Chontita Rattanakul, Yongwimon Lenbury,

Nathnarong Khajohnsaksumeth, Charin Modchang

E-ISSN: 2224-2902

177

Volume 19, 2022

abetes conditions, the pancreas has to work harder to

release sufficient amount of insulin in order to control

the glucose in the bloodstream to a satisfactory low

level. It is well-known that insulin resistance takes

place when too much glucose in the blood decreases

the cells' ability to absorb and use blood sugar for en-

ergy [14]. We will therefore introduce a term in our

model, which is based on that in [6], to take into ac-

count the reduction in the ability of cells to control

plasma glucose when its concentration becomes too

high. For a patient with prediabetes conditions, we

express the insulin resistance term as a function of

glucose concentration, so that it rises as the glucose

level increases.

Recently, Ratanakul et al.[15] attempted to incor-

porate this insulin resistance effect more explicitly, as

well as account for the diversified time scales in their

model, which was analyzed by the geometric singu-

lar perturbation method, upon assuming that the level

of β-cells remains relatively constant and thus reduc-

ing the model analysis into a 2-dimensional singular

perturbation argument. In our present work, the func-

tional βcells is tracked, and although it is considered

to vary slowly, it is not taken to be constant in our

analysis, making it necessary to carry out the geomet-

ric singular perturbation analysis in the 3-dimensional

space.

To account for comorbidity between diabetes and

COVID-19, an equation describing the rate of change

in the inflammatory response is added to the mod-

el, following that given in the work by Barbiero and

Lió [7]. In this equation, the rate of change of in-

flammatory response is taken to depend on abnormal

ACE2 activity, SARS-CoV-2 affinity with ACE2, the

re-enforcement of inflammation due to drug surplus

left inside the body, inflammation rate due to glucose

surplus, and the anti-inflammatory response rate as

suggested in [7]

We will explain in subsequent sections that, based

on values of our model system parameters, reported

in various literature, we may conclude that the level

of plasma glucose varies relatively quickly, whereas

the number of β-cells changes at a much larger

time scale, and insulin concentration, on the other

hand, can be said to change at an intermediate speed

relatively. We shall then use two small parameters to

merge such diversified time-scales. Arguments based

on singular perturbation theory is thus the appropriate

geometric method of analysis which is applied to

identify different dynamic behavior permitted by our

model.

2 Materials and Methods

Insulin, I(t), is a chemical messenger produced in one

part of the body to have an action on another. It is a

protein responsible for regulating blood glucose lev-

els, G(t), as part of the metabolism. The control of

blood sugar, or glucose, by insulin is a negative feed-

back mechanism. When blood sugar rises, from our

food intake for example, receptors in the body sense

a change. In response, the control center, namely the

pancreas, is stimulated. In the pancreas, the function-

al beta cells, βf(t), are designed to function as "fuel

sensors" triggered by glucose[16]. As glucose levels

rise in the plasma of the blood, uptake and metabolis-

m by the pancreas beta cells are enhanced, leading to

insulin secretion into the bloodstream which has the

effect of lowering plasma glucose levels. Once blood

sugar levels reach a state of equilibrium (homeosta-

sis), the pancreas stops releasing insulin.

Thus, the drop in the rate of increase of plasma

glucose sends a negative feed-back to signal the func-

tional beta cells in the pancreas to reduce the produc-

tion of insulin. In this way, the glucose-insulin con-

trol system in a healthy human's body functions has

a negative feed-back loop that regulates the glucose

level in the blood stream.

We are thus led to the following model system of

differential equations, which has been adapted from

[6],[7], and [15].

dt =R0−EG0G−SI

IRS +i, G(t0) = G0

(1)

dt =σβfGνh

α+Gνh−kI, I(t0) = I0(2)

dβf

dt =−r0+r1G−r2G2βf, βf(t0) = βf0

(3)

dIRP

dt =kSARS ∆kACE2,0+kDDs+kGG

−keff IRP , IRP (t0) = IP R0(4)

The insulin resistance term IRS , in the denominator

of the last term in (1), is expressed as the following

function of glucose concentration:

IRS =a0G2+b0G+c0(5)

following the observation mentioned earlier that in-

sulin resistance occurs when excess glucose in the

blood reduces the ability of the cells to absorb and

use blood sugar for energy [15], so that the ability of

cells to control plasma glucose declines when glucose

concentration becomes too high. Here, a0is the sec-

ond order coefficient for insulin resistance, b0is the

first order coefficient for insulin resistance, and c0is

insulin resistance at zero glucose.

WSEAS TRANSACTIONS on BIOLOGY and BIOMEDICINE

DOI: 10.37394/23208.2022.19.20

Chontita Rattanakul, Yongwimon Lenbury,

Nathnarong Khajohnsaksumeth, Charin Modchang

E-ISSN: 2224-2902

178

Volume 19, 2022

Table 1: State variables and parameters appearing in

the model equations (1) - (5) with their descriptions

Item Description Unit

Gglucose concentration G∗

in bloodstream

Iinsulin concentration I∗

in bloodstream

βfdensity of functional β∗

βcells

IRS insulin resistance I∗

IRP inflammatory response I∗

ttime t∗

R0rate of glucose production G∗(t∗)−1

at zero glucose

SIinsulin sensitivity I∗

RS (t∗I∗)−1

EG0specific rate of glucose (t∗)−1

removal at zero insulin

iinsulin resistance self- I∗

inhibition coefficient

µ i +c0,c0being insulin I∗

resistance at zero glucose

σmaximal rate of insulin I∗(β∗

ft∗)−1

secretion by βcells

αglucose concentration (G∗)νh

yielding 50% of insulin

secretion

kcombined specific rate of , (t∗)−1

insulin uptake at the liver

kidneys, and insulin

receptors

νhpower coefficient for Hill- None

shaped glycemia effect on

pancreatic insulin release,

set to 4 as done in [6]

r0death rate of functional β∗

f(t∗)−1

βcells at zero glucose

r1I-order coefficient for β∗

f(G∗t∗)−1

βcell replication

r2II-order coefficient for (G∗)−2(t∗)−1

βcell replication

kSARS inflammation rate due I∗

RP (k∗

ACE2t∗)−1

to SARS-CoV-2

kACE2,0normal ACE-catalyzed k∗

ACE2

conversion rate

∆multiple of normal None

conversion rate

kDspecific inflammation rate I∗

RP (D∗

rt∗)−1

due to ACEi surplus

Dssurplus drug concentration D∗

kGspecific inflammation rate I∗

RP (G∗t∗)−1

due to glucose surplus

keff specific anti- (t∗)−1

inflammatory response rate

a0II-order coefficient for I∗

RS (G∗)−2

insulin resistance

b0I-order coefficient for I∗

RS (G∗)−1

insulin resistance

c0insulin resistance at I∗

zero glucose

The variables and parameters appearing in (1)-(5)

and their definitions are summarized in Table 1. The

unit of each of the variables or parameters is denoted

here as the star,(·)∗, of the corresponding quantity (·),

as it depends on the choice of units used by a particu-

lar researcher to measure these quantities.

In Equation (1) for the rate of change of glucose

concentration, following that given in [15], the first

term on the right, R0, is the rate of glucose produc-

tion at zero glucose. The last term is the inhibition on

the rise of glucose exerted by the secretion of insulin,

where EG0is the specific rate of glucose removal at

zero insulin, and SIrepresents the normal insulin sen-

sitivity [7].

In Equation (2), the same as that used in [7] and

by De Gaetano et al. [6], the first term on the right

accounts for the increase in insulin level due to the

increasing level of glucose, with σbeing the maximal

rate of secretion of insulin by β-cells, and αthe glu-

cose concentration yielding 50% of insulin secretion.

Here, νhis the power coefficient for the Hill-shaped

glycemia effect on pancreatic insulin release. This is

set to 4 in [6]. The last term here is the rate of re-

moval of insulin by natural means depending on how

much insulin is present, with kbeing the removal rate

constant of variation.

Equation (3) describes the rate of change of func-

tional β-cells, following [7], where r0is the death rate

of β-cells at zero glucose, r1the first order coefficient

for β-cells replication, and r2is the second order co-

efficient for β-cells replication [7].

Equation (4) has been modified from that men-

tioned in [7]. The rate of increase of inflammato-

ry response IRP is taken to vary, in the absence of

drug residuals, glucose or insulin, as a product of

inflammation factor due to SARS-CoV-2 infectivi-

ty, with kSARS representing the SARS-CoV-2 affin-

ity with ACE2, and the abnormal activity of ACE2,

∆kACE2,0, which we express here as a multiple ∆of

the normal value kACE2,0of kACE2. According to

[17], for someone who has diabetes, taking an ACE

inhibitor or ARB can help to control high blood pres-

sure which increases the risk of problems from dia-

betes. Diabetes can damage the blood vessels in the

kidneys, and high blood pressure can also damage the

kidneys. According to [7], ACE inhibitor or ARB

treatments could increase ACE2 abundance and thus

enhance viral entry. The second term on the right of

(4) therefore accounts for the increase in inflammato-

ry response due to drug that remains in the patient's

body, kDbeing the rate constant in the inflammation

rate due to ACE inhibitor surplus. Assuming contin-

ual treatments by ACE inhibitor, the quantity DSin

(4) is set to the mean value of 0.0372 mg/ml on day

29 of Sotrovimab treatment, as given in [18].

The third term in (4) is the inflammation rate due

WSEAS TRANSACTIONS on BIOLOGY and BIOMEDICINE

DOI: 10.37394/23208.2022.19.20

Chontita Rattanakul, Yongwimon Lenbury,

Nathnarong Khajohnsaksumeth, Charin Modchang

E-ISSN: 2224-2902

179

Volume 19, 2022

to glucose surplus, with kGbeing the rate constant for

this mechanism. Finally, the last term in (4) accounts

for the anti-inflammatory response, which varies di-

rectly as the inflammation response at time t, with

constant of variation keff .

Before we carry out the singular perturbation anal-

ysis, we need the following result concerning the in-

variance property of the model system. To that end,

we let

G=R0

EG0

, G =¯

EG0

,¯

I=σ

k,(6)

βf=¯

r2G2,¯

IRP =¯

keff

,(7)

k1=R0−SI

,¯

k2=−r0+r1¯

G, (8)

k3=kSARS ∆kACE2,0+kDDS+kG¯

G(9)

Lemma 2

Noting that all model parameters are positive, we

let Bbe the set in the G, I, βf, IRP -space defined as

B={G, I, βf, IRP ∈R4:G≤G≤¯

0< I ≤¯

I, 0< βf≤¯

βf,0< IRP ≤¯

IRP }

(10)

Then, Bis positive invariant and all solutions starting

in Bare uniformly bounded, provided

min{¯

k1,¯

k2}>0(11)

Proof

Suppose G0, I0, βf0, IRP 0∈ B. First, we show

that G, I, βf,and IRP remain positive.

By contradiction, if Gwere to be negative at

some point in B, then there has to be a t0>0where

G(t0) = 0 and G′(t0)<0. However, by (1), at this

point,

dt t=t0

=R0−EG0G(t0)−SI

G(t0)I(t0)

IRS (t0) + i

=R0>0

which is a contradiction, and thus Gnever vanishes.

Similarly, if Iwere to be negative at some point in

B, then there has to be a t0>0where I(t0) = 0 and

I′(t0)<0. However, by (2), at this point,

dt t=t0

=σβfGνh(t0)

α+Gνh−kI(t0)

=σβfGνh(t0)

α+Gνh>0(12)

since we have already shown that Gis always posi-

tive. This is a contradiction, and thus Inever vanish-

es.

Next, if βfwere to be negative at some point in

B, then there has to be a t0>0where βf(t0)=0

and β′

f(t0)<0. However, by (3), at this point

dβf

dt t=t0

=−r0+r1G(t0)−r2G2(t0)βf(t0)

=−r0+r1G(t0)

>−r0+r1G > 0

due to (11). This is a contradiction, and thus βfnever

vanishes.

Finally, if IRP were to be negative at some point

in B, then there has to be a t0>0where IRP (t0) = 0

and I′

RP (t0)<0. However, by (4), at this point,

dIRP

dt t=t0

=kSARS ∆kACE2,0+kDDS+kGG(t0)

−keff IRP (t0)

> kSARS ∆kACE2,0+kDDs>0

which, again, is a contradiction, and thus IRP never

vanishes.

Now, we consider the system (1)-(4)

withG0, I0, βf0, IRP 0∈ B, and first show

that Iis bounded. From (2), we have

dt =σβfGνh

α+Gνh−kI ≤σ¯

βf−kI

Then, by the comparison theorem, one has

I≤I0e−kt +σ

k1−e−kt≤¯

Ie−kt+¯

I1−e−kt=¯

using (6). Thus, Iremains bounded above by ¯

Then, from (1), we have

dt =R0−EG0G−SI

IRS +i≤R0−EG0G

Thus, by the comparison theorem, we have

G≤G0e−EG0t+R0

EG01−e−EG0t

≤¯

Ge−EG0t+¯

G1−e−EG0t=¯

using (6). Therefore, Gremains bounded above by

G. On the other hand, (1) gives

dt ≥R0−EG0G−SI

=¯

k1−EG0G

Thus, by the comparison theorem, we have

G≥G0e−EG0t+¯

EG01−e−EG0t

≥Ge−EG0t+G1−e−EG0t=G(13)

using (6) and (11). Therefore, Gremains bounded

WSEAS TRANSACTIONS on BIOLOGY and BIOMEDICINE

DOI: 10.37394/23208.2022.19.20

Chontita Rattanakul, Yongwimon Lenbury,

Nathnarong Khajohnsaksumeth, Charin Modchang

E-ISSN: 2224-2902

180

Volume 19, 2022

below by G. From (3),

dβf

dt ≤ −r0+r1¯

G−r2G2βf=¯

k2−r2G2βf

using (9). Then, by the comparison theorem, we have

βf≤βf0e−r2G2t+¯

r2G21−e−r2G2t

≤¯

βfe−r0r2

r1t+¯

βf1−e−r0r2

r1t=¯

βf(14)

using (7). Thus βfremains bounded above by ¯

βf.

Finally, (4) gives

dIRP

dt ≤kSARS ∆kACE2,0+kDDS+kG¯

G−keff IRP

≤¯

k3−keff IRP

using (9) and (11). Then, by the comparison theory,

we have

IRP ≤IRP 0e−kef f t+¯

keff 1−e−kef f t

≤¯

IRP e−kt +¯

IRP 1−e−kt=¯

IRP

using (7). Thus, IRP remains bounded above by ¯

IRP .

Hence, we have shown that Bis positive invari-

ant, and all solutions starting in Bremain uniformly

bounded. .

With the above result, we are in the position to an-

alyze the model to identify various dynamic behavior

permitted by the system. Substituting (5) into (1), we

introduce the following change of variables and pa-

rameters, to scale out their units, and to reflect, with

the use of small scaling factors εand δ, the obser-

vation mentioned earlier that Gappears to be varying

at a very fast speed. βfvaries at an extremely slow

speed, whereas Ivaries at an intermediate speed, rel-

atively, with time. As for IRP , we see in [7] that kef f

is noticeable larger than kDor kG, it is reasonable to

assume that IRP varies at more or less the same speed

as I. Letting

G′=G

G∗, I′=I

I∗, β′

f=βf

β∗

, I′

RS =IRS

I∗

,(15)

I′

RP =IRP

I∗

, R′

0=R0t∗

G∗, E′

G0=EG0t∗,(16)

S′

I=SII∗t∗

a0(G∗)2, b′

0=a0b0

G∗, µ =a0(c0+i)

(G∗)2,(17)

σ′=σt∗β∗

εI∗, α′=α

(G∗)νh, k′=kt∗

ε,(18)

r′

0=r0t∗

εδβ∗

, r′

1=r1G∗t∗

εδβ∗

, r′

2=r2(G∗)2t∗

εδ ,(19)

k′

ACE2=kACE2

k∗

ACE2

, k′

SARS =k∗

SARS k∗

ACE2t∗

I∗

,(20)

D′

S=DS

D∗

, k′

D=kDD∗

St∗

I∗

,(21)

k′

G=kGG∗t∗

I∗

, k′

eff =kef f t∗(22)

where (·)∗denotes the unit of the corresponding quan-

tity (·), and substituting the above new variables and

parameters into (1) - (4), then dropping the primes,

we arrive at the following model system, equivalent

to (1) - (4).

dt =R0−GEG0+SI

G2+b0G+µ

≡f(G, I)(23)

dt =εσβfGνh

α+Gνh−kI

≡g(G, I, βf)(24)

dβf

dt =εδ−r0+r1G−r2G2βf

≡h(G, βf)(25)

dIRP

dt =kSARS ∆kACE2+kDDS+kGG−kef f IRP

(26)

where now the variables are dimensionless, since we

have scaled out their corresponding units. Their or-

ders of magnitude depend on our choice of the units,

G∗, I∗, . . . and so on. However, their rates of change

are of different orders of magnitude, due to the small

sizes of εand δ, as well as the sizes of the parameters

appearing in the equations.

We see that, once Gis known by solving (23)-

(25), IRP can be found by solving (26), which does

not affect Equations (23)-(25) in any way, we can es-

sentially leave (26) out of the singular perturbation

arguments for the time being and concentrate on the

system of 3 equations (23)-(25) which are now in the

appropriate form on which geometric singular pertur-

bation analysis can be applied in the next section.

3 Model Analysis

3.1 Equilibrium Manifolds

Considering (23) - (25), we see that when f, g, and h

are finite and nonzero, 

dt 

=O(1),

dt 

=O(ε),

and 

dβf

dt 

=O(εδ), so that, for small εand δ,G,

WSEAS TRANSACTIONS on BIOLOGY and BIOMEDICINE

DOI: 10.37394/23208.2022.19.20

Chontita Rattanakul, Yongwimon Lenbury,

Nathnarong Khajohnsaksumeth, Charin Modchang

E-ISSN: 2224-2902

181

Volume 19, 2022

the scaled glucose concentration, has the fastest dy-

namics, while I, the scaled insulin concentration, has

intermediate dynamics, and βfthe scaled functional

β-cells density, varies at the slowest speed amongst

the three state variables. Then, for sufficiently small

values of εand δ, under suitable regularity conditions,

the singular perturbation theory allows approximating

the solution of the system (23) - (25) with a sequence

of simple dynamic transitions occurring at differen-

t speeds. More specifically, for sufficiently small ε

and δ, the solution of the model system will stay with-

in a θ-tube around the approximating transitions for

some small θ,θ > 0, which approaches 0 as εand

δtend to zero. For more detailed description of the

singular perturbation theory and the geometric pertur-

bation technique, please see [19, 20]. This method of

analysis allows us to predict how the solution trajec-

tories will appear when the model's state variables are

changing at radically heterogenous time scales. With-

out using this method, many dynamic behavior may

be hidden and cannot be identified. In addition, the

delineating conditions that allow such dynamics may

not possibly be derived explicitly. If the state vari-

ables were treated as having the same dynamics, then

we may miss the solution trajectories that could be ob-

served only if one variable were recognized to have a

very slow dynamics in comparison with the other s-

tate variables. Specifically, supported by the values

of the physiological parameters reported in the liter-

ature [6] and [7], the functional beta cells have been

found to possess a very slow dynamics since the size

of its derivative is of order εδ, while glucose has a fast

dynamics since its derivative is O(1). Insulin, on the

other hand, has the intermediate dynamics, its rate of

change is of the order ε, where εand δare small and

<1.

To analyze the model, we first identify and de-

scribe the 3 equilibrium manifolds of interest.

The slow equilibrium manifold

The slow equilibrium manifold is obtained by equat-

ing the right hand side of (25) to zero. Solving it for

βfin terms of G, we arrive at the following equation.

We observe that, on this manifold, βf= 0 when

G=G∗≡r0

(27)

and βf<0if Gbecomes smaller than G∗. Moreover,

βf→0as G→ ∞ , and βf→ −∞ , as G→0. We

can easily find that, on this surface,

dβf

dG =−r1G+ 2r0

r2G

which means the slope of the graph will become 0

when

G=2r0

≡ˆ

G= 2G∗

βf=r1G−r0

r2G2(28)

Thus, the graph of this surface, for appropriate

parametric values, may be depicted as in Fig.1a)

Figure 1: The equilibrium manifolds: a) the slow

manifold, b) the intermediate manifold, c) the fast

manifold.

WSEAS TRANSACTIONS on BIOLOGY and BIOMEDICINE

DOI: 10.37394/23208.2022.19.20

Chontita Rattanakul, Yongwimon Lenbury,

Nathnarong Khajohnsaksumeth, Charin Modchang

E-ISSN: 2224-2902

182

Volume 19, 2022

By considering (25), we know that when we are

in the front of the surface where βfis large, here

dβf

dt <0, then the transition will take us in the

direction of decreasing βfuntil the surface dβf

dt = 0

is reached. On the other hand, if we start at the

point behind this surface where βfis small, then

dβf

dt >0there and the transition will take us in the

direction of increasing βfuntil the surface dβf

dt = 0

is reached. These transitions will be very slow since



dβf

dt =O(εδ).

The intermediate equilibrium manifold

The intermediate equilibrium manifold is obtained

by equating the right hand side of (24) to zero. Solv-

ing the resulting equation for Iin terms of G, we ar-

rive at the following equation.

I=σβfGνh

kα +kGνh(29)

We observe that, on this manifold, I→σβf

kas

G→ ∞ , and I→0as G→0when βfis constant.

On the other hand, when Gis kept constant,I→ ∞

as βf→ ∞ , and I→0as βf→0. The graph

of this surface, for appropriate parametric values,

may be depicted as in Fig.1b) Looking at (24), we

see that when we are above this surface where Iis

large, dI

dt <0here, and then the transition will take

us in the direction of decreasing Iuntil the surface

is reached where dI

dt = 0. On the other hand, if we

start at the point below this surface where Iis small,

then dI

dt >0there and the transition will take us in

the direction of increasing Iuntil the surface dI

dt = 0

is reached at which point the transition stops. These

transitions will be made at an intermediate speed

since 

dt 

=O(ε).

The fast equilibrium manifold

The fast equilibrium manifold is obtained by e-

quating the right hand of (23) to zero. Solving the

resulting equation for Iin terms of G, we arrive at

the following equation.

I=R0−EG0G)(a0G2+b0G+µ)

SIG(30)

We observe that, on this manifold, I→ −∞ as

G→ ∞ ,I→ ∞ as G→0, and I= 0 at G=R0

EG0

Moreover, since (30) does not involve βf, this surface

will be parallel to the βf- axis.

Differentiating Iwith respect to Ggives

dG =a0R0−EG0b0G2−2a0EG0G3−R0µ

SIG2

(31)

Letting

B=−a0R0−EG0b0

2a0EG0

and D=R0µ

2a0EG0

(32)

we can thus prove the following result.

Lemma 2 The curve resulting from the intersection

between the fast manifold ˙

G= 0 and the slow mani-

fold ˙

βf= 0 has 2 relative extrema at G=GM1>0

and G=GM2> GM1>0provided (11) holds,

a0R0−EG0b0>0(33)

and

D < −4

27B3(34)

Proof

From the derivative found in (31), we get that the

relative extrema will occur on the curve if the numer-

ator of (31) vanishes, that is

(a0R0−EG0b0)G2−2a0EG0G3−R0µ= 0 (35)

By (33), we have B < 0and D > 0, while the

discriminant of (35) is

∂=−4B3D−27D2

which is positive since (34) holds. Therefore, we

can conclude that (35) has 2 positive solutions by the

Descartes' Rule of Signs. This means that there exist

2 relative extrema, G=GM1and G=GM2, on this

curve with GM1,2positive. 

Let the two relative extrema on this curve on the

fast manifold be denoted by (GM1, IM1, βfM1)and

(GM2, IM2, βfM 2)where

IM1,2=(R0−EG0GM1,2)(a0G2

M1,2+b0GM1,2+µ)

SIGM1,2

(36)

We note that, by (35), GM1,2satisfy the equation

M1,2−R0

2EG0

M1,2+R0µ

2a0EG0

= 0

Thus,

R0−EG0GM1,2=EG0GM1,2+R0µ

a0G2

M1,2

WSEAS TRANSACTIONS on BIOLOGY and BIOMEDICINE

DOI: 10.37394/23208.2022.19.20

Chontita Rattanakul, Yongwimon Lenbury,

Nathnarong Khajohnsaksumeth, Charin Modchang

E-ISSN: 2224-2902

183

Volume 19, 2022

By (36), we then see that IM1,2>0, which means

both relative extrema are located in the first octant of

the (G, I, βf, IRP )-space in which all state variables

are positive. The graph of this surface, for appropri-

ate parametric values, may be depicted as in Fig. 1c).

If we start at a point to the above and to the right of

the fast manifold, according to (23), dG

dt <0here,

then the transition will take us in the direction of de-

creasing G. On the other hand, starting from the area

below and to the left of the fast manifold, we will be

in the region where dG

dt >0and hence the transition-

s will be in the direction of increasing G. Thus, as

indicated in Fig.1c), we see that the fast manifold has

two stable branches, namely the portion of the surface

where G<GM1and the one where GM2< G . It

has one unstable branch where GM1< G < GM2.

3.2 Geometric Singular Perturbation

Analysis

We now let I1,2be the values of Ion the intersection

curve between the intermediate and slow manifolds at

the points where G=GM1,2, that is

I1,2=σβf1GM1,2GM1,2νh

kα +kGM1,2νh(37)

with

βf1=r1GM1−r0

r2G2

(38)

and also let S(GS, IS, βf,S )be the critical point

where the 3 manifolds intersect, under appropriate

parametric values, where

f(GS, IS) = g(GS, IS, βf,S ) = h(GS, βf,S )

The following theorem can then be stated.

Theorem 1 The model system (23)-(25) admits a peri-

odic solution in the form of a limit cycle surrounding

the equilibrium point S(GS, IS, βf,S ), at which the

3 manifolds intersect, provided (11), (33), (34) hold,

and

GM0> G (39)

I1< IM1, IM2< I2(40)

where G=GM0is the smallest root of the following

equation:

IM2=(R0−EG0G)(a0G2+b0G+µ)

SIG

Proof

The inequality (40) insures that the graph of the

fast manifold where f= 0, for all βf, will

have, as depicted in Fig.1c), two relative extrema

GM12, IM12, βf M 2. The inequality (39), consider-

ing the definition of G, insures that the value of βf

remains positive when Ireaches the height IM2.

The 3 equilibrium manifolds in this case are shown

together in Fig.2a), where the fast manifold f= 0 is

shown in red, the intermediate manifold is in green

and the slow manifold is in blue. The transitiontions

are traced in purple, where the fast transitions are in-

dicated by triple arrows, double arrows indicate ones

with intermediate speed, and a single arrow indicates

a slow transition.

We see clearly in this figure that if inequality (40)

holds, then the green curve where the fast and inter-

mediate manifolds intersect, along which f=g= 0,

lies in the region where the manifold f= 0 is unstable

for all Gin an open interval ΩS=GS−ω, GS+ω

for some ω > 0.I1< IM1means the intermedi-

ate manifold lies below the fast manifold for Gin the

interval GS−ω, GS, whereas IM2< I2means

the intermediate manifold lies above the fast manifold

for Gin the interval GS, GS+ω. Thus, the curve

f=g= 0 lies in the region where the fast manifold

is stable for Gin ΩS.

Starting from a generic initial point marked with

Ain Fig.2a), AG0, I0, βf,0∈ B , such that G0∈

Ωs, G0> GS, βf> βf,S , since dG

dt <0here,

a fast transition will develop, in the direction of de-

creasing G, toward the point Bon the nearest stable

branch of the manifold f= 0, while Iand βfremain

frozen. The transition is made at a fast speed relative-

ly because 

dt 

=O(1) .

Once Bis reached, the intermediate system has

become active. Here, we are below the intermediate

manifold so that dI

dt >0and a transition is made, in

the direction of increasing Iat an intermediate speed,

since 

dt 

=O(ε), until the point Cis reached,

where the stability is lost, and a quick jump brings

the system to the other stable branch of the manifold

f= 0 at point D. Here, dI

dt >0and so a mo-

tion at intermediate speed develops down the surface

f= 0 until the point Eis reached where the stability

of the manifold is again lost and a quick jump brings

the trajectory to the point Fon other stable branch.

Next, for the same reason as before, a transition of

intermediate speed develops on this manifold in the

direction of increasing Iagain to Hjust missing the

point Csince βfhas been slowly decreasing. The

cycling motion repeats while βfdecreases slowly, as

WSEAS TRANSACTIONS on BIOLOGY and BIOMEDICINE

DOI: 10.37394/23208.2022.19.20

Chontita Rattanakul, Yongwimon Lenbury,

Nathnarong Khajohnsaksumeth, Charin Modchang

E-ISSN: 2224-2902

184

Volume 19, 2022



dβf

dt 

=O(ε, δ), until the manifold f=h= 0

is reached at a point J. A quick jump to Kis fol-

lowed by a upward transition to Lwhere the stability

is lost. A quick jump is then made to Mfollowed by a

downward transition at intermediate speed back to J

thereby completing a stable limit cycle JKLM J sur-

rounding the unstable steady state SGS, IS, βf,S .

Of equal interest, is the case in which the system

is stable and the state variables all tend towards their

respective steady state values. By similar arguments,

we can prove the following result.

Theorem 2 The model system (23)-(25) admits an

equilibrium solution, GS, IS, βf,S , which is local-

ly asymptotically stable provided (11), (32), (34) and

(39) hold, and the inequalities in either Case 1 or Case

2 below are satisfied.

Case 1: I1> IM1,and I2> IM2(41)

Case 2: I1< IM1,and I2< IM2(42)

Proof

The arguments for Case 1 or Case 2 are essentially

the same. We will therefore only show Case 2 here. In

this case, the shapes and positions of the 3 equilibrium

manifolds are as depicted in Fig.2b), for appropriate

parametric values. The inequalities in (42) insure that

a portion of the curve f=g= 0 intersects the curve

f=h= 0 at the point SGS, IS, βf,S that is located

on a stable portion of the fast manifold.

Again, Starting from a generic initial point marked

with Ain Fig.2b), AG0, I0, βf,0∈ B , such that

G0∈Ωs, where, say G0< GS, βf> βf,S since

dt >0here, a fast transition will develop, in the di-

rection of increasing G, toward the point Bon the n-

earest stable branch of the fast manifold f= 0, while,

for the time being, Iand βfremain frozen.

Once Bis reached, the intermediate system has

become active. Here, dI

dt <0and a transition is made

at intermediate speed in the direction of decreasing I,

until the point Cis reached, where the stability is lost.

As before, a quick jump brings the system to the other

stable branch of the manifold f= 0 at point D.

Here,dI

dt >0and so a motion at intermediate speed

develops upward along the surface f= 0 until

we arrive at the point Eon the curve where the

intermediate manifold intersects the fast manifold

and the intermediate dynamic now becomes active.

Since this is still on the stable portion of the fast

manifold, stability is maintained and a transition is

then made slowly along this curve in the direction of

increasing βf, since dβf

dt >0here. Eventually, the

point where the curve f=g= 0 intersects the slow

manifold so that the 3 manifolds intersect at the point

GS, IS, βf,sis reached, marked by Sin Fig.2b),

where transitions stop since dG

dt =dI

dt =dβf

dt = 0

here. Thus, this is the case where the steady state is

locally asymptotically stable. 

Figure 2: The fast (red), intermediate (green) and

slow (blue) manifolds are here shown in the 3-

dimensional space with transitions of fast speed

(triple arrows), intermediate speed (double arrows)

and slow speed (single arrow) in the case that a) the

model system admits a limit cycle surrounding the

steady state S, described in Theorem 1, and b) the

trajectory tends to a stable steady state in case 2 of

Theorem 2

In the next section, we show numerical simulation-

s of the model under investigation to support our the-

oretical predictions.

4 Model Simulation

We present here numerical simulations of the model

system (23) - (26) in various cases mentioned above.

The values of those parameters not involved in the sat-

isfaction of conditions that delineate different dynam-

ic behaviors are estimated or based on those utilized

in the literatures, mainly in [6],[7], and [15]. The

WSEAS TRANSACTIONS on BIOLOGY and BIOMEDICINE

DOI: 10.37394/23208.2022.19.20

Chontita Rattanakul, Yongwimon Lenbury,

Nathnarong Khajohnsaksumeth, Charin Modchang

E-ISSN: 2224-2902

185

Volume 19, 2022

units of the variables have already been scaled out by

the introduction of new variables and parameters in

(15) - (22). Thus, what units the parameters have de-

pends on the units with which you choose to measure

the state variables, ensuring that both sides of Equa-

tions (23) - (26) have the same units. To specify phys-

iological values for the model parameters, we rely as

much as possible on [6] and [7], where rigorous re-

view of the literature was carried out to identify phys-

iological values for the relevant parameters in their

model. Since information for some parameters is lim-

ited, most were based on empirical observations and

their values were derived using well established sta-

tistical methodologies [7]. Some of the parameters

are here varied to satisfy the inequalities required for

each of the cases to hold. The variations reflect vari-

ous symptoms of the patients which will be explained

accordingly.

It must be mentioned here that, according to Vice-

conti et al. [21], the utilization of integrative ap-

proaches to biomedical research creates the need to

manage large amounts of data that are radically het-

erogeneous in structure, complicated by the fact that

more often than not these data are measured or ob-

served at very different dimensional and/or temporal

scales. Their study in [21] describes a first attempt at

providing a strategical environment for homogeneous

biomedical data defined over radically different spa-

tial or temporal scales. New strategies for the man-

agement are currently discovered and developed for

the management of temporal multi-scaling or highly

heterogenous data types. It is, however, beyond the

scope of our research.

The model system (23) - (26) is simulated with

parametric values satisfying inequalities identified in

Theorem 1, in which case the model system is predict-

ed to admit limit cycles towards which the solution

trajectory tends as time progresses, as seen in Fig.3.

The corresponding time series of the 4 state variables

are shown in Fig.4, exhibiting sustained oscillation

mimicking clinical data reported in several literatures.

Here, G0= 2, I0= 2, βf0= 3, IRP 0= 1,

t0= 0, R0= 1.12, SI= 0.889, EG0= 0.054,

a0= 1, b0= 0.2347, c0= 0.000313, i = 0.000187,

σ= 0.20592, α = 0.0017, k = 6.5, r0= 0.00001,

r1= 0.0884, r2= 0.024, kD= 0.001,

kSARS = 0.15, kACE2= 0.4, kACE2,0= 0.385,

DS= 0.0372, kG= 0.1, keff = 0.689,∆ = 1,

ε= 0.001, δ = 0.8.

The inflammatory response IRP also oscillates as

time passes. This reflects a control system which is

working comparatively efficiently, in which an in-

crease in the glucose level triggers the beta cells to be-

come functional and secrete insulin in response. The

rise in insulin is then followed by the satisfactory drop

in the glucose level. It is then no longer necessary for

the beta cells to produce the hormone insulin, the lev-

el of which then drops with a little time lag following

the glucose drop.

The numerical solution of (23) - (26) in Case 2,

where inequalities identified in Case 2 of Theorem 2

are satisfied, is shown in Fig.5 and 6. Case 1 is simi-

lar, so that the simulation has been omitted.

Figure 3: Numerical solution of (23) - (26) in the case

that the inequalities in Theorem 1 are satisfied. The

solution trajectory is seen here to tend towards the sta-

ble limit cycle as time passes as theoretically predict-

ed.

Here, G0= 0.1, I0= 0.1, βf0= 0.1,

IRP 0= 0.2, t0= 0, R0= 1.6, EG0= 0.8,

SI= 0.04998, a0= 1, b0= 0.00001,

c0= 0.000413, i = 0.000087, σ = 0.00352, α = 2,

k= 11, r0= 0.001, r1= 0.0884, r2= 0.074,

kSARS = 0.15, kACE2= 0.4, kACE2,0= 0.385,

kD= 0.001, DS= 0.0372, kG= 0.1,

keff = 0.689,∆ = 1, ε = 0.31, δ = 0.28.

We can observe clearly here that there is a definite

delay after the level of glucose in the bloodstream

increases before insulin is secreted in response. The

difference from the previous case where sustained

oscillation is exhibited is that the value of R0is larger,

WSEAS TRANSACTIONS on BIOLOGY and BIOMEDICINE

DOI: 10.37394/23208.2022.19.20

Chontita Rattanakul, Yongwimon Lenbury,

Nathnarong Khajohnsaksumeth, Charin Modchang

E-ISSN: 2224-2902

186

Volume 19, 2022

Figure 4: The time courses of a) glucose concentra-

tion, b) insulin concentration, c) functional beta cells,

and d) the inflammatory response, in the case shown

in Fig.3. All 4 state variables are seen to oscillate pe-

riodically as time progresses.

which drives glucose to increase more quickly, while

a smaller SIhas the effect of slowing down the

response of insulin. Moreover, σis very much smaller

which means the increase in Gis felt to a lesser degree

by the rate of insulin production, so that it can manage

to bring down the glucose level only a little before it

settles to a rather high steady state level. More impor-

tantly, we see that the inflammatory response is able

to rise up to a sustainably high level, which is not de-

sirable in the clinical point of view of a patient with

comorbidity with COVID-19.

Figure 5: Numerical solution of (23) - (26) in the

case that the inequalities in Theorem 2 are satisfied.

The solution trajectory is seen here to tend towards

the locally asymptotically stable steady state as time

passes as theoretically predicted.

5 Discussion and Interpretation

Considering the conditions (34) and (39) in both

Theorem 1 and Theorem 2 which ensure that the

solutions remain bounded and local stability is main-

tained, we see that we must have

WSEAS TRANSACTIONS on BIOLOGY and BIOMEDICINE

DOI: 10.37394/23208.2022.19.20

Chontita Rattanakul, Yongwimon Lenbury,

Nathnarong Khajohnsaksumeth, Charin Modchang

E-ISSN: 2224-2902

187

Volume 19, 2022

Figure 6: The time courses of a) glucose concentra-

tion, b) insulin concentration, c) functional beta cells,

d) and the inflammatory response, in the case shown

in Fig.5. All 4 state variables are seen to tend toward

their respective steady state values as time progresses.

µ < a0R2

27E2

(43)

for (34) to hold, and

r0< r1GM0(44)

for (39) to hold. Inequality (43) will be satisfied if

µ, the sum of c0and i, should be sufficiently small

so that it is smaller than the quantity on the right of

the inequality (43). This means that the zero glucose

coefficient in the numerator of the last term in (23),

reflecting insulin resistance, is sufficiently small.

Essentially, insulin resistance should be sufficiently

small when negligible glucose is present.

In addition, r0, the death rate of β-cells at zero glu-

cose, must be small in comparison to r1, the first order

coefficient for β-cells replication. This means β-cells

should be slow in dying and, at the same time, they

should be able to replicate quickly even with low glu-

cose activation. If the above conditions are not met,

the system should be recognized as defective. Abnor-

mal dynamics could be expected if these 2 conditions

are violated. Moreover, considering the position of

the slow manifold given by (28) relative to the fast

and intermediate manifolds in the phase space, we can

see that trajectories, traced by the singular perturba-

tion arguments, will develop in a way that is indicative

of symptoms of serious insulin resistance if

βfg=0 =Iˆ

Gkα +kˆ

Gνh

σˆ

Gνh

>r1ˆ

G−r0

r2ˆ

G2=βfh=0 (45)

which reflects a disfunctional control system since it

appears not to be able to produce insulin, or activate

the β-cells to do so to decrease the glucose level in the

bloodstream. Condition (45) can be achieved if kis

sufficiently large or σis sufficiently small, in which

case k≫1or σ≪1insulin is not secreted sufficient-

ly in response to high glucose and consequently, in-

sulin level continues to drop necessitating exogenous

insulin supplements to bring the system under con-

trol. We show a simulated time series of the model's

solution in this scenario in Fig.7 with G0= 0.8, I0=

0.5, βf0= 0.01, IRP O = 2, t0= 0, R0= 2, EG0=

0.06, SI= 0.01, a0= 1, b0= 0.00001, c0=

0.000313, i = 0.000187, σ = 0.009, α = 0.07, k =

8, r0= 0.007, r1= 0.018, r2= 0.074, kSARS =

0.15, kACE2= 0.4, kACE2,0= 0.385, kD=

0.001, Ds= 0.0372, kG= 0.1, kff = 0.689, δ =

1, ε = 0.21, δ = 0.28.

Considering (26), the presence of high ACE2 in-

hibitor can contribute to high rate of inflammatory

WSEAS TRANSACTIONS on BIOLOGY and BIOMEDICINE

DOI: 10.37394/23208.2022.19.20

Chontita Rattanakul, Yongwimon Lenbury,

Nathnarong Khajohnsaksumeth, Charin Modchang

E-ISSN: 2224-2902

188

Volume 19, 2022

response increase. ACE inhibitors are a class of drug

known to cause inflammation and allergic reactions

[23]. Fine tuning of the drug therapy may be a good s-

trategy to compensate for worrisome constant inflam-

matory activation, in the case where the rise of glu-

cose level is unchecked and insulin level is persistent-

ly low.

Although inflammation is beneficial, since it pro-

tects against infection and injury, when that response

is constantly triggered, it can be harmful the body in-

stead of healing it in the long run. As reported in

[24], the prolonged inflammatory activation and sub-

sequent response commonly result in increased sus-

ceptibility to coronavirus infection leading to more

severe respiratory failure and end-organ damage in

patients with SARS-CoV-2.

Figure 7: The time courses of (a) glucose concen-

tration, (b) insulin concentration, (c) functional beta

cells, and (d) the inflammatory response, in the sce-

nario in which the glucose level quickly rises while

insulin drops and the inflammatory response rises to

a high level.

6 Conclusion

A model of insulin-glucose control system, incorpo-

rating inflammatory response, which, when triggered

constantly, over time, can harm the body instead of

healing it, has been studied. The inflammatory acti-

vation and subsequent response, may more common-

ly than not lead to increased susceptibility to coron-

avirus infection. Theoretical analysis by the geomet-

rical singular perturbation technique, well suited for

such systems subject to highly diversified time scales,

has been carried out, and numerical simulations dis-

cussed in terms of the delineating conditions on the

system's parameters. The state variables tracked by

our model, the glucose and insulin concentrations, the

functional βcells and the inflammatory response will

eventually tend asymptotically to their corresponding

steady state values, or oscillate in a bounded fashion,

provided that the zero glucose coefficient in the nu-

merator of the last term which represents insulin re-

sistance, in the rate equation for glucose level, is suffi-

ciently small. Essentially, insulin resistance should be

sufficiently small when negligible glucose is present.

The conditions that are indicative of a dysfunction-

al system, appearing not to be able to produce insulin,

or activate the β-cells to do so to decrease the glucose

level in the bloodstream, are manifested if kis too

large or σis too small, in which case k≫1or σ≪1.

Insulin is not secreted sufficiently in response to high

glucose and consequently, insulin level continues to

drop necessitating exogenous insulin supplements to

bring the system under control.

Further research could involve investigating the

impacts of food intakes and exercises, as well as tem-

poral drug administrations and periodic hormone sup-

plements, which can be incorporated in to the model

derived here with relative ease. Such studies can lead

to additional insights and valuable suggestions of ef-

ficient treatment protocols for patients who could be

exhibiting various symptoms of comorbidity.

Acknowledgment. This research has received

funding support from the NSRF via the Program Man-

agement Unit for Human Resource & Institutional

Development, Research and Innovation (grant num-

ber B05F640231). Also, this research is partially sup-

ported by the Centre of Excellence in Mathematics,

Ministry of Higher Education, Science, Research and

Innovation, Thailand(grant number RG-01-65-01-1).

References:

[1] S. Babashov, Predicting the Dynamics of Covid-

19 Propagation in Azerbaijan based on Time

Series Models,WSEAS TRANSACTIONS on

ENVIRONMENT and DEVELOPMENT, DOI:

10.37394/232015.2022.18.99

[2] G. Makanda, A mathematical model for the pre-

diction of the impact of coronavirus (COVID-

19) and social distancing effect, WSEAS TRANS-

ACTIONS on SYSTEMS and CONTROL, DOI:

10.37394/23203.2020.15.60

[3] Y. Riyani, S. Andriana, K. Mardiah, L. Suher-

ma, B. Riyadhi, A. Arianto, K. Khamimi,

J. Jakfar, E. Endri, Stock Market Reaction-

s before and during the COVID-19 Pandemic:

Evidence from Indonesia, WSEAS TRANSAC-

TIONS on BUSINESS and ECONOMICS DOI:

10.37394/23207.2022.19.104

WSEAS TRANSACTIONS on BIOLOGY and BIOMEDICINE

DOI: 10.37394/23208.2022.19.20

Chontita Rattanakul, Yongwimon Lenbury,

Nathnarong Khajohnsaksumeth, Charin Modchang

E-ISSN: 2224-2902

189

Volume 19, 2022

[4] C. Smith, A dangerous new coronavirus com-

plication was discovered - and it never goes

away if you get it.(accessed: 06.03.2022),

https://bgr.com/science/coronavirus-symptoms-

complications-diabetes-onset-after-covid-19/.

[5] D. Campbell, Covidl19: people with type

1 diabetes more likely to die than those

with type 2 l study, Diabetes (accessed:

06.03.2022) ,https://www.theguardian.com/

society/2020/may/20/type-1-diabetics-type-2-

coronavirus-nhs-study

[6] A. De Gaetano, T. Hardy, B. Beck, E. Abu-

Raddad, P. Palumbo,J.B. Valleskey, N. Pϕrksen,

Mathematical models of diabetes progres-

sion. Am J Physiol Endocrinol Metab. 2008

Dec;295(6):E1462-79. doi: 10.1152/ajpen-

do.90444.2008.

[7] P. Barbiero , P. Li ó,The Computation-

al Patient has Diabetes and A COVID,

medRxiv (accessed: 06.03.2022),http-

s://doi.org/10.1101/2020.06.10.20127183.

[8] Y. Lenbury, S. Ruktamatakul, S. Amorn-

samarnkul, Modeling insulin kinetics: responses

to a single oral glucose administration or

ambulatory-fed conditions. Biosystems. 2001

Jan;59(1):15-25.

[9] D.V. Giang, Y. Lenbury, A. De Gaetano, Delay

model of glucose-insulin systems: Global stabil-

ity and oscillated solutions conditional on delays,

Journal of Mathematical Analysis and Applica-

tions, Volume 343, Issue 2, 2008, Pages 996-

1006.

[10] P. Palumbo and S. Panunzi and A. De Gae-

tano, Qualitative behavior of a family of delay-

differential models of the Glucose-Insulin sys-

tem. Discrete and Continuous Dynamical System-

s. Series B. No.7, 2007, 399-424.

[11] L. Kardar , A. Fallah , S. Gharibzadeh , F. Moz-

tarzadeh, Application of fuzzy logic controller

for intensive insulin therapy in type 1 diabetic

mellitus patients by subcutaneous route. WSEAS

Transactions on Systems and Control, Vol 3, No.

9, 2008, 712-721.

[12] M. Chuedoung , W. Sarika , Y. Lenbury,

Dynamical analysis of a nonlinear model for

glucose-insulin system incorporating delays and

b-cells compartment. Nonlinear Analysis - theory

Methods & Applications - NONLINEAR ANAL-

THEOR METH APP. Vol.71, No. 12, 2009,

e1048-e1058.

[13] X. Cao, D. Liu, S. Yu Wang.,Mathematical Mod-

elling and Stability Analysis for Diabetes Predict-

ing System ,WSEAS Transaction on Mathematics,

Vol.14, 178-191, 2015.

[14] A. Felman, What to know about

insulin resistance. MedicalNewsTo-

day, medically reviewed by Deborah

Weatherspoon on March 26, 2019,http-

s://www.medicalnewstoday.com/articles/305567.

php.

[15] C. Ratanakul, Y. Lenbury, J. Suksamran, Sin-

gular Perturbation Analysis for Identification of

Dynamic Behaviour and Stability of a Nonlin-

ear Model of Long Term Progression of Diabetes

Mellitus. WSEAS Transactions on Mathematic-

s,2020, 523-530.

[16] Mechanism Of Action Of Insulin

On Blood Glucose Level, Diabetes

Talk,https://diabetestalk.net/insulin/mechanism-

of-action-of-insulin-on-blood-glucose-

evel♯:∼:text=Insulin%20has%20a%20number

%20of%20actions%20on%20the,the%20cells

%20by%20activating %20the%20sodium-

potassium%20cellular%20channels

[17] K. Romito, A. Husney, T. O'Young, D.

C.W. Lau, Learning About ACE Inhibitors

and ARBs for Diabetes, available from:

https://healthy.kaiserpermanente.org/health-

wellness/health-encyclopedia/he.learning-about-

ace-inhibitors-and-arbs-for-iabetes.abq4936

♯:∼:text=When%20you%20have%20diabetes

%2C%20taking%20an%20ACE%20inhibitor,

High%20blood%20pressure%20can%20damage

%20the%20kidneys%2C%20too.

[18] Art 5(3) - Sotrovimab (VIR-7831 -

GSK4182136) for the treatment of COVID-19

(GSK) - assessment report, European Medicines

Agency, 20 May 2021

[19] S. Muratori, S. Rinaldi.,Low- and high- frequen-

cy oscillations in three-dimensional food chain

systems. ,SIAM Journal of Applied Mathematics,

Vol.52,1992,1688-1706.

[20] S. Rinaldi, M. Scheffer, Geometric Analysis of

Ecological Models with Slow and Fast Processes.

Ecosystems 3, 2000, 507�521.

[21] M. Viceconti, G. Clapworthy, D.Testi,

F.Taddei, N.McFarlane, Multimodal fusion

of biomedical data at different temporal

and dimensional scales, Computer Method-

s and Programs in Biomedicine, Vol. 102,

WSEAS TRANSACTIONS on BIOLOGY and BIOMEDICINE

DOI: 10.37394/23208.2022.19.20

Chontita Rattanakul, Yongwimon Lenbury,

Nathnarong Khajohnsaksumeth, Charin Modchang

E-ISSN: 2224-2902

190

Volume 19, 2022

Issue 3, June 2011, Pages 227-237, http-

s://doi.org/10.1016/j.cmpb.2010.04.017)

[22] K. Hobson,Is Inflammation Bad for Y-

ou or Good for You? Health News from

NPR,https://www.npr.org/sections/health-

shots/2017/07/21/538377221/is-inflammation-

bad-for-you-or-good-for-you,2017.

[23] T. Eyre and V. Van-Hamel-Parsons and L. Mun

Wang and K. A. Hughes and Timothy J Little-

wood, An unusual cause of anaemia of chronic

disease: lisinopril-induced chronic inflammatory

state.PubMed, 2011, doi: 10.1155/2011/939080.

[24] R. M. Inciardi, S. D. Solomon, P. M. Ridker, M.

Metra, Coronavirus 2019 Disease (COVID-19),

Systemic Inflammation, and Cardiovascular Dis-

ease. Journal of the American Heart Association,

2020.

Contribution of individual authors to

the creation of a scientific article

(ghostwriting policy)

Chontita Rattanakul carried out parts of the analysis,

checked all proofs and carried out model simulations

and figures production.

Yongwimon Lenbury constructed the model equa-

tions in the form that singular perturbation analysis

could be consequently carried out, and provided the

drawings and arguments concerning the geometric

singular perturbation.

Nathnarong Khajohnsaksumeth created the computer

codes for the simulations and production of the

manuscript.

Charin Modchang provided the initial fundamental

understandings of the problem of interest, specifical-

ly the mechanisms related to comorbidity between

diabetes and COVID-19, which assisted in the model

formulation.

Sources of funding for research

presented in a scientific article or

scientific article itself

This research has received funding support from the

NSRF via the Program Management Unit for Human

Resource & Institutional Development, Research and

Innovation (grant number B05F640231). Also, this

research is partially supported by the Centre of Excel-

lence in Mathematics, Ministry of Higher Education,

Science, Research and Innovation, Thailand (grant

number RG-01-65-01-1).

Creative Commons Attribution

License 4.0 (Attribution 4.0

International , CC BY 4.0)

This article is published under the terms of the

Creative Commons Attribution License 4.0

WSEAS TRANSACTIONS on BIOLOGY and BIOMEDICINE

DOI: 10.37394/23208.2022.19.20

Chontita Rattanakul, Yongwimon Lenbury,

Nathnarong Khajohnsaksumeth, Charin Modchang

E-ISSN: 2224-2902

191

Volume 19, 2022