Mathematical Modelling Co-existence of Diabetes and COVID-19 : Deterministic and Stochastic Approach

Diabetes plays a vital role in the epidemiological cause of COVID-19 widely. So that the study of COVID-19 morbidity along with diabetes is to be elevated for the better solutions of public healthcare policies. Even though the COVID-19 infection along with diabetes was present in the human body, COVID-19 only having the potential to infect others and diabetes can remain with the particular human so that we call it a co-existence of two diseases. The infection point of view is considered only in the COVID-19 part. Diabetes will not be transmitted from one to other. Based on these ideas, we have developed a mathematical model which consists of ﬁve compartments corresponding to ﬁve population classes, namely, diabetes susceptible, diabetes patients, COVID-19 susceptible, COVID-19 infected, and COVID-19 recovered class. We have shown the non-negativity and invariant region of the system and discussed the existence of equilibrium points. Our system exhibits two equilibrium points. The stabilities of those equilibria of the model are studied. Also, by using the next-generation matrix method the basic reproduction number ( R 0 ) is calculated. Further, we perform the sensitivity analysis of R 0 to see the efﬁcacy and inﬂuence of each of the parameters involved in R 0 and observed that reduction of transmission coefﬁcient ( α 1 ) from diabetes susceptible class to diabetes class is the most critical factor to control the co-morbidity. We have attempted to ﬁt our model with the data given by the World Health Organization(WHO) [1] and it also suits well with the data. The deterministic model is extended into the stochastic model and by using numerical simulations our results of stochastic and deterministic models are compared. Our numerical ﬁndings were illustrating the fact that the chance of getting COVID-19 by diabetes patients is high if they have come into contact with COVID-S. 19 infected people. The results also illustrate the robustness of our model from the eco-epidemiological perspective. It is also obtained and highlighted that the burden of diabetes and COVID-19 co-existence and the role of the α 1 in the severity of the disease.


