Modelling the role of correctional services on gangs: insights through a mathematical model

Research has shown that gang membership increases the chances of offending, antisocial behaviour and drug use. Gang membership should be acknowledged as part of crime prevention and policy designs, and when developing interventions and preventative programmes. Correctional services are designed to rehabilitate convicted offenders. We formulate a deterministic mathematical model using nonlinear ordinary differential equations to investigate the role of correctional services on the dynamics of gangs. The recruitment into gang membership is assumed to happen through an imitation process. An epidemic threshold value, Rg, termed the gang reproduction number, is proposed and defined herein in the gangs’ context. The model is shown to exhibit the phenomenon of backward bifurcation. This means that gangs may persist in the population even if Rg is less than one. Sensitivity analysis of Rg was performed to determine the relative importance of different parameters in gang initiation. The critical efficacy ε* is evaluated and the implications of having functional correctional services are discussed.


Introduction
Correctional services in South Africa provide needs-based correctional sentence plans and interventions that are based on an assessment of the security risk and criminal profile of individuals. The corrections target all elements associated with offending behaviour and focus on the offence for which a person was sentenced to correctional supervision, remanded in a correctional centre or released on parole [1]. Correctional programmes and/or interventions can be viewed as a structured set of 2017 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.

Model formulation
We consider a population whose size is N(t), at any time t. The population is divided into four disjoint independent classes or compartments based on an individual's status and risk factors with respect to gang membership. The class S n (t) represents individuals not at risk of becoming gang members, S r (t) represents individuals at risk of becoming gang members, G(t) represents gang members and lastly, C(t) represents those in correctional services. The total population at any time t is thus given by The general population enter the susceptible population at a rate Λ. Among individuals entering the susceptible population, we have that a proportion p of these individuals are recruited into the class of susceptible individuals not at risk of joining a gang and the complementary proportion (1 − p) join susceptible individuals at risk of joining a gang. Therefore, we neglect the possible recruitment of individuals already belonging to gangs. Transition rate from no risk susceptibility into at risk susceptibility is represented by θ . Unlike in [10], the change in the risk status is not driven by interacting with gang members but simply by change of one's environmental conditions. A typical example is that of a slowing down economy that results into retrenchment. Once the economic status of an individual changes, then susceptibility to committing crimes may increase. This is particularly important in the South African context as different living environments often determine the risk.
The recruitment of individuals into gangs is assumed to follow an imitation process, described comprehensively in [19]. We propose an initiation function f (S r , G) = βG(1 + ηG)S r into gang membership that is driven by imitation, with β as the effective contact rate and η as the imitation  coefficient. Once initiated, we also assume that gang members can either revert back to compartment S r at a rate σ 2 or are sent to correctional facilities through convictions and sentencing at a rate γ . Depending on the efficacy ε (where ε ∈ [0, 1]) of correctional services, a released inmate may either join a gang again at a rate (1 − ε)σ 1 or may rejoin the community as either a susceptible at risk at a rate (1 − ν)εσ 1 or those not at risk at a rate νεσ 1 . By efficacy, we mean the measure to which a policy, programme or initiative meets its intended result with ε = 1 signifying that no individuals will revert to gangs when they leave correctional services. This represents completely effective correctional programmes. ε = 0 signifies that all individuals in correctional facilities will revert back to gangs upon their release, while 0 < ε < 1 implies that correctional programmes will be effective to some degree. In reality, ε ∈ (0, 1). A summary of the description of parameters together with their estimated values is given in table 1. Figure 1 shows the movement of humans as their status with respect to gang membership changes. Combining the parameters, assumptions and the schematic diagram, we have the following general set of nonlinear ordinary differential equations: where we assume that all the model parameters are positive. The positivity of the solutions of system (2.1) can easily be established if S n0 ≥ 0, S r0 ≥ 0, G 0 ≥ 0, C 0 ≥ 0, see for instance [20][21][22].

