Introduction

Drug abuse, encompassing illicit substances, misused prescriptions, and even over-the-counter medications, poses a significant global threat1. According to the United Nations Office on Drugs and Crime (UNODC) World Drug Report \(2023\), the issue is widespread, with an estimated \(296\) million people (\(5.8\%\) aged 15–64) using drugs in \(2021\). While cannabis is the most common drug, opioids pose the greatest threat, linked to overdoses and health problems often associated with heroin and fentanyl, while also addressing organized crime, corruption, and terrorism2. Furthermore, crime transcends borders, affecting every nation on the world stage. It leaves a deep scar on individuals, communities, and entire economies. The human cost is immense, inflicting physical and emotional trauma on victims. Financially, crime leads to devastating losses, hindering economic growth3. Furthermore, it erodes trust in the justice system and weakens the social fabric that binds communities together4.

The intricate link between drugs and crime is undeniable. While drug use can fuel criminal activity to finance the habit, social and economic hardships can also be root causes. The type of drug, frequency of use, and an individual’s mental health further complicate this relationship5. This drug-crime nexus is constantly evolving as drug use patterns, law enforcement approaches, and societal views shift6,7. Understanding these dynamics is critical for crafting effective solutions. Research suggests a clear correlation: drug users are more likely to commit crimes, arrested individuals are often under the influence when committing offenses, and drug trafficking breeds violence8,9.

A multifaceted strategy is crucial to tackling the intertwined issues of crime and drug abuse. Educational programs and community initiatives informed10 can prevent drug use by educating individuals about the dangers and promoting healthy choices. Community-based programs fostering collaboration between residents and various stakeholders11 build resilience within communities, further deterring crime and drug abuse. Law enforcement plays a vital role through proactive measures like community policing and disrupting drug trafficking networks, as outlined by the United Nations Office on Drugs and Crime12. Accessible and effective treatment programs are essential for those struggling with addiction, as research by13,14 demonstrates the success of various recovery approaches.

The world around us is brimming with complex systems and social issues like crime and drug abuse are prime examples. To understand these intricate problems, scientists rely on a powerful tool: mathematical modeling. These models act like virtual laboratories, allowing researchers to simulate real-life scenarios and analyze the dynamic interactions between different factors15. This approach is not limited to social sciences; mathematical modeling plays a crucial role in simulating and analyzing phenomena across various scientific disciplines, including economics, physics, and biological sciences16,17,18,19. In the realm of social issues specifically, researchers have employed mathematical models to tackle problems like crime dynamics, exploring the impact of interventions such as age-specific programs or rehabilitation measures20,21,22. Similarly, models have been developed to understand drug addiction and its control, illicit drug use patterns, and even analyze the cost-effectiveness of different intervention strategies23,24,25. Previous research has explored the relationship between substance abuse and crime, often focusing on isolated aspects of the problem. For instance, studies in South Africa have examined the correlation between substance abuse and related crime26,27, while others have investigated the link between drug abuse and violence using optimal control methods28. While these studies provide valuable insights, they frequently overlook the complex co-dynamics between crime and substance abuse, neglecting crucial elements such as crime dynamics, co-engaged prone populations, and cross-recidivism. Mathematical modeling offers a promising avenue for understanding these intricate relationships and informing the development of effective prevention and intervention strategies29. In the subsequent sections, we will delve into a model that explicitly addresses these limitations and captures the co-evolution of crime and substance abuse.

In the article, researchers explore a new mathematical model to understand the link between crime and drug abuse. This model treats both behaviors as competing forces within a system with predictable rules. The model considers rehab as a way to influence these dynamics, potentially reducing both crime and drug use. Another factor, Confine, might represent incarceration or other methods to curb criminal activity. The model examines how repeated criminal behavior and drug use become ingrained habits, potentially leading to a point of no return (entrenched threshold). By analyzing how different factors influence this threshold and potential tipping points (criticality), researchers hope to develop better prevention strategies. This model, while promising, could benefit from real-world data validation, the inclusion of additional social and economic factors, and potentially even a shift to a model that incorporates randomness to account for unforeseen circumstances.

The full research paper delves deeper into the specifics of the crime-drug abuse model. “Model development” section lays the groundwork by explaining the model’s assumptions, how it’s mathematically formulated, and its basic behavior. “Results and discussion” section then results and discussion, analyzes the stability of two key situations: a state with no crime or drug abuse, and a state where both are prevalent. Additionally, it explores how a simpler model with just one behavior (crime or drug abuse) reacts to changes in treatment effectiveness. Moving on, “Numerical results” section uses simulations to assess the potential impact of interventions like increased access to rehab or stricter law enforcement. Finally, “Conclusion” section wraps things up with the main conclusions about the model’s ability to represent the dynamics of crime and substance abuse. The researchers acknowledge any limitations and suggest areas for further exploration, potentially including incorporating more factors or making the model more applicable to real-world scenarios.

Model development

The model assumes a well-mixed population where everyone has an equal chance of interacting with each other. To monitor criminal or substance use behavior among the population of size \(N(t)\), we divide it into thirteen exclusive compartments. The categories include the Susceptible \((S)\) class, who are vulnerable to both crime and substance abuse. There is also the Crime-Prone \((P_c)\) class, which consists of individuals at a high risk of committing crimes. The substance Abuse-Prone \((P_a)\) class comprises people at a high risk of developing substance abuse. The Co-occurrence-Prone \((P_e)\) class has individuals with a heightened risk of developing both crime and substance abuse if exposed to either. The Criminal \((I_c)\) class comprises individuals actively committing crimes, while the substance Abuse \((I_a)\) class consists of individuals with active substance abuse. The Co-occurred \((I_e)\) class comprises individuals involved in both crime and substance abuse. Furthermore, we consider individuals for treatments \((T_c, T_a)\) for crime and substance abuse, respectively, and co-occurring conditions \((T_e)\) to prevent further diffusion. Finally, the recovery category encompasses individuals who have successfully ceased criminal activity \((R_c)\), substance abuse \((R_a)\), and engagement in co-occurring behaviors \((R_e)\).

