Impact of pangolin bootleg market on the dynamics of COVID-19 model.

In this paper we consider ant-eating pangolin as a possible source of the novel corona virus (COVID-19) and propose a new mathematical model describing the dynamics of COVID-19 pandemic. Our new model is based on the hypotheses that the pangolin and human populations are divided into measurable partitions and also incorporates pangolin bootleg market or reservoir. First we study the important mathematical properties like existence, boundedness and positivity of solution of the proposed model. After finding the threshold quantity for the underlying model, the possible stationary states are explored. We exploit linearization as well as Lyapanuv function theory to exhibit local stability analysis of the model in terms of the threshold quantity. We then discuss the global stability analyses of the newly introduced model and found conditions for its stability in terms of the basic reproduction number. It is also shown that for certain values of R 0 , our model exhibits a backward bifurcation. Numerical simulations are performed to verify and support our analytical findings.

Introduction

Corona virus disease (COVID-19) is the current outbreak, which is caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). This disease was reported for the first time in the Wuhan city of China. The respiratory track of the infected person is badly affected. Common symptoms of this disease include difficulty in breathing, dry cough, sever pain and fever. This horrible outbreak has infected more than three million people all over the world [1], [2], [3], [4]. The virus is zoonotic in nature, i.e., it can be transmitted from animals to human and cause disease. Several studies indicate that ant-eating pangolin (Manis Pentadactyla) is the common vector of CoVID-19 [1]. Similar studies have also been recorded in the past, the pandemic of SARS Virus back in 2003 which had been transmitted from Bats to a Civet [3]. Some data shows that the DNA sequence of the Corona virus discovered in the lung sample of pangolin was found to be similar to SARS CoV-2, known as CoVID-19 and both of these viruses have specific type of proteins on their surface known as spike proteins [5], [6].

Pangolins are amongst those animals that are heavily trafficked in the world, which has resulted in the extension of this animal. After the out break of Covid-19 pangolins are considered to have a great link with its initial out break in China. Though the evidences are not adequate, still Chinese government has propagated this message and is working to pass legislation for banning the use of this animal. More attention is required for the conservation of pangolins; therefore, governments need to take more serious actions against the wildlife trade. The main question that arises about Covid-19 is to know from where it emerged? Evidence shows that several animals are the source of Covid-19. Addressing this question is very essential for public health authorities and policy makers because it will help us to prevent out breaks in future. Many theories have been initiated, but potentially the most noticeable cause is pangolin. There are two reasons behind illegal trade and hunting of pangolin. The first reason is the delicacy of its meat and is mostly used in China and Vietnam. The second one is that pangolin scales is mostly used in traditional remedies of China [11]. Due to these two main factors pangolin is considered as the world’s most trafficked animal. In press conference held on 7th February 2020 by ‘South China Agriculture University, Guangzhou, it was stated that Covid-19 came from pangolin.

Pangolin is genetically similar to human’s corona virus at 99.1 per cent where it is related to a specific site known as receptor binding domain (RBD) [7]. The pangolin feeds on the excreta and saliva of the bats, which is the source of Corona virus and due to this feeding pangolin was found to be infected by the virus. It has been revealed from the whole genome sequences comparison that pangolin shows similarity of 90.3 with human corona virus. Corona virus inserts its crucial part, known as RBD, to human cell. There is 99.1 percent similarity in their RBD of the viruses [7]. Tao et al. (2020) reported that pangolin acts as a source of CoVID-19 [8]. The newly identified corona virus is interestingly found to be the most trafficked disease and would be a threat to public health if this illegal trafficking across the world is not controlled [9]. Illegally trafficked pangolin share 85.5 to 92.4 percent similarity of their DNA with virus in human being found in the frozen cell samples of the pangolin[10]. The anti-smuggling agency in Malaysia also found corona virus in cells of frozen tissues of 18 Malayan pangolin, where corona virus RNA was found is six of 43 organ samples [12]. Other studies also identified SARS-CoV-2-related viruses in Malayan pangolin, native in Southeast Asia, that were smuggled into Southern China [13], [14], [15], [16]. Several biological homology observed a very minor difference between the scaly ant-eater pangolin corona virus and SARS Cov-2. There was also similarity among Bat Corona virus, Rat G-13, and SARS Cov-2. These findings suggest that scaly ant-eater pangolin acts as intermediate hosts of SARS Cov-2 for the transmission of Bat corona virus (CoVID-19) [12].

