A Lyapunov–Schmidt method for detecting backward bifurcation in age-structured population models

Abstract

Backward bifurcation is an important property of infectious disease models. A centre manifold method has been developed by Castillo-Chavez and Song for detecting the presence of backward bifurcation and deriving a necessary and sufficient condition for its occurrence in Ordinary Differential Equations (ODE) models. In this paper, we extend this method to partial differential equation systems. First, we state a main theorem. Next we illustrate the application of the new method on a chronological age-structured Susceptible-Infected-Susceptible (SIS) model with density-dependent recovery rate, an age-since-infection structured HIV/AIDS model with standard incidence and an age-since-infection structured cholera model with vaccination.

Keywords:

• Lyapunov–Schmidt theory
• backward bifurcation
• age-structured population model
• epidemic model

AMS Subject Classification:

• 92D30

1. Introduction

Backward bifurcation plays an important role in infectious diseases and their control. Having methods for deriving necessary and sufficient conditions for backward bifurcation can delineate the reasons for its occurrence and help devise control strategies that combat its negative effects. A recent review article [10] collects methods for analysis of backward bifurcation in ODE models and age-structured partial differential equation (PDE) models. The reader may find some mechanisms that create backward bifurcations for endemic steady states in [7].

Most methods were first developed or applied in ODEs but have an analogue for age-structured PDEs. A very popular method for detecting backward bifurcation in ODE models was developed by Castillo-Chavez and Song [2]. Their method is based on the centre manifold theorem and is one of the preferred methods for modern analysis of backward bifurcation in ODE population models. As developed in [2], however, this method cannot be applied to PDEs and does not work for delay equations, age-structured population models, or diffusion models.

In this paper, our goal is to extend the method to infinite dimensional systems. In the next section, we state a theorem analogous to Castillo-Chavez and Song's theorem for infinite dimensional systems, which gives the necessary and sufficient condition for backward bifurcations to occur. In one of our remarks, we argue that the backward bifurcating solution is unstable, at least close to the bifurcating point, while the forward bifurcating solution is locally stable, also at least close to the bifurcating point. This is an expected universal property for transcritical bifurcation. In the rest of the paper, we focus on concrete epidemic examples. The first example we consider is a chronological age-structured SIS model with non-linear recovery rate. For this example, the theory applies directly because the boundary condition is linear, so the differential operator becomes a generator of a $C_0$-semigroup. We also illustrate this example with a numerical simulation that shows the backward bifurcation, once the condition for it has been met. Example 2 is an age-since-infection HIV/AIDS model with standard incidence. It is known from the literature that this model does not exhibit backward bifurcation and the result of the analysis concludes exactly that. The importance of this example is the development of the tools to handle the non-linear, non-local boundary condition associated with age-since-infection structured models. Example 3 is an age-since-infection structured cholera model considered as a test example in [10], where we derive a necessary and sufficient condition for backward bifurcation using several different methods. Here we use the newly proposed Lyapunov–Schmidt method and we reach a similar condition for backward bifurcation. Compared to the methods in [10] the Lyapunov–Schmidt method developed here is somewhat easier to use but it is based on more abstract approach.

We anticipate that the methodology introduced here can be applied to other infinite dimensional models, such as diffusion epidemic models, or delay epidemic models. We expect that this methodology will be helpful to pave the road for deriving necessary and sufficient conditions for backward bifurcation in infinite dimensional systems.

2. Lyapunov–Schmidt method for PDEs

In this section, we introduce an analogue of the Castillo-Chavez and Song Theorem [2] designed for PDEs. The proof is based on the Lyapunov–Schmidt theory.

Let X and Z be Banach spaces, and $x\in X$ and $\alpha \in \mathbb{R}$ is a parameter. We consider the following abstract differential equation (1) $\frac{\mathrm{d}x}{\mathrm{d}t}=F(x,\alpha ),F:X\times \mathbb{R}⟶Z.$(1) Without loss of generality, we assume that 0 is an equilibrium so $F(0,\alpha )=0\mathrm{for}\mathrm{all}\alpha \in \mathbb{R}.$

Theorem 2.1

Assume:

• $\mathcal{A}=D_{x}F(0,{\alpha }_0)$ is the linearization around the equilibrium 0, evaluated at a critical value of parameter ${\alpha }_0,$ such that $\mathcal{A}$ is a closed operator with a simple isolated eigenvalue zero and remaining eigenvalues having negative real part. Let ${\hat{v}}_0$ be the unique (up to a constant) positive solution of $\mathcal{A}v=0$.

• $F(x,\alpha )\in C^2(U_0\times I_0,Z)$ for some neighbourhood $U_0$ of 0 and interval $I_0$ containing ${\alpha }_0$.

• Assume $Z^{\ast }$ is the dual of Z and $⟨\cdot ,\cdot ⟩$ is the pairing between Z and $Z^{\ast }$. Assume ${\hat{v}}_0^{\ast }\in Z^{\ast }$ is the unique (up to a constant) positive vector satisfying $⟨\mathcal{A}x,{\hat{v}}_0^{\ast }⟩=0\text{for all }x\in X,$ that is $\dim (\ker {\mathcal{A}}^{\ast })=1$, where ${\mathcal{A}}^{\ast }$ is the adjoint of $\mathcal{A}$, and $\ker {\mathcal{A}}^{\ast }=\mathrm{span}\{{\hat{v}}_0^{\ast }\}$.

• Assume $⟨D_{x\alpha }^{2}F(0,{\alpha }_0){\hat{v}}_0,{\hat{v}}_0^{\ast }⟩\ne 0$ where $D_{x\alpha }F(0,{\alpha }_0)$ is the second derivative of F with respect to x and α.

Then, the direction of the bifurcation is determined by the numbers (2) $\begin{array}{rl}a& =⟨D_{xx}^{2}F(0,{\alpha }_0)[{\hat{v}}_0,{\hat{v}}_0],{\hat{v}}_0^{\ast }⟩,\\ b& =⟨D_{x\alpha }^{2}F(0,{\alpha }_0){\hat{v}}_0,{\hat{v}}_0^{\ast }⟩,\end{array}$(2) where $D_{xx}^{2}F(0,{\alpha }_0)[h_1,h_2]$ is the second derivative of F with respect to x applied to the functions $h_1$ and $h_2$. If $b>0,$ the bifurcation is backward if and only if a>0 and forward if and only if a<0.

Proof Theorem 2.1.

(Theorem 1) To determine the direction of the bifurcation, we study the solution of $F(x,\alpha )=0$ in a neighbourhood of $(0,{\alpha }_0)\in X\times \mathbb{R}$. Let $\mathcal{A}=D_{x}F(0,{\alpha }_0).$ Since $\mathcal{A}$ is closed and has a simple isolated eigenvalue, it follows from [4, p. 284] that $\mathcal{A}$ is a Fredholm operator of index zero, that is $\mathcal{A}$ satisfies:

• $\dim (\ker \mathcal{A})=1$;

• $\dim (Z/\mathrm{im}\mathcal{A})=1$;

• $\mathrm{im}\mathcal{A}$ is closed in Z

where $\ker \mathcal{A}=\{x\in X:\mathcal{A}x=0\},$ and $\mathrm{im}\mathcal{A}=\{y\in Z:y=\mathcal{A}x\text{ for some }x\in X\}.$ Because, $\dim (Z/\mathrm{im}\mathcal{A})=1$, we have $Z/\mathrm{im}\mathcal{A}=\mathrm{span}\{[{\hat{z}}_0]\}$ where $[{\hat{z}}_0]={\hat{z}}_0+\mathrm{im}\mathcal{A}$ and ${\hat{z}}_0\notin \mathrm{im}\mathcal{A}$. Let $Z_0=\mathrm{span}\{{\hat{z}}_0\}$. Then following [3] we have: $Z=Z_0\oplus \mathrm{im}\mathcal{A}.$ Similarly, if $\ker \mathcal{A}=\mathrm{span}\{{\hat{v}}_0\}$, then $X=\ker \mathcal{A}\oplus X_0,$ where $X_0$ is the topological complement of $\ker \mathcal{A}$ in X.

It follows that there exist projections $P\in \mathcal{L}(X,\ker \mathcal{A})$ and $Q\in \mathcal{L}(Z,Z_0)$ such that $\ker \mathcal{A}=\mathrm{im}P$ and $\mathrm{im}\mathcal{A}=\ker Q$. Consequently we can write $x\in X$ in the form $x=v+y$ where $v=Px\in \ker \mathcal{A}$. That implies that there exists a constant c such that $v=c{\hat{v}}_0$. In addition $y=(I-P)x\in X_0$. Using the projection Q we can decompose $F(x,\alpha )=0$ into two equations (3) $\begin{array}{rl}& (I-Q)F(v+y,\alpha )=0,\\ & QF(v+y,\alpha )=0.\end{array}$(3) Let $\mathcal{G}(v,y,\alpha ):=(I-Q)F(v+y,\alpha )$. Then the first equation in (3(3) $\begin{array}{rl}& (I-Q)F(v+y,\alpha )=0,\\ & QF(v+y,\alpha )=0.\end{array}$(3) ) is written as $\mathcal{G}(v,y,\alpha )=0$. In particular, we have $\mathcal{G}(0,0,{\alpha }_0)=(I-Q)F(0,{\alpha }_0)=0$. Furthermore, $D_y\mathcal{G}(0,0,{\alpha }_0)=(I-Q)D_{x}F(0,{\alpha }_0)=(I-Q)\mathcal{A}=\mathcal{A}.$ The restriction $\mathcal{A}{|}_{X_0}$ is a one-to-one mapping between $X_0$ and $\mathrm{im}\mathcal{A}$ and it is therefore invertible. Hence $D_y\mathcal{G}(0,0,{\alpha }_0){|}_{X_0}$ is invertible. As a result we can apply the Implicit Function Theorem. From the Implicit Function Theorem, it follows that there exists a neighbourhood $\mathrm{\Omega }\in \ker \mathcal{A}$ with $0\in \mathrm{\Omega }$ and an interval $I\in \mathrm{\Re }$ with ${\alpha }_0\in I$ and a $C^2$ map $\psi :\mathrm{\Omega }\times I⟶X_0$ with $y=\psi (v,\alpha )$ such that $\psi (0,{\alpha }_0)=0$ and $(I-Q)F(v+\psi (v,\alpha ),\alpha )\equiv 0\mathrm{for}\mathrm{all}(v,\alpha )\in \mathrm{\Omega }\times I.$ Furthermore, we define a map $\mathcal{B}:\mathrm{\Omega }\times I⟶Z_0$ by $\mathcal{B}(v,\alpha ):=QF(v+\psi (v,\alpha ),\alpha ).$ Then the second equation in (3(3) $\begin{array}{rl}& (I-Q)F(v+y,\alpha )=0,\\ & QF(v+y,\alpha )=0.\end{array}$(3) ) can be written as (4) $\mathcal{B}(v,\alpha )=0,$(4) which is called the bifurcation equation.