Main assumptions

  1. (1)

    Population Dynamics The model reflects real-world changes in the population by including new individuals (births) and those exiting the system through recovery or death.

  2. (2)

    Force of initiation We assume that individuals actively involved in both crime and drug abuse (\(I_e\)) exert an equivalent influence on others to engage in either crime (\(I_c\)) or drug abuse (\(I_a\)) as those solely involved in these respective behaviors. However, individuals undergoing treatment (\(T_j\)) are considered to have a reduced likelihood of influencing others to adopt these behaviors due to their engagement in rehabilitation programs and potential changes in social circles.

  3. (3)

    Co-occurrence An important aspect of the model is co-occurrence, which refers to the potential link between criminal activity and substance abuse. Engaging in one behavior increases the likelihood of engaging in the other. This creates a high-risk group (denoted by \(P_e\)) that is prone to experiencing both simultaneously. Individuals may move into this group due to substance use while committing crimes, or vice versa. However, for simplicity, the model assumes that the chance of starting both behaviors at the same time is negligible.

  4. (4)

    Recovery Dynamics The model acknowledges the complexities of recovery, especially for individuals struggling with both crime and substance abuse (co-occurrence). An individual in the co-occurrence class \((I_e)\) and \(T_e\) might recover from one behavior while continuing the other, leading to a transition to either the substance-only \((I_a, T_a)\) or crime-only \((I_c, T_c)\) class, respectively. However, the model assumes that once an individual has fully recovered from a specific behavior, they gain permanent immunity from relapsing into that same behavior, reflecting lessons learned from rehabilitation experiences.

  5. (5)

    Intervention Strategies The model explores the impact of intervention strategies on individuals engaged in criminal activities or substance abuse. A crucial assumption is adherence to these interventions, which can significantly influence their effectiveness. Two primary types of intervention are considered:

    1. (a)

      Outpatient Treatment This approach involves providing treatment to individuals while they continue to live within their community. Outpatient programs typically include counseling, therapy, education, and support groups. This type of treatment is less disruptive to an individual’s daily life but may require a higher level of self-discipline and support.

    2. (b)

      Inpatient Treatment This strategy involves isolating individuals who pose a high risk to themselves or others. For criminal activities, this could be represented by incarceration, while for substance abuse, it could involve residential treatment programs, detoxification, or court-mandated treatment centers. Inpatient treatment provides a controlled environment aimed at reducing exposure to negative influences and facilitating more intensive recovery efforts.

    The model uses the parameter \(\omega\) to measure the effectiveness of these interventions, with adjustments made for potential reductions in impact due to adherence challenges.

  6. (6)

    Cross-Recidivism or Relapse The model incorporates the concept of cross-recidivism or relapse, recognizing that individuals who have recovered from one behavior may still be vulnerable to developing another. Drawing on the idea of “cross-immunity” from previous research30,31, the model introduces parameters \(\tau _c\) and \(\tau _a\) (ranging from 0 to 1) to capture the degree of this susceptibility. However, consistent with the assumption of permanent immunity, the model posits that successful recovery from a specific behavior provides lasting resilience, preventing relapse into the same behavior.

The model incorporates a constant inflow of new individuals \(\Lambda\) susceptible to crime and substance abuse \((S)\). All individuals are subject to a natural death rate \((\mu)\). Susceptible individuals can transition to a pre-engaged state influenced by exposure to criminal or substance-using behaviors. This results in three categories: crime-prone (\(P_c\)), substance abuse-prone (\(P_a\)), and prone to both (\(P_e\)). Individuals in these pre-engaged states may progress to actively engaging in crime (\(I_c\)), substance abuse (\(I_a\)), or both (\(I_e\)) at rates \(\sigma _c\), \(\sigma _a\), and \(\sigma _e\), respectively. Those actively engaged can enter treatment programs for crime (\(T_c\)), substance abuse (\(T_a\)), or co-occurring issues (\(T_e\)) at rates \(\theta _c, \theta _a\), and \(\theta _e\). However, treatment programs are not entirely effective, and individuals may relapse or re-offend (\(\kappa _c, \kappa _a\), and \(\kappa _e\)) back into their respective active states. Ultimately, individuals can recover from crime (\(R_c\)), substance abuse (\(R_a\)), or both (\(R_e\)) at rates \(\gamma _c, \gamma _a\), and \(\gamma _e\), respectively. Additionally, individuals in treatment can recover and exit the system at rates \(\phi _a, \phi _c\), and \(\phi _e\). A visual representation of these dynamics is presented in Fig. 1. The exact formulas for the force of initiation of crime \((\xi _c)\) and substance abuse \((\xi _a)\) are explained in the following equations:

$$\begin{aligned} \xi _c=\dfrac{\beta _c\left( I_c+I_e+\omega (T_c+T_e)\right) }{N},\quad \xi _a =\dfrac{\beta _a\left( I_a+I_e+\omega (T_a+ T_e)\right) }{N}. \end{aligned}$$
(1)

Within the equations we just discussed (1), \(\beta _c\) and \(\beta _a\) represent the rates at which crime and substance abuse spread through the population, respectively. These rates were derived considering the model’s assumptions (mentioned earlier), the visual representation in Fig. 1, and the definitions of the parameters in Table  1. This combination allows us to establish the following system of equations:

Fig. 1
figure 1

The schematic diagram of crime and substance abuse co-dynamics.

Full size image
$$\begin{aligned} \begin{aligned} \frac{dS}{dt}&=\Lambda -\xi _c S-\xi _aS-\mu S,\\ \frac{dP_c}{dt}&=\xi _c (S+\tau _aR_a)-(\sigma _c+\mu ) P_c,\\ \frac{dP_a}{dt}&=\xi _a(S+\tau _cR_c)-(\sigma _a+\mu )P_a,\\ \frac{dI_c}{dt}&=\sigma _cP_c+\gamma _aI_e+\kappa _c T_c-(1-\gamma _a)\xi _a I_c-(\theta _c+\gamma _c +\delta _c+\mu ) I_c,\\ \frac{dI_a}{dt}&=\sigma _aP_a+\gamma _cI_e+\kappa _a T_a-(1-\gamma _c)\xi _cI_a-(\theta _a+\gamma _a+\delta _a+\mu ) I_a,\\ \frac{dP_e}{dt}&= (1-\gamma _a) \xi _a I_c+(1-\gamma _c)\xi _c I_a -(\sigma _e+\mu )P_e,\\ \frac{dI_e}{dt}&=\sigma _eP_e+\kappa _eT_e-(\theta _e+\gamma _c+\gamma _a+\gamma _e+\delta _e+\mu )I_e,\\ \frac{dT_{c}}{dt}&=\theta _cI_c+\phi _a T_e-(\kappa _{c}+\phi _c +\delta _{c}+\mu ) T_c,\\ \frac{dT_a}{dt}&=\theta _aI_a+\phi _c T_e-(\kappa _a+\phi _a +\delta _a+\mu )T_a,\\ \frac{dT_e}{dt}&=\theta _eI_e-(\kappa _e+\phi _a+\phi _c+\phi _e+\delta _e+\mu )T_e,\\ \frac{dR_{c}}{dt}&=\gamma _c I_c+\phi _cT_c-\xi _a\tau _cR_c-\mu R_c,\\ \frac{dR_a}{dt}&=\gamma _a I_a+\phi _a T_a-\xi _c\tau _a R_a-\mu R_a,\\ \frac{dR_e}{dt}&= \gamma _e I_e+\phi _e T_e-\mu R_e. \end{aligned} \end{aligned}$$
(2)

