In search of COVID-19 transmission through an infected prey.

This paper considers a nonlinear dynamical model of an ecosystem, which has been established through combining the classical Lotka-Volterra model with the classic SIR model. This nonlinear system consists of a generalist predator that subsists on two prey species in which disease is becoming endemic in one of them. The dynamical analysis methods prove that the system has a chaotic attractor and extreme multistability behavior, where there are infinitely many attractors that coexist under certain conditions. The occurrence of extreme multistability demonstrates the high sensitivity of the system for the initial conditions, which means that tiny changes in the original prey species could enlarge and be widespread, and that could confirm through studying the complexity of the time series of the system's variables. Simulation results of the sample entropy algorithm show that the changes in the system's variables expand over time. It is reasonable now to consider the endemic in the prey species of the system could evolve to be pandemic such as COVID-19. Consequently, our results could provide a foresight about the unpredictability of the COVID-19 outbreak in its original host species as well as after the transmission to other species such as humans.

Introduction

In the second week of December 2019, unknown viral infection was identified in a small local wild animal market in Wuhan city, Hubei province in China [1]. The Chinese Center for Disease Control and Prevention (CCDC) identified this infection as a novel coronavirus infection, which is announced as COVID-19 [2, 3] and spreads with several deadly waves [4, 5].

COVID-19 is the new virus in the coronavirus family and mostly spreads through transmission among human. In this family, it has been reported that most of these types of viruses have an origin in animals. The infection spreading from animal to human can also be observed as a starting points of many outbreaks, from influenza to HIV and from SARS to COVID-19 [6]. The coronavirus causing COVID-19 is highly related to those that caused the SARS (severe acute respiratory syndrome) pandemic in 2003 [7]. Despite speculation, it is early to predict whether COVID-19 will disappear or stay with us permanently like the flu [8].

However, a predator–prey model was published recently by A. Eilersen et. al [9]. In this model, there is 100% fatal disease in both prey and predator species. In fact, the dynamics of predator–prey–pathogen interactions in general have received significant attention in recent decades [10]. Most attention has been given to the interaction between predators and single-host epidemics [11-13]. Therefore, studying a dynamical model, which describes the interplay between predation and multi-host fatal disease, could be a reasonable option to understand the unpredictable behavior of COVID-19.

Numerous predator–prey models in the literature show chaotic and quasi-periodic attractors for given sets of parameters and initial conditions. For discrete-time systems, the 1D logistic map [14], which describes animal reproduction, is one of the most well-known examples of chaos. Fewer instances of chaos were found in ecological models in continuous time. However, over and above investigations of the dynamics of the nonlinear system revealed some unexpected results. It observed that a dynamical system can exhibit ambiguous behavior, which is more complex than existing chaos in the system [15-17]. The coexisting attractors or multistability is the name of this behavior [18, 19]. Multistability has reported in both continuous and discrete systems [20, 21].

In this article, we have considered the 4D dynamical model of an ecosystem, which constructs by combining the classical Lotka-Volterra model with the classic SIR model [9]. This model describes the interplay between a generalist predator subsists on two prey species in which disease is becoming endemic in one of them. Several numerical analysis methods, including phase space, Lyapunov exponents, and bifurcation diagram, are used in the research to investigate the system’s dynamics. This research, through numerical simulations, demonstrates that the predator–prey system can exhibit various types of multistability behaviors, in which the coexistence of infinity many quasi-periodic attractors, and the coexistence of infinity many chaotic attractors can be observed for the same parameters values. Moreover, the complex multistability behavior of coexisting chaotic and quasi-periodic is identified for specific parameters values. The different types of multistability behaviors reflect the high sensitivity of the system. Thus, tiny changes in the original prey species could enlarge and be widespread, and that could confirm through studying the complexity of the time series of the system’s variables. In this regard, a contour plot based on the sample entropy algorithm has been applied to estimate the complexity of each variable. The complexity results revealed that the coexisting attractors of the same set of parameters exhibit different complexity values; hence, the changes in the system’s variables expand over time. Therefore, it is reasonable now to consider the endemic in the prey species of the system could evolve to be pandemic such as COVID-19.

