A mathematical model with nonlinear relapse: conditions for a forward-backward bifurcation

ABSTRACT

We constructed a Susceptible-Addicted-Reformed model and explored the dynamics of nonlinear relapse in the Reformed population. The transition from susceptible considered at-risk is modeled using a strictly decreasing general function, mimicking an influential factor that reduces the flow into the addicted class. The basic reproductive number is computed, which determines the local asymptotically stability of the addicted-free equilibrium. Conditions for a forward-backward bifurcation were established using the basic reproductive number and other threshold quantities. A stochastic version of the model is presented, and some numerical examples are shown. Results showed that the influence of the temporarily reformed individuals is highly sensitive to the initial addicted population.

KEYWORDS:

• Nonlinear relapse
• backward bifurcation
• epidemic models
• social determinants
• addiction

1. Introduction

Mathematical models applied to infectious diseases have become common; more recently, an insurmountable number of models arose during the Covid-19 pandemic [4,5,10,15–17,27] (and references therein). In general, when studying infectious diseases, mathematical models help understand disease transmission dynamics. Furthermore, in some cases, mathematical models can provide insight to health authorities to construct and develop efficient public health policies [10].

Infectious diseases have been a burden to public health for some time. The transmission mechanisms of pathogens are mainly close contact with an infectious host, airborne, via a vector, and in some cases via contact with an infected area [2]. However, more recently, health authorities worldwide have been vigilant with the high number of mental health incidents, highlighted primarily by the Covid-19 pandemic [33]. Social factors provide a unique challenge to construct mathematical models that include social aspects not typically included in epidemic models. However, the incorporation of social determinants in these models is inherently difficult. In previous work, social determinants were introduced as ‘epidemics’ where transmission happened through social interactions, similar to infectious diseases. For example, a drinking dynamics model using a nonlinear system of differential equations [26], where the ‘infectious’ class was the drinking population, and the interaction between nondrinkers and drinkers simulated an epidemic process. We based our theoretical framework on the latter. Other models simulating social dynamics include: a bulimia model [11], drug models [1,31], a sex worker industry model [7], and smoking models [6,23,28–30] among others.

Modeling social interactions as epidemic processes can provide a helpful understanding of the phenomenon of interest. Here, we model addiction as an infectious disease where the interactions between the non-addicted and addicted individuals can cause an ‘epidemic’ process and confer an ‘infection’. Drug addiction has been a problem worldwide for many decades [3,24,32]. In particular, when the crack ‘epidemic’ of the 1980s was in full force, the derivative of cocaine, a pure and more expensive narcotic, led to a faster addiction and deterioration of individuals that consumed the drug [9,12].

Furthermore, relapse rates of addicted individuals, especially those that used potent narcotics such as crack cocaine, methamphetamine, fentanyl, and heroin, among others, are very high [18,20]. In the model constructed here, we looked at nonlinear relapse rates and the influence of those who recovered and want to provide support for non-addicted, presumed susceptible individuals. This is done via a general function that depends on the temporarily recovered population and other parameters.

Epidemic models have helped describe the transmission dynamics of infectious pathogens and derive strategies for their control, prevention, and reduction of incidence, among others. Here, we provide a theoretical framework to study social phenomena via an epidemic model and highlight the sensitivity of initial conditions.

The article is organized as follows: in Section 2, we give details of the mathematical model. In Section 3, we present the mathematical analysis. In Section 4, we provide a stochastic version of the model and provide some numerical examples. Finally, in Section 5, we provide a discussion based on our results.

2. Mathematical model

An addiction model is a mathematical framework used to describe the dynamics of substance abuse and addiction in a population. The model assumes that the total population is divided into three compartments: susceptible individuals (S), addicted individuals (A), and temporarily reformed individuals ($\tilde{S}$). The model we consider is based on [26], where authors explored the impact of nonlinear influence on drinking behaviour dynamics.

The model is similar to the SIR model used to describe the spread of infectious diseases [2,14], but with some key differences. In an addiction model, the transition from susceptible to addicted is not solely determined by contact with an addicted individual but rather by a complex interplay of individual and environmental factors, such as genetics, social norms, and substance access.

The recruitment rate, β, represents the strength of social influence on susceptible (at-risk) individuals. In this context, transmission is a collective behaviour rather than an individual consequence; i.e. recruitment is not typically the work of a single individual, but instead is a result of the collective influence of a group of individuals as a whole [8]. Moreover, $\kappa \in [0,1]$ denotes the cost of addiction, and $\nu \in [0,1]$ is the willingness of reformed individuals to deter at-risk individuals from addiction. We then consider a positive, strictly decreasing smooth function g defined by: (1) $g_{\kappa ,\nu }\left(\tilde{S}\right)=\frac{\kappa }{1+\nu \frac{\tilde{S}}{N}},$(1) which acts as a reducing factor that impacts the at-risk population by reducing the probability that an individual transitions to the addicted class.

