A Partial Differential Equation Model with Age-Structure and Nonlinear Recidivism: Conditions for a Backward Bifurcation and a General Numerical Implementation
Fabio Sanchez Juan G. Calvo Esteban Segura Zhilan Feng
math.DS q-bio.PE 91D10 (Primary), 35Q91 (Secondary)


A Partial Differential Equation Model with Age-Structure and Nonlinear Recidivism: Conditions for a Backward Bifurcation and a General Numerical Implementation
    [
    December 30, 2023
===============================================================================================================================================================










We formulate an age-structured three-staged nonlinear partial differential equation model that features nonlinear recidivism to the infected (infectious) class from the temporarily recovered class. Equilibria are computed, as well as local and global stability of the infection-free equilibrium. As a result, a backward-bifurcation exists under necessary and sufficient conditions. A generalized numerical framework is established and numerical experiments are explored for two positive solutions to exist in the infectious class.









§ INTRODUCTION
 
Ordinary differential equation models with nonlinear recidivism have been explored in previous studies <cit.>. These models exhibit a phenomena known as a backward bifurcation, a behavior that is strongly correlated to initial conditions, and more specifically, to the initial number of infectious individuals. The implications of a backward bifurcation in epidemic models are crucial in understanding how to develop control mechanisms for the disease. 

The main concept in these models is the study of the threshold quantity ℛ_0, the basic reproductive number, which represents the number of secondary infections caused by a infectious individual when introduced into a mostly susceptible population. Typically, the condition of having ℛ_0<1 is sufficient to have the disease die out in compartmental models. In other words, the infection-free equilibrium is stable. Conversely, the disease prevails and the endemic equilibrium becomes asymptotically stable when ℛ_0>1. However, if there is a backward bifurcation, then ℛ_0<1 is not a sufficient condition for the disease to die out. Furthermore, efforts to control the disease will depend on the density of infectious individuals in the population or having ℛ_0 ≪ 1. 

In ordinary differential equation models, when studying the bifurcation behavior, it is represented by the graph of the infected class versus ℛ_0. This is represented as a function of the transmission rate. This behavior translates into the homologous nonlinear partial differential equation model. The bifurcation in our model is described graphically by a surface plot represented by the infected class at a steady state distribution as a function of age and ℛ_0. The main difference is the age dependence which can lead to the development of better control efforts using the individuals age distribution.

Other studies have included age-structure in their models <cit.>. However, in the model presented here, besides including age-structure, we also include the possibility of recidivism, which leads to a backward bifurcation under some conditions. A strictly theoretical approach was studied in <cit.>. 

The rest of this paper is organized as follows. In <Ref>, we describe the homologous partial differential equation model, based on the nonlinear ordinary differential equation model presented in <cit.>. In <Ref>, we study the infection-free non-uniform steady state distribution, its stability, and we define the basic reproductive number, ℛ_0. We study the endemic non-uniform steady state distribution in <Ref>, where we also analyze the existence of two positive solutions if necessary and sufficient conditions hold, which results in a backward bifurcation. In <Ref>, a general numerical framework is presented, including an example based on the parameters used in <cit.>, as well as an example with age-dependent parameters. Finally, in <Ref> we give some conclusions of the model.





§ A MODEL WITH NONLINEAR RECIDIVISM
 
We construct an age-dependent nonlinear partial differential equation model where S(t,a), I(t,a) and R(t,a) represent susceptible, infected and temporarily recovered individuals at time t and age a, respectively. Susceptible individuals can become infected by having contact with an infectious individual at rate β(a) B(t,a), where B(t,a) is the probability that a person contacted by a susceptible of age a is infectious at time t and β(a) is the transmission rate at age a; they can also exit the system at rate μ(a), where μ(a) is the exit rate of the system at age a. Infected individuals receive treatment at rate ϕ(a), can recover at rate γ(a) or exit the system at rate μ(a). Analogously, temporarily recovered individuals can become infected (subsequent infections) by having contact with an infectious individual at rate ρ(a) B(t,a) and can exit the system at rate μ(a). In this case, the class R(t,a) acts as another susceptible" class.

The age-dependent mixing contact structure is modeled via the mixing density p(t,a,a'), which gives the proportion of contact that individuals of age a have with individuals of age a', given that they had contact with somebody at time t. We will restrict ourselves to the case of proportional mixing; this is, 
    p(t,a,a') ≡ p(t,a) =  c(a)n(t,a)/∫_0^∞ c(a) n(t,a) da,
 with c(a) the age-specific per capita contact/activity rate; see, e.g., <cit.>. We then define the force of infection 
    B(t) := ∫_0^∞I(t,a)/n(t,a) p(t,a)da,
 where n(t,a):= S(t,a)+I(t,a)+R(t,a) is the total population.

The model we just described is given by the system of partial differential equations

    ( ∂/∂ t + ∂/∂ a) S(t,a)    =  - β(a) S(t,a) B(t) -μ(a)S(t,a),
    ( ∂/∂ t + ∂/∂ a) I(t,a)    = β(a) S(t,a) B(t) -(μ(a) + ϕ(a) + γ(a)) I(t,a)+ρ(a) R(t,a) B(t),
    ( ∂/∂ t + ∂/∂ a) R(t,a)    =  (ϕ(a) + γ(a)) I(t,a)- ρ(a) R(t,a) B(t) - μ(a) R(t,a),
 
for t>0, a>0, and boundary conditions for t=0 and a=0 given by 

    S(t,0)   =Λ,     I(t,0)   =0,     R(t,0)   =0,
    
    S(0,a)   =S_0(a),     I(0,a)   =I_0(a),     R(0,a)   =R_0(a),

where Λ denotes the birth rate (assumed constant), and S_0(a), I_0(a), R_0(a) are given initial conditions. For simplicity, we will assume that β(a), μ(a), ϕ(a), γ(a), ρ(a) are continuous functions.

The total population n(t,a) satisfies the initial value problem

    ( ∂/∂ t +∂/∂ a) n(t,a)    = -μ(a)n(t,a),
    
    n(t,0)    = Λ,
    
    n(0,t)    = S_0(a)+I_0(a)+R_0(a) =: n_0(a),

with t>0 and a>0, where Λ represents the constant rate at which newborns enter the population. When solving by using the method of characteristics we obtain 

    n(t,a) = {[ n_0(a-t) ℱ(a)ℱ(a-t)            if  t<a,;               Λℱ(a)            if  t>a, ].

where 
    ℱ(a)=e^ -∫_0^a μ(h)dh.
 

We note that n(t,a) is continuous if and only if the compatibility condition n_0(0)=Λ is satisfied. At demographic steady state, 
    lim_t→∞ n(t,a)=Λℱ(a),
 and therefore 
    p_∞(a) := lim_t→∞ p(t,a)=c(a) ℱ(a)/∫_0^∞ c(a) ℱ(a) da.
  
We rescale variables by making the substitution

    s(t,a) := S(t,a)n(t,a),   i(t,a) := I(t,a)n(t,a),   r(t,a) := R(t,a)n(t,a).

Thus, 

    B(t) = ∫_0^∞ i(t,a) p(t,a) da
 
and system (<ref>) becomes
 

    ( ∂/∂ t + ∂/∂ a) s(t,a)    = - β(a) s(t,a) B(t),
    
    ( ∂/∂ t + ∂/∂ a) i(t,a)    = β(a) s(t,a) B(t) - (ϕ(a)+γ(a)) i(t,a)+ρ(a) r(t,a) B(t),
    
    ( ∂/∂ t + ∂/∂ a) r(t,a)    =  (ϕ(a)+γ(a)) i(t,a)- ρ(a) r(t,a) B(t),

along with boundary conditions

    s(t,0)   =1,     i(t,0)   =0,     r(t,0)   =0,
    
    s(0,a)   =s_0(a),     i(0,a)   =i_0(a),     r(0,a)   =r_0(a),


where s_0(a)+i_0(a)+r_0(a) = 1. 



§ INFECTION-FREE NON-UNIFORM STEADY STATE DISTRIBUTION AND BASIC REPRODUCTIVE NUMBER
 
In this section we explore the infection-free steady state distribution and study its local and global stability. We also compute the basic reproductive number ℛ_0, and the critical value ℛ_C that will determine global stability to the infection-free steady state.

For the infection-free non-uniform steady state distribution, system (<ref>) simplifies to
 

    ds^*(a)da   = 0,
    
    di^*(a)da   = 0 ,
    
    dr^*(a)da   =  0,

with a > 0 and along with boundary conditions

    s^*(0)=1,    i^*(0)=0,   r^*(0)=0.



Thus, system (<ref>) supports the infection-free non-uniform 

distribution given by

    s^*(a) = 1,   i^*(a) = 0,   r^*(a) =  0.


The local stability of (<ref>) is studied using the perturbations

    s(t,a)    = ŝ(a)e^λ t + s^*(a) = ŝ(a)e^λ t + 1,
    
    i(t,a)    = î(a)e^λ t,
    
    r(t,a)    = r̂(a)e^λ t, 
    
    B(t)     = B_0 e^λ t,

where

    B_0 = ∫_0^∞î(a)p_∞ (a)da.

We then have that

    ( ∂/∂ t + ∂/∂ a) i(t,a) = λî(a)e^λ t + d î(a)da e^λ t,

and from (<ref>), we obtain 

    λî(a) + d î(a)da   = β(a) s^*(a) B_0 -(ϕ(a)+γ(a)) î(a) + ρ(a) r^*(a) B_0 
       = β(a) B_0 - (ϕ(a)+γ(a)) î(a).

Thus, î(a) solves the first-order linear ordinary differential equation

    d î(a)da + ( λ +  ϕ(a)+γ(a)) î(a) = β(a) B_0  with  î (0)=0,









which solution is given by

    î (a) 
    
       =  B_0 ∫_0^a  β(h) e^-λ (a-h) -∫_h^a (ϕ(k) +γ(k) )dk dh.
 
Multiplying by p_∞(a) and integrating, we get from (<ref>) that

    B_0 = B_0 ∫_0^∞∫_0^a  p_∞ (a) β(h) e^-λ (a-h) - ∫_h^a (ϕ(k) + γ(k))dk  dh da.

Therefore, for B_0 ≠ 0,

    1 = ∫_0^∞∫_0^a  p_∞ (a) β(h) e^-λ (a-h) - ∫_h^a (ϕ(k) +γ(k))dk  dh da.

Then we have the corresponding Euler-Lotka equation (see, e.g., <cit.>)

    G(λ):= ∫_0^∞∫_0^a  p_∞ (a) β(h) e^-λ (a-h) - ∫_h^a (ϕ(k)+γ(k))dk  dh da =1.

We then define the control adjusted reproductive number, ℛ_0, by

    ℛ_0 :=G(0) = ∫_0^∞∫_0^a  p_∞ (a) β(h) e^- ∫_h^a (ϕ(k)+γ(k))dk  dh da.

We also consider

    ℛ_C := ∫_0^∞∫_0^a  p_∞ (a) β(h)  dh da.


We then establish the following theorem:


If ℛ_0 < 1, the infection-free non-uniform steady state distribution (<ref>) is locally asymptotically stable and unstable if ℛ_0>1. 

 
Observe that G^' (λ)<0 and G^''(λ) >0 for all λ, and that lim_λ→∞ G(λ) = 0 and lim_λ→ -∞ G(λ) = ∞. Then G(λ) is a continuous, strictly decreasing, convex function of λ that takes on all positive values. Thus (<ref>) has exactly one real root and it is a dominant root (see <cit.>).
Therefore we have that (<ref>) has a unique negative real solution if and only if ℛ_0<1 and it has a unique positive real solution if and only if ℛ_0>1.
Then, the fact that the unique real root is dominant guarantees the local asymptotic stability of the infection-free non-uniform steady state distribution, i.e., it is locally asymptotically stable if ℛ_0 <1 and unstable if ℛ_0 >1.


We also have the following result about global stability:

Assume that ℛ_C  < 1. Then, the infection-free steady state distribution of System (<ref>) is globally asymptotically stable.




Define q(t,a):= i(t,a) + r(t,a). From (<ref>), q satisfies the equation 

    ( ∂/∂ t + ∂/∂ a) q(t,a) = β(a) s(t,a) B(t).

We use then the method of characteristics to obtain that for t>a

    s(t,a)    = e^ -∫_0^a β(h) B(h+t-a) dh,
    
    q(t,a)    = ∫_0^a β(h) B(h+t-a) s(h+t-a,h) dh.


Define Q(a) := lim sup_t→∞ q(t,a), I(a) := lim sup_t→∞ I(t,a) and C := ∫_0^∞ p_∞(a) I(a) da. Thus, from (<ref>) we deduce that

    I(a)≤ Q(a)≤ C∫_0^a β(h) e^-C∫_0^h β(τ) dτ dh = 1-e^-C∫_0^a β(τ) dτ.

Multiplying by p_∞(a) and integrating with respect to a, we then have that

    C ≤∫_0^a (1-e^-C∫_0^a β(τ) dτ) p_∞(a) da.

Suppose that C> 0. We can write then

    1 ≤1C∫_0^a (1-e^-C∫_0^a β(τ) dτ) p_∞(a) da =: F(C).


It is straightforward to verify that 
    lim_C→ 0 F(C) =  ℛ_C,
 
and that 

    F'(C) = 1C^2∫_0^∞ p_∞(a) e^-C∫_0^a β(τ) dτ(1+ C ∫_0^a β(τ) dτ -e^ C∫_0^a β(τ) dτ) da < 0

since e^x > 1+x for x≠ 0. Therefore, 
    F(C)≤lim_C→ 0 F(C) = ℛ_C < 1,

which is a contradiction. We then conclude that C=0. Since i(t,a)≥ 0 and r(t,a) ≥ 0, we conclude that lim_t→∞ i(t,a) = lim_t→∞ r(t,a) = 0 and lim_t→∞ s(t,a) = 1.

 


§ ENDEMIC NON-UNIFORM STEADY STATE DISTRIBUTION
 

In this section we analyze necessary conditions for the existence of non-uniform steady state distributions. We first prove that this is the case if ℛ_0 > 1. We then compute explicit conditions in order to observe a backward bifurcation, for the particular case of constant coefficients μ, β,ϕ, γ, ρ, c.



 §.§ Existence of at least one endemic non-uniform steady state

In this section we establish conditions for which at least one positive solution exists. For our model, the existence of solutions mostly rely under the threshold condition ℛ_0.

 
There exists an endemic non-uniform steady state of system (<ref>) when ℛ_0 > 1.


A non-uniform steady state age-distribution is a solution of the nonlinear system
 

    ds^*(a)da   = - B^* β(a)  s^*(a),
    
    di^*(a)da   = B^* β(a) s^*(a) -  (ϕ(a)+γ(a)) i^*(a)+ B^*ρ(a) r^*(a) ,
    
    dr^*(a)da   =  (ϕ(a)+γ(a)) i^*(a) - B^*ρ(a) r^*(a) ,

with a > 0, along with boundary conditions

    s^*(0)=1,    i^*(0)=0,   r^*(0)=0,

where 

    B^* := ∫_0^∞ i^*(a) p_∞ (a) da.



Consider the linear system of equations with parameter B given by 
 

    ds_B^*(a)da   =  - B β(a)  s_B^*(a) ,
    
    di_B^*(a)da   = B β(a) s_B^*(a) -( ϕ(a)+γ(a)) i_B^*(a)+ Bρ(a) r_B^*(a) ,
    
    dr_B^*(a)da   =  (ϕ(a)+γ(a)) i_B^*(a) - Bρ(a) r_B^*(a) ,

along with boundary conditions

    s_B^*(0)=1,    i_B^*(0)=0,   r_B^*(0)=0.



Given the solution (s_B^*(a), i_B^*(a), r_B^*(a)) of system (<ref>), let

    H(B) := ∫_0	^∞ i_B^*(a) p_∞ (a) da.

It is clear that (s_B^*(a), i_B^*(a), r_B^*(a)) satisfies system (<ref>) if and only if B is a fixed point of H; i.e., H(B) = B. Furthermore, if B=0 then 

i_B^*(a) = 0 and H(0)=0.
Thus, in order to guarantee existence of at least one non-trivial solution to (<ref>), it is just necessary to prove that H(B) has a fixed point B>0. Consider the auxiliary function G(B) = H(B)/B for B≠ 0. Clearly G(B) is continuous in [0,1] by defining G(0) = H'(0). Since H(1)<1, we have that G(1)<1. Thus, we just need to prove that G(0)>1. If this is the case, then there exists B∈ (0,1) such that G(B)=1.

We then compute explicitly G(B). First, from (<ref>), we compute that

    s_B^*(a) = e^-B∫_0^aβ(h)dh .

Let X_B(a) = [i_B^*(a), r_B^*(a)]^T. We can write (<ref>) and (<ref>) as the linear system

    dda	X_B(a) = M(a) X_B(a)  + f_B(a),

where 

    M(a)=
      [ -(ϕ(a)+γ(a))        Bρ(a);    ϕ(a)+γ(a)       -Bρ(a);              ],  f_B(a)=
      [ Bβ(a) s_B^*(a);              0;                ].

We then observe that 

    X^(1)_B(a) :=
    [
    [  1; -1 ]] e^-∫_0^a (ϕ(h)+γ(h)+ Bρ(h))dh

is a solution of the homogeneous system dda	X_B(a) = M(a) X_B(a) with initial condition X^(1)_B(0) = [1, -1]^T. By Abel's formula, we compute a second linear independent solution X^(2)_B(a) for the homogeneous system. After straightforward computations, a linear independent solution with initial condition X^(2)_B(0) = [0,1]^T is given by

    X^(2)_B(a) :=
    [
    [ 1 - e^-∫_0^a (ϕ(h)+γ(h)+ Bρ(h))dh - ∫_0^a (ϕ(τ)+γ(τ)) e^-∫_τ^a  (ϕ(h)+γ(h)+ Bρ(h))dh dτ;       e^-∫_0^a (ϕ(h)+γ(h)+ Bρ(h))dh+∫_0^a (ϕ(τ)+γ(τ)) e^-∫_τ^a  (ϕ(h)+γ(h)+ Bρ(h))dh dτ ]].

Therefore, the fundamental matrix Φ(a) for system (<ref>) is given by 

    Φ(a) = [X^(1)_B(a)  X^(2)_B(a) ].

By variation of parameters, the solution of system (<ref>) is given by 
    X_B(a)  = Φ(a) C(a),
 where C(a) satisfies the equation 

    C'(a) = Φ^-1(a) f_B(a).

Since 

    [
    [ 0; 0 ]] = X_B(0)  = Φ(0) C(0) = 
    [
    [  1  0; -1  1 ]] C(0),

we deduce that C(0) = [0,0]^T. Then, by integrating (<ref>) we obtain that 

    C(a) = [ [ C_1(a); C_2(a) ]] = [
    [ B∫_0^a   β(τ) s_B^*(τ) ( 1+∫_0^τ (ϕ(t)+γ(t)) e^∫_0^t (ϕ(h)+γ(h)+Bρ(h))dh dt )  dτ;                                                          B∫_0^a  β(τ) s_B^*(τ) dτ ]].

After some simplifications, we conclude that r_B^*(a) and i_B^*(a) are given explicitly by

    r_B^*(a) =    ∫_0^a (ϕ(τ) +γ(τ))( 1- s_B^*(τ )) e^-∫_τ^a ( ϕ(h)+γ(h) + B ρ(h)) dh dτ, 
    
    i_B^*(a) =    1 - s_B^*(a) - r_B^*(a).

Since lim_B→ 01-s_B(a)B =  ∫_0^a β (h) dh,
we obtain that 
    G(0) = lim_B→ 0G(B) =  ∫_0^∞∫_0^a β(τ) e^-∫_τ^a ( ϕ(h)+γ(h) )dh p_∞ (a) dτ da = ℛ_0.

By hypothesis ℛ_0 >1. Therefore, there exists B∈ (0,1) such that G(B) = 1, completing our proof.


 In the case of constant coefficients, we can compute that 
    ℛ_0 = βμ+ϕ+γ,   ℛ_C = βμ;
 see (<ref>) below. 
Thus, a large value for β implies the existence of at least one endemic non-uniform steady state, while a small value for β guarantees that the disease tends to zero. In the next section we explore necessary and sufficient conditions that lead to the existence of a backward bifurcation.









 §.§ Conditions for a backward bifurcation

In this section we explore conditions for the existence of two positive steady states when ℛ_0<1. For simplifity, it is assumed that μ, β, ρ, ϕ, γ and c are constant in the characterization of the backward bifurcation. 

 We assume constant parameters as an illustrative example, however, initial conditions remain age-dependant.




By solving (<ref>) we obtain that (see (<ref>))

    s_B^*(a) := e^-aB β .

Moreover, (<ref>) becomes a linear system with constant coefficients given by

    dda	X_B(a) = M X_B(a)  + f_B(a),

where 

    M=
      [ -(ϕ+γ)     Bρ;    ϕ+γ    -Bρ;        ],  f_B(a)=
      [ Bβ s_B^*(a);           0;             ].

Observe that the eigenvalues of M are given by λ_1 = 0 and λ
_2 = - (ϕ+γ) - B ρ, with corresponding eigenvectors v_1 = [Bρ, ϕ+γ]^T, v_2 = [-1, 1]^T. Then, the solution of (<ref>) is given by

    i_B^*(a) =     Bρϕ+γ+Bρ-B(ρ-β)(ϕ+γ+B(ρ-β))· e^-Bβ a
       -B β (ϕ+γ) (ϕ+γ+Bρ)(ϕ+γ+B(ρ-β))· e^-(ϕ+γ+Bρ)a,

and 

    r_B^*(a) =    ϕ+γϕ+γ+Bρ-ϕ+γϕ+γ+B(ρ-β)· e^-Bβ a
       +B β (ϕ+γ) (ϕ+γ+Bρ)(ϕ+γ+B(ρ-β))· e^-(ϕ+γ+Bρ)a.

Note that 
    p_∞(a) = lim_t→∞ p(t,a)=ℱ(a)/∫_0^∞ℱ(a) da = μ e^-μ a.
 

After some simplifications, we obtain that 

    G(B) = B+μρ(B+μβ)( B+μ+ϕ+γρ).

In particular, the basic reproductive number is given by

    ℛ_0 = G(0) = βμ +ϕ+γ   and    ℛ_C = βμ.


It can be proven that the equation G(B)=1 has two real roots B∈ (0,1) if and only if the following conditions are satisfied:

    ℛ_0 <1 < ℛ_C,
    
       β > ρμ/μ+(√(ρ)-√(ϕ+γ))^2.





 

 

These conditions can be represented graphically as the interior of the volume depicted in <Ref> in the space (ϕ+γ/μ,β/μ,ρ/μ). In <Ref> we present a typical graph for G(B) with constant parameters (ϕ,γ,β,ρ,μ) satisfying (<ref>). It is clear that G(B) = B has two different solutions, and for each value there is a non-trivial solution (s^*(a),i^*(a),r^*(a)) as shown in <Ref>. We remark that conditions (<ref>) are necessary and sufficient in order to guarantee existence of two non-zero steady states.

 


In <Ref> we show a backward bifurcation obtained by fixing μ, ρ, ϕ, γ as given in (<ref>) and varying β. In this case, we have infected (infectious) individuals i^* as a function of age a and ℛ_0 = β/μ+ϕ+γ. This is typical behavior when we have a backward bifurcation and ℛ_0<1, which implies ℛ_0<1 is not a sufficient condition to control the disease. Hence, the condition of ℛ_0 ≪ 1 is necessary to reach a infection-free state (see Region 1 in <Ref>). When conditions (<ref>) are met, we observe two non-trivial solutions for i^*(a) (see region 2). This implies that under the scenario where there is a critical mass of infectious individuals in a population, it would make the control of the disease more difficult. This is relevant for entities responsible to develop prevention/control strategies. For ℛ_0>1, <Ref> implies the existence of a non-uniform steady state distribution (see Region 3). 



§ NUMERICAL EXPERIMENTS
 
In this section, we first describe the numerical implementation for solving system (<ref>). We then present different examples confirming our theoretical results obtained in <Ref> and <Ref>. All parameters are in units of year^-1. After analyzing the model in its general form, we implement a numerical example based on the analogous nonlinear ordinary differential equation model in <cit.>. That is, we look at the nonlinear dynamics of drinkers as an epidemicprocess where individuals who are considered problem drinkers can cause other individuals to start drinking. Eventually, these individuals can temporarily recover but as state in <cit.>, relapse rates are high and the probability of never drinking again is small.



 §.§ Numerical implementation
 
We will discretize (<ref>) with a first-order upwind finite difference scheme and will approximate the solution on the physical domain of interest given by the rectangle { (t,a)∈ [0,T]× [0,A]}. We first construct a uniform grid with equidistant points. Divide [0,T] and [0,A] into N_T and N_A subintervals, respectively. Thus, the nodes (t_j,a_k) on the rectangular mesh are given by 

    (t_j,a_k) = ( jΔ t, kΔ a),

for j∈{ 0,1,…,N_T}, k∈{ 0,1,…,N_A}, where 

    Δ t := T/N_T,  Δ a := A/N_A
 are the corresponding step sizes. We require that Δ t <Δ a in order to satisfy the CFL and stability conditions of the scheme.

For any function x and a grid point (t_j,a_k), we denote the approximation of x(t_j,a_k) by x_k^j. 
Since

    ( ∂/∂ t + ∂/∂ a) x(t_j,a_k) =    x(t_j+Δ t, a_k)-x(t_j,a_k)/Δ t + 𝒪(Δ t) 
          +  x(t_j,a_k)- x(t_j,a_k-Δ a)/Δ a + 𝒪(Δ a),

we approximate the derivatives by

    ( ∂/∂ t + ∂/∂ a) x(t_j,a_k) ≈x^j+1_k - x^j_k/Δ t + x^j_k - x^j_k-1/Δ a.

Similarly, define β_k:=β(a_k), ϕ_k:=ϕ(a_k),  γ_k:=γ(a_k), ρ_k:=ρ(a_k) and B^j:=B(t_j). Then, by evaluating at all the grid points, the discretization for system (<ref>) is given by the explicit system

    s^j+1_k - s^j_k/Δ t + s^j_k - s^j_k-1/Δ a   = - β_k s^j_k B^j,
    i^j+1_k - i^j_k/Δ t + i^j_k - i^j_k-1/Δ a   = β_k s^j_k B^j - (ϕ_k+γ_k) i^j_k+ρ_k r_k^j B^j , 
    r^j+1_k - r^j_k/Δ t + r^j_k - r^j_k-1/Δ a   = (ϕ_k+γ_k) i^j_k-ρ_k r_k^j B^j,

for 1≤ k≤ N_A and 0≤ j≤ N_T-1. We recall that the initial conditions (<ref>) provide the values s_0^j, s_k^0, i_0^j, i_k^0, r_0^j, and r_k^0. The integral 
    B^j = B(t_j) = ∫_0^∞ i(t_j,a) p(t_j,a) da
 is approximated via MATLAB's command . We note that for a fixed time t_j, B^j depends only on values of i at the same time t_j.

Solving for s^j+1_k, i^j+1_k and r^j+1_k from (<ref>), we obtain that


    s^j+1_k     = s^j_k + Δ t (- β_k s^j_k B^j  - s^j_k - s^j_k-1/Δ a),
    
    
    i^j+1_k    = i^j_k +  Δ t( β_k s^j_k B^j - (ϕ_k+γ_k) i^j_k+ρ_k r_k^j B^j -  i^j_k - i^j_k-1/Δ a),
    
     
    r^j+1_k    = r^j_k+ Δ t( (ϕ_k+γ_k) i^j_k-ρ_k r_k^j B^j  -  r^j_k - r^j_k-1/Δ a),


for 1≤ k≤ N_A and 0≤ j≤ N_T-1. Therefore, starting with the initial conditions s_0^j, s_k^0, i_0^j, i_k^0, r_0^j, r_k^0, we can compute, in successive time steps, the values of the unknowns on the grid points. The convergence of the scheme has been established in <cit.>.



 §.§ Global stability
 
In this example we confirm the result shown in <Ref>, where global stability 
to the infection-free state
is guaranteed when ℛ_C<1. We consider constant coefficients given by

    [ μ = 0.0125,   β= 0.011,      ϕ= 60,       γ=13,    ρ=76.65. ]

Thus, ℛ_0 ≈  1.5· 10^-4 and ℛ_C = 0.88. We then obtain i(t,a) as shown in <Ref>. The parameter values used in our model are based on  <cit.>.

 

For the given set of parameters all solutions for i(t,a) go to zero as t→∞. In this case, we are in the first region of the backward bifurcation (see <Ref>) where ℛ_0 ≪ 1. Here β<ρ and β is not sufficiently large to sustain enough drinkers in the system.



 §.§ Backward bifurcation
 
We now consider a particular choice of parameters that satisfy conditions (<ref>). In particular, we fix 

    [ μ = 0.0125,      β= 60,      ϕ= 60,       γ=13,    ρ=76.65. ]

For this case, 

    ℛ_0 ≈ 0.8218   and   ℛ_C = 4800.



We are now in region 2 of <Ref>, where we need ℛ_0<1<ℛ_C. We then have two non-trivial steady-state solutions for system (<ref>) depicted in <Ref>. We then obtain i(t,a) as described in <Ref> for two different initial conditions i(0,a) ; see <Ref>. 







 
 

In <Ref> we obtain the stable solution that corresponds to the steady-state distribution shown in <Ref>. For a smaller initial condition i(0,a), we confirm that the steady-state distribution shown in <Ref> is unstable, since i(t,a) tends to the (stable) zero solution; see <Ref>. In this case, albeit having β<ρ, β is sufficiently large to sustain multiple drinking steady states dependent on the initial state of the drinking population. That is, ℛ_0<1 is not a sufficient condition for the stability of the drinking-free state. 















 §.§ Endemic non-uniform steady state distribution
  
According to <Ref>, there exists a non-uniform steady state distribution if ℛ_0>1. In this example we take

    [ μ = 0.0125,     β= 120,      ϕ= 60,       γ=13,    ρ=76.65, ]

for which ℛ_0 ≈  1.6436.

We obtain i(t,a) as shown in <Ref>. It is clear that i(t,a) converges to a non-zero solution.
 

Here β>ρ and we have ℛ_0>1, which is sufficient for the drinking steady state to prevail and the initial conditions do not play a role (see <Ref>).



 §.§ Numerical experiments with age-dependent parameters
 




Lastly, we present an example where the parameters ϕ(a), γ(a) β(a), ρ(a), μ(a) and c(a) are age-dependent with particular distributions; see <Ref>. These were chosen arbitrarily, however, we pursued biological significance when choosing each parameter distribution. 




 

 
 
In <Ref>, we see that the solution i(t,a) goes to the infection-free state when ℛ_C<1, as shown in <Ref>, where global stability was attained. <Ref> illustrates an endemic state when ℛ_0>1; see <Ref>. The case ℛ_0<1<ℛ_C is shown in <Ref>; by using different initial conditions we have the possibility of multiple steady states as shown in the bifurcation diagram <Ref>.


Overall, when parameters are functions of a we can obtain the same behavior as when parameters are constant. This is more realistic but has proven to be challenging to obtain explicit conditions for each case. Numerically, albeit challenging, it is possible to explore the parameter space of the bifurcation <Ref> and produce similar results. Numerical results are highly dependent on the approximation of B(t), which depends on the unknown i(t,a). In order to compute efficiently and accurately B^j, we have used Chebyshev expansions based on 's library  <cit.>. <Ref> was computed by solving G(B) = B for several values of β, where (s^*(a),i^*(a),r^*(a)) are computed by using the command , which allows to improve running times. Further exploration of the model with variable coefficients is required. 
 


§ CONCLUSIONS
 
We have studied an age-structured epidemic model with nonlinear recidivism. The infection-free steady state distribution was computed, as well as its local and global stability. The existence of the endemic non-uniform steady state distribution is guaranteed when ℛ_0>1. Moreover, our analysis shows the existence of multiple endemic equilibria when ℛ_0<1, and necessary conditions that lead to a backward bifurcation were computed. As a numerical example we used alcohol dynamics to showcase the results of the model.

Numerical experiments were conducted to illustrate the different scenarios where there exist multiple drinking states. Typically, ℛ_0<1 is a sufficient condition for a disease" to die out. However, when incorporating social factors, in this case there exists the possibility of an endemic state even when ℛ_0<1. In our model, the basic reproductive number is a function of the treatment and recovery rate of individuals. 

We have necessary conditions for two positive endemic steady states, which highlights the importance and significance of the initial density of the drinking population. In this case, it is possible to have an endemic state when ℛ_0<1 if the initial number of drinkers is large enough, as seen in <Ref>. Short treatment periods, in turn mostly ineffective, induce more individuals into the temporarily recovered class where the likelihood of relapse is high. This creates a population of new susceptible individuals whose susceptibility and risk can be measured by the strength of environmentally induced recidivism rates. Therefore making prevention strategies far more challenging.

The implications of an age-structure model with nonlinear recidivism could lead to a better understanding of applications where age is an important factor when implementing prevention/control strategies. 





§ ACKNOWLEDGEMENTS

The authors would like to thank the Research Center in Pure and Applied Mathematics and the Mathematics Department at Universidad de Costa Rica for their support during the preparation of this manuscript. 



99
allen2008mathematical
Allen, L. J., Brauer, F., Van den Driessche, P., and Wu, J. (2008).
Mathematical Epidemiology, Vol. 1945.,
Springer.


arino2003global
Arino, J., McCluskey, C. C., and van den Driessche, P. (2003).
Global results for an epidemic model with vaccination that exhibits backward bifurcation.
SIAM J. Appl. Math. 64(1):260–276.

brauer2004
Brauer, F. (2004).
Backward bifurcations in simple vaccination models.
J Math Anal Appl. 298(2):418 – 431.

brauer2011
Brauer, F. (2011).
Backward bifurcations in simple vaccination/treatment models.
J Biol Dyn. 5(5):410–418.

ccc2002
Castillo-Chavez, C. and Huang, W. (2002).
Age-structured core group model and its impact on std dynamics.
Mathematical Approaches for Emerging and Reemerging Infectious Diseases: Models, Methods, and Theory, Editors: Castillo-Chavez, C., Blower, S., Driessche, P. van d., Kirschner, D. and Yakubu, A.-A.

ccc2004
Castillo-Chavez, C., and Song, B. (2004).
Dynamical models of tuberculosis and their applications.
Math Biosci Eng. 1(9):361–404.

Driscoll2014
Driscoll, T. A., Hale, N., and Trefethen, L. N. (2014).
Chebfun Guide.
Pafnuty Publications.

feng2000model
Feng, Z., Castillo-Chavez, C., and Capurro, A. F. (2000).
A model for tuberculosis with exogenous reinfection.
Theor Popul Biol. 57(3):235–247.

finney1999relapse
Finney, J, M. R. T. C. (1999).
The course of treated and untreated substance use disorders: remission and resolution, relapse and mortality.
Oxford University Press.

holmes1994partial
Holmes, E. E., Lewis, M. A., Banks, J., and Veit, R. (1994).
Partial differential equations in ecology: spatial interactions and population dynamics.
Ecology. 75(1):17–29.

jin1998relapse
Jin, H., Rourke, S. B., Patterson, T. L., Taylor, M. J., and Grant, I. (1998).
Predictors of relapse in long-term abstinent alcoholics.
J Stud Alcohol. 59:640–646.

kot2001elements
Kot, M. (2001).
Elements of Mathematical Ecology.
Cambridge University Press.

kribs2002
Kribs-Zaleta, C. M., and Martcheva, M. (2002).
Vaccination strategies and backward bifurcation in an age-since-infection structured model.
Math Biosci. 177-178(Supplement C):317–332.
  
martcheva2003progression
Martcheva, M., and Thieme, H. R. (2003).
Progression age enhanced backward bifurcation in an epidemic model with super-infection.
J Math Biol. 46(5):385–424.  

moss2007relapse
Moss, R. H., and Moss, B. S. (2007).
Rates and predictors of relapse after natural and treated remission from alcohol use disorders.
Addiction. 101(2):212–222.

sanchez2018x
Sanchez, F. and Calvo, J.G. (2018).
The role of short-term immigration on disease dynamics: An SIR model with age-structure.
arXiv:1811.05923
  
sanchez2007
Sanchez, F., Wang, X., Castillo-Chavez, C., Gorman, D., and Gruenewald, P. (2007).
Drinking as an Epidemic-A Simple Mathematical Model with Recovery and Relapse.
In: Witkiewitz, K.A., Marlatt, G.A., eds. Therapistś Guide to Evidence-Based Relapse Prevention
Elsevier Inc., pp. 353–368.

shim2006
Shim, E., Feng, Z., Martcheva, M., and Castillo-Chavez, C. (2006).
An age-structured epidemic model of rotavirus with vaccination.
J. Math. Biol. 53(4):719–746.

song2006
Song, B., Castillo-Garsow, M., Rios-Soto, K., Mejran, M., Henso, L., and Castillo-Chavez, C. (2006).
Raves, clubs and ecstasy: the impact of peer pressure.
Math Biosci Eng. 3(1):249–266.
  
song2013
Song, B., Du, W., and Lou, J. (2013).
Different types of backward bifurcations due to density-dependent treatments.
Math Biosci Eng. 10(5-6):1651—1668.