At a given time t, the total population size \(N\) is given by

$$\begin{aligned} N(t)= S(t) + P_c(t) + P_a(t) +P_e(t) + I_c(t) + I_a(t) + I_e(t)+ T_c(t) + T_a(t) +T_e(t)+ R_c(t) +R_a(t) +R_e(t). \end{aligned}$$

Additionally, the system (2) requires that all the starting values (initial data) for each population class must be non-negative. This is shown in the following equation.

$$\begin{aligned} \begin{aligned}&S(0)\ge 0,~P_c(0)\ge 0, ~ P_a(0)\ge 0, ~P_e(0)\ge 0,\\&I_c(0)\ge 0, ~I_a(0)\ge 0,~ I_e(0)\ge 0, ~T_c(0)\ge 0, ~T_a(0)\ge 0, \\&T_e(0)\ge 0, ~R_c(0)\ge 0,~ R_a(0)\ge 0, \text { and } R_e(0)\ge 0. \end{aligned} \end{aligned}$$
(3)
Table 1 Model parameters and their descriptions.
Full size table

Basic properties of the main model

Because we are dealing with problems related to population dynamics, for meaningful biological interpretation, all the variables must be positive and bounded for all time \(t\ge 0\).

Lemma 2.1

(Positivity of solutions) Let\(( S(t),P_c(t),P_a(t), P_e(t)\),\(I_c(t), I_a(t), I_e(t)\),\(T_c(t), T_a(t), T_e(t)\),\(R_c(t),R_a(t), R_e(t))\)be the solution of the system (2). With initial conditions (3), the solution exists and remains non-negative for all\(t \ge 0\).

Proof

Consider system (2), it follows from (3) that

$$\begin{aligned} \frac{dS}{dt}|_{S=0}&=\Lambda >0, \qquad \frac{dP_c}{dt}|_{P_{c}=0}=\xi _c S+\xi _c\tau _aT_a\ge 0,\quad \forall ~ S,T_a\ge 0, \\ \frac{dP_a}{dt}|_{P_a=0}&=\xi _aS+\xi _a\tau _cT_c\ge 0, \quad \forall ~ S,T_c\ge 0,\\ \frac{dI_c}{dt}|_{I_c=0}&=\sigma _cP_c+\gamma _aI_e+\kappa _c T_c\ge 0,\quad \forall ~ P_c,I_e, T_c\ge 0,\\ \frac{dI_a}{dt}|_{I_a=0}&=\sigma _aP_a+\gamma _cI_e+\kappa _a T_a\ge 0,\quad \forall ~P_a, I_e, T_a\ge 0,\\ \frac{dP_e}{dt}|_{P_e=0}&= (1-\gamma _a) \xi _a I_c+(1-\gamma _c)\xi _c I_a \ge 0,\quad \forall ~I_c, I_a,\ge 0, \, \text {and},\,\, 0\le \gamma _c,\gamma _a \le 1\\ \frac{dI_e}{dt}|_{I_e=0}&=\sigma _eP_e+\kappa _eT_e\ge 0, \forall ~ P_e, T_e\ge 0,\\ \frac{dT_{c}}{dt}|_{T_{c}=0}&=\theta _cI_c+\phi _a T_e\ge 0 ,\quad \forall ~ I_c, T_e\ge 0,\quad \frac{dT_a}{dt}|_{T_a=0}=\theta _aI_a+\phi _c T_e\ge 0,\quad \forall ~ I_a, T_e\ge 0,\\ \frac{dT_e}{dt}|_{T_e=0}&= \theta _eI_e\ge 0, \forall ~ I_e\ge 0,\quad \frac{dR_{c}}{dt}|_{R_{c}=0}=\gamma _cI_c+\phi _cT_c\ge 0, \forall ~ I_c, T_c\ge 0,\\ \frac{dR_a}{dt}|_{R_a=0}&=\gamma _aI_a+\phi _aT_a\ge 0, \forall ~I_a, T_a\ge 0,\, \frac{dR_e}{dt}|_{R_e=0}= \gamma _eI_e+\phi _eT_e\ge 0, \forall ~ I_e, T_e\ge 0. \end{aligned}$$

We now prove that all solutions of (2) are non-negative as \(t\rightarrow \infty\). \(\square\)

Lemma 2.2

(Boundedness of solutions) The closed set

$$\begin{aligned} \Omega&= \Big \{ (S,P_c,P_a, P_e,I_c,I_a, I_e,T_c,T_a, T_e,R_c,R_a, R_e)\in \mathbb {R}^{13}_+: S+P_c+P_a+ P_e+I_c+I_a\\&\quad +I_e+T_c+T_a+ T_e+R_c+R_a+R_e\le \dfrac{\Lambda }{\mu }\Big \} \end{aligned}$$

is positively-invariant for the model (2).

Proof

Adding all equations in the model (2) gives the rate of change in a total population over time:

$$\begin{aligned} \begin{aligned} N^\prime (t)&= \Lambda - \mu N - (\delta _c I_c + \delta _c T_c +\delta _a I_a + \delta _aT_a +\delta _eI_e + \delta _eT_e).\\&\le \Lambda - \mu N. \end{aligned} \end{aligned}$$
(4)

Thus \(N^\prime (t) + \mu N\le \Lambda .\)

Now applying the theorem of differential inequalities43, we obtain

$$\begin{aligned}0\le N(t)\le \frac{\Lambda }{\mu }+\left( N(0)-\frac{\Lambda }{\mu }\right) e^{-\mu t}.\end{aligned}$$

As \(t\rightarrow \infty\), we have

$$\begin{aligned}0\le N(t)\le \Lambda /\mu .\end{aligned}$$

Hence, the region \(\Omega\) is positive invariant and attracts all solutions in \(\mathbb {R}^{13}_+\). \(\square\)

Results and discussion

This section analyzes the mathematical model for the co-occurrence of crime and substance abuse presented in system  2.

Stability of crime and substance use-free equilibrium