Introduction
According to World Health Organization(WHO), the novel coronavirus disease 2019 or simply COVID-19 is an infectious disease caused by SARS-CoV-2 i.e Severe Acute Respiratory Syndrome CoronaVirus-2 [2,3,4,5]. Coronavirus transmission can occur from infected individuals to susceptible individuals via direct communication (contaminated hands) or through droplets spread by coughing or sneezing from an infected person or airborne spread. The symptoms of COVID-19 include fatigue, dry cough, fever, diarrhea, chest pain, loss of smell and taste senses, shortness of breath, and breathing difficulties. COVID-19 has been declared as a public health emergency on January 30, 2020 [6] and pandemic on March 11 2020 [7] by WHO, due to its exponential worldwide growth.
Mortality rates are higher in older people and patients with pre-existing medical conditions like diabetes, chronic lung disease, hypertension, malignancy, cardiovascular disease, and cerebrovascular diseases [14], [15]. Recently, epidemiological studies are focusing on the co-infection of COVID-19 and comorbidities, because comorbidities have been recognized as a risk factor for increased COVID-19 cases and higher fatality [16], [17]. Also, various studies have shown that diabetes is one of the major co-morbidity associated with severe illness and death from COVID-19 infection [18,19,20].
Diabetes is a chronic condition that occurs when there is an imbalanced level of glucose in the blood [21]. According to the report of the International Diabetes Federation (IDF), [22] and WHO [23], 463 million adults are having diabetes; nearly 760 billion USD was spent on healthcare; more than 310 million diabetes patients are living in urban areas and about 11.3% of deaths were due to diabetes.
According to the Centers for Disease Control and Prevention (CDC), [24], the rate of fatality from COVID-19 is up to 50% higher in patients with diabetes than in the individual who does not have diabetes [25]. In [26], [27], it was reported that 20-30% of patients who died from COVID-19 had diabetes. A survey conducted in China showed that the higher prevalent comorbidity in COVID-19 patients was diabetes with 8.2% [28].
However, the above attestations demonstrate the fact that co-infection of diabetes and COVID-19 have a poor prognosis and increased death rate. So, analyzing how diabetes aggravates COVID-19 outcomes can give us better insight to restrain the disease and development of effective public health policies. At a recent time, several researchers tried to study the effects of co-infection of diabetes and COVID-19 through their respective fields.
For example, in [29] the author mentioned that COVID-19 patients with diabetes are susceptible to lower extremity ischemic and other complicated diseases. Rimesh Pal and Sanjay [30] studied an unholy interaction of COVID-19 and diabetes. They found that concurrent COVID-19 makes glucose control difficult to treat due to fluctuations in blood with diabetes patients.
Furthermore, in [31] authors have presented co-infection related worries and subsequent outcomes. Here authors have suggested that improving knowledge and precaution measures are necessary for people with diabetes during the COVID-19 pandemic. Interesting research for diabetes and COVID-19 co-infection disease is carried out in [32], where authors have exhibited the fact that metformin treatment might benefit the COVID-19 patients with diabetes.
On the other hand, in mathematical modeling, there is a lack of research focus on studies dealing with co-infection of diabetes and COVID-19. Although, many mathematical models have been derived to understand the patterns of each disease (COVID-19 [33,34,35,36,37] and diabetes [38,39,40,41,42]), as of now only one study [43] has explored diabetes and COVID-19 co-infection. This paper [43] emphasizes that the negative impact of quarantine and several strategies to the lifting of the quarantine. In addition, the author worked on the optimal control techniques to provide an optimal treatment guideline for diabetes and COVID-19 co-infected patients.
Motivated by the aforementioned study, in this paper, we analyze the dynamics of the co-infection of diabetes and COVID-19. Our aim is twofold: (i) We framed a deterministic model that describes diabetes and COVID-19 coinfection. We show the existence of equilibrium and discussed its stability.
(ii) Actually, stochastic models are more realistic as compared to deterministic models. Indeed, the stochastic models are more capable of capturing randomness in epidemic disease [44,45,46,47,48]. So our proposed model is extended to a stochastic model to see the effects of stochastic perturbations.
The structure of the rest of the paper is as follows. In section 2 , we have constructed a mathematical model, which is based on some basic assumptions. The section 3 , covers the basic properties of our system. In section 4 , we have discussed the existence of the equilibrium points and in subsection 4.3 the basic reproduction number (R 0 ) is calculated. In section 5 , the sensitivity analysis of R 0 is demonstrated which helps us to find out the key parameters. In section 6 , the stability of equilibrium points is investigated. In section 7 , we tried to fit the data from WHO. In section 8 , we present our derived stochastic model extended from the determin-istic model. In section 9 , numerical simulations are provided to verify the analytical findings. Finally, in section 10 , we conclude some of the key findings and ecological implications of our research.

Model formulation
In this section we have proposed a mathematical model for co-infection of diabetes and COVID-19 disease. To derive our model equations, we divide the total population N (t) into five compartments namely Diabetes susceptible class S d (t), Diabetes class D(t), COVID-19 Susceptible class S c (t), COVID-19 Infected class I c (t) and COVID-19 Recovered class R c (t).
For developing the mathematical model, we make the following assumptions.
i. All compartments are mutually exclusive. ii. Total population is given by iii. The population of our model is homogeneously mixed. This means that each individual has an equal chance to enter into COVID-19 infected compartment. iv. A simple mass-action type interaction for the disease transmission is considered. v. Diabetes is a non-communicable disease caused by an imbalanced level of glucose. So there is no transmission of diabetes between two people. vi. If either the production rate of insulin in the pancreas is low or the body cannot use the produced insulin effectively, then diabetes susceptible individual (S d ) moves to diabetes class (D) at a rate of α 1 . vii. The COVID-19 susceptible population (S c ) is composed of individuals who have not yet been infected by the COVID-19 but may be infected through contact with COVID-19 infected people. viii. Since COVID-19 is not age-specific i.e. any adult age group can be affected if they come into contact with infected people so that we considered the transmission between S d to I c at the rate of α 2 S c to I c at the rate of β D to I c at the rate of ρ ix. γ is the rate at which infective individuals in I c recovered from the disease enter to the recovered class R c . x. The individuals are losing their immunity against COVID-19 and hence, moving to COVID-19 susceptible class. xi. The general people may stay in the recovered class itself but the diabetes patients were automatically moved to their native diabetes class. Rate of recovery of an individual due to treatment (may be with temporary immunity) δ Rate of movement of individuals from R c to D α 1 Rate of progression of individuals from S d to D α 2 Rate of disease transmission from S d to I c ρ Rate of disease transmission from D to I c Rate of movement of individuals from R c to S c due to loss of immunity µ Natural mortality rate µ 1 Mortality rate due to COVID-19 Fig. 1: Schematic diagram of the model (1) Based on these assumptions, the following mathematical model of co-infection of diabetes and COVID-19 has been formulated: with initial conditions Fig. 1 represents the flow of individuals from one class to the other. Also, the biological meanings of parameters involved in system (1) are given in Table 1 .
3 Qualitative properties of the model

