Application of the Modified Adomian Decomposition Method on a Mathematical Model of COVID-19

Justina Mulenga^*, Patrick Azere Phiri

Department of Mathematics, School of Mathematics and Natural Sciences, The Copperbelt University, Kitwe, Zambia.

DOI: 10.4236/jamp.2023.119169                

Abstract

In this study, we constructed and analysed a mathematical model of COVID-19 in order to comprehend the transmission dynamics of the disease. The reproduction number (R_C) was calculated via the next generation matrix method. We also used the Lyaponuv method to show the global stability of both the disease free and endemic equilibrium points. The results showed that the disease-free equilibrium point is globally asymptotically stable if R_C < 1 and the endemic equilibrium point is globally asymptotically stable if R_C > 1. We further used the Adomian decomposition method and the modified Adomian decomposition method to obtain the solutions of the model. Numerical analysis of the model was done using Sagemath 9.0 software.

Keywords

COVID-19, Stability Analysis, Equilibrium Points, Adomian Decomposition Method, Modified Adomian Decomposition Method, Numerical Analysis

1. Introduction

Mathematical disease models are important tools in analyzing the spread and the control of infectious diseases. These models are based on dividing the host population into compartments, each containing individuals that are identical in terms of their status with respect to the disease under consideration. For example, in the SIR model there are three compartments namely Susceptibles (S), Infections (I) and Removed (R). From the compartments, the set of equations will arise. The equations specify how the sizes of the compartment change overtime, and are linear and nonlinear. Solutions of these equations will yield, S(t) the size of susceptibles at time t, I(t) the size of infectious individuals at the time t, and R(t) the size of the removed population at a given time t. Solving nonlinear equations is more challenging than solving linear equations. Different authors use different methods to solve nonlinear equations.

In [1] , Homotopy Analysis Method (HAM) was used to solve the differntial equations of the Susptible, Exposed, Infective and Removed (SEIR) epidemic model. The results showed that the series solution converged very fast and the solutions were accurate. The Adomian Decomposition Method (ADM) was used to solve the mathematical model of malaria in [2] . The Laplace Adomian Decomposition Method (LADM) was applied to analyse the SEIR epidemic models of measles in reference [3] . In [4] , the ADM was employed to solve the Human Immunodeficiency Virus (HIV) infection model of latently infected cells. In [5] , the Laplace-Adomian decomposition method and the Homotopy perturbation techniques were used to obtain approximate solutions of the equation of the fractional dynamics of a corona virus mathematical model under a Caputo derivative.

The novel corona virus, also known as COVID-19, is a deadly disease which came into being in late 2019. It is believed to break out in Wuhan city in China [6] . As of 23^rd May 2022, it had affected 227 countries with a total of 527,804,291 cases worldwide and 6,300,525 deaths [7] . Infected individuals present with fever, coughing, and sneezing as some of the symptoms. In severe cases, patients are diagnosed with pneumonia and shortness of breath. It is transmitted from human to human by inhaling droplets exhaled through normal breathing, coughing, and sneezing from infected persons [8] . Approximate incubation period is from 2 - 14 days from the time of contact. However, it may go up to 27 days. As for the asymptomatic individuals, they do not develop any symptoms and are not aware of their infection, yet they can transmit the disease [9] . As a result, preventive measures such as lock-down policy, use of face masks in public places, washing of hands with soap, and using hand sanitizer were introduced in most affected countries.

This study focuses on the construction of a mathematical model to analyze the transmission dynamics of COVID-19 and solve it with modified Adomian decomposition method. The rest of the paper is organized as follows: in Section 2 we formulate the mathematical model of COVID-19 and analyze it. In Section 3 we solve the mathematical model of COVID-19 using Adomain Decomposition Method (ADM) and Modified Adomian Decomposition Method (MADM). Section 4 presents the numerical analysis of the model and discussion. The conclusion comes in Section 5.

2. Model Formulation

The mathematical model of COVID-19 is developed by dividing the total population N into seven classes. The variables Susceptible (S), Exposed (E), Symptomatic (I_s), Asymptomatic (I_a), Hospitalized (H), Quarantined (Q) and Recovered (R) are used to represent the classes. Individuals are recruited into the susceptible class through the rate Θ. The susceptible population is exposed to the disease through contact with symptomatic and asymptomatic infectious individuals. The parameters ${\rho }_s$ and ${\rho }_a$ represent the effective contact rates for individuals in the symptomatic and asymptomatic infectious classes, respectively. There is a fraction of individuals who use face masks in the population and it is given as $0<\kappa \le 1$ , whereas $0<\omega \le 1$ represents the expected decrease in the risk of infection as a result of using face masks. The exposed individuals progress to infectious classes at the rate $\sigma$ . A fraction $0<q\le 1$ of the exposed shows no symptoms and they proceed to asymptomatic infectious class whereas the remaining fraction $\left(1-q\right)$ shows symptoms of the disease and hence proceeds to the symptomatic class. The symptomatic infectious are hospitalized at the rate of $\pi$ and are quarantined at the rate of $\gamma$ . The asymptomatic infectious recover at the rate $\tau$ . The quarantined are hospitalized at the rate $\epsilon$ and they recover at the rate $\nu$ whereas the hospitalized recover at the rate $\beta$ . There is a natural death rate of $\delta$ for individuals in all classes. Additionally, individuals in the symptomatic, hospitalized, asymptomatic, and quarantined classes have disease-induced death rates of ${\alpha }_1$ , ${\alpha }_2$ , ${\alpha }_3$ , and ${\alpha }_4$ , respectively. The summary of the description of the variables and parameters is given in Table 1.

Using the symbols and variables described in Table 1 we draw the compartmental model that shows the progression of the disease, given in Figure 1.

Using the variables, parameters and the compartmental model we derive the system of ordinary differential model equations as follows:

$\begin{array}{l}\frac{\text{d}S}{\text{d}t}=\Theta -{\rho }_s\left(1-\kappa \omega \right)S{I}_s-{\rho }_a\left(1-\kappa \omega \right)S{I}_a-\delta S,\\ \frac{\text{d}E}{\text{d}t}={\rho }_s\left(1-\kappa \omega \right)S{I}_s+{\rho }_a\left(1-\kappa \omega \right)S{I}_a-\left(\delta +\sigma \right)E,\\ \frac{\text{d}{I}_s}{\text{d}t}=\left(1-q\right)\sigma E-\left(\delta +{\alpha }_1+\gamma +\pi \right)I_s,\\ \frac{\text{d}{I}_a}{\text{d}t}=q\sigma E-\left(\delta +{\alpha }_3+\tau \right)I_a,\\ \frac{\text{d}H}{\text{d}t}=\pi I_s+\epsilon Q-\left(\delta +\beta +{\alpha }_2\right)H,\\ \frac{\text{d}Q}{\text{d}t}=\gamma I_s-\left(\delta +\epsilon +\nu +{\alpha }_4\right)Q,\\ \frac{\text{d}R}{\text{d}t}=\beta H+\nu Q+\tau I_a-\delta R.\end{array}\}$ (1)

