Influence of incubation delays on COVID-19 transmission in diabetic and non-diabetic populations – an endemic prevalence case

Monalisa Anand; Palla Danumjaya; Ponnada Raja Sekhara Rao

doi:10.1515/cmb-2023-0115

Open Access Published by De Gruyter Open Access December 31, 2023

Influence of incubation delays on COVID-19 transmission in diabetic and non-diabetic populations – an endemic prevalence case

Monalisa Anand , Palla Danumjaya and Ponnada Raja Sekhara Rao

From the journal Computational and Mathematical Biophysics

https://doi.org/10.1515/cmb-2023-0115

Abstract

The study of dynamics of diabetic population infected by COVID-19 is of pressing concern as people with diabetes are considered to be at higher risk of severe illness from COVID-19. A three-compartment mathematical model to describe the interactions of diabetic population and non-diabetic population both infected by COVID-19 with a susceptible population is considered. Time delays in incubation periods of COVID-19 in diabetic and non-diabetic populations are introduced. Besides the basic properties of such a dynamical system, both local and global stability of endemic equilibrium, are studied. The lengths of time delays are estimated for which the stability of the system is preserved locally, while sufficient conditions on system parameters are obtained for global stability. Numerical examples are provided to establish the theory, and simulations are provided to visualize the examples. It is noted that an increase in length of time delay in either of infected populations leads to oscillations in susceptible population but has no impact on infected populations.

Keywords: COVID-19; diabetes; delay differential equation; equilibria; stability analysis; numerical experiments

MSC 2010: 92C42; 93C43; 92D25

1 Introduction

COVID-19 is a highly infectious respiratory illness caused by the novel coronavirus, SARS-CoV-2, and it was first identified in Wuhan, China, in December 2019. The transmission of this disease is via contaminated air droplets or small airborne particles containing the virus. General symptoms include fever, cough, headache, fatigue, and loss of taste and smell, among many others. Individuals of all ages are at risk of contracting this infection; however, patients with age greater than 60 years with underlying medical condition have an increased risk of developing serious complications [6]. The most prevalent comorbidity of COVID-19 is hypertension (42.46%), followed by diabetes (39.72%), old k-chest (20.54%), bronchial asthma (16.43%), coronary artery disease (13.69%), chronic kidney disease (13.69%), and valvular heart disease (6.84%) [4]. Studies [5,8,11,29] have highlighted that COVID-19 patients with diabetes, cardiovascular diseases, hypertension, and other comorbidities are likely to develop a more severe course and progression of the disease, which could develop into a life threatening situation leading to death.

COVID-19 patients with diabetes experience significantly higher rates of complications and mortality [16,19,28]. Joshi and Pozzilli [16] emphasized on the COVID-19 induced diabetes, which is a novel phenomena observed in critically ill patients. COVID-19 patients with diabetes frequently experience uncontrolled hyperglycaemia and episodes of acute hyperglycaemic crisis, requiring exceptionally high doses of insulin [30]. COVID-19 has been associated with the development of new-onset hyperglycemia and worsening of glycemic control in pre-existing diabetes, due to direct pancreatic damage by the virus, body’s stress response to infection, and use of diabetogenic drugs such as corticosteroids in the treatment of severe COVID-19 [35]. Diabetes is one of the major risk factors for fatal outcomes from COVID-19, and a 64% greater risk of diabetes was found in patients with COVID-19 [22]. Patients with diabetes are vulnerable to infection because of hyperglycemia and impaired immune function [15]. During the COVID-19 pandemic, diabetes care was also affected and had a negative impact on people with diabetes, as the healthcare providers were overloaded by COVID-19 cases [12].

Several studies [1,18,20,26,33] suggest a mathematical model to study the dynamics of the spread of COVID-19 in diabetic patients. Mathematical models have been used in literature since decades to examine, analyze, and predict future events and to bring the system to a desirable and controllable state. The real impact of mathematical modeling on public health came with the need for evaluating intervention strategies for newly emerging and re-emerging pathogens [21]. Hence, the study of the interaction between COVID-19 infection and diabetes to understand the spread of epidemic, as well as identification and control, is of utmost importance in current scenario. Kouidere et al. [20] suggested a mathematical modelling of the spread of the COVID-19 pandemic, highlighting the negative impact of quarantine on diabetic population. Khan et al. [18] proposed a first-order nonlinear mathematical model for hypertensive or diabetic patients open to COVID-19. Okyere and Ackora-Prah [24] suggested that more attention should be given to COVID-19 patients with diabetes as they have increased chance of death. Various mathematical modelling of COVID-19 disease has been done extensively in literature with or without comorbidities [2,23,25,31,32].

Our primary focus in this article revolves around examining the impact of COVID-19 on individuals with and without diabetes, incorporating incubation time delays in both populations. In our recent study [1], we developed a mathematical model that initially did not account for time delays. However, we acknowledge that in real-life scenarios, there is a period of time before symptoms of the infection manifest in the host’s body, indicating that it is not an instantaneous process. Therefore, mathematical models utilizing delay differential equations are employed to address this limitation of ordinary differential equation (ODE) models. Delay differential equations involve time derivatives that depend on the solution at previous time points. Over the past few decades, the literature has seen an increased focus on the analysis and control of continuous-time delay systems. The main objective of the model presented in this article is to comprehend the dynamics of COVID-19 infection among individuals with and without diabetes, specifically examining the tolerable level of incubation time delay that allows the system to reach a desirable state.

The outline of the article is as follows: Section 2 introduces our mathematical model incorporating time delays, which represents the susceptible population, as well as the COVID-19-infected populations of both diabetic and non-diabetic individuals. We define the model parameters and derive the characteristic equation. Moving to Section 3, we explore the conditions under which the endemic equilibrium point will be locally stable, considering various scenarios of time delays. Section 4 focuses on the global stability of the endemic equilibrium point. In Section 5, we validate the theoretical results obtained in the previous sections through numerical experiments. Finally, we conclude with a discussion and provide insight into future prospects in the last section.

2 The mathematical model

We present a nonlinear delay mathematical model to study the transmission of COVID-19 infection in diabetic and non-diabetic population. We denote the total susceptible population by x ( t ) ; the infected diabetic population by y ( t ) , and z ( t ) denotes the infected non-diabetic population at time t . Here, we have introduced two different time delays in our model to study their consequences on the dynamic behaviour of the model. Below, we briefly discuss about the biological significance of these time delays.

τ 1 : The time delay τ 1 ≥ 0 has been incorporated to account for the presence of the incubation period of the COVID-19 virus specifically in the diabetic population. Thus, τ 1 represents the duration between infection and the activation or onset of COVID-19 symptoms in individuals with diabetes.
τ 2 : The delay τ 2 ≥ 0 takes into consideration the time delay in the onset of symptoms of COVID-19 in the non-diabetic population.

Figure 1 shows a pictorial representation of our model.

Figure 1

Schematic diagram of the proposed model.

We consider the following system of delay differential equations:

(2.1) x ′ ( t ) = f ( x ) − d x − a f 1 ( x ) y − b f 2 ( x ) z + p α 1 y + q α 2 z , y ′ ( t ) = a 1 α f 1 ( x ( t − τ 1 ) ) y + b 1 β f 2 ( x ( t − τ 2 ) ) z − d 1 y − α 1 y , z ′ ( t ) = a 1 ( 1 − α ) f 1 ( x ( t − τ 1 ) ) y + b 1 ( 1 − β ) f 2 ( x ( t − τ 2 ) ) z − d 2 z − α 2 z ,

with initial populations

x ( t ) = χ ( t ) > 0 , y ( t ) = ζ ( t ) > 0 , z ( t ) = η ( t ) > 0 , t ∈ [ − τ , 0 ] , τ = max { τ 1 , τ 2 } .

In equation (2.1), f ( x ) is the growth rate function of the susceptible population x , whereas f 1 and f 2 are the conversion functions of y and z , respectively. They are converting the susceptible population into infected ones after contact. We assume these functions to be non-negative, as they define the growth rate of the populations. Several examples of these types of functions are available in literature, such as a x 1 + x (Holling Type II) and a x k 1 + x k (Holling Type III). The model consists of three main variables, two time delays, and several parameters. These variables and parameters are important to capture the dynamics occurring in the transmission of COVID-19 among diabetic and non-diabetic groups in a certain population. Here, all the parameters are non-negative owing to the biological significance. The parameters d , d 1 , and d 2 represent the specific death rates of the susceptible, infected diabetic, and infected non-diabetic populations. The death could be natural due to ageing, diseases, etc. for susceptible population, while it could be due to COVID-19 and its comorbidities for infected populations. Rest of the parameters are given in Table 1 for a clear understanding.

Table 1

Notations of the model parameters with their biological significance

Parameters	Biological significance
a	Contact (exposure) rate of y with x
b	Contact (exposure) rate of z with x
a 1	Rate of infection of exposed y
b 1	Rate of infection of exposed z
α	Part of diabetic population infected by y
β	Part of diabetic population infected by z
α 1	Quarantine + treatment rate of y
α 2	Quarantine + treatment rate of z
p	Part of recovered population of y becoming susceptible again
q	Part of recovered population of z becoming susceptible again
d	Specific mortality rate of x
d 1	Specific mortality rate of y
d 2	Specific mortality rate of z

