Mathematical model of COVID-19 with comorbidity and controlling using non-pharmaceutical interventions and vaccination.

Pandemic is an unprecedented public health situation, especially for human beings with comorbidity. Vaccination and non-pharmaceutical interventions only remain extensive measures carrying a significant socioeconomic impact to defeating pandemic. Here, we formulate a mathematical model with comorbidity to study the transmission dynamics as well as an optimal control-based framework to diminish COVID-19. This encompasses modeling the dynamics of invaded population, parameter estimation of the model, study of qualitative dynamics, and optimal control problem for non-pharmaceutical interventions (NPIs) and vaccination events such that the cost of the combined measure is minimized. The investigation reveals that disease persists with the increase in exposed individuals having comorbidity in society. The extensive computational efforts show that mean fluctuations in the force of infection increase with corresponding entropy. This is a piece of evidence that the outbreak has reached a significant portion of the population. However, optimal control strategies with combined measures provide an assurance of effectively protecting our population from COVID-19 by minimizing social and economic costs.

Introduction

The World health organization (WHO) declared COVID-19 as a worldwide health hazard at beginning of 2020. The first scenario was observed in China (Wuhan) [1]. In this connection, COVID-19 has been transmitted in whole world. In order to conquer the pandemic, various countries took a measure of a lockdown scenario as well as successively unlock process for political impacts and affecting to socioeconomic condition [2-5]. Further, the mortality and morbidity rate has been different across different countries in the world [6, 7].

COVID-19 is transmitted by people through their respiratory droplets from direct contact with contaminated objects and infected people [8, 9]. However, vaccination is in the beginning process to defeat COVID-19. Different studies can be available in stages of comorbidity (like heart disease, diabetes, lung disease, etc.) within the infected human being. All human beings with comorbidity or medical illness are at high risk to get an infection than normal or healthy people [10-12]. From the statistical data, this is reported that the patients with COVID-19 having comorbidity like kidney disease and type-2 diabetes were transferred to hospital with intensive or critical care unit [13-15].

As the number of COVID-19 follows an increasing–decreasing trend, entropy leads us to explain the heterogeneity in diffusion or force of infection [16]. Different methods like structural entropy [17] and evolutionary entropy [18] are used to study the dynamics of infectious disease, and to forecast the pandemic [19, 20]. Here, the basis of Shannon entropy is usually applied and considered for the spread of pandemic due to significant analogy which follows Boltzmann’s classical thermodynamics [21].

In the COVID-19 pandemic, various countries have implemented different non-pharmaceutical interventions (NPIs) like wearing masks, social distancing and proper sanitation, etc., to mitigate disease spread. This process helps to slow but not to stop the spread. The literature can be found to study the optimal control strategies [22-27]. Different mathematical models of epidemiology have been investigated with their corresponding transmission dynamics [28, 29]. Lie et al. [34] studied distributed delay-based COVID-19 model. In COVID-19 model, Khyar et al. [35] investigated the global dynamics of multi-strain. Moreover, various authors [36, 37] have applied optimal control theory for mitigation of COVID-19. Different authors studied the vaccination process in COVID-19, and incorporating NPIs in the vaccination process has not been yet studied according to our knowledge [38, 39]. Upadhyay et al. [40] have analyzed the effectiveness of social distancing measures which has been implemented in the early stages of infection. Moreover, different model formulations have been done in COVID-19 [41-44].

However, in order to defeat COVID-19 like a pandemic, suitable vaccination strategies are essential to be implemented with NPIs. For this purpose, combined fundamental strategies are inevitably required under optimal control problem to defeat COVID-19 like, an infectious disease with a minimum cost of the vaccination process. Our main aim is here to investigate the mathematical model of COVID-19 based on bifurcation theory and also to address the optimal control strategy to defeat COVID-19 under NPIs and vaccination [45-47]. Data-driven modeling is an essential tool for investigating transmission dynamics in the landscape of global health which is employed in our study [48-50]. Further, incorporating comorbidity in the COVID-19 model is one of ongoing research trends [51]. Moreover, it can be noted that the transmission dynamic of COVID-19 with comorbidity and its optimal control strategies under the combination of NPIs and vaccination campaign have not been yet explored to the best of our knowledge.

The subsequent part of the paper are as follows: In Sect. 2, a mathematical model of COVID-19 with comorbidity is proposed and explored. The qualitative dynamics are studied by employing bifurcation analysis in Sect. 3. In Sect. 4, optimal control problem is formulated with non-pharmaceutical interventions (NPIs) and vaccination. In Section 5, the model curve fitting and also sensitivity analysis are performed. Moreover, the trend in force of infection is investigated by Shannon entropy. Moreover, numerical simulations are present for optimal control strategies. Finally, in Sect. 6, we present our conclusion from our study.