Model analysis
3.1. Invariant region It follows from system (2.1) that dN/dt = Λ − μN. Then, sup t→∞ N(t) ≤ Λ/μ. We can thus study (2.1) in following feasible region which is positively invariant with respect to system (2.1). This means that our system is well posed and all solutions of system (2.1) with (S n0 , S r0 , G 0 , C 0 ) ∈ R 4 + remain in Γ for all t > 0.

Gang-free equilibrium and the gang reproduction number
The model has a gang-free equilibrium given by a scenario depicting a gang-free state in the community or society. The gang reproduction number R g of the model, is defined herein in the gang membership context as the average number of people that each single gang member will initiate to a gang during his/her membership in a wholly susceptible population. This threshold quantity is analogous to the basic reproduction number in mathematical epidemiology described in [23,24]. Usually, R g < 1 implies that gangs will decline, whereas R g > 1 implies that gangs will persist within a community and R g = 1 requires further investigation. The determination of R g is done using the next generation matrix approach [24]. This method has been explored in many papers [25][26][27][28][29]. Driessche & Watmough [24] describe the following method to determine the reproduction number: with each x i ≥ 0, be the number of individuals in each compartment. Denote m to be the number of compartments corresponding to infected individuals where the epidemiological interpretation of the model determines between infected and uninfected compartments. More than one interpretation is possible for some models. Define X s to be the set of all disease free states given by be the rate of transfer of individuals into compartment i by all other means, and V − i (x) be the rate of transfer of individuals out of compartment i. It is assumed that each function is continuously differentiable at least twice in each variable. Consider the disease transmission model with non-negative initial conditions given by If x 0 is a disease free equilibrium of (3.1) and f i (x) satisfies assumptions (A1)-(A5) given in Driessche & Watmough [24], then the reproduction number of (3.1) is the spectral radius of the next generation matrix FV −1 where where F is non-negative and V is a non-singular M-matrix. Using this method we have The gang members' compartments are G and C, giving m = 2. Then

Sensitivity analysis
We examine which model parameter has the greatest effect on the value of the gang reproduction number R g . Determining these parameters is useful in reducing the recruitment of new gang members given that R g is directly related to gang initiation. Following Chitnis et al. [30], we calculate the sensitivity indices of the gang reproduction number R g , to the parameters in the model. These indices indicate how sensitive R g is to a change in each parameter, in other words, this tells us how crucial each parameter is to gang initiation. Since there are usually errors in data collection and presumed parameter values, sensitivity analysis is commonly used to determine the robustness of model predictions to parameter values [30]. Sensitivity indices allow us to measure the relative change in a state variable when a parameter changes. The normalized forward sensitivity index (NFSI) of the gang reproduction number R g to a parameter is the relative change in the variable R g to the relative change in a given parameter. A directly proportional normalized sensitivity index indicates that an increase/decrease in the parameter value brings about an increase/decrease, respectively, in the value of R g , whereas, an inversely proportional normalized sensitivity index indicates that an increase in the parameter value brings about a decrease in the value of R g . When R g is a differentiable function with respect to each of its parameters, then the sensitivity index may be alternatively defined using partial derivatives as follows.
Then, for every parameter q ∈ V, the NFSI of R g is defined as: Using an explicit formula for R g (3.2) and definition 3.1, the sensitivity indices of R g with respect to each of its parameters are calculated. Recall that μ is the natural death rate. Thus, the sensitivity index of R g with respect to μ has been omitted because it is clear that increase in this rate is neither ethical nor practical. , .
From the calculations here we see that R g is most sensitive to changes in the values of β and Λ. An increase in either of these results in an increase of the same proportion in R g and a decrease in either of these will bring about an equivalent decrease in the value of R g ; they are directly proportional. Also, R g has a direct proportional relationship with parameters σ 1 and θ, however with a proportionally smaller increase or decrease. Parameters σ 2 , p, γ and ε have an inversely proportional relationship with R g ; an increase in any of them will bring about a decrease in R g . This is a reflection that increasing the efficiency To further understand the model reproduction number in the context of gangs, we can deduce the threshold efficacy by setting R g = 1. It can easily be established that , (3.4) is the basic reproduction number, the value of R g in the absence of correctional services i.e ε = 0. In the absence of correctional services (obtained by setting ε = 0), the model assumes that there is no rehabilitative correction and individuals released from correctional facilities go back into gangs. An efficacy of ε = 0 depicts totally dysfunctional correctional services, while ε = 1 signifies that correctional services will be 100% effective. A high value of the efficacy of correctional services in any given population impacts the reproduction number over time. The question then is: what is the threshold efficacy necessary for the reduction of the reproduction number to below one? Absence of correctional services here means that jails do not act as rehabilitation and correctional facilities. So gangs can be contained or eradicated if the efficacy of correctional services is maintained above ε * . This clearly shows the need to have restorative and corrective prisons for gang members. Some corrective interventions include skilling, counselling and education of inmates. Research has shown that offenders who undergo programmes such as the provision of education, employment and other correctional programmes (e.g. substance abuse, violence prevention, sexual offending prevention, family violence prevention), at the most appropriate time in the offender's sentence, contributes to safe transition to the community. Education programmes in custodial settings are known to decrease recidivism [3].

