Mathematical modeling and analysis of the co-dynamics of crime and drug abuse

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 method⁴⁴ 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}and\({\mathscr {R}_a}\}, 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 and \(\mathscr {R}_a, 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-Chavez⁴⁵, 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-Chavez⁴⁵. 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}

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 and \({\mathfrak {R}}_c 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 research^46,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)⁴⁸. 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

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

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$.