key: cord-0699600-ewuz1xjk
authors: Martsenyuk, Vasyl; Augustynek, Krzysztof; Urbas, Andrzej
title: On qualitative analysis of the nonstationary delayed model of coexistence of two-strain virus: Stability, bifurcation, and transition to chaos
date: 2020-10-16
journal: Int J Non Linear Mech
DOI: 10.1016/j.ijnonlinmec.2020.103630
sha: 4a9f28c84cfd45936107ce0d91565d203f511d81
doc_id: 699600
cord_uid: ewuz1xjk

The model of interaction of two strains of the virus is considered in the paper. The model is based on a nonstationary system of differential equations with delays and takes into account populations of susceptible, first-time and re-infected individuals across two strains. For small values of the delays, the conditions of global asymptotic stability are obtained with the help of Lyapunov functionals technique. In some special cases the exponential estimates are constructed. On the basis of numerical modeling complex chaotic solutions of the model are obtained. Their investigation is performed with help of nonlinear characteristics, namely bifurcation maps based on Poincare sections and the maximal Lyapunov exponent were obtained.

Current epidemiological diseases are characterized by significant diversity, fast virus mutations, weakness of natural immunity, the possibility of re-infectiousness. In turn, it demands the development of the model of coexistence of different virus strains. In [1] such a model was developed in the case of two serotypes of dengue virus. Namely, the general flowchart is presented in Figure 1 . Here we consider the following sub-populations: S, sub-population of persons susceptible to both virus strains; I k , k = 1, 2, sub-population of persons which are infected with the kth virus strain; R k , k = 1, 2, sub-population of persons which are recovered after kth virus strain (we suppose that the corresponding immunity has been obtained), but they are susceptible to the jth virus strain, j = 1, 2, j = k; Y k , k = 1, 2, sub-population of persons, which are recovered as a result of jth virus strain, j = k, currently they are infected with the kth virus strain; R, sub-population of persons which are recovered after being infected with both virus strains. Let N be the total population size. Since S + k=1,2 I k + R k + Y k + R = N , we may omit the sub-population R.

Here we give brief explanations for the parameters.

• There is a recruitment rate µN into the susceptible class.

• The natural death rate of all sub-populations is µ.

• The rate of having primary infectiousness with the kth virus strain is β k I k , k = 1, 2.

• The recovered rates for sub-populations I k or Y k are α k , k = 1, 2.

• The rate of having secondary infectiousness with the jth virus strain is σjβjR k Ij, k, j = 1, 2, k = j.

• The parameter σj describes an increase or decrease of susceptibility to strain j as of secondary infection due to possible immune altering as a result of primary infections. In [2] it was investigated the simplification of the model of coexistence of two virus strains. The model is assigned for describing a spread of various virus strains (for example, pandemic and seasonal influenza). In the model there were made the assumptions: 1) compartments of latent persons are not considered; 2) influenza spread is assumed to be necessarily accompanied by the symptoms, i.e., the absence of compartments of the asymptotically infected; 3) a general size of population N is considered constant. Hence, when considering the diagram of transition states we get the following system of ordinary differential equations (ODE) S = µ(N − S) − (β1I1 + β2I2)S,

(1) where k, j = 1, 2, k = j.