Many people who are working in various disciplines of science and technology, struggling continuously to curb or at least control this global killer. In the thick of these wide spread efforts mathematical modeling also play a big role. Mathematical modeling has attracted the attention of considerable number of researchers in the past few decades, see, e.g., [29], [30], [31]. Various infectious diseases have been studied through mathematical models, e.g., see [32], [33], [34], [35], [36], [37], [38]. With the aid of mathematical models, researchers first try to understand the spread of an infectious disease and then develop strategies to eradicate or control it from a community, see, e.g., [39], [40], [41], [42].

Like other infectious diseases, various models can be found in literature which describe the transmission dynamics of the current pandemic [43], [44], [45], [46], [47], [48]. Several factors have been taken into account for the possible spread of the COVID-19. To the best of our knowledge, the impact of pangolin bootleg market on the dynamics of COVID-19 has not been reported yet.

In this article, we take this possibility into account and develop a novel mathematical model which may be helpful to understand the transmission of the diseases between pangolin and humans. It will also play a role in estimating the impact of controls at both national and global level and to integrating and interpreting scrapie monitoring results. The paper is organized as follows. In Section, we present our model and describe its basic mathematical properties. We then, in Section “Dynamical aspects”, discuss the two possible non-trivial equilibrium states and compute the threshold quantity for the model under consideration. In the same section, we also study the stability analysis of both states of the model and discuss the backward bifurcation. We perform numerical simulations in Section “Numerical simulations and discussion” and discuss our derived results in detail. Finally, we conclude our work in Section “Conclusions”.

Formulation of the model

Mathematical model used in this research paper is an attempt to discuss the transmission of the novel Corona virus through ant-eating pangolin. How does this pandemic disease impact the population? The World Health Organization (WHO) has declared the COVID-19 as a public emergency and is searching for the root cause of its spread. The current study is an attempt to answer some of the basic questions regarding the current pandemic. We consider ant-eating pangolin as a vector to the current virus.

The infected pangolin with SARS-CoV-2 virus infects human population in bootleg market. These suspected humans, when interact with other susceptible humans in the community, spread the disease. Thus in our model, we intend to focus our concentration on two types of interacted population, namely, the human and ant-eating pangolin.

All the parameters used in the model are suitably designed to be non-negative and represent specific rates and periods. The description of the parameters involved in our model is given in Table 1 . In the scenario of Table 1, mathematical model for the current problem can be casted in the form of the following autonomous differential equations:

Table 1

Variables and parameters description of the proposed PH-model described by (1).

Variables (Parameters)	Description
P_s(t)H_s(t)	Class of susceptible pangolin [human]
P_e(t)H_e(t)	Class of exposes pangolin [human]
P_i(t)H_i(t)	Class of infected pangolin [human]
P_R(t)H_R(t)	Class of removed pangolin [human]
H_A(t)	Class of asymptomatic human
M(t)	Pangolin bootleg market or reservoir
n_pn_h	Birth rate of pangolin (human)
N_pN_h	Total number of pangolin (human)
m_pm_h	The natural mortality rate of pangolin (human)
γ_pγ_h	Transmission rate from P_iH_i to P_sH_s
ρ_p=n_pN_p,ρ_h=n_hN_h
T_p	Transmission rate of P_e to P_i
1T_p	The incubation period of pangolin
Σ_p	Transmission rate of P_i to P_R
1Σ_p	The infectious period of pangolin
0⩽κ≤1	The transmissibility of H_A that of H_i
β_p′	The transmission rate from black market M to H_s
T_h=η_hf_h	where
η_h	The symptomatic infection rate proportion of human and
1f_h	is the incubation period of human
T_h′=η_h′f_h′	where
η_h′	The asymptomatic infection rate proportion of human and
1fh′	is the period between infection to symptoms known as latent period of human
1∑_h	The symptomatic infection period of human
1∑∑_h	The asymptomatic infection period of human
C_h	The exfoliating coefficient from H_i to bootleg market M
C_h′	The exfoliating coefficient from H_A to bootleg market M
β_p	The peddle rate of pangolin in black market
ε	The wastage rate of pangolin in M compartment
1ε	The life span of infected virus in M compartment

We will investigate our model (1) under the following biologically feasible initial conditions and

It is worth noticing that the proposed system (1) satisfies the following conservative laws

Note that and , where and denote the total population size of pangolin and human being respectively at time t. We observe that Eqs. (3), (4) are exact differential equations with respective solutions given by,

In the following we discuss some basic properties of our model. The first result is stated and proved in the following theorem.