The rest of this paper is organized as follows. In Sect. 2, we briefly describe the dynamical model of the combination of Lotka–Volterra model with the classic SIR model, and we then investigate its stability and parameters-related chaotic regions. Section 3 provides an overview of the dynamics changes in the system when initials-boosted, and this can illustrate the apparent unpredictability of epidemic outbreaks in prey species and humans. In Sect. 4, we conduct numerical experiment based on sample entropy algorithm to illustrate the complexity for the time series of predator and prey species.

The infected prey model

This section studies a 4D dynamical model of an ecosystem [9], which consists of a generalist predator subsists on two prey species in which a disease is becoming endemic in one of them. Generally, a disease can be spillover from an infected prey to human when the infection is widespread in the prey species [6]. If we consider that one of the coronaviruses, such as COVID-19, is widespread in the prey species of the 4D ecosystem. Consequently, the mode of transmission from the infected prey to human might be occurred in two scenarios, the consumption of infected prey as a source of food, and direct contact or contact through an intermediate host.

However, the 4D ecosystem model is established by combining classical Lotka–Volterra model with the classic SIR model, which can be defined as follows:The above dynamical system is a predator–prey model which consists of a generalist predator subsists on two prey species in which a disease is becoming endemic in one of them. In this model, a, b, c, d, R are the positive constant parameters, and , , , are the state variables defined as follows: is the healthy populations of the susceptible prey species; is infected populations of the prey species; is the populations of the immune prey species; and is the population of predators.

Equilibria and stability

By setting the left-hand side of differential equations in the system (1) equal to zero, one can obtain five different equilibria as follows:Linearizing the system (1) at an arbitrary equilibrium , the following Jacobian matrix is obtained:where . To determine the stability of each equilibrium point, one has to compute the corresponding eigenvalues.

Proposition 1

The following statements are true.

is unstable saddle point.

is unstable point if and/or .

is unstable point when .

are unstable points.

Proof

It is easy to obtain the eigenvalues of the system at by solving the following equation:

Therefore, the eigenvalues of the system at are given by Hence, the equilibrium point is unstable saddle.

The eigenvalues of the system at are given by It is well known that an equilibrium point is unstable only when there is at least one eigenvalue with positive real part. Thus, is unstable point if and/or .

The eigenvalues of the system at are given by Since all the parameters of the system are greater then zero, there is only one eigenvalue can have positive real part, which is . As a result, is unstable point only when .

The eigenvalues of the system at are given by Since , then is unstable equilibrium point. Finally, the characteristic equation of the system at has mixed signs. By Routh–Hurwitz criterion, the real parts of the corresponding eigenvalues have also positive and negative signs, which means that is unstable equilibrium point.

Chaotic regions

As can be seen in the previous subsection, the 4D infected prey model (1) has five different equilibria in which there are always at least three unstable equilibria. Consequently, two important matters can be concluded: (1) it is quite possible to have chaos in this system; (2) in case of presence chaotic attractor, this attractor is called self-excited.

Maximum Lyapunov Exponents (MLE) is one of the fundamental algorithms used to identify chaos in nonlinear dynamical systems. Generally, the MLE of a dynamics system is calculated versus one parameter varying. In this paper, we calculate the MLE of the system (1) with two parameters varying at the same time, which can provide a wider vision about the chaotic regions of the system. Therefore, four sets of parameters are chosen for simulating four cases of MLE with step size of 0.01, as shown in Fig.  1. Based on this figure and with the support of time series of each case, we observe that the system shows chaotic behavior only when MLE . Otherwise, the system shows no chaos.

Fig. 1

The chaotic regions of the system for the initial conditions (1, 1, 1, 1) versus two parameters varying in which the color bars refer to the maximum values of the Lyapunov exponent, so the positive values show the chaotic regions and negative values display the periodic regions: a , , ; b , , ; c , , ; d

However, it crucial to clarify that the chaotic regions of the 4D epidemic outbreaks system (1) in Fig.  1 are valid only when the initial conditions are chosen to be (1, 1, 1, 1) .