The total population is denoted by N and given as,

$N=S+E+I_a+I_s+H+Q+R$ (2)

and $\frac{\text{d}N}{\text{d}t}$ is given by

$\frac{\text{d}N}{\text{d}t}=\Theta -\delta N-\left({\alpha }_1{I}_s+{\alpha }_{2}H+{\alpha }_3{I}_a+{\alpha }_{4}Q\right).$ (3)

Figure 1. Compartmental model of the transmission of COVID-19.

Table 1. Symbols and description of parameters.

2.1. Positivity and Boundedness of Solutions

2.1.1. Positivity of Solutions

Since the variables of model equation system (1) represent humans, it is important to show that they are positive.

Theorem 1 If the initial values are given by $S\left(0\right)>0$ , $E\left(0\right)>0$ , $I_a\left(0\right)>0$ , $I_s\left(0\right)>0$ , $H\left(0\right)>0$ , $Q\left(0\right)>0$ , $R\left(0\right)>0$ , then the solutions $S\left(t\right)>0$ , $E\left(t\right)>0$ , $I_a\left(t\right)>0$ , $I_s\left(t\right)>0$ , $H\left(t\right)>0$ , $Q\left(t\right)>0$ , $R\left(t\right)>0$ of the system Equation (1) are positive for all $t>0$ .

Proof: Using the first equation of equation system (1), we have the following:

$\begin{array}{c}\frac{\text{d}S}{\text{d}t}=\Theta -\left[{\rho }_s\left(1-\kappa \omega \right)I_s+{\rho }_a\left(1-\kappa \omega \right)I_a+\delta \right]S\\ \ge -\left[{\rho }_s\left(1-\kappa \omega \right)I_s+{\rho }_a\left(1-\kappa \omega \right)I_a+\delta \right]S,\end{array}$

$\frac{\text{d}S}{\text{d}t}\ge -\left[{\rho }_s\left(1-\kappa \omega \right)I_s+{\rho }_a\left(1-\kappa \omega \right)I_a+\delta \right]S,$

${\displaystyle {\int }_0^t}\frac{\text{d}S}{S}\text{d}S\ge -{\displaystyle \int }\left[{\rho }_s\left(1-\kappa \omega \right)I_s+{\rho }_a\left(1-\kappa \omega \right)I_a+\delta \right]\text{d}t,$

$\ln S\left(t\right)-\ln S\left(0\right)\ge -{\displaystyle \int }\left[{\rho }_s\left(1-\kappa \omega \right)I_s+{\rho }_a\left(1-\kappa \omega \right)I_a+\delta \right]\text{d}t,$

$S\left(t\right)\ge S\left(0\right){\text{e}}^{-{\displaystyle \int }\left[{\rho }_s\left(1-\kappa \omega \right)I_s+{\rho }_a\left(1-\kappa \omega \right)I_a+\delta \right]\text{d}t}.$

Thus S is positive since $S\left(0\right)$ is positive and the exponential function ${\text{e}}^{-{\displaystyle \int }\left[{\rho }_s\left(1-\kappa \omega \right)I_s+{\rho }_a\left(1-\kappa \omega \right)I_a+\delta \right]\text{d}t}$ is always positive. Using the same method, we can prove the rest of the equations of system Equation (1) and show that $E\left(t\right)>0$ , $I_a\left(t\right)>0$ , $I_s\left(t\right)>0$ , $H\left(t\right)>0$ , $Q\left(t\right)>0$ , $R\left(t\right)>0$ .

2.1.2. Boundedness of the Solutions

Theorem 2 All positive solutions presented in Theorem (1) are bounded.

Proof: In the absence of the disease, i.e. ${\alpha }_1={\alpha }_2={\alpha }_3={\alpha }_4=0$ , Equation (3) becomes $\frac{\text{d}N}{\text{d}t}=\Theta -\delta N$ , which can be written as $N=\frac{\Theta }{\delta }-\left[\frac{\Theta -\Theta N_0}{\Theta }\right]{\text{e}}^{-\delta t}$ , and as $t\to \infty$ , $N\to \frac{\Theta }{\delta }$ . Thus, we conclude that

$N=\frac{\Theta }{\delta }\mathrm{.}$ (4)

Hence, the positive solutions of model Equation (1) are bounded.

2.2. Existence and Stability of Equilibrium Points

The equilibrium points are calculated by equating the left hand side of equation system (1) to zero. This leads to the following system of equations,

$\begin{array}{l}\Theta -{\rho }_s\left(1-\kappa \omega \right)S{I}_s-{\rho }_a\left(1-\kappa \omega \right)S{I}_a-\delta S=0,\\ {\rho }_s\left(1-\kappa \omega \right)S{I}_s+{\rho }_a\left(1-\kappa \omega \right)S{I}_a-(\delta +\sigma )E=0,\\ \left(1-q\right)\sigma E-\left(\delta +{\alpha }_1+\xi +\pi +\gamma \right)I_s=0,\\ q\sigma E-\left(\delta +{\alpha }_3+\tau \right)I_a=0,\\ \pi I_s+\epsilon Q-\left(\beta +\delta +{\alpha }_2\right)H=0,\\ \gamma I_s-\left(\epsilon +\nu +\delta +{\alpha }_4\right)Q=0,\\ \xi I_s+\nu Q+\tau I_a-\delta R=0.\end{array}\}$ (5)

2.2.1. The Disease Free Equilibrium Point ( $d_{fep}$ )

In the absence of infection $I_a=I_s=H=Q=0$ . Solving the equations of system 5, we obtain $d_{fep}=\left(\frac{\Theta }{\delta }\mathrm{,0,0,0,0,0,0}\right)$ .

2.2.2. The Reproduction Number of the Model

We compute the basic reproduction number $R_C$ of the model system (1) by following [10] . The basic reproduction number measures the average number of new infections generated by a single infected person in a completely susceptible population. Using the next generation matrix method, we consider the classes that are actively transmitting the disease. These are:

$\begin{array}{l}\frac{\text{d}E}{\text{d}t}={\rho }_s\left(1-\kappa \omega \right)S{I}_s+{\rho }_a\left(1-\kappa \omega \right)S{I}_a-\left(\delta +\sigma \right)E,\\ \frac{\text{d}{I}_s}{\text{d}t}=\left(1-q\right)\sigma E-\left(\delta +{\alpha }_1+\xi +\pi +\gamma \right)I_s,\\ \frac{\text{d}{I}_a}{\text{d}t}=q\sigma E-\left(\delta +{\alpha }_3+\tau \right)I_a.\end{array}\}$ (6)

Next, we find the matrices $F_i$ and $V_i$ which are the rate of appearance of new infections in compartment i and transfer of individuals into and out of compartment i by all other means, respectively, given by