We assume that the functions f , f 1 , and f 2 satisfy the local Lipschitz conditions, and hence, equation (2.1) has unique solutions that are continuable in their maximal intervals of existence for appropriate initial conditions [13]. Furthermore, following arguments as in the study by Rao and Kumar [27], one can see that the solutions of equation (2.1) are also non-negative in their domains of definition to be suitable to represent a biological situation.

2.1 Equilibrium solution and characteristic equation

In general, studies on influence of time delays on biological models consider one or more of five interesting cases: stability invariance, instability invariance, stability change, instability change, and instability switching of equilibria [3,10]. Studies on mathematical models of infectious diseases such as equation (2.1) considers two types of equilibria, viz, a disease-free equilibrium and an epidemic equilibrium. However, in our present study, as we are interested in the environment where the disease, i.e. the COVID-19 infection, persists under the influence of time delays, we consider only the existence of a disease equilibrium and establish conditions for stability of it, implying the prevalence of a disease environment. We assume that the endemic equilibrium solution of the system (2.1), denoted by ( x * , y * , z * ) , exists and is unique, recalling the conditions on parameters and functional relations in the study by Anand et al. [1], as time delays have no impact on the existence of equilibria. For simplicity, hereafter, we assume the growth function f ( x ) as a constant, say f ( x ) = k , throughout the text.

In order to study the local stability of this endemic equilibrium solution, we first linearise the system (2.1) around the endemic equilibrium point ( x * , y * , z * ) . For simplicity, we will denote x ( t − τ 1 ) as x τ 1 and x ( t − τ 2 ) as x τ 2 . The linearised system with respect to the endemic equilibrium point is given as follows:

(2.2) x ′ ( t ) y ′ ( t ) z ′ ( t ) = H A 1 A 2 0 B 1 B 2 0 C 1 C 2 x ( t ) y ( t ) z ( t ) + 0 0 0 α D 1 0 0 ( 1 − α ) D 1 0 0 x τ 1 0 0 + 0 0 0 β D 2 0 0 ( 1 − β ) D 2 0 0 x τ 2 0 0 ,

where H = − d − a f 1 ′ ( x ) y − b f 2 ′ ( x ) z , A 1 = p α 1 − a f 1 ( x ) , A 2 = q α 2 − b f 2 ( x ) , B 1 = a 1 α f 1 ( x τ 1 ) − d 1 − α 1 , B 2 = b 1 β f 2 ( x τ 2 ) , C 1 = a 1 ( 1 − α ) f 1 ( x τ 1 ) , C 2 = b 1 ( 1 − β ) f 2 ( x τ 2 ) − d 2 − α 2 , D 1 = a 1 f 1 ′ ( x τ 1 ) y , and D 2 = b 1 f 2 ′ ( x τ 2 ) z . Let us assume

S = H A 1 A 2 0 B 1 B 2 0 C 1 C 2 , T = 0 0 0 α D 1 0 0 ( 1 − α ) D 1 0 0 , U = 0 0 0 β D 2 0 0 ( 1 − β ) D 2 0 0 .

Then, the characteristic equation corresponding to the endemic equilibrium point ( x * , y * , z * ) is given as follows:

F ( m ) = det ( m I − S − T − U ) = 0 ,

where I is the identity matrix. Hence, the characteristic equation is

(2.3) F ( m ) = m 3 + l 1 m 2 + l 2 m + l 3 + l 4 e − m τ 1 + l 5 e − m τ 2 + l 6 m e − m τ 1 + l 7 m e − m τ 2 ,

where l 1 = − ( H + B 1 + C 2 ) , l 2 = B 1 C 2 − B 2 C 1 + ( B 1 + C 2 ) H , l 3 = − H ( B 1 C 2 − B 2 C 1 ) , l 4 = α A 1 C 2 D 1 − ( 1 − α ) A 1 B 2 D 1 − α A 2 C 1 D 1 + ( 1 − α ) A 2 B 1 D 1 , l 5 = β A 1 C 2 D 2 − ( 1 − β ) A 1 B 2 D 2 − β A 2 C 1 D 2 + ( 1 − β ) A 2 B 1 D 2 , l 6 = − α A 1 D 1 − ( 1 − α ) A 2 D 1 , and l 7 = − β A 1 D 2 − ( 1 − β ) A 2 D 2 .

We shall now proceed to find the stability of the endemic equilibrium point, which enables us to understand the importance of the delay parameters τ 1 and τ 2 and the role they are playing in ensuring the stability of the system in hand.

3 Local stability analysis

In this section, we study the local stability of the endemic equilibrium point ( x * , y * , z * ) of the system (2.1). We will consider the following different cases:

τ 1 > 0 , τ 2 = 0 : Here, we consider that there is a latency or incubation period of the COVID-19 virus in the diabetic population; however, in the non-diabetic population, the symptoms arrive as soon as they are infected.
τ 1 = 0 , τ 2 > 0 : This is the opposite of the above case. That is, there is no latency or incubation period of the COVID-19 virus in the diabetic population, but the non-diabetic population takes some time to show the symptoms.
τ 1 > 0 , τ 2 > 0 : This is the most realistic real-life scenario. Both the populations (diabetic and non-diabetic) are taking time to show the symptoms after getting infected by COVID-19.

Here, we will consider these cases and find the range of the time delays, for which the endemic equilibrium point will be locally stable.

3.1 Delay in incubation of diabetic population ( τ 1 > 0 , τ 2 = 0 )

This is the case where there is no latency or incubation period of the infection in the non-diabetic population, but the diabetic population takes time τ 1 > 0 to show first sign of symptoms.

The characteristic equation (2.3) in this case becomes

(3.1) F ( m , τ 1 ) = m 3 + l 1 m 2 + ( l 2 + l 7 ) m + l 3 + l 5 + ( l 4 + l 6 m ) e − m τ 1 .

We rewrite the characteristic equation (3.1) as F ( m , τ 1 ) = P ( m ) + Q ( m ) e − m τ 1 , where,

(3.2) P ( m ) = m 3 + l 1 m 2 + ( l 2 + l 7 ) m + l 3 + l 5 and Q ( m ) = l 4 + l 6 m .

We will denote the real and imaginary parts of the root m by ℜ ( m ) and ℑ ( m ) , respectively. Let us consider the following result, which will help us estimate delays in such cases.

Lemma 3.1

A system with characteristic equation F ( m , τ 1 ) = P ( m ) + Q ( m ) e − m τ 1 is asymptotically stable if and only if the following conditions are satisfied [14,27]:

F 1 ( m ) = P ( m ) + Q ( m ) ≠ 0 , ℜ ( m ) ≥ 0 ,
F 2 ( m ) = P ( m ) − Q ( m ) ≠ 0 , ℜ ( m ) = 0 , m ≠ 0 ,
F 3 ( m ) = P ( m ) + 1 − m R 1 + m R Q ( m ) ≠ 0 , ℜ ( m ) = 0 , m ≠ 0 , 0 < R < ∞ ,

where R denotes the pseudo-delay and is related to τ 1 as τ 1 = 2 ν ( tan − 1 ( ν R ) + k π ) , ∀ k ∈ Z .

The condition ( i ) in Lemma 3.1 ensures stability for τ 1 = 0 ; hence, by continuity of solutions, this condition ensures stability even for sufficiently small delay τ 1 > 0 ; ( i i ) confirms that no pure imaginary zeros exist in the limiting case; and condition ( i i i ) establishes the absence of imaginary zeros. F 1 , F 2 , and F 3 are obtained by letting e − m τ 1 = 1 , − 1 , and 1 − m R 1 + m R , respectively. Lemma 3.1 is a necessary and sufficient condition, and violation of any of the conditions ( i ) , ( i i ) , or ( i i i ) reflects corresponding change in stability of characteristic equation F ( m , τ 1 ) = P ( m ) + Q ( m ) e − m τ 1 . For more details on the applicability of this lemma, one may refer to [14].

Using Lemma 3.1, we will find the conditions for the preservation of stability in the system. We have, F ( m , τ 1 ) = P ( m ) + Q ( m ) e − m τ 1 . Putting e − m τ 1 = 1 , we obtain

(3.3) F 1 ( m ) = P ( m ) + Q ( m ) = m 3 + l 1 m 2 + ( l 2 + l 7 ) m + l 3 + l 5 + l 4 + l 6 m ≠ 0

for ℜ ( m ) ≥ 0 .

Next, let e − m τ 1 = − 1 , then

(3.4) F 2 ( m ) = P ( m ) − Q ( m ) = m 3 + l 1 m 2 + ( l 2 − l 6 + l 7 ) m + l 3 + l 5 − l 4 .

Let us now put, m = i ν in equation (3.4). Then,

(3.5) F 2 ( i ν ) = − i ν 3 − l 1 ν 2 + i ν ( l 2 − l 6 + l 7 ) + l 3 − l 4 + l 5 .

Clearly, F 2 ( m ) = 0 if and only if ν 2 = l 2 − l 6 + l 7 and ν 2 = l 3 − l 4 + l 5 l 1 . Hence, we can say that, F 2 ( m ) ≠ 0 for ℜ ( m ) = 0 , m ≠ 0 , if the following condition holds ensuring that no purely imaginary zeros exist in the limiting case:

(3.6) ( l 2 − l 6 + l 7 ) l 1 ≠ l 3 − l 4 + l 5 .