The model admits a crime- and substance abuse-free equilibrium \(\textbf{E}^0\), where no individuals are involved in either behavior. This equilibrium is given by: \(\textbf{E}^0 = (S^0, P_c^0, P_a^0, P_{e}^0\), \(I_c^0, I_a^0, I_{e}^0, T_c^0, T_a^0, T_{e}^0\), \(R_c^0, R_a^0, R_{e}^0) = \left( \dfrac{\Lambda }{\mu },0,0,0,0,0,0,0,0,0,0,0,0\right)\). Before stability of \(\textbf{E}^0\), we employ the next-generation matrix method44 to compute the control reproduction number. The next-generation matrices \(F\) and \(V\) for the model are defined as follows:

$$\begin{aligned}F=\left( \begin{array}{lllllllll} 0 & 0 & 0 & \beta _c & 0 & \beta _c & \omega \beta _c & 0 & \omega \beta _c \\ 0 & 0 & 0 & 0 & \beta _a & \beta _a & 0 & \omega \beta _a & \omega \beta _a \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \end{array} \right) , \quad V=\left( \begin{array}{lllllllll} \textbf{g}_c & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & \textbf{g}_a & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & \textbf{g}_e & 0 & 0 & 0 & 0 & 0 & 0 \\ -\sigma _c & 0 & 0 & \textbf{b}_c & 0 & -\gamma _a & -\kappa _c & 0 & 0 \\ 0 & -\sigma _a & 0 & 0 & \textbf{b}_a & -\gamma _c & 0 & -\kappa _a & 0 \\ 0 & 0 & -\sigma _e & 0 & 0 & \textbf{b}_e & 0 & 0 & -\kappa _e \\ 0 & 0 & 0 & -\theta _c & 0 & 0 & \textbf{d}_c & 0 & 0 \\ 0 & 0 & 0 & 0 & -\theta _a & 0 & 0 & \textbf{d}_a & 0 \\ 0 & 0 & 0 & 0 & 0 & -\theta _e & 0 & 0 & \textbf{d}_e \\ \end{array} \right) ,\end{aligned}$$

where \(g_i=\sigma _i+\mu\), \(b_i=\theta _i+\gamma _i+\delta _i+\mu\), \(d_i=\kappa _i+\phi _i+\delta _i+\mu\), for \(i=c,a\), \(g_e=\sigma _e+\mu\), \(b_e=\theta _e+\gamma _e+\delta _e+\mu +\gamma _c+\gamma _a\), and \(d_e=\kappa _e+\phi _e+\delta _e+\mu +\phi _c+\phi _a\).

The spectral radius \(\rho (FV^{-1})\) is the maximum of \(\mathscr {R}_c\) and \(\mathscr {R}_a\), where:

$$\begin{aligned} \mathscr {R}_c= \dfrac{\beta _c \sigma _c\textbf{d}_c+\omega \beta _c \sigma _c \theta _c}{\textbf{g}_c(\textbf{r}_c\textbf{d}_c+\textbf{q}_c\theta _c)},\qquad \mathscr {R}_a=\dfrac{\beta _a \sigma _a\textbf{d}_a+\omega \beta _a \sigma _a \theta _a}{\textbf{g}_a(\textbf{r}_a\textbf{d}_a+\textbf{q}_a\theta _a)}, \end{aligned}$$

where \(\textbf{r}_i=\gamma _i+\delta _i+\mu ,\quad \textbf{q}_i=\phi _i+\delta _i+\mu\) for \(i=c,a\). The co-occurrence reproduction number \({\mathscr {R}_{T}}\) is defined as the maximum of \({\mathscr {R}_{c}}\) and \({\mathscr {R}_a}\):

$$\begin{aligned} {\mathscr {R}_{T}} =\max \{{\mathscr {R}_c}, {\mathscr {R}_a}\}. \end{aligned}$$
(5)

Noted that, \(\mathscr {R}_c\) represents the crime-based reproduction number, and \(\mathscr {R}_a\) represents the substance abuse-based reproduction number.

Theorem 3.1

If\(\{{\mathscr {R}_c}<1\)and\({\mathscr {R}_a}\}<1\), then\(\textbf{E}^0\)is locally asymptotically stable (LAS). Otherwise, it is unstable.

Proof

The Jacobin matrix of the system (2) at the equilibrium \(\textbf{E}^0\) is

$$\begin{aligned}J(E^0)=\left( \begin{array}{lllllllllllll} -\mu & 0 & 0 & -\beta _c & -\beta _a & 0 & -\beta _a-\beta _c & -\beta _c \omega & -\beta _a \omega & -\beta _a \omega -\beta _c \omega & 0 & 0 & 0 \\ 0 & -\textbf{g}_c & 0 & \beta _c & 0 & 0 & \beta _c & \beta _c \omega & 0 & \beta _c \omega & 0 & 0 & 0 \\ 0 & 0 & -\textbf{g}_a & 0 & \beta _a & 0 & \beta _a & 0 & \beta _a \omega & \beta _a \omega & 0 & 0 & 0 \\ 0 & \sigma _c & 0 & -\textbf{b}_c & 0 & 0 & \gamma _a & \kappa _c & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & \sigma _a & 0 & -\textbf{b}_a & 0 & \gamma _c & 0 & \kappa _a & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & -\textbf{g}_e & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & \sigma _e & -\textbf{b}_e & 0 & 0 & \kappa _e & 0 & 0 & 0 \\ 0 & 0 & 0 & \theta _c & 0 & 0 & 0 & -\textbf{d}_c & 0 & \phi _a & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & \theta _a & 0 & 0 & 0 & -\textbf{d}_a & \phi _c & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & \theta _e & 0 & 0 & -\textbf{d}_e & 0 & 0 & 0 \\ 0 & 0 & 0 & \gamma _c & 0 & 0 & 0 & \phi _c & 0 & 0 & -\mu & 0 & 0 \\ 0 & 0 & 0 & 0 & \gamma _a & 0 & 0 & 0 & \phi _a & 0 & 0 & -\mu & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & \gamma _e & 0 & 0 & \phi _e & 0 & 0 & -\mu \\ \end{array} \right) .\end{aligned}$$

The eigenvalues of \(J(E^0)\) are \(\lambda _{1,2,3,4}=-\mu ,\, \lambda _5=-g_e\), and the remaining eigenvalues are determined by the following characteristic polynomials:

$$\begin{aligned} P_e(\lambda )&=\lambda ^2+(\textbf{b}_e+\textbf{d}_e)\lambda +\textbf{b}_e \textbf{d}_e-\theta _e \kappa _e,\\ P_a(\lambda )&=\lambda ^3+ (\textbf{b}_a+\textbf{d}_a+\textbf{g}_a)\lambda ^2 + (\textbf{r}_a\textbf{d}_a+\textbf{q}_a\theta _a+(\textbf{b}_a +\textbf{d}_a) \textbf{g}_a -\beta _a \sigma _a)\lambda +\textbf{g}_a(\textbf{r}_a\textbf{d}_a+\textbf{q}_a\theta _a)(1-\mathscr {R}_a),\\ P_c(\lambda )&=\lambda ^3+(\textbf{b}_c+\textbf{d}_c+\textbf{g}_c)\lambda ^2 +(\textbf{r}_c\textbf{d}_c+\textbf{q}_c\theta _c+(\textbf{b}_c +\textbf{d}_c) \textbf{g}_c-\beta _c \sigma _c)\lambda + \textbf{g}_c(\textbf{r}_c\textbf{d}_c+\textbf{q}_c\theta _c)(1-\mathscr {R}_c), \end{aligned}$$

The eigenvalues of \(P_e(\lambda )\) have negative real parts. By the Routh–Hurwitz criterion, the eigenvalues of \(P_c(\lambda )\) and \(P_a(\lambda )\) have negative real parts if \(\mathscr {R}_c<1\) and \(\mathscr {R}_a<1\), respectively. \(\square\)

Global stability of crime-substance abuse free equilibrium

To investigate the global stability of the crime and substance abuse co-occurrence model using Castillo-Chavez45, we will reformulate the system presented in the model (2) as follows: we denote a vector of non-engaged individuals by \(\textbf{x}=(S, R_c, R_a, R_e)^T\in \mathbb {R}^4\) and that of engaged individuals by \(\textbf{I}=(P_c, P_a, P_e, I_c, I_a, I_e, T_c, T_a, T_e)^T\in \mathbb {R}^9\) such that

$$\begin{aligned} \begin{aligned} \dfrac{d\textbf{x}}{dt}=F(\textbf{x}, \textbf{I}),~~ \dfrac{d\textbf{I}}{dt}=G(\textbf{x}, \textbf{I}). \end{aligned} \end{aligned}$$
(6)

We also denote the crime and substance abuse free equilibrium, \({\textbf{E}}^0=(\textbf{x}^*,0,0,0,0,0,0,0,0,0)\) where \(\textbf{x}^*\) define the crime and substance abuse free equilibrium of the system \({d\textbf{x}}/{dt}\). Moreover, we state the following conditions:

H1

For \(F(\textbf{x}, \textbf{I})|_{\textbf{x}^*}\), \(\textbf{x}^*\) is globally asymptotically stable GAS,

H2

\(G(\textbf{x}, \textbf{I})=A\textbf{I}-\hat{G}(\textbf{x}, \textbf{I})\), \(\hat{G}(\textbf{x}, \textbf{I})\ge 0\) for \((\textbf{x}, \textbf{I})\in \Omega\),

where \(A=G(\textbf{x}^*, 0)\) is an M-matrix. The following claim is true when the given two conditions hold for the system (2):

Theorem 3.2

The crime and substance abuse free equilibrium\({\textbf{E}}^0=(\textbf{x}^*,0,0,0,0,0,0,0,0,0)\)of system (6), equivalent to (2) is globally asymptotically stable (GAS) if assumptionsH1andH2holds.

Proof

Let us now consider the prove for assumptions H1 and H2 using Carlos Castillo-Chavez45. From (6), we have

$$\begin{aligned} \begin{array}{l} F(\textbf{x}, 0)= \begin{pmatrix} \Lambda -\mu S\\ 0 \\ 0 \end{pmatrix},\, \hat{G}(\textbf{x}, \textbf{I})=\begin{pmatrix} \hat{G}_{1c}\\ \hat{G}_{2a}\\ \hat{G}_{3c}+\hat{G}_{3a}\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0 \end{pmatrix}, \, A= \begin{pmatrix} -\textbf{g}_c & 0 & 0 & \beta _c & 0 & \beta _c & \beta _c & 0 & \beta _c \\ 0 & -\textbf{g}_a & 0 & 0 & \beta _a & \beta _a & 0 & \beta _a & \beta _a \\ 0 & 0 & -\textbf{g}_e & \beta _c & \beta _a & \beta _c+\beta _a & \beta _c & \beta _a & \beta _c+\beta _a \\ \sigma _c & 0 & 0 & -\textbf{b}_c & 0 & \gamma _a & \kappa _c & 0 & 0 \\ 0 & \sigma _a & 0 & 0 & -\textbf{b}_a & \gamma _c & 0 & \kappa _a & 0 \\ 0 & 0 & \sigma _e & 0 & 0 & -\textbf{b}_e & 0 & 0 & \kappa _e \\ 0 & 0 & 0 & \theta _c & 0 & 0 & -\textbf{d}_c & 0 & 0 \\ 0 & 0 & 0 & 0 & \theta _a & 0 & 0 & -\textbf{d}_a & 0 \\ 0 & 0 & 0 & 0 & 0 & \theta _e & 0 & 0 & -\textbf{d}_e \\ \end{pmatrix} . \end{array} \end{aligned}$$

where

$$\begin{aligned} \hat{G}_{1c}&=\beta _c\left[ \left( 1-\frac{(S+\tau _a R_a)}{N}\right) (I_c+I_e)+\left( 1-\frac{\omega (S+\tau _a R_a)}{N}\right) (T_c+T_e)\right] ,\\ \hat{G}_{2a}&=\beta _a\left[ \left( 1-\frac{(S+\tau _c R_c)}{N}\right) (I_a+I_e)+\left( 1-\frac{\omega (S+\tau _c R_c)}{N}\right) (T_a+T_e)\right] ,\\ \hat{G}_{3c}&=\beta _c\left[ \left( 1-\frac{(1-\gamma _c)I_a}{N}\right) \left( I_c+I_e\right) +\left( 1-\frac{(1-\gamma _c)\omega I_a}{N}\right) \left( T_c+T_e\right) \right] ,\\ \hat{G}_{3a}&=\beta _a\left[ \left( 1-\frac{(1-\gamma _a)I_c}{N}\right) \left( I_a+I_e\right) +\left( 1-\frac{(1-\gamma _a)\omega I_c}{N}\right) \left( T_a+T_e\right) \right] . \end{aligned}$$

It follows that \(\omega\),\(\gamma _c\),\(\gamma _a\),\(\tau _c\),\(\tau _a \in (0,1]\), and \(0\le (S+\tau _{c} R_{c})\),\(\omega (S+\tau _c R_c)\),\(\le S+R_c\le N\), \(0\le (S+\tau _{a} R_{a})\),\(\omega (S+\tau _a R_a)\),\(\le S+R_a\le N\), \(0\le (1-\gamma _a)I_c\),\(\omega (1-\gamma _a)I_c\le I_c\le N\), \(0\le (1-\gamma _c)I_a\),\(\omega (1-\gamma _c)I_a\le I_a\le N\) then \(\hat{G}(\textbf{x}, \textbf{I})\ge 0\) (H2 holds). Moreover,