$\begin{array}{l}{F}_i=\left(\begin{array}{c}\left(1-\kappa \omega \right)I_{s}S{\rho }_s+\left(1-\kappa \omega \right)I_{a}S{\rho }_a\\ 0\\ 0\end{array}\right);\\ V_i=\left(\begin{array}{c}\left(\delta +\sigma \right)E\\ \left({\alpha }_1+\delta +\gamma +\pi +\xi \right)I_s-\left(1-q\right)\sigma E\\ \left({\alpha }_3+\delta +\tau \right)I_a-q\sigma E\end{array}\right).\end{array}$

Then taking partial derivatives of both $F_i$ and $V_i$ on the disease-free equilibrium point we get,

$F=\left(\begin{array}{ccc}0& \left(1-\kappa \omega \right)S{\rho }_s& \left(1-\kappa \omega \right)S{\rho }_a\\ 0& 0& 0\\ 0& 0& 0\end{array}\right)$

and

$V=\left(\begin{array}{ccc}\delta +\sigma & 0& 0\\ -\left(1-q\right)\sigma & {\alpha }_1+\delta +\gamma +\pi +\xi & 0\\ -q\sigma & 0& {\alpha }_3+\delta +\tau \end{array}\right).$

Next we calculate the inverse of V which is given by:

$V^{-1}=\left(\begin{array}{ccc}\frac{1}{\delta +\sigma }& 0& 0\\ \frac{\left(1-q\right)\sigma }{\left({\alpha }_1+\delta +\xi +\pi +\gamma \right)\left(\delta +\sigma \right)}& \frac{1}{{\alpha }_1+\delta +\xi +\pi +\gamma }& 0\\ \frac{q\sigma }{\left(\delta +\sigma \right)\left({\alpha }_3+\delta +\tau \right)}& 0& \frac{1}{{\alpha }_3+\delta +\gamma }\end{array}\right).$

The next generation matrix is given as the following product:

$F{V}^{-1}=\left(\begin{array}{ccc}\frac{{\rho }_s\left(1-\kappa \omega \right)S\left(1-q\right)\sigma }{\left(\delta +\sigma \right)\left({\alpha }_1+\delta +\xi +\gamma +\pi \right)}+\frac{{\rho }_a\left(1-\kappa \omega \right)Sq\sigma }{\left(\delta +\sigma \right)\left({\alpha }_3+\delta +\pi \right)}& \frac{{\rho }_s\left(1-\kappa \omega \right)S}{{\alpha }_1+\delta +\xi +\gamma +\pi }& \frac{{\rho }_a\left(1-\kappa \omega \right)}{{\alpha }_3+\delta +\pi }\\ 0& 0& 0\\ 0& 0& 0\end{array}\right)\mathrm{,}$

and we calculate the eigenvalues of $F{V}^{-1}$ as follows:

$\left|\begin{array}{ccc}\left(\frac{{\rho }_s\left(1-\kappa \omega \right)S\left(1-q\right)\sigma }{\left(\delta +\sigma \right)\left({\alpha }_1+\delta +\xi +\gamma +\pi \right)}+\frac{{\rho }_a\left(1-\kappa \omega \right)Sq\sigma }{\left(\delta +\sigma \right)\left({\alpha }_3+\delta +\pi \right)}\right)-{\lambda }_1& \frac{{\rho }_s\left(1-\kappa \omega \right)S}{{\alpha }_1+\delta +\xi +\gamma +\pi }& \frac{{\rho }_a\left(1-\kappa \omega \right)}{{\alpha }_3+\delta +\pi }\\ 0& 0-{\lambda }_2& 0\\ 0& 0& 0-{\lambda }_3\end{array}\right|\mathrm{<mo>\; =\; </mo>0.}$

Thus ${\lambda }_1=\frac{{\rho }_s\left(1-\kappa \omega \right)S\left(1-q\right)\sigma }{\left(\delta +\sigma \right)\left({\alpha }_1+\delta +\xi +\gamma +\pi \right)}+\frac{{\rho }_a\left(1-\kappa \omega \right)Sq\sigma }{\left(\delta +\sigma \right)\left({\alpha }_3+\delta +\pi \right)}$ , ${\lambda }_2=0$ , ${\lambda }_3=0$ .

The reproduction number is given by the largest eigenvalue of the determinant of the matrix $F{V}^{-1}$ and so

$R_C=\frac{{\rho }_s\left(1-\kappa \omega \right)S\left(1-q\right)\sigma }{\left(\delta +\sigma \right)\left({\alpha }_1+\delta +\gamma +\pi +\xi \right)}+\frac{{\rho }_a\left(1-\kappa \omega \right)Sq\sigma }{\left(\delta +\sigma \right)\left({\alpha }_3+\delta +\tau \right)}.$ (7)

Equation (7) consists of two reproduction numbers. The first reproduction number $\frac{{\rho }_s\left(1-\kappa \omega \right)S\left(1-q\right)\sigma }{\left(\delta +\sigma \right)\left({\alpha }_1+\delta +\gamma +\pi +\xi \right)}$ defines the number of new COVID-19 cases generated from the symptomatic infected individuals in class $I_s$ . The second reproduction number is $\frac{{\rho }_a\left(1-\kappa \omega \right)Sq\sigma }{\left(\delta +\sigma \right)\left({\alpha }_3+\delta +\tau \right)}$ which defines the number

of new COVID-19 cases generated from the asymptomatic infectious individuals in class $I_a$ . Hence the reproduction number is written as,

$R_C=R_{I_s}+R_{I_a}.$ (8)

2.2.3. Global Stability of the Disease-Free Equilibrium Point

The global stability of the disease-free equilibrium $D_{fep}$ is established by the following theorem.

Theorem 3 If $R_C<1$ , the disease-free equilibrium is globally asymptotically stable in $\Omega$ and unstable if $R_C>1$

Proof: Let the Lyapunov function be

$V=d_{1}E+d_2{I}_s+d_3{I}_a,$ (9)

where,

$d_1=\left(\delta +{\alpha }_3+\tau \right)\left(\delta +{\alpha }_1+\gamma +\pi \right),$

$d_2=\left(\delta +{\alpha }_3+\tau \right){\rho }_s\left(1-\kappa \omega \right)S,$

$d_3=\left(\delta +{\alpha }_1+\gamma +\pi \right){\rho }_a\left(1-\kappa \omega \right)S.$

From Equation (9), the derivative is given as

$\frac{\partial V}{\partial t}=d_1\frac{\text{d}E}{\text{d}t}+d_2\frac{\text{d}{I}_s}{\text{d}t}+d_3\frac{\text{d}{I}_a}{\text{d}t}.$ (10)

But $\frac{\text{d}E}{\text{d}t},\frac{\text{d}{I}_s}{\text{d}t},\frac{\text{d}{I}_a}{\text{d}t}$ are given in equation system (1). Therefore,