Extreme multistability behaviors

In the previous section, we have demonstrated that the system (1) shows chaotic attractors as well as non-chaotic attractors. These two different behaviors are appeared with only one specific set of initial conditions. In fact, initiating from a one set of initial values may not give a full picture about the dynamics of a chaotic system. Therefore, this section investigates the dynamics of the system (1) through choosing randomly several sets of initial conditions in which the chaotic and non-chaotic regions of the system are examined with these sets of initial conditions.

When the parameters of the system (1) are selected as , , , , and R to be a control parameter for the range , the bifurcation diagram of the state variable and MLE are depicted in Fig. 2a, b, respectively. In this figure, we have chosen nine different sets of initial conditions to examine the non-chaotic and chaotic regions of the system. As can be seen in Fig. 2, the system exhibits extreme multistability behaviors in which the coexistence of nine different quasi-periodic attractors is primarily occurred within the range . Meanwhile, the coexistence of several chaotic and quasi-periodic attractors is appeared for . To further illustrate the case of coexisting quasi-periodic attractors, we depict the phase portraits and the time series of the system for , as shown in Fig. 3. The coexistence of nine quasi-periodic attractors is quite clear in the time series of the state variable as well as two different projections. However, it is important to mention that choosing many other sets of initial conditions can produce infinitely many other coexisting quasi-periodic attractors.

Fig. 2

Dynamics of the system (1) with respect to R for , , , with nine different initial conditions where , 0.05, 0.1, 0.15, 0.2, 0.25, 0.35, 0.45: a bifurcation diagram; b maximum Lyapunov exponents

Fig. 3

Quasi-periodic attractors of the system (1) for , , , , and with nine different initial conditions , where , 0.05, 0.1, 0.15, 0.2, 0.25, 0.35, 0.45: a and b phase portraits of coexisting infinitely many quasi-periodic attractors; c the time series for the state variable

Additionally, Fig. 4a, b depicts the coexisting bifurcation models of the state variable and MLE, respectively, for the parameters , , , , and . Five different sets of initial conditions were selected in this figure to check the dynamics of the system, in which the corresponding orbits begin with , where equal to 0.1 (blue), 0.2 (green), 0.3 (red), 0.4 (black), and 0.5 (cyan). The evolution of these orbits begin with coexisting quasi-periodic attractors and then spread over a larger region of the equilibrium points to produce coexistence of chaotic attractors, as illustrated in Fig. 5. As can be seen in Fig. 5, five different chaotic attractors are generated by the system via selecting a certain set of parameters and different sets of initial conditions. In fact, except for these coexisting chaotic attractors, there are infinitely many others, which can be generated by selecting appropriate sets of initial conditions.

Fig. 4

Dynamics of the system (1) with respect to R for , , , with five different initial conditions where , 0.2, 0.3, 0.4, 0.5: a bifurcation diagram; b maximum Lyapunov exponents

Fig. 5

Phase portraits and the time series of infinitely many chaotic attractors for , , , , and with five different initial conditions where , 0.2, 0.3, 0.4, 0.5: a plane; b plane; c the time series of the state variable

The most complicated behavior in nonlinear dynamical systems is the transit from chaotic to non-chaotic behaviors (and vice versa) for a certain sets of system parameters. This is due to inefficient techniques for analyzing the dynamical systems with different sets of initial conditions, especially when the systems exhibit infinitely many coexisting attractors even in a small region of these sets. To check whether the system (1) shows this complicated transitions, we depict the coexisting bifurcation models of the state variable and MLE for the parameters , , , , and , as shown in Fig. 6a, b, respectively. The orbits in Fig. 6 colored in blue, green, and purple begin with the initial conditions (2.5, 0.2, 0.2, 0.2) , (2, 0.2, 0.2, 0.2) , and (0.5, 0.2, 0.2, 0.2) , respectively. It can be seen that there are several transitions from chaotic to quasi-periodic attractors for all three orbits. For instance, the orbit of purple color drops from chaos to quasi-periodic when , as illustrated in Fig. 7a–c. Therefore, it can be observed that the behavior of system (1) is unpredictable in which choosing only three sets of initial conditions can generate coexistence between chaotic and non-chaotic attractors.