For non-negative initial conditions (2) , all the feasible solutions of (1) are bounded and are not identically trivial on any interval and enter the region, .

Assume the set be any solution of the system (1) subject to the non negative initial conditions  ,  as well as and for . Similarly for any non negative initial values of , and the solution is bounded below by zero. By adding equations of the system involving pangolin classes in (1) we haveSolution of this equation has the form with , the initial pangolin population. From this equation, it is easy to deduce that when the time grows without bound.

Moving on the same line and considering classes involving human population in the proposed model one can obtain

If represents the initial human population then solution of the last equation looks like . Again, it is concluded from this relation that as . Therefore, all solutions of the system (1) enter the feasible region,

, where , and . Thus, one may conclude that is the largest positive invariant set. □

Dynamical aspects

There are several aspects of a dynamical system which provide information about the complex nature of a system. Among these features equilibrium points are the most import. Equilibrium points help us to define states corresponding to constant operating conditions in a dynamical system. Depending upon nature of a dynamical system, there may exist zero, one or more than one equilibrium points. This section is devoted to focus our attention on determining stationary states of the model (1) and to perform the stability analysis of the two equilibrium states. First we consider the situation when there is no infection of the disease in a community, i.e., the so-called disease free equilibrium. In the following we denote such equilibrium by . The second possible state is that one when the disease exists in a community. This state is known as the endemic state and will be denoted by .

Inserting , and in the given model and solving the system of autonomous differential Eq. (1) for and we obtain

Before performing the stability analysis of the both the disease free and endemic equilibria of the model, we consider it necessary to compute the so-called basic reproduction number, which plays a key role in the stability analysis.

Spectral radius theory

The basic reproductive number or spectral radius is calculated by exploiting the traditional next generating matrix method [17]. The threshold quantity is obtained by considering the subsystems, the state-of-infection , and the state-of-infectiousness , of (1) following [18]. In the matrix form the dynamical system (1) can be written as , where , andand

The Jacobians of the above matrices F and V are, respectively calculated as

Note that the matrix V is non-singular and its multiplicative inverse has the form

With the help of the matrices (9), (11), the next generating matrix of the proposed system is calculated as follows,

The partial reproductive numbers of ant-eating pangolin and human population are defined by, and .

Thus, the next generation matrix isHere . Following [19] the basic reproductive number is the product of the two partial reproductive numbers. Thus we have

The spectral radius of the matrix gives the threshold quantity [19]. Following [19], we conclude that whenever then each eigenvalue of the next generation matrix has negative real part and the disease free equilibrium is locally asymptotically stable. see Fig. 1 .

Fig. 1

Schematic diagram (Flow Chart) of the proposed model.

Plots of the basic reproduction number in terms of some parameters of the model under consideration has been depicted in Fig. 2 .

Fig. 2

Profiles of the threshold quantity , in terms of various parameters involved in model (1) under consideration.

Local stability at the disease free equilibrium state

The purpose of this subsection is to show the local stability analysis of the underlying model. In fact if we can show that the system (1) is locally stable at disease free as well as at disease present equilibrium points, then the commutative epidemiological situations are not too much dissimilar from that particular equilibrium point that will transform into that particular equilibrium point. It will also give us information that the equilibrium points posses the ability to be stable in response to small perturbations: i.e., to say that if the situation is moved away from the equilibrium point, then the system will resume its previous position. For the underlying model (1) we provide the result related to local stability at the disease free equilibrium state in the form of the following theorem.

The disease free equilibrium of the model (1) is locally asymptotically stable for , otherwise unstable.

The Jacobian matrix for the local stability of the disease free equilibrium state can be examined by linearizing the model (1) around the disease free state as

We perform the traditional row operations over the above matrix to obtainwhere , , , and .

The characteristics equation of in (16) has the formwhere ,  , and .

There are ten eigenvalues of (17), among them two eigenvalues are and with multiplicity 2. The eigenvalues and , being negative, have negative real parts. The other three eigenvalues of (17) can be obtained by solving

Following Routh-Hurwitz Criterion, [20], the three eigenvalues will have negative real parts provided for with and .

For we obtain for . Thus all the eigenvalues of the characteristic equation have negative real parts if and only if . This clearly shows that disease free equilibrium is locally asymptotically stable.

For or equivalently , we have and limit , when , there exists such that . From here one can deduce that the disease free equilibrium state is unstable.

Global stability of the disease free equilibrium

