Demographic population cycles and ℛ₀ in discrete-time epidemic models

ABSTRACT

We use a general autonomous discrete-time infectious disease model to extend the next generation matrix approach for calculating the basic reproduction number, ${\mathcal{R}}_0$, to account for populations with locally asymptotically stable period k cycles in the disease-free systems, where $k\ge 1$. When ${\mathcal{R}}_0<1$ and the demographic equation (in the absence of the disease) has a locally asymptotically stable period k population cycle, we prove the local asymptotic stability of the disease-free period k cycle. That is, the disease goes extinct whenever ${\mathcal{R}}_0<1$. Under the same period k demographic assumption but with ${\mathcal{R}}_0>1$, we prove that the disease-free period k population cycle is unstable and the disease persists. Using the Ricker recruitment function, we apply our results to discrete-time infectious disease models that are formulated for Susceptible-Infectious-Recovered (SIR) infections with and without vaccination, and Infectious Salmon Anemia Virus (ISA$\mathit{v}$) infections in a salmon population. When ${\mathcal{R}}_0>1$, our simulations show that the disease-free period k cycle dynamics drives the SIR disease dynamics, but not the ISAv disease dynamics.

KEYWORDS:

• Next generation matrix method
• Period k cycles
• Ricker recruitment
• ${\mathcal{R}}_0$

1. Introduction

The next generation matrix approach for calculating the basic reproduction number, ${\mathcal{R}}_0$, has been very useful in determining biologically meaningful thresholds for guaranteeing local stability (disease extinction) or instability (disease persistence) of the disease-free equilibrium of infectious disease epidemic models [2,7,8,10,18,25,26]. For autonomous infectious disease models (without seasonality or periodicity), the next generation matrix approach assumes the existence of a locally asymptotically stable disease-free equilibrium in the demographic equation. That is, in the absence of the disease, the population is assumed to be at a static steady state [2,7,10]. However, density dependent birth or recruitment functions, such as the Ricker or Quadratic discrete-time models, are known to generate locally asymptotically periodic (cyclic) or erratic (chaotic) fluctuations in autonomous disease-free demographic equations, where the static equilibrium point is unstable [20]. Such prominent and persistent cyclic fluctuations have been observed in the abundance of some sockeye salmon populations [17], and the Ricker model was first used to study population cycles in fish [22]. Some mammals in high northern latitudes, for example snowshoe hares (Lepus americanus), show regular population cycles. Herbivores with high-amplitude population cycles, such as voles, lemmings, snowshoe hares and forest Lepidoptera, form the heart of terrestrial food web dynamics. Also, red grouse populations in Britain exhibit cyclic fluctuations in abundance, with periodic crashes [14,23].

In this paper, we use a general autonomous discrete-time epidemic model to extend the next generation matrix approach for calculating the basic reproduction number, ${\mathcal{R}}_0$, to account for populations with periodic demographic dynamics. In these models, the basic reproduction number determines local stability $({\mathcal{R}}_0<1,\mathrm{disease}\mathrm{extinction})$ or instability of the extinction periodic cycle $({\mathcal{R}}_0>1,\mathrm{disease}\mathrm{persistence})$, where the disease-free dynamics is locally asymptotically periodic. We apply the extended next generation matrix approach for calculating ${\mathcal{R}}_0$ to two discrete-time epidemic models with Ricker recruitment functions, a Susceptible - Infectious - Recovered (SIR) model with and without vaccination, and an Infectious Salmon Anemia Virus (ISA$\mathit{v}$) model.

In the SIR model, the dynamics of the Ricker demographic equation in the absence of the disease is qualitatively equivalent to the dynamics of the total population, $(S+I+R)$. However, in the ISA$\mathit{v}$ model, the dynamics of the salmon population equation in the absence of the disease is qualitatively different from that of the total salmon population, $(S+I)$. As a result, when ${\mathcal{R}}_0>1$, the period of the disease-free susceptible population cycle in the SIR model is numerically found to determine the period of the infectious population cycle. However, in the ISA$\mathit{v}$ model, the period of the infectious salmon population can be different from that of the disease-free salmon population. In particular, the ISA$\mathit{v}$ disease is capable of forcing population fluctuations in the salmon population (see Figure 5).

Others have computed ${\mathcal{R}}_0$ for time-periodic (nonautonomous) discrete-time and continuous-time disease epidemic models (for example, see [3] and the references therein). The discrete-time autonomous epidemic models for this study are not seasonal or time-periodic. Our model's demographic population cycles are forced by the density-dependent one-hump Ricker recruitment function via period-doubling bifurcations route to chaos [6,11,12,22,26,29].

In Section 2, we derive the next generation matrix method for calculating ${\mathcal{R}}_0$ in autonomous discrete-time epidemic models with demographic population cycles. To demonstrate how the results of Section 2 can be applied to study specific autonomous discrete-time diseases with demographic population cycles, we introduce SIR and ISA$\mathit{v}$ disease models in Sections 3 and 4, respectively. For  each of these models, we compute ${\mathcal{R}}_0$ and use it to establish disease persistence and extinction. In Section 3, we compute a vaccination reproduction number, ${\mathcal{R}}_v,$ and use it to determine how adopting a vaccination control policy in the cyclic environment lowers the values of ${\mathcal{R}}_0$. Our simulations, Figure 2 in Section 3 and Figure 5 in Section 4, illustrate the relationship between ${\mathcal{R}}_0$ and the infectious population cycles. We summarize our results in the Concluding Remarks section.

2. Next generation matrix method and demographic population cycles

A next generation matrix method for computing the basic reproduction number, ${\mathcal{R}}_0,$ for discrete-time compartmental epidemic models with a stable fixed point was developed by Allen and van den Driessche in Theorem 2.1 of [2]; see also Theorem 3.3 in Li and Schneider [18] and Theorem 3 in Cushing and Yicang [8]. As in [2], to describe the method for ${\mathcal{R}}_0$ when the disease-free dynamics is a period k population cycle, for $t=0,1,2,\dots ,$ and $k\ge 1$, we write a general discrete-time compartmental infectious disease model in the form (1) $x(t+1)=G\left(x\left(t\right)\right)$(1) where $\begin{array}{rl}& x\left(t\right)={(x_1\left(t\right),x_2\left(t\right),\dots ,x_n\left(t\right))}^T\in {\mathbb{R}}_{+}^n,\\ & G:{\mathbb{R}}_{+}^n\to {\mathbb{R}}_{+}^n\mathrm{is}\mathrm{in}{C}^1\left({\mathbb{R}}_{+}^n\right),\end{array}$ and ${\mathbb{R}}_{+}^n=\{Y=(y_1,y_2,\dots ,y_n)\mid y_i\ge 0\mathrm{for}\mathrm{each}i\in \{1,2,\dots ,n\left\}\right\}$. Thus, Model (1(1) $x(t+1)=G\left(x\left(t\right)\right)$(1) ) describes the disease dynamics of the states of the population over discrete-time intervals. Suppose that m<n and the states are ordered so that at generation $t\ge 0$, the first m states, denoted as $X_0\left(t\right)=\left(x_1\right(t),x_2(t),\dots ,x_m(t){)}^T$, are the infectious states, and the remaining n−m states are the healthy noninfectious states, denoted as $X_1\left(t\right)=\left(x_{m+1}\right(t),x_{m+2}(t),\dots ,x_n(t){)}^T$. Hence, as in [2], Model (1(1) $x(t+1)=G\left(x\left(t\right)\right)$(1) ) can be written as the following system of difference equations. (2) $\begin{array}{rl}{X}_0(t+1)& =G_0(X_0\left(t\right),X_1\left(t\right))\\ X_1(t+1)& =G_1(X_0\left(t\right),X_1\left(t\right)).\end{array}$(2) Using knowledge of vectors of population sizes in the disease and non-disease compartments, $X_0\left(t\right)$ and $X_1\left(t\right)$, at time t, Model (2(2) $\begin{array}{rl}{X}_0(t+1)& =G_0(X_0\left(t\right),X_1\left(t\right))\\ X_1(t+1)& =G_1(X_0\left(t\right),X_1\left(t\right)).\end{array}$(2) ) predicts the vectors of population sizes in the disease and non-disease compartments, $X_0(t+1)$ and $X_1(t+1)$, at time $(t+1)$, where $t=0,1,2,\dots$. The specific model's unit of time depend on the disease application. For example, the unit of time could be a convenient time for a follow up census. The basic reproduction number, ${\mathcal{R}}_0$, the lifetime production of infections produced per infectious individual, is independent of the model's time scale. We will use Model (2(2) $\begin{array}{rl}{X}_0(t+1)& =G_0(X_0\left(t\right),X_1\left(t\right))\\ X_1(t+1)& =G_1(X_0\left(t\right),X_1\left(t\right)).\end{array}$(2) ) to compute ${\mathcal{R}}_0$ when the demographic population dynamics is oscillatory with period k population cycle.

We assume that the disease-free system has a unique positive disease-free period k population cycle at $z=\{z_{1\mathrm{\infty }},z_{2\mathrm{\infty }},\dots ,z_{k\mathrm{\infty }}\}$, where $k\ge 1$, and for $t=0,1,2,\dots ,$ $X_0\left(t\right)=0$ and $X_1\left(t\right)\ge 0$. In addition, we assume that linearization of Model (1(1) $x(t+1)=G\left(x\left(t\right)\right)$(1) ) about the disease-free period k population cycle yields the linearized system (3) $Y(t+1)=(DG(0,z_{k\mathrm{\infty }})\cdot DG(0,z_{(k-1)\mathrm{\infty }})\cdot \cdots DG(0,z_{1\mathrm{\infty }}))Y\left(t\right),$(3) where for each $j\in \{1,2,\dots ,k\}$, $DG(0,z_{j\mathrm{\infty }})$ is the $n\times n$ Jacobian matrix of Model (1(1) $x(t+1)=G\left(x\left(t\right)\right)$(1) ) evaluated at the disease-free period k cycle. For each $j\in \{1,2,\dots ,k\}$, $DG(0,z_{j\mathrm{\infty }})$ has the following form: (4) $DG(0,z_{j\mathrm{\infty }})=\left(\begin{array}{cc}{F}_j+T_j& O\\ A_j& C_j\end{array}\right),$(4) where the $m\times m$ submatrices $F_j$ and $T_j$ are nonnegative and O is the zero matrix. Each matrix, $F_j$ and $T_j$ is obtained from differentiation with respect to states $X_0$ and evaluated at the disease-free k cycle, $(0,z_{j\mathrm{\infty }})$. As in [2], by letting $G_0\left(X_0\right(t),X_1(t\left)\right)=\mathcal{F}\left(t\right)+$ $\mathcal{T}\left(t\right)$, where $\mathcal{F}$ is the vector of new infections that survive the time interval and $\mathcal{T}$ is the vector of all the other transitions, we respectively obtain $F_j$ and $T_j$. Thus, for each $j\in \{1,2,\dots ,k\}$, $F_j$ is the matrix of new infections and $T_j$ is the transition matrix. Since some of the population may die either from the disease infection or natural causes during the k−cycle, we assume that the spectral radius of the product of the k transition matrices, $\rho (T_k{T}_{k-1}\cdots T_1)$, is less than 1.