Fig. 6

Dynamics of the system (1) with respect to R for , , , with the initial conditions where , 2, 0.5: a bifurcation diagram; b maximum Lyapunov exponents

Fig. 7

Phase portraits and the time series of the coexisting two chaotic attractors with one quasi-periodic attractor for the parameters: , , , , with the initial conditions where , 2, 0.5: a plane; b plane; c the time series of the state variable

Different chaotic systems have different degrees of sensitivity. Based on our simulation results, it can be concluded that the system (1) is very sensitive to its parameters and initial conditions, in which a tiny alteration in the initial values can drive the system to exhibit completely different behavior. In other words, each set of system parameters can produce infinity many coexisting attractors by choosing appropriate sets of initial conditions. As a result, the high-sensitivity performance of the system (1) provides a clear explanation about the unpredictability of the epidemic among infected preys. However, the next section would describe the complexity of the attractors that generate by the system (1).

Complexity of coexisting attractors

This section investigates the complexity performance of each variable of system (1) using one of the fundamental algorithms, which is Sample Entropy (SamEn) [22]. Generally, the SamEn measure is used to estimate how much extra information is required to predict the th output of a trajectory using its previous (t) outputs. Obviously, larger positive SamEn values indicate that the trajectories need a more additional information to be predicted. In our simulations, the parameters of SamEn measure are set as the embedding dimension and tolerance parameter Standard Deviation.

As aforementioned in the previous section, the system (1) shows extreme multistability behaviors, which means that tiny changes in the original prey species could enlarge and be widespread, and that could confirm through studying the complexity of the time series of the system’s variables. Therefore, let us display the complexity of the variables of the system (1) in case of coexistence of different chaotic attractors. This can be achieved by contour plots based on SamEn algorithm, in which x-axis is a varied initial condition and y-axis is a varied control parameter, as shown in Fig. 8a–d. The parameters of this case are set as , , , , and ; meanwhile, the initial conditions are selected to be , where . Second, Fig. 8e–h depicts the complexity of the variables of the system (1) in case coexisting chaotic and quasi-periodic attractors. Here, the parameters of the system are chosen to be , , , , and , whereas the initial conditions are selected to be , where . Based on the above two cases, several matters can be observed as follows: In conclusion, the complexity results of the system’s variables confirm the high sensitivity to initial conditions; hence, the changes in the system’s variables expand over time. Therefore, it is reasonable now to consider the endemic in the prey species of the system could evolve to be pandemic such as COVID-19.

Fig. 8

The complexity results of the system (1) where the color bars refer to the values of the SamEn algorithm: a–d for the parameters , , , and with the initial conditions , where ; e–h for the parameters , , , and with the initial conditions , where

For all values of the parameter R, the system with varying a single initial condition exhibits different complexity values, especially for the variables , , and , in which the complexity values of the first case increase with the increasing of the initial condition ; meanwhile, the complexity values of the second case decrease with the increasing of the initial condition ;

The variable of healthy prey is the most complicated variable in the system; meanwhile, it is obvious that the variable of infected prey is more complex than the variable of immune prey .

Conclusions

This paper investigated the 4D dynamical model of an ecosystem, which constructs by combining the classical Lotka–Volterra model with the classic SIR model. This model describes the interplay between a generalist predator subsists on two prey species in which disease is becoming endemic in one of them. Several numerical analysis methods find that this model shows chaotic behavior and extreme multistability behaviors. The occurrence of extreme multistability comes from the high sensitivity of the system. That means a tiny change in the system’s variable could enlarge over time. This might interpret the widespread disease among prey species till becoming endemic. Since disease spillover from animals to humans could occur when the infection is widespread in the prey species, it is sensible to consider that the endemic in the infected prey of the 4D model can jump to humans. However, the complexity of the disease in infected prey species can be evaluated by the Sample entropy algorithm. Contour plot results have demonstrated the unpredictability of disease in the prey species. The results of this paper could provide a demonstration of the randomness of the COVID-19 outbreak among prey species and humans.