Analysis of stochastic dynamics in a multistable logistic-type epidemiological model.

Motivated by the important problem of analyzing and predicting the spread of epidemics, we propose and study a discrete susceptible-infected model. This logistic-type model accounts such significant parameters as the rate of infection spread due to contacts, mortality caused by disease, and the rate of recovery. We present results of the bifurcation analysis of regular and chaotic survival regimes for interacting susceptible and infected subpopulations. Parametric zones of multistability are found and basins of coexisting attractors are determined. We also discuss the particular role of specific transients. In-phase and anti-phase synchronization in the oscillations of the susceptible and infected parts of the population is studied. An impact of inevitably present random disturbances is studied numerically and by the analytical method of confidence domains. Various mechanisms of noise-induced extinction in this epidemiological model are discussed.

Introduction

The rapidly evolving COVID-19 pandemic has made the traditionally academic epidemiology a priority problem for the global community. The oppressive uncertainty in the prospects makes us look for answers to questions about possible scenarios for the development of epidemics and about the main mechanisms that control the spread of infection. Even a huge amount of information on the development of the epidemic in different countries does not make it possible to understand the internal nonlinear mechanisms of the observed changes. In the study of possible scenarios for the development of epidemics, mathematical modeling plays an important role [1, 2].

Currently, mathematical epidemiology is an actively developing scientific field. Exploring mathematical models, one can study an influence of various biological factors on the epidemiological dynamics. Among these factors, the most significant are the rate of infection spread due to contacts, mortality caused by disease, and the rate of recovery.

In the theoretical study of epidemics, mathematical models with discrete time are effectively used [3-6]. In population dynamics, an important role is played by the simple conceptual models with logistic- and Ricker-type maps [4, 7, 8]. Due to the strong nonlinearity, these models exhibit complex regular and chaotic regimes. Multistability in such systems significantly complicates the analysis and requires taking into account not only the values of the parameters, but also the initial data. The theoretical base in this study is the methods of nonlinear dynamics [9-11] and the mathematical theory of bifurcations [12-14]. To understand the features of the behavior of nonlinear systems, it is necessary to study attractors, their stability, basins of attraction and separatrices between them.

In the study of population systems, it is generally accepted to take into account the impact of unavoidable random perturbations [8, 15–20]. Even small stochastic disturbances in biological systems can dramatically change its dynamics [21-23]. In particular, random disturbances in population systems can lead to the noise-induced extinction [24-27]. Currently, stochastic effects in epidemiological models attract the attention of many researchers [28-33].

In the present paper, we study epidemic spreading processes on the base of a discrete nonlinear two-compartment mathematical model of the interacting susceptible and infected individuals. The processes of infection of susceptible individuals and recovery of the infected ones determine the mutual flows in this system. It is assumed that dynamics of each isolated subpopulation is governed by the logistic map. Logistic-type models are widely used in the study of population dynamics. The logistic map combines simplicity and representativeness, allowing even in the one-dimensional case to simulate equilibrium regimes, as well as regular and chaotic oscillations [8, 34, 35].

In Section 2, we introduce the model and study variability of regimes of deterministic dynamics under variation of key system parameters. We study bifurcation scenarios, attractors, and their basins for different values of the recovery rate, and the rate of spread of infection because of contacts between infected individuals and susceptible ones. Here, along with the conditions of extinction, types of synchronization in the regular and chaotic oscillatory modes of survival are discussed.

In Section 3, we study how stochastic fluctuations in the parameter of the spread of infection affect the behavior of this susceptible-infected model. Various scenarios of noise-induced extinction are determined and investigated by statistics extracted from direct numerical simulation. For the parametric analysis of the noise-induced phenomena, we use semi-analytical method of confidence domains allowing one to take into account mutual arrangement of attractors, their basins, and stochastic sensitivity.

Deterministic model

We consider a system where the total population is divided into two compartments, namely susceptible individuals x and infected individuals y. In mathematical analysis, we use the logistic-type model:Here, is the parameter of growth rate, is the rate of spread of infection because of contacts of infected individuals with susceptible ones, l is the additional death rate of the infected due to the disease, and r is the recovery rate. All biological parameters are non-negative. Since for both subpopulations in (1) the carrying capacity equals one, we consider .

