Modelling fractional-order dynamics of COVID-19 with environmental transmission and vaccination: A case study of Indonesia

1.   Introduction

Coronavirus disease (COVID-19) is a fatal illness caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), a novel coronavirus ^[1]. The disease first originated from the Wuhan City, Hubei Province in China in 2019 and has spread throughout every part of the world. As a result of no cure to this epidemic, on 31 December 2019, World Health Organisation (WHO) received a report of the disease. On 30 January 2020, WHO declared the COVID-19 outbreak a global health emergency ^[1]. The disease was declared a global pandemic by the WHO on March 11, 2020 ^[2]. In Indonesia, the reported COVID-19 of confirmed cases on June 1, 2021, is 1, 826, 527 with 50723 deaths, 101325 active cases, and 1, 674, 479 recovered cases ^[3].

The COVID-19 disease spreads from person to person through infected air droplets due to sneezing or coughing, can also be transmitted through human contacts. Initially, there was no vaccine or treatment for SAR-CoV-2. Still, some preventive measures such as hand washing, lock-downs, quarantine, use of face masks, sanitizers are used by different countries to reduce the rate of disease transmission. Most countries implemented the closure of schools and air travel to help reduce the importation of new cases, which primarily drives/leads to local transmission of the virus in the susceptible population.

Vaccination programs play an essential role in controlling the spread of infectious diseases, including COVID-19. The first mass COVID-19 vaccination program began in early December 2020. According to WHO, 13 different vaccines have been administered ^[4]. The implementation of COVID-19 vaccines such as AstraZeneca, Moderna, Pfizer, Johnson, and Johnson has helped reduce the severity of the illness once susceptible humans contract the virus.

Numerous researchers have used mathematical models that use deterministic, stochastic, and fractional order to explain the COVID-19 transmission. These models explain various strategies which can be used to control or mitigate the spread of the disease in different countries by incorporating lockdown measures, the pathogen in the environment, optimal control, non-pharmaceutical interventions, cost-effectiveness analysis, role of the media campaign, treatment and vaccination (see, for instance, ^{[5,6,7,8,9,10]}). The authors in ^[8] studied the impact of early detection and vaccination strategy in the COVID-19 eradication program in Jakarta, Indonesia, by using an $SVEIQAR$ deterministic model. Results from sensitivity analysis of their model show that rapid testing is crucial in reducing the basic reproduction number when COVID-19 is endemic in the population rather than contact tracing. Further, their results indicated that vaccination strategy could relax social distancing rules while maintaining the basic reproduction number at the minimum level and eradicating COVID-19 from the population with a higher vaccination rate. The authors in ^[10], used an $SIQRD$ model for non-age and age structures to assess the impact of COVID-19 vaccination in several provinces in Indonesia. Their study suggested that proper vaccine implementation is in the early stage of the pandemic will help suppress the number of active cases and reduce the total number of deaths. On the other hand, vaccines implementation should be consistent and not only for one or two months as this will not reduce the number of infected persons. In contrast, the prioritisation of the active and older adults (above 50) for vaccination will significantly reduce total Indonesia COVID-19 deaths.

To capture the memory effect and the non-local in nature, many authors have utilized fractional models to study infectious diseases such as Syphilis ^[11], Dengue fever ^[12] and COVID-19 ^{[13,14,15,16,17,18,19,20,21,22,23,24,25]}. For more on fractional models on COVID-19 infection interested readers is referred to the following papers ^{[26,27,28,29,30]}. The study in ^[26] employed a set of six fractional differential equations to model the transmission of COVID-19 by considering the disease transmission through the dead individuals and the pathogen in the environment for Saudi Arabia. Analyses of the fractional model were presented, as well as numerical simulation carried out in their paper. Their results suggest that a decrease in the pathogen viral load due to symptomatic and asymptomatic COVID-19 infected individuals and asymptomatic infected individuals' viral contribution in the environment decreases the disease burden significantly.

In the present study, we investigate the potential impact of the environmental pathogen which drives the COVID-19 spread in the human population, thus causing the resurgent of the infection. Vaccination is currently used to curb the spread of the disease to some extent. Still, due to the mutation of the original virus, which is more harmful, a fraction of those vaccinated can also become exposed and may contract the disease at a low rate. Thus, we introduce a modification parameter to ascertain this phenomenon. We note that the environment plays a vital role in the COVID-19 transmission, as it can be transported in the air or exist in our physical environment. Hence, we extend an integer and fractional epidemic model of COVID-19 for Indonesia by incorporating vaccination and environmental transmission of the pathogen.

