A mathematical model of the evolution and spread of pathogenic coronaviruses from natural host to human host.

Coronaviruses are highly transmissible and are pathogenic viruses of the 21st century worldwide. In general, these viruses are originated in bats or rodents. At the same time, the transmission of the infection to the human host is caused by domestic animals that represent in the habitat the intermediate host. In this study, we review the currently collected information about coronaviruses and establish a model of differential equations with piecewise constant arguments to discuss the spread of the infection from the natural host to the intermediate, and from them to the human host, while we focus on the potential spillover of bat-borne coronaviruses. The local stability of the positive equilibrium point of the model is considered via the Linearized Stability Theorem. Besides, we discuss global stability by employing an appropriate Lyapunov function. To analyze the outbreak in early detection, we incorporate the Allee effect at time t and obtain stability conditions for the dynamical behavior. Furthermore, it is shown that the model demonstrates the Neimark-Sacker Bifurcation. Finally, we conduct numerical simulations to support the theoretical findings.

Introduction

According to the International Committee on Taxonomy of Viruses (ICTV), coronaviruses are members of the sub-family Coronavirinae in the family Coronaviridae and the order Nidovirales. Coronavirinae consists of four genera groups; α-coronavirus, β-coronavirus, γ-coronavirus, and δ-coronavirus [1,2]. Recently, these groups are divided in terms of phylogenetic clustering while before they were sorted based on serology. All viruses of Nidovirales order are enveloped, non-segmented positive-sense RNA viruses, where within this, the coronaviridae has the most significant identified RNA genomes, containing approximately 30 kgbases (kb) genomes. On the other hand, all coronaviruses have animal origin [2,3]. Table 1 below shows that the four genera of coronaviruses originated from animals.

Table. 1

Genera of CoV and the pathogenic class.

Coronavirinae Genera	α-CoV	β-CoV	γ-CoV	δ-CoV
Pathogenic Class	mammals	mammals	both non-mammal and mammals	both non-mammal and mammals

The HCoV-NL63 and HCoV-229E are coronaviruses that cause mild infections in humans. On the other hand, SADS-CoV, which has swine as an intermediate host role, does not show any evidence of infections in humans. HCoV-OC43 and HCoV-HKU1 are both coronaviruses and are also mostly harmless to the human body that has a rodent-borne origin. HCoV-229E and HCoV-OC43 were isolated nearly 50 years ago, while HCoV-NL63 and HCoV-HKU1 were identified in 2003. Coronaviruses have not considered as highly pathogenic until the outbreak of SARS-CoV in 2003 and MERS-CoV in 2012. The spread of SARS-CoV in China (Guangdong) indicated that a coronavirus was transmitted from bats to an intermediate host like market civets, and from there to the human host, while the outbreak of MERS-CoV in the middle east countries also came from bats to dromedary camels as an intermediate host, and then, it was transmitted to human [4], [5], [6], [7], [8]. These viruses cause respiratory and intestinal infections, including fever, dizziness, and cough.

On the 12th of December 2019, a new virus form of Coronaviridae was reported in China (Wuhan). The outbreak was associated again with intermediate hosts like reptilians, while the natural host was assumed as bats. This novel virus was designated at first as WH—Human 1 and was referred after that as COVID-19 by the WHO. COVID-19 was characterized by two members of β-coronavirus; human-origin coronavirus (SARS-CoV Tor2) and bat-origin coronavirus (bat-SL-CoVZC45), while intensive studies show that it was most closely related to the bat-origin coronavirus [9]. As a result of the above discussion, the primary assumption was set as the natural host of COVID-19 is the bat, which infects the human population with a domestic intermediate host. Fig. 1 shows the animal origins of the human coronaviruses SARS-CoV, MERS-CoV, and COVID-19.

Fig. 1

Natural-Intermediate and Human Host transmission.

