A black swan and canard cascades in an SIR infectious disease model

1.   Introduction

Typical mathematical models for the spread of an infectious disease usually implicitly employ an assumption that is so common that in the majority of cases it is used without remarks regarding its justification. This assumption is that the population size is constant in time. However, while for the humans this assumption is more or less justified, as for the majority of human infections the epidemic processes are considerably faster than the demographic processes, for animal infectious diseases the assumptions may be invalid and, in any case, must be used with a certain degree of care. In a model, formulated in the form of a system of ordinary differential equations, the presence of processes progressing at very different time scales usually means that this particular system is a singularly perturbed system, and, hence, it should be treated like such.

The theory and applications of singularly perturbed systems of differential equations, traditionally connected with the problems of fluid dynamics and nonlinear mechanics, has been developed intensively and its methods are actively applied to a wide range of problems in all areas of natural science. This can be explained by the fact that the singularly perturbed systems appear naturally in modelling various real-life processes that are characterized by slow and fast motions simultaneously present. In the majority of applications, it is necessary to consider the behaviour of the studied system as a whole rather than its separate trajectories. This implies investigating the system dynamics by means of a qualitative analysis rather than straightforward numerical simulations.

Typically, to investigate a singularly perturbed system, a combination of asymptotic and geometrical techniques of analysis should be applied ^[1,2,3]. The essence of this approach consists in separating out the slow motions in the system under investigation. Then the order of the differential system can be decreased. However, when the slow integral manifold is attracting, the resulting reduced system of a lower order still inherits the essential elements of the qualitative behaviour of the original system in the corresponding domain and reflects the behaviour of the original models to a high order of accuracy. A mathematical justification of this method can be given by means of the theory of integral (or invariant) manifolds for singularly perturbed systems ^[2]. The foundations of the integral manifolds theory were laid by N. Bogolyubov ^[4] and significant impact on the development of the method was provided by N. Bogolyubov and Yu Mitropolskii ^[5,6] and J. Hale ^[7,8,9]. As to the singularly perturbed systems, pioneering papers were published by K. Zadiraka ^[10,11]. Asymptotic expansions of slow integral manifolds were used firstly in ^[12,13,14]. The geometric theory of singular perturbations was further developed in ^{[1,15,16,17,18]}. A more detailed account of the history of the issue can be found in ^[2].

Interesting results based on the application of this theory to the study of biological models were obtained in ^[19,20,21].

In this paper, we deal with such specific objects of the geometric theory of singular perturbations as slow invariant manifold with a change of stability, and in particular canards, canard cascades, and black swans.

The term "canard" were originally given by French mathematicians to the intermediate periodic trajectories of the van der Pol equation between the small and the large orbits due to their special shapes ^[22]. Later such trajectories were found and investigated in other types of singularly perturbed systems including those of higher dimensions, where the phase variables are vectors (see, for example, ^{[1,2,18,23,24]} and references therein). From the geometrical theory of singular perturbations viewpoint, a canard may be considered as a result of gluing stable and unstable slow invariant manifolds at one point of the breakdown surface due to the availability of an additional scalar parameter in the differential system.

Another example of invariant surfaces of variable stability is a canard cascade ^[3,25], which is a result of gluing stable (attractive) and unstable (repulsive) slow invariant manifolds at several points of the breakdown surface. To allow the possibility of such gluing, an additional function of a slow variable and several additional parameters are required. If we take an additional function of a vector variable parameterizing the breakdown surface, we can glue the stable and unstable slow invariant manifolds at all points of the breakdown surface simultaneously. As a result, we obtain the continuous stable/unstable invariant surface or black swan ^[1,3,18,23]. Such surfaces may be considered as a multidimensional analogue of the notion of a canard.

In this paper, we consider all these objects in a simple Susceptible-Infectious-Removed (or SIR) infectious disease model, where we assume a variable population size. A combination of demographic and epidemic processes in a single model naturally leads to a singularly perturbed model.

2.   Singularly perturbed systems: mathematical background for time scales separation technique

Let us consider an autonomous system of ordinary differential equations

where $X \in R^{m}$, $Y \in R^{n}$, $t \in R$; $\varepsilon$ is a small positive parameter, and $\alpha$ is an additional parameter. Such systems are called singularly perturbed systems, since for $\varepsilon = 0$ the ability to specify an arbitrary initial condition for $Y(t)$ is lost. The usual approach to the qualitative study of system (2.1), (2.2) is to consider, firstly, the so-called limiting system

or the quasi-steady-state assumption (QSSA) system