Finally, let e − m τ 1 = 1 − m R 1 + m R , for some R > 0 . Then, F 3 ( m ) becomes

F 3 ( m ) = P ( m ) + 1 − m R 1 + m R Q ( m ) .

Substituting P ( m ) and Q ( m ) in the above equation, we obtain

(3.7) F 3 ( m ) = m 4 R + m 3 ( 1 + l 1 R ) + m 2 ( l 1 + R ( l 2 + l 7 − l 6 ) ) 1 + m R + m ( l 2 + l 6 + l 7 + R ( l 3 + l 5 − l 4 ) ) + l 3 + l 4 + l 5 1 + m R .

Now, if we can show that F 3 ( m ) ≠ 0 for ℜ ( m ) = 0 , m ≠ 0 , 0 < R < ∞ , then equation (3.1) has no pure imaginary zeros. In addition, F 3 ( m ) = 0 only if numerator is equal to 0. Hence, we will assume that l 3 + l 4 + l 5 ≠ 0 , otherwise m = 0 . Putting m = i ν in numerator of F 3 ( m ) , we obtain

(3.8) F 3 ( i ν ) = ν 4 R − i ν 3 ( 1 + l 1 R ) − ν 2 ( l 1 + R ( l 2 + l 7 − l 6 ) ) + i ν ( l 2 + l 6 + l 7 + R ( l 3 + l 5 − l 4 ) ) + l 3 + l 4 + l 5 .

Separating real and imaginary parts, we obtain

ℜ ( F 3 ) = ν 4 R − ν 2 ( l 1 + R ( l 2 + l 7 − l 6 ) ) + l 3 + l 4 + l 5 , ℑ ( F 3 ) = − ν 3 ( 1 + l 1 R ) + ν ( l 2 + l 6 + l 7 + R ( l 3 + l 5 − l 4 ) ) .

From the ℑ ( F 3 ) part, we have ν ≠ 0 , hence ν 2 = l 2 + l 6 + l 7 + R ( l 3 − l 4 + l 5 ) 1 + l 1 R . Thus, no real ν exists if l 2 + l 6 + l 7 + R ( l 3 − l 4 + l 5 ) < 0 for any R > 0 , establishing that there is no change in stability for any τ 1 ≥ 0 . Next, using this value of ν 2 and putting in ℜ ( F 3 ) , we obtain

(3.9) M 1 R 3 + N 1 R 2 + P 1 R + Q 1 = 0 ,

where M 1 , N 1 , P 1 , and Q 1 are given as follows:

M 1 = ( l 3 − l 4 + l 5 ) 2 − l 1 ( l 3 − l 4 + l 5 ) ( l 2 − l 6 + l 7 ) ,
N 1 = 2 ( l 2 + l 6 + l 7 ) ( l 3 − l 4 + l 5 ) + l 1 2 ( l 3 + l 4 + l 5 ) + l 1 l 6 2 − l 1 ( l 2 + l 7 ) 2 − ( l 2 − l 6 + l 7 ) ( l 3 − l 4 + l 5 ) − l 1 2 ( l 3 − l 4 + l 5 ) ,
P 1 = ( l 2 + l 6 + l 7 ) 2 + 2 l 1 ( l 3 + l 4 + l 5 ) − ( l 2 + l 7 ) 2 + l 6 2 − l 1 ( l 3 − l 4 + l 5 ) − l 1 2 ( l 2 + l 7 + l 6 ) ,
Q 1 = l 3 + l 4 + l 5 − l 1 ( l 2 + l 6 + l 7 ) .

Hence, we have the following theorem.

Theorem 3.1

Assume the endemic equilibrium solution of equation (2.1) is locally stable for τ 1 = τ 2 = 0 . Then, for any length of delay τ 1 > 0 , τ 2 = 0 , the endemic equilibrium solution is locally stable if the following conditions hold:

( l 2 − l 6 + l 7 ) l 1 ≠ l 3 − l 4 + l 5 ,
M 1 , N 1 , P 1 , Q 1 > 0 ,
N 1 P 1 > M 1 Q 1 .

A change of stability occurs at τ 1 = τ 1 + > 0 if any of the coefficient M 1 , N 1 , P 1 , or Q 1 is negative. Here, τ 1 + = 2 ν + tan − 1 ( ν + R + ) , where

ν + = l 2 + l 6 + l 7 + R ( l 3 − l 4 + l 5 ) 1 + l 1 R

and R + is solution of equation (3.9).

3.2 Delay in incubation of non-diabetic population ( τ 1 = 0 , τ 2 > 0 )

This case corresponds to the scenario where there is no latency or incubation period of the infection in the diabetic population, but the non-diabetic population takes time τ 2 > 0 to show first sign of symptoms.

We have τ 1 = 0 , τ 2 > 0 . Then, equation (2.3) becomes

(3.10) F ( m , τ 2 ) = m 3 + l 1 m 2 + ( l 2 + l 6 ) m + l 3 + l 4 + e − m τ 2 ( l 5 + l 7 m ) .

Let us write the characteristic equation (3.10) as F ( m , τ 2 ) = P ( m ) + Q ( m ) e − m τ 2 , where

(3.11) P ( m ) = m 3 + l 1 m 2 + ( l 2 + l 6 ) m + l 3 + l 4 , Q ( m ) = l 5 + l 7 m .

Similar as the above case, we have F ( m , τ 2 ) = P ( m ) + Q ( m ) e − m τ 2 , and we will use Lemma 3.1 to study the stability of this system. Putting e − m τ 2 = 1 , we obtain

(3.12) F 1 ( m ) = P ( m ) + Q ( m ) = m 3 + l 1 m 2 + ( l 2 + l 6 ) m + l 3 + l 4 + l 5 + l 7 m ≠ 0

for ℜ ( m ) ≥ 0 , ensuring stability for τ 2 = 0 .

Next, let e − m τ 2 = − 1 , then

(3.13) F 2 ( m ) = P ( m ) − Q ( m ) = m 3 + l 1 m 2 + ( l 2 + l 6 − l 7 ) m + l 3 + l 4 − l 5 .

Let us now put, m = i ν . Then,

(3.14) F 2 ( i ν ) = − i ν 3 − l 1 v 2 + i ν ( l 2 + l 6 − l 7 ) + l 3 + l 4 − l 5 .

Clearly, F 2 ( m ) = 0 if ν 2 = l 2 + l 6 − l 7 and ν 2 = l 3 + l 4 − l 5 l 1 . Hence, we can say that, F 2 ( m ) ≠ 0 , and there exists no pure imaginary zeros in limiting case if

(3.15) l 1 ( l 2 + l 6 − l 7 ) ≠ l 3 + l 4 − l 5 .

Finally, let e − m τ 2 = 1 − m R 1 + m R for some R > 0 . Then, F 3 ( m ) becomes

F 3 ( m ) = P ( m ) + 1 − m R 1 + m R Q ( m ) .

Substituting P ( m ) and Q ( m ) in the above equation, we obtain

(3.16) F 3 ( m ) = m 4 R + m 3 ( 1 + l 1 R ) + m 2 ( l 1 + R ( l 2 + l 6 − l 7 ) ) 1 + m R + m ( l 2 + l 6 + l 7 + R ( l 3 + l 4 − l 5 ) ) + l 3 + l 4 + l 5 1 + m R .

Therefore, F 3 ( m ) = 0 only if numerator is equal to 0. Hence, l 3 + l 4 + l 5 ≠ 0 ; otherwise, m = 0 . For pure imaginary zeros, putting m = i ν in numerator of F 3 ( m ) , we obtain

(3.17) F 3 ( i ν ) = ν 4 R − i ν 3 ( 1 + l 1 R ) − ν 2 ( l 1 + R ( l 2 + l 6 − l 7 ) ) + i ν ( l 2 + l 6 + l 7 + R ( l 3 + l 4 − l 5 ) ) + l 3 + l 4 + l 5 .

Separating real and imaginary parts, we obtain

ℜ ( F 3 ) = ν 4 R − ν 2 ( l 1 + R ( l 2 + l 6 − l 7 ) ) + l 3 + l 4 + l 5 , ℑ ( F 3 ) = − ν 3 ( 1 + l 1 R ) + ν ( l 2 + l 6 + l 7 + R ( l 3 + l 4 − l 5 ) ) .

From the ℑ ( F 3 ) part, we have ν ≠ 0 , hence ν 2 = l 2 + l 6 + l 7 + R ( l 3 + l 4 − l 5 ) 1 + l 1 R . Thus, no real ν exists if l 2 + l 6 + l 7 + R ( l 3 + l 4 − l 5 ) < 0 for R > 0 , establishing that stability is preserved for any length of delay τ 2 > 0 . Next, using this value of ν 2 and putting in ℜ ( F 3 ) , we obtain

(3.18) M 2 R 3 + N 2 R 2 + P 2 R + Q 2 = 0 ,

where

M 2 = ( l 3 + l 4 − l 5 ) 2 − l 1 ( l 3 + l 4 − l 5 ) ( l 2 + l 6 − l 7 ) ,
N 2 = 2 ( l 2 + l 6 + l 7 ) ( l 3 + l 4 − l 5 ) + l 1 l 7 2 − l 1 ( l 2 + l 6 ) 2 + l 1 2 ( l 3 + l 4 + l 5 ) − ( l 2 + l 6 − l 7 ) ( l 3 + l 4 − l 5 ) − l 1 2 ( l 3 + l 4 − l 5 ) ,
P 2 = ( l 2 + l 6 + l 7 ) 2 + 2 l 1 ( l 3 + l 4 + l 5 ) − ( l 2 + l 6 ) 2 + l 7 2 − l 1 ( l 3 + l 4 − l 5 ) − l 1 2 ( l 2 + l 7 + l 6 ) ,
Q 2 = l 3 + l 4 + l 5 − l 1 ( l 2 + l 6 + l 7 ) .