To further explain the function g, it is important to note that the function is a decreasing function of $\tilde{S}$. As the proportion of reformed individuals in the population increases (i.e. as $\tilde{S}$ increases), the reducing factor provided by g decreases, reducing the probability that susceptible individuals transition to the addicted class.

Moreover, the reducing effect of g on the addicted population is greater for larger values of ν, thus a greater willingness of reformed individuals to deter susceptible (at-risk) individuals from addiction lead to a stronger reduction in the addicted population. Similarly, larger values of κ correspond to a greater cost associated with addiction.

Overall, the function g captures the impact of reformed individuals on addiction dynamics by providing a mechanism through which they can influence the at-risk population and reduce the probability of transitioning to the addicted class.

The model is governed by a set of ordinary differential equations, which describe the rates of change of the number of individuals in each compartment over time. These equations are: (2) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mu N-\text{β}g\left(\tilde{S}\right)S\frac{A}{N}-\mu S,\\ \frac{\mathrm{d}A}{\mathrm{d}t}& =\text{β}g\left(\tilde{S}\right)S\frac{A}{N}+\varphi \tilde{S}\frac{A}{N}-(\mu +\gamma )A,\\ \frac{\mathrm{d}\tilde{S}}{\mathrm{d}t}& =\gamma A-\varphi \tilde{S}\frac{A}{N}-\mu \tilde{S},\end{array}$(2) where $N=S+A+\tilde{S}$ is the total (assumed constant) population. For simplicity, we re-scale system (2(2) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mu N-\text{β}g\left(\tilde{S}\right)S\frac{A}{N}-\mu S,\\ \frac{\mathrm{d}A}{\mathrm{d}t}& =\text{β}g\left(\tilde{S}\right)S\frac{A}{N}+\varphi \tilde{S}\frac{A}{N}-(\mu +\gamma )A,\\ \frac{\mathrm{d}\tilde{S}}{\mathrm{d}t}& =\gamma A-\varphi \tilde{S}\frac{A}{N}-\mu \tilde{S},\end{array}$(2) ) by substituting $s=\frac{S}{N}$, $a=\frac{A}{N}$, $\tilde{s}=\frac{\tilde{S}}{N}$, in order to use percentages instead of absolute numbers. We then obtain the equivalent system of equations: (3a) $\begin{array}{rl}\frac{\mathrm{d}s}{\mathrm{d}t}& =\mu -\text{β}g\left(\tilde{s}\right)sa-\mu s,\end{array}$(3a) (3b) $\begin{array}{rl}\frac{\mathrm{d}a}{\mathrm{d}t}& =\text{β}g\left(\tilde{s}\right)sa+\varphi \tilde{s}a-(\mu +\gamma )a,\end{array}$(3b) (3c) $\begin{array}{rl}\frac{\mathrm{d}\tilde{s}}{\mathrm{d}t}& =\gamma a-\varphi \tilde{s}a-\mu \tilde{s}.\end{array}$(3c) It is clear that $s+a+\tilde{s}=1$ and the reducing factor (1(1) $g_{\kappa ,\nu }\left(\tilde{S}\right)=\frac{\kappa }{1+\nu \frac{\tilde{S}}{N}},$(1) ) is therefore given by (4) $g\left(\tilde{s}\right)=\frac{\kappa }{1+\nu \tilde{s}}\in [0,1].$(4) The first equation represents the rate of change of the susceptible population, which is proportional to the number of susceptible individuals who come into contact with addicted individuals. The second equation represents the rate of change of the addicted population, which is proportional to the number of susceptible individuals who become addicted minus the number of addicted individuals who temporarily recover. The third equation represents the reformed population's change rate, which is proportional to the number of addicted individuals who recover.

Examples of possible applications to the model include drugs like heroin, fentanyl, methamphetamine, and others where relapse is common. The model can provide insights into the dynamics of addiction and the impact of interventions such as harm reduction strategies and substance abuse treatment.

3. Mathematical analysis