Explicitly stating, if we analyze the environmental origin of COVID-19, then we adopt the hypothesis that says that infected bats may spread the disease. These bats of genus Rhinolophus are mainly in the area of Shatan River Valley. Domestic animals like snakes in that area are hunted for the food market in Wuhan, which had an intermediate host role in the transmission. Finally, this virus spillover from the intermediate hosts to cause several diseases in human. A virus that started with an endemic pathogenic behavior in China (Wuhan) reaches somehow to a pandemic point worldwide.

In this study, we give brief information about coronaviruses and introduce the genera of Covid-19. We explain the transmission from the natural host to the intermediate one, and from there to the human host. The mathematical model in Section 2 shows the spread from animal to human and human to human. We noticed in our model that the transmission risk from human to human is higher and more dangerous than the effect from animal to human. In the end, the study reveals that the infected class who do not know they are infected (because of late recognized symptoms) is the major problem of a pandemic spread worldwide.

The model

Many research areas in biology or medicine are attractive topics for scholars engaging in applied mathematics since mathematical modeling has an essential role in understanding the dynamics of many diseases and biological phenomena. Over the years, biological models have been formulated mathematically [10], [11], [12], [13]. Mainly, studies are restricted to integer-order differential equations. However, it is seen that many problems in biology, as well as in other fields like engineering, finance, and economics, can be formulated successfully by differential equations with piecewise constant arguments [14], [15], [16], [17]. In mathematical modeling, for an overlapping species population, it is convenient to use differential equations; on the contrary, for a non-overlapping species one, it is suitable to use difference equations [14]. However, there are some dynamics in the environment, which combine the tools of both differential-difference equations concerning time t. For such biological events, it is preferred to construct models of differential equations with both continuous and discrete-time. Relevant studies that have considered the discrete and continuous-time effects as differential equations with piecewise constant arguments are given in [17], [18], [19], [20], [21].

The question of whether the coronavirus is used as a biological weapon or not is out of our primary objective in the paper. However, this question was under consideration for the plague disease [22].

We consider here the pandemic infection that occurs when the virus is transmitted to the human body from the intermediate host and continues to spread from human-to-human. The first three equations on the system show an SI (susceptible-infected) model to explain the transmission from human-to-human, where S is the susceptible class, C 1 is the infected class, which does not know that they are infected because of the late occurred symptoms of COVID-19 and C2 represents the infected class that knows they are infected. The spillover from the intermediate infected class M to the human host S denotes a predator-prey mathematical model, while for the transmission from the natural host N, which is the bat population, to intermediate host M is a host-parasite model of Holling Type II. Thus, the mathematical model of this biological phenomena is modeled as follows;whererepresent the Holling type II functions. All the parameters in (1) belong to ℝ and is the integer part of t ∈ [0,  ∞).

The susceptible S is composed of individuals that have not contacted the infection but can get infected through contact with the humans that do not know they are infected and from the intermediate hosts. r1 is the population growth rate of the susceptible population and α 1 denotes the logistic rate. p is a rate of the susceptible population per year. The susceptible lost their class following contacts with infectives C 1 and the intermediate host M at a rate β 1 and β 2, respectively. σ 1 shows the parameter of the interaction between the hunted M class and the predator S population.

The C 1class does not know that they have COVID-19. In this equation, r2 is the population growth rate of the class, while α 2 is the logistic rate. The population of this class decreases after screening at a rate θ and be aware of the infection. Another possibility is that after the S-C1 contact, the symptoms occur in early stages so that both classes noticed that they are infected, which is given with the rate ɛ1. The intermediate host infected group could also show early symptoms to be aware of the infection, which is given by a rate of ɛ2. The logistic rate of C2 is denoted as α 3.

M is the domestic animal as an intermediate class in the corona transmission spread. r3 is the intrinsic growth rate of the population, while α 4 is the logistic rate. σ 2 shows the effect on the hunted M during the interaction between the intermediate host and susceptible class. γ denotes the predation rate in the host-parasite scheme.