In [2] there were investigated the steady states of (1), namely we have three equilibrium states E0 = (N, 0, 0, 0, 0, 0, 0), E1 = (S 1 , I 1 , 0, R 1 , 0, 0, 0), E2 = (S 2 , 0,

µ + α2 are the basic reproductive rates. From the viewpoint of epidemiology the equilibrium states have the following meanings: E0 means the absence of infectiousness; E1 means the presence of the first strain alone; E2 -the second strain alone.

The main stability result for (1) states [2] Theorem 1. E1 is the locally asymptotically stable state of the system (1) if R1 > max(1, σ1 − 1) and the inequality R2 < R1 holds 1 .

The time-delay effects play important roles in the models to capture epidemic dynamics. Here we only mention the most recent study by Liu et al. [3] , which presents the first epidemic dynamic model of COVID-19, and which is based on the distributed time delays. Using the worldwide experimental data on coronavirus infection, they showed the advantages upon conventional SEIR models without delay, when predicting and controlling of COVID-19.

Being inspired by the idea of application of time delays and taking into account the nonstationarity of epidemic dynamics, here we modify the system (1) as,

where k, j = 1, 2, k = j.

Here τ k > 0, k = 1, 2 are the incubation periods, the time during which the infectious agent develops in the host, and only after that time the infected person becomes infectious itself. Therefore, the number of actively infected persons at time t is arising from the contacts of actual population of susceptible and infected at time t − τ k , where τ k are the discrete time delays.

Further in (2) we omit the differential equation for R(t), the number of re-recovered person, since it can be checked that

The biologically significant area for (2) is

We mention that there are the other developments related to the problem of the virus coexistence [4] - [6] . In [7] the model was designed with the help of the system with continuously distributed time delay. Although the problem of coexistence of different strains of the virus has been investigated in a number of studies, the impact of the incubation period, which corresponds to the delay value in the differential equation systems, has not yet been properly investigated. After all, different strains of the virus may be characterized by different values of the incubation period, and the question of how the relationship between these values affects the dynamic behavior of the model is of considerable interest in qualitative research.

We will use the following basic notations S I 1 • R+, the real positive numbers,

• let for any function f :

• 

where τmax = max(τ k , k = 1, 2).

Theorem 2. Consider the system, which is given by the equations (2) and initial conditions (3) . Then there is unique positive solution S(t), I k (t), R k (t), Y k (t), k = 1, 2 satisfying (2) and (3) at t > 0.

Proof. Given anyŜ(t),Î k (t),R k (t),Ŷ k (t) ∈ C + [−τmax, 0], k = 1, 2, since the righthand sides of (2) imply Lipschitz condition, there exists a unique trajectory of (2) starting from (3) [8] .

Henceforth, we will focus our attention on positiveness of the solution of the system (2) .

Since N (t) ≥ S(t), from the first equation of (2) we get

We prove that I k (t) > 0, t > 0. If not, then there exists t k,1 > 0 such as I k (t k,1 ) = 0 and I k (t) ≥ 0, t ∈ [−τmax, t k,1 ]. Furthermore,

Hence,

From the third equation of system (2), we have

From the third equation of system (2), we have

J o u r n a l P r e -p r o o f

When investigating the global asymptotic stability of equilibrium point, they define it as a stable solution which is a global attractor also [9] (see Definition 1.1). The natural desire is to introduce the notion of a globally asymptotically stable system, not focusing on the equilibrium point, which can be done when considering two arbitrary different solutions attracting each other along time. We introduce the following definition of stability. (3) is said to be globally asymptotically stable if, for any two solutions S(t),

The proof of the main stability result will be based on the auxiliary lemma.

Theorem 3. Suppose, there exist the positive values hS, hI k , hR k , hY k , k = 1, 2 such that the functions BS, BI k , BR k , BY k , k = 1, 2 are non-negative and for any sequences ξ l l=1,∞ and η l l=1,

Then the system (2) given the initial conditions (3) is globally asymptotically stable.

Proof. Let S(t), I k (t), R k (t), Y k (t) and s(t), i k (t), r k (t), y k (t)), k = 1, 2 be any two solutions of the system (2) which start from the corresponding initial conditions satisfying (3). Define

The right-upper derivatives of VS(t), VI k (t), VR k (t), VY k (t) along the solutions of system (2) with initial conditions (3) are given by 6 J o u r n a l P r e -p r o o f

Define

Let for values hS, hI k , hR k , hY k , k = 1, 2

J o u r n a l P r e -p r o o f

Using equations (5), the right-upper derivative of V (t) is given by

Applying integration of (6) from 0 to t, we get

When using Lemma 1 to sub-integral function of (7) (which is uniformly continuous), we get (4), which completes the proof.

The following auxiliary results will be used for ordinary differential equations. Lemma 2. [11] Suppose that x ∈ C 1 ([0, ∞), R) satisfies the following inequality

