Numerical analysis of a bi-modal covid-19 SITR model

This study presents a structure preserving nonstandard finite difference scheme to analyze a susceptible-infected-treatment-recovered (SITR) dynamical model of coronavirus 2019 (covid-19) with bimodal virus transmission in susceptible population. The underlying model incorporates the possible treatment measures as the emerging scenario of covid-19 vaccines. Keeping in view the fact that the real time data for covid-19 is updated at discrete time steps, we propose a new structure preserving numerical scheme for the proposed model. The proposed numerical scheme produces realistic solutions of the complex bi-modal SITR nonlinear model, converges unconditionally to steady states and reflects dynamical consistency with continuous sense of the model. The analysis of the model reveals that the model remains stable at the steady state points. The basic reproduction number Rcovid falls less than 1 when treatment rate is increased and disease will die out. On the other hand, it predicts that human population may face devastating effects of pandemic if the treatment measures are not strictly implemented.


Introduction
Novel coronavirus (covid- 19) belongs to a large class of fatal viruses that has infected millions of human across the planet and has posed serious challenges not only to individual life-styles but also to economies and GDPs of even developed countries. The infection caused by coronaviruses is identified by symptoms ranging from common cold illness to critical respiratory problems. The covid-19 was first identified as core cause of several cases of pneumonia in December 2019, when the very first respiratory infection case was recorded in Chinese city, Wuhan [1]. Generally, the coronavirus transmission takes place when an exposed individual catches infected droplets released by sneezing, exhaling or coughing of infected individuals. Normally, 80-85% of the infected individuals get recovered from the covid-19 pandemic without any special treatment [2,3]. The worldwide confirmed covid-19 infected cases have risen up to nearly 110 million, and over 2.4 million disease induced deaths have occurred. As of February 18, 2021, according to Worldometers [4] USA, India and Brazil have faced maximum death tolls amounting 0.5, 0.16 and 0.24 million respectively.
The rise in number of covid-19 cases in February 2020 led to emphasize more on modernizing the research on covid-19 predictions by exploiting latest data updates in more effective and complex covid-19 models. Like several other research issues related to covid-19 disease, the dependable forecasting of its transmission dynamics is also an essential component of the investigations. The study of epidemiological dynamics is an extensively transited research issue in mathematical modelling and simulation. There exists a variety of models that have been proposed by making specific assumptions related to modes of virus propagation. For example, susceptibleinfected-susceptible (SIS) models were used to formulate nonlinear occurrence rates in paired epidemic assumptions [5], the susceptible-infected-recovered (SIR) models employ recovery rates from the epidemic [6,7], susceptible-exposed-infectedrecovered (SEIR) type models predict the proliferation of epidemic to constitute exposed class [8], the concept of quarantined individuals is modeled through susceptible-exposedinfected-quarantined-recovered (SEIQR) integer and fractional order dynamical systems [9,10]. The fractional order modelling of dynamics of covid-19 disease were primarily explored in studies presented by Khan and Atangana [11], Khan et al. [12], Khan, Ullah & Kumar [13], Atangana & _ IG retaraz, [14] and Fehaid et al. [15]. Atangana [16] investigated the dynamics of covid-19 with new fractal-fractionalonly lockdown strategies before the emergence of any vaccine. Ullah and Khan [17] investigated the influence of non-pharmaceutical intrusions on the changing aspects of covid-19 by incorporating optimal control measures.
An important point to be noted is that a majority of existing deterministic models consist of ordinary differential equations (ODEs) by assuming an invariant diffusion process in the population. However, several other approaches exist which employ diffusion based partial differential equations (PDEs) to model the imperfect propagation [18][19][20][21]. Over and above the ODEs or PDEs based usual models, some other fruitful reports exist in literature that present individual-based [22], agent-based [23], network based [24], evolutionary computation based [25,26], stochastic [27] and statistical [28] models to investigate the proliferation of epidemics in complex scenarios of the epidemics. The reliable simulations of such models have also proved to be a basis for relevant applications [29]. The reliability of discretization based numerical methodologies for ODEs/PDEs based models is its stability, consistency and convergence. However, due to discrete nature of data related to subpopulations of the model, the finite-difference schemes preserving all essential characteristics like boundedness and/ or positivity of the state variables of the model are focused [30][31][32][33].
Biologists, mathematical analysts and actuarial scientists are also keenly following the transmission dynamics of corona virus to accelerate the diagnostic measures and preparation of effective vaccines. In this connection, a study on the symptoms of covid-19 victims in ICUs was conducted by Cao et al. [34,35]. Nesteruk [36] suggested a susceptible-infected-recovered (SIR) covid-19 model and mined some useful strategies by conducting a statistics based detailed analysis of parameters employed in the model. Later on, a revised SIR model [37] was proposed to highlight the procedure for predicting an actual number of infections that may exert additional burden on isolation units and ICUs. Zeb et al. [38] proposed an SEIIR (susceptible-exposed-infected-isolated-recovered) model to extract some meaningful values of isolation rates of infected individuals from covid-19 pandemic to isolation centers. The final volume of covid-19 pandemic was predicted by Batista [39] through a logistic growth model.
Reportedly, after an extensive research the scientists have succeeded in preparing vaccine for prevention or treatment of covid-19 infection. However, the preparation, analysis and supply of vaccine to each and every individual need a longer processing time. The dense and highly populated places of developing countries, awaiting vaccines, are vulnerable to face multiple waves of pandemic. In the absence of vaccines, other treatment measures including social distancing, quarantine, lockdowns, wearing masks/gloves and use of sanitizer were suggested by the health authorities to prevent susceptible population from becoming exposed to virus.
In this study, we assume that the covid-19 spreads in populations with constant diffusion. For this purpose, we propose a variant of a bi-modal SITR model in which the proliferation of covid-19 directly takes place from two susceptible classes and the infected individuals are admitted for treatment. The proposed model consists of coupled nonlinear ODEs involving different real parameters. The motivation behind proposing such epidemiological is the current trends of saving, dissemination and understanding of covid-19 data include the count of susceptible, infected, treated and recovered human from covid-19 in publically published almost all countries [40]. By imposing specific assumptions on the transmission modes we develop a coupled system of nonlinear ODEs. We will compute and analyze the steady states of the model along with their stability analyses regarding the computed basic reproductive number of the model. The crux of our proposed scheme is that the continuous model undergoes a non-local discretization, and we will show that the resulting discrete model preserves various important structural properties like positivity and boundedness which sort the proposed scheme as an efficient simulation tool of the SITR model. This aspect is extremely important because the state variables are sizes of subpopulations [41,42].
The remainder of present study is structured as follows. Section 2 we introduces the underlying bi-modal SITR mathematical model. The description will also include the assumptions made for derivation and physically relevant parameters to make it maximally realistic. In Section 3, we will compute the steady state points of the model. The basic reproductive number R covid will be computed as well as investigated for sensitivity and stability results for the steady state points will be established in the same section. The numerical simulations of the proposed bi-modal SITR model using an NSFD numerical scheme based on Mickens theory [43] have been presented in Section 4. The closing part of the manuscript consists of some concluding remarks and future directions.