The linearized system reduces to the following form: (5) $Y(t+1)=\left(\begin{array}{cc}(F_k+T_k)(F_{k-1}+T_{k-1})\cdots (F_1+T_1)& O\\ \mathbf{A}& C_k{C}_{k-1}\cdots C_1\end{array}\right)Y\left(t\right).$(5) In the disease-free state, we assume that disease-free period k cycle is locally asymptotically stable. That is, the spectral radius of the product matrix, $C_k{C}_{k-1}\cdots C_1$, is less than one. Therefore, we assume (6) $\rho (C_k{C}_{k-1}\cdots C_1),\rho (T_k{T}_{k-1}\cdots T_1)<1.$(6) Notice that the matrix in (5(5) $Y(t+1)=\left(\begin{array}{cc}(F_k+T_k)(F_{k-1}+T_{k-1})\cdots (F_1+T_1)& O\\ \mathbf{A}& C_k{C}_{k-1}\cdots C_1\end{array}\right)Y\left(t\right).$(5) ) is a block matrix with $\rho (C_k{C}_{k-1}\cdots C_1)<1$. Consequently, the stability of the linearized system, (5(5) $Y(t+1)=\left(\begin{array}{cc}(F_k+T_k)(F_{k-1}+T_{k-1})\cdots (F_1+T_1)& O\\ \mathbf{A}& C_k{C}_{k-1}\cdots C_1\end{array}\right)Y\left(t\right).$(5) ), depends on the eigenvalues of the product matrix $(F_k+T_k)(F_{k-1}+T_{k-1})\cdots (F_1+T_1)$, and is independent of the matrix $\mathbf{A}$. To provide a definition for the basic reproduction, ${\mathcal{R}}_0,$ we let $\mathbf{T}\mathbf{=}{T}_k{T}_{k-1}\cdots T_1$ and $\mathbf{F}\mathbf{=}(F_k+T_k)(F_{k-1}+T_{k-1})\cdots (F_1+T_1)-\mathbf{T}.$ Then $\mathbf{F}\mathbf{+}\mathbf{T}\mathbf{=}(F_k+T_k)(F_{k-1}+T_{k-1})\cdots (F_1+T_1).$ The product matrices, $\mathbf{F}$ and $\mathbf{T}$ are entrywise non-negative. That is, $\mathbf{F}\ge 0$ and $\mathbf{T}\ge 0$. We assume that, $\mathbf{F}\mathbf{+}\mathbf{T}$ is irreducible. Furthermore, along the period k population cycle, $\mathbf{F}$ is the effective matrix of new infections and $\mathbf{T}$ is the transition matrix. Since $\mathbf{T}\ge \mathbf{0}$ and $\rho \left(\mathbf{T}\right)<1$, it follows that $(Id-\mathbf{T}{)}^{-1}=Id+{\mathbf{T}\mathbf{+}\mathbf{T}}^2+\cdots +{\mathbf{T}}^n+\cdots ,$ where Id is an $n\times n$ identity matrix. Thus, $(Id-\mathbf{T}{)}^{-1}$ is entrywise non-negative. Let $\psi \left(0\right)$ be the number of initially infected individuals. Then $\mathbf{F}(Id-\mathbf{T}{)}^{-1}\psi (0)$ is an entrywise non-negative vector that gives the expected number of new infections over the period k population cycle. The matrix $\mathbf{F}(Id-\mathbf{T}{)}^{-1}$ has $(i,j)$ entry equal to the expected number of secondary infections in compartment i produced by an infected individual introduced in compartment j. Consequently, when the disease-free system has a unique positive locally asymptotically stable disease-free period k population cycle, the matrix $Q=\mathbf{F}(Id-\mathbf{T}{)}^{-1}$ is the next generation matrix, and the basic reproduction number is ${\mathcal{R}}_0=\rho \left(\mathbf{F}\right(Id-\mathbf{T}{)}^{-1})=\rho \left(Q\right).$ Using the above periodic next generation matrix method notation, we summarize the relationship between the unique positive disease-free period k population cycle at $z_{\mathrm{\infty }}=\{z_{1\mathrm{\infty }},z_{2\mathrm{\infty }},\dots ,z_{k\mathrm{\infty }}\}$, and ${\mathcal{R}}_0$ in the following result.

Theorem 2.1

If $(0,z_{\mathrm{\infty }})$ is a disease-free period k population cycle of the system $x(t+1)=G\left(x\left(t\right)\right),$ then $(0,z_{\mathrm{\infty }})$ is locally asymptotically stable if ${\mathcal{R}}_0=\rho \left(\mathbf{F}\right(Id-\mathbf{T}{)}^{-1})<1,$ but unstable if ${\mathcal{R}}_0>1,$ where using the above notation $\mathbf{F}\ge 0,$ $\mathbf{T}\ge 0,$ $\mathbf{F}\mathbf{+}\mathbf{T}$ is irreducible and $\rho \left(\mathbf{T}\right)<1$.

The proof of Theorem 2.1 is similar to that of Theorem 2.1 in [2] and is omitted. When the disease-free system has a unique locally asymptotically stable disease-free equilibrium (fixed point), then k=1 and $Q=\mathbf{F}(Id-\mathbf{T}{)}^{-1}$ reduces to $Q=F_1(Id-T_1{)}^{-1}$ [2]. In a recent paper, van den Driessche and Yakubu [26], proved that when ${\mathcal{R}}_0<1$ and the demographic dynamics is asymptotically constant or under geometric growth (non-oscillatory), then the disease-free equilibrium point (fixed point) is globally asymptotically stable. Can a disease-free population cycle with period $k\ge 2$ be globally asymptotically stable? This question was answered in 2002 by Elaydi and Yakubu who showed that this is not possible in autonomous discrete-time models such as Model (1(1) $x(t+1)=G\left(x\left(t\right)\right)$(1) ) ([12], Theorem 3).

3. SIR model with demographic population cycles

Susceptible - Infectious - Recovered (SIR) disease epidemic models have been used to study childhood diseases including chicken pox (varicella) and for the host population in the study of vector borne diseases (for example, malaria) [4,5,19,25]. To introduce a SIR disease epidemic model with periodic demographic dynamics, we assume that at each generation $t\in \{0,1,2,\dots \}$, each member of a population is either susceptible ($S_t$), infectious (infected with the disease, $I_t$) or recovered from the disease with life-long immunity ($R_t$), where the total population is $N_t=S_t+I_t+R_t.$ That is, we let $S_t$, $I_t$, $R_t$ and $N_t$, respectively denote the population density of susceptible, infectious, recovered and total population of individuals at time $t\in \{0,1,2,\dots \}$.

Following Yakubu [29] and the references therein, we assume that susceptible individuals who interact with infectious individuals become infectious with probability $\hat{\varphi }\left(I_t\right)=(1-\varphi (I_t\left)\right)$ and remain susceptible with probability $\varphi \left(I_t\right)$ per unit time interval, where the ‘escape’ function $\varphi :{\mathbb{R}}_{+}\to [0,1]$ is a nonlinear decreasing smooth concave up function with $\varphi \left(0\right)=1$. That is, ${\varphi }^{\mathrm{\prime }}\left(I_t\right)<0$ and ${\varphi }^{\mathrm{\prime \prime }}\left(I_t\right)>0$ for all $I_t\ge 0$. For example, when infections are modelled as Poisson processes, then $\varphi \left(I_t\right)=\exp (-\beta I_t)$ with the transmission constant $\beta >0$.

Furthermore, we assume that infectious individuals recover from the infection with constant probability $\gamma \in (0,1)$ and remain infectious with constant probability $\hat{\gamma }=(1-\gamma )$ per generation. We assume that the probability of natural death for each individual in the population is $d\in (0,1)$ per generation, and probability of staying alive per generation is $\hat{d}=(1-d)$. We assume that a susceptible individual has to be in the infectious class before recovering, and recovered individuals have life-long immunity. That is, we consider a disease that is not fatal so we ignore death due to the disease.

To include periodic demography dynamics, we let $g:{\mathbb{R}}_{+}\to {\mathbb{R}}_{+}$ denote the nonlinear differentiable recruitment (birth or immigration) function of individuals to the susceptible class per unit time interval, where $g\left(0\right)=\underset{t\to \mathrm{\infty }}{lim}g\left(N_t\right)=0$. The Ricker recruitment function, $g\left(N_t\right)=r{N}_t{e}^{-b{N}_t},$ is an example of the recruitment function for this study, with the intrinsic growth rate r>0 and the scaling parameter b>0.

Our discrete-time SIR model with mass action type incidence, $\varphi \left(I_t\right)$, implicitly assumes three distinct temporal phases. At the end of each generation, susceptibles become infectious while infectious recover; a fraction of each class is removed; then susceptibles, infectious and recovered reproduce into the susceptible class. These important assumptions distinguish our discrete-time SIR epidemic model from a similar continuous-time differential equation model. Typically, continuous-time differential equation models with similar well-defined distinct temporal phases are non-autonomous. Taking into account the temporal ordering of events, as in [26], we derive our SIR model in the following three steps.

1. Disease Transmission and Recovery: $\begin{array}{rl}{S}_{\left(1\right)}& =S_t\varphi \left(I_t\right)\\ I_{\left(1\right)}& =S_t(1-\varphi \left(I_t\right))+(1-\gamma )I_t\\ R_{\left(1\right)}& =\gamma I_t+R_t\end{array}$ That is, after disease transmission and recovery, $S_{\left(1\right)}$, $I_{\left(1\right)}$ and $R_{\left(1\right)}$ denote densities of susceptibles, infectious and recovered individuals, respectively.

2. Natural Death (Survival): $\begin{array}{rl}{S}_{\left(2\right)}& =(1-d)S_{\left(1\right)}\\ I_{\left(2\right)}& =(1-d)I_{\left(1\right)}\\ R_{\left(2\right)}& =(1-d)R_{\left(1\right)}\end{array}$ That is, after disease transmission, recovery and natural death, $S_{\left(2\right)}$, $I_{\left(2\right)}$ and $R_{\left(2\right)}$ denote densities of susceptibles, infectious and recovered individuals, respectively.

3. Reproduction or recruitment (S, E, I and R into S): $\begin{array}{rl}{S}_{\left(3\right)}& =g\left(N_t\right)+S_{\left(2\right)}\\ I_{\left(3\right)}& =I_{\left(2\right)}\\ R_{\left(3\right)}& =R_{\left(2\right)}\end{array}$ That is, after disease transmission, recovery, natural death and reproduction, $S_{\left(3\right)}$, $I_{\left(3\right)}$ and $R_{\left(3\right)}$ denote densities of susceptibles, infectious and recovered individuals, respectively.

These assumptions and notation lead to the following discrete-time SIR epidemic model. (7) $\begin{array}{rl}{S}_{t+1}& =g\left(N_t\right)+(1-d)S_t\varphi \left(I_t\right)\\ I_{t+1}& =(1-d)S_t(1-\varphi \left(I_t\right))+(1-\gamma )(1-d)I_t\\ R_{t+1}& =\gamma (1-d)I_t+(1-d)R_t,\end{array}$(7) where $t=0,1,2,\dots$. We study Model (7(7) $\begin{array}{rl}{S}_{t+1}& =g\left(N_t\right)+(1-d)S_t\varphi \left(I_t\right)\\ I_{t+1}& =(1-d)S_t(1-\varphi \left(I_t\right))+(1-\gamma )(1-d)I_t\\ R_{t+1}& =\gamma (1-d)I_t+(1-d)R_t,\end{array}$(7) ) with initial conditions $(S_0,I_0,R_0)\in {\mathbb{R}}_{+}^3$. Note that in Model (7(7) $\begin{array}{rl}{S}_{t+1}& =g\left(N_t\right)+(1-d)S_t\varphi \left(I_t\right)\\ I_{t+1}& =(1-d)S_t(1-\varphi \left(I_t\right))+(1-\gamma )(1-d)I_t\\ R_{t+1}& =\gamma (1-d)I_t+(1-d)R_t,\end{array}$(7) ), we assume that events happen in the following order: disease transmission and recovery, survival (death) and reproduction or recruitment. However, in real biological systems, these three events may happen in different orders.

Theorem 3.1

In Model (7(7) $\begin{array}{rl}{S}_{t+1}& =g\left(N_t\right)+(1-d)S_t\varphi \left(I_t\right)\\ I_{t+1}& =(1-d)S_t(1-\varphi \left(I_t\right))+(1-\gamma )(1-d)I_t\\ R_{t+1}& =\gamma (1-d)I_t+(1-d)R_t,\end{array}$(7) ), for each time $t\in \{1,2,\dots \},(S_t,I_t,R_t)\in {\mathbb{R}}_{+}^3$ and there is no population explosion whenever $(S_0,I_0,R_0)\in {\mathbb{R}}_{+}^3$.

Proof.