For any parameters, the system (1) has the equilibrium corresponding to the total extinction of both susceptible and infected individuals. This equilibrium is stable for . For , the system possesses the equilibrium which marks extinction of the infected only. Note that in the system (1), extinction of the susceptible with survival of the infected can occur only if and . In this case, the equilibrium appears.

Besides the degenerate equilibria , , and , the system (1) can have a non-degenerate equilibrium with . The coordinate of the equilibrium is the solution of the quadratic equationThe coordinate of the equilibrium is found explicitly:Note that for , the equilibrium has coordinatesThe non-degenerate equilibrium describes the stationary regime of coexistence of susceptible and infected subpopulations.

The study of variability in the behavior of the system under consideration is based on the bifurcation analysis of the model (1), description of its attractors and their basins of attraction.

System dynamics without recovery ()

Let us consider the model (1) for the case when infected individuals do not recover (). In this case, a behavior of the population is determined by parameters , and . Here, we fix and change and .

First, consider corresponding to the equilibrium mode in the classical logistic map. In Fig. 1, it is shown how equilibria of the system (1) change with an increase in the rate of the spread of the infection. Because of additional mortality of the infected individuals due to disease, at the value of is less than . As can be seen, increasing results in a monotonic decrease in the susceptible population until extinction for . The critical value is found from the condition , where is determined in (3). The infected subpopulation changes non-monotonically but stabilizes to when becomes extinct.

Fig. 1

Coordinates of the equilibrium of the system (1) versus for

It should be emphasized that variability of the system (1) dynamics is determined not only by the system parameters but also by the choice of initial states. Indeed, the system (1) is multistable: four equilibria and coexist. In Fig. 2, for different values of , basins of attraction of these equilibria along with sample phase trajectories tending to them are shown by different colors. The arrows show the direction of iterations. Here, the basin of corresponding to the total extinction is shown by blue. The basin of the equilibrium which marks extinction of the infected individuals is shown by green. The basin of the equilibrium which marks extinction of the susceptible individuals is shown by red. The basin of the equilibrium corresponding to coexistence of both subpopulations is shown by white. In this white basin, a sample trajectory tending to the equilibrium is plotted by pink.

Fig. 2

Basins of attraction of equilibria of the system (1) with . Here, basin of is shown by blue, the basin of is plotted by green, and the basin of is shown by red. By white, we show basin of the non-degenerate equilibrium . Broken lines with arrows show iterations of the system (1) solutions starting from the different points belonging to different basins

As can be seen, with increase of , the basin of the non-degenerate equilibrium decreases (see Fig. 2a–d), and disappears (Fig. 2e). Here, the persistence zone of this equilibrium regime is Generally, with increase of , a geometry of the basins gets more complicated: simply connected basins transform into multiply connected. This leads to the fact that trajectories starting from close initial values can converge to different equilibria corresponding to fundamentally different modes of coexistence–extinction. Such sensitivity to the choice of initial data significantly complicates predicting the behavior of the system, even in considered case where attractors of the system are equilibria only. Moreover, these figures clearly demonstrate how such a division of the basins into parts complicates the transients themselves.

Consider now how dynamic regimes of the system (1) change under increasing for other values of the parameter .

Let . In this case, isolated subsystems (for ) exhibit stable two-periodic regimes. Details of transformations of attractors under increasing can be seen in Fig. 3a. With increasing , two-periodic regime transforms into quasi-periodic oscillations with attractors in a form of close invariant curves. In the -interval of close invariant curves, windows of periodic attractors are seen. Further, through Neimark–Sacker bifurcation, the close invariant curves transform into equilibria. Examples of attractors for several values of are presented in Fig. 3b. Here, for we have 2-cycle, for and we show closed invariant curves, for the attractor is 9-cycle, and we have a stable equilibrium for .

Fig. 3

System (1) with : a bifurcation diagrams versus , b attractors: 2-cycle for , closed invariant curves for and , 9-cycle for , and the stable equilibrium for

It can be resumed that the strengthening of the spread of infection deforms initial oscillatory regimes from periodic to quasiperiodic, and further suppresses oscillations to equilibrium mode. With a further increase in the parameter , these equilibria behave similar to case . However, here so the persistence zone is slightly enlarged compared to .