We have the following theorem.

Theorem 3.2

Assume the endemic equilibrium solution to equation (2.1) is locally stable for τ 1 = τ 2 = 0 . Then, for any length of delay τ 1 = 0 , τ 2 > 0 , the endemic equilibrium solution is locally stable if the following conditions hold:

l 1 ( l 2 + l 6 − l 7 ) ≠ l 3 + l 4 − l 5 ,
M 2 , N 2 , P 2 , Q 2 ≥ 0 ,
N 2 P 2 > M 2 Q 2 .

A change of stability occurs at τ 2 = τ 2 + > 0 if any of the coefficient M 2 , N 2 , P 2 , or Q 2 is negative. Here, τ 2 + = 2 ν + tan − 1 ( ν + R + ) , where

ν + = l 2 + l 6 + l 7 + R ( l 3 + l 4 − l 5 ) 1 + l 1 R

and R + is solution to equation (3.18).

3.3 Delays in incubation of both diabetic and non-diabetic populations ( τ 1 , τ 2 > 0 )

This is the most important and relevant scenario, where the incubation periods τ 1 and τ 2 of both the populations are greater than zero.

Letting m = μ + i ν in equation (2.3), then μ and ν are real solutions of

ℜ ( F ( μ , ν ) ) = μ 3 − 3 μ ν 2 + l 1 ( μ 2 − ν 2 ) + l 2 μ + l 3 + l 4 cos ( ν τ 1 ) e − τ 1 μ + l 5 cos ( ν τ 2 ) e − τ 2 μ + l 6 ( μ cos ( ν τ 1 ) + ν sin ( ν τ 1 ) ) e − τ 1 μ + l 7 ( μ cos ( ν τ 2 ) + ν sin ( ν τ 2 ) ) e − τ 2 μ , ℑ ( F ( μ , ν ) ) = 3 μ 2 ν − ν 3 + 2 l 1 μ ν + l 2 ν − l 4 sin ( ν τ 1 ) e − τ 1 μ − l 5 sin ( ν τ 2 ) e − τ 2 μ + l 6 ( ν cos ( ν τ 1 ) − μ sin ( ν τ 1 ) ) e − τ 1 μ + l 7 ( ν cos ( ν τ 2 ) − μ sin ( ν τ 2 ) ) e − τ 2 μ .

To find the condition for change of stability, we substitute μ = 0 in the above equations to obtain

(3.19) ℜ ( F ( i ν ) ) = − l 1 ν 2 + l 3 + l 4 cos ( ν τ 1 ) + l 5 cos ( ν τ 2 ) + l 6 ν sin ( ν τ 1 ) + l 7 ν sin ( ν τ 2 ) ,

(3.20) ℑ ( F ( i ν ) ) = − ν 3 + l 2 ν − l 4 sin ( ν τ 1 ) − l 5 sin ( ν τ 2 ) + l 6 ν cos ( ν τ 1 ) + l 7 ν cos ( ν τ 2 ) .

For the preservation of the stability of the delay system, we use the following conditions of the Nyquist criteria ℑ ( F ( i ν 0 ) ) > 0 and ℜ ( F ( i ν 0 ) ) = 0 , where ν 0 is the smallest positive root of ℜ ( F ( i ν 0 ) ) = 0 . For more details, we refer to [7,9].

Let ν = ν 0 be the smallest positive root of ℜ ( F ( i ν ) ) = 0 . We write ℜ ( F ( i ν ) ) = 0 as

(3.21) l 1 ν 2 − l 3 = l 4 cos ( ν τ 1 ) + l 5 cos ( ν τ 2 ) + l 6 ν sin ( ν τ 1 ) + l 7 ν sin ( ν τ 2 ) .

After replacing sin and cos by their maximum values, we rewrite the above equation (3.21) as follows:

l 1 ν 2 − l 3 < ∣ l 4 ∣ + ∣ l 5 ∣ + ∣ ν ∣ ∣ l 6 ∣ + ∣ ν ∣ ∣ l 7 ∣ .

Let us consider the equation

(3.22) l 1 ν 0 2 − ν 0 ( ∣ l 6 ∣ + ∣ l 7 ∣ ) − l 3 − ∣ l 4 ∣ − ∣ l 5 ∣ = 0 .

The quadratic equation (3.22) has two roots, and let ν + be the bigger root, which is given as follows:

(3.23) ν + = ∣ l 6 ∣ + ∣ l 7 ∣ + ∣ l 6 ∣ 2 + ∣ l 7 ∣ 2 + 2 ∣ l 6 ∣ ∣ l 7 ∣ + 4 l 1 ( l 3 + ∣ l 4 ∣ + ∣ l 5 ∣ ) 2 l 1 .

We now use ℑ ( F ( i ν 0 ) ) > 0 and obtain

(3.24) − ν 3 + l 2 ν + l 6 ν cos ( ν τ 1 ) + l 7 ν cos ( ν τ 2 ) > l 4 sin ( ν τ 1 ) + l 5 sin ( ν τ 2 ) .

Let, ϕ = − ν 3 + l 2 ν + l 6 ν cos ( ν τ 1 ) + l 7 ν cos ( ν τ 2 ) and ψ = l 4 sin ( ν τ 1 ) + l 5 sin ( ν τ 2 ) . Then, clearly ϕ > ψ .

We now need to find the range of the delay for which the system will be stable. We observe that ψ < ( ∣ l 4 ∣ + ∣ l 5 ∣ ) ν + ( τ 1 + τ 2 ) . Hence,

ψ ( τ 1 + τ 2 ) ν + ≤ ∣ l 4 ∣ + ∣ l 5 ∣ .

In addition,

ϕ ν + ≤ ν + 2 + ∣ l 2 ∣ + ∣ l 6 ∣ cos ( ν + τ 1 ) + ∣ l 7 ∣ cos ( ν + τ 2 ) < ν + 2 + ∣ l 2 ∣ + ∣ l 6 ∣ + ∣ l 7 ∣ .

Therefore,

ϕ ( τ 1 + τ 2 ) ν + ≤ ν + 2 + ∣ l 2 ∣ + ∣ l 6 ∣ + ∣ l 7 ∣ τ 1 + τ 2 .

Hence, we can write

∣ l 4 ∣ + ∣ l 5 ∣ < ν + 2 + ∣ l 2 ∣ + ∣ l 6 ∣ + ∣ l 7 ∣ τ 1 + τ 2 .

Thus, we obtain

(3.25) ( τ 1 + τ 2 ) < ν + 2 + ∣ l 2 ∣ + ∣ l 6 ∣ + ∣ l 7 ∣ ∣ l 4 ∣ + ∣ l 5 ∣ .

Hence, we can formulate the following theorem.

Theorem 3.3

Assume that the endemic equilibrium point of equation (2.1) is locally stable for τ 1 = τ 2 = 0 . Then, the stability of the endemic equilibrium point is preserved for τ 1 , τ 2 > 0 if

( τ 1 + τ 2 ) < ν + 2 + ∣ l 2 ∣ + ∣ l 6 ∣ + ∣ l 7 ∣ ∣ l 4 ∣ + ∣ l 5 ∣ .

Remark

Note that condition (i) of Lemma 3.1 discusses the case of delay-free system τ 1 = τ 2 = 0 . Thus, from equation (3.3), we have

(3.26) F ( m ) = m 3 + l 1 m 2 + ( l 2 + l 6 + l 7 ) m + l 3 + l 4 + l 5 .

Hence, if the parameters of the system satisfy the conditions

l 1 > 0 , l 2 + l 6 + l 7 > 0 , l 3 + l 4 + l 5 > 0 , l 1 ( l 2 + l 6 + l 7 ) > l 3 + l 4 + l 5 ,

then equation (3.26) will have all the roots with negative real parts. Therefore, the system (2.1) is asymptotically stable for τ 1 = τ 2 = 0 if the above conditions hold true. We observe that our results are matching with those of Theorem 3.3 of [1].

We now proceed to find the global stability of the endemic equilibrium point of equation (2.1).

4 Global stability analysis

Global stability is an important study in disease dynamics, which specifies the stability of the equilibrium in the broader sense – whether the equilibrium point is always stable irrespective of the changes in the initial conditions. In this section, we will deal with the global stability aspect of the endemic equilibrium point. We make the following assumptions on the conversion functions f 1 and f 2 before we establish our results.

(4.1) 0 ≤ ℳ 1 ≤ f 1 ( x ) ≤ N 1 , 0 ≤ ℳ 2 ≤ f 2 ( x ) ≤ N 2 ,

for all x . These assumptions are realistic in view of examples of conversion functions that are commonly used in such biological scenarios mentioned in Section 2.

4.1 Global stability analysis for endemic equilibrium point

Let ( x * , y * , z * ) be the endemic equilibrium point. Then, they satisfy the following equations:

(4.2) 0 = k − d x * − a f 1 ( x * ) y * − b f 2 ( x * ) z * + p α 1 y * + q α 2 z * ,