and then to draw conclusions on the qualitative behavior of the full system (2.1), (2.2) for a sufficiently small non-zero $\varepsilon$. Please note that these two systems, namely, the limiting system and the QSSA system, coincide if functions $f$ and $g$ are independent from $\varepsilon$; that is, if $f = f(X, Y, \alpha)$ and $g = g(X, Y, \alpha)$ hold. For the majority of applied real-life problems, where there is no singularities such as bifurcations, the use of the limiting system or QSSA system instead of the full system gives acceptable results. However, in some cases the approximations (2.3), (2.4) or (2.5), (2.6) are too rough and lack some important properties of the original system. This typically occurs in the systems with bifurcations. For such models, a method, which relies on the theory of integral manifolds, can be used to obtain a more accurate approximation. (An example of this approach can be seen in ^[26] where it was applied to a model of marine bacteriophage dynamics.) The objective is to essentially replace the original system by a system of a lower order on an integral manifold whose dimension is equal to that of the slow subsystem, which, at the same time, preserves the essential properties of the original system (that is, it is asymptotically equivalent to the original system ^[26]).

The degenerate equation (2.4) describes a slow surface $S$ (or slow curve, in the case of a scalar variable $X$). The subset $\Gamma$ of $S$ given by

is referred to as a breakdown surface (breakdown curve, breakdown or jump points in the scalar case) ^[27]. The surface $\Gamma$ divides the slow surface $S$ into foliations where

By the implicit function theorem, a leaf of the slow surface can be represented by a well-defined vector-function $Y = \varphi(X, \alpha)$. A slow surface can comprise several foliations determined by different functions $Y = \varphi_i(X, \alpha)$, domains of which can intersect, depending on a structure of the slow surface. It should be noted that the conditions of the implicit function theorem do not hold on the leaf boundary defined by (2.7). The leaf $S^s$ of the slow surface is stable (or attractive), if the spectrum of the Jacobian matrix

is located in the left open complex half-plane. If there is at least one eigenvalue of the Jacobian matrix with a positive real part, then the corresponding leaf of $S$, say $S^u$, is unstable (or repulsive). The case of complex eigenvalues is not considered in this paper.

In an $\varepsilon$–neighborhood of $S^s$ ($S^u$) there exists a stable (unstable) slow invariant manifold. The slow invariant manifold is defined as an invariant surface of slow motions.

Definition 1. A smooth surface $S_{\varepsilon}$ in $R^m \times R^n$ is called an invariant manifold of the system (2.1), (2.2), if any trajectory of the system that has at least one point in common with $S_{\varepsilon}$ lies entirely on $S_{\varepsilon}$; i.e., if $\big(X(t_0), Y(t_0)\big) \in S$ for any $t_0$, then the integral curve $\big(X(t, \varepsilon), Y(t, \varepsilon)\big) \in S$ for all $t$.

The only invariant manifolds of (2.1), (2.2) discussed here are these of dimension $m$ (the dimension of the slow variable $X$) that can be represented as graphs of vector-valued functions $Y = h(X, \alpha, )_\varepsilon$, where $h(X, \alpha, \varepsilon)$ is assumed to be a sufficiently smooth function of $\varepsilon$. Such invariant manifolds are called manifolds of slow motions; the origin of this term lies in nonlinear mechanics.

If a slow invariant manifold of the system is attractive, then any trajectory of (2.1), (2.2) can be represented as a trajectory on the attractive slow invariant manifold plus asymptotically negligible terms. The flow on the slow invariant manifold is described by equation

The presence of the additional scalar parameter $\alpha$ provides a possibility of gluing the stable and unstable invariant manifolds at the jump point to form a single trajectory that is referred to as the canard. As we mentioned, the term "canard" (or "duck–trajectory") has been originally associated with periodic trajectories only. However, in this paper, by the term "canard" we call a one-dimensional slow invariant manifold of variable stability ^[18].

Definition 2. Trajectories, which initially move along a stable slow invariant manifold and then continue for a while along an unstable slow invariant manifold, are called canards or duck-trajectories.

Remark. Note that a canard is always a one-dimensional slow invariant manifold of variable stability, but the reverse is not always true. For example, in the system

the straight line $y = 0$ represents an invariant manifold of variable stability, but it is not a canard trajectory because it is not a trajectory; it consists of three trajectories: $x < 0, y = 0$; $x = y = 0$, and $x > 0, y = 0$. For more details, the reader can refer to Chapter 8 in ^[2].