We first analyze the addiction-free equilibrium, $(s_0^{\ast },a_0^{\ast },{\tilde{s}}_0^{\ast })=(1,0,0)$, that can be used to determine the basic reproductive number, $R_0$. In epidemiology, the basic reproductive number represents the number of secondary infections produced by an average infected individual; the disease typically dies out when this number is less than one. When it is greater than one, there will be an epidemic [2]. In this context, we consider $R_0$ to measure the strength of the social influence of addicted individuals to recruit individuals into a vice. As we will demonstrate, $R_0>1$ implies establishing an infectious agent, and $R_0<1$ typically implies that the number of addicted individuals decreases and goes to zero. Our model can sustain an addiction when $R_0<1$ under particular initial conditions.

3.1. Basic reproductive number and addiction-free equilibrium

We use the next generation operator method [14] to compute $R_0$. Let $\mathcal{F}=\text{β}g\left(\tilde{s}\right)sa+\varphi \tilde{s}a\mathrm{a}\mathrm{n}\mathrm{d}\mathcal{V}=(\mu +\gamma )a,$where $\mathcal{F}$ and $\mathcal{V}$ contains all terms flowing into a and flowing out of a, respectively. It holds that $F=\frac{\mathrm{\partial }\mathcal{F}}{\mathrm{\partial }a}{|}_{(s_0^{\ast },a_0^{\ast },{\tilde{s}}_0^{\ast })}=\text{β}g\left(0\right)\mathrm{a}\mathrm{n}\mathrm{d}V=\frac{\mathrm{\partial }\mathcal{V}}{\mathrm{\partial }a}{|}_{(s_0^{\ast },a_0^{\ast },{\tilde{s}}_0^{\ast })}=\mu +\gamma .$The basic reproductive number is then: $R_0=F{V}^{-1}=\frac{\text{β}g\left(0\right)}{\mu +\gamma }=\frac{\text{β}\kappa }{\mu +\gamma },$where $\frac{1}{\mu +\gamma }$ represents the average amount of time spent in the addicted class. We then have the following result:

Theorem 3.1

The addiction-free equilibrium is stable if and only if $R_0<1$.

Proof.

The Jacobian of system (3a(3a) $\begin{array}{rl}\frac{\mathrm{d}s}{\mathrm{d}t}& =\mu -\text{β}g\left(\tilde{s}\right)sa-\mu s,\end{array}$(3a) ) is given by: $J(s,a,\tilde{s})=\left[\begin{array}{ccc}{\displaystyle \frac{-\text{β}\kappa }{1+\nu \tilde{s}}}a-\mu & {\displaystyle \frac{-\text{β}\kappa }{1+\nu \tilde{s}}}s& {\displaystyle \frac{\text{β}\kappa }{(1+\nu \tilde{s}{)}^2}}sa\nu \\ {\displaystyle \frac{\text{β}\kappa }{1+\nu \tilde{s}}}a& {\displaystyle \frac{\text{β}\kappa }{1+\nu \tilde{s}}}s+\varphi \tilde{s}-(\mu +\gamma )& {\displaystyle \frac{-\text{β}\kappa }{(1+\nu \tilde{s}{)}^2}}sa\nu +\varphi a\\ 0& \gamma -\varphi \tilde{s}& -\varphi a-\mu \end{array}\right],$and evaluating the addicted-free equilibrium yields: $J(1,0,0)=\left[\begin{array}{ccc}-\mu & -\text{β}\kappa & 0\\ 0& \text{β}\kappa -(\mu +\gamma )& 0\\ 0& \gamma & -\mu \end{array}\right].$The eigenvalues of this matrix are: $\begin{array}{rl}{\lambda }_1,{\lambda }_2& =-\mu <0,\\ {\lambda }_3& =\text{β}\kappa -(\mu +\gamma ).\end{array}$Note that ${\lambda }_3<0⟺R_0<1$, and the result holds since the equilibrium is stable if all the eigenvalues are negative.

3.2. Endemic equilibria

To study the prevalence of addiction and endemic equilibria in our model, we need to define an analogous basic reproductive number for the reformed class, denoted by $R_{\varphi }$, that represents the strength of social influence of addicted individuals to recruit reformed individuals back into addiction. It is defined by $R_{\varphi }=\frac{\varphi }{\mu +\gamma }.$Since reformed individuals have already been exposed to a previous addiction, we assume that the rate at which a reformed individual relapse is higher than the rate at which an at-risk individual becomes addicted, i.e. $R_{\varphi }>R_0$. As we will show, endemic equilibria exist whenever $R_0>1$, and under certain initial conditions when $R_0<1$ and $R_{\varphi }>1$.

We study two cases: $\nu =0$ (absence of willingness factor) and $0<\nu \le 1$. In the first case, the social influence of reformed individuals on the at-risk population is absent; this is, $g\left(\tilde{s}\right)=\kappa$. These two cases allow us to explore reformed individuals' impact on population dynamics.