The Lyapunov–Schmidt reduction says that there exists a neighbourhood U of $0\in \mathrm{\Omega }$ such that each solution of $\mathcal{B}(v,\alpha )=0$ corresponds one-to-one to some solution of $F(x,\alpha )=0$. Furthermore, it is shown in [8, Corollary I.2.4] that (5) $\mathcal{B}(0,{\alpha }_0)=0,{\mathcal{B}}_v(0,{\alpha }_0)=0.$(5) From $F(0,\alpha )=0$ for all $\alpha \in \mathrm{\Re }$, we have that $\psi (0,\alpha )=0\mathrm{for}\mathrm{all}\alpha \in I.$ Furthermore, $D_{\alpha }\psi (0,\alpha )=0\mathrm{for}\mathrm{all}\alpha \in I.$ Consequently, $\mathcal{B}(0,\alpha )=0\mathrm{for}\mathrm{all}\alpha \in I.$ This corresponds to the trivial solution. This last equality means that we can write $\mathcal{B}$ in the following way [8]: $\mathcal{B}(v,\alpha )={\int }_0^1\frac{\mathrm{d}}{\mathrm{d}t}\mathcal{B}(tv,\alpha )\mathrm{d}t.$ Alternatively, $\mathcal{B}(v,\alpha )={\int }_0^1{D}_v\mathcal{B}(tv,\alpha )v\mathrm{d}t\mathrm{for}\mathrm{all}(v,\alpha )\in \mathrm{\Omega }\times I.$ For $v\in \mathrm{\Omega }$, we set $v=s{\hat{v}}_0$, $s\in (-\delta ,\delta )$ and we obtain the nontrivial solutions for $s\ne 0$ by solving $\tilde{\mathcal{B}}(s,\alpha )\equiv {\int }_0^1{D}_v\mathcal{B}(ts{\hat{v}}_0,\alpha ){\hat{v}}_0\mathrm{d}t=0\mathrm{for}s\in (-\delta ,\delta ),s\ne 0.$ The assumptions on the smoothness of F give that $\tilde{\mathcal{B}}\in C^1((-\delta ,\delta )\times I,Z_0)$. Furthermore, by the properties (5(5) $\mathcal{B}(0,{\alpha }_0)=0,{\mathcal{B}}_v(0,{\alpha }_0)=0.$(5) ) we have $\tilde{\mathcal{B}}(0,{\alpha }_0)=0$. We would like to apply the Implicit Function Theorem and express $\alpha =\varphi (s)$. To do that, we need $D_{\alpha }\tilde{\mathcal{B}}(0,{\alpha }_0)\ne 0$. Computation similar to [8, p. 19] give that (6) $\begin{array}{rl}{D}_{\alpha }(D_v\mathcal{B}(v,\alpha ){\hat{v}}_0)& =D_{\alpha }(Q{D}_{x}F(v+\psi (v,\alpha ),\alpha )({\hat{v}}_0+D_v\psi (v,\alpha ){\hat{v}}_0)\\ & =Q{D}_{xx}^{2}F(v+\psi (v,\alpha )\alpha ,\alpha )[{\hat{v}}_0+D_v\psi (v,\alpha ){\hat{v}}_0,D_{\alpha }\psi (v,\alpha )]\\ & +Q{D}_{x}F(v+\psi (v,\alpha ),\alpha )D_{v\alpha }^2\psi (v,\alpha ){\hat{v}}_0\\ & +Q{D}_{x\alpha }^{2}F(v+\psi (v,\alpha ),\alpha )({\hat{v}}_0+D_v\psi (v,\alpha ){\hat{v}}_0).\end{array}$(6) Setting $(v,\alpha )=(0,{\alpha }_0)$ we find that the first term is equal to zero from the properties of ψ and the second term is also zero by the definition of the projection Q. Hence, we obtain (7) $D_{\alpha }\tilde{\mathcal{B}}(0,{\alpha }_0)=Q{D}_{x\alpha }^{2}F(0,{\alpha }_0){\hat{v}}_0,$(7) and $Q{D}_{x\alpha }^{2}F(0,{\alpha }_0){\hat{v}}_0\ne 0$ because of assumptions A3 and A4. In particular, $D_{x\alpha }^{2}F(0,{\alpha }_0){\hat{v}}_0\ne 0$ and $D_{x\alpha }^{2}F(0,{\alpha }_0){\hat{v}}_0\notin \mathrm{im}\mathcal{A}$. The Implicit Function Theorem now gives the existence of $\varphi (s)$ such that $\varphi :(-\delta ,\delta )⟶I$ with $\varphi (0)={\alpha }_0$ and $\tilde{\mathcal{B}}(s,\varphi (s))=0\mathrm{for}\mathrm{all}s\in (-\delta ,\delta ).$ Then $\mathcal{B}(s{\hat{v}}_0,\varphi (s))=s\tilde{\mathcal{B}}(s,\varphi (s))=0$ for $s\in (-\delta ,\delta )$. Consequently, $F(x(s),\alpha (s))=0\mathrm{for}s\in (-\delta ,\delta )$ where (8) $\begin{array}{rl}x(s)& =s{\hat{v}}_0+\psi (s{\hat{v}}_0,\varphi (s)),\\ \alpha (s)& =\varphi (s).\end{array}$(8) This is the essence of the Crandall-Rabinowitz Theorem [8]. Next we claim that the tangent vector to the nontrivial solution curve (8(8) $\begin{array}{rl}x(s)& =s{\hat{v}}_0+\psi (s{\hat{v}}_0,\varphi (s)),\\ \alpha (s)& =\varphi (s).\end{array}$(8) ) at the bifurcation point $(0,{\alpha }_0)$ is given by $({\hat{v}}_0,{\alpha }^{\mathrm{\prime }}(0))\in X\times \mathbb{R}.$ Differentiating $x(s)$ in (8(8) $\begin{array}{rl}x(s)& =s{\hat{v}}_0+\psi (s{\hat{v}}_0,\varphi (s)),\\ \alpha (s)& =\varphi (s).\end{array}$(8) ) with respect to s gives (9) $\frac{\mathrm{d}}{\mathrm{d}s}x(s){|}_{s=0}={\hat{v}}_0+D_v\psi (0,{\alpha }_0){\hat{v}}_0+D_{\alpha }\psi (0,{\alpha }_0)\dot{\alpha }(0)={\hat{v}}_0.$(9) The bifurcation at the critical value $(0,{\alpha }_0)$ is backward if and only if ${\alpha }^{\mathrm{\prime }}(0)<0$. To obtain formulas to compute ${\alpha }^{\mathrm{\prime }}(0)$, we recall that $\alpha (s)=\varphi (s)$ and therefore $\tilde{\mathcal{B}}(s,\alpha (s))=0$ for all $s\in (-\delta ,\delta )$. Thus, (10) $\frac{\mathrm{d}}{\mathrm{d}s}\tilde{\mathcal{B}}(s,\alpha (s)){|}_{s=0}=D_s\tilde{\mathcal{B}}(0,{\alpha }_0)+D_{\alpha }\tilde{\mathcal{B}}(0,{\alpha }_0)\dot{\alpha }(0)=0,$(10) where $D_{\alpha }\tilde{\mathcal{B}}(0,{\alpha }_0)\ne 0$ is given in (7(7) $D_{\alpha }\tilde{\mathcal{B}}(0,{\alpha }_0)=Q{D}_{x\alpha }^{2}F(0,{\alpha }_0){\hat{v}}_0,$(7) ). Furthermore, (11) $\begin{array}{rl}{D}_s\tilde{\mathcal{B}}(0,{\alpha }_0)& ={\int }_0^1{D}_{vv}^2\mathcal{B}(0,{\alpha }_0)[{\hat{v}}_0,t{\hat{v}}_0]\mathrm{d}t\\ & =Q{D}_{xx}^{2}F(0,{\alpha }_0)[{\hat{v}}_0,{\hat{v}}_0]{\int }_0^{1}t\mathrm{d}t=\frac{1}{2}Q{D}_{xx}^{2}F(0,{\alpha }_0)[{\hat{v}}_0,{\hat{v}}_0].\end{array}$(11) Using the dual pairing and from Equation (10(10) $\frac{\mathrm{d}}{\mathrm{d}s}\tilde{\mathcal{B}}(s,\alpha (s)){|}_{s=0}=D_s\tilde{\mathcal{B}}(0,{\alpha }_0)+D_{\alpha }\tilde{\mathcal{B}}(0,{\alpha }_0)\dot{\alpha }(0)=0,$(10) ) we have $\frac{1}{2}⟨Q{D}_{xx}^{2}F(0,{\alpha }_0)[{\hat{v}}_0,{\hat{v}}_0],{\hat{v}}_0^{\ast }⟩+⟨Q{D}_{x\alpha }^{2}F(0,{\alpha }_0)[{\hat{v}}_0,{\hat{v}}_0],{\hat{v}}_0^{\ast }⟩\dot{\alpha }(0)=0.$ Notice that (12) $\begin{array}{rl}& ⟨Q{D}_{xx}^{2}F(0,{\alpha }_0)[{\hat{v}}_0,{\hat{v}}_0],{\hat{v}}_0^{\ast }⟩=⟨D_{xx}^{2}F(0,{\alpha }_0)[{\hat{v}}_0,{\hat{v}}_0],{\hat{v}}_0^{\ast }⟩,\\ & ⟨Q{D}_{x\alpha }^{2}F(0,{\alpha }_0){\hat{v}}_0,{\hat{v}}_0^{\ast }⟩=⟨D_{x\alpha }^{2}F(0,{\alpha }_0){\hat{v}}_0,{\hat{v}}_0^{\ast }⟩.\end{array}$(12) These equalities hold because $⟨z,{\hat{v}}_0^{\ast }⟩=0$ for all $z\in \mathrm{im}\mathcal{A}$. Hence, we have ${\alpha }^{\mathrm{\prime }}(0)=-\frac{1}{2}\frac{⟨D_{xx}^{2}F(0,{\alpha }_0)[{\hat{v}}_0,{\hat{v}}_0],{\hat{v}}_0^{\ast }⟩}{⟨D_{x\alpha }^{2}F(0,{\alpha }_0){\hat{v}}_0,{\hat{v}}_0^{\ast }⟩}=-\frac{1}{2}\frac{a}{b}.$ If we assume that b>0, the bifurcation curve bifurcates backward if ${\alpha }^{\mathrm{\prime }}(0)<0$, that is, if a>0. If ${\alpha }^{\mathrm{\prime }}(0)>0$, that is a<0, then the bifurcation is forwards.

Remark 2.1

The condition b>0 means that ${\mathcal{R}}_0$ increases with the parameter α. If b<0, that means that ${\mathcal{R}}_0$ decreases as α increases and here we will not consider this case as a backward bifurcation case.

Remark 2.2

As is usually observed in the ODE case, stability of the bifurcating solutions can be determined such that the backwardly bifurcated solution is unstable, while the forwardly bifurcated solution is locally stable as long as the parameter s is small enough. According to Britton [1], let us show the Factorization Theorem. Let $(x(s),\alpha (s))$ be the bifurcating solutions as (13) $F(x(s),\alpha (s))=0,s\in (-\delta ,\delta ),(x(0),\alpha (0))=(0,{\alpha }_0),$(13) and $\frac{\mathrm{d}}{\mathrm{d}s}x(s){|}_{s=0}={\hat{v}}_0.$ Define $\mathcal{A}(s):=D_{x}F(x(s),\alpha (s)).$ Then $\mathcal{A}(0)=\mathcal{A}$. The dominant eigenvalue $\sigma (s)$ and corresponding eigenvector are introduced as follows: (14) $\mathcal{A}(s){\hat{v}}_0(s)=\hat{\sigma }(s){\hat{v}}_0(s),$(14) where $\hat{\sigma }(0)=0,{\hat{v}}_0(0)={\hat{v}}_0.$ Define the dual eigenvector ${\hat{v}}_0^{\ast }(s)$ as ${\mathcal{A}}^{\ast }(s){\hat{v}}_0^{\ast }(s)=\hat{\sigma }(s){\hat{v}}_0^{\ast }(s),$ where they are normalized as follows: $⟨{\hat{v}}_0(s),{\hat{v}}_0^{\ast }(s)⟩=1.$ Especially, we write as ${\hat{v}}_0^{\ast }(0)={\hat{v}}_0^{\ast }$.

Differentiating (13(13) $F(x(s),\alpha (s))=0,s\in (-\delta ,\delta ),(x(0),\alpha (0))=(0,{\alpha }_0),$(13) ) with respect to s, we have (15) $D_{x}F(x,\alpha )x^{\prime }(s)+D_{\alpha }F(x,\alpha ){\alpha }^{\prime }(s)=0$(15) From (15(15) $D_{x}F(x,\alpha )x^{\prime }(s)+D_{\alpha }F(x,\alpha ){\alpha }^{\prime }(s)=0$(15) ), taking the duality pairing with ${\hat{v}}_0^{\ast }(s)$, it follows that (16) $\begin{array}{rl}⟨D_{x}F(x,\alpha )x^{\prime }(s),{\hat{v}}_0^{\ast }(s)⟩& =⟨\mathcal{A}(s)x^{\prime }(s),{\hat{v}}_0^{\ast }(s)⟩\\ & =⟨x^{\prime }(s),{\mathcal{A}}^{\ast }(s){\hat{v}}_0^{\ast }(s)⟩\\ & =\hat{\sigma }(s)⟨x^{\prime }(s),{\hat{v}}_0^{\ast }(s)⟩\\ & =-⟨D_{\alpha }F(x,\alpha ){\alpha }^{\prime }(s),{\hat{v}}_0^{\ast }(s)⟩\\ & =-{\alpha }^{\prime }(s)⟨D_{\alpha }F(x,\alpha ),{\hat{v}}_0^{\ast }(s)⟩.\end{array}$(16) Then we arrive at an expression (17) $\hat{\sigma }(s)=-{\alpha }^{\prime }(s)\frac{⟨D_{\alpha }F(x,\alpha ),{\hat{v}}_0^{\ast }(s)⟩}{⟨x^{\prime }(s),{\hat{v}}_0^{\ast }(s)⟩}$(17) Using the Taylor expansion around s = 0, we have $\begin{array}{rl}& D_{\alpha }F(x,\alpha )=D_x{D}_{\alpha }F(0,{\alpha }_0)x^{\prime }(s)s+O(s^2),\\ & \alpha (s)={\alpha }_0+{\alpha }^{\prime }(0)s+O(s^2),\\ & ⟨x^{\prime }(s),{\hat{v}}_0^{\ast }(s)⟩=1+O(s),\end{array}$ where we used $D_{\alpha }F(0,{\alpha }_0)=0$. Then it holds that (18) $\begin{array}{rl}\hat{\sigma }(s)& =-(\alpha (s)-{\alpha }_0)⟨D_x{D}_{\alpha }F(0,{\alpha }_0){\hat{v}}_0,{\hat{v}}_0^{\ast }⟩+O(s^2)\\ & =-(\alpha (s)-{\alpha }_0)b+O(s^2).\end{array}$(18) Therefore, if we assume that b>0, for small s>0, the backward bifurcation (${\alpha }^{\prime }(0)<0$) implies the bifurcated solution is unstable ($\hat{\sigma }(s)>0$), while the forward bifurcation (${\alpha }^{\prime }(0)>0$) implies the bifurcated solution is locally stable ($\hat{\sigma }(s)<0$) as long as we could assume that $\sigma (s)$ is the dominant eigenvalue. In many concrete examples, we could show the dominance of $\sigma (s)$ for small s using complex analysis arguments applied to the characteristic equation.

In the following sections, we illustrate how the above theory can be applied to concrete examples from epidemic modeling.

3. Example 1: age-structured SIS epidemic model

3.1. An SIS epidemic model with density-dependent recovery rate

First we consider the following SIS model with age-structure and non-linear recovery rate: (19) $\begin{array}{rl}{s}_t+s_a& =-\kappa (a)s(a,t){\int }_0^{\omega }\beta (\sigma )i(\sigma ,t)\mathrm{d}\sigma +\frac{\gamma (a)i(a,t)}{A+{\int }_0^{\omega }i(a,t)\mathrm{d}a}\\ i_t+i_a& =\kappa (a)s(a,t){\int }_0^{\omega }\beta (\sigma )i(\sigma ,t)\mathrm{d}\sigma -\frac{\gamma (a)i(a,t)}{A+{\int }_0^{\omega }i(a,t)\mathrm{d}a}\end{array}$(19) where a is chronological age, t is time, and ω is the maximal age, assumed finite. Furthermore, $s(a,t)$ and $i(a,t)$ are the age-density functions of susceptibles and invectives, respectively, $\kappa (a)$ is the susceptibility of susceptible individuals, $\beta (a)$ is the age-specific infectivity of infective individuals, $\gamma (a)$ is the standard recovery rate, and A is a non-zero constant. That is, it is assumed that the effective recovery rate is density-dependent, it is decreasing as the size of infecteds increases. Although this is a toy model to illustrate the bifurcation calculation, it would not necessarily be unrealistic, if the recovery needs medical supports and the medical resource is limited. Here we assume that the time scale of epidemic is so short that the natural death rate of host population can be neglected¹.

We assume that $\kappa (a),\beta (a),\gamma (a)\in L^{\mathrm{\infty }}(0,\omega ).$ Although model (19(19) $\begin{array}{rl}{s}_t+s_a& =-\kappa (a)s(a,t){\int }_0^{\omega }\beta (\sigma )i(\sigma ,t)\mathrm{d}\sigma +\frac{\gamma (a)i(a,t)}{A+{\int }_0^{\omega }i(a,t)\mathrm{d}a}\\ i_t+i_a& =\kappa (a)s(a,t){\int }_0^{\omega }\beta (\sigma )i(\sigma ,t)\mathrm{d}\sigma -\frac{\gamma (a)i(a,t)}{A+{\int }_0^{\omega }i(a,t)\mathrm{d}a}\end{array}$(19) ) is with separable transmission rate, this approach can be applied in more general situations.

The model has the following initial and boundary conditions: (20) $\begin{array}{rl}s(0,t)& =B\\ i(0,t)& =0\\ s(a,0)& =s_0(a)\\ i(a,0)& =i_0(a)\end{array}$(20) where B is the birth rate, assumed constant. Furthermore, we assume that $s_0(a)+i_0(a)=B$ for all $a\in [0,\omega ]$. The equation of the total population size is $n_t+n_a=0$. Hence, the total population is constant $n(a,t)=B$. Let, $s(a,t)=B-i(a,t)$. Thus, we reduce the system to a unique equation: (21) $\begin{array}{rl}& i_t+i_a=\kappa (a)(B-i(a,t)){\int }_0^{\omega }\beta (\sigma )i(\sigma ,t)\mathrm{d}\sigma -\frac{\gamma (a)i(a,t)}{A+{\int }_0^{\omega }i(a,t)\mathrm{d}a},\\ & i(0,t)=0.\end{array}$(21)

3.2. Application of the Lyapunov–Schmidt method

Now, we consider the equation for the equilibria which is the base equation to which we will apply the Lyapunov–Schmidt backward bifurcation theorem. (22) $-\frac{\mathrm{d}i}{\mathrm{d}a}+\kappa (a)(B-i(a)){\int }_0^{\omega }\beta (\sigma )i(\sigma )\mathrm{d}\sigma -\frac{\gamma (a)i(a)}{A+{\int }_0^{\omega }i(a)\mathrm{d}a}=0,i(0)=0.$(22) This equation clearly has the disease-free equilibrium $i^0(a)=0$. Let us introduce a bifurcation parameter $\overline{\beta }$ as $\beta (a)=\overline{\beta }{\beta }_0(a),$ where a normalized function ${\beta }_0\in L_{+}^1(0,\omega )\cap L_{+}^{\mathrm{\infty }}(0,\omega )$ will be determined below. Thus, we have a problem $F(i,\overline{\beta })=0,$ where $F(i,\overline{\beta }):=(Ai)(a)+\kappa (a)(B-i(a))\overline{\beta }{\int }_0^{\omega }{\beta }_0(\sigma )i(\sigma )\mathrm{d}\sigma -\frac{\gamma (a)i(a)}{A+{\int }_0^{\omega }i(a)\mathrm{d}a},$ where A is the differential operator given by $(Af)(a)=-f^{\prime }(a)$ with domain $\mathcal{D}(A)=\{f\in AC(0,\omega ):f(0)=0\},$ where $AC(0,\omega )$ is the set of absolutely continuous functions on the interval $(0,\omega )$. The operator F is defined from $\mathcal{D}(A)\times \mathbb{R}$ to $L^1(0,\omega )$.

To compute $D_{i}F(0,\overline{\beta })$ we consider a parameter s and a function $h(a)$ with properties to be specified later and we compose $\phi (s)=F(i+sh,\overline{\beta })$. Taking ${\phi }^{\prime }(s)$ and setting s = 0 we obtain (23) $\begin{array}{rl}{D}_{i}F(i,\overline{\beta })h& =(Ah)(a)-\kappa (a)h(a)\overline{\beta }{\int }_0^{\omega }{\beta }_0(\sigma )i(\sigma )\mathrm{d}\sigma +\overline{\beta }\kappa (a)(B-i){\int }_0^{\omega }{\beta }_0(\sigma )h(\sigma )\mathrm{d}\sigma \\ {\displaystyle }& -\frac{\gamma (a)h(a)}{A+{\int }_0^{\omega }i(a)\mathrm{d}a}+\frac{\gamma (a)i(a)}{(A+{\int }_0^{\omega }i(a)\mathrm{d}a{)}^2}{\int }_0^{\omega }h(\sigma )\mathrm{d}\sigma .\end{array}$(23) Setting i = 0, we have $D_{i}F(0,\overline{\beta })h=(Ah)(a)+\overline{\beta }\kappa (a)B{\int }_0^{\omega }{\beta }_0(\sigma )h(\sigma )\mathrm{d}\sigma -\frac{\gamma (a)h(a)}{A}.$ This is not exactly the operator $\mathcal{A}$ since the derivative is evaluated at an arbitrary value of the parameter $\overline{\beta }$.

To see the stability of the disease-free equilibrium, we look for eigenvalues of the operator $D_{i}F(0,\overline{\beta })h=\lambda h$. This leads to the linearized equation of the disease-free equilibrium: (24) $h_a+\lambda h+\frac{\gamma (a)}{A}h=\kappa (a)B\mathcal{H},h(0)=0,$(24) where $\mathcal{H}={\int }_0^{\omega }\beta (\sigma )h(\sigma )\mathrm{d}\sigma .$ Solving the above differential equation we have $h(a)=B\mathcal{H}{\int }_0^a{\mathrm{e}}^{-\lambda (a-s)-{\int }_s^a(\gamma (\eta )/A)\mathrm{d}\eta }\kappa (s)\mathrm{d}s.$ Substituting $h(a)$ in $\mathcal{H}$ and cancelling since we are looking for non-zero h, we obtain the following characteristic equation: $1=B{\int }_0^{\omega }\beta (a){\int }_0^a{\mathrm{e}}^{-\lambda (a-s)-{\int }_s^a(\gamma (\eta )/A)\mathrm{d}\eta }\kappa (s)\mathrm{d}s\mathrm{d}a.$ From here we can define the basic reproduction number ${\mathcal{R}}_0=\overline{\beta }B{\int }_0^{\omega }{\beta }_0(a){\int }_0^a{\mathrm{e}}^{-{\int }_s^a(\gamma (\eta )/A)\mathrm{d}\eta }\kappa (s)\mathrm{d}s\mathrm{d}a.$ The normalized function ${\beta }_0$ is assumed to satisfy (25) $1=B{\int }_0^{\omega }{\beta }_0(a){\int }_0^a{\mathrm{e}}^{-{\int }_s^a(\gamma (\eta )/A)\mathrm{d}\eta }\kappa (s)\mathrm{d}s\mathrm{d}a.$(25) Then we have ${\mathcal{R}}_0=\overline{\beta }$ and ${\lambda }^{\ast }=0$ is a simple isolated eigenvalue of the characteristic equation when $\overline{\beta }=1$. It is not hard to see that all remaining eigenvalues have negative real part.

We define the operator $\mathcal{A}$: $\mathcal{A}=D_{i}F(0,1),$ which has the same domain as $\mathcal{D}(A)$. Looking for $\ker \mathcal{A}$, we need to solve $\mathcal{A}v=0.$ Following a procedure similar to obtaining $h(a)$, we obtain ${\hat{v}}_0(a)=B{\int }_0^a{\mathrm{e}}^{-{\int }_s^a(\gamma (\eta )/A)\mathrm{d}\eta }\kappa (s)\mathrm{d}s.$ We notice that from (25(25) $1=B{\int }_0^{\omega }{\beta }_0(a){\int }_0^a{\mathrm{e}}^{-{\int }_s^a(\gamma (\eta )/A)\mathrm{d}\eta }\kappa (s)\mathrm{d}s\mathrm{d}a.$(25) ) we have ${\int }_0^{\omega }{\beta }_0(a){\hat{v}}_0(a)\mathrm{d}a=1.$ To obtain ${\hat{v}}_0^{\ast }$ we need to formally find the adjoint operator, and obtain its kernel. Since $Z=L^1(0,\omega )$, the dual space is $Z^{\ast }=L^{\mathrm{\infty }}(0,\omega )$. Let ${\mathcal{A}}^{\ast }$ be the adjoint with $\mathcal{D}({\mathcal{A}}^{\ast })=\{\psi (a)\in AC(0,\omega ):\psi (\omega )=0\}.$ For $\psi \in \mathcal{D}({\mathcal{A}}^{\ast })$, we have after integration by parts (26) $\begin{array}{rl}& {\int }_0^{\omega }[-h_a+\kappa (a)B{\int }_0^{\omega }{\beta }_0(s)h(s)\mathrm{d}s-\frac{\gamma (a)}{A}h(a)]\psi (a)\mathrm{d}a\\ & ={\int }_0^{\omega }[{\psi }_a+B{\int }_0^{\omega }\kappa (s)\psi (s)\mathrm{d}s{\beta }_0(a)-\frac{\gamma (a)}{A}\psi (a)]h(a)\mathrm{d}a.\end{array}$(26) Hence, ${\mathcal{A}}^{\ast }\psi ={\psi }_a+B{\int }_0^{\omega }\kappa (s)\psi (s)\mathrm{d}s{\beta }_0(a)-\frac{\gamma (a)}{A}\psi (a)$ with $\psi (\omega )=0$. To find ${\hat{v}}_0^{\ast }$ we need to determine the kernel of the adjoint operator. For that purpose, we have to solve the problem: (27) ${\psi }_a+B{\beta }_0(a)\mathcal{K}-\frac{\gamma (a)}{A}\psi =0,\psi (\omega )=0,$(27) where $\mathcal{K}={\int }_0^{\omega }\kappa (s)\psi (s)\mathrm{d}s.$ Solving this system as before, we obtain ${\hat{v}}_0^{\ast }(a)=B{\int }_a^{\omega }{\beta }_0(a){\mathrm{e}}^{{\int }_s^a(\gamma (\eta )/A)\mathrm{d}\eta }\mathrm{d}s.$ It is not hard to see from integration by parts that ${\int }_0^{\omega }\kappa (a){\hat{v}}_0^{\ast }(a)\mathrm{d}a=1.$ Next, we compute a and b. First we compute b. We have, $D_{i\overline{\beta }}^{2}F(0,1){\hat{v}}_0=B\kappa (a){\int }_0^{\omega }{\beta }_0(s){\hat{v}}_0(s)\mathrm{d}s.$ Hence, we have (28) $b=B{\int }_0^{\omega }\kappa (a){\hat{v}}_0^{\ast }(a)\mathrm{d}a{\int }_0^{\omega }{\beta }_0(s){\hat{v}}_0(s)\mathrm{d}s=B>0.$(28) Next, we compute a. First, we need to compute $D_{ii}^{2}F(0,\overline{\beta })[h,h_2]$. To do that we consider (29) $\begin{array}{rl}{\phi }_2(s)& =D_{i}F(i+s{h}_2,\overline{\beta })h\\ & =-h_a-\overline{\beta }\kappa (a)h{\int }_0^{\omega }{\beta }_0(a)(i+s{h}_2)\mathrm{d}a+\overline{\beta }\kappa (a)(B-i-s{h}_2){\int }_0^{\omega }{\beta }_0(a)h(a)\mathrm{d}a\\ & -\frac{\gamma (a)h(a)}{A+{\int }_0^{\omega }(i+s{h}_2)\mathrm{d}a}+\frac{\gamma (a)(i+s{h}_2(a))}{(A+{\int }_0^{\omega }(i+s{h}_2)\mathrm{d}a{)}^2}{\int }_0^{\omega }h(\eta )\mathrm{d}\eta .\end{array}$(29) Differentiating with respect to s and setting s = 0 and i = 0, we obtain (30) $\begin{array}{rl}{D}_{ii}^{2}F(0,\overline{\beta })[h,h_2]& =-\overline{\beta }\kappa (a)h(a){\int }_0^{\omega }{\beta }_0(a)h_2(a)\mathrm{d}a-\overline{\beta }\kappa (a)h_2(a){\int }_0^{\omega }{\beta }_0(a)h(a)\mathrm{d}a\\ & +\frac{\gamma (a)h(a)}{A^2}{\int }_0^{\omega }{h}_2(a)\mathrm{d}a+\frac{\gamma (a)h_2(a)}{A^2}{\int }_0^{\omega }h(a)\mathrm{d}a.\end{array}$(30) We notice that this form is symmetric with respect to h and $h_2$. Thus, we have $D_{ii}^{2}F(0,1)[{\hat{v}}_0,{\hat{v}}_0]=-2\kappa (a){\hat{v}}_0(a){\int }_0^{\omega }{\beta }_0(a){\hat{v}}_0(a)\mathrm{d}a+2\frac{\gamma (a){\hat{v}}_0(a)}{A^2}{\int }_0^{\omega }{\hat{v}}_0(a)\mathrm{d}a.$ If we recall ${\int }_0^{\omega }{\beta }_0(a){\hat{v}}_0(a)da=1$, multiplying with ${\hat{v}}_0^{\ast }$ and integrating, we obtain (31) $a=-2{\int }_0^{\omega }\kappa (a){\hat{v}}_0(a){\hat{v}}_0^{\ast }(a)\mathrm{d}a+\frac{2}{A^2}{\int }_0^{\omega }\gamma (a){\hat{v}}_0(a){\hat{v}}_0^{\ast }(a)\mathrm{d}a{\int }_0^{\omega }{\hat{v}}_0(a)\mathrm{d}a.$(31) Hence, the following is a necessary and sufficient condition for backward bifurcation (a>0) (32) ${\int }_0^{\omega }\kappa (a){\hat{v}}_0(a){\hat{v}}_0^{\ast }(a)\mathrm{d}a<{\int }_0^{\omega }\frac{\gamma (a){\hat{v}}_0(a){\hat{v}}_0^{\ast }}{A^2}\mathrm{d}a{\int }_0^{\omega }{\hat{v}}_0(a)\mathrm{d}a.$(32)

4. Example 2: infection-age structured HIV/AIDS epidemic model

The example in this section is an HIV/AIDS epidemic model developed in [11]. Although it is already known that this model does not possess backward bifurcation, it is instructive to apply the technique on this rather simple model under formal setting before we apply it to more complex models. In fact, if we extend this model to take into account chronological age structure, we can show that the backward bifurcation of endemic steady states can occur [5,6], and many age-structured epidemic model can be formulated as the same kind of the Cauchy problem with non-densely defined generator and non-local boundary condition.

Consider the infection-age (age-since-infection) structured HIV/AIDS model with the standard incidence: (33) $\begin{array}{rl}& S^{\prime }=\mathrm{\Lambda }-\overline{\beta }\frac{S}{N}{\int }_0^{\mathrm{\infty }}{\beta }_0(\tau )i(\tau ,t)\mathrm{d}\tau -\mu S,\\ & i_t+i_{\tau }=-(\mu +\alpha (\tau ))i,\\ & i(0,t)=\overline{\beta }\frac{S(t)}{N(t)}{\int }_0^{\mathrm{\infty }}{\beta }_0(\tau )i(\tau ,t)\mathrm{d}\tau ,\end{array}$(33) where $S(t)$ is the number of susceptibles, $i(\tau ,t)$ is the infection-age density of infectives and $N(t)=S(t)+{\int }_0^{\mathrm{\infty }}i(\tau ,t)\mathrm{d}\tau .$ is the total size of sexually active population, which is not constant because of the disease-induced death rate $\alpha (\tau )$. We assume that ${\beta }_0(\tau )$ is normalized as (34) ${\int }_0^{\mathrm{\infty }}{\beta }_0(\tau )\pi (\tau )\mathrm{d}\tau =1,$(34) where $\pi (\tau ):=\exp (-{\int }_0^{\tau }(\mu +\alpha (\sigma ))\mathrm{d}\sigma ),$ denotes the survival rate of infecteds, so we can use $\overline{\beta }$ as a bifurcation parameter, and $\overline{\beta }$ is no other than the basic reproduction number.

According to [9], in order to take into account the boundary condition, we define the extended state space $X:=\mathbb{R}\times \hat{X}$ associated with $\hat{X}:=L^1(0,\mathrm{\infty })\times \mathbb{R}$. We also define the subspace by $X_0:=\mathbb{R}\times {\hat{X}}_0$, where ${\hat{X}}_0:=L^1(0,\mathrm{\infty })\times \{0\}$. Then $X_0$ can be identified with the usual state space $\mathbb{R}\times L^1(0,\mathrm{\infty })$. Let $E_0=(\mathrm{\Lambda }/\mu ,\hat{0})\in X$ be the disease-free steady state, where $\hat{0}=(0,0{)}^{\mathrm{T}}\in {\hat{X}}_0$. Define a linear operator $\hat{A}:\mathcal{D}(\hat{A})\subset {\hat{X}}_0\to \hat{X}$ by $\hat{A}\left(\begin{array}{c}u\\ 0\end{array}\right)=\left(\begin{array}{c}-u^{\prime }-(\mu +\alpha )u\\ -u(0)\end{array}\right),$ for $u\in W^{1,1}(0,\mathrm{\infty })$ with $\mathcal{D}(\hat{A})=W^{1,1}(0,\mathrm{\infty })\times \{0\}\subset {\hat{X}}_0$.

Let $\hat{i}(t)=(i(t,\cdot ),0{)}^{\mathrm{T}}\in {\hat{X}}_0$. Then the basic problem (36(36) $\frac{\mathrm{d}x}{\mathrm{d}t}=Ax+G(x),x(0)\in X_0,t\ge 0,$(36) ) can be formulated as a semilinear Cauchy problem in X: (35) $\begin{array}{rl}\frac{\mathrm{d}S(t)}{\mathrm{d}t}& =-\mu S(t)+G_1(S(t),\hat{i}(t)),\\ \frac{\mathrm{d}\hat{i}(t)}{\mathrm{d}t}& =\hat{A}\hat{i}(t)+G_2(S(t),\hat{i}(t)),\end{array}$(35) where for $x=(x_1,{\hat{x}}_2{)}^{\mathrm{T}}\in X_0$, $G_1:X_0\to \mathbb{R}$ and $G_2:X_0\to \hat{X}$ are defined by $\begin{array}{rl}{G}_1(x):& =\mathrm{\Lambda }-\overline{\beta }\frac{x_1}{N(x)}{\int }_0^{\mathrm{\infty }}{\beta }_0(\tau )x_2(\tau )\mathrm{d}\tau ,\\ G_2(x):& =\left(\begin{array}{c}0\\ \overline{\beta }\frac{x_1}{N(x)}{\int }_0^{\mathrm{\infty }}{\beta }_0(\tau )x_2(\tau )\mathrm{d}\tau \end{array}\right)\end{array}$ and $N(x)$ is a positive functional of $x=(x_1,{\hat{x}}_2{)}^{\mathrm{T}}\in X_0$ defined by $N(x):=x_1+{\int }_0^{\mathrm{\infty }}{x}_2(\tau )\mathrm{d}\tau .$ Moreover, define a linear operator $A:\mathcal{D}(A)\subset X_0\to X$ and a nonlinear operator $G:X_0\to X$ by $A\left(\begin{array}{c}{x}_1\\ {\hat{x}}_2\end{array}\right)=\left(\begin{array}{c}-\mu x_1\\ \hat{A}{\hat{x}}_2\end{array}\right),G(x):=\left(\begin{array}{c}{G}_1(x)\\ G_2(x)\end{array}\right)$ with $\mathcal{D}(A)=\mathbb{R}\times \mathcal{D}(\hat{A})$. Then $\overline{\mathcal{D}(A)}=X_0$ is not dense in X, but A satisfies the Hille-Yosida estimate. Then (35(35) $\begin{array}{rl}\frac{\mathrm{d}S(t)}{\mathrm{d}t}& =-\mu S(t)+G_1(S(t),\hat{i}(t)),\\ \frac{\mathrm{d}\hat{i}(t)}{\mathrm{d}t}& =\hat{A}\hat{i}(t)+G_2(S(t),\hat{i}(t)),\end{array}$(35) ) can be formulated as a well-posed abstract Cauchy problem in X as follows: (36) $\frac{\mathrm{d}x}{\mathrm{d}t}=Ax+G(x),x(0)\in X_0,t\ge 0,$(36) where $x=(x_1,{\hat{x}}_2{)}^{\mathrm{T}}\in X_0$.

In order to examine the bifurcation at the disease-free steady state of (39(39) ${\psi }_2(\tau )={\psi }_2(0){\int }_{\tau }^{\mathrm{\infty }}{\beta }_0(s)\frac{\pi (s)}{\pi (\tau )}\mathrm{d}s.$(39) ), let $F(x,\overline{\beta }):=Ax+G(x),$ for $x=(x_1,{\hat{x}}_2)\in X_0$. Then $F(x_0,\overline{\beta })=0$ if $x_0=(\mathrm{\Lambda }/\mu ,\hat{0})$ is the disease-free steady state. Then the linearized operator $\mathcal{A}:=D_{x}F(x_0,\overline{\beta })$ acting on X is calculated as follows: (37) $\mathcal{A}x:=\left(\begin{array}{c}-\mu x_1-\overline{\beta }{\displaystyle {\int }_0^{\mathrm{\infty }}{\beta }_0(\tau )x_2(\tau )\mathrm{d}\tau }\\ -x_2^{\prime }(\tau )-(\mu +\alpha (\tau ))x_2(\tau )\\ -x_2(0)+\overline{\beta }{\displaystyle {\int }_0^{\mathrm{\infty }}{\beta }_0(\tau )x_2(\tau )\mathrm{d}\tau }\end{array}\right),$(37) where its domain is given by $\mathcal{D}(\mathcal{A})=\mathcal{D}(A)=\mathbb{R}\times W^{1,1}(0,\mathrm{\infty })\times \{0\}.$ Let us consider the eigenvalues of $\mathcal{A}$. Solving the differential equation from $\mathcal{A}x=\lambda x$ for $x\in X_0$, we have $x_2(\tau )=x_2(0)\pi (\tau ){\mathrm{e}}^{-\lambda \tau }.$ Inserting the above expression into the third equation in $\mathcal{A}x=\lambda x$, we have $-x_2(0)+x_2(0)\overline{\beta }{\int }_0^{\mathrm{\infty }}{\beta }_0(\tau )\pi (\tau ){\mathrm{e}}^{-\lambda \tau }\mathrm{d}\tau =0.$ Then we obtain the characteristic equation (38) $1=\overline{\beta }{\int }_0^{\mathrm{\infty }}{\beta }_0(\tau )\pi (\tau ){\mathrm{e}}^{-\lambda \tau }\mathrm{d}\tau ,$(38) which has a simple eigenvalue of zero when $\overline{\beta }=1$. To find the eigenvector ${\hat{v}}_0$ associated with eigenvalue zero, we need to solve $\mathcal{A}x=0$. Then we can choose $x_2(\tau )=\pi (\tau )$ and the boundary condition (the third equation of $\mathcal{A}x=0$) is trivially satisfied at $\overline{\beta }=1$. From the first equation of $\mathcal{A}x=0$, we have $x_1=-(1/\mu )$ and ${\hat{v}}_0=(-\frac{1}{\mu },\pi (\tau ),1).$ To find the adjoint operator ${\mathcal{A}}^{\ast }$, let $x=(x_1,{\hat{x}}_2{)}^{\mathrm{T}}\in \mathcal{D}(\mathcal{A})\subset X_0$, we multiply $\mathcal{A}x$ with the vector $\psi =({\psi }_1,{\psi }_2,{\psi }_3{)}^{\mathrm{T}}\in X^{\ast }:=\mathbb{R}\times L^{\mathrm{\infty }}(0,\mathrm{\infty })\times \mathbb{R}$. Then we have $\begin{array}{rl}⟨\psi ,\mathcal{A}x⟩& =[-\mu x_1-{\psi }_1\overline{\beta }{\int }_0^{\mathrm{\infty }}{\beta }_0(\tau )x_2(\tau )\mathrm{d}\tau ]{\psi }_1\\ & +{\int }_0^{\mathrm{\infty }}(-x_2^{\prime }(\tau )-(\mu +\alpha (\tau ))x_2(\tau )){\psi }_2(\tau )\mathrm{d}\tau \\ & +[-x_2(0)+\overline{\beta }{\int }_0^{\mathrm{\infty }}{\beta }_0(\tau )x_2(\tau )\mathrm{d}\tau ]{\psi }_3.\end{array}$ Observe that $\begin{array}{rl}& {\int }_0^{\mathrm{\infty }}(-x_2^{\prime }(\tau )-(\mu +\alpha (\tau ))x_2(\tau )){\psi }_2(\tau )\mathrm{d}\tau \\ & =-x_2(\tau ){\psi }_2(\tau ){|}_0^{\mathrm{\infty }}+{\int }_0^{\mathrm{\infty }}{x}_2(\tau ){\psi }_2^{\prime }(\tau {)}_{\tau }\mathrm{d}\tau -{\int }_0^{\mathrm{\infty }}(\mu +\alpha (\tau ))x_2(\tau ){\psi }_2(\tau )\mathrm{d}\tau \\ & =x_2(0){\psi }_2(0)+{\int }_0^{\mathrm{\infty }}[{\psi }_2^{\prime }(\tau )-(\mu +\alpha (\tau )){\psi }_2(\tau )]x_2(\tau )\mathrm{d}\tau ,\end{array}$ where we assume that ${\psi }_2(\mathrm{\infty })=0$. Hence, we have $\begin{array}{rl}⟨\psi ,\mathcal{A}x⟩& =-\mu x_1{\psi }_1+{\int }_0^{\mathrm{\infty }}(\overline{\beta }{\beta }_0(\tau )(-{\psi }_1+{\psi }_2(0))+{\psi }_2^{\prime }(\tau )\\ & -(\mu +\alpha (\tau )){\psi }_2(\tau ))x_2(\tau )\mathrm{d}\tau \\ & +x_2(0)(-{\psi }_3+{\psi }_2(0))\\ & =⟨{\mathcal{A}}^{\ast }\psi ,x⟩,\end{array}$ which should hold for $\psi \in \mathcal{D}({\mathcal{A}}^{\ast })$ and $x\in \mathcal{D}(\mathcal{A})\subset X_0$.

Thus we choose the domain of ${\mathcal{A}}^{\ast }$ as $\mathcal{D}({\mathcal{A}}^{\ast })=\{\psi =({\psi }_1,{\psi }_2,{\psi }_3)\in \mathbb{R}\times W^{1,\mathrm{\infty }}\times \mathbb{R}:{\psi }_2(\mathrm{\infty })=0,{\psi }_3={\psi }_2(0)\}.$ Then the adjoint operator ${\mathcal{A}}^{\ast }$ is defined on $\mathcal{D}({\mathcal{A}}^{\ast })\subset X^{\ast }=\mathbb{R}\times L^{\mathrm{\infty }}\times \mathbb{R}$ as follows: ${\mathcal{A}}^{\ast }\psi =\left(\begin{array}{c}-\mu {\psi }_1\\ \overline{\beta }{\beta }_0(-{\psi }_1+{\psi }_2(0))+{\psi }_2^{\prime }(\tau )-(\mu +\alpha (\tau )){\psi }_2\\ {\psi }_2(0)\end{array}\right),\psi \in D({\mathcal{A}}^{\ast }).$ To find ${\hat{v}}_0^{\ast }$, we solve ${\mathcal{A}}^{\ast }\psi =0$ with $\overline{\beta }=1$. We have ${\psi }_1=0$. Solving the differential equation in ${\mathcal{A}}^{\ast }\psi =0$, we have (39) ${\psi }_2(\tau )={\psi }_2(0){\int }_{\tau }^{\mathrm{\infty }}{\beta }_0(s)\frac{\pi (s)}{\pi (\tau )}\mathrm{d}s.$(39) Hence, we have the eigenfunction ${\hat{v}}_0^{\ast }=(0,{\psi }_2(\tau ),{\psi }_2(0))$, where ${\psi }_2(0)>0$.

Computing the second derivative for F, we have $D_x{D}_{\overline{\beta }}F(x_0,1){\hat{v}}_0=\left(\begin{array}{c}-{\displaystyle {\int }_0^{\mathrm{\infty }}{\beta }_0(\tau )\pi (\tau )\mathrm{d}\tau }\\ 0\\ {\displaystyle {\int }_0^{\mathrm{\infty }}{\beta }_0(\tau )\pi (\tau )\mathrm{d}\tau }\end{array}\right).$ Hence, we obtain (40) $b=⟨D_x{D}_{\overline{\beta }}F(x_0,1){\hat{v}}_0,{\hat{v}}_0^{\ast }⟩={\psi }_2(0){\int }_0^{\mathrm{\infty }}{\beta }_0(\tau )\pi (\tau )\mathrm{d}\tau >0.$(40) On the other hand, $D_{xx}^{2}F(x_0,1)[{\hat{v}}_0,{\hat{v}}_0]$ can be computed as follows: $D_{xx}^{2}F(x_0,1)[{\hat{v}}_0,{\hat{v}}_0]={\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }h\mathrm{\partial }k}F(x_0+(h+k){\hat{v}}_0,1)|}_{(h,k)=(0,0)}=\left(\begin{array}{c}\frac{2\mu }{\mathrm{\Lambda }}{\displaystyle {\int }_0^{\mathrm{\infty }}\pi (\tau )\mathrm{d}\tau }\\ 0\\ -{\displaystyle \frac{2\mu }{\mathrm{\Lambda }}{\int }_0^{\mathrm{\infty }}\pi (\tau )\mathrm{d}\tau }\end{array}\right).$ Thus, we have (41) $a=⟨D_{xx}^{2}F(x_0,1)[{\hat{v}}_0,{\hat{v}}_0],{\hat{v}}_0^{\ast }⟩=-{\psi }_2(0)\frac{2\mu }{\mathrm{\Lambda }}{\int }_0^{\mathrm{\infty }}\pi (\tau )\mathrm{d}\tau <0.$(41) Hence, it follows from (41(41) $a=⟨D_{xx}^{2}F(x_0,1)[{\hat{v}}_0,{\hat{v}}_0],{\hat{v}}_0^{\ast }⟩=-{\psi }_2(0)\frac{2\mu }{\mathrm{\Lambda }}{\int }_0^{\mathrm{\infty }}\pi (\tau )\mathrm{d}\tau <0.$(41) ) and (40(40) $b=⟨D_x{D}_{\overline{\beta }}F(x_0,1){\hat{v}}_0,{\hat{v}}_0^{\ast }⟩={\psi }_2(0){\int }_0^{\mathrm{\infty }}{\beta }_0(\tau )\pi (\tau )\mathrm{d}\tau >0.$(40) ) that $-\frac{1}{2}a/b>0$, so the bifurcation at $(x,\overline{\beta })=(x_0,1)$ is a forward one, which is a conclusion that we expected.