(4.3) 0 = a 1 α f 1 ( x * ) y * + b 1 β f 2 ( x * ) z * − ( d 1 + α 1 ) y * ,

(4.4) 0 = a 1 ( 1 − α ) f 1 ( x * ) y * + b 1 ( 1 − β ) f 2 ( x * ) z * − ( d 2 + α 2 ) z * .

Subtracting equations (4.2)–(4.4) from equation (2.1) and rearranging, we obtain

(4.5) ( x − x * ) ′ = − d ( x − x * ) − a f 1 ( x ) ( y − y * ) − a y * ( f 1 ( x ) − f 1 ( x * ) ) − b f 2 ( x ) ( z − z * ) − b z * ( f 2 ( x ) − f 2 ( x * ) ) + p α 1 ( y − y * ) + q α 2 ( z − z * ) ,

(4.6) ( y − y * ) ′ = a 1 α f 1 ( x τ 1 ) ( y − y * ) + a 1 α y * ( f 1 ( x τ 1 ) − f 1 ( x * ) ) + b 1 β f 2 ( x τ 2 ) ( z − z * ) + b 1 β z * ( f 2 ( x τ 2 ) − f 2 ( x * ) ) − d 1 ( y − y * ) − α 1 ( y − y * ) ,

(4.7) ( z − z * ) ′ = a 1 ( 1 − α ) f 1 ( x τ 1 ) ( y − y * ) + a 1 ( 1 − α ) y * ( f 1 ( x τ 1 ) − f 1 ( x * ) ) + b 1 ( 1 − β ) f 2 ( x τ 2 ) ( z − z * ) + b 1 ( 1 − β ) z * ( f 2 ( x τ 2 ) − f 2 ( x * ) ) − d 2 ( z − z * ) − α 2 ( z − z * ) .

Let us assume that there exist positive real constants S 1 and S 2 such that

(4.8) ∣ f 1 ( x ) − f 1 ( x * ) ∣ ≤ S 1 ∣ x − x * ∣ , ∣ f 2 ( x ) − f 2 ( x * ) ∣ ≤ S 2 ∣ x − x * ∣ ,

where x ≠ x * . We are now in a position to establish global stability of the endemic equilibrium. We employ a Lyapunov-Krasovskii-type functional to cancel out the delay terms and manipulate the delay-free terms using inequality analysis [17]. We note that the functional employed in [27] is also applicable here. Now we have the following theorem.

Theorem 4.1

The endemic equilibrium point ( x * , y * , z * ) of the system (2.1)is globally asymptotically stable if the following conditions hold:

d − ( a − a 1 ) y * S 1 − ( b − b 1 ) z * S 2 > 0 ,
d 1 + ( 1 − p ) α 1 + a ℳ 1 − a 1 N 1 > 0 , and
d 2 + ( 1 − q ) α 2 + b ℳ 2 − b 1 N 2 > 0 ,

where ℳ 1 , N 1 , and ℳ 2 , N 2 are positive real constants defined above.

Proof

We consider the functional

(4.9) V ( t ) ≡ V ( x ( t ) , y ( t ) , z ( t ) ) = ∣ x − x * ∣ + ∣ y − y * ∣ + ∣ z − z * ∣ + a 1 y * ∫ t − τ 1 t ∣ f 1 ( x ( u ) ) − f 1 ( x * ) ∣ d u + b 1 z * ∫ t − τ 2 t ∣ f 2 ( x ( u ) ) − f 2 ( x * ) ∣ d u .

Clearly, V ( x * , y * , z * ) = 0 and V ( t ) ≥ 0 . The upper Dini derivative of V along the solutions of the system (2.1) is given as follows:

D + V ≤ − d ∣ x − x * ∣ − a f 1 ( x ) ∣ y − y * ∣ − a y * ∣ f 1 ( x ) − f 1 ( x * ) ∣ − b f 2 ( x ) ∣ z − z * ∣ − b z * ∣ f 2 ( x ) − f 2 ( x * ) ∣ + p α 1 ∣ y − y * ∣ + q α 2 ∣ z − z * ∣ + a 1 α f 1 ( x τ 1 ) ∣ y − y * ∣ − d 1 ∣ y − y * ∣ + a 1 α y * ∣ f 1 ( x τ 1 ) − f 1 ( x * ) ∣ + b 1 β f 2 ( x τ 2 ) ∣ z − z * ∣ + b 1 β z * ∣ f 2 ( x τ 2 ) − f 2 ( x * ) ∣ − α 1 ∣ y − y * ∣ + a 1 ( 1 − α ) f 1 ( x τ 1 ) ∣ y − y * ∣ + a 1 ( 1 − α ) y * ∣ f 1 ( x τ 1 ) − f 1 ( x * ) ∣ + b 1 ( 1 − β ) f 2 ( x τ 2 ) ∣ z − z * ∣ + b 1 ( 1 − β ) z * ∣ f 2 ( x τ 2 ) − f 2 ( x * ) ∣ − d 2 ∣ z − z * ∣ − α 2 ∣ z − z * ∣ + a 1 y * ∫ t − τ 1 t ∣ f 1 ( x ( u ) ) − f 1 ( x * ) ∣ d u + b 1 z * ∫ t − τ 2 t ∣ f 2 ( x ( u ) ) − f 2 ( x * ) ∣ d u .

After simplifying, we obtain

(4.10) D + V ≤ − d ∣ x − x * ∣ − ( d 1 + ( 1 − p ) α 1 ) ∣ y − y * ∣ − ( d 2 + ( 1 − q ) α 2 ) ∣ z − z * ∣ − ∣ y − y * ∣ ∣ a f 1 ( x ) − a 1 f 1 ( x τ 1 ) ∣ − ∣ z − z * ∣ ∣ b f 2 ( x ) − b 1 f 2 ( x τ 2 ) ∣ − ( a − a 1 ) y * ∣ f 1 ( x ) − f 1 ( x * ) ∣ − ( b − b 1 ) z * ∣ f 2 ( x ) − f 2 ( x * ) ∣ .

We rewrite the above equation (4.10) by using the bounds (4.1) and (4.8), and we arrive at

D + V < − ∣ x − x * ∣ ( d − ( a − a 1 ) y * S 1 − ( b − b 1 ) z * S 2 ) − ( d 1 + ( 1 − p ) α 1 + a f 1 ( x ) − a 1 N 1 ) ∣ y − y * ∣ − ( d 2 + ( 1 − q ) α 2 + b f 2 ( x ) − b 1 N 2 ) ∣ z − z * ∣ ≤ − ( d − ( a − a 1 ) y * S 1 − ( b − b 1 ) z * S 2 ) ∣ x − x * ∣ − ( d 1 + ( 1 − p ) α 1 + a ℳ 1 − a 1 N 1 ) ∣ y − y * ∣ − ( d 2 + ( 1 − q ) α 2 + b ℳ 2 − b 1 N 2 ) ∣ z − z * ∣ < 0

by assumptions (i)–(iii). Thus, V ( t ) is negative definite. Furthermore, by definition V ( t ) ≥ 0 and V ( t ) = V ( x , y , z ) = 0 if and only if x = x * , y = y * , z = z * . Thus, V ( t ) is a suitable Lyapunov functional. This completes the proof.□

We shall now present some illustrative examples in order to verify the results.

5 Numerical experiments

In this section, we perform some numerical experiments to validate our theoretical results using MATLAB software dde23. Since, we could not obtain a single source of data to reflect all parameters, we have assumed the parameters to illustrate the results. We perform simulations for various lengths of delays in all three cases, i.e. (i) τ 1 = τ 2 = 0 , (ii) τ 1 > 0 , τ 2 = 0 , (iii) τ 1 = 0 , τ 2 > 0 , and (iv) τ 1 > 0 , τ 2 > 0 . The motive of this study is to find the incubation time of the COVID-19 virus in the diabetic and non-diabetic population, depending on the different parameters. We will assume the units of the parameters to be as per day.

Example 1

We consider the system (2.1) with high exposure rates a and b and low conversion rates a 1 and b 1 . We also take the treatment parameters α 1 and α 2 to be nearly equal to the conversion rates, and the mortality rates ( d , d 1 , d 2 ) are taken as low values. Let f 1 ( x ) = f 2 ( x ) = x 1 + x and α = β . Let the parameter values be k = 300 , d = 0.0001 , a = 2 , b = 3 , p = 0.5 , q = 0.2 , α 1 = 0.7 , α 2 = 0.5 , α = 0.5 , a 1 = 0.75 , b 1 = 0.5 , d 1 = 0.01 , and d 2 = 0.008 . The initial populations are x ( t ) = 100 , y ( t ) = 20 , z ( t ) = 80 , t ∈ [ − τ , 0 ] , and τ = max { τ 1 , τ 2 } . The endemic equilibrium point of the system (2.1) is given by ( x * , y * , z * ) = ( 49 , 54 , 75 ) .

Delay free. Let τ 1 = τ 2 = 0 . Then, we see that all the conditions for delay-free system are satisfied. The solution profile of the delay-free case is given in Figure 2.

$Figure 2 Solution profile for τ 1 = τ 2 = 0 {\tau }_{1}={\tau }_{2}=0 of Example 1.$

Figure 2

Solution profile for τ 1 = τ 2 = 0 of Example 1.