The model

We formulate here SEICR model by introducing comorbidity in infected individuals for the transmission dynamics of COVID-19. The mathematical model contains five compartments, susceptible individuals (S), exposed individuals (E), infected without comorbidity individuals (I), infected with comorbidity individuals (C), and recovered individuals (R). It is assumed that . We investigate COVID-19 model with and without comorbidity. Based on biological viewpoint, the nonlinear mathematical model of COVID-19 in a short time window is:The initial values are assumed in model (1)Here, it is considered that , where represents the primary date of disease spreading of COVID-19 model in (1). The external recruitment increases susceptible individuals through by immigration or birth to the society. The number of the individuals reduces after getting an infection with a rate , through from susceptible individuals, infected with and without comorbidity. These interaction produces exposed individual. It is also considered that the recovery rate for infected with and without comorbidities is and , respectively. The model assumes a normal death rate for each subpopulations. Here, the period of disease incubation is , and exposed individuals become infected with and without comorbidity with a fraction and , respectively. The scenario of the model is given in Fig. 1. The description of parameters with values is provided in Table 1.

Fig. 1

Pictorial scenario of SEICR model. This shows various interactions: susceptible (S), exposed (E), infected without comorbidity (I), infected with commodity (C), and recovered (R). The inward and outward arrows represent the increasing and decreasing of population, respectively

Table 1

The values of the parameters used in the SEICR model (1)

Parameter	Description	Value	Reference
Π_s=μ×N	Average recruitment rate	4×10^4	Assumed
β_s	Transmission rate	0.4945	Estimated
1μ	Average life expectancy at birth	70.4 years	[55]
1α_e	COVID-19 incubation period	5.2 days	[56]
ξ	Fraction of exposed individuals	0.21	Estimated
γ_i	Recovery rate of infected with no comorbidity	0.1245	Estimated
γ_c	Recovery rate of infected with comorbidity	0.1241	Estimated

SEICR model analysis

Now, we investigate a few basic qualitative properties of the SEICR model (1) at biologically feasible equilibria with . We have assumed that following lemmas hold.

Lemma 3.1

The solutions of SEICR model (1) with initial values (2) lie in for all .

Proof

In order to proof the positivity of (1), we consider that any solution in remains positive for all In the system (1), we see that Here, this is the positivity of solutions in the region . is established as positive invariant set of SEICR model.

Lemma 3.2

The solutions of SEICR model (1) with initial values (2) are uniformly bounded within the region

In order to prove boundedness, we sum up all equations in the model (1), which provides Differentiating both sides, we getwhich implies Without loss of generality, we can express all equations in model (1) as . So, a bounded set can be definedwhich is further a positive invariant set of the SEICR model (1). Therefore, all solution trajectories initiating from interior of always remain within the domain . This assures that the growth of all individuals cannot be unbounded or exponential for the time window.

In the following, numerical results show that the solutions of model (1) are always positive and bounded by changing the various parameter values.

Basic reproduction number

The basic reproduction number is a dimensionless number and measures the expected values of decreasing or increasing outbreak. The infection-free equilibrium is considered as The components E, I, C in system (1) are explicitly appended with a disease transmission. The matrices , represent respective new infection and transition, given byThe dominant eigenvalue of matrix represents basic reproduction number, where , and . So, we obtain

Stability analysis of Infection free equilibrium

Theorem 1

The equilibrium becomes unstable and locally asymptotically stable for

Proof

We get the matrix J of (1) asFor our convenience, we considerThe two eigenvalues of matrix at are . The roots of following equation provide the rest of eigenvalues:which can be expressed asDenoteNow put , if , thenConsequently, , which gives . For , we get negative real parts in equation . Thus, all eigenvalues are negative for . Hence, becomes local asymptotically stable for Again, for ,then, for . The matrix contains . Thus, for , becomes unstable.

Theorem 2

The infection-free equilibrium is globally asymptotically stable for the system (1) if in the bounded region

We can express the model (1) aswhere , shows the number of susceptible or uninfected individual compartments, shows infected individual compartments. is infection-free equilibrium of the model (1). Now, G(X, U), and H(X, U) are given as follows:Now, it can be noted that In order to show global asymptotically stability, the following conditions must hold. From the condition (A), the compartments of the model (1) can be expressed asSolving the above system analytically, we get and Now, and as Thus, is globally asymptotically stable for . Hence, the condition (A) holds.

For , is globally asymptotically stable,

for where is a Metzler matrix (where non-diagonal components are non-negative) in the region