Remark 4.1

In this example, just as in the ODE case [2], we see that components of the eigenvector ${\hat{v}}_0$, that correspond to positive entries in the disease-free equilibrium, may be negative. This is in general the case with the PDEs as well.

5. Example 3: age-since-infection structured cholera model

In this section, we consider a cholera model as an age-since-infection example with backward bifurcation. The model was first introduced in [10]. Here we use the model and formal calculations along Theorem 1 to derive the necessary and sufficient condition for backward bifurcation associated with this model. In particular, if τ is the age-since-infection of infectious individuals, and θ is the age of the bacteria in the environment, we denote by $i(\tau ,t)$ the density of the infectives and by $u(\theta ,t)$ the density of the bacteria in the environment.

Let also S be the number of susceptibles, V the number of vaccinated, and R the number of recovered individuals. With these notations, the age-structured model becomes: (42) $\begin{array}{rl}{\displaystyle S^{\prime }}& =\mathrm{\Lambda }-\frac{S}{B+D}{\int }_0^{\mathrm{\infty }}\beta (\theta )u(\theta ,t)\mathrm{d}\theta -(\mu +\psi )S+wR,\\ V^{\prime }& =\psi S-\sigma \frac{V}{B+D}{\int }_0^{\mathrm{\infty }}\beta (\theta )u(\theta ,t)\mathrm{d}\theta -\mu V,\\ i_t+i_{\tau }& =-(\mu +\gamma (\tau ))i,\\ i(0,t)& =\frac{S+\sigma V}{B+D}{\int }_0^{\mathrm{\infty }}\beta (\theta )u(\theta ,t)\mathrm{d}\theta ,\\ R^{\prime }& ={\int }_0^{\mathrm{\infty }}\gamma (\tau )i(\tau ,t)\mathrm{d}\tau -(\mu +w)R,\\ u_{\theta }+u_t& =-\delta (\theta )u,\\ u(0,t)& ={\int }_0^{\mathrm{\infty }}\eta (\tau )i(\tau ,t)\mathrm{d}\tau .\end{array}$(42) Again the basic system can be formulated as a Cauchy problem with non-densely defined generators, so here we skip repeating the abstract formulation by using the extended state space.