Since $(S_0,I_0,R_0)\in {\mathbb{R}}_{+}^3$, each equation of Model (7(7) $\begin{array}{rl}{S}_{t+1}& =g\left(N_t\right)+(1-d)S_t\varphi \left(I_t\right)\\ I_{t+1}& =(1-d)S_t(1-\varphi \left(I_t\right))+(1-\gamma )(1-d)I_t\\ R_{t+1}& =\gamma (1-d)I_t+(1-d)R_t,\end{array}$(7) ) is a sum of nonnegative terms. Hence, $(S_t,I_t,R_t)\in {\mathbb{R}}_{+}^3$ for all $t\ge 0$. Next, we prove boundedness of orbits. Summing all the three equations of Model (7(7) $\begin{array}{rl}{S}_{t+1}& =g\left(N_t\right)+(1-d)S_t\varphi \left(I_t\right)\\ I_{t+1}& =(1-d)S_t(1-\varphi \left(I_t\right))+(1-\gamma )(1-d)I_t\\ R_{t+1}& =\gamma (1-d)I_t+(1-d)R_t,\end{array}$(7) ), at each time $t\in \{0,1,2,\dots \}$, the total population is governed by the demographic equation $N_{t+1}=g\left(N_t\right)+(1-d)N_t,$ and $N_t=S_t+I_t+R_t\ge 0$. Since, $g\left(N_t\right)\ge 0$ is smooth and $\underset{t\to \mathrm{\infty }}{lim}g\left(N_t\right)=0$ there exists K>0 and $\epsilon \in (0,1)$ such that $g\left(N_t\right)<\epsilon$ for all $t\ge K$. Thus, $\underset{t\to \mathrm{\infty }}{lim}{N}_t<\mathrm{\infty }$, the set of iterates of all non-negative initial population sizes are bounded and there is no population explosion in Model (7(7) $\begin{array}{rl}{S}_{t+1}& =g\left(N_t\right)+(1-d)S_t\varphi \left(I_t\right)\\ I_{t+1}& =(1-d)S_t(1-\varphi \left(I_t\right))+(1-\gamma )(1-d)I_t\\ R_{t+1}& =\gamma (1-d)I_t+(1-d)R_t,\end{array}$(7) ).

At each time $t\in \{0,1,2,\dots \}$, when $I_t=R_t=0$, Model (7(7) $\begin{array}{rl}{S}_{t+1}& =g\left(N_t\right)+(1-d)S_t\varphi \left(I_t\right)\\ I_{t+1}& =(1-d)S_t(1-\varphi \left(I_t\right))+(1-\gamma )(1-d)I_t\\ R_{t+1}& =\gamma (1-d)I_t+(1-d)R_t,\end{array}$(7) ) reduces to the following disease-free equation. (8) $S_{t+1}=g\left(S_t\right)+(1-d)S_t.$(8) Since $g\left(0\right)=0$, Model (7(7) $\begin{array}{rl}{S}_{t+1}& =g\left(N_t\right)+(1-d)S_t\varphi \left(I_t\right)\\ I_{t+1}& =(1-d)S_t(1-\varphi \left(I_t\right))+(1-\gamma )(1-d)I_t\\ R_{t+1}& =\gamma (1-d)I_t+(1-d)R_t,\end{array}$(7) ) and Equation (8(8) $S_{t+1}=g\left(S_t\right)+(1-d)S_t.$(8) ) respectively have the trivial equilibrium, $P_{-}=(S,I,R)=(0,0,0)$ and $P_{-}=\left(0\right)$. To include periodic demography dynamics in Model (7(7) $\begin{array}{rl}{S}_{t+1}& =g\left(N_t\right)+(1-d)S_t\varphi \left(I_t\right)\\ I_{t+1}& =(1-d)S_t(1-\varphi \left(I_t\right))+(1-\gamma )(1-d)I_t\\ R_{t+1}& =\gamma (1-d)I_t+(1-d)R_t,\end{array}$(7) ), we assume that the demographic Equation (8(8) $S_{t+1}=g\left(S_t\right)+(1-d)S_t.$(8) ) has a unique positive locally asymptotically stable period k population cycle at $z_{\mathrm{\infty }}=\{z_{1\mathrm{\infty }},z_{2\mathrm{\infty }},\dots ,z_{k\mathrm{\infty }}\}$, where $k\ge 1$. That is, in Model (8(8) $S_{t+1}=g\left(S_t\right)+(1-d)S_t.$(8) ) we assume that $z_{\mathrm{\infty }}$ exists and (9) $-1<{\mathrm{\Pi }}_{s=1}^k(g^{\mathrm{\prime }}\left(z_{s\mathrm{\infty }}\right)+(1-d))<1,$(9) where $k\ge 1$.

3.1. ${\mathcal{R}}_0$ for SIR model with demographic population cycles

To compute the basic reproduction number, ${\mathcal{R}}_0$, for Model (7(7) $\begin{array}{rl}{S}_{t+1}& =g\left(N_t\right)+(1-d)S_t\varphi \left(I_t\right)\\ I_{t+1}& =(1-d)S_t(1-\varphi \left(I_t\right))+(1-\gamma )(1-d)I_t\\ R_{t+1}& =\gamma (1-d)I_t+(1-d)R_t,\end{array}$(7) ) by the next generation matrix method, for each $j\in \{1,2,\dots ,k\}$, $F_j=[-\hat{d}{z}_{j\mathrm{\infty }}{\varphi }^{\mathrm{\prime }}\left(0\right)]\mathrm{and}{T}_j=\left[\hat{\gamma }\hat{d}\right],$ where the demographic Equation (8(8) $S_{t+1}=g\left(S_t\right)+(1-d)S_t.$(8) ) is assumed to have a unique positive locally asymptotically stable period k population cycle. Hence, the transition matrix is $\mathbf{T}\mathbf{=}{T}_k{T}_{k-1}\cdots T_1={\left(\hat{\gamma }\hat{d}\right)}^k$ and the matrix of new infections is $\begin{array}{rl}\mathbf{F}& =(F_k+T_k)(F_{k-1}+T_{k-1})\cdots (F_1+T_1)\mathbf{-}\mathbf{T}\\ & ={\left(\hat{d}\right)}^k(\left({\mathrm{\Pi }}_{j=1}^k[-z_{j\mathrm{\infty }}{\varphi }^{\mathrm{\prime }}\left(0\right)+\hat{\gamma }]\right)-{\left(\hat{\gamma }\right)}^k).\end{array}$ For this model, the matrices $\mathbf{T}$ and $\mathbf{F}$ are positive scalars, and $\rho \left(\mathbf{T}\right)<1$. Hence, ${\mathcal{R}}_0=\rho \left(\mathbf{F}{(Id-\mathbf{T})}^{-1}\right)=\frac{{\left(\hat{d}\right)}^k(\left({\mathrm{\Pi }}_{j=1}^k[-z_{j\mathrm{\infty }}{\varphi }^{\mathrm{\prime }}\left(0\right)+\hat{\gamma }]\right)-{\left(\hat{\gamma }\right)}^k)}{1-{\left(\hat{\gamma }\hat{d}\right)}^k}.$ By Theorem 2.1, the disease-free period k population cycle in Model (7(7) $\begin{array}{rl}{S}_{t+1}& =g\left(N_t\right)+(1-d)S_t\varphi \left(I_t\right)\\ I_{t+1}& =(1-d)S_t(1-\varphi \left(I_t\right))+(1-\gamma )(1-d)I_t\\ R_{t+1}& =\gamma (1-d)I_t+(1-d)R_t,\end{array}$(7) ) is locally asymptotically stable when ${\mathcal{R}}_0<1$ and unstable when ${\mathcal{R}}_0>1$. When d=0 (no death) and k=1 in Model (7(7) $\begin{array}{rl}{S}_{t+1}& =g\left(N_t\right)+(1-d)S_t\varphi \left(I_t\right)\\ I_{t+1}& =(1-d)S_t(1-\varphi \left(I_t\right))+(1-\gamma )(1-d)I_t\\ R_{t+1}& =\gamma (1-d)I_t+(1-d)R_t,\end{array}$(7) ), it becomes a simple SIR model with fixed point (non-oscillatory) dynamics, and ${\mathcal{R}}_0$ reduces to ${\mathcal{R}}_0=-z_{\mathrm{\infty }}{\varphi }^{\mathrm{\prime }}\left(0\right)/\gamma$.

Next, we consider an extension of Model (7(7) $\begin{array}{rl}{S}_{t+1}& =g\left(N_t\right)+(1-d)S_t\varphi \left(I_t\right)\\ I_{t+1}& =(1-d)S_t(1-\varphi \left(I_t\right))+(1-\gamma )(1-d)I_t\\ R_{t+1}& =\gamma (1-d)I_t+(1-d)R_t,\end{array}$(7) ). Suppose that a fraction p of individuals are vaccinated at recruitment into the susceptible population. Furthermore, suppose that the vaccine is perfectly effective, so everyone receiving the vaccine is protected from the infection. We assume that $(1-p)g\left(N_t\right)$ enter the S compartment and the remaining $pg\left(N_t\right)$ enter the recovered compartment. Then Model (7(7) $\begin{array}{rl}{S}_{t+1}& =g\left(N_t\right)+(1-d)S_t\varphi \left(I_t\right)\\ I_{t+1}& =(1-d)S_t(1-\varphi \left(I_t\right))+(1-\gamma )(1-d)I_t\\ R_{t+1}& =\gamma (1-d)I_t+(1-d)R_t,\end{array}$(7) ) becomes (10) $\begin{array}{rl}{S}_{t+1}& =(1-p)g\left(N_t\right)+(1-d)S_t\varphi \left(I_t\right)\\ I_{t+1}& =(1-d)S_t(1-\varphi \left(I_t\right))+(1-\gamma )(1-d)I_t\\ R_{t+1}& =pg\left(N_t\right)+\gamma (1-d)I_t+(1-d)R_t.\end{array}$(10) The following system is the disease-free system for Model (10(10) $\begin{array}{rl}{S}_{t+1}& =(1-p)g\left(N_t\right)+(1-d)S_t\varphi \left(I_t\right)\\ I_{t+1}& =(1-d)S_t(1-\varphi \left(I_t\right))+(1-\gamma )(1-d)I_t\\ R_{t+1}& =pg\left(N_t\right)+\gamma (1-d)I_t+(1-d)R_t.\end{array}$(10) ). (11) $\begin{array}{rl}{S}_{t+1}& =(1-p)g(S_t+R_t)+(1-d)S_t\\ R_{t+1}& =pg(S_t+R_t)+(1-d)R_t.\end{array}$(11) Model (11(11) $\begin{array}{rl}{S}_{t+1}& =(1-p)g(S_t+R_t)+(1-d)S_t\\ R_{t+1}& =pg(S_t+R_t)+(1-d)R_t.\end{array}$(11) ) has a trivial DFE, $P_{-}=(S,R)\equiv (0,0)$. Summing the equations for $S_{t+1}$ and $R_{t+1}$ in Model (11(11) $\begin{array}{rl}{S}_{t+1}& =(1-p)g(S_t+R_t)+(1-d)S_t\\ R_{t+1}& =pg(S_t+R_t)+(1-d)R_t.\end{array}$(11) ) results in the following one-dimensional equation in $(S_t+R_t)$. (12) $S_{t+1}+R_{t+1}=g(S_t+R_t)+(1-d)(S_t+R_t).$(12) Equation (12(12) $S_{t+1}+R_{t+1}=g(S_t+R_t)+(1-d)(S_t+R_t).$(12) ) is similar to (8(8) $S_{t+1}=g\left(S_t\right)+(1-d)S_t.$(8) ). By assumption, Equations (8(8) $S_{t+1}=g\left(S_t\right)+(1-d)S_t.$(8) ) and (12(12) $S_{t+1}+R_{t+1}=g(S_t+R_t)+(1-d)(S_t+R_t).$(12) ) have the unique positive locally asymptotically stable period k population cycle, $z_{\mathrm{\infty }}=\{z_{1\mathrm{\infty }},z_{2\mathrm{\infty }},\dots ,z_{k\mathrm{\infty }}\}$. Consequently, the disease-free system with vaccination, Model (11(11) $\begin{array}{rl}{S}_{t+1}& =(1-p)g(S_t+R_t)+(1-d)S_t\\ R_{t+1}& =pg(S_t+R_t)+(1-d)R_t.\end{array}$(11) ), has the unique positive period k population cycle, $z_{p\mathrm{\infty }}=\{((1-p)z_{1\mathrm{\infty }},p{z}_{1\mathrm{\infty }}),((1-p)z_{2\mathrm{\infty }},p{z}_{2\mathrm{\infty }}),\dots ,((1-p)z_{k\mathrm{\infty }},p{z}_{k\mathrm{\infty }})\}.$ In Model (11(11) $\begin{array}{rl}{S}_{t+1}& =(1-p)g(S_t+R_t)+(1-d)S_t\\ R_{t+1}& =pg(S_t+R_t)+(1-d)R_t.\end{array}$(11) ) we assume that $z_{p\mathrm{\infty }}$ is locally asymptotically stable. That is, (13) $\rho (J_k{J}_{k-1}\cdots J_1)<1,$(13) where $J_s=\left[\begin{array}{cc}(1-p)g^{\mathrm{\prime }}\left(z_{s\mathrm{\infty }}\right)+(1-d)& (1-p)g^{\mathrm{\prime }}\left(z_{s\mathrm{\infty }}\right)\\ p{g}^{\mathrm{\prime }}\left(z_{s\mathrm{\infty }}\right)& p{g}^{\mathrm{\prime }}\left(z_{s\mathrm{\infty }}\right)+(1-d)\end{array}\right]$ for each $s\in \{1,2,3,\dots ,k\}$.