Description of the bi-modal SITR model
The considered covid-19 dynamical model consists of five nonoverlapping compartments named as not infected susceptible, not infected but aged or serious-disease-suffering susceptible, infected, under treatment and the recovered subpopulations respectively. At any time instant t 2 ½0; 1Þ, the real valued differentiable functions S 1 t ð Þ, S 2 t ð Þ, and IðtÞ, denote the number of susceptible individuals that are not infected from the infection, the numbers of susceptible individuals that are aged or suffering from some other serious disease and seriously infected individuals from coronavirus respectively. Reportedly, the covid-19 vaccine has been prepared in technologically advanced countries and is in process for its provision to the common man. Since the latent period of the virus is small, so there is possibility of multiple covid-19 waves until each and every human is treated. However, a large proportion of the human population can be treated by implementing strict precautionary measures like social distancing, wearing masks and smart lockdowns. Therefore, a non-negative and differentiable real valued function TðtÞ is used to denote the number of individuals undergoing some treatment at any time t. Finally, the function R t ð Þ is used to denote the number of recovered individuals from the covid-19 disease as a result of using treatment measures. This work assumes that the two susceptible classes S 1 and S 2 are constantly saturated at the rates K 1 and K 2 respectively. At the same time, a constant natural death rate a is considered for all subpopulations. The constants b 1 and b 2 present the infection rates of two susceptible classes due to close contacts with the infected individuals which results in recruitment of new infected cases to the infected compartment.
The assumption made here is that the treatment provided to the infected individuals is at a constant rate l. The parameter q is used to represent the recovery rate in the result of treatment.
Using the above stated assumptions, the transmission dynamics of covid-19 in human population in the presence of treatment measures are depicted by the compartmental model exhibited in Fig. 1. The model consists of following system of five governing equations.
The system is constrained to the following system of initial conditions: The parameters presenting the transmission rates of covid-19 transmission model and the initial states are presented in Table 1.
The model (1) has following feasible region: LettingK ¼ K 1 þ K 2 the dynamics of size N t ð Þ of entire population is given by: Integrating over ½0; t we get: It shows that the system (1) possesses all solutions that fall in feasible space C. The positive invariance of the feasible region is guaranteed by its positivity and boundedness.