3.2.1. Case $\nu =0$: absence of willingness factor

Solving for the endemic equilibria of system (3) when $\nu =0$ yields to: $s_1^{\ast }=\frac{\mu }{\text{β}\kappa a+\mu },{\tilde{s}}_1^{\ast }=\frac{\gamma a}{\varphi a+\mu }.$Substituting these values into (3b(3b) $\begin{array}{rl}\frac{\mathrm{d}a}{\mathrm{d}t}& =\text{β}g\left(\tilde{s}\right)sa+\varphi \tilde{s}a-(\mu +\gamma )a,\end{array}$(3b) ) yields the quadratic equation $x_2{a}^2+x_{1}a+x_0=0$, where the coefficients are given by: $\begin{array}{rl}{x}_2& =R_{\varphi }{R}_0,\\ x_1& =R_0(1-R_{\varphi })+R_{\mu }{R}_{\varphi },\\ x_0& =R_{\mu }(1-R_0),\end{array}$where $R_{\mu }=\mu /(\mu +\gamma )\in (0,1)$. In this case, our system is very similar to the drinking model studied in [26]. We can construct a bifurcation diagram to analyze the stability of the endemic equilibria as a function of $R_0$ (as κ varies). Stability depends on the value of $R_0$ and the initial addicted population size. In Figure 1, we show a typical backward bifurcation curve, typically occurring in systems with nonlinear relapse rates [13,26,34]. Furthermore, the system exhibits hysteresis; i.e. it is highly sensitive to initial conditions [14]. After straightforward computations, the quadratic equation has a double root when $R_0$ is equal to (5) $R_c=R_{\mu }{R}_{\varphi }\frac{1+R_{\varphi }+2\sqrt{R_{\varphi }(1-R_{\mu })}}{(R_{\varphi }-1{)}^2+4R_{\mu }{R}_{\varphi }}.$(5) There are no positive endemic equilibria if $R_0<R_c$, and it can be shown that there are two positive endemic equilibria when $R_c<R_0<1$; see [26].

Figure 1. Backward bifurcation with parameters $\mu =0.00015$, $\text{β}=0.009$, $\gamma =0.0027$, $\nu =0$, $\varphi =0.005$, and κ varies. The dotted vertical line represents the critical value $R_c\approx 0.52$ for which there are no positive endemic equilibria if $R_0<R_c$. There are two positive endemic equilibria when $R_c<R_0<1$. Parameter values were taken from [26].

3.2.2. Case $0<\nu \le 1$: presence of willingness factor

Analyzing the system when $0<\nu \le 1$ allows us to explore the impact of the social influence of reformed individuals on the at-risk population. Solving for the endemic equilibria of system (3a(3a) $\begin{array}{rl}\frac{\mathrm{d}s}{\mathrm{d}t}& =\mu -\text{β}g\left(\tilde{s}\right)sa-\mu s,\end{array}$(3a) ) when $0<\nu \le 1$ leads to the following: $s_2^{\ast }=\frac{\mu D}{\text{β}\kappa a+\mu D},{\tilde{s}}_2^{\ast }=\frac{\gamma a}{\varphi a+\mu },$where $D=1+\frac{\nu \gamma a}{\varphi a+\mu }$. Substituting these values into Equation (3b(3b) $\begin{array}{rl}\frac{\mathrm{d}a}{\mathrm{d}t}& =\text{β}g\left(\tilde{s}\right)sa+\varphi \tilde{s}a-(\mu +\gamma )a,\end{array}$(3b) ) yields the cubic equation $f\left(x\right)=0$, where (6) $f\left(x\right)=x_3{a}^3+x_2{a}^2+x_{1}a+x_0,$(6) and the coefficients are given by: $\begin{array}{rl}{x}_3& =R_{\varphi }^2{R}_0,\\ x_2& =R_{\varphi }[R_0(1-R_{\varphi })+R_{\mu }(R_0+R_{\varphi })+\nu R_{\mu }(1-R_{\mu }\left)\right],\\ x_1& =R_{\mu }\left[\nu \right(1-R_{\mu })+R_0(1-R_{\varphi })+R_{\varphi }(1-R_0)+R_{\mu }{R}_{\varphi }],\\ x_0& =R_{\mu }^2(1-R_0).\end{array}$Figure 2 shows a typical bifurcation exhibiting both a forward and a backward behaviour for model (3) as a function of $R_0=R_0\left(\kappa \right)$. Let $R_c$ and $R_0^{\ast }$ (with $R_c<R_0^{\ast }$) be the values of $R_0$ for which (6(6) $f\left(x\right)=x_3{a}^3+x_2{a}^2+x_{1}a+x_0,$(6) ) has two double roots, similarly as (5(5) $R_c=R_{\mu }{R}_{\varphi }\frac{1+R_{\varphi }+2\sqrt{R_{\varphi }(1-R_{\mu })}}{(R_{\varphi }-1{)}^2+4R_{\mu }{R}_{\varphi }}.$(5) ). These two constants are thresholds determining the endemic equilibria for a given value of $R_0$. We remark that:

1. If $R_0>1$, there exists at least one positive equilibrium state, since $f\left(0\right)=x_0<0$ and $f\left(1\right)=R_0(1-R_{\mu })(R_{\mu }+R_{\varphi })+R_{\mu }(1+R_{\varphi })(R_{\mu }+R_{\varphi }+\nu (1-R_{\mu }\left)\right)>0$.

2. If $R_0<1$ and $R_{\varphi }<1$, the coefficients of f are all positive, and therefore there is no positive equilibrium state.

Figure 2. Forward-backward bifurcation with parameters $\mu =0.00015$, $\gamma =0.0027$, $\text{β}=0.009$, $\nu =0.8$, $\varphi =0.0044$ and κ varies. The dotted lines ($R_0=R_c$, $R_0=1$ and $R_0=R_0^{\ast }$ from left to right) separate the domain into four regions: Region 1 ($R_0<R_c$) with no positive equilibria, Region 2 ($R_c<R_0<1$) with two positive equilibria, Region 3 ($1<R_0<R_0^{\ast }$) with three positive equilibria, and Region 4 ($R_0^{\ast }<R_0$) with one positive equilibrium. Parameter values were taken from [26].

We now demonstrate conditions for the number of endemic equilibria as a function of $R_0$ and $R_{\varphi }$, according to each region of interest as depicted in Figure 2. In Region 2, two endemic equilibria exist, and the system exhibits backward behaviour. In this case, the initial addicted population size determines if the addicted population can establish itself or decrease to zero:

Theorem 3.2

If $\nu >0$, a necessary condition for two positive equilibria is $0<R_c<R_0<1$ and $R_{\varphi }>1$.

Proof.

If the system has two positive equilibria, the polynomial f in Equation (6(6) $f\left(x\right)=x_3{a}^3+x_2{a}^2+x_{1}a+x_0,$(6) ) has three real roots, which implies $R_0>R_c$. In addition, since $x_3>0$, only one of these equilibria has to be negative, which implies that $f\left(0\right)={\mu }^2(1-R_0)>0$, which implies $R_0<1$.

In Region 3, three endemic equilibria exist (two stable and one unstable). This implies that two end-states can occur and that the long-term population can establish itself at either a large or small size, depending on the initial addicted population size:

Theorem 3.3

If $\nu >0$, a necessary condition for three positive equilibria is $1<R_0<R_0^{\ast }$ and $R_{\varphi }>1$.

Proof.

If the system has three positive equilibria, then the polynomial f in Equation (6(6) $f\left(x\right)=x_3{a}^3+x_2{a}^2+x_{1}a+x_0,$(6) ) has three real equilibria, which implies $R_c<R_0<R_0^{\ast }$. If these equilibria are all positive, in particular, this implies $f\left(0\right)={\mu }^2(1-R_0)<0$, which implies $R_0>1$.

In Region 4, a unique endemic equilibrium exists, and the addicted population will establish itself, regardless of the initial addicted population size:

Theorem 3.4

If $\nu >0$, a sufficient condition for a unique positive equilibrium is $1<{\mathcal{R}}_0^{\ast }<{\mathcal{R}}_0$ and ${\mathcal{R}}_{\varphi }>1$.

Proof.

If ${\mathcal{R}}_0>{\mathcal{R}}_0^{\ast }$, then the polynomial f in Equation (6(6) $f\left(x\right)=x_3{a}^3+x_2{a}^2+x_{1}a+x_0,$(6) ) has a single real root $a_{+}$. Since ${\mathcal{R}}_0>1$, this implies that $f\left(0\right)={\mu }^2(1-{\mathcal{R}}_0)<1$, which implies $a_{+}>0$.

3.3. Effect of reducing the relapse rate and the willingness factor

We first explore the impact of relapse rate, ϕ, on the model. Reducing the relapse rate while still maintaining ${\mathcal{R}}_{\varphi }>1$ results in a change in the behaviour of the system as shown in Figure 3. A lower relapse rate causes the bifurcation to shift to the right, where ${\mathcal{R}}_0<1$ guarantees an infectious-free equilibrium. This highlights the crucial role that relapse plays in addicted population dynamics. If relapse can be lowered below a critical threshold, the infectious population may be managed by just controlling $R_0$.