For delayed differential equations the following result will be used. In the work [13] the notion of integral power function of the order a was introduced and investigated.

Definition 2. Integral power function of the order a is called the special function

where for the positive t we mean the integral in the sense of principle value, i.e.

We will use the auxiliary result. 

From the second equation of (2) we get

We assume that

Then, due to Halanay inequality

where λ k > 0 is the root of

From the third equation of (2) we get

Applying Lemma 2 we have

J o u r n a l P r e -p r o o f

With the help of Lemma 4, we get

Let's briefly recall that the key point of the paper [2] and Section 3 is stability research of deterministic behavior of the systems (1) and (2).

In order to demonstrate the application of the results from the previous sections related to global asymptotic stability, we consider the following values for the parameters of (2) Basing on Theorem 2, we conclude that there is a unique positive solution of the system (2) under the parameters' values (12) and (13) . When applying the Theorem 3, we let hS = hI k = hR k = hY k = 1, k = 1, 2. Then the conditions of the Theorem 3 hold and the system (2) is globally asymptotically stable. Moreover, since the parameters of the system satisfy conditions (8), we can construct the exponential estimates (9) and (11) for I k and R k , k = 1, 2. The phase trajectories demonstrating the tending to the equilibrium state are shown in Figure 2 .

The results of numerical integration shows us "deterministic" solutions tending to steady states or limit circles. In practice, experimental research display more complex non-predictable behavior of epidemies. It turns out that the non-stationary model (2) that demonstrates chaotic trajectories.

As was shown in the previous research, the value of delay is the thing that affects significantly the qualitative behavior of the delayed system [7] , [14] , [15] . Here we also For further studies, which are not restricted with the globally asymptotically stability, we let the following values for the parameters of (2) 

We consider the values of parameters µ = 0.001, β1 = 0.000001,β2 = 0.000001, α1 = 0.01,α2 = 0.01 (15) for non-stationary case of (2). On the other hand, if

then the system (2) is stationary one. Such values of the parameters were inspired by the work [2] . Further we will recall the system (2) given by the values (14) and (15) as non-stationary system, whereas the system given by (14) and (16) 

Primary numerical data to make a decision on qualitative behavior of the system may be obtained from the bifurcation plots. Here we construct them on the basis of Poincare sections [16] . Firstly we have constructed phase trajectories and Poincare sections for different values of delays τ1 = τ2 ∈ [0, 30] and visualised them in 3D-space given by the coordinates S, I1, I2.

For any solution (S(t), I1(t), I2(t), R1(t), R2(t), Y1(t), Y2(t)) ∈ Ω of (2), (3) we consider Poincare section with respect to the hyperplane H = (S, I1, I2, R1, R2, Y1, Y2) | I1 = I 1 ,

. At the primary stage of analysis Poincare sections described above were used. In case of stationary system the 3D phase plots are depicted in Figure 3 . We observed at different values of the delay trajectories tending to stable focus, limit focus or "strange" attractor (see Figure 3 ). Corresponding Poincare sections are in Figure 6 .

Nonstationary system (2), (3) demonstrates more complex behavior, starting from limit cycles and tending to unstable manifolds as depicted in Figure 5 Here we have used Poincare sections with the aim of construction of bifurcation plots (see Figure 7 ). As it can be seen from the bifurcation plots, even for small values of time delays nonstationary system tends to periodic solutions whereas stationary one converges to stable focus. 2D bifurcation plots allows us to investigate dynamics only for some special cases of τ1 and τ2. General analysis of common influence of τ1 and τ2 on qualitative behavior requires more adequate visual presentation, which will be offered in the next subsection.

Complete research of qualitative behavior of the model can be executed with the help of investigation of the corresponding time series based on characteristics of nonlinear dynamics. Here we will follow to the general computational approach which was developed and applied earlier in the work [17] .

For the purpose of research of nonlinear dynamics, we have used nonlinearTseries package. The package is a universal software for analyzing signals and time series of any nature. Possible applications of this software include investigation of biological and physiological systems, mechanical oscillations, electrical signals, population dynamics, time series in financial branches, economics, and many others. Here, we also have the kit of standard signal processing functions [18] . The nonlinear statistics and algorithms, including deterministic chaos characteristics, are among them.