N represents the natural host (bat population) of COVID-19 in this dynamic system. r4 is the intrinsic growth rate and α 5 is the logistic rate of the population. δ shows the conversion factor of the natural host. e is the attack rate of the bat population to infect the M, while ω (0 < ω ≤ 1) represents the fraction of the potential infectivity of the natural host. h is the rate of average time spent on infecting the domestic intermediate class, which is also known as the handling time.

Local and global stability analysis

Herein, we investigate the local and global stability of the system (1). Before proceeding to the main result, we need some preparations.

Integration of system (1) on an interval of leads towhich means that for positive initial conditions, the solutions of (3) are positive as well. Moreover, on an interval of we can write system (3) aswhich is a system of Bernoulli equations. Integrating both sides of (4) concerning t on [n, t) and taking we get a difference equation system such asTherefore, we conclude that any global analysis of (5) represents the behavior of (1), since (5) is a solution of (1) for In this section, we want to consider the local and global stability of the co-existing equilibrium point that represents the positive equilibrium point of system (5). The Jacobian matrix of (5) around Λ is given bywhere and Thus, the characteristic equation of (6) iswhereWe need the following theorem to prove the local stability.

Theorem 1. Linearized Stability Theorem [23]

Letwhere for p, q ∈  ℝ+ the characteristic equation isand the initial conditions are  ∈  ℝ+.

If |λ 1, 2| < 1 then the equilibrium of Eq. (9) is locally asymptotically stable. Thus, |λ 1, 2| < 1 if and only ifTheorem 2. Let Λ be the positive equilibrium point of system (5). Assume that the basic reproduction numbers are R 01 < 1 and R 02 < 1. Furthermore, let,, ɛ1 < 1 and

. If andwhere , and then the positive equilibrium point of system (5) is locally asymptotically stable.

Proof. Let us consider at first (8), where we obtainsince . Thus, we need to consider the quadratic equations given byandThe characteristic Eq. (13) can be rewritten in a form ofwhich implieswhereis the basic reproduction number, that shows the transmission potential of the class. In applying the Linearized Stability Theorem to (15), we obtainFromwe getwhere R 01 < 1, and for ɛ1 < 1. On the other side, consideringwe havewhere.

From (18) and (19), we obtainwhere , andAdditionally, the characteristic Eq. (14) can also be rewritten such aswhich implieswhereis the basic reproduction number of the intermediate-natural host classes. From the Linearized Stability Theorem, we want to consider the conditions for the given inequalityFromwe havewhere R 02 < 1, and Moreover, from (24), we obtainwhich holds for the inequalitiesandand Considering (26) together with (27) and (28), we obtainandwhere and . This completes the proof. □

Remark 1. Theorem 2 shows that among the human hosts, those who do not know they are infected, are the control class in the spread. In contrast, between the animal hosts, the intermediate class plays a dominant role, since that one has the essential role in transmitting the virus from animal to human. The transmission potential for both and are R 01 < 1 and R 02 < 1. Moreover, the susceptible class and the C 1 class is stable based on two parameters, which are the awareness of the symptoms and the screening rate. We noticed that class C 1 should be more aware of the symptoms that might become from the susceptible class as well as from the intermediate class, than the S class to stop the outbreak. For the susceptible class, it is more important to keep the population rate per year non-infected. The transmission of the virus to the offspring would reach an uncontrollable phenomenon worldwide.

Theorem 3. Let be a positive solution to the system (5). Then the following statements are true.then is increasing monotonically.then is decreasing monotonically.

If

If

Proof.This completes the proof.

Let be a positive solution to system (5). From (31), we obtain

Assume that be a positive solution to the system (5). From (32), we get

To prove the global stability in Theorem 2, we use the following notations to simplify the computations: Theorem 4. Let Λ be the positive equilibrium point of system (5) and assume that the conditions in Theorem 2 and Theorem 3/(i) hold. If andthen the positive equilibrium point is globally asymptotically stable and where

denotes the positive solution of system (5).

Proof. Let V be an appropriate Lyapunov function defined bywhere and .