It is possible for $z_{p\mathrm{\infty }}$ to be a locally asymptotically stable period k population cycle of Model (11(11) $\begin{array}{rl}{S}_{t+1}& =(1-p)g(S_t+R_t)+(1-d)S_t\\ R_{t+1}& =pg(S_t+R_t)+(1-d)R_t.\end{array}$(11) ), where $z_{\mathrm{\infty }}$ is a locally asymptotically stable period k population cycle of Model (8(8) $S_{t+1}=g\left(S_t\right)+(1-d)S_t.$(8) ) or (12(12) $S_{t+1}+R_{t+1}=g(S_t+R_t)+(1-d)(S_t+R_t).$(12) ). In the following result, we illustrate this possibility when k=1.

Proposition 3.2

Let $z_{\mathrm{\infty }}=\left\{z_{1\mathrm{\infty }}\right\}$ be a positive fixed point of Model (8(8) $S_{t+1}=g\left(S_t\right)+(1-d)S_t.$(8) ) or (12(12) $S_{t+1}+R_{t+1}=g(S_t+R_t)+(1-d)(S_t+R_t).$(12) ), where $-1<g^{\mathrm{\prime }}\left(z_{1\mathrm{\infty }}\right)+(1-d)<-1+d.$ Then $z_{\mathrm{\infty }}=\left\{z_{1\mathrm{\infty }}\right\}$ is a locally asymptotically stable fixed point of Model (8(8) $S_{t+1}=g\left(S_t\right)+(1-d)S_t.$(8) ) or (12(12) $S_{t+1}+R_{t+1}=g(S_t+R_t)+(1-d)(S_t+R_t).$(12) ), and $z_{p\mathrm{\infty }}=\left\{\right((1-p)z_{1\mathrm{\infty }},p{z}_{1\mathrm{\infty }}\left)\right\}$ is a locally asymptotically stable positive fixed point of Model (11(11) $\begin{array}{rl}{S}_{t+1}& =(1-p)g(S_t+R_t)+(1-d)S_t\\ R_{t+1}& =pg(S_t+R_t)+(1-d)R_t.\end{array}$(11) ).

Proof.

Since $-1<g^{\mathrm{\prime }}\left(z_{1\mathrm{\infty }}\right)+(1-d)<-1+d<0$, $z_{\mathrm{\infty }}=\left\{z_{1\mathrm{\infty }}\right\}$ is a locally asymptotically fixed point of Model (8(8) $S_{t+1}=g\left(S_t\right)+(1-d)S_t.$(8) ) or (12(12) $S_{t+1}+R_{t+1}=g(S_t+R_t)+(1-d)(S_t+R_t).$(12) ). Clearly, $z_{p\mathrm{\infty }}=\left\{\right((1-p)z_{1\mathrm{\infty }},p{z}_{1\mathrm{\infty }}\left)\right\}$ is a positive fixed point of Model (11(11) $\begin{array}{rl}{S}_{t+1}& =(1-p)g(S_t+R_t)+(1-d)S_t\\ R_{t+1}& =pg(S_t+R_t)+(1-d)R_t.\end{array}$(11) ). Next, we use a local stability analysis to determine the stability of $z_{p\mathrm{\infty }}$. The Jacobian matrix of Model (11(11) $\begin{array}{rl}{S}_{t+1}& =(1-p)g(S_t+R_t)+(1-d)S_t\\ R_{t+1}& =pg(S_t+R_t)+(1-d)R_t.\end{array}$(11) ) evaluated at $z_{p\mathrm{\infty }}$ is $J_1=\left[\begin{array}{cc}(1-p)g^{\mathrm{\prime }}\left(z_{1\mathrm{\infty }}\right)+(1-d)& (1-p)g^{\mathrm{\prime }}\left(z_{1\mathrm{\infty }}\right)\\ p{g}^{\mathrm{\prime }}\left(z_{1\mathrm{\infty }}\right)& p{g}^{\mathrm{\prime }}\left(z_{1\mathrm{\infty }}\right)+(1-d)\end{array}\right].$ Thus, $\begin{array}{rl}Trace\left(J_1\right)& =g^{\mathrm{\prime }}\left(z_{1\mathrm{\infty }}\right)+2(1-d)\\ \mathit{D}\mathit{e}\mathit{t}\mathit{e}\mathit{r}\mathit{m}\mathit{i}\mathit{n}\mathit{a}\mathit{n}\mathit{t}\left(J_1\right)& =(1-d)(g^{\mathrm{\prime }}\left(z_{1\mathrm{\infty }}\right)+(1-d\left)\right).\end{array}$ By assumption, $-1<g^{\mathrm{\prime }}\left(z_{1\mathrm{\infty }}\right)+(1-d)<-1+d$. Hence, $-d<$ $Trace\left(J_1\right)<0$. That is, $\left|Trace\right(J_1\left)\right|<d$. Also, $-1<g^{\mathrm{\prime }}\left(z_{1\mathrm{\infty }}\right)+(1-d)<-1+d<0$ implies that $d<$ Determinant$\left(J_1\right)+1<2$. Consequently, $\left|Trace\left(J_1\right)\right|<\mathit{D}\mathit{e}\mathit{t}\mathit{e}\mathit{r}\mathit{m}\mathit{i}\mathit{n}\mathit{a}\mathit{n}\mathit{t}\left(J_1\right)+1<2,$ and by the Jury test, $z_{p\mathrm{\infty }}=\left\{\right((1-p)z_{1\mathrm{\infty }},p{z}_{1\mathrm{\infty }}\left)\right\}$ is a locally asymptotically stable positive fixed point of Model (11(11) $\begin{array}{rl}{S}_{t+1}& =(1-p)g(S_t+R_t)+(1-d)S_t\\ R_{t+1}& =pg(S_t+R_t)+(1-d)R_t.\end{array}$(11) ).

The vaccination reproduction number, also called a control reproduction number, denoted by ${\mathcal{R}}_v$, is given by (14) ${\mathcal{R}}_v=\frac{{\left(\hat{d}\right)}^k(\left({\mathrm{\Pi }}_{j=1}^k[-(1-p)z_{j\mathrm{\infty }}{\varphi }^{{}^{\mathrm{\prime }}}\left(0\right)+\hat{\gamma }]\right)-{\left(\hat{\gamma }\right)}^k)}{1-{\left(\hat{\gamma }\hat{d}\right)}^k}$(14) Notice that ${\mathcal{R}}_v\le {\mathcal{R}}_0$, and adopting a vaccination policy decreases the reproduction number. As in Model (7(7) $\begin{array}{rl}{S}_{t+1}& =g\left(N_t\right)+(1-d)S_t\varphi \left(I_t\right)\\ I_{t+1}& =(1-d)S_t(1-\varphi \left(I_t\right))+(1-\gamma )(1-d)I_t\\ R_{t+1}& =\gamma (1-d)I_t+(1-d)R_t,\end{array}$(7) ), the disease-free population cycle of Model (10(10) $\begin{array}{rl}{S}_{t+1}& =(1-p)g\left(N_t\right)+(1-d)S_t\varphi \left(I_t\right)\\ I_{t+1}& =(1-d)S_t(1-\varphi \left(I_t\right))+(1-\gamma )(1-d)I_t\\ R_{t+1}& =pg\left(N_t\right)+\gamma (1-d)I_t+(1-d)R_t.\end{array}$(10) ) with vaccination, $z_{p\mathrm{\infty }}^0=\{((1-p)z_{1\mathrm{\infty }},0,p{z}_{1\mathrm{\infty }}),((1-p)z_{2\mathrm{\infty }},0,p{z}_{2\mathrm{\infty }}),\dots ,((1-p)z_{k\mathrm{\infty }},0,p{z}_{k\mathrm{\infty }})\},$ is locally asymptotically stable and the disease goes extinct when ${\mathcal{R}}_v<1$. However, $z_{p\mathrm{\infty }}^0$ is unstable and the disease persists when ${\mathcal{R}}_v>1$.

When the disease-free system without vaccination, (8(8) $S_{t+1}=g\left(S_t\right)+(1-d)S_t.$(8) ), has the unique positive equilibrium point, $z_{\mathrm{\infty }}=\left\{z_{1\mathrm{\infty }}\right\}$, and $-1<g^{\mathrm{\prime }}\left(z_{1\mathrm{\infty }}\right)+(1-d)<-1+d$, then by Proposition 3.2, $z_{\mathrm{\infty }}=\left\{z_{1\mathrm{\infty }}\right\}$ and $z_{p\mathrm{\infty }}=\left\{\right(S,R)=((1-p)z_{1\mathrm{\infty }},p{z}_{1\mathrm{\infty }}\left)\right\}$ are locally asymptotically stable equilibrium points of Models (8(8) $S_{t+1}=g\left(S_t\right)+(1-d)S_t.$(8) ) and (11(11) $\begin{array}{rl}{S}_{t+1}& =(1-p)g(S_t+R_t)+(1-d)S_t\\ R_{t+1}& =pg(S_t+R_t)+(1-d)R_t.\end{array}$(11) ) respectively. In this case, ${\mathcal{R}}_0=\frac{-\hat{d}{z}_{1\mathrm{\infty }}{\varphi }^{\mathrm{\prime }}\left(0\right)}{1-\left(\hat{\gamma }\hat{d}\right)}\mathrm{and}{\mathcal{R}}_v=(1-p){\mathcal{R}}_0,$ To bring ${\mathcal{R}}_v$ below the threshold value 1, the fraction that needs to be vaccinated to give herd immunity is known to be $p>1-\frac{1}{{\mathcal{R}}_0}$. In the next section, we will use the Ricker recruitment function and demographic population cycles in Model (10(10) $\begin{array}{rl}{S}_{t+1}& =(1-p)g\left(N_t\right)+(1-d)S_t\varphi \left(I_t\right)\\ I_{t+1}& =(1-d)S_t(1-\varphi \left(I_t\right))+(1-\gamma )(1-d)I_t\\ R_{t+1}& =pg\left(N_t\right)+\gamma (1-d)I_t+(1-d)R_t.\end{array}$(10) ) to determine how to bring ${\mathcal{R}}_v$ below the threshold value 1, where ${\mathcal{R}}_0>1$.

3.2. Ricker recruitment function