The outline of the manuscript is as follows: In Section 2, we give mathematical preliminaries. Section 3 presents a mathematical model that describes the spread of COVID-19 in Indonesia with the incorporation of vaccination and environmental pathogen. Formulation of fractional COVID-19 model is given in Section 4 and the model analyses is presented in Section 5. Section 6 presents the numerical simulations using the Caputo operator. The paper is concluded in Section 7 discussing the results obtained from the study.

2.   Fractional calculus preliminaries

Here, we give a few mathematical definitions for fractional derivatives, which shall be used to study the COVID-19 model with vaccination.

Definition 1. Let $\chi\in\mathcal{C}^{k}$ be a function, then the classical Caputo derivative with a fractional order $\tau\in(k-1, k)$ for $k\in\mathbb{N}$ is thus defined as follows ^[31]:

Thus, we see clearly that $^\mathcal{C}\mathcal{D}_{t}^{\tau}(\chi(t))\to\chi(t)^{\prime}$ if $\tau\to 1.$

Definition 2. The integral operator for Definition 1 with a fractional order $\tau\ge0$ for the function $\chi:\mathbb{R}\to\mathbb{B}$ is given by

for $0 < \tau < 1, \; \; t > 0.$

Next, we define the equilibrium for any Caputo fractional-order ordinary differential model. Hence, it is give as follows:

Definition 3. If ($\chi^{*}$) is the steady state, then ^[32]

for $\tau\in(0, 1)$ if and only if $f(t, \chi^{*}) = 0.$

3.   Model formulation

In this section, we formulate an $SVEIRP$ model, which consists of the human population and the pathogens in the environment. The human compartments consist of the susceptible ($S(t)$), vaccinated ($V(t)$), exposed ($E(t)$), infectious ($I(t)$), and recovered individuals at any time $t$ ($R(t)$). While the pathogen/virus in the environment at any time $t$ is denoted by $P(t).$ Therefore, the total human population over time is given by $N(t) = S(t)+V(t)+E(t)+I(t)+R(t).$

Humans are recruited into the susceptible class at a rate $\pi.$ These susceptible humans are exposed to the COVID-19 disease at a rate of force of infection denoted by

where $\beta_{1}$ and $\beta_{2}$ are the effective contact rate for the infectious individuals and pathogens in the environment, respectively. The exposed individuals can become infected at a rate $\gamma.$ The infectious individuals can then recover at a rate $\sigma.$ In this model, we have assumed that all human compartments have a natural mortality rate $\mu$, and the infectious individuals assume an induced death rate $\delta.$ Further, the susceptible individuals can be vaccinated at a rate $\gamma$ while those in the vaccination compartment can contract the infectious after being exposed at rate $\lambda\omega.$ Here $\omega$ is the modification parameter which is much much less than the rate of the initial infection of COVID-19 infection. On the other hand, we have assumed that humans shade the pathogen into the environment at a per capita rate $\xi.$ The model parameters are summarised in Table 1 below and the schematic diagram describing the interactions between the compartments given in Figure 1.

Combining the model assumption and descriptions with parameters in Table 1, we have the following system.

subject to the initial conditions

Parameter estimation and model fitting

In this subsection, we estimate the biological parameters of the COVID-19 model (3.1) using the reported active cases (undergoing treatment) of COVID-19 in Indonesia for the given period starting from 1 June 2021 to 20 July 2021 ^[33]. We utilize the well-known least square fitting technique to estimate the parameters. The initial value of the total population of Indonesia based on the data is $N(0) = 270, 200, 000$ ^[34] and the average life span is $\dfrac{1}{71.47}$ years ^[35]. The initial populations of vaccinated, exposed, infectious, and recovered as stated in ^[3]. The data are given by $V(0) = 10714300, E(0) = 61108, I(0) = 101325$, and $R(0) = 1674479$. The initial of the susceptible population is taken as $S(0) = 257, 648, 788.$ We assume the initial population of $P(0)$ is given by $P(0) = 100000$. The estimated and fitted parameters are displayed in Table 2. The fitting result to the reported data through our model is compared in Figure 2. Hence, it can be seen that the simulation model is quite in accordance with the actual data.