Non-negativity of solutions
Theorem 1 Every solution of system (1) with initial conditions (2) exists in the in- Proof Since the right hand side of system (1) is completely continuous and locally Lipschitzian on C (space of continuous functions), the solution (S d (t), D(t), S c (t), I c (t), R c (t)) of (1) with initial conditions (2) exists and is unique on [0, ξ), where 0 < ξ ≤ +∞.
From the fourth equation of (1), we have Also, from the fifth equation of (1), we get It follows from the third equation of the system (1) that, Hence, Similarly, it follows from the first equation of the system (1) that, Finally, from the second equation of (1), we have Therefore, we can see that This completes the proof.

Invariant region
Theorem 2 The feasible region Ω defined by Proof Adding the equations of the system (1) we obtain The solution N (t) of the differential equation (3) has the following property, where χ = Λ+Φ and N (0) represents the sum of the initial values of the variables. As enters Ω or approach it asymptotically. Hence it is positively invariant under the flow induced by the system (1). Thus in Ω, the model (1) is well-posed epidemiologically and mathematically. Hence it is sufficient to study the dynamics of the model in Ω.

Existence of equilibrium
In this section, we analyze the existence of equilibrium points of the system (1). System (1) has two types of equilibrium points:

Diseases Free Equilibrium
In Diseases Free Equilibrium (DFE)

Endemic Equilibrium
At an Endemic Equilibrium (EE) , disease is present and the followings hold: Solving the equations of system (1) at the equilibrium state we get Now, putting the values of S * d , D * , S * c , R * c into the fourth equation of (1) and simplifying we obtain where Obviously A 1 is negative, However, the signs of A 2 , A 3 , A 4 are not obvious. Now, applying Descartes's rule of signs in equation (4) we obtain some of the results. These results are summarized in the Table 2 . From Table 2, it is easy to observe that if A 4 is positive, then there exists at least one positive value for I * c . i.e, at least one non-trivial endemic equilibrium. Summarizing the above discussions we arrive to the following result: (4) is positive, then there exists at least one non-trivial endemic equilibrium of the system (1).
The basic reproduction number The basic reproduction number (R 0 ) for this proposed model (1) is defined as the average number of new infective individuals generated by a single infective individual during his/her infected period when introduced into a susceptible population. One of the theoretical methods for the estimation of the R 0 is by the construction of the next-generation matrix for the system (1) as follows [49], [50]: The basic reproduction number (R 0 ) of model (1) is given by the largest eigenvalue of F V −1 , and is obtained as,