The basic reproductive number
The most important measure in the covid-19 pandemic is the basic reproduction number that is denoted by R covid here. If  R covid < 1 then the disease will die out and if R covid > 1 then the disease will persist in population. We calculate R covid by applying next generation matrix approach considering the compartments of infectious population and the treatment subclass [6,8,9]. Let X ¼ I; T ½ t then The Jacobeans of f X ð Þ and v X ð Þ are denoted by F andÑ respectively and are given as under.
Substituting virus free point C 0 in the product FÑ À1 and computing spectral radius of resulting matrix we obtain following value of R covid .

Sensitivity of basic reproductive number
Let us investigate the behavior of R covid due to variation in each of the model parameters and compute the relevant effect on R covid and the possible potential for control of the covid-19 disease. For this purpose we compute the normalized forward sensitivity index (NFSI) [44] of a parameter p, denoted by NFSI (p), for R covid using following relation.
The NFSI's of parameters K 1 ,K 2 , b 1 , b 2 and l are given respectively as under.
The above calculations clearly show that R covid increases when any of the parameters K 1 ,K 2 , b 1 , and b 2 increases whereas an increase in treatment rate l results in a decrease in R covid . One can observe from Fig. 2 that R covid is the most sensitive to the parameters K 1 and b 1 . The parts (a) and (b) of Fig. 3 demonstrate the behavior of R covid under variation of transmission rates b 1 and b 2 and the treatment rate l respectively.

Equilibrium analysis
There exist two equilibrium points of the model that are found by solving the equations obtained by equating time derivatives in system (1) to zero. Clearly, there exists a disease free equilibrium (DFE) point which is: The endemic equilibrium (EE) point is: with

Stability analysis
This subsection presents the stability of the proposed model at the calculated equilibrium points.
Theorem 1. The disease free equilibrium point C 0 of system (1) is locally asymptotically stable if and only if R covid < 1.
Proof. Let us express the right hand sides of governing equations by following functions F i : 1 i 5.
The Jacobean of the system (10) at C 0 is given by: Solving the associated characteristic equation we get the following eigenvalues of Þ =a À ða þ lÞ, k 4 ¼ Àa À q < 0. Now Hence all the eigenvalues of J 0 at DFE point are negative. h Theorem 2. The endemic equilibrium point C Ã is locally asymptotically stable if and only if R covid > 1.
Proof. The Jacobean of the system (10) at C Ã is given by: It can be further modified as under: The associated characteristic equation is: Clearly, there are two negative eigenvalues of J Ã that are k 1 ¼ Àa and k 2 ¼ Àa À q. Rest of the eigenvalues are the roots of following depressed equations.
Since the algebraic form of solution of the above equation is quite complicated, therefore we solve it numerically. Using the values of I Ã ; S Ã 1 ; S Ã 2 and values of parameters from Table 1 Fig . 3 Dynamics of R covid under variations of (a) b 1 and b 2 (b) l.

Simulation results
This section is devoted to analyze the simulation results of the developed model. For this purpose we construct a non-standard finite difference scheme for the underlying model. For any natural number n, we approximate the continuous derivative f 0 ðt n Þ at time step t n by the following discrete time derivative to all the state variables [45].
Representing the values of state variables at time t n by S n 1 , S n 2 , I n , T n and R n and employing Eq. (11) we get the NSFD scheme in the following form: The initial data S 0 1 , S 0 2 , I 0 , T 0 and R 0 for the scheme (12) is presented in Table 1.
We    states cannot be achieved by during coevolution of treated population T relative to S 1 and R whenever b 1 exceeds its optimally adjusted value presented in Table 1.

Conclusion
In present study, we have considered a bi-modal SITR system that models the transmission of covid-19 by dividing the normalized susceptible population into two non-overlapping subpopulations. Employing non-negative initial states, constant parameters and available covid-19 data, our proposed method considered the uninfected susceptible, un-infected but having some other serious disease or aged susceptible, infected, treated and recovered subpopulations. The steady state points of the model were computed in the absence as well as the presence of the pandemic. The basic reproductive number R covid was computed and analyzed for its sensitivity. The sensitivity analysis revealed that the proposed model can guide in drawing optimal values of the model parameters in order to eliminate the disease. The local stability of steady state points was proved by employing both of the conditions R covid < 1 and R covid > 1. Additionally, a nonlocal nonstandard finite-difference scheme based on Mickens' theory was constructed for simulating the proposed covid-19 model addressing the degree of availability of treatment measures. The simulation results favor the true dynamics of the model in the sense that the optimal combination of parameters can be achieved by implementing effective treatment measures. As a future direction the proposed bimodal SITR model can be extended to fractional order, stochastic and delayed dynamical senses.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.