Figure 3. Forward-backward bifurcation with parameters $\mu =0.00015$, $\gamma =0.0027$, $\text{β}=0.009$, $\nu =0.8$, κ varies, and (a) $\varphi =0.0044$, ${\mathcal{R}}_{\varphi }=1.54$ ; (b) $\varphi =0.004$, ${\mathcal{R}}_{\varphi }=1.4035$. Parameter values were taken from [26].

Varying the willingness factor, ν, yields significant changes in the behaviour of the model; see Figure 4 where we show the bifurcation diagrams for four different values of ν. Low values of ν (indicative of low interaction between the reformed and at-risk classes) yield a backward bifurcation similar to Figure 1. As ν increases, the system moves through the state shown in Figure 3(b) and continues to shift to the right. When $\nu =1$, the state qualitatively resembles Figure 3(b), for which ${\mathcal{R}}_0<1$ guarantees stability for the addiction-free equilibrium. This implies that reformed individuals helping at-risk individuals have the potential to significantly impact the long-term addicted population size, despite high relapse rates (${\mathcal{R}}_{\varphi }>1$).

Figure 4. Bifurcation diagrams for varying values of ν with parameters $\mu =0.00015$, $\gamma =0.0027$, $\text{β}=0.009$, $\varphi =0.004$, κ varies and ν as shown on each graph. Parameter values were taken from [26].

4. Numerical results

Besides the deterministic model (2(2) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mu N-\text{β}g\left(\tilde{S}\right)S\frac{A}{N}-\mu S,\\ \frac{\mathrm{d}A}{\mathrm{d}t}& =\text{β}g\left(\tilde{S}\right)S\frac{A}{N}+\varphi \tilde{S}\frac{A}{N}-(\mu +\gamma )A,\\ \frac{\mathrm{d}\tilde{S}}{\mathrm{d}t}& =\gamma A-\varphi \tilde{S}\frac{A}{N}-\mu \tilde{S},\end{array}$(2) ), we also consider a discrete stochastic model to compare the behaviour of both models and the dependence on the initial conditions. For the stochastic model, probability rates between states are given in Table 1, which are straightforwardly obtained from (2(2) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mu N-\text{β}g\left(\tilde{S}\right)S\frac{A}{N}-\mu S,\\ \frac{\mathrm{d}A}{\mathrm{d}t}& =\text{β}g\left(\tilde{S}\right)S\frac{A}{N}+\varphi \tilde{S}\frac{A}{N}-(\mu +\gamma )A,\\ \frac{\mathrm{d}\tilde{S}}{\mathrm{d}t}& =\gamma A-\varphi \tilde{S}\frac{A}{N}-\mu \tilde{S},\end{array}$(2) ). We only replace the recruitment rate by a constant Λ, chosen initially as $\mathrm{\Lambda }=N\left(0\right)\mu$, to essentially have a constant population. For simplicity, we take $t\in \{0,\mathrm{\Delta }t,2\mathrm{\Delta }t,3\mathrm{\Delta }t,\dots \}$ and the number of events on each time step $\mathrm{\Delta }t$ is assumed to follow a Poisson distribution with mean equal to the rates shown in Table 1.

Table 1. Transition rates for the stochastic model similar to (2(2) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mu N-\text{β}g\left(\tilde{S}\right)S\frac{A}{N}-\mu S,\\ \frac{\mathrm{d}A}{\mathrm{d}t}& =\text{β}g\left(\tilde{S}\right)S\frac{A}{N}+\varphi \tilde{S}\frac{A}{N}-(\mu +\gamma )A,\\ \frac{\mathrm{d}\tilde{S}}{\mathrm{d}t}& =\gamma A-\varphi \tilde{S}\frac{A}{N}-\mu \tilde{S},\end{array}$(2) )

Event	Effect	Rate
Addiction	$A\leftarrow A+1$, $S\leftarrow S-1$	$\text{β}g\left(\tilde{S}\right)SA/N$
Relapse	$A\leftarrow A+1$, $\tilde{S}\leftarrow S-1$	$\varphi A\tilde{S}/N$
Recovery	$\tilde{S}\leftarrow \tilde{S}+1$, $A\leftarrow A-1$	$\gamma A$
Recruitment	$S\leftarrow S+1$	Λ
Exit from S	$S\leftarrow S-1$	$\mu S$
Exit from A	$A\leftarrow A-1$	$\mu A$
Exit from $\tilde{S}$	$\tilde{S}\leftarrow \tilde{S}-1$	$\mu \tilde{S}$