Some regimes of non-degenerate coexistence of the susceptible and infected are illustrated by time series for in Fig. 4. In Fig. 4a for , the synchronization of oscillations of the susceptible and infected is shown. Here, in the transient zone, one can see a stage of in-phase synchronization. Further, this regime transforms into anti-phase synchronization. In Fig. 4b for , there is a complex quasi-periodic regime with essential changes in amplitudes and alternation of intervals with in-phase and anti-phase synchronization. In Fig. 4c for , we show 9-cycle with the phase shift by one. In Fig. 4d, one can see how for oscillations are suppressed to the equilibrium mode.

Fig. 4

Time series for x- and y-coordinates of solutions of the system (1) with for a , b , c , d

In Fig. 5 for , it is shown how basins of coexisting attractors change with increase of the parameter . The basins of the degenerate equilibria , , and the non-degenerate regimes are shown by the same colors as for . Here, peculiarities of the complicating geometry of the basins are clearly seen.

Fig. 5

Basins of attraction and phase trajectories of the system (1) with for a , b , c , d

Consider now a case when the isolated subsystem of susceptible individuals is chaotic, whereas the isolated subsystem of infected individuals operates in a two-periodic mode. For example, such a case will be for the fixed . In Fig. 6, for , we show attractors of the system (1) plotted versus parameter .

Fig. 6

Bifurcation diagram of the system (1) with versus . Here, x- and y-coordinates of coexisting attractors are shown in blue and red

Here, three new key points should be mentioned. First of them is connected with the appearance of -zone where the non-degenerate regime exists in a form of two coexisting attractors shown by red and blue. Second, multiple transitions between order and chaos take place. Third point is in the change of the type of extinction: as the bifurcation parameter passes , the non-degenerate chaotic regime of coexistence transforms into chaotic oscillations in the subpopulation of the susceptible individuals and total extinction of the infected ones. Note that the persistence -zone for the whole population essentially contracts to

Details of mutual arrangement of attractors and their basins for different values of can be seen in Fig. 7. For small , non-degenerate regimes are presented by two coexisting chaotic attractors (black and pink in Fig. 7a). The black attractor corresponds to the regime of anti-phase chaotic synchronization of the susceptible and infected parts of the population, and the pink attractor presents regime of in-phase chaotic synchronization.

Fig. 7

Attractors and their basins for the system (1) with for a , b , c , b , e , f

In Fig. 7b for , it is shown how anti-phase synchronization has a form of 4-cycle (black). Here, the system (1) exhibits coexistence of regular and chaotic synchronization. In Fig. 7c, for , both in-phase and anti-phase regimes are regular: 4-cycle (black) and 12-cycle (pink) coexist.

Figure 7d shows that the system (1) for exhibits a single non-degenerate regime with in-phase synchronization in a form of the 4-cycle. In Fig. 7e, for , this non-degenerate regime is chaotic, and basin (green) of the equilibrium consists of many “islands”. Figure 7f shows that there is no white basin; hence, the system (1) for does not have any non-degenerate dynamical regimes. Here, in the phase plane, green color dominates, and the system mostly tends to .

Such sharp destruction of chaotic non-degenerate dynamical regime is explained by the crisis bifurcation. Indeed, as passes the value from left to right, the chaotic attractor touches the border of its basin, and disappears because of crisis bifurcation [36].

Figure 8 illustrates details of time series for several non-degenerate oscillatory regimes. In Fig. 8a, for , where the system (1) has two coexisting non-degenerate oscillatory attractors, we show chaotic in-phase (left) and regular anti-phase (right) synchronization. In Fig. 8b, for , regular in-phase synchronization in a form of 4-cycle can be seen. Chaotic synchronization is illustrated in Fig. 8c for . Here, zones of temporal destruction of the synchronization can be seen.

Fig. 8

Time series of solutions of the system (1) with : a chaotic in-phase (left) and regular anti-phase (right) synchronization for ; b regular in-phase synchronization in a form of 4-cycle for ; c chaotic synchronization for

To summarize, even for , variation of the parameters and implies complicated transitions between regular and chaotic dynamic regimes with different forms of synchronization and non-trivial dependence on initial states.

System dynamics with recovery ()

Consider now how variation of the parameter r changes scenarios of system (1) dynamics. In this study, we restrict ourselves by the “simplest” variant when subsystems are in the stable equilibrium modes. Different non-degenerate dynamic regimes of the system (1) under the change of parameters r and are shown by bifurcation diagrams in Fig. 9 for the fixed values and . Here, we plot x- and y-coordinates of attractors versus parameter r for three values of .