The Ricker model was first used to model fish populations. However, it is used in theoretical ecology and epidemiology to model all types of population dynamics and recruitment functions [6,11,12,22,26,29]. When the recruitment function is the Ricker Model, $g\left(N_t\right)=r{N}_t{e}^{-b{N}_t}$, the disease-free Equation (8(8) $S_{t+1}=g\left(S_t\right)+(1-d)S_t.$(8) ) without vaccination reduces to (15) $S_{t+1}=(r{e}^{-b{S}_t}+(1-d))S_t.$(15) In this case, $S_{t+1}\le (r+(1-d\left)\right)S_t$. Let ${\mathcal{R}}_d=\frac{r}{d}.$ Then, in the absence of the disease the susceptible population goes extinct when ${\mathcal{R}}_d<1$. However, when ${\mathcal{R}}_d>1$ the disease-free population persists on a fixed point or cyclic or chaotic attractor. When $1<{\mathcal{R}}_d<e^{2/d}$, the disease-free system persists on the fixed point attractor, $z_{\mathrm{\infty }}=\ln {\mathcal{R}}_d/b$. When ${\mathcal{R}}_d>e^{2/d}$, then the disease-free Equation (8(8) $S_{t+1}=g\left(S_t\right)+(1-d)S_t.$(8) ) is known to undergo period-doubling bifurcations route to chaos [6,11,12,26]. Typically, populations that exhibit periodic cycles take on their highest and lowest population sizes at regular time intervals. To illustrate population cycles in (8(8) $S_{t+1}=g\left(S_t\right)+(1-d)S_t.$(8) ), we choose the probability of natural death to be equal to the probability of staying alive per generation so that d=0.5. To scale the population sizes, we choose the scaling parameter, b, to equal 0.1. In Figure 1, we illustrate period-doubling bifurcations in the disease-free Equation (8(8) $S_{t+1}=g\left(S_t\right)+(1-d)S_t.$(8) ) as the intrinsic growth rate, r, is varied between 20 and 140, where d=0.5 and b=0.1. With our choice of parameters, (15(15) $S_{t+1}=(r{e}^{-b{S}_t}+(1-d))S_t.$(15) ) has an asymptotically stable fixed point when $0.5<r<0.5e^4\approx 27.299$, and the fixed point undergoes period-doubling bifurcation at $r=0.5e^4\approx 27.299$. Figure 1 shows windows of population cycles of period $1,2,4,8,\dots$, and the range of the values of the parameter r for the local asymptotic stability of the periodic population cycles.

Figure 1. Period-doubling bifurcations in the disease-free Equation (15(15) $S_{t+1}=(r{e}^{-b{S}_t}+(1-d))S_t.$(15) ), where d=0.5, b=0.1, on the horizontal axis $20\le r\le 140$ and on the vertical axis $20\le S\le 120.$

Figure 2. Bifurcation diagram for SIR model: As ${\mathcal{R}}_0$ increases from values less than 1 to values greater than 1, the model dynamics changes from disease extinction to disease persistence on a locally asymptotically stable period 2 population cycle.

From Figure 1, the disease-free Equation (15(15) $S_{t+1}=(r{e}^{-b{S}_t}+(1-d))S_t.$(15) ) exhibits a locally asymptotically stable period 2 population cycle at (16) $z_{\mathrm{\infty }}=\{z_{1\mathrm{\infty }},z_{2\mathrm{\infty }}\}=\{37.883,65.757\}$(16) when $r=e^4\approx 54.598$, and a locally asymptotically stable period 4 population cycle at (17) $z_{\mathrm{\infty }}=\{z_{1\mathrm{\infty }},z_{2\mathrm{\infty }},z_{3\mathrm{\infty }},z_{4\mathrm{\infty }}\}=\{41.473,85.957,44.560,73.743\}$(17) when $r=e^{4.6}\approx 99.484$.

3.3. SIR model with Ricker recruitment function

Now, we consider Model (7(7) $\begin{array}{rl}{S}_{t+1}& =g\left(N_t\right)+(1-d)S_t\varphi \left(I_t\right)\\ I_{t+1}& =(1-d)S_t(1-\varphi \left(I_t\right))+(1-\gamma )(1-d)I_t\\ R_{t+1}& =\gamma (1-d)I_t+(1-d)R_t,\end{array}$(7) ) with $g\left(N_t\right)=r{N}_t{e}^{-b{N}_t}$, where the model's unit of time is a year. To compute ${\mathcal{R}}_0$ for this example, as in [6] and [26], we assume that the SIR disease infections are modelled as Poisson processes with $\varphi \left(I_t\right)=\exp (-\beta I_t)$ where $\beta >0$. When the disease-free system has a locally asymptotically stable period 2 population cycle, $z_{\mathrm{\infty }}=\{z_{1\mathrm{\infty }},z_{2\mathrm{\infty }}\}$, then (18) ${\mathcal{R}}_0=\frac{{\left(\hat{d}\right)}^2(z_{1\mathrm{\infty }}{z}_{2\mathrm{\infty }}{\beta }^2+\hat{\gamma }(z_{1\mathrm{\infty }}+z_{2\mathrm{\infty }})\beta )}{1-{\left(\hat{\gamma }\hat{d}\right)}^2}.$(18) To illustrate the impact of increasing ${\mathcal{R}}_0$ from values less than 1 to values greater than 1 on Model (7(7) $\begin{array}{rl}{S}_{t+1}& =g\left(N_t\right)+(1-d)S_t\varphi \left(I_t\right)\\ I_{t+1}& =(1-d)S_t(1-\varphi \left(I_t\right))+(1-\gamma )(1-d)I_t\\ R_{t+1}& =\gamma (1-d)I_t+(1-d)R_t,\end{array}$(7) ) with $g\left(N_t\right)=r{N}_t{e}^{-b{N}_t}$, let $b=0.1,\gamma =0.5,d=0.5\mathrm{and}r=e^4$ while β varies between 0 and 0.1. From (16(16) $z_{\mathrm{\infty }}=\{z_{1\mathrm{\infty }},z_{2\mathrm{\infty }}\}=\{37.883,65.757\}$(16) ), with our choice of parameters, the disease-free Equation (8(8) $S_{t+1}=g\left(S_t\right)+(1-d)S_t.$(8) ) has a locally asymptotically stable period 2 population cycle at $z_{\mathrm{\infty }}=\{z_{1\mathrm{\infty }},z_{2\mathrm{\infty }}\}=\{37.883,65.757\}$. Furthermore, as the infection rate β increases from zero to values greater than $\beta \approx 0.03$, ${\mathcal{R}}_0$ increases from values less than 1 (disease extinction, see Figure 2) to values greater than 1 (disease persistence on a period 2 population cycle, see Figure 2). As predicted by Theorem 2.1, Figure 2 shows that the disease-free period 2 population cycle is locally asymptotically stable when ${\mathcal{R}}_0<1$ and unstable when ${\mathcal{R}}_0>1$.

For example, when all parameters are kept fixed at their current values and $\beta =0.05$, then ${\mathcal{R}}_0\approx 2.352>1$ and Figures 2 and 3 show that Model (7(7) $\begin{array}{rl}{S}_{t+1}& =g\left(N_t\right)+(1-d)S_t\varphi \left(I_t\right)\\ I_{t+1}& =(1-d)S_t(1-\varphi \left(I_t\right))+(1-\gamma )(1-d)I_t\\ R_{t+1}& =\gamma (1-d)I_t+(1-d)R_t,\end{array}$(7) ) with $g\left(N_t\right)=r{N}_t{e}^{-b{N}_t}$ has a locally asymptotically stable endemic period 2 population cycle at $\left\{\right(S,I,R)=(23.910,9.930,4.042)\mathrm{and}(S,I,R)=(54.092,7.161,4.504)\}.$

Figure 3. Model (7(7) $\begin{array}{rl}{S}_{t+1}& =g\left(N_t\right)+(1-d)S_t\varphi \left(I_t\right)\\ I_{t+1}& =(1-d)S_t(1-\varphi \left(I_t\right))+(1-\gamma )(1-d)I_t\\ R_{t+1}& =\gamma (1-d)I_t+(1-d)R_t,\end{array}$(7) ) with $g\left(N_t\right)=r{N}_t{e}^{-b{N}_t}$ has a locally asymptotically stable endemic period-2 population cycle when $b=0.1,\gamma =d=0.5,\beta =0.05\mathrm{and}r=e^4$.

3.4. Vaccination in SIR model with ricker recruitment function

In this section, we consider Model (10(10) $\begin{array}{rl}{S}_{t+1}& =(1-p)g\left(N_t\right)+(1-d)S_t\varphi \left(I_t\right)\\ I_{t+1}& =(1-d)S_t(1-\varphi \left(I_t\right))+(1-\gamma )(1-d)I_t\\ R_{t+1}& =pg\left(N_t\right)+\gamma (1-d)I_t+(1-d)R_t.\end{array}$(10) ) with the Ricker recruitment function $g\left(N_t\right)=r{N}_t{e}^{-b{N}_t}$. The disease-free Model (11(11) $\begin{array}{rl}{S}_{t+1}& =(1-p)g(S_t+R_t)+(1-d)S_t\\ R_{t+1}& =pg(S_t+R_t)+(1-d)R_t.\end{array}$(11) ) becomes (19) $\begin{array}{rl}{S}_{t+1}& =(1-p)r(S_t+R_t)e^{-b(S_t+R_t)}+(1-d)S_t\\ R_{t+1}& =pr(S_t+R_t)e^{-b(S_t+R_t)}+(1-d)R_t.\end{array}$(19) The population dynamics of Model (19(19) $\begin{array}{rl}{S}_{t+1}& =(1-p)r(S_t+R_t)e^{-b(S_t+R_t)}+(1-d)S_t\\ R_{t+1}& =pr(S_t+R_t)e^{-b(S_t+R_t)}+(1-d)R_t.\end{array}$(19) ) is similar to that of Equation (15(15) $S_{t+1}=(r{e}^{-b{S}_t}+(1-d))S_t.$(15) ). As in Equation (15(15) $S_{t+1}=(r{e}^{-b{S}_t}+(1-d))S_t.$(15) ), the trivial disease-free equilibrium of Model (19(19) $\begin{array}{rl}{S}_{t+1}& =(1-p)r(S_t+R_t)e^{-b(S_t+R_t)}+(1-d)S_t\\ R_{t+1}& =pr(S_t+R_t)e^{-b(S_t+R_t)}+(1-d)R_t.\end{array}$(19) ), $P_{-}=(S,R)=(0,0)$, is globally asymptotically stable and the population goes extinct when ${\mathcal{R}}_d<1$. However, $P_{-}$ is unstable, and the disease-free population persists when ${\mathcal{R}}_d>1$.  When $1<{\mathcal{R}}_d<e^{2/d}$, as in the one-dimensional disease-free system without vaccination, the disease-free population of Model (19(19) $\begin{array}{rl}{S}_{t+1}& =(1-p)r(S_t+R_t)e^{-b(S_t+R_t)}+(1-d)S_t\\ R_{t+1}& =pr(S_t+R_t)e^{-b(S_t+R_t)}+(1-d)R_t.\end{array}$(19) ) persists on the fixed point attractor, $z_{p\mathrm{\infty }}=\left\{\right(S,R)=((1-p)\left(\ln {\mathcal{R}}_d/b\right),p\left(\ln {\mathcal{R}}_d/b\right)\left)\right\}$. When ${\mathcal{R}}_d>e^{2/d}$, as in Equation (15(15) $S_{t+1}=(r{e}^{-b{S}_t}+(1-d))S_t.$(15) ), the positive fixed point of the 2-dimensional Model (19(19) $\begin{array}{rl}{S}_{t+1}& =(1-p)r(S_t+R_t)e^{-b(S_t+R_t)}+(1-d)S_t\\ R_{t+1}& =pr(S_t+R_t)e^{-b(S_t+R_t)}+(1-d)R_t.\end{array}$(19) ) undergoes period-doubling bifurcations route to chaos [6,11,12,26]. To illustrate this Equation (15(15) $S_{t+1}=(r{e}^{-b{S}_t}+(1-d))S_t.$(15) ) induced period-doubling bifurcations in Model (19(19) $\begin{array}{rl}{S}_{t+1}& =(1-p)r(S_t+R_t)e^{-b(S_t+R_t)}+(1-d)S_t\\ R_{t+1}& =pr(S_t+R_t)e^{-b(S_t+R_t)}+(1-d)R_t.\end{array}$(19) ), set p=0.78 and keep all the other parameters fixed at their current values in Figure 1. As in Figure 1, Figure 4 shows that, under a vaccination programme, both the susceptible and recovered disease-free populations undergo period-doubling bifurcations as r is varied between 20 and 140, where d=0.5 and b=0.1. That is, the population dynamics of the one-dimensional disease-free model without vaccination determines qualitatively the dynamics of the two-dimensional disease-free system with vaccination (bifurcation points of the two models are the same but S values are lower with the vaccination at p=0.78).