$\begin{array}{c}\frac{\partial V}{\partial t}=\left(\delta +{\alpha }_3+\tau \right)\left(\delta +{\alpha }_1+\gamma +\pi \right)\left[{\rho }_s\left(1-\kappa \omega \right)S{I}_s+{\rho }_a\left(1-\kappa \omega \right)S{I}_a-\left(\delta +\sigma \right)E\right]\\ +\left(\delta +{\alpha }_3+\tau \right){\rho }_s\left(1-\kappa \omega \right)S\left[\left(1-q\right)\sigma E-\left(\delta +{\alpha }_1+\pi +\gamma \right)I_s\right]\\ +\left(\delta +{\alpha }_1+\gamma +\pi \right){\rho }_a\left(1-\kappa \omega \right)S\left[q\sigma E-\left(\delta +{\alpha }_3+\tau \right)I_a\right].\end{array}$ (11)

The expansion of Equation (11) yields,

$\begin{array}{c}\frac{\partial V}{\partial t}=\left(\delta +{\alpha }_3+\tau \right)\left(\delta +{\alpha }_1+\gamma +\pi \right){\rho }_s\left(1-\kappa \omega \right)S{I}_s\\ +\left(\delta +{\alpha }_3+\tau \right)\left(\delta +{\alpha }_1+\gamma +\pi \right){\rho }_a\left(1-\kappa \omega \right)S{I}_a\\ -\left(\delta +{\alpha }_3+\tau \right)\left(\delta +{\alpha }_1+\gamma +\pi \right)\left(\delta +\sigma \right)E\\ +\left(\delta +{\alpha }_3+\tau \right){\rho }_s\left(1-\kappa \omega \right)S\left(1-q\right)\sigma E\end{array}$

$\begin{array}{c}-\left(\delta +{\alpha }_3+\tau \right){\rho }_s\left(1-\kappa \omega \right)S\left(\delta +{\alpha }_1+\gamma +\pi \right)I_s\\ +\left(\delta +{\alpha }_1+\gamma +\pi \right){\rho }_a\left(1-\kappa \omega \right)Sq\sigma E\\ -\left(\delta +{\alpha }_1+\gamma +\pi \right){\rho }_a\left(1-\kappa \omega \right)S\left(\delta +{\alpha }_3+\tau \right)I_a.\end{array}$ (12)

Equation (12) can be simplified to

$\begin{array}{c}\frac{\partial V}{\partial t}\le \left(\delta +{\alpha }_3+\tau \right){\rho }_s\left(1-\kappa \omega \right)S\left(1-q\right)\sigma E+\left(\delta +{\alpha }_1+\gamma +\pi \right){\rho }_a\left(1-\kappa \omega \right)Sq\sigma E\\ -\left(\delta +{\alpha }_3+\tau \right)\left(\delta +{\alpha }_1+\gamma +\pi \right)\left(\delta +\sigma \right)E,\\ \le \left(\delta +{\alpha }_3+\tau \right)\left(\delta +{\alpha }_1+\gamma +\pi \right)\left(\delta +\sigma \right)\\ \cdot \left[\frac{{\rho }_s\left(1-\kappa \omega \right)S\left(1-q\right)\sigma }{\left(\delta +{\alpha }_1+\gamma +\pi \right)\left(\delta +\sigma \right)}+\frac{{\rho }_a\left(1-\kappa \omega \right)Sq\sigma }{\left(\delta +{\alpha }_3+\tau \right)\left(\delta +\sigma \right)}-1\right]E,\\ \le \left(\delta +{\alpha }_3+\tau \right)\left(\delta +{\alpha }_1+\gamma +\pi \right)\left(\delta +\sigma \right)\cdot \left[R_C-1\right]E.\end{array}$

Hence, $\frac{\partial V}{\partial t}\le 0$ if $R_C<1$ and $\frac{\partial V}{\partial t}=0$ if $E=0$ . By LaSalle’s Invariance

Principle [11] , we conclude that the disease-free equilibrium ( $D_{fep}$ ) of the model of COVID-19 is globally asymptotically stable in $\Omega$ whenever $R_C<1$ .

2.2.4. The Endemic Equilibrium Point

The endemic equilibrium point is denoted by $E_{ep}=\left\{S^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{E}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{I}_s^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{I}_a^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{H}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{Q}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{R}^{\mathrm{<mo>\; *\; </mo>}}\right\}$ and is calculated as,

$\begin{array}{l}{S}^{*}=\frac{\Theta }{{\rho }_s\left(1-\kappa \omega \right)I_s+{\rho }_a\left(1-\kappa \omega \right)I_a+\delta },\\ E^{*}=\frac{{\rho }_s\left(1-\kappa \omega \right)I_s+{\rho }_a\left(1-\kappa \omega \right)I_a}{\delta +\sigma },\\ I_s^{*}=\frac{\left(1-q\right)\sigma E}{{\alpha }_1+\delta +\pi +\gamma },\\ I_a^{*}=\frac{q\sigma E}{{\alpha }_3+\tau +\delta },\\ H^{*}=\frac{\pi I_s+\epsilon Q}{{\alpha }_2+\delta +\beta },\\ Q^{*}=\frac{\gamma I_s}{{\alpha }_4+\delta +\nu +\epsilon },\\ R^{*}=\frac{\beta H+\nu Q+\tau I_a}{\delta }.\end{array}\}$ (13)

2.2.5. Global Stability Analysis of the Endemic Equilibrium Point

The Lyapunov asymptotic stability theorem is used to prove the global asymptotic stability of the endemic equilibrium point. Using the method for constructing the Lyapunov function discussed in [12] , we formulate the Lyapunov function for model Equation (1).

Theorem 4 If $R_C>1$ , then the endemic equilibrium point $E_{ep}$ of model Equation (1) is globally asymptotically stable in the region Ω.

Proof First, we define $L\mathrm{<mo>\; :\; </mo>}\left\{\left(S,E,I_s,I_a,Q,H,R\right)\in \Omega :S,E,I_s,I_a,Q,H,R>0\right\}\to \Re$ . Consider the function below:

$\begin{array}{c}L\left(S,E,I_s,I_a,Q,H,R\right)=\ln (\left(S-S^{*}\right)+\left(E-E^{*}\right)+\left(I_s-I_s^{*}\right)+\left(I_a-I_a^{*}\right)\\ +\left(Q-Q^{\mathrm{<mo>\; *\; </mo>}}\right)+\left(H-H^{\mathrm{<mo>\; *\; </mo>}}\right)+\left(R-R^{\mathrm{<mo>\; *\; </mo>}}\right)+1)\end{array}$ (14)

The derivative of L along the solutions of the model in Equation (1) is given by the expression:

$\begin{array}{c}\stackrel{\dot{}}{L}=\frac{\partial L}{\partial S}\frac{\text{d}S}{\text{d}t}+\frac{\partial L}{\partial E}\frac{\text{d}E}{\text{d}t}+\frac{\partial L}{\partial I_s}\frac{\text{d}{I}_s}{\text{d}t}+\frac{\partial L}{\partial I_a}\frac{\text{d}{I}_a}{\text{d}t}+\frac{\partial L}{\partial Q}\frac{\text{d}Q}{\text{d}t}+\frac{\partial L}{\partial H}\frac{\text{d}H}{\text{d}t}+\frac{\partial L}{\partial R}\frac{\text{d}R}{\text{d}t}\\ =\frac{\frac{\text{d}S}{\text{d}t}+\frac{\text{d}E}{\text{d}t}+\frac{\text{d}{I}_s}{\text{d}t}+\frac{\text{d}{I}_a}{\text{d}t}+\frac{\text{d}Q}{\text{d}t}+\frac{\text{d}H}{\text{d}t}+\frac{\text{d}R}{\text{d}t}}{\left(S-S^{\mathrm{<mo>\; *\; </mo>}}\right)+\left(E-E^{\mathrm{<mo>\; *\; </mo>}}\right)+\left(I_s-I_s^{\mathrm{<mo>\; *\; </mo>}}\right)+\left(I_a-I_a^{\mathrm{<mo>\; *\; </mo>}}\right)+\left(Q-Q^{\mathrm{<mo>\; *\; </mo>}}\right)+\left(H-H^{\mathrm{<mo>\; *\; </mo>}}\right)+\left(R-R^{\mathrm{<mo>\; *\; </mo>}}\right)+1}\end{array}$ (15)

From Equation (4), all solutions of Equation (13) satisfy the equality

$N^{*}=S^{*}+E^{*}+I_s^{*}+I_a^{*}+Q^{*},H^{*}+R^{*}=\frac{\Theta }{\delta }.$ (16)

In addition $N={\text{e}}^{\delta t+k}+\frac{\Theta }{\delta }$ where, k is the value that satisfies the condition $N_0=\frac{\Theta }{\delta }$ . Thus $L=\ln \left(N-N^{*}+1\right)\ge 0$ . Therefore,

$\stackrel{\dot{}}{L}=\frac{1}{N-N^{*}}\frac{\text{d}N}{\text{d}t}=\frac{\delta }{N-\frac{\Theta }{\delta }+1}\left(\frac{\Theta }{\delta }-N\right)\le 0.$

$L=0$ and $\stackrel{\dot{}}{L}=0$ are satisfied if and only if $N=\frac{\Theta }{\delta }$ .

L is positive definite and $\stackrel{\dot{}}{L}$ is negative definite, hence the function L is the Lyapunov function for model Equation (1) and the endemic equilibrium $E_{ep}$ is globally asymptotically stable by the Lyapunov asymptotic stability analysis [13] . Hence the proof.

3. Adomian Decomposition Method

In this section we apply the Adomian decomposition method proposed by George Adomian in the mid 1980’s, [14] to solve the system Equation (1).

3.1. Basic Concepts of the Adomian Decomposition Method

Let us consider the initial value problem expressed as

$Lu+Nu+Ru=g,$ (17)

where L is the linear operator, N is the nonlinear operator and R is the remaining linear part. By defining the inverse operator of L as $L^{-1}$ , we introduce it on both sides of Equation (17) to get,

$L^{-1}Lu=L^{-1}g-L^{-1}\left[Ru+Nu\right]\mathrm{.}$ (18)

Solving for u in Equation (18) leads to,

$u=\phi \left(0\right)+L^{-1}g-L^{-1}\left[Ru+Nu\right],$ (19)

where,