$$\begin{aligned} \lim \limits _{t\rightarrow \infty }F(\textbf{x}(t), \textbf{I}(t))|_{\textbf{x}^*}=\lim \limits _{t\rightarrow \infty }F(\textbf{x}(t), 0)=\left( \frac{\Lambda }{\mu }, 0,0\right) =\textbf{x}^*, \text {{\textbf { H1}} holds}. \end{aligned}$$

\(\square\)

Proposition 3.1

The crime and substance abuse free equilibrium , \({\textbf{E}}^0\) of system (2) is globally asymptotically stable (GAS) if \(\max \{{\mathscr {R}}_c,\mathscr {R}_a\}\le 1\).

Entrenched equilibrium

Now, we introduce a new equilibrium state for the model, called the co-occurrence entrenched equilibrium and denoted by \(\textbf{E}^{*}\). This equilibrium represents a scenario where both crime and substance abuse are prevalent and persist within the population. It’s mathematically defined by the vector:

$$\begin{aligned} \textbf{E}^{*}= \Big [{S^{*}},{P_c^{*}}, {P_a^{*}}, {P_{e}^{*}},{I_c^{*}}, {I_a^{*}}, {I_{e}^{*}},{T_c^{*}}, {T_a^{*}}, {T_{e}^{*}},{R_a^{*}}, {R_a^{*}}, {R_{e}^{*}}\Big ]. \end{aligned}$$

Due to the complexity of analytical computations, we will simulate numerically the existence and stability of \(\textbf{E}^{*}\) in Sect. 4.

Invasion reproduction number

We use a next-generation approach to derive the crime invasion reproduction number, \(_c{\mathscr {R}}_a\), when substance abuse behavior is present. The classes involved are \(P_c\), \(P_e\), \(I_c\), \(I_e\), \(T_c\), and \(T_e\). The procedures outlined in Sect. 3.1 are repeated to obtain \(_c{\mathscr {R}}_a\).

$$\begin{aligned} _c{\mathscr {R}}_a&=S^{*}{\mathscr {R}}_c +R_a^{*}\eta _{a}{\mathscr {R}}_c+ I_a^{*}{\mathfrak {R}}_c. \end{aligned}$$
(7)

Analogously, the reproductive number for substance abuse invasion under the assumption that the criminal is a resident is defined as

$$\begin{aligned} \begin{aligned} _a{\mathscr {R}}_c&= S^{*}{\mathscr {R}}_a + R_c^{*}\eta _{c}{\mathscr {R}}_a + I_c^{*}{\mathfrak {R}}_a. \end{aligned} \end{aligned}$$
(8)

Equation (7) shows that the \(S^{*}\mathscr {R}_c\) term represents the secondary cases of susceptible individuals that can arise from one individual engaging in criminal activities in the population. On the other hand, the \({R_a^*}\eta _{a}{\mathscr {R}}_c\) term indicates the secondary cases that can arise from one person involved in criminal activities among individuals who have recovered from substance abuse. The other term, \(I_a^{*}{\mathfrak {R}}_c\), gives the co-occurrence reproduction number generated by criminals in the sub-population of substance abuse engaged (\({I_a^{*}}\)). The interpretation of Eq. (8) is treated similarly to (7). The quantities \({\mathfrak {R}}_c\) and \({\mathfrak {R}}_a\) are given by

$$\begin{aligned} {\mathfrak {R}}_c =\dfrac{\omega \beta _{c}\sigma _{e}\theta _{e}+\textbf{d}_{e}\beta _{c}\sigma _{e}}{\textbf{g}_{e}({\textbf{r}}_{e}\textbf{d}_{e}+{\textbf{q}}_{e}\theta _{e})}, \, \text {and} \quad {\mathfrak {R}}_a =\dfrac{\omega \beta _{a}\sigma _{e}\theta _{e}+\textbf{d}_{e}\beta _{a}\sigma _{e}}{\textbf{g}_{e}({\textbf{r}}_{e}\textbf{d}_{e}+{\textbf{q}}_{e}\theta _{e})}. \end{aligned}$$

If \({\mathfrak {R}}_a<1\) and \({\mathfrak {R}}_c<1\) then in the long-run \(I_c^{*}{\mathfrak {R}}_a\rightarrow 0\) and \(I_a^{*}{\mathfrak {R}}_c\rightarrow 0\); this implies that co-engaged sub-population will de-escalate. With these results, we establish the following theorem.

Theorem 3.3

Assume that\({\mathscr {R}_c}>1\)and\({\mathscr {R}_a}>1\). The system (2) has a co-existence entrenched equilibrium provided that\(_c{\mathscr {R}}_a>1\)or\(_a{\mathscr {R}}_c>1\). Moreover, the equilibrium will have co-engaged individuals provided that\({\mathfrak {R}}_a >1\)and or\({\mathfrak {R}}_c>1\).

Sensitivity analysis

Understanding which factors within a population model most influence crime and substance abuse spread is crucial for designing effective control strategies and prioritizing future research46,47. Sensitivity analysis on the model’s Control Reproduction Number (CRN), \(\mathscr {R}_j\), can be conducted using Partial Rank Correlation Coefficients (PRCC) and Latin Hypercube Sampling (LHS)48. PRCC assesses CRN’s sensitivity to parameter changes, even with non-linear relationships. LHS efficiently explores the entire parameter space, providing a comprehensive understanding of how parameter variations affect CRN. By identifying the most influential factors, researchers gain valuable insights for targeted interventions and future research improvements.

The PRCC analysis in Fig. 2 unveils how different aspects of the model influence the ease with which crime or substance abuse spreads, captured by the Control Reproduction Number, \(\mathscr {R}_j\). Parameters like diffusion rate (\(\beta _j\)) and progression rate (\(\sigma _j\)) act like accelerators, positively correlated with CRN. Higher values mean people are more likely to be exposed and become infected, raising the spread potential. Similarly, a leaky rehabilitation system (\(\omega\)) and high relapse rates (\(\kappa _j\)) fuel the spread by making it harder for individuals to recover permanently. Conversely, some parameters act like brakes. Higher rehabilitation/incarceration rates (\(\theta _j\)) effectively remove individuals from the active pool, reducing the CRN and hindering spread. Faster recovery rates, either from treatment (\(\phi _j\)) or natural processes (\(\gamma _j\)), also slam on the brakes by speeding up the removal of infected individuals, making it harder for the spread to take hold. By pinpointing these key factors, the PRCC analysis provides a roadmap for designing effective control strategies. By focusing on interventions that target parameters like diffusion rate, rehabilitation effectiveness, and relapse rates, policymakers can maximize their impact in curbing the spread of crime and substance abuse.