Definition 3. Trajectories which initially move along an unstable slow invariant manifold and then continue for a while along a stable slow invariant manifold are called false canard trajectories.

If it is necessary to glue stable and unstable slow invariant manifolds at several breakdown points, then one needs several additional parameters and, as a result, one obtains a slow invariant manifold with multiple change of stability. Such invariant manifolds are referred to as a canard cascade ^[2,24,25].

Definition 4. A continuous one-dimensional slow invariant manifold of system (2.1), (2.2), which contains at least two canards or false canards, is called a canard cascade.

In the case of a planar system, if we take an additional function whose arguments are a vector parameter and a slow variable, we can glue the stable (attractive) and unstable (repulsive) slow invariant manifolds at all breakdown points at the same time. The outcome of this procedure is a canard cascade ^[2,24,25]. If we take an additional function of a vector variable parameterizing the breakdown surface, then we can glue the stable and unstable slow invariant manifolds at all points of the breakdown surface (or curve) at the same time. As a result we obtain a continuous stable/unstable (attractive/repulsive) invariant surface or a black swan ^[1,18,23]. Such surfaces are considered as a multidimensional analogue of a canard.

3.   An SIR model with variable population size

To formulate a Susceptible–Infected–Removed (or SIR) model with variable population size, we assume that a population of size $N(t)$ is divide into the sub-populations of susceptible, uninfected and removed (i.e., recovered and, hence, immune) individuals, of sizes $S(t)$, $I(t)$ and $R(t)$, respectively. That is, equality

holds. We also make the following assumptions:

● in the absence of the infectious disease, the population growth is governed by the logistic growth model, with the carrying capacity $K$ and per capita population growth rate $r$, that is

● all the newborn individuals are susceptible (that is, there is no "vertical transmission");

● the susceptibles are infected by the infectious individuals at rate $\beta S I$;

● immediately after infection, the infected individual moves from the $S$ compartment into the $I$ compartment;

● apart from the natural mortality, the infected individuals are also dying of the disease at per capita rate $\sigma$;

● the infected individuals recover at per capita rate $\delta$;

● the recovered (or removed) individuals have a life-long immunity.

Please note that according to these assumptions, all the newborns go only to the susceptible compartment, whereas the natural (that is, not associated with the disease) deaths occur from all three compartments. Therefore, we have to separate the births and the deaths in the logistic equation. We assume that all the three subpopulations $S$, $I$ and $R$ have the same per capita natural mortality rate and denote the per capita birth and death rates as $\alpha$ and $\mu$, respectively. Then,

where, $\alpha$ and $\mu$ are the birth and death rates, and $g$ and $h$ the weights reflecting the impacts of food shortage on the birth and death rates. These parameters are related to the population growth rate $r$ by the equalities

and

For a specific species, $\alpha$ and $\mu$ can be more or less directly measured (and, in fact, it appears that it is easier to determine these parameters than the growth rate $r$). A specific choice of the weights $h$ and $g$ is not so straightforward and may depend on many factors. For the majority of species, $\alpha, \mu > r$. (In the majority of cases, $\alpha$ and $\mu$ are significantly larger than $r$.)

Taking into consideration the disease-induced mortality, the total population size $N(t)$ satisfies the equation

Hence, under the above-stated assumptions, the dynamics of the three subpopulation of the $SIR$ model with variable population size are described by the following differential equations:

It is easy to see that if the population size $N$ is assumed to be constant, then system (3.5) reduces to the "standard" SIR model

where $B = \alpha N \left(1-h\frac{N}{K}\right)$ and $\bar\mu = \mu \left(1+g\frac{N}{K}\right)$.

The constant population size assumption also allows to omit one of the equations of system (3.6) (usually it is the equation for $R$) reducing the system to two equations (usually, to the first two equations); the third variable then can be found from the equality $N = S+I+R$. For model (3.5) with variable population size such a straightforward reduction is impossible. However, it may be convenient to use Eq (3.4) for $N(t)$ instead of one of the equations of system (3.5). In particular, it can be convenient to replace the equation for the removed subpopulation $R$ by this equation. This yields the following system, which we analyze further in this paper:

4.   Results

Let us consider system (3.7) where, for simplicity, we assume that $g = 0$. (This assumption implies that the shortage of food inhibits the reproduction but not affect the lifespan of the adults; this may correspond to a situation where all, or the most of, mortalities related to the food shortage occur at a very early stage of existence.) Moreover, let us assume that the system phase space is the set

(The later assumption is needed as the positive octant

is not an invariant set of the system, whereas the set $\Omega$ is a positive invariant set.) Let us use re-scaling