Numerical simulations allow us to examine the addicted population dynamics over time for each region depicted in Figure 2; see Figures 5–8. We consider three different populations, labeled low ($N\approx 10,000$), medium ($N\approx 100,000$) and high ($N\approx 1,000,000$). Results for the stochastic model are scaled by $N\left(t\right)$, to compare results with the deterministic solution of the scaled model (3a(3a) $\begin{array}{rl}\frac{\mathrm{d}s}{\mathrm{d}t}& =\mu -\text{β}g\left(\tilde{s}\right)sa-\mu s,\end{array}$(3a) ). We include the deterministic solution $a\left(t\right)$ to compare the time series for both models.

Figure 5. Addicted population time series for Region 1 of the forward-backward bifurcation, with parameters $\mu =0.00015$, $\text{β}=0.009$, $\gamma =0.0027$, $\kappa =0.2$, $\nu =0.8$, $\varphi =0.0044$, $R_0=0.6315$, $R_{\varphi }=1.5439$, $I\left(0\right)/N\left(0\right)=0.15$. We present the mean $I\left(t\right)/N\left(t\right)$ for the stochastic model (100 simulations) for (a) low, (b) medium, and (c) high populations. Gray-shaded regions correspond to the 5th and 95th percentiles. Blue dots correspond to the deterministic solution.

We remark that the discrete stochastic model allows more realistic modeling of the inherent uncertainty in real-world situations and can provide insights on possible extreme events. Even though there are no significant differences in the long-term behaviour, a stochastic model can provide additional insights into the short-term dynamics of the system depending on the population size and initial conditions. For example, in the context of addiction modeling, there is often a high degree of variability in individual behaviour and response to treatment, and a stochastic model might capture this variability. Moreover, we observe cases where the addicted population is extinct in the stochastic model, even though this is not the case for the deterministic model.

Figure 5 (region 1) shows that the addicted population will, over time, decrease to the addicted-free equilibrium for both stochastic and deterministic models, independent of the total population size. We observe more variability among simulations with smaller populations, as expected.

Figure 6 (region 2) shows that the addicted population for the deterministic model will either establish itself or decrease to the addicted-free equilibrium, depending on the initial addicted population size. For the stochastic model, a similar behaviour occurs when the initial addicted population is small (1%). Nevertheless, if the total population is small enough, it is possible to obtain an addicted-free state. We observe a large variation among simulations when there is an initial addicted compartment with 10% of the total population in the low population case.

Figure 6. Addicted population time series for Region 2 of the forward-backward bifurcation, with parameters $\mu =0.00015$, $\text{β}=0.009$, $\gamma =0.0027$, $\kappa =0.3111$, $\nu =0.8$, and $\varphi =0.0044$, $R_0=0.9824$, $R_{\varphi }=1.5439$. We present the mean $I\left(t\right)/N\left(t\right)$ for the stochastic model (100 simulations) for (a) low, (b) medium, and (c) high populations, with $I\left(0\right)/N\left(0\right)=0.01$ (top) and $I\left(0\right)/N\left(0\right)=0.10$ (bottom). Gray-shaded regions correspond to the 5th and 95th percentiles. Blue lines correspond to the deterministic solution.

Figure 7 (region 3) shows that the addicted population will establish itself at either a large or small size, once again depending on the initial addicted population size. Figures 6 and 7 illustrate sensitivity to initial conditions for the deterministic model; whether the addicted community is established, and how large the community is, is dependent on how pervasive the initial population of addicted individuals is. For medium and large populations, the stochastic models preserve the qualitative behaviour of the deterministic curves. Nevertheless, in the low population cases, there are cases where an addicted-free state is reached.

Figure 7. Addicted population time series for Region 3 of the forward-backward bifurcation, with parameters $\mu =0.00015$, $\text{β}=0.009$, $\gamma =0.0027$, $\kappa =0.3243$, $\nu =0.8$, $\varphi =0.0042$, $R_0=1.0241$, $R_{\varphi }=1.4736$. We present the mean $I\left(t\right)/N\left(t\right)$ for the stochastic model (100 simulations) for (a) low, (b) medium, and (c) high populations, with $I\left(0\right)/N\left(0\right)=0.01$ (top) and $I\left(0\right)/N\left(0\right)=0.10$ (bottom). Gray-shaded regions correspond to the 5th and 95th percentiles. Blue lines correspond to the deterministic solution.