Local stability of the gang-free steady state
We shall now prove the local stability of the gang-free equilibrium point G 0 whenever the gang reproduction number R g is less than unity.
Theorem 3.2. The gang-free equilibrium point G 0 is locally asymptotically stable if R g < 1 and unstable otherwise.
Proof. The Jacobian matrix evaluated at G 0 is The eigenvalues are given by λ 1 = −(μ + θ ), λ 2 = −μ and the solution of This gives We note that when R g < 1, then the remaining eigenvalues will be both negative. This completes the proof.

Gang-persistent equilibrium
In this section, we determine the gang-persistent equilibrium point denoted by G * = (S * n , S * r , G * , C * ). The gang-persistent equilibrium always satisfies From the last equation of (3.5), we have that Substituting equation (3.6) into the first and third equation of (3.5) leads to Substituting equations (3.6) and (3.7) into the second equation of (3.5) leads to the following quadratic equation in terms of G * aG * 2 + bG * + c = 0, and .
For the gang-persistent equilibrium to exist, the solutions of (3.8) must be real and positive. We note Since a = 0, equation (3.8) is a quadratic equation with respect to G * . Let the discriminant of (3.8) be denoted by , so that = R * g + R g − 1. (3.10) Solving (3.10) for = 0 in terms of R g , we get We clearly note the following relations: Equation (3.8) has real roots provided ≥ 0. We thus have the following results on existence of the gang-persistent equilibrium. (H1) Let η = 0. Then, system (2.1) has a unique gang-persistent equilibrium when R g > 1.
Epidemiologically, a backward bifurcation entails that it is not enough to only reduce the basic reproductive number to less than 1 to eliminate a disease. On most part there are two distinct bifurcations at R 0 = 1 namely; forward (supercritical) and backward (subcritical). A backward bifurcation happens when R 0 is less than unity, a small positive unstable equilibrium appears while the disease-free equilibrium and a larger positive equilibrium are locally asymptotically stable. On the other hand, a forward bifurcation happens when R 0 crosses unity from below, a small positive asymptotically stable equilibrium appears and the disease-free equilibrium looses its stability [31]. The phenomenon of backward bifurcation was first found in epidemiological models by Huang et al. [32]. Studies supporting backward bifurcations include those in [33][34][35][36][37].
From theorem 3.3, we observe that if η = 0, then system (2.1) has a unique gang-persistent equilibrium when R g > 1. For this case, the bifurcation at R g = 1 is forward. However, if η > 0, system (2.1) has two different gang-persistent equilibria when η > η * and 1 − R * g < R g < 1. Hence, system (2.1) has a backward bifurcation at R g = 1 from the gang-free equilibrium to two gang-persistent equilibria. To conclude, we now show existence of backward bifurcation.