and notations

Under these assumptions, taking into account (3.3) and the relation $N = S+I+R$, system (3.7) can be rewritten as

If an average duration of the infectious period is considerably shorter then an average lifespan (that is $\sigma + \delta \gg r$), which condition holds for the majority of infectious diseases, then $\varepsilon$ is a small positive dimensionless parameter. Hence, the geometric theory of singular perturbations can be applied to the analysis of system (4.1).

Please note that in model (4.1) the parameter $R_0$ is the basic reproduction number of the infection; that is, $R_{0}$ is the expected number of secondary infections produced by a single infectious introduced into the entirely susceptible population.

The limiting system for (4.1) is

In ^[28], where a similar model was considered, it is noted that the last equation in (4.2) can be solved for $y$, and then a quasi-steady-state approximation cannot be used for order reduction of system (4.1). However, it is easy to see that this equation has an obvious solution $y\equiv 0$, which is a slow surface of (1). Moreover, $y\equiv 0$ is the exact slow invariant manifold of the system (4.1), such that $y\equiv 0$ is stable (attractive) for $s < 1/R_{0}$ and unstable (repelling) for $s > 1/R_{0}$. That is, $y\equiv 0$ is a black swan.

The flow on the black swan is described by the system

It is a straightforward exercise to check that the exact solution of this system is given by

where $C_1$ and $C_2$ are constants. Constant $C_2$ should be negative in order to ensure that $s < n$ and a corresponding solution is biologically feasible. However, to show a wider picture, at this stage we prefer to not limit ourselves to the biologically feasible solutions only.

System (4.3) has two equilibrium states; namely, a saddle $O(0, 0)$ and a stable node $P_1(1, 1)$. The trajectories of system (4.3) are described by the expression

Figures 1(a, b) show the phase portraits of this system for $\phi > 1$ and $\phi < 1$, respectively. In these Figures, one can see that there are trajectories on the slow invariant manifold that intersect the breakdown line $s = 1/R_{0}$. These trajectories are either canards, or false canards, since they are one-dimensional stable/unstable slow invariant manifolds. Moreover, in Figure 1(b) it can be seen that there are trajectories that intersect the breakdown line twice, and, for $\phi > 1$ and $R_{0} < 1$, even three times. This means that such trajectories are canard cascades. Due to the invariance of the surface $y\equiv 0$, the corresponding canards exist in the full system (4.1); see Figure 2. From the simple analysis above, it also immediately follows that in the full system (4.1) the equilibrium state $(1, 1, 0)$ is asymptotically stable (a node) when $R_0 \leq 1$ (in this case it belongs to the stable part of the black swan $y\equiv 0$) and unstable (a saddle) when $R_0 > 1$ (in this case it belongs to the unstable part of the black swan $y\equiv 0$).

Note that in the full system (4.1), in addition to the equilibria $\bar O (0, 0, 0)$ and $\bar P_1 (1, 1, 0)$, there is one more positive equilibrium $\bar P_2$ that does not lie on the slow invariant manifold $y \equiv 0$ and has the coordinate $s = 1/b$. Depending on the relationship between the parameter values, this singular point can be either stable or unstable. The most interesting case is when this point, being a stable focus (see Figure 3), loses stability when parameters are varied. In this case a supercritical Hopf bifurcation occurs in the system, and a stable limit cycle appears in the system phase space (see Figure 4). This result corresponds to the paradox of enrichment in the predator-prey systems ^[29]. In contrast to simpler 2-dim predator-prey models, in 3-dim and higher dimensional systems the same paradox may occur even for the bilinear incidence rate ^[26].

Trajectories of the system (4.1), starting from initial points in the basin of attraction of the stable part of the plane $P\equiv 0$, are the spirals containing a stable/unstable part of the slow motion along $P\equiv 0$, i.e., they are canards. The similar situation one can observe in the standard 2-dim SIR model.

The comparatively simple model (4.1) considered in this paper exhibits a rather complex non-linear dynamics that includes canard cascades and Hopf bifurcation. The authors have to stress the the majority of the nonlinear phenomena exhibited by this simple model do not exist in the standard 2-dim SIR model. This fact is particularly striking because the assumptions that we applied to this model were specifically suppoused to make this model close to the standard SIR model.

Acknowledgments

Elena Shchepakina and Vladimir Sobolev are supported by the Ministry of Education and Science of the Russian Federation under the Competitiveness Enhancement Program of Samara University (2013–2020).

Conflict of interest

The authors declare there is no conflict of interest.