Fig. 2
figure 2

The values of (PRCC) on the outcome of \(\mathscr {R}_j\).

Full size image

Figure 3 utilizes contour plots to unveil critical insights into controlling the spread of crime and substance abuse behaviors. These plots delve into how two key parameters influence the control reproduction number(\(\mathscr {R}_j\)), which reflects the transmission potential within the population. The first plot explores the interplay between diffusion rate (\(\beta _j\)) and rehabilitation rate \((\theta _j)\). Here, we observe that stronger rehabilitation programs (higher \(\theta _j\)) generally lead to a decrease in \(\mathscr {R}_j\). This suggests their effectiveness in curbing the spread of these behaviors. However, the plot also reveals a counteracting force: a higher diffusion rate (\(\beta\)), indicating greater interaction between individuals, tends to push \(\mathscr {R}_j\) upwards. While not explicitly shown here, the concept extends to other factors. Additional plots examining \(\beta _j\) against factors like recidivism rate (\(\kappa _j\)) or treatment-based recovery rate (\(\phi _j\)) would likely reveal further dynamics. A higher \(\kappa _j\) (individuals returning to criminal/substance abuse) might counteract the initial rise in \(\mathscr {R}_j\) caused by \(\beta _j\). Similarly, a higher \(\phi _j\) (treatment effectiveness) would likely lead to a significant decrease in \(\mathscr {R}_j\).

Fig. 3
figure 3

Contour plots of the conrol reproduction number (CRN), as a function of diffusion rate (\(\beta\)) and rehabilitation rate (\(\theta\)), recidivism rate (\(\kappa\)), and treatment based recovery rate (\(\phi\)).

Full size image

Numerical results

We simulated the model (2) using the baseline parameter values tabulated in Table 1 (unless otherwise stated), to assess the impact of the various control and mitigation strategies against the diffusion of crime and substance abuse and visualized the theoretical analysis. Simulations of the model  2 were conducted for scenarios with reproduction numbers both below and above \(1\). The results support the theoretical predictions about the stability of the crime and substance -abuse-free equilibrium, as well as the entrenched equilibrium.

Figure 4 act as a bridge between the theoretical predictions from the Theorems 3.1, 3.2 and the actual behavior of the model under different conditions. The simulations showcase how the model’s long-term outcome depends on the values of the reproduction numbers less than one. This translates to a situation where interventions or natural processes are effectively reducing the spread of both crime and substance abuse. The simulation reflects this by showing the system converging towards a state with no crime or substance abuse (equilibrium \(\mathbf {E^0}\)) as time goes to infinity. This aligns perfectly with Theorem 3.1, which both predict this stable outcome when the reproduction numbers are all below \(1\). In contrast, Fig. 5 depicts a scenario where all reproduction numbers are greater than one. Here, the model predicts a more concerning outcome. The simulation shows the system settling into a persistent state with both crime and substance abuse present (equilibrium \(\mathbf {E^*}\)). This behavior matches the prediction of Theorem 3.3. When the reproduction numbers are all above unity, the model suggests an entrenched state where both crime and substance abuse become established and persist within the population.

While the theoretical analysis of the co-existence equilibrium in our model (Eq. 2) can be mathematically challenging, it provides valuable insights. However, obtaining exact analytical solutions can be tedious. To address this, Fig. 6 presents solutions derived numerically from the model equations. These numerical solutions effectively complement the theoretical analysis and provide a clearer understanding of the co-existence equilibrium behavior.

Fig. 4
figure 4

Dynamic stability of \(\textbf{E}^0\), where \(\mathscr {R}_c= 0.1226<1\), \(\mathscr {R}_a =0.5142<1\), \({}_c{\mathscr{R}}_a= 0.12261<1\), \({}_a{\mathscr {R}}_c=0.51411<1\).

Full size image
Fig. 5
figure 5

Dynamics stability of effect of \(\textbf{E}^*\) whenever the reproduction numbers are \(\mathscr {R}_c= 7.1071>1\), \(\mathscr {R}_a =7.381>1\), \({}_c{\mathscr {R}}_a= 2.3104>1\), \({}_a{\mathscr {R}}_c=2.1974>1\).

Full size image
Fig. 6
figure 6

Co-existence of crime and substance abuse with co-engagement.

Full size image

Figure 6 presents the results of numerical simulations on the co-existence equilibrium of our model (Eq. 2). These simulations reveal a key finding: a stable positive co-existence equilibrium. This means the model predicts a scenario where both crime and substance abuse can persist within the population at a long-term stable level. The simulations also highlight a crucial factor influencing this behavior: the reproduction number. A large reproduction number suggests a high potential for both crime and substance abuse to spread. This aligns with the theoretical analysis, which emphasizes the importance of reproduction numbers in determining the model’s long-term outcome. Furthermore, the figure suggests that co-occurrence, the presence of both crime and substance abuse within individuals, is sustained as long as the invasion reproduction numbers are greater than one. Invasion reproduction numbers, as discussed earlier, quantify the potential for spread when one behavior (crime or substance abuse) is already established. So, this result implies that as long as the conditions favor the invasion of either crime or substance abuse into a population where the other behavior is already present, the co-occurrence of both will be maintained.

Following the insights gained from Fig. 6 on the co-existence equilibrium and the critical role of reproduction numbers, let’s now shift our focus to how interventions can impact the prevalence of individuals actively engaged in crime or substance abuse (denoted by \(I_c\) and \(I_a\)).

Figure 7 provides a deeper look within the crime-substance abuse co-occurrence model by examining how rehabilitation (\(\theta _a\)) and incarceration (\(\theta _c\)) impact specific subclasss. Simulations were conducted using baseline parameter values from Table 1 with varying combined rates for these interventions. The figure reveals clear trends: increasing incarceration rates (\(\theta _c\)) lead to a decrease in the criminal population (\(I_c\)), suggesting that interventions effectively remove individuals from active criminal behavior. Similarly, the prevalence of substance abusers (\(I_a\)) decreases with a higher rehabilitation rate (\(\theta _a\)), highlighting the effectiveness of strong rehabilitation programs in reducing active substance abuse. Additionally, the figure shows that co-occurrence (\(I_e\)), the presence of both behaviors within individuals, also declines with a higher rehabilitation rate (\(\theta _e\)). This suggests that by effectively addressing substance abuse, rehabilitation programs can indirectly reduce the co-occurrence of crime and substance abuse within the population. By dissecting these subclass effects, Fig. 7 underscores the targeted impact of rehabilitation and incarceration strategies. These interventions not only hold promise for reducing the overall co-occurrence of crime and substance abuse but also for decreasing the populations of criminals and substance abusers individually.