Now the matrix V and for the model (1) areIt can be noted that for the region Therefore, infection-free equilibrium of the system (1) is globally asymptotically stable in region for

Existence of infection present or endemic equilibrium

Now, we consider the force of infection at endemic equilibrium (say) asBy simplifying the (1) at equilibrium state, we obtain , , and .

Replacing into (5), we get the equation in whereSince , it can be noted that the model (1) has a unique endemic equilibrium point for and no positive endemic equilibrium for

Now, we study the local stability of the unique endemic equilibrium .

Bifurcation analysis

Theorem 3

The SEICR model (1) is locally asymptotically stable around the endemic equilibrium for . Moreover, the system undergoes transcritical bifurcation at

For , assuming as bifurcation parameter, center manifold theorem is used for stability analysis at .

Now, the right eigenvector having zero eigenvalue to the variational matrix of system (1) at is given by , where In similar way, the left eigenvector having zero eigenvalue to variational matrix at is given by , whereWe incorporate a few symbols for SEICR model system as follows: and where Now, we calculate at and getThe rest of derivatives at becomes zero. Again, we determine the coefficient a and b based on well-established Theorem 4.1 in Castillo-Chavez et al. [52]Now, we replace all to find the coefficient a and b at threshold , we getandHere, the values of a and b are negative and positive, respectively. The system (1) undergoes forward bifurcation at . The endemic equilibrium is locally asymptotically stable for

Optimal control problem

Control measures play a significant role in the mitigation of COVID-19 transmission. Now, we focus on optimal control strategies for employing non-pharmaceutical interventions along with the vaccination process in India. As the vaccination process is applied to categorized individuals like front-liner workers: doctors, nurses, clinicians, polices, etc., a suitable implementation of non-pharmaceutical interventions is thus essential to flatten the disease transmission. An optimal control approach can help to minimize the combined cost of control [53]. The optimal control problem is given by introducing time-varying controls u(t), which represents the implementation of non-pharmaceutical interventions, and , which represents a fraction of the vaccination process.with satisfying the initial condition (2). Here, our intention is to mitigate the number of infected individual with and without comorbidity in the presence of non-pharmaceutical interventions with ongoing vaccination process. We construct the objective functional as:subject to SEICR model (9). Here, and represent the weights to balance the factors. and measure the relative cost of strategic interventions and vaccination over time are considered for nonlinearity due to implication of severe strategy. According to our purpose, we only study on a finite time period . The optimal control is defined asHere is admissible control class defined on i.e.,The control pair is bounded between 0 and 1. Here, the controls become zero implying no extra measures and vaccination are implemented. On the other hand, the controls become the maximum value; one implying interventions and vaccination are . Indeed, perfectly implementation is impossible in reality and thus, we consider .

Proposition 4.1

The solutions of the optimal control (9) exist and bounded for a finite time period .

Now, Pontryagin’s maximum principle [54] is applied to evaluate the necessary conditions for optimality . In order to find the optimality in the model (9), we have constructed the augmented Hamiltonian function as the cost functionalFor optimality, we drive the adjoint equations as follows:where adjoint equations satisfy transversality conditions . For optimality conditions, the Hamiltonian is to be optimized with respect to control pair . We haveThe corresponding control pair is given as follows:

Numerical simulation and results validation

Now, we validate the analytical findings. To do this, we fit the parameter values of the SECIR model.

Model curve-fitting

We consider April 25, 2020, to January 31, 2021, for model curve fitting. For this purpose, we take daily cases of COVID-19 in India [57]. We fit the model (1) for daily cases of COVID-19 in India. We enlist the SEICR model parameters in Table 1 by estimating from the data. By fitting the model to the basis of daily report, four parameters of the SEICR model have been fitted, which are (a) transmission rate , (b) fraction of exposed individuals to become infected, (c) recovery rate of infected with no comorbidity , (d) recovery rate of infected with no comorbidity (). The initial values of the model (1) are assumed as . In MATLAB, inbuilt solver is applied to fit daily cases based on the nonlinear least square. In Table 1, the estimated parameters are given by the above technique. The daily new cases of India are fitted and shown in Fig. 2. The value of is computed by taking the parameter values of the model from Table 1.

Fig. 2

The fitted SEICR model (1) to confirmed daily data in India. Reported confirmed cases are seen in black dots, and the red line shows the fitted line

To identify the influential parameters to infected with and without comorbidity individual, we next execute sensitivity analysis.

Sensitivity analysis