Figure 4. Susceptible and recovered disease-free populations of Model (19(19) $\begin{array}{rl}{S}_{t+1}& =(1-p)r(S_t+R_t)e^{-b(S_t+R_t)}+(1-d)S_t\\ R_{t+1}& =pr(S_t+R_t)e^{-b(S_t+R_t)}+(1-d)R_t.\end{array}$(19) ) undergo period-doubling bifurcations as r is varied between 20 and 140, where p=0.78, d=0.5 and b=0.1.

When $b=0.1,\gamma =0.5,d=0.5,\beta =0.05\mathrm{and}r=e^4$ in Model (7(7) $\begin{array}{rl}{S}_{t+1}& =g\left(N_t\right)+(1-d)S_t\varphi \left(I_t\right)\\ I_{t+1}& =(1-d)S_t(1-\varphi \left(I_t\right))+(1-\gamma )(1-d)I_t\\ R_{t+1}& =\gamma (1-d)I_t+(1-d)R_t,\end{array}$(7) ) with $g\left(N_t\right)=r{N}_t{e}^{-b{N}_t}$ and no vaccination $(p=0)$, then ${\mathcal{R}}_0=2.352>1$ where the period of the disease-free population is 2 (see Equation (18(18) ${\mathcal{R}}_0=\frac{{\left(\hat{d}\right)}^2(z_{1\mathrm{\infty }}{z}_{2\mathrm{\infty }}{\beta }^2+\hat{\gamma }(z_{1\mathrm{\infty }}+z_{2\mathrm{\infty }})\beta )}{1-{\left(\hat{\gamma }\hat{d}\right)}^2}.$(18) ) in Section 3.3 for ${\mathcal{R}}_0$ formula). However, keeping all the parameters fixed at their current values, and introducing vaccination at p=0.78 decreases the reproduction number to ${\mathcal{R}}_v=0.232<1$, from Equation (14(14) ${\mathcal{R}}_v=\frac{{\left(\hat{d}\right)}^k(\left({\mathrm{\Pi }}_{j=1}^k[-(1-p)z_{j\mathrm{\infty }}{\varphi }^{{}^{\mathrm{\prime }}}\left(0\right)+\hat{\gamma }]\right)-{\left(\hat{\gamma }\right)}^k)}{1-{\left(\hat{\gamma }\hat{d}\right)}^k}$(14) ) in Section 3.1. That is, as in populations with non periodic demographic dynamics, vaccination programmes can play a role in controlling infectious diseases in populations with periodic demographic dynamics.

4. ISAv model with demographic population cycles

Infectious salmon anemia (ISA), a finfish disease caused by a virus that belongs to a family of viruses called Orthomyxoviridae, has caused significant mortality among salmon farms in Northern Europe, Canada, Maine, and Chile. Worldwide, economic losses due to ISA have been in the billions of dollars. Due to the severity of ISA, the European Union includes infectious salmon anemia in its list of the most dangerous diseases of fish, and it is one of just 10 virus infections of finfish that is reportable to the World Organization for Animal Health [9,15,21,28]. Currently, there are no treatment options available for ISA. Typically, ISAv disease outbreaks have been managed through a combination of regulatory measures and husbandry practices, including restricted movements of fish between farms, enforced slaughtering, use of all-in–all-out programmes at farms, and disinfection of slaughterhouses and processing plants. ISAv vaccine preparations provide partial protection, but their use is limited by difficulties associated with vaccine delivery to large numbers of fish [13,16,21,27].

In [21], Milliken and Pilyugin introduced a two-patch (migration-linked farm and wild salmon populations) ordinary differential equations model of infectious salmon anemia virus with logistic growth. Here, we introduce a single patch discrete-time ISAv model with Ricker stock recruitment, which has been used to describe stock and recruitment in fisheries [22]. As in [21], we assume that at each time $t\in \{0,1,2,\dots \}$, each live salmon is either susceptible, $S_t$ or infectious (infected with ISA disease), $I_t$. That is, we let $S_t$, $I_t$, and $N_t$ = $S_t$ + $I_t$ respectively denote the population size of susceptible, infectious, and total population of live salmon at each time t. Once infected, salmon do not recover from the ISAv disease. Thus, we use an SI epidemic model with no recovery class to describe the salmon population [1]. At each time $t\in \{0,1,2,\dots \}$, we denote the virus population size by $V_t$.

We assume that the susceptible salmon population experiences Ricker growth. As in [21], we assume that infected salmon cannot reproduce. The intrinsic birth rate of salmon per time interval is denoted by r>0. Consequently, assuming all newborn salmon are susceptible, $g\left(S_t\right)=r{S}_t{e}^{-b{S}_t},$ where the scaling parameter b>0. Also, per unit time interval, $d\in (0,1)$  is the fraction of salmon that die ‘naturally’, $\hat{d}=(1-d)$ is the fraction of salmon that survive, $\mu \in (0,1)$ is the constant ISA induced mortality, $\hat{\mu }=(1-\mu )$ is the fraction of ISA infectious salmon that survive, $d_V\in (0,1)$ is the constant fraction of virus that is cleared and ${\hat{d}}_V=(1-d_V)$ is the fraction of virus that survive.

At each time $t\in \{0,1,2,\dots \},$ we assume that a fraction of susceptible salmon, $\theta \in (0,1)$, become infected from direct contact with infectious salmon with probability ${\hat{\varphi }}_I\left(I_t\right)=(1-{\varphi }_I(I_t\left)\right)$, and the remaining susceptible salmon, $\hat{\theta }=(1-\theta )$, become infected via contact with the ISA virus with probability ${\hat{\varphi }}_V\left(V_t\right)=(1-{\varphi }_V(V_t\left)\right)$, where the ‘escape’ functions ${\varphi }_I,{\varphi }_V:{\mathbb{R}}_{+}\to [0,1]$ are nonlinear decreasing smooth concave up functions with ${\varphi }_I\left(0\right)={\varphi }_V\left(0\right)=1$, ${\varphi }_I^{\mathrm{\prime }}<0$, ${\varphi }_V^{\mathrm{\prime }}<0$, ${\varphi }_I^{\mathrm{\prime \prime }}>0$ and ${\varphi }_V^{\mathrm{\prime \prime }}>0$. For example, when ISAv infections are modelled as Poisson processes, then ${\varphi }_I\left(I_t\right)=\exp (-{\beta }_I{I}_t)$ and ${\varphi }_V\left(V_t\right)=\exp (-{\beta }_V{V}_t)$ with the transmission constants ${\beta }_I$, ${\beta }_V>0$. As in [21], we account for virus shedding in our model. During each time interval, $\delta I_t$ is the population of virus shed, where $\delta >0$.

As in Model (7(7) $\begin{array}{rl}{S}_{t+1}& =g\left(N_t\right)+(1-d)S_t\varphi \left(I_t\right)\\ I_{t+1}& =(1-d)S_t(1-\varphi \left(I_t\right))+(1-\gamma )(1-d)I_t\\ R_{t+1}& =\gamma (1-d)I_t+(1-d)R_t,\end{array}$(7) ), our discrete-time ISAv model implicitly assumes three distinct temporal phases. At the end of each time interval, susceptible salmon populations become infectious, a fraction of infectious salmon die, shedding occurs; a fraction of each salmon class is removed (natural death and virus clearing); and susceptible salmon populations reproduce into the susceptible class. Taking into account the temporal ordering of events, we derive our ISAv model in the following three steps.

1. Disease transmission, disease induced death, and shedding: $\begin{array}{rl}{S}_{\left(1\right)}& =S_t(\theta {\varphi }_I\left(I_t\right)+\hat{\theta }{\varphi }_V\left(V_t\right))\\ I_{\left(1\right)}& =S_t(\theta {\hat{\varphi }}_I\left(I_t\right)+\hat{\theta }{\hat{\varphi }}_V\left(V_t\right))+\hat{\mu }{I}_t\\ V_{\left(1\right)}& =V_t+\delta I_t\end{array}$ That is, after disease transmission, disease induced death, and shedding, $S_{\left(1\right)}$, $I_{\left(1\right)}$ and $V_{\left(1\right)}$ denote sizes of the susceptible salmon, infectious salmon and virus populations respectively.

2. Natural Death (salmon survival and virus clearing): $\begin{array}{rl}{S}_{\left(2\right)}& =\hat{d}{S}_{\left(1\right)}\\ I_{\left(2\right)}& =\hat{d}{I}_{\left(1\right)}\\ V_{\left(2\right)}& ={\hat{d}}_V{V}_{\left(1\right)}\end{array}$ That is, after disease transmission, disease induced death, shedding, and natural death, $S_{\left(2\right)}$, $I_{\left(2\right)}$ and $V_{\left(2\right)}$ denote sizes of the susceptible salmon, infectious salmon and virus populations respectively.

3. Reproduction or Ricker recruitment (S into S): $\begin{array}{rl}{S}_{\left(3\right)}& =r{S}_t{e}^{-b{S}_t}+S_{\left(2\right)}\\ I_{\left(3\right)}& =I_{\left(2\right)}\\ V_{\left(3\right)}& =V_{\left(2\right)}\end{array}$ That is, after disease transmission, disease induced death, shedding, natural death and susceptible salmon reproduction or recruitment, $S_{\left(3\right)}$, $I_{\left(3\right)}$ and $V_{\left(3\right)}$ denote sizes of the susceptible salmon, infectious salmon and virus populations, respectively.

These assumptions and notation lead to the following discrete-time ISAv epidemic model. (20) $\begin{array}{rl}{S}_{t+1}& =r{S}_t{e}^{-b{S}_t}+\hat{d}{S}_t(\theta {\varphi }_I\left(I_t\right)+\hat{\theta }{\varphi }_V\left(V_t\right))\\ I_{t+1}& =\hat{d}(S_t(\theta {\hat{\varphi }}_I\left(I_t\right)+\hat{\theta }{\hat{\varphi }}_V\left(V_t\right))+\hat{\mu }{I}_t)\\ V_{t+1}& ={\hat{d}}_V(V_t+\delta I_t),\end{array}$(20) where $t=0,1,2,\dots$. We study Model (20(20) $\begin{array}{rl}{S}_{t+1}& =r{S}_t{e}^{-b{S}_t}+\hat{d}{S}_t(\theta {\varphi }_I\left(I_t\right)+\hat{\theta }{\varphi }_V\left(V_t\right))\\ I_{t+1}& =\hat{d}(S_t(\theta {\hat{\varphi }}_I\left(I_t\right)+\hat{\theta }{\hat{\varphi }}_V\left(V_t\right))+\hat{\mu }{I}_t)\\ V_{t+1}& ={\hat{d}}_V(V_t+\delta I_t),\end{array}$(20) ) with initial conditions $(S_0,I_0,V_0)\in {\mathbb{R}}_{+}^3$.

Theorem 4.1