Case 1 ( τ 1 > 0 , τ 2 = 0 ). We calculate M 1 = 0.00009 , N 1 = 0.016 , P 1 = − 0.0046 , and Q 1 = − 0.066 from Theorem 3.1. We see that there is one sign change in the equation M 1 R 3 + N 1 R 2 + P 1 R + Q 1 = 0 , and hence from Routh-Hurwitz criteria, there exists one root R with positive real part, which is coming to be R = 2.17706 . With this value of R , we obtain ν + = 0.236 and τ 1 + = 4.022 from Theorem 3.1. Hence, theoretically, the system will be stable till τ 1 + ≤ 4.022 for τ 2 = 0 and will be unstable when τ 1 + > 4.022 . Numerically, we can see from Figure 3, the solutions are approaching the equilibrium solution ( x * , y * , z * ) for τ 1 = 4 , whereas for τ 1 = 4.2 , the solutions does not converge. Hence, stability change is clearly evident as per theory.

$Figure 3 Solution profiles for τ 1 > 0 , τ 2 = 0 {\tau }_{1}\gt 0,{\tau }_{2}=0 .$

Figure 3

Solution profiles for τ 1 > 0 , τ 2 = 0 .

Case 2 ( τ 1 = 0 , τ 2 > 0 ). Here, M 2 = − 0.00009 , N 2 = 0.014 , P 2 = − 0.0075 , and Q 2 = − 0.066 from Theorem 3.2. Similar to the above case, we see that there are two sign changes in M 2 R 3 + N 2 R 2 + P 2 R + Q 2 = 0 , and hence there exists two R with positive real parts. The positive roots are R = 2.4528 and R = 155.58 . With R = 155.58 , we obtain τ 2 + = 50.4 . This high value of delay τ 2 may represent the case where the non-diabetic population is a carrier (strong immunity and/or does not show symptoms easily) but not identified as infected for a long time. Since, we are more interested in finding the incubation period after which the population will show symptom, we will not proceed with this case. With R = 2.4528 , we obtain ν + = 0.2328 and τ 2 + = 4.4574 from Theorem 3.2. Hence, the system is stable for τ 2 + ≤ 4.4574 when τ 1 = 0 . We see that numerically, from Figure 4, a change in τ 2 from 4.45 to 4.7 changes the stability of the system. Hence, our theoretical results are validated.

$Figure 4 Solution profiles for τ 1 = 0 , τ 2 > 0 {\tau }_{1}=0,{\tau }_{2}\gt 0 .$

Figure 4

Solution profiles for τ 1 = 0 , τ 2 > 0 .

Case 3 ( τ 1 , τ 2 > 0 ). From Theorem 3.3, we find τ 1 + τ 2 < 5.38 . Hence, the system should be stable till τ 1 + τ 2 < 5.38 and unstable when τ 1 + τ 2 ≥ 5.38 . From Figure 5, we see that the system is stable till τ 1 + τ 2 ≤ 3.6 .

$Figure 5 Solution Profiles for τ 1 , τ 2 > 0 {\tau }_{1},{\tau }_{2}\gt 0 .$

Figure 5

Solution Profiles for τ 1 , τ 2 > 0 .

Discussion. From our chosen set of parameter values, we see that the exposure rates a and b are much higher than the conversion rates a 1 and b 1 , respectively. We took the treatment parameters α 1 and α 2 nearly equal to the conversion rates and the death rates are comparatively small. With these parameter values, we see that the system is becoming stable after long time, and if the delays or incubation time is crossing their threshold value, then the system starts oscillation. The delays or incubation periods estimated in this Example 1 are close to the real-life scenario [36]. We also note from the Figures 3–5 that the oscillation due to the delays are only in the x , i.e. the susceptible population, and there is no oscillation in the y and z populations. This is due to the fact that we have introduced the delay parameters only in the susceptible population and none in the infected diabetic and non-diabetic populations. As we are increasing the delay range, oscillations with higher amplitude are observed in the x population – a chaotic situation continues. In addition, the condition of global stability in Theorem 4.1 is not satisfied, and hence, the system is not globally stable around the equilibrium point ( x * , y * , z * ) = ( 49 , 54 , 75 ) . That is, a small change in the initial conditions might give rise to instability around the equilibrium point ( x * , y * , z * ) = ( 49 , 54 , 75 ) .

Example 2

In this example, we alter the conversion parameters (i.e. reduce a 1 and increase b 1 ) while increasing the proportions of recovered turning susceptible diabetic and non-diabetic populations ( p and q ). We let f 1 ( x ) = f 2 ( x ) = x 1 + x and α = β = 0.5 , k = 300 , d = 0.0001 , a = 2 , b = 3 , p = 0.7 , q = 0.3 , α 1 = 0.7 , α 2 = 0.5 , a 1 = 0.6 , b 1 = 0.6 , d 1 = 0.01 , and d 2 = 0.008 . The initial populations are x ( t ) = 100 , y ( t ) = 20 , z ( t ) = 80 , t ∈ [ − τ , 0 ] , and τ = max { τ 1 , τ 2 } . The endemic equilibrium point is given by ( x * , y * , z * ) = ( 76 , 55 , 77 ) .

Delay free. Let τ 1 = τ 2 = 0 . Then, all the conditions for delay-free system are satisfied. The solution profile of the delay-free case is given in Figure 6.

$Figure 6 Solution profile for τ 1 = τ 2 = 0 {\tau }_{1}={\tau }_{2}=0 of Example 2.$

Figure 6

Solution profile for τ 1 = τ 2 = 0 of Example 2.

Case 1 ( τ 1 > 0 , τ 2 = 0 ). We calculate M 1 = − 0.000075 , N 1 = 0.0056 , P 1 = − 0.0056 , and Q 1 = − 0.0257 from Theorem 3.1. We see that there are two sign changes in the equation M 1 R 3 + N 1 R 2 + P 1 R + Q 1 = 0 , and hence, from Routh-Hurwitz criteria, there exists two roots R with positive real parts. The roots are R = 2.7666 and R = 73.58 . With R = 73.58 , we obtain very high values of τ 1 ; hence, as discussed in Case 2 of the previous example, we will not proceed with this scenario. With R = 2.7666 , we obtain ν + = 0.16 and τ 1 + = 5.2135 from Theorem 3.1. Hence, the system will be stable till τ 1 + ≤ 5.2135 given τ 2 = 0 and will be unstable when τ 1 + > 5.2135 . Numerically, we see from Figure 7 that a change in τ 1 from 5 to 5.5 changes the stability of the system.

$Figure 7 Solution profiles for τ 1 > 0 , τ 2 = 0 {\tau }_{1}\gt 0,{\tau }_{2}=0 .$

Figure 7

Solution profiles for τ 1 > 0 , τ 2 = 0 .

Case 2 ( τ 1 = 0 , τ 2 > 0 ). Here, M 2 = 0.0000738 , N 2 = 0.0083 , P 2 = − 0.00087 , and Q 2 = − 0.025 from Theorem 3.2. Similar to the above example, we see that there is one sign change in the coefficients of M 2 R 3 + N 2 R 2 + P 2 R + Q 2 = 0 , and hence there exists one R with positive real part, which is R = 1.7984 . With this R , we obtain ν + = 0.166 and τ 2 + = 3.5 from Theorem 3.2. Hence, the system is stable for τ 2 + ≤ 3.5 given τ 1 = 0 . From Figure 8, we observe that a change in τ 2 from 3.3 to 3.6 changes the stability of the system.

$Figure 8 Solution profiles for τ 1 = 0 , τ 2 > 0 {\tau }_{1}=0,{\tau }_{2}\gt 0 .$

Figure 8

Solution profiles for τ 1 = 0 , τ 2 > 0 .

Case 3 ( τ 1 , τ 2 > 0 ). From Theorem 3.3, we obtain τ 1 + τ 2 < 5.413 . Hence, the system is stable till τ 1 + τ 2 < 5.413 and should be unstable for τ 1 + τ 2 ≥ 5.413 . From Figure 9, we see that the system is stable till τ 1 + τ 2 ≤ 3.7 .

$Figure 9 Solution Profiles for τ 1 , τ 2 > 0 {\tau }_{1},{\tau }_{2}\gt 0 .$

Figure 9

Solution Profiles for τ 1 , τ 2 > 0 .

Discussion. In this example, we changed the conversion parameters slightly, i.e. reduced a 1 (conversion rate corresponding to the diabetic population) and increased b 1 (conversion rate corresponding to the non-diabetic population). Compared to Example 1, there have been changes in the delays (with τ 1 increasing and τ 2 slightly decreasing). These changes indicate that more non-diabetic individuals exposed to COVID-19 are transitioning into the infected population, implying a weaker immune response and a shorter incubation period ( τ 2 ) compared to Example 1. On the other hand, the diabetic population exposed to the COVID-19 virus shows a lower rate of transition into the infected population compared to Example 1, suggesting better immunity due to medication. Consequently, the diabetic population has a longer incubation period compared to their non-diabetic counterparts. Similar to Example 1, we see from the Figures 7–9 that the oscillations are only in x , and none in y and z , and oscillations with larger amplitude are observed in the x population as we increase the delay values leading to a chaotic situation. The explanation of this phenomena is given in the discussion of Example 1. In addition, in this example, the condition of global stability in Theorem 4.1 is not satisfied, and hence, the system is not globally stable around the equilibrium point ( x * , y * , z * ) = ( 76 , 55 , 77 ) .