We perform uncertainty analysis by Latin hypercube sampling (LHS) method as well as sensitivity analysis by Partial rank correlation coefficient (PRCC) method. The whole technique is available in Marino et al. [58]. The sensitivity of the parameters, namely , , , , , and , is based on uniform distribution with confidence domain . We assign baseline values by these parameter values. The outputs values and scenarios are reported in Table 2 and seen by bar diagrams in Fig. 3a, b. From Fig. 3a, it is observed that , , and are sensitive parameters to the infected individuals without comorbidity. On the other hand, Fig. 3b shows that , , , and are sensitive parameters to the infected individual with comorbidity. We further investigate the sensitivity indices corresponding to . From Fig. 3c, it can be seen that and are negatively correlated, and , and are positively correlated to .

Table 2

Sensitivity of the parameters of SEICR model (1) to I, C and , where i= 100, 150, 200th day

Description	β_s	μ	α_e	ξ	γ_c	γ_i
I_100	0.8124	- 0.1042	- 0.9018	- 0.1251	- 0.4871	0.1851
I_150	0.8751	- 0.4167	- 0.8501	- 0.0916	- 0.7192	0.1473
I_200	0.7110	- 0.3876	- 0.7912	- 0.0981	- 0.6875	0.0916
C_100	0.4619	0.2179	0.7901	0.9128	- 0.1927	- 0.5012
C_150	0.6972	0.1986	0.2892	0.8715	- 0.5102	- 0.6179
C_200	0.7125	0.1275	0.1917	0.2571	- 0.6171	0.7251
R_0	0.8912	0.1779	0.7126	- 0.5267	- 0.9512	0.4576

Fig. 3

PRCC to a infected individual without comorbidity (I), b infected individual with comorbidity and c basic reproduction number (). PRCC values of different parameters with significance level 0.05. LHS approach and uniform distribution with size of sample (500) are considered

Basic reproduction number measures the expectation of variation in the epidemic evolution, and we investigate the influence of parameter variation on .

Effects of parametric variation

The impact of parameter variations is further investigated under , , in Fig. 4. It is also seen that only the increasing value of can shift to in Fig. 4a. Further, the simultaneously increasing value increase the value of in Fig. 4b, c.

Fig. 4

Matrix plots showing the changing nature in basic reproduction number () of SEICR model under parametric variations: a vs , b vs , c vs

Bifurcation diagram

In order to investigate, the existence of endemic equilibrium in the SEICR model (1), forward bifurcation occurs with , given in Fig. 5a, where unique positive endemic equilibrium exists for . From an epidemiological perspective, the scenario of forward bifurcation indicates that the disease persists. It indicates that susceptible individuals must have the awareness to the disease spreading, like social distancing, using a mask, maintaining proper sanitation, etc.

Fig. 5

a vs plot indicating forward bifurcation of SEICR model with . b Effects of variation in on forward bifurcation with : the figure shows that the boundary of forward bifurcation region increases gradually with the increasing . Keeping other parameters value fixed as in Table 1. Here, DFE represents disease (infection)-free equilibrium

Moreover, it can be observed numerically that a forward bifurcation region increases gradually with increases of the value of exposed individuals gets infected, with , given in Fig. 5b. Indeed, the region of the curve is essential to be smaller to reduce the infection like the value needs near to one.

As infection intensity varies, the average of force of infection is considered for our study.

Effects on force of infection by varying parameter

Now, the variation in force of infection is studied by taking , where N represents the length of series . Figure 6a–c shows the fluctuations of mean force of infection under , and , respectively. Figure 6a shows respective increasing trend, and Fig. 6b, c shows decreasing trend in . Moreover, it can be seen that mean force of infection increases with the increase of , average force of infection decreases with the increase of and in the model (1).

Fig. 6

a–c represent vs plot (with ), vs plot (with ) and vs plot (with ). d–f: represent box plot where five measures are lower adjacent (LA), first quartile , median (M), third quartile , upper adjacent (UA) and outlier (OL). Red (+) sign indicates outlier events leading to least predictable

We further study the distribution of the force of infection by box plot in Fig. 6d–f. The box plot usually provides more information on dispersion or variability of the force of infection. It shows an impression about the disease transmission of COVID-19. This contains five measures, namely lower adjacent (LA), first quartile , median (M), third quartile , upper adjacent (UA), and outlier (OL). Red (+) sign indicates outlier events leading to least predictable.

Further, we study the combined effect under parametric plane and on ; the corresponding matrix plots are given in Fig. 7a, b. In Fig. 7a, b, it can be noted that the mean force of infection depends on and for any value of , mentioning by various color regions of increasing value of .

Fig. 7

a, b represent over matrix plot, where and over matrix plot, where . The color bar represents values of . The other parameter’s values remained the same as given in Table 1