In Model (20(20) $\begin{array}{rl}{S}_{t+1}& =r{S}_t{e}^{-b{S}_t}+\hat{d}{S}_t(\theta {\varphi }_I\left(I_t\right)+\hat{\theta }{\varphi }_V\left(V_t\right))\\ I_{t+1}& =\hat{d}(S_t(\theta {\hat{\varphi }}_I\left(I_t\right)+\hat{\theta }{\hat{\varphi }}_V\left(V_t\right))+\hat{\mu }{I}_t)\\ V_{t+1}& ={\hat{d}}_V(V_t+\delta I_t),\end{array}$(20) ), for each time $t\in \{1,2,\dots \},$ $(S_t,I_t,V_t)\in {\mathbb{R}}_{+}^3$ and there is no population explosion whenever $(S_0,I_0,V_0)\in {\mathbb{R}}_{+}^3.$

Proof.

Clearly, $(S_0,I_0,V_0)\in {\mathbb{R}}_{+}^3$ implies $(S_t,I_t,V_t)\in {\mathbb{R}}_{+}^3$ for each time $t\ge 0$. From the equation of $S_{t+1},$ $S_{t+1}\le (r{e}^{-b{S}_t}+(1-d))S_t,$ and the susceptible salmon population is bounded. Furthermore, from the equations of $I_{t+1},$ $I_{t+1}\le (1-d)S_t+(1-d)(1-\mu )I_t.$ Using the boundedness of the susceptible salmon population, we obtain that the infectious salmon population is also bounded. Hence, there is no population explosion in the salmon population of Model (20(20) $\begin{array}{rl}{S}_{t+1}& =r{S}_t{e}^{-b{S}_t}+\hat{d}{S}_t(\theta {\varphi }_I\left(I_t\right)+\hat{\theta }{\varphi }_V\left(V_t\right))\\ I_{t+1}& =\hat{d}(S_t(\theta {\hat{\varphi }}_I\left(I_t\right)+\hat{\theta }{\hat{\varphi }}_V\left(V_t\right))+\hat{\mu }{I}_t)\\ V_{t+1}& ={\hat{d}}_V(V_t+\delta I_t),\end{array}$(20) ). From the equation of $V_{t+1}$, the virus population is also bounded. Hence, the set of iterates of all non-negative initial population sizes are bounded and there is no population explosion in Model (20(20) $\begin{array}{rl}{S}_{t+1}& =r{S}_t{e}^{-b{S}_t}+\hat{d}{S}_t(\theta {\varphi }_I\left(I_t\right)+\hat{\theta }{\varphi }_V\left(V_t\right))\\ I_{t+1}& =\hat{d}(S_t(\theta {\hat{\varphi }}_I\left(I_t\right)+\hat{\theta }{\hat{\varphi }}_V\left(V_t\right))+\hat{\mu }{I}_t)\\ V_{t+1}& ={\hat{d}}_V(V_t+\delta I_t),\end{array}$(20) ).

At each time $t\in \{0,1,2,\dots \}$, when there is no ISAv disease and $I_t=V_t=0$, Model (20(20) $\begin{array}{rl}{S}_{t+1}& =r{S}_t{e}^{-b{S}_t}+\hat{d}{S}_t(\theta {\varphi }_I\left(I_t\right)+\hat{\theta }{\varphi }_V\left(V_t\right))\\ I_{t+1}& =\hat{d}(S_t(\theta {\hat{\varphi }}_I\left(I_t\right)+\hat{\theta }{\hat{\varphi }}_V\left(V_t\right))+\hat{\mu }{I}_t)\\ V_{t+1}& ={\hat{d}}_V(V_t+\delta I_t),\end{array}$(20) ) reduces to Model (15(15) $S_{t+1}=(r{e}^{-b{S}_t}+(1-d))S_t.$(15) ). Notice that when the salmon population of Model (20(20) $\begin{array}{rl}{S}_{t+1}& =r{S}_t{e}^{-b{S}_t}+\hat{d}{S}_t(\theta {\varphi }_I\left(I_t\right)+\hat{\theta }{\varphi }_V\left(V_t\right))\\ I_{t+1}& =\hat{d}(S_t(\theta {\hat{\varphi }}_I\left(I_t\right)+\hat{\theta }{\hat{\varphi }}_V\left(V_t\right))+\hat{\mu }{I}_t)\\ V_{t+1}& ={\hat{d}}_V(V_t+\delta I_t),\end{array}$(20) ) is missing, the virus population decreases to zero. Consequently, ${\mathcal{R}}_d=r/d<1$ implies extinction of the salmon and virus populations in Model (20(20) $\begin{array}{rl}{S}_{t+1}& =r{S}_t{e}^{-b{S}_t}+\hat{d}{S}_t(\theta {\varphi }_I\left(I_t\right)+\hat{\theta }{\varphi }_V\left(V_t\right))\\ I_{t+1}& =\hat{d}(S_t(\theta {\hat{\varphi }}_I\left(I_t\right)+\hat{\theta }{\hat{\varphi }}_V\left(V_t\right))+\hat{\mu }{I}_t)\\ V_{t+1}& ={\hat{d}}_V(V_t+\delta I_t),\end{array}$(20) ). However, ${\mathcal{R}}_d>1$ implies persistence of the salmon population, and the demographic Equation (15(15) $S_{t+1}=(r{e}^{-b{S}_t}+(1-d))S_t.$(15) ) has a unique positive locally asymptotically stable period k salmon population cycle $z_{\mathrm{\infty }}=\{z_{1\mathrm{\infty }},z_{2\mathrm{\infty }},\dots ,z_{k\mathrm{\infty }}\}$, where $k\ge 1.$

4.1. ${\mathcal{R}}_0$ for ISAv model

Assuming virus shedding is not a new infection, we compute the basic reproduction number, ${\mathcal{R}}_0$, for Model (20(20) $\begin{array}{rl}{S}_{t+1}& =r{S}_t{e}^{-b{S}_t}+\hat{d}{S}_t(\theta {\varphi }_I\left(I_t\right)+\hat{\theta }{\varphi }_V\left(V_t\right))\\ I_{t+1}& =\hat{d}(S_t(\theta {\hat{\varphi }}_I\left(I_t\right)+\hat{\theta }{\hat{\varphi }}_V\left(V_t\right))+\hat{\mu }{I}_t)\\ V_{t+1}& ={\hat{d}}_V(V_t+\delta I_t),\end{array}$(20) ) by the next generation matrix method, as in Section 3.1. For each $j\in \{1,2,\dots ,k\}$, $\begin{array}{rl}{F}_j& =\left[\begin{array}{cc}-\hat{d}{z}_{j\mathrm{\infty }}\theta {\varphi }_I^{\mathrm{\prime }}\left(0\right)& -\hat{d}{z}_{j\mathrm{\infty }}\hat{\theta }{\varphi }_V^{\mathrm{\prime }}\left(0\right)\\ 0& 0\end{array}\right]\\ T_j& =\left[\begin{array}{cc}\hat{d}\hat{\mu }& 0\\ {\hat{d}}_V\delta & {\hat{d}}_V\end{array}\right],\end{array}$ and $F_j+T_j=\left[\begin{array}{cc}-\hat{d}{z}_{j\mathrm{\infty }}\theta {\varphi }_I^{\mathrm{\prime }}\left(0\right)+\hat{d}\hat{\mu }& -\hat{d}{z}_{j\mathrm{\infty }}\hat{\theta }{\varphi }_V^{\mathrm{\prime }}\left(0\right)\\ {\hat{d}}_V\delta & {\hat{d}}_V\end{array}\right],$ where the salmon demographic threshold ${\mathcal{R}}_d>1$ and Equation (15(15) $S_{t+1}=(r{e}^{-b{S}_t}+(1-d))S_t.$(15) ) has a unique positive locally asymptotically stable period k salmon population cycle. Hence, the nonnegative transition matrix is $\mathbf{T}\mathbf{=}{T}_k{T}_{k-1}\cdots T_1=\left[\begin{array}{cc}{\left(\hat{d}\right)}^k{\left(\hat{\mu }\right)}^k& 0\\ \delta \sum _{i=0}^{k-1}{\left(\hat{d}\right)}^i{\left(\hat{\mu }\right)}^i{\left({\hat{d}}_V\right)}^{k-i}& {\left({\hat{d}}_V\right)}^k\end{array}\right]$ with $\rho \left(\mathbf{T}\right)<1$, and the nonnegative matrix of new infections is $\mathbf{F}\mathbf{=}(F_k+T_k)(F_{k-1}+T_{k-1})\cdots (F_1+T_1)-\mathbf{T},$ where $\mathbf{F}\mathbf{+}\mathbf{T}$ is irreducible. Consequently, the basic reproduction number for Model (20(20) $\begin{array}{rl}{S}_{t+1}& =r{S}_t{e}^{-b{S}_t}+\hat{d}{S}_t(\theta {\varphi }_I\left(I_t\right)+\hat{\theta }{\varphi }_V\left(V_t\right))\\ I_{t+1}& =\hat{d}(S_t(\theta {\hat{\varphi }}_I\left(I_t\right)+\hat{\theta }{\hat{\varphi }}_V\left(V_t\right))+\hat{\mu }{I}_t)\\ V_{t+1}& ={\hat{d}}_V(V_t+\delta I_t),\end{array}$(20) ) is ${\mathcal{R}}_0=\rho \left(\mathbf{F}\right(Id-\mathbf{T}{)}^{-1}),$ when ${\mathcal{R}}_d>1$ and Equation (15(15) $S_{t+1}=(r{e}^{-b{S}_t}+(1-d))S_t.$(15) ) has a unique positive locally asymptotically stable period k salmon population cycle. By Theorem 2.1, when ${\mathcal{R}}_d>1$ the positive disease-free period k salmon population cycle in Model (20(20) $\begin{array}{rl}{S}_{t+1}& =r{S}_t{e}^{-b{S}_t}+\hat{d}{S}_t(\theta {\varphi }_I\left(I_t\right)+\hat{\theta }{\varphi }_V\left(V_t\right))\\ I_{t+1}& =\hat{d}(S_t(\theta {\hat{\varphi }}_I\left(I_t\right)+\hat{\theta }{\hat{\varphi }}_V\left(V_t\right))+\hat{\mu }{I}_t)\\ V_{t+1}& ={\hat{d}}_V(V_t+\delta I_t),\end{array}$(20) ) is locally asymptotically stable when ${\mathcal{R}}_0<1$ and unstable when ${\mathcal{R}}_0>1$.

When $1<{\mathcal{R}}_d<e^{2/d}$, the disease-free system persists on a fixed point attractor $($non-oscillatory$)$, and Equation (15(15) $S_{t+1}=(r{e}^{-b{S}_t}+(1-d))S_t.$(15) ) has a positive locally asymptotically stable fixed point salmon population (k=1), $z_{\mathrm{\infty }}=\ln {\mathcal{R}}_d/b$. Assuming that virus shedding is not a new infection, and proceeding exactly as before we obtain that at the DFE, $(S,I,V)=(\ln {\mathcal{R}}_d/b,0,0)$, the matrix of new infections that survive the time interval is $\mathbf{F}=\left[\begin{array}{cc}-\frac{\ln {\mathcal{R}}_d}{b}\hat{d}\theta {\varphi }_I^{\mathrm{\prime }}\left(0\right)& -\frac{\ln {\mathcal{R}}_d}{b}\hat{d}\hat{\theta }{\varphi }_V^{\mathrm{\prime }}\left(0\right)\\ 0& 0\end{array}\right]$ the transition matrix is $\mathbf{T}=\left[\begin{array}{cc}\hat{d}\hat{\mu }& 0\\ {\hat{d}}_V\delta & {\hat{d}}_V\end{array}\right],$ and ${\mathcal{R}}_0=\rho \left(\mathbf{F}{(Id-\mathbf{T})}^{-1}\right)={\mathcal{R}}_{0I}+{\mathcal{R}}_{0V}.$ where ${\mathcal{R}}_{0I}=\frac{-\hat{d}\left(\ln {\mathcal{R}}_d\right)\theta {\varphi }_I^{\mathrm{\prime }}\left(0\right)}{b(1-\hat{d}\hat{\mu })}\mathrm{and}{\mathcal{R}}_{0V}=\frac{-\hat{d}{\hat{d}}_V\delta \left(\ln {\mathcal{R}}_d\right)\hat{\theta }{\varphi }_V^{\mathrm{\prime }}\left(0\right)}{b{d}_V(1-\hat{d}\hat{\mu })}.$ As in cholera models, the direct and indirect routes of ISAv transmission enter ${\mathcal{R}}_0$ in an additive way (see, for example, [21,24–26]).