Fig. 9

Bifurcation diagrams of the system (1) with versus parameter r for a , b , c

In Fig. 9a, for , one can see how an increase of the recovery parameter r transforms the equilibrium regime into the oscillatory mode in a form of two-cycle. Interestingly, the amplitude of recovery-induced x-oscillations is much larger than for y-oscillations.

In Fig. 9b for , there is an r-zone with the additional quasiperiodic regime in a form of 2-torus. More complex behavior with periodic, quasiperiodic and chaotic attractors is shown in Fig. 9c for . Here, a narrow zone of bistability with coexistence of the equilibrium and 2-cycle is demonstrated in the enlarged fragment.

Further, we focus on the most intriguing case . Attractors and basins for the system (1) for are shown in Fig. 10 for several values of the parameter r.

Fig. 10

Attractors and their basins for the system (1) with : a the equilibrium for , b coexistence of the equilibrium (red point) and 2-cycle for , c 2-cycle for , d 2-torus for , e 6-piece chaos for , f 1-piece chaos for . Basin of attraction of is shown by blue, for is shown by green

Remember that for , the equilibrium which characterizes extinction of the susceptible part of the population with survival of the infected part does not exists. Therefore, there is no red domains in Fig. 10.

For small recovery rate r, the non-degenerate regime is the equilibrium (see Fig. 10a for ). Under increase of r, a stable 2-cycle appears additionally, and the system (1) becomes bistable (see Fig. 10b for ). With the further increase of r, the equilibrium loses its stability, and 2-cycle becomes a single attractor. This cycle is shown in Fig. 10c for . Further, due to Neimark–Sacker bifurcation, this 2-cycle transforms into quasi-periodic regime in a form of 2-torus (see Fig. 10d for ). After destruction of this torus, the system (1) undergoes a chain of bifurcations with alternation of regular and chaotic modes. In Fig. 10e, 6-piece chaos is shown for . For , the system exhibits non-degenerate oscillatory regime in a form of one-piece chaos (see Fig. 10f for ).

In the following Section, we will show how noise can cause transformations of the deterministic dynamic regimes.

Stochastic system

To study noise-induced phenomena in the system of susceptible-infected population, we will use the following stochastic model:This model takes into account stochastic fluctuations in the parameter of the spread of infection: . Here, are uncorrelated white Gaussian noises with parameters , and is the noise intensity.

First, consider how random noise deforms non-degenerate regime of interaction of susceptible and infected subpopulations in the equilibrium mode. In the deterministic mode, such a regime is observed, for example, for (see Fig. 9c).

In Fig. 11a, for this set of parameters, we plot time series of the stochastic system (4) solutions starting from the stable equilibrium for two values of noise intensity: (blue), (red). For smaller noise, random trajectories exhibit fluctuations near the deterministic equilibrium. For larger noise, after some transient process, solutions show a sharp change: y-coordinate vanishes, and x-coordinate stabilizes to constant. This means that subpopulation of infected individuals becomes extinct, and the susceptible part, because of absence of infection, transits to the corresponding deterministic equilibrium mode.

Fig. 11

Noise-induced extinction of the infected subpopulation in the stochastic system (4) with : a time series for (blue), (red), b random states versus noise intensity , c phase trajectories and basins, d confidence ellipses ()

Details of such transformation can be seen in Fig. 11b where random states are plotted versus noise intensity . For weak noise, the increase of the noise intensity causes a natural growth of dispersion. As the parameter passes some threshold, a behavior of the stochastic system transforms into the deterministic regime with the equilibrium . To clarify a mechanism of the such transformation, one has to consider the basin (white) of attraction of the non-degenerate equilibrium and the basin (green) of attraction of the degenerate equilibrium in Fig. 11c.

For , random states (blue) of the system (4) totally belong to the basin of . For larger noise with intensity , random states (red) fall into the basin of and further reside on the axis where approach the equilibrium .

For the parametric study of such noise-induced extinction of the infected subpopulation, one can apply the stochastic sensitivity analysis and geometrical method of confidence ellipses (see mathematical details in [37] and some applications in [19, 38, 39]).