4.   The fractional model

We transform the integer model (3.1) into a fractional differential model such that yields

Using the Caputo operator, we convert system (3.1) to the following system of fractional order differential equation

subject to initial conditions as in (3.2).

5.   Model analyses

This model analyses assumes that the recovered individuals do not contribute to COVID-19 infections and transmissions. Next, we establish some basic properties out model (4.2).

5.1.   The model basic properties

We prove that the model system (4.2) is positive, bounded, and biologically meaningful over the modelling time for $t\geq0.$

Proof. Adding the total change in the human population gives

The solutions of (5.1) yields

We have that $N(t)\leq \dfrac{\pi}{\mu}$ as $t\to\infty.$ Thus, the COVID-19 model feasible region defined as $\bar{\Omega}$ which consists of the human population pathogen in the environment is given by

Hence, the COVID-19 model (4.2) which incorporates vaccination is biological feasible and bounded for all time $t\geq0,$ which completes the proof.

5.2.   Existence and uniqueness of solutions of the fractional model

For simplicity in our mathematical analyses we will set in (4.2) as follows:

Here, we investigate and prove that the solutions to Eq (4.2) exist and are unique over the modelling time using the fixed point theory.

Let $\mathcal{B}(\mathcal{G})$ denote the Banach space which consists of a real valued continuous function over the interval $\mathcal{G} = [0, d]$ whose norm are defined as follows for our model compartments:

and $||P|| = \sup_{t\in \mathcal{G}}|P(t)|.$ Further, we also define the norm

Proof. We convert model (4.2) using the Caputo integral and obtain

Applying Definition 2, we have that

with the following kernels:

The expressions $\mathscr{K}_{i}$ in Eq (5.3) for $i = 1, 2, \cdots, 5$ satisfies the Lipschitz for the model respective compartments $(S(t), V(t), E(t), I(t), P(t))$ have an upper bound. Firstly, we present the proof for the susceptible ($S$) compartment. Let $S$ and $\bar{S}$ be two different functions, we thus have that

Setting $\ell_{1} = (||\alpha_{2}V+(\lambda+Q_{0})||)$, we get

Similarly for the remaining state variables, we have that

where $\ell_{1}, \ell_{2}, \ell_{3}, \ell_{4}$ and $\ell_{5}$ represents the Lipschitz constants for each kernels $\mathscr{K}_{i}$ for $i = 1, 2, \cdots, 5$ respectively. Thus, the Lipschitz condition is satisfied. To complete the uniqueness proof by applying the Banach theory, we also show that system (5.2) is recursive in nature. Therefore

After some algebraic calculations, the difference between the successive terms together with the model initial conditions, we obtain

where in (5.5) $S_{k}(t) = \sum_{i = 0}^{k}\Psi_{S, i}(t), V_{k}(t) = \sum_{i = 0}^{k}\Psi_{V, i}(t), E_{k}(t) = \sum_{i = 0}^{k}\Psi_{E, i}(t), I_{k}(t) = \sum_{i = 0}^{k}\Psi_{I, i}(t)$ $\text{and}\; P_{k}(t) = \sum_{i = 0}^{k}\Psi_{P, i}(t).$ Furthermore, let

then we have that

Thus, since the conditions for a Lipschitz function are all satisfied by $S(t), V(t), E(t), I(t), P(t)$ and the functions $\mathcal{K}_{i}$ for $i = 1, \cdots, 5,$ we conclude that it is bounded. The proof ends here.

5.3.   The disease free equilibrium and the basic reproduction number

At disease free equilibrium (DFE) we have a whole susceptible population without an infection, therefore $E^{*} = I^{*} = P^{*} = 0.$ By applying Definition 3, the DFE denoted by $\mathscr{E}_{0}$ is thus given by

We employ the method of Next-generation matrix to compute the basic reproduction number ($\mathscr{R}_{0}$) of COVID-19 infection spread as in ^[36]. The rate of generation of new infections and transition matrices results from the compartments $E, \; I$ and $P.$ Hence, evaluation at the DFE we have that