When ${\mathcal{R}}_d>e^{2/d}$ and Equation (15(15) $S_{t+1}=(r{e}^{-b{S}_t}+(1-d))S_t.$(15) ) has a positive locally asymptotically stable period 2 salmon population cycle (k=2), $z_{\mathrm{\infty }}=\{z_{1\mathrm{\infty }},z_{2\mathrm{\infty }}\}$, $\begin{array}{rl}(Id-\mathbf{T}{)}^{-1}& =\left[\begin{array}{cc}\frac{1}{1-{\left(\hat{d}\right)}^2{\left(\hat{\mu }\right)}^2}& 0\\ \frac{\delta \hat{d}\hat{\mu }{\hat{d}}_V+{\left({\hat{d}}_V\right)}^2\delta }{(1-{\left(\hat{d}\right)}^2{\left(\hat{\mu }\right)}^2)(1-{\left({\hat{d}}_V\right)}^2)}& \frac{1}{1-{\left({\hat{d}}_V\right)}^2}\end{array}\right],\\ \mathbf{F}& =\left[\begin{array}{cc}{f}_{11}& f_{12}\\ -\hat{d}{\hat{d}}_V\delta \theta {\varphi }_I^{\mathrm{\prime }}\left(0\right)z_{1\mathrm{\infty }}& -\hat{d}{\hat{d}}_V\delta \hat{\theta }{\varphi }_V^{\mathrm{\prime }}\left(0\right)z_{1\mathrm{\infty }}\end{array}\right],\end{array}$ where $f_{11}={\left(\hat{d}\right)}^2{\theta }^2{\varphi }_I^{\mathrm{\prime }}{\left(0\right)}^2{z}_{1\mathrm{\infty }}{z}_{2\mathrm{\infty }}-{\left(\hat{d}\right)}^2\hat{\mu }\theta {\varphi }_I^{\mathrm{\prime }}\left(0\right)(z_{1\mathrm{\infty }}+z_{2\mathrm{\infty }})-\hat{d}{\hat{d}}_V\delta \hat{\theta }{\varphi }_V^{\mathrm{\prime }}\left(0\right)z_{2\mathrm{\infty }}$ and $f_{12}={\left(\hat{d}\right)}^2\theta {\varphi }_I^{\mathrm{\prime }}\left(0\right)\hat{\theta }{\varphi }_V^{\mathrm{\prime }}\left(0\right)z_{1\mathrm{\infty }}{z}_{2\mathrm{\infty }}-{\left(\hat{d}\right)}^2\hat{\mu }\hat{\theta }{\varphi }_V^{\mathrm{\prime }}\left(0\right)z_{1\mathrm{\infty }}-\hat{d}{\hat{d}}_V\hat{\theta }{\varphi }_V^{\mathrm{\prime }}\left(0\right)z_{2\mathrm{\infty }}.$

4.2. ISAv model and demographic period 2 cycle

In this section, we consider Model (20(20) $\begin{array}{rl}{S}_{t+1}& =r{S}_t{e}^{-b{S}_t}+\hat{d}{S}_t(\theta {\varphi }_I\left(I_t\right)+\hat{\theta }{\varphi }_V\left(V_t\right))\\ I_{t+1}& =\hat{d}(S_t(\theta {\hat{\varphi }}_I\left(I_t\right)+\hat{\theta }{\hat{\varphi }}_V\left(V_t\right))+\hat{\mu }{I}_t)\\ V_{t+1}& ={\hat{d}}_V(V_t+\delta I_t),\end{array}$(20) ) with b=0.1, d=0.5 and $r=e^4\approx 54.598$, where the model's unit of time is a year. The disease-free equation is Model (15(15) $S_{t+1}=(r{e}^{-b{S}_t}+(1-d))S_t.$(15) ). From Figure 1, the disease-free system has a locally asymptotically stable period 2 population cycle, $z_{\mathrm{\infty }}=\{z_{1\mathrm{\infty }},z_{2\mathrm{\infty }}\}=\{37.883,65.757\}$. To compute ${\mathcal{R}}_0$, we assume that the ISAv disease infections are modelled as Poisson processes with ${\varphi }_I\left(I_t\right)=\exp (-{\beta }_I{I}_t)$ and ${\varphi }_V\left(V_t\right)=\exp (-{\beta }_V{V}_t)$ with varying ${\beta }_I$, and set the following parameter values in Model (20(20) $\begin{array}{rl}{S}_{t+1}& =r{S}_t{e}^{-b{S}_t}+\hat{d}{S}_t(\theta {\varphi }_I\left(I_t\right)+\hat{\theta }{\varphi }_V\left(V_t\right))\\ I_{t+1}& =\hat{d}(S_t(\theta {\hat{\varphi }}_I\left(I_t\right)+\hat{\theta }{\hat{\varphi }}_V\left(V_t\right))+\hat{\mu }{I}_t)\\ V_{t+1}& ={\hat{d}}_V(V_t+\delta I_t),\end{array}$(20) ) $\begin{array}{rl}& b=0.1,d=0.5,{\beta }_V=0.01,d_V=0.2,r=e^4,\delta =0.2,\\ & \mu =0.7\mathrm{and}\theta =\mathrm{0.6.}\end{array}$ As the infection rate of infectious salmon, ${\beta }_I$ increases from zero to values greater than ${\beta }_I\approx 0.046$, ${\mathcal{R}}_0$ increases from values less than 1 (disease extinction, see Figure 5) to values greater than 1 (disease persistence on a period or chaotic population cycle, see Figure 5). As predicted by Theorem 2.1, Figure 5 shows that the disease-free period 2 population cycle is locally asymptotically stable when ${\mathcal{R}}_0<1$ and unstable when ${\mathcal{R}}_0>1$.

Figure 5 shows ISAv disease induced period-doubling bifurcations route to chaos in the infectious salmon population, where the disease-free susceptible salmon population is on a locally asymptotically stable period 2 cycle. For example, when ${\beta }_I=0.056$ and all the other parameters are kept fixed at their current values in Figure 5, then ${\mathcal{R}}_0=1.178>1$ and Figure 6 shows that Model (7(7) $\begin{array}{rl}{S}_{t+1}& =g\left(N_t\right)+(1-d)S_t\varphi \left(I_t\right)\\ I_{t+1}& =(1-d)S_t(1-\varphi \left(I_t\right))+(1-\gamma )(1-d)I_t\\ R_{t+1}& =\gamma (1-d)I_t+(1-d)R_t,\end{array}$(7) ) has a locally asymptotically stable endemic period 4 population cycle at $\begin{array}{rl}& \{(37.060,10.969,12.818)\to (62.263,7.640,13.764)\to \\ & (29.748,9.249,13.456)\to (93.442,5.744,13.735)\},\end{array}$ where the period of the disease-free susceptible salmon population cycle is 2. That is, unlike in Model (7(7) $\begin{array}{rl}{S}_{t+1}& =g\left(N_t\right)+(1-d)S_t\varphi \left(I_t\right)\\ I_{t+1}& =(1-d)S_t(1-\varphi \left(I_t\right))+(1-\gamma )(1-d)I_t\\ R_{t+1}& =\gamma (1-d)I_t+(1-d)R_t,\end{array}$(7) ), in Model (20(20) $\begin{array}{rl}{S}_{t+1}& =r{S}_t{e}^{-b{S}_t}+\hat{d}{S}_t(\theta {\varphi }_I\left(I_t\right)+\hat{\theta }{\varphi }_V\left(V_t\right))\\ I_{t+1}& =\hat{d}(S_t(\theta {\hat{\varphi }}_I\left(I_t\right)+\hat{\theta }{\hat{\varphi }}_V\left(V_t\right))+\hat{\mu }{I}_t)\\ V_{t+1}& ={\hat{d}}_V(V_t+\delta I_t),\end{array}$(20) ) the equation for the total salmon population, $S_{t+1}+I_{t+1}$, is different from the disease-free equation, and the period of the disease-free susceptible salmon population cycle does not determine the period of the infectious salmon population cycle when ${\mathcal{R}}_0>1$.

Figure 5. Bifurcation diagram for ISAv Model (20(20) $\begin{array}{rl}{S}_{t+1}& =r{S}_t{e}^{-b{S}_t}+\hat{d}{S}_t(\theta {\varphi }_I\left(I_t\right)+\hat{\theta }{\varphi }_V\left(V_t\right))\\ I_{t+1}& =\hat{d}(S_t(\theta {\hat{\varphi }}_I\left(I_t\right)+\hat{\theta }{\hat{\varphi }}_V\left(V_t\right))+\hat{\mu }{I}_t)\\ V_{t+1}& ={\hat{d}}_V(V_t+\delta I_t),\end{array}$(20) ): As ${\mathcal{R}}_0$ increases from values less than 1 to values greater than 1, the model dynamics changes from disease extinction to disease persistence on a periodic or chaotic population cycle.

Figure 6. Model (20(20) $\begin{array}{rl}{S}_{t+1}& =r{S}_t{e}^{-b{S}_t}+\hat{d}{S}_t(\theta {\varphi }_I\left(I_t\right)+\hat{\theta }{\varphi }_V\left(V_t\right))\\ I_{t+1}& =\hat{d}(S_t(\theta {\hat{\varphi }}_I\left(I_t\right)+\hat{\theta }{\hat{\varphi }}_V\left(V_t\right))+\hat{\mu }{I}_t)\\ V_{t+1}& ={\hat{d}}_V(V_t+\delta I_t),\end{array}$(20) ) with $g\left(S_t\right)=r{S}_t{e}^{-b{S}_t}$ has a locally asymptotically stable period 4 population cycle when ${\beta }_I=0.056$ and all the other parameter values are kept fixed at their current values in Figure 5.

5. Concluding remarks

Populations do not grow indefinitely over time, and density dependence or other factors tend to drive populations toward their carrying capacity. However, it is possible for real populations to experience some population fluctuations around the carrying capacity. These fluctuations can be periodic (cyclic) or erratic (chaotic). In this paper, we use a general autonomous discrete-time infectious disease model to extend the next generation matrix approach for calculating the basic reproduction number, ${\mathcal{R}}_0$, to account for populations with cyclic demographic dynamics. When ${\mathcal{R}}_0<1$ and the demographic equation (in the absence of the disease) has a locally asymptotically stable period k population cycle, we prove the local asymptotic stability of the disease-free period k cycle where $k\ge 1$. That is, the disease dies out whenever ${\mathcal{R}}_0<1$. Also, under the same period k demographic assumption, we prove that the disease-free period k population cycle is unstable and the disease persists when ${\mathcal{R}}_0>1$ (Theorem 2.1). An extension of the next generation matrix approach for calculating the basic reproduction number, ${\mathcal{R}}_0$, to account for populations with chaotic and non-periodic demographic dynamics is an open problem.

We apply our theoretical results to discrete-time infectious disease models with a Ricker recruitment function that are formulated for SIR infections with and without vaccination, and ISAv infections in a salmon population. When ${\mathcal{R}}_0>1$, in the SIR model, the period of the disease-free S-cycle is numerically shown to be the same as the period of the I-cycle. That is, the demographic dynamics drives the SIR disease dynamics (see Figure 2). However, we illustrate, in Figure 5, that the demographic dynamics does not drive the ISAv disease dynamics. Questions on establishing the possible connections between ISAv disease infections and observed population cycles in some salmon populations are open.

Acknowledgments

We thank the referees for their useful comments and suggestions. Preliminary ideas for this research were discussed during ICMA-VI Conference at University of Arizona in Tucson.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This research was partially supported by NSERC, through a Discovery Grant (P.vdD.). A.-A.Y. was partially supported by DHS Center Of Excellence for Command, Control and Interoperability at Rutgers University, NSF Computational Sustainability Grant # CCF - 1522054, and NSF Award # DMS-1743144.