As force of infection has disorder trend, Shannon entropy [21] can be applied to measure the disorder due to non-uniform distribution of disease spread.

Entropy within force of infection

Shannon entropy (H) iswhere N represents the event’s length and is the probability of the event , i.e., occurrence-based frequency on non-recurrent-based event . The disorder of on , and are seen in Fig. 8a–c. Figure 8a shows the increasing trend of with increase of . decrease with and in Fig. 8b, c. From Figs. 6a–c and 8a–c, it can be mentioned that the trend of is positively correlated with the variation in force of infection . Further, we also study the combined dependency of on . In the similar way, comparing Figs. 7a, b and 8d, e, we see the alike pattern of enhancing entropy of and . This assures the dynamically disorder in force of infection in (1) for the same parametric plane.

Fig. 8

Variation of by varying a , b , c d and e . The color bar represents values of . The rest of parameters are fixed as shown in Table 1

Optimal strategies over time

Now, we study the effect of combined measures under non-pharmaceutical interventions(NPIs) and vaccination under optimal control strategy. We consider to observe the control scenarios in Fig. 9. Figure 9a shows that susceptible individual becomes vaccinated under NPIs near bout 100 days. As whole population goes through NPIs with vaccination event, the fraction of exposed individual decreases gradually with increase of vaccinated susceptible in Fig. 9b. In similar way, it can be observed that infected individual with and without comorbidity acquire less infection due to the same reason during vaccination process in Fig. 9c, d. Figure 9e shows that the number of recovery individuals increases with the decreasing number of the susceptible individuals along with the infected individuals under combined measures. How non-pharmaceutical interventions as well as vaccination event can be implemented, seen in Fig. 9f, g. Further, it can be seen that although NPIs can be imposed lightly after few days, indeed needs to be continued till almost all individual becomes vaccinated. Figure 9h assures that a combined control strategy decreases the cost of intervention as well as vaccination event.

Fig. 9

Dynamics of optimality system (9) resulting from combination of NPIs and vaccination process. The individual counts-based trajectories of susceptible (S), exposed (E), infected without (I) and with (C) comorbidity and recovered (R) are in (a)–(e). The control trajectories of combined strategies are in (f) and (g) by the variation of v and u, respectively. The trajectory of cost functional is in (h). Here, amount of vaccine and weights are

Conclusion

In this article, we have formulated and studied the disease dynamics of COVID-19 through a mathematical model considering infected subpopulations as with and without commodities due to the complexity of infection transmission. The parameter of the model has been fitted to new daily COVID-19 cases of India. The qualitative dynamics have been investigated, and the basic reproduction number has been derived by using the next-generation matrix method. The model is asymptotically stable at infection-free equilibrium for . Based on the center manifold theorem, the model experiences forward bifurcation due to an increase of fraction in the exposed individuals having commodity. This assures that infection presents in society for greater than unity for , basic reproduction number. From an epidemiological perspective, comorbidity individuals get gradually infection due to lack of precautions and surveillance and like social distancing, proper sanitation wearing masks, etc. In this situations, susceptible individuals become infected and turned to exposed individuals. Indeed, exposed individuals can prevent COVID-19 infection due to strong immunity in this connection.

On the study of transmission dynamics, the average fluctuation in the force of infection has an evidence of increasing or decreasing trend with transmission rate (), the recovery rate of infected without and with comorbidity , , respectively. Further, an increasing trend to the mean force of infection has been indicated through the composite effect. Higher Shannon entropy production, i.e., more disorder in mean force of infection, indicated more strengthening the force of infection according to dynamical perspective. This might cause a dangerous situation in the population. Finally, optimal control strategies have shown that a combination of non-pharmaceutical interventions and vaccination events can be effective measures to diminish COVID-19 by minimizing the social and economic cost.

Our study reveals that the production of higher entropy means more disorder in disease transmission which can be helpful evidence to understand the state of the force of infection to computational biologists. Moreover, proposed strategy in controlling COVID-19 through proper implementations of NPIs and vaccination events can be fruitful to our society. Moreover, a study of pulsed vaccination strategies remains as our future research scope.

Table 3

Different measures of box plot corresponding to , and

Box Plot	LA	Q_1	M	Q_3	UA	OL
β_s	0.5×10^-3	1.9×10^-3	2.7×10^-3	5.3×10^-3	7.3×10^-3	–
γ_i	3.1×10^-3	3.3×10^-3	3.7×10^-3	4.4×10^-3	5.8×10^-3	5.9×10^-3
γ_c	2.5×10^-3	3.2×10^-3	5.1×10^-3	7.9×10^-3	13.5×10^-3	–