$\phi \left(0\right)=(\begin{array}{ll}u\left(0\right)\hfill & \text{if}L=\frac{\text{d}}{\text{d}x},\hfill \\ u\left(0\right)+x{u}^{\prime }\left(0\right)\hfill & \text{if}L=\frac{{\text{d}}^2}{\text{d}{x}^2},\hfill \\ u\left(0\right)+x{u}^{\prime }\left(0\right)+\frac{x^2}{2!}{u^{\prime }}^{\prime }\left(0\right)\hfill & \text{if}L=\frac{{\text{d}}^3}{\text{d}{x}^3}.\hfill \\ ⋮\hfill & \hfill \\ u\left(0\right)+x{u}^{\prime }\left(0\right)+\frac{x^2}{2!}{u^{\prime }}^{\prime }\left(0\right)+\cdots +\frac{x^n}{n!}{u}^n\left(0\right)\hfill & \text{if}L=\frac{{\text{d}}^{n+1}}{\text{d}{x}^{n+1}}.\hfill \end{array}$

The Adomian Decomposition Method assumes that the unknown function $u$ can be expressed by an infinite series of the form,

$u={\displaystyle \underset{n=0}{\overset{\infty }{\sum }}}{u}_n.$ (20)

where the component $u_n$ will be determined recursively. This method also defines the nonlinear term by the Adomian polynomials. More precisely, the ADM assumes that the nonlinear operator can be decomposed by an infinite series of polynomials given by,

$N\left(u\right)={\displaystyle \underset{n=0}{\overset{\infty }{\sum }}}{A}_n,$ (21)

where $A_n$ ’s are the Adomian’s polynomials defined as,

$A_n=A_n\left(u_0,u_1,u_2,\cdots ,u_n\right)$ (22)

and are calculated using the formular in [15]

$A_n=\frac{1}{n!}\frac{{\text{d}}^n}{\text{d}{\zeta }^n}{\left[F\left({\displaystyle \underset{k=0}{\overset{\infty }{\sum }}}{u}_k{\zeta }^k\right)\right]}_{\zeta =0},n=0,1,2,\cdots .$ (23)

Substituting Equation (20) and Equation (21) into Equation (19) and using the fact that R is a linear operator we obtain,

${\displaystyle \underset{n=0}{\overset{\infty }{\sum }}}{u}_n=\phi \left(0\right)+L^{-1}g-L^{-1}\left({\displaystyle \underset{n=0}{\overset{\infty }{\sum }}}R\left(u_n\right)\right)-L^{-1}\left({\displaystyle \underset{n=0}{\overset{\infty }{\sum }}}{A}_n\right)$ (24)

Therefore the formal recurrence algorithm could be defined by,

$\begin{array}{l}{u}_0=\phi \left(0\right)+L^{-1}g\\ ⋮\\ u_{n+1}=-L^{-1}\left(R\left(u_n\right)\right)-L^{-1}\left(A_n\right)\end{array}$ (25)

The n-term approximation of the solution is given by,

${\Psi }_n\left(x\right)={\displaystyle \underset{k=0}{\overset{n-1}{\sum }}}{u}_k\left(x\right).$ (26)

The advantage of this method is that it solves problems in a direct way and in uncomplicated manner without linearization, perturbation or any unpreferable assumptions that may change the physical behaviour of the problem under discussion.

Solutions of the Mathematical Model of COVID-19 Using Adomian Decomposition Method

In order to explicitly construct approximate solutions of the system described by Equation (1), we employ the Adomain decomposition method. To apply the Adomian decomposition method we choose: $S\left(0\right)$ , $E\left(0\right)$ , $I_s\left(0\right)$ , $I_a\left(0\right)$ , $H\left(0\right)$ , $Q\left(0\right)$ , $R\left(0\right)$ as initial approximations of $S\left(t\right)$ , $E\left(t\right)$ , $I_s\left(t\right)$ , $I_a\left(t\right)$ , $H\left(t\right)$ , $Q\left(t\right)$ , $R\left(t\right)$ . Using Equation (12) for calculating the Adomian polynomials for $S{I}_s$ and $S{I}_a$ as $A_n$ and $B_n$ respectively, we apply Equation (24) to each of the equations in model system (1) to obtain the recursive algorithm for each equation as follows:

$S_0=S\left(0\right)+L^{-1}\Theta$ ; $S_{n+1}=-{\rho }_s\left(1-\kappa \omega \right)L^{-1}\left(A_n\right)-{\rho }_a\left(1-\kappa \omega \right)L^{-1}\left(B_n\right)-\delta L^{-1}\left(S_n\right)$ , $n\ge 0$

$E_0=E\left(0\right)$ ; $E_{n+1}={\rho }_s\left(1-\kappa \omega \right)L^{-1}\left(A_n\right)+{\rho }_a\left(1-\kappa \omega \right)L^{-1}\left(B_n\right)-\left(\delta +\sigma \right)L^{-1}\left(E_n\right)$ , $n\ge 0$

$I_{s0}=I_s\left(0\right)$ ; $I_{s\left(n+1\right)}=\left(1-q\right)\sigma L^{-1}\left(E_n\right)-\left(\gamma +{\alpha }_1+\delta +\pi \right)L^{-1}\left(I_{sn}\right)$ , $n\ge 0$

$I_{a0}=I_a\left(0\right)$ ; $I_{a\left(n+1\right)}=q\sigma L^{-1}\left(E_n\right)-\left(\tau +{\alpha }_3+\delta \right)L^{-1}\left(I_{an}\right)$ , $n\ge 0$

$H_0=H\left(0\right)$ ; $H_{n+1}=\pi L^{-1}\left(I_{sn}\right)+\epsilon L^{-1}\left(Q_n\right)-\left(\beta +{\alpha }_2+\delta \right)L^{-1}\left(H_n\right)$ , $n\ge 0$

$Q_0=Q\left(0\right)$ ; $Q_{n+1}=\gamma L^{-1}\left(I_{sn}\right)-\left(\epsilon +\nu +{\alpha }_4+\delta \right)L^{-1}\left(Q_n\right)$ , $n\ge 0$

$R_0=R\left(0\right)$ ; $R_{n+1}=\beta L^{-1}\left(H_n\right)+\nu L^{-1}\left(Q_n\right)+\tau L^{-1}\left(I_{an}\right)-\delta L^{-1}\left(R_n\right)$ , $n\ge 0$ .

Using Equation (26) the solutons to Equation (1) are obtained.

3.2. The Modified Adomian Decomposition Method

Here, we present the Modified Adomian Decomposition Method (MADM) which we use to solve the system of equations arising from the mathematical model of COVID-19. In [16] , the authors proposed the modification of the ADM by introducing the terms ${\displaystyle {\sum }_{n=0}^{\infty }}{a}_n{t}^n-p{\displaystyle {\sum }_{n=0}^{\infty }}{a}_n{t}^n$ into the calculations of the standard ADM. To understand the procedure of MADM we consider Equation (23) and insert ${\displaystyle {\sum }_{n=0}^{\infty }}{a}_n{t}^n-p{\displaystyle {\sum }_{n=0}^{\infty }}{a}_n{t}^n$ to obtain:

${\displaystyle \underset{n=0}{\overset{\infty }{\sum }}}{u}_n=\phi \left(0\right)+L^{-1}g+L^{-1}{\displaystyle \underset{n=0}{\overset{\infty }{\sum }}}{a}_n{t}^n-p{L}^{-1}{\displaystyle \underset{n=0}{\overset{\infty }{\sum }}}{a}_n{t}^n-L^{-1}\left({\displaystyle \underset{n=0}{\overset{\infty }{\sum }}}R\left(u_n\right)\right)-L^{-1}\left({\displaystyle \underset{n=0}{\overset{\infty }{\sum }}}{A}_n\right)$ (27)

So the recursive algorithm for MADM is defined as,

$\begin{array}{l}{u}_0=\phi \left(0\right)+L^{-1}{\displaystyle \underset{n=0}{\overset{\infty }{\sum }}}{a}_n{t}^n\\ u_1=-p{L}^{-1}{\displaystyle \underset{n=0}{\overset{\infty }{\sum }}}{a}_n{t}^n-L^{-1}R\left(u_0\right)-L^{-1}\left(A_0\right)+L^{-1}g\\ ⋮\\ u_{n+1}=-L^{-1}\left(R\left(u_n\right)\right)-L^{-1}\left(A_n\right)\end{array}$ (28)

In this method we set $u_1=0$ and $p=1$ so that we solve for the coefficients $a_n$ ’s for $n=0,1,2,\cdots$ . The approximation of the solution is found by replacing the coefficients in the solution equation given by:

$u\left(x\right)=\phi \left(0\right)+L^{-1}{\displaystyle \underset{n=0}{\overset{\infty }{\sum }}}{a}_n{t}^n.$ (29)

The advantage of using MADM is that it reduces the computational size of the problem being solved because it involves the calculation of $A_0$ only.

Solution of the Mathematical Model of COVID-19 Using the Modified Adomian Decomposition Method

Applying Equation (27) on each equation of system (1), we obtain the recursive relationship for each equation as follows:

$S_0=S\left(0\right)+L^{-1}\left({\displaystyle {\sum }_{n=0}^{\infty }}{a}_n{t}^n\right)$ ,

$\begin{array}{c}{S}_1=-p{L}^{-1}\left({\displaystyle {\sum }_{n=0}^{\infty }}{a}_n{t}^n\right)-{\rho }_s\left(1-\kappa \omega \right)L^{-1}\left(A_0\right)+L^{-1}\left(\Theta \right)\\ -{\rho }_a\left(1-\kappa \omega \right)L^{-1}\left(B_0\right)-\delta L^{-1}\left(S_0\right),\end{array}$

$S_{n+1}=-{\rho }_s\left(1-\kappa \omega \right)L^{-1}\left(A_n\right)-{\rho }_a\left(1-\kappa \omega \right)L^{-1}\left(B_n\right)-\delta L^{-1}\left(S_n\right)$ , $n\ge 1$ .

$E_0=E\left(0\right)+L^{-1}\left({\displaystyle {\sum }_{n=0}^{\infty }}{a}_n{t}^n\right)$ ,

$\begin{array}{c}{E}_1=-p{L}^{-1}\left({\displaystyle {\sum }_{n=0}^{\infty }}{a}_n{t}^n\right)+{\rho }_s\left(1-\kappa \omega \right)L^{-1}\left(A_0\right)+{\rho }_a\left(1-\kappa \omega \right)L^{-1}\left(B_0\right)\\ -\left(\delta +\sigma \right)L^{-1}\left(E_0\right),\end{array}$

$E_{n+1}={\rho }_s\left(1-\kappa \omega \right)L^{-1}\left(A_n\right)+{\rho }_a\left(1-\kappa \omega \right)L^{-1}\left(B_n\right)-\left(\delta +\sigma \right)L^{-1}\left(E_n\right)$ , $n\ge 1$ .

$I_{s0}=I_s\left(0\right)+L^{-1}\left({\displaystyle {\sum }_{n=0}^{\infty }}{a}_n{t}^n\right)$ ,

$I_{s1}=-p{L}^{-1}\left({\displaystyle {\sum }_{n=0}^{\infty }}{a}_n{t}^n\right)+\left(1-q\right)\sigma L^{-1}\left(E_0\right)-\left(\delta +{\alpha }_1+\gamma +\pi \right)L^{-1}{I}_{s0},$

$I_{s\left(n+1\right)}=\left(1-q\right)\sigma L^{-1}\left(E_n\right)-\left(\delta +{\alpha }_1+\gamma +\pi \right)L^{-1}{I}_{sn}$ , $n\ge 1$ .

$I_{a0}=I_a\left(0\right)+L^{-1}\left({\displaystyle {\sum }_{n=0}^{\infty }}{a}_n{t}^n\right)$ ,

$I_{a1}=-p{L}^{-1}\left({\displaystyle {\sum }_{n=0}^{\infty }}{a}_n{t}^n\right)+q\sigma L^{-1}\left(E_0\right)-\left(\delta +{\alpha }_3+\tau \right)L^{-1}\left(I_{a0}\right)$ ,

$I_{a\left(n+1\right)}=q\sigma L^{-1}\left(E_n\right)-\left(\delta +{\alpha }_3+\tau \right)L^{-1}\left(I_{an}\right)$ , $n\ge 1$

$H_0=H\left(0\right)+L^{-1}\left({\displaystyle {\sum }_{n=0}^{\infty }}{a}_n{t}^n\right)$ ,

$H_1=-p{L}^{-1}\left({\displaystyle {\sum }_{n=0}^{\infty }}{a}_n{t}^n\right)+\pi L^{-1}\left(I_{s0}\right)+\epsilon L^{-1}\left(Q_0\right)-\left({\alpha }_2+\delta +\beta \right)L^{-1}\left(H_0\right)$ ,

$H_{n+1}=\pi L^{-1}\left(I_{sn}\right)+\epsilon L^{-1}\left(Q_n\right)-\left({\alpha }_2+\delta +\beta \right)L^{-1}\left(H_n\right)$ , $n\ge 1$ .

$Q_0=Q\left(0\right)+L^{-1}\left({\displaystyle {\sum }_{n=0}^{\infty }}{a}_n{t}^n\right)$ ,

$Q_1=-p{L}^{-1}\left({\displaystyle {\sum }_{n=0}^{\infty }}{a}_n{t}^n\right)+\gamma L^{-1}\left(I_{s0}\right)-\left({\alpha }_4+\delta +\nu +\epsilon \right)L^{-1}\left(Q_0\right)$ ,

$Q_{n+1}=\gamma L^{-1}\left(I_{sn}\right)-\left({\alpha }_4+\delta +\nu +\epsilon \right)L^{-1}\left(Q_n\right)$ , $n\ge 1$ .

$R_0=R\left(0\right)+L^{-1}\left({\displaystyle {\sum }_{n=0}^{\infty }}{a}_n{t}^n\right)$ ,

$R_1=-p{L}^{-1}\left({\displaystyle {\sum }_{n=0}^{\infty }}{a}_n{t}^n\right)+\beta L^{-1}\left(H_0\right)+\nu L^{-1}\left(Q_0\right)+\tau L^{-1}\left(I_{a0}\right)-\delta L^{-1}\left(R_0\right)$ ,

$R_{n+1}=\beta L^{-1}\left(H_n\right)+\nu L^{-1}\left(Q_n\right)+\tau L^{-1}\left(I_{an}\right)-\delta L^{-1}\left(R_n\right)$ , $n\ge 1$ .

Letting $S_1=E_1=I_{s1}=I_{a1}=H_1=Q_1=R_1=0$ and setting $p=1$ we find the $a_n$ ’s for $n=1,2,3,\cdots$ . and replace in Equation (29) to write the solution for equation system (1).

4. Results and Discussion

The results obtained from solving system Equation (1) using ADM are compared with the results obtained using MADM. We will use COVID-19 data for Zambia [17] as initial values and are given as $S\left(0\right)=1715$ , $E\left(0\right)=54$ , $I_s\left(0\right)=20$ , $I_a\left(0\right)=34$ , $H\left(0\right)=9$ , $Q\left(0\right)=786$ , $R\left(0\right)=85$ . The natural death rate, $\delta$ is calculated by taking the reciprocal of the average life expectancy (in

months). In Zambia the life expectancy is 64.70 years [18] , hence $\delta =\frac{1}{64.70\times 12}$ .

The population of Zambia is $N\left(0\right)=19470000$ [18] . Thus the recruitment Θ is

estimated as $\frac{\Theta }{\delta }=N\left(0\right)$ . The rest of the parameters are estimated from the

literature as given in Table 2. Using the initial values and the parameter values we draw the graphs using SageMaths 9.0 software. The results are shown in Figures 2(a)-(g).

Table 2. Parameter values.

Figure 2. Numerical analysis of MADM and ADM solutions of the model of COVID-19

We can see that the graphs in Figure 2(a) and Figure 2(b), MADM and ADM solutions are close to each other. In Figures 2(c)-(g), MADM solutions are exactly the same as the ADM solutions by coincident. In Figure 2(a), Figure 2(c), Figure 2(d) and Figure 2(f), the population is reducing whereas in Figure 2(b), Figure 2(e) and Figure 2(g), the population is increasing.

5. Conclusion

The model of COVID-19 has been formulated and its dynamic behaviour investigated. We showed that the population classes are non-negative. Using the next generation matrix, we calculated the basic reproduction number, $R_C$ , which is useful in guiding control strategies. By constructing Lyapunov functions, we proved global stability of disease free and endemic equilibrium points. If the basic reproduction number is less than 1, all solutions converge to the disease free equilibrium point and the disease dies out from the population. When the basic reproduction number is greater than 1, the endemic equilibruim point is globally stable, meaning that the disease will persist and the number of infected individuals tends to be a positive constant. Lastly, it is well known that analytical solutions of nonlinear ordinary differential equations are difficult to find. In this paper, we solved nonlinear and linear Ordinary Differential Equations (ODEs) using the ADM and MADM. It can be seen from numerical solutions in Figure 2 that we demonstrated the ability of ADM and MADM in solving ODEs. The series solutions converge very rapidly and from the graphs in Figure 2, we can conclude that the ADM and MADM are very efficient and accurate methods in solving nonlinear mathematical models.

Data Availability

The data used to support the findings of this study is indicated within the article.

Conflicts of Interest

No potential conflict was reported by the authors.

References

[1]	Inyama, S.C., Ekeamadi, G.U., Uwagboe, O.M., Omame, A., Mbachu, H.I. and Uwakwe, J.I. (2019) Application of Homotopy Analysis Method for Solving an SEIRS Epidemic Model. Mathematical Modelling and Applications, 4, 36-48. https://doi.org/10.11648/j.mma.20190403.11
[2]	Abioyei, A.I., Peter, O.J., Ayoade, A.A., Uwaheren, O.A. and Ibrahim, O. (2020) Application of Adomian Decomposition Method on a Mathematical Model of Malaria. Advances in Mathematics: Scientific Journal, 9, 417-435. https://doi.org/10.37418/amsj.9.1.33
[3]	Muhammad, F., Muhammad, U.S., Aqeel, A. and Ahmad, M.O. (2018) Analysis and Numerical Solution of SEIR Epidemic Model of Measles with Non-Integer Time Fractional Derivatives by Using Laplace Adomian Decomposition Method. Ain Shams Engineering Journal, 9, 93391-3397. https://doi.org/10.1016/j.asej.2017.11.010
[4]	Nigar, A., Saeed, A., Sartaj, A. and Gul, Z. (2019) The Adomian Decomposition Method for Solving HIV Infection Model of Latently Infected Cells. Matrix Science Mathematic, 3, 05-08. https://doi.org/10.26480/msmk.01.2019.05.08
[5]	Ali, A.A., Rahmamn, M., Shah, Z. and Kuman, P. (2022) Investigation of a Time-Fractional COVID-19 Mathematical Model with Singular Kernel. Advances in Continuous and Discrete Models, 2022, Article No. 34. https://doi.org/10.1186/s13662-022-03701-z
[6]	Enahoro, I., Oluwaseun, O.S., Ngonghala, C. and Gumel, A.B. (2020) Mathematical Modeling and Analysis of COVID-19 Pandamec in Nigeria. Mathematical Biosciences and Engineering, 17, 7192-7220. https://doi.org/10.3934/mbe.2020369
[7]	Worldometer. https://www.worldometers.info/
[8]	Kounakep, Y.T., Tchoumi, S.Y., Fotsa, D.J., Kamba, F.G.T., Ngounou, D., Mboula, E., Kamla, V. and Kamgang, J.C. (2021) Modelling the Anti-COVID-19 Individual or Collective Containment Strategies in Cameroon. Applied Mathematical Sciences, 15, 63-78. https://doi.org/10.12988/ams.2021.914395
[9]	Worldometer (2020) Coronavirus Incubation Period. https://www.worldometers.info/coronavirus
[10]	Brauer, F., van den Driessche, P. and Wu, J. (2008) Mathematical Epidemiology. Springer-Verlag, Berlin. https://doi.org/10.1007/978-3-540-78911-6
[11]	LaSalle, J.P. (1976) The Stability of Dynamical Systems. Regional Conference Series in Applied Mathematics. SIAM, Philadephia.
[12]	Xu, D.G., Xu, X.Y., Yang, C, H. and Gui, W.H. (2015) Global Stability of a Variation Epidemic Spreading Model on Complex Networks. Mathematical Problems in Engineering, 2015, Article ID: 365049. https://doi.org/10.1155/2015/365049
[13]	Lyapunov, A.M. (1992) The General Problem of the Stability of Motion. Taylor and Francis, London.
[14]	Adomian, G. (1983) Stochastic Systems. Academic Press, New York.
[15]	Adomian, G. and Rach, R. (1983) Inversion of Nonlinear Stochastic Operators. Journal of Mathematical Analysis and Applications, 91, 39-46. https://doi.org/10.1016/0022-247X(83)90090-2
[16]	Jafar, B. and Hosseini, K. (2016) A Modified Adomian Decomposition Method for Singular Initial Value Emden-Fowler Type Equations. International Journal of Applied Mathematical Research, 5, 69-72. https://doi.org/10.14419/ijamr.v5i1.5666
[17]	Ministry of Health and Zambia National Public Health Institute (2022) Zambia COVID-19 Statistics Daily Status Update.
[18]	https://www.macrotrends.net>ZMB.
[19]	Mnganga, J.M. and Zachariah, N.S. (2020) Mathematical Model of COVID-19 Transmission Dynamics and Control Strategies. Journal of Applied & Computational Mathematics, 9, 1-3.
[20]	Camilla, R., Mirjam, S., Peter, S., Gisela, B., Guenter, F., Claudia, W., Thorbjörn, Z., Verena, T., Christian, J., Wolfgang, G., Michael, S., Christian, D., Patrick, V., Katrin, Z., Sabine, Z., Roman, W. and Michael, H. (2019) Transmission of 2019-nCoV Infection from an Asymptomatic Contact in Germany: Case Report. The New England Journal of Medicine, 382, 970-971. https://doi.org/10.1056/NEJMc2001468
[21]	Matiur, R., Saeed, A., Matoog, R.T., Nawal, A.A. and Tahir, K. (2021) Study on the Mathematical Modelling of COVID-19 with Caputo-Frabrizio Operator. Chaos, Solitons & Fractals, 150, Article ID: 111121. https://doi.org/10.1016/j.chaos.2021.111121
[22]	Duah, D., Iddi, S., Adu, B., Aheto, M.J., Kojo, M.S., Fobil, J. and Bosomprah, S. (2021) Mathematical Modelling of COVID-19 Infection Dynamics in Ghana: Impact Evaluation of Integrated Government and Individual Level Interventions. Infectious Disease Modelling, 6, 381-397. https://doi.org/10.1016/j.idm.2021.01.008
[23]	Pakwan, R., Sherif, E.S. and Irthit, I. (2021) A Mathematical Model of COVID-19 Pandemic: A Case Study of Bangkok, Thailand. Computational and Mathematical Methods in Medicine, 2021, Article ID: 6664483. https://doi.org/10.1155/2021/6664483
[24]	Chernet, T.D. and Gemechis, F.D. (2020) Modeling and Optimal Control Analysis of Transmission Dynamics of COVID-19: A Case of Ethopia. Alexandria Engineering Journal, 60, 719-732. https://doi.org/10.1016/j.aej.2020.10.004
[25]	Mushanyu, J., Chazuka, Z., Mudzingwa, F. and Ogbogbo, C. (2021) Modelling the Impact of Detention on COVID-19 Transmission Dynamics in Ghana. RMS: Research in Mathematics & Statistics, 8, Article ID: 1953722. https://doi.org/10.1080/27658449.2021.1953722