The confidence ellipse approximates a dispersion of random states around the stable equilibrium. To construct the confidence ellipse, one has to find the stochastic sensitivity matrix W for the stable equilibrium . This matrix is a unique solution of the following equation:where J is the Jacobi matrix at the equilibrium, and the matrix Q characterizes random disturbances. The confidence ellipse around the equilibrium point an be written aswhere are eigenvalues of W, and are orthonormal eigenvectors. Here, and are coordinates of the ellipse in the basis of and P is the fiducial probability.

In Fig. 11d, along with basins of attraction, we plot confidence ellipses for (blue) and (red). For , the ellipse belongs to the basin of . This fact predicts localization of random states near the equilibrium . For , the expanded ellipse partially occupies points of the basin of . Such geometrical arrangement predicts onset of noise-induced transitions to the basin of and successive extinction of the infected subpopulation y.

Consider now how noise impacts the non-degenerate chaotic regime. Such deterministic regime is observed in the system (1) for (see chaotic attractor in Fig. 10f).

In Fig. 12a, it is shown how dispersion of random states of the stochastic system (4) changes with the increase on . Here, we observe an -zone where noisy chaos transforms into regular regime of stabilization at the equilibrium . Comparing with the case (Fig. 11b), one can see that such noise-induced extinction of the infected occurs at a slightly lower noise intensity . Details of the statistical analysis of this phenomenon are presented in Fig. 12b where plots of probability of noise-induced extinction of the infected subpopulation are shown for different values of r. Here, it is well seen how the -zone for onset of noise-induced extinction depends on the recovery parameter r: the larger r the smaller noise causes the extinction.

Fig. 12

Noise-induced extinction of the infected subpopulation in the stochastic system (4) with : a random states for versus noise intensity , b probability of noise-induced extinction of the infected subpopulation for different r

Consider now a case when both populations become extinct because of random noise. In Fig. 13a, we show an example of transition from the non-degenerate chaotic attractor to the equilibrium . Here, time series of the stochastic system (4) with are plotted for .

Fig. 13

Noise-induced extinction of both subpopulations in the stochastic system (4) with : a time series for , b probability of noise-induced extinction for different r

Results of the statistical analysis of such total extinction are presented in Fig. 13b. Here, probability of the extinction of both susceptible and infected subpopulations is shown for several values of the parameter r. Note that in the considered range of the parameter , this probability does not reach one.

Comparing Figs. 12b and 13b, one can assume that the subpopulation of the susceptible is more resistant to noise than the subpopulation of the infected. For example, for , the infected become extinct with probability close to one, whereas the total extinction has probability close to 0.2. This means that probability of the survival of the susceptible individuals is about 0.8.

To summarize, noise can destroy the persistent coexistence of the susceptible and infected parts of the population and cause extinction of one subpopulation or the whole population system.

Conclusion

This paper is devoted to the actual problem of analyzing possible epidemiological scenarios using mathematical models of the spread of infections. We are studying this problem on the basis of a discrete susceptible-infected model, taking into account such significant parameters as the rate of infection spread due to contacts, mortality caused by disease, and the rate of recovery. As shown by the bifurcation analysis, even in the deterministic case, this simple conceptual logistic-type model demonstrates a wide variety of dynamic modes, in particular, regular and chaotic survival regimes for both subpopulations, as well as various extinction scenarios. It was found that even an insignificant variation in the parameters can radically change the geometry of the basins of attractors corresponding to these regimes. The paper investigates the influence of the random disturbances in the parameter of the infection spread due to contacts between susceptible and infected individuals. The phenomenon of noise-induced transformation of the survival mode into the extinction mode was discovered and investigated. A parametric study of this phenomenon was carried out using statistics extracted from direct numerical modeling, as well as an analytical method of confidence regions, which takes into account the relative position of the basins of coexisting attractors and the stochastic sensitivity of these attractors. This work sheds light on the conditions and internal mechanisms of very complex and difficult to predict modes of behavior of epidemiological systems. When studying such systems, the key questions are the variety of equilibrium and oscillatory modes of behavior, transitions between them, leading to qualitative changes in the processes of infection spread, including extinction mechanisms. These biological questions can be answered using the mathematical theory of bifurcations, by studying attractors and their basins, and by analyzing noise-induced transitions. In the context of the biological problems under consideration, an important result of our study is the demonstration of the constructive abilities of the confidence region method using the stochastic sensitivity technique.