The change along the solutions of the system isBy considering the first equation of system (5), we haveFrom (33), we obtainMoreover,ifThen, this implies that ΔV 1(n) < 0 and, thus we have Similarly, we can obtain the conditions andThus, , , and , which completes the proof. □

IV. Spread of Coronavirus with Control Parameters for an Infection in Early Detection

In [24], Verhulst considered the logistic growth function to explain mono-species growth. If x represents the population size at time t, then the logistic growth equation has the formwhere r and K are positive numbers. However, in biological phenomena, many situations require modifications in the main model to explain the growth of the population in low density-size, which is well known as the Allee effect [25]. The Allee effect can be divided into two main types: (i) strong Allee effect and (ii) weak Allee effect. A population with a strong Allee effect will have a critical population size, which is the threshold of the population, and any size that is less than the threshold will go to extinction without any further aid. On the other hand, a population with a weak Allee effect will reduce the per capita growth rate at lower population density or size [26], [27], [28].

Let us incorporate an Allee function to the C 1(t) class at discrete time t such aswhereare functions of Holling type II.

Integrating both sides of (45) on [n, t) and taking we get a difference equation system such asLetwhere we obtain , if the conditions of Theorem 3/(i) hold andThus, screening in discrete time is an essential control parameter to support the Allee function in stabilizing the effect of the spread.

Remark 2 The susceptible class and the classes who do not know they are infected are the main populations that affect the Allee function in stabilizing the spread of transmission. While it is essential to keep human non-infected, the other essential aim is to detect the infected class before the symptoms occur.

For a strong Allee effect, let us assume that the Allee function is given bywhere K 0 represents the Allee threshold of the infected class, which does not know they are infected.

Theorem 5. If, then Proof. On the contrary, assume that From , any ϖ > 0 withthere exists t ϖ > 0 such thatFor t > t ϖ and we have as t → ∞, where we obtain a contradiction. □

To avoid redundancy with Section 3, we stated the following theorems without proof.

Theorem 6. Let Λ be the positive equilibrium point of system (46). Assume that the basic reproduction numbers are R 01 < 1 and R 02 < 1. Thus, if conditions for Theorem 2 holds forthen the equilibrium point Λ of system (46) is locally asymptotically stable. □

In applying a weak Allee effect on system (46), we assume that the Alle function is given bywhere E 1 is the Allee coefficient of the population class, that does not know they are infected.

Theorem 7. Let Λ be the positive equilibrium point of system (46). Assume that the basic reproduction numbers are R 01 < 1 and R 02 < 1. Thus, if conditions for Theorem 2 holds forthen the equilibrium point of system (46) is locally asymptotically stable. □

V. Neimark-Sacker Bifurcation Analysis

In this section, we analyze the conditions for a Neimark-Sacker bifurcation for system (5). The following theorem is essential.

Theorem 8. [29] For a quadratic polynomiala pair of complex conjugate roots of (51) lie on the unit circle if and only ifTheorem 9. Assume that , 1 < 1, and. If where , and , and the basic reproductive numbers are 01 < 1 and 02 < 1, then both and classes show Neimark-Sacker bifurcation.

Proof. The class: Because of the characteristic Eq. (15), we haveFrom (i) we havewhere. It is evident that (ii) and (iii) hold, since ℓ1 < 0 andℓ0 > 0. The condition (iv), shows thatwhere 01 < 1, and for 1 < 1. Considering both (53) and (54), we obtainandwhere .

The class: Considering the characteristic Eq. (22), we haveFrom the conditions of Theorem 9, we obtain where and . This completes the proof. □

From Theorem 9, the characteristic equation of the class iswhere . Thus, the complex eigenvalues areAdditionally, the complex eigenvalues of the are where the characteristic equation is given as

Simulation results and conclusion

Numerical simulations