In this subsection of the manuscript, we study the global properties of the underlying Covid-19 model. The existence of locally stability of the system (1) depends upon the initial conditions but the global stability of a system is regardless of the subsidiary initial conditions. The local stability at a particular point suggests that situation will exist for a particular short interval of time, while the global stability of a system infers the occurrence of that particular situation forever regardless of any subsidiary initial or boundary conditions. Disease prevails if we are able to show that the system is globally stable at endemic equilibrium point. Moreover, we can get the information that the a disease dies out if the system is globally stable at disease free equilibrium point. Therefore, in this part of the paper we provide results about the global stability analysis of model (1). The global stability at disease free equilibrium state of the proposed model (1) can be seen from the following theorem.

The disease free equilibrium of the system (1) is globally asymptotically stable in the region if , otherwise unstable.

In order to establish global stability of the model (1) at the disease free equilibrium point , we construct the following Lyapunov functionwhere for are non negative constants that will be determined later on. Calculating the time rate derivative of (19) along with solutions of the provided system (1), we have

Using and in (19), one may arrive at

By choosing and , the last equation can be seen in the form,

Thus, if and only if , and . Also is negative whenever , and . In other words the function is negative semi-definite. Thus following the Lascalle’s invariant rule [21], the largest invariant compact set is a singleton set of disease free equilibrium point. Thus, is globally asymptotically stable in . □

Endemic equilibrium point and backward bifurcation

The system (1) has endemic equilibrium point where the infected individuals in each class are non-zero. We need to take the following necessary steps.

Assume that be an arbitrary endemic equilibrium state of the model under investigation (1). Solving the autonomous equations of the model(1) simultaneously at steady state gives,

Assume , and substituting , in the seventh equation of the endemic equilibrium point of the system (1) at the static state. After some algebra we obtain the following quadratic equationwhere

It can be noticed that the coefficient u is positive. Since, , then w can either be positive or negative for certain values of parameters involved in w. Therefore, Eq. (22) can have a positive root. This implies existence of the endemic equilibrium point. Finally at the left of there exists an interval such that

Now if and or , Eq. (22) has no positive solution and there exists no endemic equilibrium state at all. This whole colloquy shows that the equilibria of Eq. (22) continuously depends upon the value of . Therefore, we establish the following results for the different values of coefficients of Eq. (22).

The model (1) satisfies the following facts:

The PH model (1) has unique endemic equilibrium, if and .

The PH model (1) has two endemic equilibria on , if and .

PH model (1) has no endemic equilibria on , in other case.

Now to explore the phenomenon of backward bifurcation [22], [23], we set the discriminant of the quadratic relation in equal to zero i.e . The critical point for , has the form

It follows that backward bifurcation occurs when or equivalently . The epidemiological significance of backward bifurcation needs a basic requirement of . This condition is necessary but not sufficient to eradicate the pandemic of COVID-19. In some scenario,the pandemic elimination needs the initial state parameters and some time initial size of the population of a system. The phenomenon of backward bifurcation in model (1) recommends biologic feasibility of control of the pandemic whenever the threshold quantity and also depends upon the initial population of the model (1).

The PH model (1) observe backward bifurcation if .

The PH model (1) goes through backward bifurcation when if and only if .

For the sufficient condition, consider the graph of . For one may set for specific real values of parameters in w. Thus the graph of passes through for . Further, if , we have has a root namely . Now when we increase w to a non zero positive value and , which guarantees the exists of an open interval containing w where has two real and positive roots. In short words, this shows there must exists two endemic equilibria if . If , Eq. (22) has no positive real roots if . □

Local stability of endemic equilibrium state

Following theorem is enough to show that the system(1) performs local stability at the endemic equilibrium point .

The unique endemic equilibrium state of system (1) is locally asymptotically stable if , and is unstable if the threshold quantity is less than unity.

The Jacobian matrix for the unique endemic equilibrium state of the PH model (1) is given bywhere . To discuss nature of the eigenvalues of the matrix (26), we make elementary row operations to obtain

In the above matrix (27) we have assumed

The matrix in (27) has the following eigenvalues

The eigenvalue if and only if the condition holds.

Inserting values of in and , the given condition takes the form and arrive at

We note that all the coefficients of Eq. (28) are positive if . Thus from the above stated theorem it is concluded that the system (1) has local asymptotic stability at endemic equilibrium point [24], [25]. □

Global stability of endemic equilibrium state

In this section we explore the global stability of the endemic equilibrium point of the said model (1) in terms of the basic reproduction number .

The endemic equilibrium state of PH model (1) is globally asymptotically stable on , if and

