![MathJax Logo](/templates/jsp/_style2/_tandf/pb2/images/math-jax.gif)
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.
AMS Subject Classification:
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 [Citation10] 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 [Citation7].
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 [Citation2]. 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 [Citation2], 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 -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 [Citation10], 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 [Citation10] 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 [Citation2] designed for PDEs. The proof is based on the Lyapunov–Schmidt theory.
Let X and Z be Banach spaces, and and
is a parameter. We consider the following abstract differential equation
(1)
(1) Without loss of generality, we assume that 0 is an equilibrium so
Theorem 2.1
Assume:
is the linearization around the equilibrium 0, evaluated at a critical value of parameter
such that
is a closed operator with a simple isolated eigenvalue zero and remaining eigenvalues having negative real part. Let
be the unique (up to a constant) positive solution of
.
for some neighbourhood
of 0 and interval
containing
.
Assume
is the dual of Z and
is the pairing between Z and
. Assume
is the unique (up to a constant) positive vector satisfying
that is
, where
is the adjoint of
, and
.
Assume
where
is the second derivative of F with respect to x and α.
Then, the direction of the bifurcation is determined by the numbers
(2)
(2) where
is the second derivative of F with respect to x applied to the functions
and
. If
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 in a neighbourhood of
. Let
Since
is closed and has a simple isolated eigenvalue, it follows from [Citation4, p. 284] that
is a Fredholm operator of index zero, that is
satisfies:
;
;
is closed in Z
where
and
Because,
, we have
where
and
. Let
. Then following [Citation3] we have:
Similarly, if
, then
where
is the topological complement of
in X.
It follows that there exist projections and
such that
and
. Consequently we can write
in the form
where
. That implies that there exists a constant c such that
. In addition
. Using the projection Q we can decompose
into two equations
(3)
(3)
Let
. Then the first equation in (Equation3
(3)
(3) ) is written as
. In particular, we have
. Furthermore,
The restriction
is a one-to-one mapping between
and
and it is therefore invertible. Hence
is invertible. As a result we can apply the Implicit Function Theorem. From the Implicit Function Theorem, it follows that there exists a neighbourhood
with
and an interval
with
and a
map
with
such that
and
Furthermore, we define a map
by
Then the second equation in (Equation3
(3)
(3) ) can be written as
(4)
(4) which is called the bifurcation equation.
The Lyapunov–Schmidt reduction says that there exists a neighbourhood U of such that each solution of
corresponds one-to-one to some solution of
. Furthermore, it is shown in [Citation8, Corollary I.2.4] that
(5)
(5) From
for all
, we have that
Furthermore,
Consequently,
This corresponds to the trivial solution. This last equality means that we can write
in the following way [Citation8]:
Alternatively,
For
, we set
,
and we obtain the nontrivial solutions for
by solving
The assumptions on the smoothness of F give that
. Furthermore, by the properties (Equation5
(5)
(5) ) we have
. We would like to apply the Implicit Function Theorem and express
. To do that, we need
. Computation similar to [Citation8, p. 19] give that
(6)
(6) Setting
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)
(7) and
because of assumptions A3 and A4. In particular,
and
. The Implicit Function Theorem now gives the existence of
such that
with
and
Then
for
. Consequently,
where
(8)
(8)
This is the essence of the Crandall-Rabinowitz Theorem [Citation8]. Next we claim that the tangent vector to the nontrivial solution curve (Equation8
(8)
(8) ) at the bifurcation point
is given by
Differentiating
in (Equation8
(8)
(8) ) with respect to s gives
(9)
(9) The bifurcation at the critical value
is backward if and only if
. To obtain formulas to compute
, we recall that
and therefore
for all
. Thus,
(10)
(10) where
is given in (Equation7
(7)
(7) ). Furthermore,
(11)
(11) Using the dual pairing and from Equation (Equation10
(10)
(10) ) we have
Notice that
(12)
(12) These equalities hold because
for all
. Hence, we have
If we assume that b>0, the bifurcation curve bifurcates backward if
, that is, if a>0. If
, that is a<0, then the bifurcation is forwards.
Remark 2.1
The condition b>0 means that increases with the parameter α. If b<0, that means that
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 [Citation1], let us show the Factorization Theorem. Let be the bifurcating solutions as
(13)
(13) and
Define
Then
. The dominant eigenvalue
and corresponding eigenvector are introduced as follows:
(14)
(14) where
Define the dual eigenvector
as
where they are normalized as follows:
Especially, we write as
.
Differentiating (Equation13(13)
(13) ) with respect to s, we have
(15)
(15) From (Equation15
(15)
(15) ), taking the duality pairing with
, it follows that
(16)
(16) Then we arrive at an expression
(17)
(17) Using the Taylor expansion around s = 0, we have
where we used
. Then it holds that
(18)
(18) Therefore, if we assume that b>0, for small s>0, the backward bifurcation (
) implies the bifurcated solution is unstable (
), while the forward bifurcation (
) implies the bifurcated solution is locally stable (
) as long as we could assume that
is the dominant eigenvalue. In many concrete examples, we could show the dominance of
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)
(19)
where a is chronological age, t is time, and ω is the maximal age, assumed finite. Furthermore,
and
are the age-density functions of susceptibles and invectives, respectively,
is the susceptibility of susceptible individuals,
is the age-specific infectivity of infective individuals,
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 neglectedFootnote1.
We assume that
Although model (Equation19
(19)
(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)
(20) where B is the birth rate, assumed constant. Furthermore, we assume that
for all
. The equation of the total population size is
. Hence, the total population is constant
. Let,
. Thus, we reduce the system to a unique equation:
(21)
(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)
(22) This equation clearly has the disease-free equilibrium
. Let us introduce a bifurcation parameter
as
where a normalized function
will be determined below. Thus, we have a problem
where
where A is the differential operator given by
with domain
where
is the set of absolutely continuous functions on the interval
. The operator F is defined from
to
.
To compute we consider a parameter s and a function
with properties to be specified later and we compose
. Taking
and setting s = 0 we obtain
(23)
(23) Setting i = 0, we have
This is not exactly the operator
since the derivative is evaluated at an arbitrary value of the parameter
.
To see the stability of the disease-free equilibrium, we look for eigenvalues of the operator . This leads to the linearized equation of the disease-free equilibrium:
(24)
(24) where
Solving the above differential equation we have
Substituting
in
and cancelling since we are looking for non-zero h, we obtain the following characteristic equation:
From here we can define the basic reproduction number
The normalized function
is assumed to satisfy
(25)
(25) Then we have
and
is a simple isolated eigenvalue of the characteristic equation when
. It is not hard to see that all remaining eigenvalues have negative real part.
We define the operator :
which has the same domain as
. Looking for
, we need to solve
Following a procedure similar to obtaining
, we obtain
We notice that from (Equation25
(25)
(25) ) we have
To obtain
we need to formally find the adjoint operator, and obtain its kernel. Since
, the dual space is
. Let
be the adjoint with
For
, we have after integration by parts
(26)
(26) Hence,
with
. To find
we need to determine the kernel of the adjoint operator. For that purpose, we have to solve the problem:
(27)
(27) where
Solving this system as before, we obtain
It is not hard to see from integration by parts that
Next, we compute a and b. First we compute b. We have,
Hence, we have
(28)
(28) Next, we compute a. First, we need to compute
. To do that we consider
(29)
(29) Differentiating with respect to s and setting s = 0 and i = 0, we obtain
(30)
(30) We notice that this form is symmetric with respect to h and
. Thus, we have
If we recall
, multiplying with
and integrating, we obtain
(31)
(31) Hence, the following is a necessary and sufficient condition for backward bifurcation (a>0)
(32)
(32)
4. Example 2: infection-age structured HIV/AIDS epidemic model
The example in this section is an HIV/AIDS epidemic model developed in [Citation11]. 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 [Citation5,Citation6], 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)
(33)
where
is the number of susceptibles,
is the infection-age density of infectives and
is the total size of sexually active population, which is not constant because of the disease-induced death rate
. We assume that
is normalized as
(34)
(34) where
denotes the survival rate of infecteds, so we can use
as a bifurcation parameter, and
is no other than the basic reproduction number.
According to [Citation9], in order to take into account the boundary condition, we define the extended state space associated with
. We also define the subspace by
, where
. Then
can be identified with the usual state space
. Let
be the disease-free steady state, where
. Define a linear operator
by
for
with
.
Let . Then the basic problem (Equation36
(36)
(36) ) can be formulated as a semilinear Cauchy problem in X:
(35)
(35)
where for
,
and
are defined by
and
is a positive functional of
defined by
Moreover, define a linear operator
and a nonlinear operator
by
with
. Then
is not dense in X, but A satisfies the Hille-Yosida estimate. Then (Equation35
(35)
(35) ) can be formulated as a well-posed abstract Cauchy problem in X as follows:
(36)
(36) where
.
In order to examine the bifurcation at the disease-free steady state of (Equation39(39)
(39) ), let
for
. Then
if
is the disease-free steady state. Then the linearized operator
acting on X is calculated as follows:
(37)
(37) where its domain is given by
Let us consider the eigenvalues of
. Solving the differential equation from
for
, we have
Inserting the above expression into the third equation in
, we have
Then we obtain the characteristic equation
(38)
(38) which has a simple eigenvalue of zero when
. To find the eigenvector
associated with eigenvalue zero, we need to solve
. Then we can choose
and the boundary condition (the third equation of
) is trivially satisfied at
. From the first equation of
, we have
and
To find the adjoint operator
, let
, we multiply
with the vector
. Then we have
Observe that
where we assume that
. Hence, we have
which should hold for
and
.
Thus we choose the domain of as
Then the adjoint operator
is defined on
as follows:
To find
, we solve
with
. We have
. Solving the differential equation in
, we have
(39)
(39) Hence, we have the eigenfunction
, where
.
Computing the second derivative for F, we have
Hence, we obtain
(40)
(40) On the other hand,
can be computed as follows:
Thus, we have
(41)
(41) Hence, it follows from (Equation41
(41)
(41) ) and (Equation40
(40)
(40) ) that
, so the bifurcation at
is a forward one, which is a conclusion that we expected.
Remark 4.1
In this example, just as in the ODE case [Citation2], we see that components of the eigenvector , 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 [Citation10]. 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 the density of the infectives and by
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)
(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 [Citation10]; β and δ are functions of θ, while γ and η are structured by τ. B in this model denotes the concentration of bacteria in the environment
This system has a disease-free equilibrium
where
We define the probability of survival of infected individual in the infectious class and the probability of survival of the bacteria in the environment:
The control reproduction number is given by (derivation can be found in [Citation10]):
(43)
(43) To have a bifurcation parameter, we decompose
as
where
. We define
and
. For an element
with
we define the non-linear operator
Next, we define the derivative of the operator F applied to the element
. We have
where
and
Computing the derivative of F at the DFE, we have:
We define the operator
where
is the value of
that makes
. We solve the system
. From the differential equations we obtain
We take
. Then
and
This is trivially satisfied. Thus, we obtain
From the equation for x we have
Further, from the equation for v we have
Finally, from the equation for r we have
Thus, the eigenvector with x and v as above is
To compute the adjoint operator, we take
and consider
Integrating the differential equations by parts, we have the following adjoint operator:
To find eigenvector that generates the kernel, we set
. From the second line, we have
. From the first line we have
. From the fifth line, we have
. Solving the equation for
we have
Then we have
Further,
Solving the equation for
we have
Integrating, we obtain:
We may take
. We take
and
. Then the eigenvector is
To compute a and b we compute the second derivatives: Next, we define the derivative of the operator F applied to the element
. We have
where
and
Thus,
(44)
(44) where
Denote by
and
Then, if
is positive, the condition for backward bifurcation becomes
(45)
(45) where
This is the same condition obtained in [Citation10] as a condition for backward bifurcation of this model. Finally we check the positivity of b. First we compute the derivative of
with respect to
. We have
Then we have,
(46)
(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 (Equation19(19)
(19) ) is that s and i are prevalencies, B = 1 and the survival rate is Type 1, that is, all deaths occur at
.
References
- N.F. Britton, Reaction-Diffusion Equations and Their Applications to Biology, Academic Press, London, 1986.
- C. Castillo-Chavez and B Song, Dynamical models of tuberculosis and their applications, Math Biosciences Eng 1(2) (2004), pp. 361–404. doi: 10.3934/mbe.2004.1.361
- S. Guo and J. Wu, Bifurcation Theory of Functional Differential Equations, Springer, New York, 2013.
- M. Haragus and G. Iooss, Local Bifurcations, Center Manifolds and Normal Forms in Infinite-Dimensional Dynamical Systems, Springer, New York, 2011.
- H. Inaba, Backward bifurcation in a HIV/AIDS epidemic model with age structure II: the Case of General Transmission Rate, RIMS Kokyuroku Vol. 1337, Mathematical Economics, Resarch Institute for Mathematical Sciences, Kyoto University, Kyoto, 2003, p. 103–111.
- H. Inaba, Endemic threshold results for age-duration-structured population model for HIV infection, Math Biosci 201 (2006), pp. 15–47. doi: 10.1016/j.mbs.2005.12.017
- H. Inaba, Age-Structured Population Dynamics in Demography and Epidemiology, Springer, Singapore, 2017.
- H. Kielhöfer, Bifurcation Theory: An Introduction with Applications to Partial Differential Equations, 2nd ed., Springer, New York, 2012.
- P. Magal, C.C. McCluskey, and G.F. Webb, Lyapunov functional and global asymptotic stability for an infection-age model, Applicable Anal 89(7) (July 2010), pp. 1109–1140. doi: 10.1080/00036810903208122
- M. Martcheva, Methods for deriving necessary and sufficient conditions for backward bifurcation, JBD 13(1) (2019), pp. 538–566. doi: 10.1080/17513758.2019.1647359
- H. Thieme and C. Castillo-Chavez, How may infection-age-dependent infectivity affect the dynamics of HIV/AIDS?, SIAM J Appl Math 53(5) (1993), pp. 1447–1479. doi: 10.1137/0153068