The parameters are described in [10]; β and δ are functions of θ, while γ and η are structured by τ. B in this model denotes the concentration of bacteria in the environment $B(t)={\int }_0^{\mathrm{\infty }}u(\theta ,t)\mathrm{d}\theta .$ This system has a disease-free equilibrium ${\mathcal{E}}_0=(S^0,V^0,0,0,0)$ where $S^0=\frac{\mathrm{\Lambda }}{\mu +\psi }{V}^0=\frac{\psi }{\mu }{S}^0.$ We define the probability of survival of infected individual in the infectious class and the probability of survival of the bacteria in the environment: $\pi (\tau )={\mathrm{e}}^{-\mu \tau +{\int }_0^{\tau }\gamma (s)\mathrm{d}s}{\pi }_{\mathrm{d}}(\theta )={\mathrm{e}}^{-{\int }_0^{\theta }\delta (s)\mathrm{d}s}.$ The control reproduction number is given by (derivation can be found in [10]): (43) $\mathcal{R}(\psi )=\frac{(S^0+\sigma V^0)}{D}{\int }_0^{\mathrm{\infty }}\beta (\theta ){\pi }_{\mathrm{d}}(\theta )\mathrm{d}\theta {\int }_0^{\mathrm{\infty }}\eta (\tau )\pi (\tau )\mathrm{d}\tau .$(43) To have a bifurcation parameter, we decompose $\beta (\tau )$ as $\overline{\beta }{\beta }_0(\tau )$ where $max{\beta }_0(\tau )=1$. We define $X=\mathbb{R}\times \mathbb{R}\times L^1(0,\mathrm{\infty })\times \mathbb{R}\times \mathbb{R}\times L^1(0,\mathrm{\infty })\times \mathbb{R}$ and $X^0=X=\mathbb{R}\times \mathbb{R}\times L^1(0,\mathrm{\infty })\times \{0\}\times \mathbb{R}\times L^1(0,\mathrm{\infty })\times \{0\}$. For an element $\varphi \in X^0$ with $\varphi =(S,V,i,0,R,u,0)$ we define the non-linear operator $F(\varphi ,\overline{\beta })=\left(\begin{array}{c}{\displaystyle \mathrm{\Lambda }-\overline{\beta }\frac{S}{B+D}{\int }_0^{\mathrm{\infty }}{\beta }_0(\theta )u(\theta )\mathrm{d}\theta -(\mu +\psi )S+wR}\\ {\displaystyle \psi S-\sigma \overline{\beta }\frac{V}{B+D}{\int }_0^{\mathrm{\infty }}{\beta }_0(\theta )u(\theta )\mathrm{d}\theta -\mu V}\\ {\displaystyle -i_{\tau }-(\mu +\gamma (\tau ))i}\\ {\displaystyle \overline{\beta }\frac{S+\sigma V}{B+D}{\int }_0^{\mathrm{\infty }}{\beta }_0(\theta )u(\theta )\mathrm{d}\theta -i(0)}\\ {\displaystyle {\int }_0^{\mathrm{\infty }}\gamma (\tau )i(\tau )\mathrm{d}\tau -(\mu +w)R}\\ {\displaystyle -u_{\theta }-\delta (\theta )u}\\ {\displaystyle {\int }_0^{\mathrm{\infty }}\eta (\tau )i(\tau )\mathrm{d}\tau -u(0)}\end{array}\right).$ Next, we define the derivative of the operator F applied to the element $h=(x,v,y,0,r,z,0)$. We have $F_{\varphi }(\varphi ,\overline{\beta })h=\left(\begin{array}{c}{\displaystyle Q_1-(\mu +\psi )x+wr}\\ {\displaystyle \psi x-Q_2-\mu v}\\ {\displaystyle -y_{\tau }-(\mu +\gamma (\tau ))y}\\ {\displaystyle Q_3-y(0)}\\ {\displaystyle {\int }_0^{\mathrm{\infty }}\gamma (\tau )y(\tau )\mathrm{d}\tau -(\mu +w)r}\\ {\displaystyle -z_{\theta }-\delta (\theta )z}\\ {\displaystyle {\int }_0^{\mathrm{\infty }}\eta (\tau )y(\tau )\mathrm{d}\tau -z(0)}\end{array}\right).$ where $\begin{array}{rl}{Q}_1& =-\overline{\beta }\frac{S}{B+D}{\int }_0^{\mathrm{\infty }}{\beta }_0(\theta )z(\theta )\mathrm{d}\theta -\overline{\beta }\frac{x}{B+D}{\int }_0^{\mathrm{\infty }}{\beta }_0(\theta )u(\theta )\mathrm{d}\theta \\ & +\overline{\beta }\frac{Sb}{(B+D{)}^2}{\int }_0^{\mathrm{\infty }}{\beta }_0(\theta )u(\theta )\mathrm{d}\theta ,\\ Q_2& =\sigma \overline{\beta }\frac{V}{B+D}{\int }_0^{\mathrm{\infty }}{\beta }_0(\theta )z(\theta )\mathrm{d}\theta +\sigma \overline{\beta }\frac{v}{B+D}{\int }_0^{\mathrm{\infty }}{\beta }_0(\theta )u(\theta )\mathrm{d}\theta \\ & -\sigma \frac{Vb}{(B+D{)}^2}{\int }_0^{\mathrm{\infty }}\beta (\theta )u(\theta )\mathrm{d}\theta ,\\ Q_3& =\overline{\beta }\frac{S+\sigma V}{B+D}{\int }_0^{\mathrm{\infty }}{\beta }_0(\theta )z(\theta )\mathrm{d}\theta +\overline{\beta }\frac{x+\sigma v}{B+D}{\int }_0^{\mathrm{\infty }}{\beta }_0(\theta )u(\theta )\mathrm{d}\theta \\ & -\overline{\beta }\frac{(S+\sigma V)b}{(B+D{)}^2}{\int }_0^{\mathrm{\infty }}{\beta }_0(\theta )u(\theta )\mathrm{d}\theta ,\end{array}$ and $b={\int }_0^{\mathrm{\infty }}z(\theta )\mathrm{d}\theta .$ Computing the derivative of F at the DFE, we have: $F_{\varphi }(\varphi ,\overline{\beta })h=\left(\begin{array}{c}{\displaystyle -\overline{\beta }\frac{S^0}{D}{\int }_0^{\mathrm{\infty }}{\beta }_0(\tau )z(\tau )\mathrm{d}\tau -(\mu +\psi )x+wr}\\ {\displaystyle \psi S-\overline{\beta }\frac{\sigma V^0}{D}{\int }_0^{\mathrm{\infty }}{\beta }_0(\tau )z(\tau )\mathrm{d}\tau -\mu V}\\ {\displaystyle -i_{\tau }-(\mu +\gamma (\tau ))i}\\ {\displaystyle \overline{\beta }\frac{S^0+\sigma V^0}{D}{\int }_0^{\mathrm{\infty }}{\beta }_0(\tau )z(\tau )\mathrm{d}\tau -y(0)}\\ {\displaystyle {\int }_0^{\mathrm{\infty }}\gamma (\tau )y(\tau )\mathrm{d}\tau -(\mu +w)r}\\ {\displaystyle -z_{\theta }-\delta (\theta )z}\\ {\displaystyle {\int }_0^{\mathrm{\infty }}\eta (\tau )y(\tau )\mathrm{d}\tau -z(0)}\end{array}\right).$ We define the operator $\mathcal{A}=F_{\varphi }({\varphi }^0,{\overline{\beta }}^{\ast })$ where ${\overline{\beta }}^{\ast }$ is the value of $\overline{\beta }$ that makes $\mathcal{R}(\psi )=1$. We solve the system $\mathcal{A}h=0$. From the differential equations we obtain $y(\tau )=y(0)\pi (\tau ),z(\tau )=z(0){\pi }_{\mathrm{d}}(\tau ).$ We take $y(0)=1$. Then $z(0)=y(0){\int }_0^{\mathrm{\infty }}\eta (\tau )\pi (\tau )\mathrm{d}\tau =\mathrm{\Sigma },$ and $y(0)={\overline{\beta }}^{\ast }\frac{S^0+\sigma V^0}{D}\mathrm{\Sigma }{\int }_0^{\mathrm{\infty }}{\beta }_0(\theta ){\pi }_{\mathrm{d}}(\theta )\mathrm{d}\theta =1.$ This is trivially satisfied. Thus, we obtain $y(\tau )=\pi (\tau ),z(\tau )=\mathrm{\Sigma }{\pi }_{\mathrm{d}}(\tau ).$ From the equation for x we have $x=\frac{-\frac{S^0}{D}\mathrm{\Sigma }{\int }_0^{\mathrm{\infty }}\beta (\tau ){\pi }_{\mathrm{d}}(\tau )\mathrm{d}\tau )+\frac{w}{\mu +w}{\int }_0^{\mathrm{\infty }}\gamma (\tau )\pi (\tau )\mathrm{d}\tau }{\mu +\psi }.$ Further, from the equation for v we have $v=\frac{\psi x-\frac{\sigma V^0}{D}\mathrm{\Sigma }{\int }_0^{\mathrm{\infty }}\beta (\theta ){\pi }_{\mathrm{d}}(\theta )\mathrm{d}\theta }{\mu }.$ Finally, from the equation for r we have $r=\frac{{\int }_0^{\mathrm{\infty }}\gamma (\tau )\pi (\tau )\mathrm{d}\tau }{\mu +w}=\frac{\mathrm{\Gamma }}{\mu +w}.$ Thus, the eigenvector with x and v as above is ${\hat{v}}_0=(x,v,\pi (\tau ),0,\frac{\mathrm{\Gamma }}{\mu +w},\mathrm{\Sigma }{\pi }_d,0).$ To compute the adjoint operator, we take $\xi \in Z$ and consider $\begin{array}{rl}⟨\mathcal{A}h,\xi ⟩& =-\frac{S^0}{D}{\int }_0^{\mathrm{\infty }}\beta (\theta )z(\theta )\mathrm{d}\theta {\xi }_1\\ & -(\mu +\psi )x{\xi }_1+wr{\xi }_1+\psi x{\xi }_2-\frac{\sigma V^0}{D}{\int }_0^{\mathrm{\infty }}\beta (\theta )z(\theta )\mathrm{d}\theta {\xi }_2\\ & -\mu v{\xi }_2-{\int }_0^{\mathrm{\infty }}(y_{\tau }+(\mu +\gamma (\tau )y(\tau )){\xi }_3(\tau )\mathrm{d}\tau \\ & +\frac{S^0+\sigma V^0}{D}{\int }_0^{\mathrm{\infty }}\beta (\theta )z(\theta )\mathrm{d}\theta {\xi }_4-y(0){\xi }_4\\ & +{\int }_0^{\mathrm{\infty }}\gamma (\tau )y(\tau )\mathrm{d}\tau {\xi }_5-(\mu +w)r{\xi }_5-{\int }_0^{\mathrm{\infty }}{z}_{\theta }+\delta (\theta z(\theta ){\xi }_6(\theta )\mathrm{d}\theta \\ & +{\xi }_7{\int }_0^{\mathrm{\infty }}\eta (\tau )y(\tau )\mathrm{d}\tau -{\xi }_{7}z(0).\end{array}$ Integrating the differential equations by parts, we have the following adjoint operator: ${\mathcal{A}}^{\ast }\xi =\left(\begin{array}{c}-(\mu +\psi ){\xi }_1+\psi {\xi }_2\\ -\mu {\xi }_2\\ ({\xi }_3{)}_{\tau }-(\mu +\gamma (\tau )){\xi }_3(\tau )+\gamma (\tau ){\xi }_5+\eta (\tau ){\xi }_6(0)\\ {\xi }_3(0)\\ w{\xi }_1-(\mu +w){\xi }_5\\ ({\xi }_6{)}_{\theta }-\delta (\theta ){\xi }_6(\theta )-{\displaystyle \frac{S^0}{D}\beta (\theta ){\xi }_1-\frac{\sigma V^0}{D}\beta (\theta ){\xi }_2+\frac{S^0+\sigma V^0}{D}\beta (\theta ){\xi }_3(0)}\\ {\xi }_6(0)\end{array}\right).$ To find eigenvector that generates the kernel, we set ${\mathcal{A}}^{\ast }\xi =0$. From the second line, we have ${\xi }_2=0$. From the first line we have ${\xi }_1=0$. From the fifth line, we have ${\xi }_5=0$. Solving the equation for ${\xi }_3$ we have $({\xi }_3{)}_{\tau }-(\mu +\gamma (\tau )){\xi }_3(\tau )+\eta (\tau ){\xi }_6(0)=0.$ Then we have ${\xi }_3(\tau )={\xi }_6(0){\int }_{\tau }^{\mathrm{\infty }}\eta (s){\mathrm{e}}^{-{\int }_{\tau }^s(\mu +\gamma (p))\mathrm{d}p}\mathrm{d}s.$ Further, ${\xi }_3(0)={\xi }_6(0){\int }_0^{\mathrm{\infty }}\eta (s){\mathrm{e}}^{-{\int }_0^s(\mu +\gamma (p))\mathrm{d}p}\mathrm{d}s={\xi }_6(0)\mathrm{\Sigma }.$ Solving the equation for ${\xi }_6(\theta )$ we have $({\xi }_6{)}_{\theta }-\delta (\theta ){\xi }_6(\theta )+\frac{S^0+\sigma V^0}{D}\beta (\theta ){\xi }_6(0)\mathrm{\Sigma }=0.$ Integrating, we obtain: ${\xi }_6(\theta )=\frac{S^0+\sigma V^0}{D}\beta (\theta ){\xi }_6(0)\mathrm{\Sigma }{\int }_{\theta }^{\mathrm{\infty }}\beta (\theta ){\mathrm{e}}^{-{\int }_{\theta }^s\delta (p)\mathrm{d}p}\mathrm{d}s.$ We may take ${\xi }_6(0)=1$. We take ${\xi }_4={\xi }_3(0)$ and ${\xi }_7={\xi }_6(0)$. Then the eigenvector is ${\hat{v}}_0^{\ast }=(0,0,{\xi }_3(\tau ),\mathrm{\Sigma },0,{\xi }_5(\theta ),1).$ To compute a and b we compute the second derivatives: Next, we define the derivative of the operator F applied to the element $h_2=(x_2,v_2,y_2,0,r_2,z_2,0)$. We have $F_{\varphi \varphi }({\varphi }^0,\overline{\beta })[h,h_2]=\left(\begin{array}{c}{\displaystyle Q_1[h,h_2]}\\ {\displaystyle Q_2[h,h_2]}\\ {\displaystyle 0}\\ {\displaystyle Q_3[h,h_2]}\\ {\displaystyle 0}\\ {\displaystyle 0}\\ {\displaystyle 0}\end{array}\right).$ where $\begin{array}{rl}{Q}_1& =-\overline{\beta }\frac{x_{2}D-S{b}_2}{D^2}{\int }_0^{\mathrm{\infty }}{\beta }_0(\theta )z(\theta )\mathrm{d}\theta -\overline{\beta }\frac{x}{D}{\int }_0^{\mathrm{\infty }}{\beta }_0(\theta )z_2(\theta )\mathrm{d}\theta \\ & +\overline{\beta }\frac{Sb}{(D{)}^2}{\int }_0^{\mathrm{\infty }}{\beta }_0(\theta )z_2(\theta )\mathrm{d}\theta ,\\ Q_2& =-\sigma \overline{\beta }\frac{v_{2}D-V{b}_2}{D^2}{\int }_0^{\mathrm{\infty }}{\beta }_0(\theta )z(\theta )\mathrm{d}\theta -\sigma \overline{\beta }\frac{v}{D}{\int }_0^{\mathrm{\infty }}{\beta }_0(\theta )z_2(\theta )\mathrm{d}\theta \\ & +\sigma \frac{Vb}{(D{)}^2}{\int }_0^{\mathrm{\infty }}\beta (\theta )z_2(\theta )\mathrm{d}\theta ,\\ Q_3& =\overline{\beta }\frac{(x_2+\sigma v_2)D-(S+\sigma V)b_2}{D^2}{\int }_0^{\mathrm{\infty }}{\beta }_0(\theta )z(\theta )\mathrm{d}\theta +\overline{\beta }\frac{x+\sigma v}{D}{\int }_0^{\mathrm{\infty }}{\beta }_0(\theta )z_2(\theta )\mathrm{d}\theta \\ & -\overline{\beta }\frac{(S+\sigma V)b}{(D{)}^2}{\int }_0^{\mathrm{\infty }}{\beta }_0(\theta )z_2(\theta )\mathrm{d}\theta .\end{array}$ and $\begin{array}{rl}& b_2={\int }_0^{\mathrm{\infty }}{z}_2(\theta )\mathrm{d}\theta .\\ & F_{\varphi \varphi }({\varphi }^0,\overline{\beta })[h,h]=\left(\begin{array}{c}{\displaystyle -2\overline{\beta }\frac{xD-S^{0}b}{D^2}{\int }_0^{\mathrm{\infty }}{\beta }_0(\theta )z(\theta )\mathrm{d}\theta }\\ {\displaystyle -2\sigma \overline{\beta }\frac{vD-V^{0}b}{D^2}{\int }_0^{\mathrm{\infty }}{\beta }_0(\theta )z(\theta )\mathrm{d}\theta }\\ {\displaystyle 0}\\ {\displaystyle 2\overline{\beta }\frac{(x+\sigma v)D-(S^0+\sigma V^0)b}{D^2}{\int }_0^{\mathrm{\infty }}{\beta }_0(\theta )z(\theta )\mathrm{d}\theta }\\ {\displaystyle 0}\\ {\displaystyle 0}\\ {\displaystyle 0}\end{array}\right).\end{array}$ Thus, (44) $\begin{array}{rl}a& =⟨F_{\varphi \varphi }({\varphi }^0,{\overline{\beta }}^{\ast })[{\hat{v}}_0{\hat{v}}_0],{\hat{v}}_0^{\ast }⟩\\ & =C_1((1+\frac{\sigma \psi }{\mu })\frac{-\frac{S^0}{D}\mathrm{\Sigma }{\int }_0^{\mathrm{\infty }}\beta (\tau ){\pi }_{\mathrm{d}}(\tau )\mathrm{d}\tau )+\frac{w}{\mu +w}{\int }_0^{\mathrm{\infty }}\gamma (\tau )\pi (\tau )\mathrm{d}\tau }{\mu +\psi }\\ & +\sigma \frac{-\frac{\sigma V^0}{D}\mathrm{\Sigma }{\int }_0^{\mathrm{\infty }}\beta (\theta ){\pi }_{\mathrm{d}}(\theta )\mathrm{d}\theta }{\mu }-\frac{(S^0+\sigma V^0)\mathrm{\Sigma }\mathrm{\Pi }}{D}).\end{array}$(44) where $C_1=2\frac{{\mathrm{\Sigma }}^2}{D}{\int }_0^{\mathrm{\infty }}{\beta }^{\ast }(\theta ){\pi }_{\mathrm{d}}(\theta )\mathrm{d}\theta .$ Denote by $\mathcal{B}={\int }_0^{\mathrm{\infty }}\beta (\tau ){\pi }_{\mathrm{d}}(\tau )\mathrm{d}\tau .$ and $\mathrm{\Gamma }={\int }_0^{\mathrm{\infty }}\gamma (\tau )\pi (\tau )\mathrm{d}\tau .$ Then, if $b:=⟨F_{\varphi \overline{\beta }}({\varphi }^0,{\overline{\beta }}^{\ast }){\hat{v}}_0,{\hat{v}}_0^{\ast }⟩$ is positive, the condition for backward bifurcation becomes (45) $Q\frac{w}{(\mu +w)(\mu +\psi )}\mathrm{\Gamma }>Q\frac{S^0}{D(\mu +\psi )}\mathrm{\Sigma }\mathcal{B}+\frac{{\sigma }^2{V}^0}{D\mu }\mathrm{\Sigma }\mathcal{B}+\frac{(S^0+\sigma V^0)\mathrm{\Sigma }\mathrm{\Pi }}{D},$(45) where $Q=(1+\frac{\sigma \psi }{\mu }).$ This is the same condition obtained in [10] as a condition for backward bifurcation of this model. Finally we check the positivity of b. First we compute the derivative of $F_{\varphi }$ with respect to $\overline{\beta }$. We have $F_{\varphi \overline{\beta }}({\varphi }^0,\overline{\beta })h=\left(\begin{array}{c}{\displaystyle -\frac{S^0}{D}{\int }_0^{\mathrm{\infty }}{\beta }_0(\theta )z(\theta )\mathrm{d}\theta }\\ {\displaystyle -\sigma \frac{V^0}{D}{\int }_0^{\mathrm{\infty }}{\beta }_0(\theta )z(\theta )\mathrm{d}\theta }\\ {\displaystyle 0,}\\ {\displaystyle \frac{S^0+\sigma V^0}{D}{\int }_0^{\mathrm{\infty }}{\beta }_0(\theta )z(\theta )\mathrm{d}\theta }\\ {\displaystyle 0}\\ {\displaystyle 0}\\ {\displaystyle 0}\end{array}\right).$ Then we have, (46) $b=⟨F_{\varphi \overline{\beta }}({\varphi }^0,{\overline{\beta }}^{\ast }){\hat{v}}_0,{\hat{v}}_0^{\ast }⟩=\frac{S^0+\sigma V^0}{D}{\mathrm{\Sigma }}^2{\int }_0^{\mathrm{\infty }}{\beta }_0(\theta ){\pi }_{\mathrm{d}}(\theta )\mathrm{d}\theta >0.$(46) This concludes the example.

Acknowledgements

The authors would like to thank the editors and two anonymous reviewers for their helpful comments, HI was supported by the japan society for promotion of Science Grant-in-Aid for Scientific Research (c) (No. 19K03614).

Disclosure statement

No potential conflict of interest was reported by the author(s).

Notes

1 Another possible interpretation for (19(19) $\begin{array}{rl}{s}_t+s_a& =-\kappa (a)s(a,t){\int }_0^{\omega }\beta (\sigma )i(\sigma ,t)\mathrm{d}\sigma +\frac{\gamma (a)i(a,t)}{A+{\int }_0^{\omega }i(a,t)\mathrm{d}a}\\ i_t+i_a& =\kappa (a)s(a,t){\int }_0^{\omega }\beta (\sigma )i(\sigma ,t)\mathrm{d}\sigma -\frac{\gamma (a)i(a,t)}{A+{\int }_0^{\omega }i(a,t)\mathrm{d}a}\end{array}$(19) ) is that s and i are prevalencies, B = 1 and the survival rate is Type 1, that is, all deaths occur at $a=\omega$.