In this sub-section, we present numerical simulations that are consistent with the theoretical results. Table 2 shows a description of the parameters that are given in system (5). We assume the initial conditions of system (5) as and . The main objective here is to demonstrate the changes in the control parameters; and , where, is the screening rate and is the rate of recognition. We emphasize that any increase in the screening rate might stop the pandemic spread. While at the same time, it is also essential to realize that the recognition of this infection depends on the continuation of the updated information regarding the novel coronavirus Covid-19. It is an essential point to realize that civilians are not necessarily knowledgeable about the infections of the coronavirus. Therefore, they should be guided about fundamental ‘health care’ applications as well as the severe phenomena worldwide through the WHO, media, health institutes.

Table. 2

Description of the parameters.

Parameter	Symbol	Rate
The growth rate of S(t)	r1	0.12
The growth rate of C1(t)	r2	0.12
The growth rate of M(t)	r3	0.18
The growth rate of N(t)	r4	0.1
Logistic rate of S(t)	α1	0.05
Logistic rate of C1(t)	α2	0.1
Logistic rate of C2(t)	α3	0.15
Logistic rate of M(t)	α4	0.01
Logistic rate of N(t)	α5	0.01
Rate of the S(t) population per year	p	1.6
Parametric lost from class S(t) to C1(t)	β1, β2	0.00134, 0.00044
Rate of interaction between S(t)−M(t)	σ1, σ2	0.0001
Predation rate	γ	0.0045
Rate of screening	θ	[0.01,0.05]
Recognition of infection	ɛ1, ɛ2	0.3
A conversion factor of N(t)	δ	0.0044
The attack rate of N(t) to M(t)	e	0.15
Rate of average time on infecting M(t)	h	0.15
Potential infectivity of N(t)	ω	ω ∈ (0, 1]

In Fig. 2 , the blue graph denotes the susceptible class and the red graph shows 1 who do not know they are infected. Fig. 2-(a) represents the transmission of the infection that occurs in epidemic form in some areas. However, it spreads intensively to pandemic phenomena worldwide and covers almost the susceptible class. Here we assume that the screening rate in the hospitals (before the symptoms appear) is around%1. Fig. 2-(b) shows the graph when the symptoms appear late so that the awareness of the infection is low. In this case, the endemic spread starts earlier and might be uncontrolled.

Fig. 2

(a) Spread of the C1 class and effect (b) Spread of the C1 class and effect on the on the susceptible S class, where susceptible S class, where ɛ1=.

Fig. 3 (a) shows that increasing the screening rate up to%5 decreases the spread of the infection, and remains the virus in epidemic form. In Fig. 3(b), we fix the screening rate to%5 but consider the awareness of the symptoms as 1= to compare the difference between Fig. 2(b) and Fig. 3(b). It is seen that to be aware of the symptoms in the early stages is an essential parameter that affects the speed of the transmission. In section IV, the rate of screening was discussed intensively, and it was shown that is one of the essential control parameters. and 2 classes, where . the susceptible class, where the susceptible class, where .

Fig. 3

(a) Spread of the C1 class and effect on (b) Spread of the C1 class and effect on.

In Fig. 4 , we considered the human-to-human infection cases. The blue graph denotes the susceptible class , the red graph the 1 class who do not know they are infected and the green graph denotes the 2 that knows they are infected. It is seen that the infected class that do not know they are infected is higher than the class that is determined as which is the infected class who is tested as positive. This means that the spread of transmission from human to human occurs mainly from the 1, which should be controlled with the parameters 1, 2 and .

Fig. 4

The spread of transmission of S, C1.

Fig. 5 shows a diagram of the population classes of (5). It is seen that the natural host and the intermediate host has a stable dynamical system in the habitat. In contrast, they have only a role as hosts in the transmission of the coronavirus. The pandemic spread undergoes when the infection is transmitted from human to human. The intermediate host (animal) shows only an endemic spread, which should be considered as a minor role in this dynamical structure.

Fig. 5

Diagram of all population classes of (5).

Fig. 6 (a)-(d) show the relation of the susceptible class () and the 1() class, who do not know they are infected. We increase the screening rate in each graph to%1,%5,%10, and%20, respectively. It is noticed that the effect of 1() relative to () decreases.

Fig. 6