in which $\phi_{2} = \Big(\dfrac{\mu+\alpha\omega}{\mu+\alpha}\Big)\cdot$ The spectral radius is thus the maximum eigenvalue of the second generation matrix $\mathscr{R}_{0} = \tau (FV^{-1})$ and is given by

Thus,

where

In the above expression, $R_{0I}$ denotes the contribution from the infected population while $R_{0P}$ represents the contribution from the pathogens in the environment. However, we note that $\dfrac{1}{Q_{1}Q_{2}}$ and $\dfrac{1}{Q_{1}Q_{2}\theta}$ connotes the duration of stay in the $I(t)$ and $P(t)$ compartments respectively.

5.4.   Local stability of the disease free equilibrium

The Jacobian matrix evaluated at the DFE gives

The eigenvalues from matrix (5.6) are $-\mu, \; -Q_{0}$ and the solutions to the determinant from the matrix

Equation (5.7) yields the following characteristic equation

where

We note that $\mathscr{R}_{0I} < \mathscr{R}_{0}$ and hence $a_{0} > 0$ and $a_{1} > 0$ if $\mathscr{R}_{0} < 1.$ Also, it is clear that $a_{1}a_{2} > a_{0}$. Therefore, their eigenvalues have negative real parts using the principles of the Routh-Hurwitz criterion. The argument of the roots of Eq (5.8) are all greater than $\frac{\tau \pi}{2}$ if $\mathscr{R}_{0} < 1.$ Hence $\mathscr{E}_{0}$ is locally asymptotically stable whenever $\mathscr{R}_{0} < 1.$

5.5.   Endemic equilibrium

The endemic equilibrium (EE) is the condition that there is a COVID-19 patient and virus in the environment, as well as the spread of the disease occurs in population. The EE denoted by $\mathscr{E}_{1} = (S^{*}, V^{*}, E^{*}, I^{*}, P^{*})$ is obtained when $S^{*}, V^{*}, E^{*}, I^{*}, P^{*}$ is not equal to zero. Equation system (4.2) can be solved using the condition of the force of infection at equilibrium point ($\lambda^{*}$), with

Setting the right-hand sides of the model (4.2) to zero and noting $\lambda = \lambda^{*}$ at equilibrium gives

where $\phi_{0} = \dfrac{\gamma }{Q_{2}}\; \text{and}\; \phi_{1} = \dfrac{\xi}{\theta}\cdot$ Substituting (5.10) into (5.9) and after some algebraic simplification we have that the endemic equilibrium of the model satisfies $\lambda^{*} = 0$ or

where

We note that $\lambda^{*} = 0$ corresponds to the DFE. Also, the coefficient $\nu_{3}$ in Eq (5.11) is always positive while the value of the coefficient $\nu_{1}$ is positive or negative depending on the value of $\mathscr{R}_0.$ If $\mathscr{R}_0 < 1$ then $\nu_{1}$ is positive and if $\mathscr{R}_0 > 1$ then $\nu_{1}$ is negative. Endemic equilibrium will exist when the solution of the Eq (5.11) is positive $(\lambda > 0).$ Thus, we establish the following results on the existence of the endemic steady states of the model (4.2):

Theorem 1. Model (4.2) has

(i) Precisely one unique endemic equilibrium if $\nu_{1} < 0$ or $\mathscr{R}_{0} > 1;$

(ii) Precisely one unique endemic equilibrium if $\nu_{2} < 0$ and either $\nu_{1} = 0$ or $\nu_{2}^{2}-4\nu_{1}\nu_{3} = 0;$

(iii) Precisely two unique endemic equilibrium if $\nu_{2} < 0,$ $\nu_{1} > 0$ and $\nu_{2}^{2}-4\nu_{1}\nu_{3} > 0;$

(iv) No endemic equilibrium otherwise.

Using the Descartes's rule of signs change, we identify the possible number of EE exhibited by Eq (5.11) as summarised in Table 3 below.

We now present the bifurcation analysis in the next subsection.

5.6.   Bifurcation analysis