Sensitivity Analysis
The basic reproduction number (R 0 ) of system (1) depends on eight parameters, namely, recruitment rate of COVID-19 and diabetes susceptible (Λ and Φ respectively), disease progression rate from S c to I c and R c to D (β and δ respectively), transmission rate by an infective from S c to I c (α 1 ), rate of movement from S d to D due to imbalanced glucose level (α 2 ), recovery rate of COVID-19 individuals with temporary immunity (γ), natural mortality rate (µ) and mortality rate due to COVID-19 (µ 1 ). Among those parameters, we can not control some of the parameters like µ and µ 1 . The partial derivatives shown below illustrate the sensitivity of the R 0 to changes in these parameters: Therefore, to study the sensitivity of R 0 of the parameters Λ, Φ, β, γ, α 1 and α 2 , following Arriola and Hyman [51], the normalized forward sensitivity index [52] with respect to each of those parameters are calculated as follows:   From the above expressions, it is clear that Θ Λ , Θ Φ , Θ β and Θ α2 are positive and Θ γ , Θ α1 are negative. This implies an increase in the parameters Λ, Φ, β, and α 2 leads to an increase in the value of the basic reproduction number, whereas the increase in the parameters γ and α 1 causes the decrease in basic reproduction number. As we see that Θ Φ = Θ α2 which implies small changes in any of these parameters in Φ and α 2 will have the same impact on reproduction number R 0 . Similarly, Θ Λ = Θ β , we conclude that Λ and β have equal impact on basic reproduction number.
Here, Fig. 2 -7, are demonstrating the impact of the parameters on reproduction number R 0 . Λ and Φ have an directly proportional relationship with R 0 , i.e, the increase in any of them will cause an increase in R 0 . This fact is demonstrated in Fig.  2. On the other hand, the parameter Φ which is the recruitment rate of diabetes has more effect on the basic reproduction number R 0 as compared to the parameter α 2 which is shown in Fig. 3.
Likewise, in Fig. 4 shows α 2 has higher impact than α 1 on R 0 . In Fig. 5 we see the relation between R 0 , Φ and γ. As expected increase in recruitment rate Λ causes increase in basic reproduction number.
However, α 1 have an inversely proportional relationship with basic reproduction number, i.e, increase of α 1 will lead to a decrease in R 0 and if α 1 decreases R 0 will increase which is reflected in Fig. 6. Next, in Fig. 7 the value of R 0 decreases with the increment of γ. As a result of our findings, we conclude that recovery(γ) helps to eradicate the disease by reducing the R 0 .

Stability Analysis
The stability results of different equilibria are summarized in the following theorems.
Proof The variational matrix evaluated at the equilibrium point E (0) is given by . Therefore, eigenvalues of the above matrix are Here, λ 1 , λ 2 , λ 3 and λ 5 are clearly real and negative. Also as R 0 < 1, then and therefore λ 4 is also real and negative. Hence, the system (1) shows local asymptotic stability at E (0) .
Theorem 5 The disease-free equilibrium E (0) is globally asymptotically stable in R 5 if R 0 < 1.
Proof We construct a Lyapunov function L S d , D, S c , I c , R c : Substituting the equation dIc dt , we get where, The characteristic equation corresponding to this variational matrix is given by the following polynomial equation is

Data fitting
There may be some real situations to adopt a mathematical model that can be found in society. The same kind of adoption we are also trying to progress and corroborate with our results along with real data. On this aspect, we have absorbed the data of COVID-19 confirmed cases as given in [1] which is worldwide during the time period May-2020 to May-2021. Even though there may be some studies about diabetes patients having COVID-19, there is no availability of data about the number of COVID-19 patients who have diabetes. It has been predicted that 15%-30% of COVID-19 patients have diabetes in some of the literature [53], [54], [55]. A report [56] from the Washington State demonstrates the fact that 33% of COVID-19 patients had known diabetes. In Satman et. al. [57], the authors conducted a study based on the mortality rates of COVID-19 patients with and without type-2 diabetes in Istanbul, Turkey. Their analysis reveals that nearly 22% of COVID-19 patients have pre-existing diabetes. However, the findings of Chen et. al. [58] showed that 15% had diabetes among 904 patients who have been hospitalized with a clinical diagnosis of infection.
Based on the above facts and data, we are also trying to follow that and consequently, we assume that 20% of the COVID-19 affected population reported by WHO have diabetes. Here we are using the reported cases (May-2020 to May-2021) of coinfection in six countries to calibrate the model (1). These data are taken from WHO [1].    For our study, we choose America [8], India [9], Brazil [10], France [11], Russia [12] and Turkey [13] as these are highly affected countries by COVID-19. Out of all these six countries, we can see that the highest number of cases are in America and less number of cases are in Turkey. Fitted data with a model solution is shown in Fig.  8-9.