Primarily, we can apply Taken's embedding theorem [19] for a dynamical system which states that the estimation of the dimension of the attractor of the initial chaotic dynamical system is based on the use for this purpose of the attractor of this system, obtained by means of a pseudophase reconstruction, i.e., a mapping that associates a point in the time series with a point (here, we denote by square brackets [i] ith element of time series; instant of time t is denoted by round brackets (t))

where i is the discrete time, h is the time delay 2 (in time discrete units), and m is the dimension of the embedding space. Studies by Takens [19] prove that using only one coordinate of a dynamical system, namely, {X[i], X[i + h], . . . , X[i + mh]}, i, m, h ∈ N, we can reconstruct the original attractor in the space of points with delays so that it will retain the most important topological properties and dynamics of the original attractor.

The dimension of the embedding space m is determined by the formula

where d is the fractal dimension of the attractor. Thus, first of all, to perform a full-fledged analysis of chaotic processes, it is necessary to determine the embedding parameters of the dynamic system necessary for the maximum predictability of the chaotic process, namely, the appropriate time delay of the time series and the dimension of the embedding. Various methods are used to select the time delay of the time series. Here we apply the autocorrelation function method, which means that, for each time delay value, the autocorrelation function calculates the coefficient of correlations between the initial time series and its modification, determined using the time delay in steps. The most suitable time delay is selected in accordance with the first zero or close to zero value of the autocorrelation function

For the studied time series, the first approximation of the graph of the autocorrelation function to zero at point h acf = 7; therefore, the most suitable time delay in accordance with this method is 7. The simplicity of the calculation of this method, which does not require large computations, makes it more accessible and is used by many researchers.

Furthermore, using the time-lag, we calculate the embedding dimension with the help of Cao's algorithm [20] . In case of the epidemiologic system (2), in the practical most of cases, Cao's algorithm estimates the embedding dimension as equal to 4. Rarely it was estimated as 5.

In practice, maximal Lyapunov exponents (MLEs) are the most often used [21] . Briefly, if MLE is positive, then the system is chaotic. On the other hand, if it is zero, then we have a periodic or quasi-periodic solution. Finally, for negative MLE, we have a fixed point as a stable attractor. The calculation of Lyapunov exponents was based on the procedure, which can be described in the following manner [22] , [23] . We compute approximations for MLE (describing scaling in phase space) for a sequence of embedding dimensions m = 4, 5, 6, 7. The estimation of MLE uses the regions of the linear dependences of approximations on embedding dimensions. Then, applying linear regression, estimates of MLE can be found. Here, we have used the function maxLyapunov for estimating the divergence rates. MLE is found when applying linear regression to them. In Figure 8 the process of finding such a linear regression is demonstrated.

The values of the MLE were firstly found in the special cases τ1 = τ2 (see Figure 9 ). Here for both stationary and nonstationary cases we observe negative MLEs for small τ . Positive values of MLEs appears in both cases as impulse jumping. For the nonstationary system we also observe continuous interval of positive MLE for approximately τ ∈ [12.5; 22.5]. Lyapunov exponents were calculated on the whole interval τ1, τ2 ∈ [0, 30] (see Figure 10 ). 2D bifurcation maps shows us mutual influence of τ1 and τ2 on the system dynamics. Namely, we see that the delay τ1 has much more impact on qualitative behavior. The reason is that at the values of parameters (14) the trajectory tends to the equilibrium E1.

Thus, we have offered conditions for the global asymptotical stability of (2). As it was shown such stable behavior can be observed only for small values of the time delays, meaning the incubation period of infectiousness.

Further we see a complex nonlinear behavior of the system trajectory. Such nonlinear behavior is caused by a change in a number of model parameters. The Figures  10 show the effect of changes in the incubation periods τ1 and τ2 on the trajectory and solutions of the equation system. It was evidenced with the help of the MLE characterizing the degree of exponential divergence of close trajectories, a positive value of which means that any two close trajectories quickly diverge over time; therefore, the system is sensitive to the values of the initial conditions, which allows us to identify a dynamic system in terms of the presence of chaotic behavior in it.