According to Kuznetsov ^[37], a bifurcation is a change that occurs if the topological structure of a system changes its parameters pass through a critical value. Hence, the system's behaviour changes as the basic reproduction number passes the value 1, that is at $\mathcal{R}_{0} = 1$ being the critical point. However, the occurrence of a backward bifurcation for any system has critical biological implications as it is not enough to $\mathcal{R}_{0}$ below unity to eliminate the disease. Thus, the reduction of $\mathcal{R}_{0}$ below the critical points denoted here as $\mathcal{R}_{0}^{c}$ is necessary to eliminate the disease.

We use the Center Manifold Theorem (CMT) as in ^[38], to show that model (4.2) exhibits a bifurcation. To employ the Center Manifold theory, we make a change of variables as follows $S = x_1$, $V = x_2$, $E = x_3$, $I = x_4,$ $P = x_5$ and using the vector notation $\textbf{X} = (x_1, x_2, x_3, x_4, x_5)^{T}$ for our model can be written in the form

given below

Suppose $\bar{\beta}$ to be the bifurcation parameter such that $\bar{\beta} = \beta_{1}$ and solving $\mathscr{R}_0 = 1,$ we have that

We compute the left and right eigenvectors of the Jacobian matrix at $\mathcal{R}_0 = 1$ and hence they are given by $W = (w_1, w_2, w_3, w_4, w_5)$ and $V = (v_1, v_2, v_3, v_4, v_5).$ The Jacobian for model system (3.1) evaluated at the DFE is found to be

The matrix $\mathscr{J}(\mathscr{E}_0\mid \mathscr{R}_0 = 1)$ has zero eigenvalue as one of its eigenvalues, thus CMT can be applied to establish the existence of bifurcation. The right and left eigenvectors are found to be

and

Using the result from Theorem 4.1 of Castillo-Chavez and Song ^[38], we calculate the coefficients $a$ as

and that of $b$ to be

where $\phi_3 = \dfrac{\mu}{\alpha+\mu}$ and $\phi_4 = \dfrac{\alpha\omega}{\alpha+\mu}.$ Hence, we have that $a$ can either be positive or negative and $b > 0.$ This implies that our model exhibits a forward bifurcation exists $a < 0$ and $b > 0$ as shown in Figure 3(a) and otherwise a backward bifurcation exists when $a > 0$ and $b > 0$ as illustrated in Figure 3(b) respectively. We observe changes in the qualitative behaviour of the model (3.1) whenever $\mathcal{R}_{0} = 1$, with $\mathcal{R}_{0}$ being the bifurcation parameter. For values of $\mathcal{R}_{0}$ greater than unity, we have a forward bifurcation, which implies that the disease will persist in the population, and decreasing $\mathcal{R}_{0}$ to values below one is not a sufficient condition to eradicate the disease. The existence of a backward bifurcation makes disease control difficult because of the co-existence of the DFE and the EE for values of $\mathcal{R}_{0}$ between $\mathcal{R}_{0}^{c}$ and 1. Therefore increase in the removal rate of the pathogen from the environment, $\theta,$ is not sufficient to introduce COVID-19 disease eradication. Implementing other control measures such as educational campaigns, effective contact tracing, quarantine and isolation of the infected, and proper hygiene is necessary to bring the population to a DFE state.

The dashed curve in Figure 4 is generated from Figure 3(b) by increasing the value of $\theta$ while decreasing the numerical value of $\beta_{1}.$ This implies that as more people practice proper hygiene, including hand washing, effective mask usage, and sanitising, it leads to human behaviour change, decreasing the effective transmission rate of COVID-19. These changes in the values of $\theta$ and $\beta_{1}$ leads to a change from backward bifurcation to a forward bifurcation, therefore making it easier to control and contain the disease. We note that most epidemiological models exhibit backward bifurcation and are influenced by a single model parameter. However, our model presents a unique and interesting scenario whereby the backward bifurcation is driven by two vital epidemiological model parameters, which are $\beta_{1}$ and $\theta.$ We now present the numerical simulations in the next section.

5.7.   Effects of parameters $\beta_{2}$ and $\alpha$ on $\mathscr{R}_{0}$

