Skip to content
BY 4.0 license 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 EMAIL logo , Palla Danumjaya and Ponnada Raja Sekhara Rao

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.

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.
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. τ 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.

  2. τ 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.

  3. τ 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]:

  1. F 1 ( m ) = P ( m ) + Q ( m ) 0 , ( m ) 0 ,

  2. F 2 ( m ) = P ( m ) Q ( m ) 0 , ( m ) = 0 , m 0 ,

  3. 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:

  1. M 1 = ( l 3 l 4 + l 5 ) 2 l 1 ( l 3 l 4 + l 5 ) ( l 2 l 6 + l 7 ) ,

  2. 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 ) ,

  3. 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 ) ,

  4. 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:

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

  2. M 1 , N 1 , P 1 , Q 1 > 0 ,

  3. 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

  1. M 2 = ( l 3 + l 4 l 5 ) 2 l 1 ( l 3 + l 4 l 5 ) ( l 2 + l 6 l 7 ) ,

  2. 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 ) ,

  3. 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 ) ,

  4. 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:

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

  2. M 2 , N 2 , P 2 , Q 2 0 ,

  3. 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:

  1. d ( a a 1 ) y * S 1 ( b b 1 ) z * S 2 > 0 ,

  2. d 1 + ( 1 p ) α 1 + a 1 a 1 N 1 > 0 , and

  3. 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 +