Figure 8. Addicted population time series for Region 4 of the forward-backward bifurcation, with parameters $\mu =0.00015$, $\text{β}=0.009$, $\gamma =0.0027$, $\kappa =0.3333$, $\nu =0.8$, $\varphi =0.0044$, $R_0=1.0525$, $R_{\varphi }=1.5439$, $I\left(0\right)/N\left(0\right)=0.01$. We present the mean $I\left(t\right)/N\left(t\right)$ for the stochastic model (100 simulations) for (a) low, (b) medium, and (c) high populations. Gray-shaded regions correspond to the 5th and 95th percentiles. Blue dots correspond to the deterministic solution.

Finally, Figure 8 (region 4) shows that the addicted population will establish itself at a (relatively) large population size despite a very small initial addicted population for the deterministic model. This occurs largely because the effects of relapse outweigh the attempts to discourage addicted involvement. Again, an addicted-free state is possible if the total population is small enough, even though the mean of 100 simulations is close to the deterministic curve.

5. Discussion

We explored an epidemic-type model that includes nonlinear relapse in the temporarily recovered population. Results showed high sensitivity to the initial addicted population. If the initial addicted population in an at-risk environment is large enough, in other words, already established, the addicted population is more likely to establish itself at an endemic equilibrium. Hence, well-established addicted populations play prominent social, political, and economic roles in the community and are thus much more difficult to control [21,25].

Our model indicates that reformed individuals play a crucial role in addicted population dynamics when relapse rates are low. A large proportion of reformed individuals serving as mentors has the potential to significantly reduce the addicted population [19,22], given that relapse rates remain under control. A lack of opportunities could lead reformed individuals to return to addiction. Rehabilitation programs that aim to reintegrate reformed individuals into productive individuals in society must offer consistent supervised rehabilitation.

While focusing on reformed individuals is essential in reducing the addicted population, recruitment into the addicted population is a major player in reducing the addicted population. The cost of becoming addicted plays a role in the basic reproductive number, ${\mathcal{R}}_0$, which implies that changes in κ have significant implications in the transmission dynamics of the system. Furthermore, informing susceptible individuals about other lifestyles and opportunities, such as education, can discourage individuals from getting involved in risky environments where they can ultimately be pulled into addiction.

Moreover, if the social influence of the reformed individuals on the at-risk susceptible population is strong, the long-term addicted population becomes manageable. For a specific cost (κ), a small region exists where $1<{\mathcal{R}}_0<{\mathcal{R}}_0^{\ast }$ and ${\mathcal{R}}_{\varphi }>1$ with multiple stable addicted populations that are highly dependent on the initial addicted population size. If the initial addicted population in an at-risk environment is large enough, in other words, already established, the addicted population establishes itself at the higher endemic equilibrium. Our model shows the influence of established problem communities, highlighting the importance of prevention programs and relapse rates.

While reformed individuals may impact addicted population dynamics, other factors, such as cost and relapse rate, play a role in the effectiveness of reformed individuals in population control. When relapse rates are low, reformed individuals are crucial to addicted population dynamics. A high value of ν can shift the forward-backward bifurcation to the point where ${\mathcal{R}}_0<1$ produces an addicted-free equilibrium. This highlights the importance of controlling relapse rates and encouraging reformed individuals to become involved with addiction prevention programs.

From our model, we also found that the cost of addiction significantly impacts the addicted population dynamics. Cost is a factor in the basic reproductive number ${\mathcal{R}}_0$, which means that changes in this value may have significant implications for the addicted population. If the costs are low to get into addiction, there is little that reformed individuals can do to decrease the growth in the addicted population. An alternative to lowering the cost of addiction is to educate at-risk individuals about the costs of addiction. This, in turn, may encourage individuals to look at these costs as a deterrent, which can ultimately help decrease the addicted population. Finally, more considerable efforts are needed not only to encourage reformed individuals to mentor individuals in an at-risk environment but also to reduce the relapse rate and help to educate the at-risk population to help contain individuals from getting into addiction.

Furthermore, incorporating more complex dynamics into the model, such as the effect of stigma, the impact of policy interventions, and the role of social networks in the spread of addiction, is of interest to public health officials. Additionally, there is potential to use addiction models to study the impact of environmental factors, such as climate change and economic instability, on substance abuse patterns. Another area of future research is developing models that consider individual differences in susceptibility to addiction and response to treatment. Overall, future work in addiction modeling can improve our understanding of the underlying mechanisms of substance abuse and inform the development of effective prevention and treatment strategies.

Acknowledgments

The authors would like to thank the support from the Research Center in Pure and Applied Mathematics and the Department of Mathematics at Universidad de Costa Rica.

Disclosure statement

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