Next, let us consider an example where the endemic equilibrium point ( x * , y * , z * ) is globally asymptotically stable.

Example 3

In this example, we take the parameter values in the system (2.1) as k = 300 , d = 0.8 , a = 1.2534 , b = 2.1431 , p = 0.5 , q = 0.7 , α 1 = 5 , α 2 = 1 , α = 0.5 , a 1 = 1.253 , b 1 = 2.142 , d 1 = 0.5 , d 2 = 0.15 , f 1 ( x ) = f 2 ( x ) = x 1 + x , and α = β .

The endemic equilibrium point is given by ( x * , y * , z * ) = ( 22 , 55 , 262 ) . We have chosen the parameters in such a way that the global stability of the endemic equilibrium point ( x * , y * , z * ) = ( 22 , 55 , 262 ) is satisfied, i.e. whenever the exposure ( a , b ) and conversion ( a 1 , b 1 ) rates are nearly equal, the treatment rates ( α 1 , α 2 ) are high compared to the mortality rates ( d 1 , d 2 ). We see that all the conditions of Theorem 4.1 are satisfied, and hence, our system is globally asymptotically stable.

We can see from Figure 10 that the stability of the delay-free system is independent of the initial populations, implying global stability.

$Figure 10 Solution profiles for τ 1 = τ 2 = 0 {\tau }_{1}={\tau }_{2}=0 at different initial conditions.$

Figure 10

Solution profiles for τ 1 = τ 2 = 0 at different initial conditions.

Next, let us consider the system where both the delays exist, i.e. τ 1 , τ 2 > 0 . From Theorem 3.3, we obtain τ 1 + τ 2 < 2.99 . Hence, the system is stable till τ 1 + τ 2 < 2.99 . Below, we have given the following plots in Figure 11 for different initial populations for τ 1 , τ 2 > 0 .

$Figure 11 Corresponding solution profiles for τ 1 , τ 2 > 0 {\tau }_{1},{\tau }_{2}\gt 0 for different initial states.$

Figure 11

Corresponding solution profiles for τ 1 , τ 2 > 0 for different initial states.

Discussion. From Figures 10 and 11, we see that system behaves similarly irrespective of the presence of delays and change in initial conditions. We observe that for the case where τ 1 , τ 2 > 0 , the system is globally stable whenever τ 1 + τ 2 < 2.99 .

Remark

From Examples 1 and 2, we observe how the parameters play an important role in influencing the time delays (incubation periods, in our case) of the system. In Example 1, our conversion rates from susceptibles to COVID-19-infected diabetic and non-diabetic populations were a 1 = 0.75 and b 1 = 0.5 , whereas in Example 2, we had a 1 = 0.6 and b 1 = 0.6 . The exposure to the COVID-19 infection was same for both the populations in both the examples ( a = 2 , b = 3 ). In Example 1, we observed that the theoretical value of delay τ 1 (incubation period of COVID-19 infection in diabetic population) is 4.022, whereas in Example 2, we obtain τ 1 = 5.2 . This is due to the fact that in Example 2, less diabetic people are getting converted to the infected population ( a 1 = 0.6 as compared to a 1 = 0.75 in Example 1), and hence, their immunity is comparatively better, due to which the onset of symptoms is taking longer as compared to Example 1. Similarly, the delay τ 2 = 3.5 in Example 2, as compared to τ 2 = 4.45 in Example 1. This is again due to the fact that our conversion parameter b 1 in Example 1 is smaller (hence less people getting converted to infected population) as compared to Example 2. In Example 3, we showcased a scenario where the system will be globally stable for all initial populations irrespective of time delays, which is only possible if the treatment parameters ( α 1 and α 2 ) are very high as compared to exposure rates ( a and b ) and conversion rates ( a 1 and b 1 ). This is relevant to real-life scenario also. Example 3 illustrates the case where the system is globally stable, i.e. independent of initial states. Furthermore, these conditions are violated in Examples 1 and 2, implying that global and local stability need not go hand-in-hand usually. On the whole, our examples illustrate the results established and explain the contribution of parameters in stabilising or destabilising the system under the influence of time delays.

6 Conclusions and future work

In this article, we have focused on discussing the COVID-19 virus, which has emerged as the most prevalent and dangerous viral infection of this decade. Specifically, we have examined its impact on individuals who are diabetic, a chronic health condition that is currently widespread. To capture the real-life dynamics of the virus within the host’s body, we have proposed a three-compartment model incorporating two time delays. These time delays have been introduced to enhance the model and represent the incubation period of the COVID-19 virus. We emphasize that ODE compartmental models do not adequately represent the transition of a susceptible individual to the infected compartment, as it is not an instantaneous process. Therefore, we have developed a delay compartmental model to better capture this aspect. In our study, we have conducted a local stability analysis of the model, where we derived upper bounds for the delay parameters τ 1 and τ 2 . These bounds have been numerically verified in Section 5. Furthermore, we have performed a global stability analysis using a Lyapunov function. Finally, we have validated our theoretical findings through numerical experiments. We chose our model parameters in such a way that the incubation periods derived should be realistic. Numerous meta-analyses [34,36] estimate the incubation time of the COVID-19 virus to be around 5 days, with a range of 1.80–18.87 days. In particular, Wu et al. [36] provided incubation periods for different variants of the COVID-19 virus. From their meta-analysis of 141 articles, they found the incubation periods for the Alpha, Beta, Delta, and Omicron variants to be 5.00, 4.50, 4.41, and 3.42 days, respectively. These values closely align with the incubation periods obtained in our examples using synthetic data. Therefore, our model, with the chosen parameters, accurately reflects the real-life scenario.

We have emphasised on the incubation periods of the virus in populations; however, the treatment period is also of significant importance. Hence, a model with recovery or treatment delay is left for future endeavor. Also, introduction of control/preventive measures such as vaccination or face mask usage in equation (2.1) could be an interesting study in the scenario of endemic prevalent conditions established here. Hence, we can further venture in this direction in our next studies. It is interesting to note that an increase in length of time delay in either of infected populations led to oscillations in susceptible population but has no impact on infected populations. This phenomenon needs to be explored further and verified with real-time environment. Authors noted that increase of lengths of time delays beyond the estimates obtained in the results leads only to disturb the stability of the endemic equilibrium leading a chaotic state. This is due to the fact that we have restricted our study to conditions of Lemma 3.1, which focuses only on non-changing of stability of the system, and parameters and functions are chosen to fit these conditions only. However, the existence of multiple bifurcation points and change of stability or instability status of the system are also interesting phenomena to explore for dynamical systems such as equation (2.1). We leave this for a future study.

Acknowledgements

The authors are thankful to the anonymous reviewers for their careful reading and valuable suggestions that led to improvements in the manuscript.

Funding information: This research received no specific grant from any funding agency, commercial or nonprofit sectors.
Conflict of interest: The authors state that there is no conflict of interest.
Ethical approval: This research did not require ethical approval.

References

[1] Anand, M., Danumjaya, P., & Rao, P. R. S. (2023). A nonlinear mathematical model on the COVID-19 transmission pattern among diabetic and non-diabetic population. Mathematics and Computers in Simulation, 210, 346–369. 10.1016/j.matcom.2023.03.016Search in Google Scholar PubMed PubMed Central

[2] Arshad, S., Siddique, I., Nawaz, F., Shaheen, A., & Khurshid, H. (2023). Dynamics of a fractional order mathematical model for COVID-19 epidemic transmission. Physica A: Statistical Mechanics and its Applications, 609, 128383. 10.1016/j.physa.2022.128383Search in Google Scholar PubMed PubMed Central

[3] Barman, B., & Ghosh, B. (2019). Explicit impacts of harvesting in delayed predator-prey models. Chaos Solitons Fractals, 122, 213–228. 10.1016/j.chaos.2019.03.002Search in Google Scholar

[4] Bhandari, S., Singh, A., Sharma, R., et al. (2020). Characteristics, treatment outcomes and role of hydroxychloroquine among 522 COVID-19 hospitalized patients in Jaipur City: An epidemio-clinical study. Journal of the Association of Physicians of India, 68(6), 13–19. Search in Google Scholar

[5] Cariou, B., Wargny, M., Boureau, A. S., et al. (2022). Impact of diabetes on COVID-19 prognosis beyond comorbidity burden: the CORONADO initiative. Diabetologia, 65, 1436–1449. 10.1007/s00125-022-05734-1Search in Google Scholar PubMed PubMed Central

[6] Cascella, M., Rajnik, M., Aleem, A., et al. Features, evaluation, and treatment of coronavirus (COVID-19) [Updated 2023 Jan 9]. Search in Google Scholar

[7] Desoer, C., & Wang, Y. T. (1980). On the generalized Nyquist stability criterion. IEEE Transactions on Automatic Control, 25(2), 187–196. 10.1109/TAC.1980.1102280Search in Google Scholar

[8] Ejaz, H., Alsrhani, A., Zafar, A., et al. (2020).COVID-19 and comorbidities: Deleterious impact on infected patients. Journal of Infection and Public Health, 13(12), 1833–1839. 10.1016/j.jiph.2020.07.014Search in Google Scholar PubMed PubMed Central