Figure 5 shows the transmission contact rate for the environment and the vaccination rate of the susceptible individual as a function of the basic reproduction number. We plotted these parameters against $\mathscr{R}_{0}$ to show their impact when $\mathscr{R}_{0} = 1$ and $\mathscr{R}_{0} = 2.$ This is because the $\mathscr{R}_{0}$ for COVID-9 is greater than 2 from model predictions for this pandemic in Indonesia (see ^[39] for instance). It can be observed that an increase in $\beta_{2}$ increases the $\mathscr{R}_{0}$ while an increase in $\alpha$ decreases the value of $\mathscr{R}_{0}.$ Further, at $\mathscr{R}_{0} = 1$ the disease will be eradicated, as also shown using mathematical analyses.

6.   Numerical simulations

In this section, we carry out numerical simulations of our model and in the following subsection we present the epidemiological parameter estimates obtained under real data from the COVID-19 active case for Indonesia. The data were used to calibrate our model and the results from fitting the data using the least square method presented hitherto.

Numerical results of the fractional model

Here, we present numerical simulations based on the fractional model using the Caputo operator. The fractional model (4.2) is then solved numerically by using the method as described in ^[40] and the biological parameter values is given in Table 2. The approximate value of the reproduction number using fitted parameters from Table 2 is $\mathscr{R}_0 \approx 8.4711$.

Figure 6 depicts the model versus data fitting for the case when $\tau = 1, 0.99$. It is clear from the real data case of COVID-19 that the fractional model has good resemblance fitting than the integer model.

Next, we varied some parameters and the Caputo operator ($\tau$), whose choice was considered according to this paper's objective to analyze theirs on the dynamical behaviour of COVID-19 incidence/resurgent in Indonesia. The simulation results of varying $\beta_{2}, \; \omega$, and $\tau$ on the state variables over time are shown graphically in Figures 7–9 respectively. In Figure 7, it can be seen that an increase in $\beta_2$ and $\tau$ leads to a decrease in the numbers of infected individuals who are infectious. While an increase in the rate of exposure of vaccinated individuals also results in a decrease in the number of infectives, as can be observed in Figure 8. It can also be observed that in both Figures 7 and 8, there is a comparative faster decrease in the number infected population at the value of $\tau = 0.95$ than when $\tau = 1$ as time progresses. Hence, the biological implication of increasing the transmission rate from the pathogen and amplification parameter, $\omega,$ will help reduce the COVID-19 transmission.

On the other hand, Figure 9 shows the potential impact of varying $\tau$ on the susceptible, exposed, infected, vaccinated, and the pathogen in the environment. Moreover, in Figure 9(a), (e) it is noticed that a decrease in the value of $\tau$ from 1 to 0.9 results to a decrease in the healthy population and the pathogen, even though the latter has a lower peak with an increase in time. In addition, from Figure 9(c), (d), we observe a trend where the increase $\tau$ results in an increase in the infected and the exposed population. However, that of the exposed has a higher peak over time when compared to that of the infectious population. However, the variation of $\tau$ on the vaccinate gives an interesting result. We see in 9b that varying $\tau$ over time leads to an increase and a decrease in the number of vaccinated individuals before and after day 200, respectively. The significance of this result can be attributed to human behaviour change during the pandemic.

7.   Conclusions

In this paper, we studied the effect of environmental transmission and vaccination on the dynamics of the COVID-19 pandemic in Indonesia. Mathematical analyses, which envisaged the integer and fractional model, were analysed for the stability and properties investigated. The disease-free equilibrium of the fractional model is locally asymptotically stable if the reproduction number is less than one. In addition, the model was found to exhibit bifurcation. We adopted the least square method to estimate and fit the model parameters using the reported active cases of COVID-19 in Indonesia. Further, the real data for the fractional model versus the integer-order model are compared. According to the simulation result, it can be seen that the fractional model has a close resemblance to the integer model. The impact of some of the key model parameters on the disease dynamics and its elimination is then shown graphically for various values of fractional-order of the Caputo derivative. The simulations show that an increase in the transmission rate of the virus from our environment increases the basic reproduction number, which contributes to more infection. Although the implementation of the vaccine plays an important role in controlling the spread of COVID-19, the remaining challenge is the efficacy of the vaccine used to deal with the new variant of COVID-19. Moreover, we conclude that the use of the fractional epidemic model provides a better understanding and biologically more insights into the disease dynamics.

Acknowledgements

This work has been supported by the Ministry of Research, Technology and Higher Education, Republic of Indonesia through the PDUPT 2021.

Conflict of interest

The authors declare that there is no conflict of interest.