We see that increasing the incubation period for the two strains of virus under consideration affects the complexity of the trajectories obtained. It should be noted that for small values of the incubation period, the obtained solutions tend to a certain melt value called the endemic solution. At the same time, an increase in the incubation period leads to a periodic solution of the system. Such phenomena in the theory of dynamical systems have been called bifurcation, which occurs when the values of the system parameters change and affect the change in the qualitative behavior of the whole model.

In this case, we transform from a steady endemic focus to a limit cycle. Such a limit cycle corresponds to the situation of periodic epidemics. The impact of other parameters on the change in the qualitative behavior of system trajectories should also be investigated. A further change in the incubation period affects the complexity of such a periodic solution. From a certain value of the incubation periods corresponding to different strains of the virus, there is a doubling of the period, then the period increases 4 times, 8 times and so on. Therefore, the model of coexistence of two strains of viruses was investigated. Such a model can be used to investigate the spread of infectious diseases. Of great importance in the model are the subpopulations of individuals susceptible to the virus, given its two strains.

Note that the nonstationary model shows much more complex behavior as compared with stationary one. Namely, as it can be seen from Figures 5 and 7, we observe the open domains in the plane (τ1, τ2) of chaotic attractors, which must be studied properly.

It is clear that the model can be developed for the cases of three strains, four, etc. In this case, a system of seven ordinary differential equations was proposed as a mathematical object. At the same time, more sophisticated models based on delayed differential equations, stochastic differential equations, partial differential derivative equations can be used in the study of the spatial spread of the epidemic. Of great importance in all these cases is a qualitative study of the nonlinear behavior of the model. We see from numerical studies that at certain values of the parameters of the solution, large values of periods are obtained. Such solutions are called quasiperiodic and correspond to a situation called in the theory of dynamical systems as "deterministic chaos".

The obtained solution trajectories of the proposed model in Figure 5 also indicate the complexity of epidemic prediction. Even in the simplest case of describing a model based on deterministic equations, we get chaotic solutions. This is due to the complexity of nonlinear interaction between subpopulations of the epidemiological model. It is undoubted that further studies should address the use of a seasonal spread of epidemiologically relevant disease that is consistent with the use of non-stationary dynamic models. Also of great importance is the inclusion in the model of populations of symptomatically and asymptotically infected persons.

Coexistence of different serotypes of dengue virus

On conditions of asymptotic stability in SIR-models of mathematical epidemiology

COVID-19: Data-driven dynamics, statistical and distributed delay models, and observations

Coexistence of a cross-diffusive dengue fever model in a heterogeneous environment

Disease persistence and serotype coexistence: An expected feature of human mobility

Coexistence of two dengue virus serotypes and forecasting for madeira island

Stability analysis of delay differential equation models of HIV-1 therapy for fighting a virus with another virus

Introduction to functional differential equations

Local and global stability for lotka-volterra systems with distributed delays and instantaneous negative feedbacks

Systèmes d'équations différentielles d'oscillations non linéaires

V. lakshmikantham, s. leela, a. a. martynyuk: Stability analysis of nonlinear systems. marcel dekker inc

Sharp estimation for the solutions of delay differential and halanay type inequalities

On the problem of chemotherapy scheme search based on control theory

Stability, bifurcation and transition to chaos in a model of immunosensor based on lattice differential equations with delay

A bifurcation path to chaos in a time-delay fisheries predator-prey model with prey consumption by immature and mature predators

Properties of impact events in the model of forced impacting oscillator: Experimental and numerical investigations

Global asymptotic stability and nonlinear analysis of the model of the square immunopixels array based on delay lattice differential equations

Nonlinear time series analysis

Detecting strange attractors in turbulence

Practical method for determining the minimum embedding dimension of a scalar time series

Distributed control and the lyapunov characteristic exponents in the model of infectious diseases

Liapunov exponents from time series

A practical method for calculating largest lyapunov exponents from small data sets