Parameter
The comparison of our model (1) and real data, during the period from May-2020 to May-2021 presented below. Fig. 8 represents a comparison of real data and our model for co-infection cases from May-2020 to May-2021. The simulated results are close to the real data in co-infection cases as depicted in Fig. 8. And in Fig. 9, fitted cumulative cases for America, India, Brazil, France, Russia, and Turkey, are shown respectively.
From the data of cumulative confirmed cases, we calculate the month-wise data. The monthly confirmed cases of co-infection worldwide during the period of study are provided in Fig. 10. From Fig. 11 we can see the scenario of monthly co-infection cases for different countries namely, America, India, Brazil, France, Russia, and Turkey, respectively.

Stochastic model
Here, we extend our deterministic model (1) to stochastic model as stochastic models are more capable of capturing random variations in biological dynamics of the problem under consideration. The derivation of the stochastic model is based on the method developed by Yuan et al. [59] . Let  Table 4 . Let us consider the case when one susceptible human becomes infected human due to the imbalanced level of glucose. In this case the state change ∆Z is denoted by ∆Z = (−1, 1, 0, 0, 0) its probability of the change is given by One can easily determine the expectation change E(∆Z) and its covariance matrix V (∆Z) associated with ∆Z by neglecting the terms higher than o(∆t). The expectation of ∆Z is given by Here it can be noted that the expectation vector and the function f are in the same form as those in deterministic model (1). Since the covariance matrix it can be approximated with diffusion matrix Ω times ∆t by neglecting the term of (∆t) 2 such that V (∆Z) ≈ E((∆Z)(∆Z) H ). Hence where the above diffusion matrix is symmetric, positive-definite and each component of this 5 × 5 diffusion matrix can be obtained as follows: V 11 = P 1 + P 2 + P 3 + P 4 = Φ + α 1 Z 1 + α 2 Z 1 Z 4 + µZ 1 , V 22 = P 2 + P 5 + P 6 + P 7 = α 1 Z 1 + δZ 5 + ρZ 2 Z 4 + µZ 2 , V 33 = P 8 + P 9 + P 10 + P 11 = Λ + Z 5 + βZ 3 Z 4 + µZ 3 , V 44 = P 3 + P 6 + P 10 + P 12 + P 13 + P 14 = βZ 3 Z 4 + ρZ 2 Z 4 + α 2 Z 1 Z 4 + γZ 4 + µZ 4 + µ 1 Z 4 , V 55 = P 5 + P 9 + P 12 + P 15 = γZ 4 + Z 5 + δZ 5 + µZ 5 , We follow the approach discussed in [59] and construct a matrix Q such that Ω = QQ H , where Q is a 5 × 12 matrix.

Numerical simulation
Analytical studies can never be definitive without numerical simulation. In this section, we have performed a numerical simulation to validate the analytical results obtained in the above sections. First, we consider the following set of parameters which corresponds to the DFE point. All parameters are in per year. Furthermore, we simulate the stochastic model (6) by using Euler-Maruyama method for the parameter set P 2 . The simulation results for both models are shown in Fig. 14 . We compare the mean of 100 runs of stochastic model simulation with the results of the corresponding deterministic model. In Fig. 15 , the time series of all the variables of the model are plotted where the deterministic simulation is represented by a black dotted curve and the stochastic model is represented by colored curves. Here, it is clear that all classes in the stochastic simulation are very close to deterministic simulation. The impact of α 1 and γ on equilibrium levels of S d , D, I c and R c , are demonstrated in the Fig. 16 -19 . In Fig. 16 variation of S d with time is shown for different values of α 1 by keeping all other parameters as in P 2 . It is observed that with the increase in the value of α 1 the equilibrium level of S d decreases.
The effect of the parameters α 1 and γ which corresponds to the transmission rate and recovery rate of diabetes population is demonstrated in Fig. 17a and 17b . From these figures, it is clear that with the increase in these parameters the diabetes population increases.   In Fig. 18a , we see that COVID-19 infected population (I c ) is directly proportional with the transmission rate (α 1 ). But I c is inversely proportional to the recovery rate of COVID-19 (γ) is demonstrated in Fig. 18b .
In Fig. 19 , the effects of COVID-19 recovered population (R c ) is shown for different pairs of α 1 and γ by keeping all other parameters constant. From this plot, it is observed that the equilibrium levels of R c increase with the increase in the rates of α 1 and γ.