In order to study the global stability of the PH model(1) under analysis we define the following Lyapunov function,

Differentiating with respect to t and using values from (1), we arrive at

Using the pre-defined values and making rearrangement, we arrive at

Since the geometric mean is always less than or equal to the arithmetic mean, we have

Thus, the defined set of parameters ensures that for all

Also , if . Then the equilibrium state is the only positive invariant set for the set of equations contained entirely in

We conclude that the non-negative endemic equilibrium state is globally asymptotically stable on [21]. □

Numerical simulations and discussion

In this section, we perform numerical simulations to verify the analytical findings in the previous section. We exploit the RK4 method to obtain the numerical results. Values of the parameters involved in the model under consideration are as follows [26], [27], [28].

, , One can easily verify that, for the above defined parametric values, the value of the threshold parameter is greater than unity, showing the existence of the disease in a community.

The first figure in Fig. 3 shows dynamics of the susceptible pangolin compartment against the time variable for different initial conditions. One can observe that this compartment is always positive, showing the existence of susceptible pangolin. At the beginning there is a sharp decay in this compartment and finally observes stable as time progresses. In second panel of the same figure, the exposed class of pangolin is depicted against time. There is a rapid increase in the exposed population and again the solution curves become stable as t increases. In the third panel, the infected class of pangolin is plotted versus time. One can observe increase this class initially and stability afterwards. The removed class of pangolin is depicted in the final panel of the same figure. This removal is either due to the disease related death or due to recovery from the disease. Fig. 4 shows the transmission dynamics of pangolin bootleg market and asymptomatic human for various initial conditions of the concerned classes. Stability of the integral curves for each compartment can be observed from the figure.

Fig. 3

Transmission dynamics of the susceptible, exposed, infected and removed pangolin compartments for various initial conditions. Values of the parameters involved in these simulations are given in Section “Numerical simulations and discussion”.

Fig. 4

Profiles of the pangolin bootleg market and Asymptomatic human classes for various initial conditions of the concerned compartment. The detail of the parametric values is given in Section “Numerical simulations and discussion”.

Fig. 5 shows the profiles of the susceptible, exposed, infected and removed human classes for different initial conditions. As before, we observe the biological feasibility of the solution curves. All the integral curves go to the stability point confirming the global stability of each compartment for the disease endemic situation.

Fig. 5

Time dynamics of the susceptible, exposed, infected and removed human compartments for various initial conditions. As before, the parametric values used in the numerical simulations are given Section “Numerical simulations and discussion”.

Conclusions

On 31st December, 2019 a global killer, Covid-19 emerged from Wuhan city of China and spread all over the world like tsunami. There is no clear cause and neither vaccine/treatment is discovered. Evidences showed that the disease can be transmitted from animals to human as well as it can be horizontally transmitted from one person to another. To date now it is discovered that the symptoms include, fever, huge coughing and difficulties in breathing. Apparently, it is found that the appearance of symptoms lies in the range 2–10 days. Segregating the pangolin and human population into mutually disjoint classes as well as incorporating pangolin blotting market or reserver, we formulate a novel Covid-19 disease model. Initially, it is shown that the solution set to the underlying system of ordinary differential equation governing Covid-19 phenomena, lies in some feasible region i.e to say that it lies below a limit and also it is above some extent. After expounding different equilibrium points we have shown that the underlying system is locally stable at disease free and disease present equilibrium points. These results suggest that if the system is disturbed in a particular interval, then the system have ability to resume its actual position. Also, it is shown that our model is globally stable under some auxiliary condition. These results can be interpreted that under some hypotheses and auxiliary conditions the disease can die out and also shown that when the disease will prevail in a society. It is worth noticing that if we use fewer hypothesis in theorem (3.8), then it needs to see whether the global stability can be achieved or not.

Authors contribution

“All authors contributed equally”.

CRediT authorship contribution statement

Abd Ullah: Formal analysis, Investigation, Methodology, Project administration, Visualization, Supervision, Writing - original draft, Writing - review & editing. Saeed Ahmad: Formal analysis, Investigation, Methodology, Software, Writing - original draft, Writing - review & editing. Ghaus ur Rahman: Formal analysis, Investigation, Methodology, Supervision, Writing - original draft, Writing - review & editing. M.M. Alqarni: Conceptualization, Formal analysis, Investigation, Methodology, Resources, Writing - original draft, Writing - review & editing. Emad E. Mahmoud: Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Validation, Visualization, Supervision, Writing - original draft, Writing - review & editing.

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.