(a) Dynamical behavior of (S(t), C1(t)), (b) Dynamical behavior of (S(t), C1(t)), where where (c) Dynamical behavior of (S(t), C1(t)), (d) Dynamical behavior of (S(t), C1(t)), where where

Finally, Fig. 7 (a) and 7 (b) shows the rate of recognition of the dynamical behavior related to the susceptible class () and the 1() class who do not know they are infected. We found that the recognition through health organizations and media are highly operative points to stop the pandemic spread and return it to its endemic form. At first, we considered the rate of recognition as and after that for which are shown in Fig. 7(a) and (b), respectively.

Fig. 7

(a) Dynamical behavior of (S(t), C1(t)), (b) Dynamical behavior of (S(t), C1(t)), where where

Conclusion

In this paper, we first review the reasons for the spread of coronaviruses from the natural host to the human host. After that, we established a model of the novel coronavirus, which is known as COVID-19, described by differential equations with piecewise constant arguments. The model is constructed in alignment with important biological and medical reasons. We divided the model into five sub-classes;

the susceptible class ,

the infected class 1, that does not know they are infected since specific symptoms do not appear,

the infected class 2 that knows they are infected because of some symptoms such as respiratory and intestinal infections, including fever, dizziness, and cough, appeared.

the intermediate domestic host , that has a transmission role from the natural host to the human host

the natural host , that are bats of genus Rhinolophus.

We considered in this study the pandemic infection case; animal to human and human to human. Therefore, the first three equations in the model show human to human transmission, while the spillover from the intermediate infected class to the human host denotes a predator-prey mathematical model, and the transmission from the natural host to intermediate host is a host-parasite model of Holling Type II.

The main results are then stated and proved. In Section 3, we analyzed the local and global stability of the co-existing equilibrium point via the Linearized Stability Theorem and a Lyapunov function, respectively. Theorem 2 and Theorem 3 show the stability results when the natural host population is under control, but the screening of C 1 is not high. We deduced that the necessary reproduction numbers R 01 < 1 and R 02 < 1, that shows the transmission potential of the and the classes, respectively. We concluded that among the human hosts, those who do not know they are infected, are the control class in the spread. In contrast, between the animal hosts, the intermediate class plays a dominant role since that class has an essential role in transmitting the disease from animal to human. We noticed that C 1 can decrease if there might be a periodic screening and awareness of information transmitted through media. For the susceptible class it is more important to keep the population rate per year non-infected. The transmission of the virus to the offspring would reach an uncontrollable phenomenon worldwide.

Later in Section IV, we incorporated the Allee function at a discrete-time t. We analyzed both weak and strong Allee effect and obtained that screening for possible inflectional cases in discrete time is an essential control parameter to support both Allee functions in stabilizing the effect of the spread. We emphasized that the susceptible class and the classes who do not know they are infected are the central populations that affect the Allee function in stabilizing the spread of transmission. While it is the priority to keep human non-infected, the other essential aim is to detect the infected class before the symptoms appear.

In Section V, we obtained that the system demonstrates a Neimark-Sacker bifurcation under specific conditions. It is seen that the basic reproduction number R 01, and the natural host has an essential role in the mentioned bifurcation.

In the end, numerical simulations, along with graphical illustrations, are presented to examine the validity of our theoretical findings. We focused on two control parameters, which are θ, the screening rate and , the rate of recognition. We obtained that if the screening percentage stays low, the spread of infection reaches to a pandemic form since the group who do not know they are infected is the significant risk group in transmission. The rate of recognition shows the behavioral act of the civilians considering the daily information from the health organizations. It is seen that any discrete-time of ‘health care’ protections would expand the pandemic spread over time.

The results of this paper studied a biomedical model that describes the character of coronavirus. The analysis of the model, as well as specific qualitative properties, are discussed throughout the paper. Our study is based on mathematical interpretations and consistent with biological and medical assumptions. We believe that our results are essential and of great significance for further investigations.

Declaration of Competing Interest

The authors declare that they have no known funding agency or personal relationships that have appeared to influence the work reported in this paper.