Result discussion and conclusion
In this research, we have proposed a mathematical model for the co-infection of diabetes and COVID-19. Here, the population is divided into five compartments namely diabetes susceptible (S d ), diabetes (D), COVID-19 susceptible (S c ), COVID-19 infected (I c ) and COVID-19 recovered population(R c ). First, the basic properties and existence of equilibrium points of the model (1) are discussed in detail. There exist two equillibria for our model: Disease Free Equilibrium(DFE) and Endemic Equilibrium(EE).
Then we computed the basic reproduction number (R 0 ) by using the next-generation matrix method. Moreover, sensitivity analysis of R 0 is performed to find the significant parameters for the disease prevalence. It is observed that the transmission coefficient (α 1 ) from the diabetes susceptible class to diabetes class and recovery rate (γ) plays a vital role in the reduction of reproduction number R 0 . This fact is shown in Fig. 6 -7.
Further, the stability of the different equilibria of the system is analyzed by using the stability theory of ordinary differential equations. The DFE point is globally asymptotically stable whenever R 0 < 1 and EE point is locally asymptotically stable under some restriction on parameters. Apart from these, the efficiency of our system is depicted through data fitting for the real data given by WHO [1].
For our study, we choose America, India, Brazil, France, Russia, and Turkey as these are most severely affected countries in the world. In Fig. 10 -11 we can see the scenario of monthly diabetes and COVID-19 co-infection cases for worldwide and different countries namely, America, India, Brazil, France, Russia, and Turkey, respectively. And in Fig. 8 -9 , fitted cumulative cases are illustrated. Our proposed model is fitted well with these real data sets. This fact clearly exhibits that our model suggesting a good amount of decrease compared with the actual projection by WHO.
Additionally, we extend our deterministic model (1) to stochastic model (6) and the simulation results of both the models are compared. It is observed that the stochastic simulation is almost close to the simulation results of the deterministic model (Fig.  14).
Finally, to verify our analytical results, some numerical experimentation is provided. The effect of transmission coefficient (α 1 ) and recovery rate (γ) on the equilibrium levels of S d , D, I c and R c are illustrated in Fig. 16 -19. We have noted that the α 1 and γ are the key parameters and a small change in these parameters increases or decreases the COVID-19 infected population.
By observing all these results and analysis of our model the following facts have arrived.
1. Comorbidity disease increases the severity of the COVID-19 (i.e., increases the number of COVID-19 patients) and poses a significant threat to human life. 2. From the disease transmission rate and its effectiveness, we found that the progression towards COVID-19 is much higher for diabetes patients when compared with normal people. 3. By introducing our recovery rate due to COVID-19 medication reduces/eradicates the COVID-19 disease alone from individuals but diabetes patients will remain with diabetes as usual and continue the treatment/medication for diabetes.
In the future, we will try to incorporate the effects of delay, genetics, and treatment for COIVD-19 and diabetes.

Data availability
The data set we have used are just parameters and initial conditions which we are mentioned in numerical simulation. For data fitting, we have used the data available in [1]. There is no special dataset other than that we have used.
Author's contributions S. Anusha : Conceptualization of the idea, Development of the model framework, Analysis and simulation, Writing. S. Athithan: Development of the model framework, Supervision, Validation and Submission process.

Funding
There is no funding available for this article.