Fig. 7
figure 7

This figure shows the simulated impact of rehabilitation and incarceration rates on the prevalence of engaged individuals in the crime-substance abuse co-occurrence model. The simulations were conducted using the baseline parameter values from Table 1 and varying the combined rehabilitation and incarceration rate (\(\theta _c = \theta _a=\theta _e\)) across different values: \(0, 0.25, 0.5, 0.75\) and \(1\).

Full size image

Now, we investigate the influence of crime recidivism rate (represented by \(\kappa _c\)), substance abuse relapse rate (represented by \(\kappa _a\)), and the combined relapse rate for co-occurrence (represented by \(\kappa _e\)) on the prevalence of criminals (\(I_c\)), substance abusers (\(I_a\)), and individuals exhibiting co-occurrence (\(I_e\)), respectively.

Figure 8 confronts us with a critical roadblock in curbing crime and substance abuse: recidivism rates. Simulations within the figure paint a troubling picture. As relapse rates for crime (\(\kappa _c\)), substance abuse recidivism rate (\(\kappa _a\)), and co-occurrence re-offend (\(\kappa _e\), where co-occurrence refers to both criminal and substance abuse behaviors) rise, so too does the prevalence of active criminals (\(I_c\)), active substance abusers (\(I_a\)), and individuals exhibiting co-occurrence (\(I_e\)) within the population. This trend underscores the urgency of tackling recidivism for any significant intervention strategy. High recidivism rates suggest our current rehabilitation programs and support systems may be falling short. To break this cycle, a multi-pronged approach is needed. This could involve strengthening rehabilitation programs with evidence-based practices, implementing robust transitional support after release, and addressing broader social and economic factors that contribute to crime and substance abuse. By prioritizing these measures, we can potentially decrease recidivism rates and ultimately reduce the prevalence of crime and substance abuse within the population. Figure 8 thus serves as a call to action, urging us to prioritize strategies that not only address criminal behavior and substance abuse but also focus on preventing individuals from returning to these activities.

Fig. 8
figure 8

Numerical simulation of the left side impact of relapse and recidivism rate on actively engaged sub-classes, where \(\kappa _c=\kappa _a=\kappa _e=0.0, 0.2, 0.4, 0.6, 0.8, 1\).

Full size image

Furthermore, we investigate the influence of cross-recidivism on the populations engaging in criminal and substance abuse activities. Cross-recidivism refers to the phenomenon where individuals who recover from one type of behavior (crime or substance abuse) relapse into the other type of behavior.

Figure 9 sheds light on the concerning impact of cross-recidivism through numerical simulations. The figure shows that as the rate of individuals relapsing from substance abuse into crime increases (\(\tau _a\)), the criminal population also increases. Similarly, a rise in the rate of individuals relapsing from crime into substance abuse (\(\tau _c\)) leads to a higher prevalence of substance abuse. This emphasizes the importance of addressing cross-recidivism in rehabilitation programs. These findings suggest that a holistic approach is crucial for policymakers and rehabilitation experts. Treatment programs should not only address criminal behavior but also consider the potential for relapse into substance abuse, and vice versa. By incorporating strategies to mitigate cross-recidivism, we can potentially create more effective interventions that keep individuals from returning to destructive behaviors.

Fig. 9
figure 9

Numerical simulation of the impact of cross-recidivism. This simulation conducted with the baseline value in Table  1 and \(\tau _a=\tau _c=0, 0.25, 0.5, 0.75, 1\).

Full size image

Conclusion

This research successfully developed and analyzed a well-posed mathematical model to study the spread of crime and substance abuse within a population. The model incorporates social interaction, rehabilitation efforts, recovery rates, and relapse rates, providing a comprehensive picture of these behaviors’ complex dynamics. Key findings bridge the gap between theoretical predictions and practical applications.

The Control Reproduction Number (\(\mathscr {R}_j\)), computed using the next-generation matrix approach for single and co-occurrence systems, is crucial for understanding the spread. A CRN exceeding unity signifies a persistent behavior, while a CRN below unity suggests interventions or natural processes that can lead to eradication. Simulations corroborated theoretical predictions about the stability of crime/substance-free and co-existence states. The model predicts a stable crime/substance-free state when interventions are successful ( \(\mathscr {R}_j< 1\)) and a persistent co-existence state when the spread is not adequately controlled ( \(\mathscr {R}_j> 1\)). Numerical simulations incorporating the invasion reproduction number revealed a stable co-existence equilibrium where both crime and substance abuse can persist long-term. This number quantifies the spread potential when one behavior is already established. The invasion reproduction number and the CRN play a critical role in determining the existence and stability of this co-existence equilibrium, as identified through theoretical analysis.

Sensitivity analysis using techniques like PRCC and LHS identified key parameters influencing the CRN, such as diffusion rate, progression rate, rehabilitation/incarceration rate, recovery rates, and relapse rates. Numerical simulations allowed us to assess the impact of these parameters. For example, increasing incarceration rates reduced the criminal population, while higher rehabilitation rates lowered substance abusers. The analysis of interventions like rehabilitation and incarceration demonstrated their effectiveness but also highlighted challenges like recidivism rates. Simulations emphasized the need for stronger rehabilitation programs, improved support systems, and interventions that address the underlying social and economic factors contributing to these behaviors.

This research offers valuable insights for policymakers. By focusing on interventions that target key parameters identified through simulations and aligned with theoretical analysis (e.g., reducing diffusion rates, improving rehabilitation effectiveness, and preventing relapse), policymakers can maximize their impact in curbing the spread of crime and substance abuse. Furthermore, the model’s well-posed framework can be adapted to study the spread of other infectious or behavioral phenomena within a population.

While this research provided valuable insights on the spread of crime and substance abuse, future work can be directed towards optimal control analysis to identify the most significant interventions. This would involve strategically allocating resources for rehabilitation, social support, and law enforcement based on factors like risk profiles and long-term impacts. Additionally, the model can be refined to consider spatial variations, media influence, and crime/substance types. Data quality and limitations like social-economic factors and population heterogeneity should also be addressed in future research. By incorporating optimal control and selective interventions, this research can significantly contribute to developing targeted strategies for curbing crime and substance abuse.