Backward bifurcation
Conditions for the existence of backward bifurcation follow from Theorem 4.1 proved in [31]. We deliberately avoid rewriting the theorem and focus on its application. Let us make the following change of variables: S n = x 1 , S r = x 2 , G = x 3 , C = x 4 , so that N = 4 n=1 x n . We now use the vector notation X = (x 1 , x 2 , x 3 , x 4 ) T . Then, model system (2.1) can be written in the form dX/dt = F(t, Let β be the bifurcation parameter, R g = 1 corresponds to The Jacobian matrix of model system (2.1) at G 0 when β = β * is given by Model system (3.12), with β = β * has a simple eigenvalue, hence the centre manifold theory can be used to analyse the dynamics of model system (2.1) near β = β * . It can be shown that J * (G 0 ), has a right eigenvector given by w = (w 1 , w 2 , w 3 , w 4 ) T , where and Further, the left eigenvector of J * (G 0 ), associated with the zero eigenvalue at β = β * is given by The computations of a and b are necessary in order to apply Theorem 4.1 in [31]. For system (3.12), the associated non-zero partial derivatives of F at the gang-free equilibrium are as follows: , , It thus follows that where η * is given in equation (3.9). Note that if η > η * then a > 0 and a < 0 if η < η * . Lastly, We thus have the following result.
From the results obtained above, we note that a backward bifurcation occurs at R g = 1 if and only if η > η * is satisfied. From this, we can deduce that when the imitation coefficient, η exceeds the critical threshold η * , then the gang population remains high leading to a backward bifurcation (figure 2). We show the existence of a backward bifurcation through numerical example by creating a bifurcation diagram around R g = 1 (figure 2). To draw a bifurcation curve (the graph of G * as a function of R g ), we fix Λ = 0.047; μ = 0.02; β = 0.3; η = 10.0; p = 0.4; ν = 0.7; θ = 0.13; γ = 0.5; σ 1 = 0.1; σ 2 = 0.5; ε = 0.5. For this case, we have that η * = 2.7595 and R g = 0.8823. Generally speaking, in many epidemic models the basic reproduction number, R 0 , which is the key concept in epidemiology can be decreased below unity to eradicate the disease. However, in our model, this classical R 0 -threshold is not the key to control the spread of gangs within a population. In fact, the existence of backward bifurcation entails that gangs may persist in the population even with values of R g less than unity. Our findings suggest that keeping the imitation coefficient η below a certain threshold η * is an effective way to avoid backward bifurcation.

Numerical simulations and results
The estimation of parameters in any model validation process is a challenging task. We make some hypothetical assumptions for the purpose of illustrating the usefulness of our model in tracking the dynamics of gangs passing through correctional services. Demographic parameters are the easiest to estimate in this instance. For the per capita death rate μ, we assume that the life expectancy of the human population is 60 years. This value has been the approximation of the life expectancy in South Africa [38]. This translates into μ = 0.0166 per year. The recruitment of individuals in the community is linked to the birth rate. The birth rate of South Africa is on average 0.028 [39]. We thus choose a value for Λ = 0.028. The parameters ε, ν and p all lie in the interval (0, 1). The remaining parameters are estimated since most of them are not available in the literature and are given in table 1.
We begin by illustrating the analytic results in which the gang-free equilibrium G 0 is locally asymptotically stable when R g < 1. The results are illustrated in figure 3.
We investigate the impact of the efficacy parameter ε on the population levels of gang members. Figure 4 shows the effects of increasing ε on the number of gang members. We hypothetically start at ε = 0.6 and observe that increasing ε lowers the number of gang members. One can quantify the percentage decrease in the number of gang members when ε is increased by 0.1. For instance, an increase of ε for 0.6 to 0.7 reduces the number of gang members by approximately 14%. Of particular importance in the fight against gangsterism is the number of convictions on committed crimes that results in the placement of gang members in correctional services. This has been an issue of considerable concern in South Africa given the existing large variation between the number of committed crimes and the number of convictions, with conviction rates being as low as 10%. In figure 5, we begin by hypothetically setting γ = 0.25 and observe that increasing conviction rates with a functional correctional system can lead to a reduction in the number of gang members. We observe that an increase of γ from 0.25 to 0.3 results in an approximate decrease of gang members by 17%.
Once an individual belongs to a gang, one has a choice of remaining a gang member and risk arrest as a result of gang related crimes or quitting all together. To investigate the choice a gang member has to undertake, we make a contour plot (figure 6) to show how the parameters σ 2 (the rate  of voluntary quitting of gang member) and γ (the rate of convictions and placement in correctional services) affect R g . The results show that increasing σ 2 coupled with decreasing γ leads to a decrease in R g . This is of particular importance as it alludes to interventions that are not correctional. In fact, if one quits being a gang member before a conviction, that is, before crimes are committed, then its beneficial to the individual and the community since resources allocated to correctional services are saved.

Conclusion
We developed a simple mathematical model to investigate the role of correctional services on gangs. Principles drawn from the literature of mathematical epidemiology have been used to model how individuals are recruited into gangs and their possible recovery. Initiation into gangs has been assumed to happen through an imitation process in which peer influence is central to joining gangs. The growth and decrease of gang members was driven by the gang reproduction number, R g , as in the case of epidemic models. However, in our model, this classical R g -threshold is not the key to control gangs within communities. In fact gangs may persist in the population even with subthreshold values of R g . It was shown to happen, in particular when the value of the imitation coefficient is high enough such that the relation η > η * is satisfied. In the absence of the imitation coefficient, that is, when η = 0, the model in this study will have a unique gang-persistent equilibrium. However, the introduction of the imitation coefficient leads to multiple equilibria and seem to be responsible for interesting dynamical aspects such as the occurrence of a backward bifurcation. This means gangs may persist in the population even with subthreshold values of R g . Thus, awareness programmes and/or specific health programmes may be employed to reduce η or, at least, to increase the value of η * . Our results put into evidence the importance to identify those social processes, as the imitation mechanism, that may facilitate or counteract the spread of gangs within a community of individuals. Some precise knowledge of these mechanisms is in fact essential to develop effective policies that will impede the spread of gangs within a community.
Sensitivity analysis have been performed by evaluating the sensitivity indices of the gang reproduction number, R g , to model parameters. Since R g is a measure of initial gang membership, these sensitivity indices allow us to determine the relative importance of different parameters in gang initiation. It was observed that R g has a direct proportional relationship with the parameters β, Λ, σ 1 and θ , whereas parameters σ 2 , p, γ and ε have an inversely proportional relationship with R g . To further understand the role of correctional services on gangs, the threshold efficacy ε * was established. It was observed that gangs can be contained or eradicated if the efficacy of correctional services is maintained above ε * . Thus, it is important to have efficient correctional services in the fight against gangs. We also investigated the role of the efficacy parameter together with other important parameters such as conviction rates and self recovery through graphical plots.
Standard statistical techniques for collecting data on gangs such as household surveys are expensive and should, at best, be carried out every three to five years. Also, reliable gang-related data is elusive. Therefore, mathematical models become useful tools as they allow the extent of the phenomenon to be estimated. While the model presented in this study is theoretical in nature, it presents very useful and practical results that can be of help to policy makers in fighting against gangsterism, gang violence and its related crimes that have ravaged communities. Like in any model development, our model is not without limitations. The model presented in this paper assumes homogeneous mixing which is practically impossible in communities with gangs. In practice, susceptibility varies. This is because of differences in behavioural, social and environmental factors. An individual-based model could be used