[9] Freedman, H. I., & Sree Hari Rao, V. (1986). Stability criteria for a system involving two time delays. SIAM Journal on Applied Mathematics, 46(4), 552–60. 10.1137/0146037Search in Google Scholar

[10] Ghosh, B., Barman, B., & Saha, M. (2023). Multiple dynamics in a delayed predator-Řprey model with asymmetric functional and numerical responses. Mathematical Methods in the Applied Sciences, 46(5), 5187–5207. 10.1002/mma.8825Search in Google Scholar

[11] Gold, M. S., Sehayek, D., Gabrielli, S., et al. (2020). COVID-19 and comorbidities: A systematic review and meta-analysis. Postgraduate Medicine, 132(8), 749–755. 10.1080/00325481.2020.1786964Search in Google Scholar PubMed

[12] Grondelle, E. V., Bruggen, S. V., Rauh, S. P., et al. (2023). The impact of the covid-19 pandemic on diabetes care: the perspective of healthcare providers across Europe. Primary Care Diabetes, 17(2), 141–147. 10.1016/j.pcd.2023.02.002Search in Google Scholar PubMed PubMed Central

[13] Hale, J. K. (1977). Theory of Functional Differential Equations, 2nd Edition, New York: Springer. 10.1007/978-1-4612-9892-2Search in Google Scholar

[14] Hertz, D., Jury, E. I., & Zeheb, E. (1984). Simplified analytic stability test for systems with commensurate time delays. IEE Proceedings - Control Theory and Applications, 131(1), 52–56. 10.1049/ip-d.1984.0008Search in Google Scholar

[15] Jeong, I. K., Yoon, K. H., & Lee, M. K. (2020). Diabetes and COVID-19: Global and regional perspectives. Diabetes Research and Clinical Practice, 166, 108303. 10.1016/j.diabres.2020.108303Search in Google Scholar PubMed PubMed Central

[16] Joshi, S. C., & Pozzilli, P. (2022). COVID-19 induced Diabetes: A novel presentation. Diabetes Research and Clinical Practice, 191, 110034. 10.1016/j.diabres.2022.110034Search in Google Scholar PubMed PubMed Central

[17] Kalmanovskii, V. B., & Nosov, V. R. (1986). Stability of Functional Differential Equations, London: Academic Press. Search in Google Scholar

[18] Khan, M. S., Samreen, M., Ozair, M., et al. (2021). Bifurcation analysis of a discrete-time compartmental model for hypertensive or diabetic patients exposed to COVID-19. The European Physical Journal, 136(8), 853. 10.1140/epjp/s13360-021-01862-6Search in Google Scholar PubMed PubMed Central

[19] Khunti, K., Aroda, V. R., Aschner, P., Chan, J. C. N., et al. (2022). The impact of the COVID-19 pandemic on diabetes services: planning for a global recovery. The Lancet Diabetes & Endocrinology, 10(12), 890–900. 10.1016/S2213-8587(22)00278-9Search in Google Scholar PubMed PubMed Central

[20] Kouidere, A., Youssoufi, L. E., Ferjouchia, H., et al. (2021). Optimal Control of Mathematical modeling of the spread of the COVID-19 pandemic with highlighting the negative impact of quarantine on diabetics people with Cost-effectiveness. Chaos, Solitons & Fractals, 145, 110777. 10.1016/j.chaos.2021.110777Search in Google Scholar PubMed PubMed Central

[21] Kretzschmar, M., & Wallinga, J. (2009). Mathematical models in infectious disease epidemiology. Modern Infectious Disease Epidemiology, 28, 209–21. 10.1007/978-0-387-93835-6_12Search in Google Scholar

[22] Lai, H., Yang, M., Sun, M., et al. (2022). Risk of incident diabetes after COVID-19 infection: A systematic review and meta-analysis. Metabolism, 137, 155330. 10.1016/j.metabol.2022.155330Search in Google Scholar PubMed PubMed Central

[23] Ng, K. Y., & Gui, M. M. (2020). COVID-19: Development of a robust mathematical model and simulation package with consideration for ageing population and time delay for control action and resusceptibility. Physica D: Nonlinear Phenomena, 411, 132599. 10.1016/j.physd.2020.132599Search in Google Scholar PubMed PubMed Central

[24] Okyere, S., & Ackora-Prah, J. (2022). A mathematical model of transmission dynamics of SARS-CoV-2 (COVID-19) with an underlying condition of diabetes. International Journal of Mathematics and Mathematical Sciences, 2022, 7984818. 10.1155/2022/7984818Search in Google Scholar

[25] Overton, C. E., Stage, H. B., Ahmad, S., et al. (2020). Using statistics and mathematical modelling to understand infectious disease outbreaks: COVID-19 as an example, Infectious Disease Modelling, 5, 409–441. 10.1016/j.idm.2020.06.008Search in Google Scholar PubMed PubMed Central

[26] Özköse, F., & Yavuz, M. (2022). Investigation of interactions between COVID-19 and diabetes with hereditary traits using real data: A case study in Turkey. Computers in Biology and Medicine, 141, 105044. 10.1016/j.compbiomed.2021.105044Search in Google Scholar PubMed

[27] Rao, P. R. S., & Kumar, M. N. (2015). A dynamic model for infectious diseases: The role of vaccination and treatment. Chaos Solitons & Fractals, 75, 34–49. 10.1016/j.chaos.2015.02.004Search in Google Scholar PubMed PubMed Central

[28] Rizvi, A. A., Kathuria, A., Mahmeed, W. A., Al-Rasadi, K., et al. (2022). Post-COVID syndrome, inflammation, and diabetes. Journal of Diabetes and its Complications, 36(11), 108336. 10.1016/j.jdiacomp.2022.108336Search in Google Scholar PubMed PubMed Central

[29] Sanyaolu, A., Okorie, C., Marinkovic., A., et al. (2020). Comorbidity and its Impact on Patients with COVID-19. SN Comprehensive Clinical Medicine, 2, 1069–1076. 10.1007/s42399-020-00363-4Search in Google Scholar PubMed PubMed Central

[30] Sathish, T., Kapoor, N., Cao, Y., Tapp, R. J., & Zimmet, P. (2021). Proportion of newly diagnosed diabetes in COVID-19 patients: A systematic review and meta-analysis. Diabetes, Obesity and Metabolism, 23(3), 870–874. 10.1111/dom.14269Search in Google Scholar PubMed PubMed Central

[31] Singh, H. P., Bhatia, S. K., Bahri, Y., & Jain, R. (2022). Optimal control strategies to combat COVID-19 transmission: A mathematical model with incubation time delay. Results in Control and Optimization, 9, 100176. 10.1016/j.rico.2022.100176Search in Google Scholar

[32] Song, H., Wang, R., Liu, S., Jin, Z., & He, D. (2022). Global stability and optimal control for a COVID-19 model with vaccination and isolation delays. Results in Physics, 42, 106011. 10.1016/j.rinp.2022.106011Search in Google Scholar PubMed PubMed Central

[33] Ssebuliba, J., Nakakawa, J. N., Ssematimba, A., & Mugisha, J. Y. T. (2022). Mathematical modelling of COVID-19 transmission dynamics in a partially comorbid community. Partial Differential Equations in Applied Mathematics, 5, 100212. 10.1016/j.padiff.2021.100212Search in Google Scholar

[34] Stephen, A. L., Kyra, H. G., Qifang, B., et al. (2020). The incubation period of coronavirus disease 2019 (COVID-19) from publicly reported confirmed cases: estimation and application. Annals of Internal Medicine, 172(9), 577–582. 10.7326/M20-0504Search in Google Scholar PubMed PubMed Central

[35] Unnikrishnan, R., & Misra, A. (2021). Diabetes and COVID19: a bidirectional relationship. Nutrition & Diabetes, 11(11), 21. 10.1038/s41387-021-00163-2Search in Google Scholar PubMed PubMed Central

[36] Wu, Y., Kang, L., Guo, Z., Liu, J., et al. (2022). Incubation period of COVID-19 caused by unique SARS-CoV-2 strains: A systematic review and meta-analysis. JAMA Network Open, 5(8), e2228008. 10.1001/jamanetworkopen.2022.28008Search in Google Scholar PubMed PubMed Central

[37] WHO Coronavirus (COVID-19) Dashboard, https://covid19.who.int/. Search in Google Scholar

Received: 2023-09-28

Revised: 2023-12-12

Accepted: 2023-12-16

Published Online: 2023-12-31

This work is licensed under the Creative Commons Attribution 4.0 International License.

Influence of incubation delays on COVID-19 transmission in diabetic and non-diabetic populations – an endemic prevalence case

Abstract

1 Introduction

2 The mathematical model

2.1 Equilibrium solution and characteristic equation

3 Local stability analysis

3.1 Delay in incubation of diabetic population ( τ 1 > 0 , τ 2 = 0 )

Lemma 3.1

Theorem 3.1

3.2 Delay in incubation of non-diabetic population ( τ 1 = 0 , τ 2 > 0 )

Theorem 3.2

3.3 Delays in incubation of both diabetic and non-diabetic populations ( τ 1 , τ 2 > 0 )

Theorem 3.3

Remark

4 Global stability analysis

4.1 Global stability analysis for endemic equilibrium point

Theorem 4.1

Proof

5 Numerical experiments

Example 1

Example 2

Example 3

Remark

6 Conclusions and future work

Acknowledgements

References

Journal and Issue

Articles in the same Issue