Incorporation of awareness programs into a model of the spread of HIV/AIDS among people who inject drugs

Mathematical modeling techniques have been used extensively during the human immunodeficiency virus (HIV) epidemic. Drug injection causes increased HIV spread in most countries globally. The media is crucial in spreading health awareness by changing mixing behavior. The published studies show some of the ways that differential equation models can be employed to explain how media awareness programs influence the spread and containment of disease (Greenhalgh et al. Appl Math Comput. 2015;251:539–563). Here we build a differential equation model which shows how disease awareness programs can alter the HIV prevalence in a group of people who inject drugs (PWIDs). This builds on previous work by Greenhalgh and Hay (1997) and Liang et al. (2016). We have constructed a mathematical model to describe the improved model that reduces the spread of the diseases through the effect of awareness of disease on sharing needles and syringes among the PWID population. The model supposes that PWIDs clean their needles before use rather than after. We carry out a steady state analysis and examine local stability. Our discussion has been focused on two ways of studying the influence of awareness of infection levels in epidemic modeling. The key biological parameter of our model is the basic reproductive number R0$$ {R}_0 $$ . R0$$ {R}_0 $$ is a crucial number which determines the behavior of the infection. We find that if R0$$ {R}_0 $$ is less than one then the disease‐free steady state is the unique steady state and moreover whatever the initial fraction of infected individuals then the disease will die out as time becomes large. If R0$$ {R}_0 $$ exceeds one there is the disease‐free steady state and a unique steady state with disease present. We also showed that the disease‐free steady state is locally asymptotically stable if R0$$ {R}_0 $$ is less than one, neutrally stable if R0$$ {R}_0 $$ is equal to one and unstable if R0$$ {R}_0 $$ exceeds one. In the last case, when R0$$ {R}_0 $$ is greater than one the endemic steady state was locally asymptotically stable. Our analytical results are confirmed by using simulation with realistic parameter values. In nontechnical terms, the number R0$$ {R}_0 $$ is a critical value describing how the disease will spread. If R0$$ {R}_0 $$ is less than or equal to one then the disease will always die out but if R0$$ {R}_0 $$ exceeds one and disease is present the disease will sustain itself and moreover the numbers of PWIDs with disease will tend to a unique nonzero value.

some of the ways that differential equation models can be employed to explain how media awareness programs influence the spread and containment of disease (Greenhalgh et al. Appl Math Comput. 2015;251:539-563). Here we build a differential equation model which shows how disease awareness programs can alter the HIV prevalence in a group of people who inject drugs (PWIDs). This builds on previous work by Greenhalgh and Hay (1997) and Liang et al. (2016).
We have constructed a mathematical model to describe the improved model that reduces the spread of the diseases through the effect of awareness of disease on sharing needles and syringes among the PWID population. The model supposes that PWIDs clean their needles before use rather than after. We carry out a steady state analysis and examine local stability. Our discussion has been focused on two ways of studying the influence of awareness of infection levels in epidemic modeling. The key biological parameter of our model is the basic reproductive number R 0 . R 0 is a crucial number which determines the behavior of the infection. We find that if R 0 is less than one then the disease-free steady state is the unique steady state and moreover whatever the initial fraction of infected individuals then the disease will die out as time becomes large. If R 0 exceeds one there is the disease-free steady state and a unique steady state with disease present. We also showed that the disease-free steady state is locally asymptotically stable if R 0 is less than one, neutrally stable if R 0 is equal to one and unstable if R 0 exceeds one. In the last case, when R 0 is greater than one the endemic steady state was locally asymptotically stable. Our analytical results are confirmed by using simulation with realistic parameter values. In nontechnical terms, the number R 0 is a critical value describing how the disease will spread. exceeds one and disease is present the disease will sustain itself and moreover the numbers of PWIDs with disease will tend to a unique nonzero value.

K E Y W O R D S
human immunodeficiency virus, local and global stability, mathematical modeling, reproduction number

J E L C L A S S I F I C A T I O N
Applied science for engineering INTRODUCTION

Models of HIV/AIDS
Human immunodeficiency virus (HIV) was discovered in the early 1980s. It is an immunodeficiency virus that weakens the immune system and leads to greater susceptibility to infection and eventually leads to acquired immune deficiency syndrome (AIDS). HIV is transmitted by various methods including homosexual transmission and heterosexual transmission. However intravenous drug use among people who inject drugs (PWIDs) is another common way that HIV can spread. The majority of infected individuals have a short flu-like illness which lasts for a couple of weeks. 1 But after that the individual may be asymptomatic for many years. Untreated individuals move from acutely infected, to asymptomatic, to pre-AIDS symptoms, to full-blown AIDS. There are tests which can detect HIV infection. As HIV/AIDS is still an important problem it is crucial to develop mathematical models which can help us to do this. The first person to study mathematical models of how HIV spreads among PWIDs was Kaplan. 2 Kaplan models how HIV is transmitted in groups of PWIDs who share needles. The places they do this are called shooting galleries. Kaplan models the spread of HIV using two differential equations, one for the PWIDs and one for needles. Caulkins and Kaplan 3 use Kaplan's original model as a basis for assessing the effect of HIV/AIDS on the size of the sharing, injecting population. In a groundbreaking paper, Kaplan and O'Keefe 4 used a simple differential equation model to examine a test needle exchange program implemented by the New Haven Health Department in Connecticut, USA. As a result of this study which definitely showed that prevalence of HIV could be decreased by up to a third by needle exchanges, needle exchanges were legalized in various states and cities in the USA.
Greenhalgh and Hay 5 examined an extension of Kaplan's original model to remove some of its deficiencies. They extended Kaplan's original model in order to permit infected PWIDs to leave a needle virus free after use and used a joint probability distribution of the chance the syringe is left virus-free after use and the PWID is infected. They also permitted PWIDs who discover that they are HIV positive to reduce their rate of sharing. Moreover Greenhalgh and Hay performed a mathematical analysis of this model.
Lewis and Greenhalgh extended Greenhalgh and Hay's model to take account of the fact that infected PWIDs go through three stages of infectiousness before developing AIDS, an initial acute infectious stage, an asymptomatic stage and a pre-AIDS stage. Lewis and Greenhalgh produced several papers with needles having only one infectious stage and PWIDs having three infectious stages (the simple model 6 ) and needles and PWIDs having three infectious stages (the Optimistic Model, 7 the Pessimistic Model, 8

and the General Model 9 ).
There is much other work on the spread of HIV among PWIDs, for example, Kretzschmar and Wiessing 10 who model HIV in social networks of PWIDs, Strathdee et al. 11 who review what causes HIV in PWIDs, and Liang et al., 12 who study a stochastic model based on Greenhalgh and Hay's model. More recent work on the spread of HIV among PWIDs includes Babaei et al. 13 who construct a mathematical model for transmission of HIV/AIDS among PWIDs in Iranian prisons. The authors model only PWIDs not syringes but use a large number of PWID classes. If no treatment is available they look at the stability of the model. They then incorporate treatment of PWIDs, look at its effects and compare the basic reproduction numbers with and without treatment. Des Jarlais et al. 14 looked at the spread of HIV among PWIDs in New York City, in the period 2012-2019, when there was little HIV and during the COVID pandemic which interfered with preventive measures. Using an agent-based model they modeled sharing syringes and equipment to measure New York HIV incidence. They found that they key drivers for the very low HIV levels were that very few PWIDs were likely to spread the disease and most sharing was between individuals of the same serostatus. It would be challenging to contain an epidemic if the pandemic continued for a long time. It was imperative that HIV prevention measures be restored to their pre-COVID levels as rapidly as possible.
Zhang et al. 15 highlight the problem of simultaneous disease transmission in multiple subpopulations. They address this by developing a social network with disease spreading in several subpopulations at the same time. This is intended to capture the multimode transmission between female sex workers and the PWIDs population. They develop a new random search method to find the best subset of individuals in the population on which to focus interventions. Fraser et al. 16 developed a dynamic HIV transmission model among key populations (KPs) to determine the extent to which their unmet prevention and treatment needs is driving HIV transmission. They estimated the proportion of new infections acquired in each KP, and the proportion due to their unprotected risky practices. They suggest 93.8% of new infections over 2020-2029 will be due to unprotected sex between men who have sex with men (MSM) or injection by PWIDs. Among Chinese young people most cases of HIV happen in MSM. Jing et al. 17 use a metapopulation approach to examine how pre-exposure prophylaxis (PrEP) will affect the spread of a theoretical novel HIV genotype in this subpopulation in Guangdong, China. This article is relevant to forecasting the spatial HIV spread in MSM from diverse regions and looks at the role of prefecture-level cities in containing HIV in the Guangdong region.
In 2009, there was an undetected hepatitis C virus (HCV) outbreak in Athens, Greece. This was closely followed by an HIV outbreak in 2010, which was detected in 2011 and in 2012 HIV control measures were introduced. These managed to lower the HIV cases. Gountas et al. 18 show that if the HCV epidemic had been identified and dealt with early the HIV epidemic could have been prevented and 35.2-53.2 million euros could have been saved. Flountzi et al. 19 later model interventions during the 2012-2013 Athens HIV outbreak.
Salman 20 discusses a model for transmission of HIV including treatment. The model does not specify the mode of transmission and considers only individuals, not individuals and syringes. The effect of treatment is to move individuals from the symptomatic class back to the asymptomatic class. A memory effect is introduced by using fractional differential equations. A disease awareness program similar to the one used by Li et al. 21 is introduced to model the effect of media awareness. A steady state and stability analysis is undertaken. Finally, the consequences of a time delay in the treatment are discussed.

Disease awareness programs
So we have established that the spread of HIV and AIDS among PWIDs is an important problem that needs attention.
However recently there has also been increased interest in the effect of disease awareness programs. Infected individuals may adjust their behavior to reduce potential contacts in the presence of high levels of disease. There are two ways to model this. The first and simplest way, and the one which we shall adopt in this article, is to reduce the disease transmission term by a factor (0 ≤ < 1) to take account of the behavioral modifications individuals make because of their knowledge of current disease levels. In the second the amount of media awareness is modeled as a separate variable and the individuals (usually the susceptible class) are split into aware and unaware individuals. 22,23

Disease awareness programs using a multiplicative factor
We shall first look at disease awareness models that reduce the disease awareness function by a factor (I) between 0 and 1. Xiao and Ruan 24 study an SIR model where the disease transmission function in the absence of an awareness program is SI, where S is the number of uninfected individuals, I is the infected individuals, and is a constant. This is multiplied by a disease awareness function Here is a positive constant. We see clearly that (I) is a strictly positive monotone decreasing function between zero and one.
Li et al. 21 look at an SIS epidemic model with constant and impulsive vaccination and where there is media awareness and the disease transmission term, again fundamentally SI, is reduced by a factor Here a and b are positive constants and a ≤ 1. Note that again (I) is a monotone decreasing function of I. This approach and the same form of media awareness function is also used by Tuenche et al. 25 in a mathematical model of the spread of influenza. Liu 26 investigates the spread of disease in an SIRS model using a similar awareness function although stochasticity is introduced. Salman 20 uses the same disease awareness function in his model. The key thing to note for our purposes is that the disease transmission term is reduced by the same disease awareness program factor as in Li et al. 21 because of awareness of infected individuals. The same function will be used to model the reduction in disease transmission in some of our numerical examples. Cui et al. 27 consider an SEI model where the disease transmission function (again effectively SI with no disease awareness) is reduced by a multiplicative factor (I) = e −mI due to the effect of disease awareness. Again we will use a similar function in our numerical examples. Liu et al. 28 consider a model for an SEIH (susceptible-exposed-infectious-hospitalized) epidemic. The fundamental disease transmission term is again SI, which is reduced by a factor e −a 1 E−a 2 I−a 3 H , where a 1 , a 2 , and a 3 are constants and E and H are the number of exposed and hospitalized people. Strictly speaking this is not the same type of disease awareness function as we have been discussing as it depends on both of the exposed and hospitalized individuals in addition to infected, but it is based on the same idea.
In a different paper, Cui et al. 29 look at an SIS model in which the basic disease transmission function with no awareness is With behavioral modification due to knowledge of disease levels this is decreased by a multiplicative term (I) = 1 − kf (I) where k < 1. Here f (I) is a positive monotone increasing function with f (0) = 0 and lim I→∞ f (I) = 1. This disease awareness program function is a generalization of the one used by Li et al.. 21 Sun et al. 30 study the effect of media-induced social-distancing on how disease spreads in a setting with two patches using a similar modification of the disease awareness function.

Disease awareness programs modeling unaware and aware individuals
Now we turn to models which use the other approach that is they divide the population into aware and unaware individuals and model the amount of disease awareness as a separate variable in some way. These models are necessarily more complex as they have more classes, but nowadays are used more often to model disease awareness programs. There are many models that use this type of disease awareness function, we can give only a small selection here. Misra et al. 22 consider a simple SIS model with aware and unaware individuals. X denotes the unaware susceptible classes, X m the aware susceptibles, Y the number of infected individuals, and M the cumulative density of media programs. Unaware uninfected individuals catch disease at rate XY and aware susceptibles at rate XM, where and are constants. The cumulative density of media awareness is modeled as where and 0 are constants. At the end of their infectious period infected individuals return to the aware susceptible class. Samanta et al. 31 consider a more complex SIS model. It is built on the model of Misra et al. 22 but allows unaware susceptibles to become infected and also at the end of their infectious period, infected individuals may become either unaware susceptibles a fraction 1 − p of the time or aware susceptibles, a fraction p of the time. Moreover individuals can move out of the aware uninfected group to the aware infected group. Greenhalgh et al. 23 further build on the model of Samanta et al.. 31 Instead of the unaware susceptibles X − becoming aware at rate X − M, they become aware at rate where k is a constant. Similarly aware susceptibles X + (t) become infected at rate where 1 is a constant. Disease awareness programs can have applications in other areas too. For example, Ma et al. 32 modeled alcoholism using a mathematical model with a time delay and awareness, using two types of individuals, aware and unaware and modeling the media awareness as a separate variable. Lastly the advent of COVID-19 has focused our attention on how people modify their behavior when there is a threat from infectious disease Musa et al. 33 suggest an epidemic model using disease awareness programs which split the population into aware and unaware individuals for COVID-19 transmission in Nigeria. They fit the model to Nigerian COVID-19 data and assess the impact of disease awareness programs on disease transmission. They explain the effect of awareness programs with regards to the basic reproduction number.
Most modern papers tend to use the more sophisticated approach of dividing the population into aware and unaware individuals. However the simpler approach is historically important and is still a crude way to model media awareness.

Research gap and contribution of the present article
However, there appears to be a research gap here because as far as we are aware no-one has yet attempted to apply disease awareness programs to models of HIV and AIDS transmission between PWIDs. This article is an attempt to fill this gap. So in this article, we use the first simpler type of disease awareness program applied to a modified version of Greenhalgh and Hay's model. However because Greenhalgh and Hay unrealistically assumed that PWIDs clean their needles after use we make a minor adjustment so that, more realistically, PWIDs clean their needles before use and then introduce the disease awareness program and examine the behavior of the disease in the resulting model.

Model description
We modify the mathematical model for HIV transmission which has been described by Liang et al., 12 multiplying the disease transmission term by the factor ( ) to represent reduction in the spread of HIV due to awareness programs. The biological parameters are as described in Table 1 adapted from Greenhalgh and Hay. 5 Note that P 1 , P 2 , P 3 , and P 4 ≥ 0 also P 1 + P 2 + P 3 + P 4 = 1. Define Let (t) be the fraction of HIV infected PWIDs at time t and let (t) be the fraction of needles infected at time t. So we introduce the model as follows:

TA B L E 1 Description of parameters
Parameter Definition 1 How fast PWIDs who are susceptible and PWIDs who have disease but are not aware of this fact visit locations where PWIDs share needles. 2 How fast infected PWIDs who are aware that they have the disease visit places where PWIDs share needles.
Chance that the PWID catches disease but the syringe remains uninfected when an initially susceptible PWID injects with an initially uninfected needle.
Chance that the PWID does not catch disease and the needle remains uninfected when an initially susceptible PWID injects with an initially infected needle.
Chance that the PWID catches the disease and the needle remains infected when an initially susceptible PWID injects with an initially infected needle.
The chance that the PWID does not catch the disease and the needle stays infected when an initially susceptible PWID injects with an initially infected needle. 1 Chance that an infected PWID leaves uninfected an initially uninfected syringe. 1 Chance that a PWID with disease leaves uninfected a needle that contained the virus before injection.
Proportion of PWIDs who that clean syringes after using them.
Gallery ratio, where = n∕m, and n is the total number of PWIDs and m is the total number of shared needles. p The chance that PWIDs with disease are aware of being infected.
Rate per PWID at which PWIDS either stop sharing needles or develop full-blown AIDS.
In general, we shall assume that is a positive monotone decreasing function with (0) = 1. We reduce the dimensions of the model in (2) and (3), by supposing that equation (3 ) is at steady state as a similar technique is used in models for HIV among PWIDs as discussed by Liang et al. 12 We do that then give the basic analytical results and simulations. We get that Next, we are going to make the model more realistic by modifying Greenhalgh and Hay's model.

Development of Greenhalgh and Hay's model
We use the notation that if f (x) and g(x) are two functions, then g(x) → 0 as x → 0. The term (n − I) 1 ( )(P 1 + P 3 )(1 − ) Δt is because there are n − I uninfected PWIDs each of whom injects at rate 1 ( ), chooses an infected needle with probability , does not clean the needle before use with probability 1 − and is infected at each injection with probability P 1 + P 3 . Rearranging Letting Δt → 0 Dividing by n, the number of PWIDs, Now we turn to the needle equations. We are going to construct and examine the differential equations for (t) the proportion of PWIDs with disease and (t), the proportion of needles with disease. We construct and examine the differential equations for these quantities. Consider the number of infected needles at time t + Δt. (i) n(1 − p ) PWIDs arrive at shooting galleries at rate 1 ( ). Also np PWIDs visit at a rate 2 ( ). Each PWID chooses one of the m shooting galleries randomly. So at one given shooting gallery PWIDs arrive at rate Here = n∕m is the number of PWIDs divided by the number of needles. So Syringes are infectious at time t and no PWIDs use them in [t, t + Δt).
(ii) For a given shooting gallery infected PWIDs enter it at rate If a PWID who is infected uses a syringe, this syringe will be infectious after being used and cleaned with probability units of injection equipment will be left infectious subsequent to use by a PWID with disease.
(iii) If we take a given syringe uninfected PWIDs come to it at a rate 1 ( )(1 − ). If a susceptible PWID uses an infectious syringe, afterwards the equipment will be capable of transmitting infection with probability (1 − P 1 − P 2 )(1 − ). So the number of syringes infectious after use by an uninfected PWID in [t, t + Δt] is Hence Subtracting i(t) from both sides and dividing by Δt Letting Dividing by m Hence we deduce that Equations (2) and (3) hold where , , , and are redefined as Note that using numbers not fractions of needles and PWIDs with disease Equations (2) and (3) become Note that > . These equations are explained by the flow diagram in Figure 1.
In Equation (15), is the rate at which susceptible PWIDs arrive at infected needles and become infectious. On the other hand I is how fast infected PWIDs stop sharing needles. Equation (16) is more complicated. The term (15) and (16) is the rate at which infected PWIDs arrive at uninfected and infect the needle, that is, the rate at which new infected needles occur. Of the two terms on the right-hand side is the rate at which uninfected PWIDs visit infected needles and either clean the needle before use or flush the needle, in other words the rate at which uninfected PWIDs visit infected needles and leave the needle uninfected. The other term on the right-hand side is and (17) is the rate at which infected PWIDs visit infected needles and leave the needles uninfected. This completes our interpretation of Equations (15) and (16). Now we move on to compute the basic reproductive number of equation model (4).

2.3
The basic reproductive number R 0 R 0 is important in epidemic models. Usually disease becomes extinct if R 0 < 1 and takes off if R 0 > 1 so to derive R 0 , we consider a new infected PWID coming into a steady state population with no disease. The basic reproduction number is the average number of PWIDs who catch the disease via only one infectious syringe. The definition of R 0 used here is similar to the definition used in Macdonald, 34 Massad et al., 35 Sanches and Massad, 36 and Van den Driessche. 37 It takes the basic infectious unit to be an infectious human. This method and the corresponding definition of R 0 are different than the next generation matrix approach as discussed by Diekmann et al., 38 Van den Driessche, 37 Van den Driessche and Watmough, 39,40

and Roberts and
Heesterbeek. 41 The next generation matrix method treats PWIDs and syringes each as separate infectious entities. Consequently, the value of R 0 derived by this method is the square root of the one which we obtained. However as each passes through one at the same time they give equivalent qualitative results if used as a threshold value.
R 0 is given as This expression for R 0 can be derived by considering a single PWID who is infected with HIV in a completely disease-free environment, when no other PWIDs have HIV and all needles are clean. From Equations (2) and (3), we have As we are near the disease-free equilibrium (DFE) we neglect second order terms in small quantities to obtain A newly infected PWID remains in the sharing injecting population for time 1  ( ) = secondary infections in PWIDs, so R 0 = . R 0 can also be derived by considering the expected number of syringes infected via only one infected PWID caused by a single syringe which has just been infected entering a steady-state population with no disease. This syringe causes infectious PWIDs who each in turn infect infectious needles. So again R 0 = ( ) ⋅ ( ) = . We will see that R 0 is a critical parameter which will determine if the disease can sustain itself.
This concludes our analysis of the basic reproduction number. In the next section we shall analytically study the one dimensional model Equation (4).

ANALYSIS OF THE ONE-DIMENSIONAL MODEL
We determined the dynamical behavior of the model in (4) depending on the basic reproductive number we started off by applying the concept of a Lipschitz continuous function and the Picard-Lindelöf theorem to the model (4) to state that there is one and only one non-negative solution, then we discussed the existence of equilibrium solutions and their stability.

Unique non-negative solution
To show that there is one and only one non-negative model solution of the model (4), we require to apply the concept of a Lipschitz continuous functions and the Picard-Lindelöf theorem.
1. Lipschitz continuous functions: 42,43 Definition: Let (X, dX) and (Y , dY ) be two metric spaces as described in Choudhary 44 where dX denotes the metric on the set X and dY is the metric on set Y , a function f ∶ X → Y is called a Lipschitz continuous function if there exists a real constant K ≥ 0 such that for all x 1 and K is called a Lipschitz constant for the function f . In particular, for a real-valued function define Y as the set of real numbers of R with the metric dY (y 1 , y 2 ) = |y 1 − y 2 |, and X might be a subset of R with the same metric. 2. The Picard-Lindelöf theorem: The Picard-Lindelöf existence theorem is an important theorem in the study of differential equations, which indicates existence and uniqueness of solutions to first-order equations with given initial conditions. Consider the initial value problem Suppose that f ∶ R × R → R is uniformly continuous in y. This means that the Lipschitz constant K is independent of t. Then for some >0 there exists a unique solution for y(t) to the initial value problem in the interval  Proof. See the Appendix. ▪ Next, we shall examine how many equilibria there are for the above model (4).

Existence of equilibria
We shall show that if R 0 is less than one then there is only the steady state with no disease whereas if R 0 exceeds one there is one and only one steady state with disease present. If R 0 ≤ 1 then there is only the DFE whereas if R 0 > 1 then there is a unique endemic equilibrium as well as the DFE. We shall first look at the case where is strictly decreases with , and then the case where it is just (possibly not strictly) monotone decreasing.

Theorem 2. Suppose that
is strictly monotone decreasing and R 0 ≤ 1 then Equation (4) will have one and only one steady-state solution where eventually there is no disease present where the disease dies out in PWIDs, * = 0. This is the only equilibrium. If R 0 > 1 there exists exactly one nonzero steady-state solution * > 0 in (0,1] as well as the DFE.
Proof. The trivial equilibrium is * = 0, and any other equilibrium must satisfy the equation Rearranging (21), we deduce that Define g 1 ( ) = 1 and g 2 ( ) = 1 There are several situations to consider.
For ∈ [ * * 1 , * * 2 ], ( ) = ( * * 1 ), g 1 ( ) < ∞ and g 2 ( ) = g 2 ( * * 1 ) = ∞. On the other hand for ∈ ( * * 2 , 1], g 1 ≥ 1 and ( ) < , so g 2 ( ) < 1. So there are no roots of Equation (22)  We have shown that if R 0 is less than or equal to one then there is only the steady state with no disease present whereas if R 0 exceeds one then there is the DFE and a unique steady state with disease present (denoted the endemic equilibrium (EE)). We shall now explore the local stability of the equilibrium. Proof. To study local stability of the equilibrium we consider whether if is slightly displaced from the equilibrium point * it will return to it or move away. We can write

Local stability of equilibrium
Now note that Hence if R 0 < 1 the DFE is locally asymptotically stable. If R 0 = 1 then the DFE is neutrally stable. If R 0 > 1 then the DFE is unstable.
As both terms are negative, df d | | | = * < 0 and the unique EE is locally asymptotically stable when it exists. This completes the proof of Theorem 3. We shall now proceed to look at the global behavior of the system. ▪

Global stability of equilibria
Theorem 4. Suppose that ( ) is monotone decreasing in . If (0) = 0 then (t) = 0 for all time. If R 0 ≤ 1 then the disease will always die out whatever the initial fraction of PWIDs infected, so we have global stability of the DFE. If R 0 > 1 and disease is initially present then over a long time the solution to Equation (4) approaches the unique endemic equilibrium.
Proof. See the Appendix. ▪ Hence we have shown that if R 0 is less than or equal to one disease will become extinct whatever the starting value. If there is no disease initially then there will never be any disease. If there is initially disease and R 0 > 1 then the solutions will tend to the unique steady state with disease present for large time. So in particular limit cycle solutions cannot exist.

Natural interpretation of results
The above analysis has been rather technical. In layman's terms R 0 is a critical value that determines if infected individuals are able to infect enough other individuals to sustain the disease. If R 0 ≤ 1 this is not the case and the infection will go extinct. Conversely if R 0 > 1 each of the number of PWIDs with disease and the number of needles infected approach a unique strictly positive level for a large time. So disease will ultimately always be present in both needles and PWIDs at these levels. Next, we are going to show some numerical simulations and confirm our theoretical analysis results.

SIMULATIONS
We support our analytical results given in Theorems 1-4 and Corollary 1 by numerical simulations. Our simulations were performed using MATLAB and the numerical ordinary differential equation solver (ode45). Our computer program was tested using comprehensive output from a large number of runs. Throughout various simulations we have used realistic parameter values for HIV among PWIDs but our main objective is to verify the analytic results which estimate the spread of HIV among PWIDs for model (4) with two disease awareness programs. We showed that if R 0 is less than or equal to one then the infection will go extinct and if infection is there initially and R 0 is greater than one then disease will approach the one and only one steady state with disease present. Motivated by Li et al. 21 and Cui et al., 27 we take two functional forms for ( ). The first one is ( ) = , where a and b are positive constants with 0 ≤ a ≤ 1, and the second is ( ) = e −m 0 n where m 0 is constant and n represents the number of the PWIDs population. We shall make similar assumptions as in Liang et al.. 12 We shall take p = 0 and assume that 1 = 2 so that PWIDs who know that they have the disease share needles at the same rate as those who are uninfected and those who are infected but do not know it. Also we take 1 , the probability that a PWID who had disease at the start clears virus from a syringe that was uninfected at the start, and 1 , the probability that a PWID who had disease at the start leaves uninfected a needle that contained virus at the start, to be zero as these probabilities are very small and in simple models of spread of HIV among PWIDs these probabilities are normally taken as zero.
We studied the behavior of the model of Equation (4) through altering R 0 by choosing different values of . The starting value was initially (0) = 1 in all cases. Figure 3 shows the plots of six simulations with disease awareness program ( ) = 1 − a b+ (taken from Li et al. 21 ), with different values of the constants a and b are shown in Figure 3. In Figure 3A,C,E, where R 0 > 1, we choose = 0.0, then from Equations (14) and (18) Figure 3, we have that R 0 > 1 for Figure 4A,C,E and R 0 < 1 for Figure 4B,D,F. So in this case as we can see from the figures on the left-hand side of Figure 4 the number of each of infected PWIDs and infected needles will both tend to zero.
In Figure 5, by using both the disease awareness program functions with the same parameters in Figures 3A,B and  4A,B, we considered five different initial values (0) ∈[0,1] of the infected PWID population who do not clean their needles before use. The cases with = 0.0 are given in Figure 5A,C with the same results as previously, we observed that if the PWID population do not clean their needles before use ( = 0.0) then R 0 > 1 and this means that over a long time the fraction of PWID population which was HIV infected tended to the unique endemic equilibrium. For the other two cases (given in in the Figure 5B,D) with = 0.7, the PWIDs often cleaned their needles successfully before use, so R 0 < 1, then the HIV virus died out after a long period of time in both PWIDs and needles. Similar to the results of Figure 4, we have that R 0 > 1 for Figure 5A,C. So in this case again the numbers of each of infected PWIDs and infected needles tends to a unique steady state level. For Figure 5B,D, R 0 < 1 and the levels of infection in both PWIDs and needles tend to zero.
We did other simulations with a variety of other starting values and a variety of other model parameters, and in each case the results of Theorems 1-4 and Corollary 1 were verified. Some of these simulations are presented in Figures A1 and  A2 in the Appendix. For R 0 ≤ 1 the disease always dies out whatever the starting values, whereas for R 0 exceeding one and infection initially in the system the infection tends to the unique steady state. So again the numerical results confirm the analytical result that limit cycle solutions cannot exist.
It is of interest to explore how the other parameters (other than R 0 ) intervene in the transmission dynamics in order to contribute to the discussion of the transmission control. But really this is a very complex question and we do not have space to do this properly here. Also we suspect that at different points in the parameter space, different parameters will have different relative effects on the transmission. We did some preliminary experiments by just varying one of the parameters of the model and keeping the other parameters fixed. For the parameter values used it looked as though only 1 , 1 and made a significant difference to the transmission dynamics and equilibrium level of HIV among PWIDs (see Figure A3 in the Appendix for the disease awareness function ( ) = e −m 0 n ). Other parameters such as 2 , , 1 , and p appeared to make little difference.

CONCLUSION
In our work, we have constructed a differential equation model of the effect of the awareness of levels of infection on the transmission of HIV between PWIDs. This work used as a basis the work of Greenhalgh and Hay 5 and Liang et al. 12 We deduced the equations to describe the improved model that reduces the spread of the diseases through the effect of awareness of disease on sharing needles and syringes among the PWID population. We performed a steady state and stability analysis for our model. Our discussion has been focused on the effect of PWIDs adjusting their behavior to known levels of disease. The key biological parameter of our model is the basic reproductive number R 0 . We deduce that the basic reproduction number is all that we need to know to decide what the model will do. If R 0 is less than or equal to one and is decreasing in then the system has a unique disease-free steady state which the system will tend to over time. If there is no disease initially present then there will never be any disease. If R 0 exceeds one then there is the DFE and additionally a unique endemic equilibrium. If there is disease initially present and R 0 exceeds one then the level of disease tends to the unique steady state. We also showed that the DFE is locally asymptotically stable if R 0 is less than one, neutrally stable if R 0 = 1 and unstable if R 0 exceeds one. For the latter case we also showed local asymptotic stability of the steady state with disease present as well.
We performed numerical simulation on the Equation (4), describing the effect of awareness programs on reducing the spread of HIV among PWIDs. We started off with realistic parameters taken from a literature review, and we assumed that the visiting rate of the shooting gallery is the same ( 1 = 2 ) for both susceptible PWIDs and the infected PWIDs whether or not they know that they are infected. The simulations were divided to simulate two disease awareness programs by changing the constants in these awareness programs for each one. Also we calculated the result of the basic reproductive number and we simulated the total of proportion of PWID population infected over time. This was calculated for different values of the basic reproduction number R 0 which was changed by altering ( ), the proportion of all PWIDs who successfully clean needles before injecting. We repeated this simulation for both awareness programs equations with different initial values of the fraction of PWIDs.
To briefly summarize our results in simple terms R 0 is a parameter that will tell us whether the disease is able to sustain itself in the population and the level of disease in both PWIDs and needles dies out. If R 0 > 1 then it is possible for the disease to sustain itself in the population and the level of disease in both PWIDs and needles tends to a unique nonzero level.
This concludes our analysis of the one-dimensional system given by (14). However recall that Equation (14) was obtained as an approximation of more realistic two-dimensional model (2) and (3) by realistically assuming that PWIDs inject quickly compared with the epidemiological changes which are slower. Hence a topic for future research is to introduce a disease awareness program into the more realistic two dimensional model. We note that as well as reducing disease transmission, some interventions may take place. Another area for future research is how to incorporate interventions into the model. This could be done by dividing the individuals with HIV into asymptomatic and symptomatic individuals and treating symptomatic individuals as in Salman 20 or dividing these individuals into more classes, including those on PrEP and antiretroviral therapy as in Jing et al., 17 in order to model these treatments. There are many other papers modeling HIV interventions, including Zhang et al., 15 Babaei et al., 13 Gountas et al., 18 Des Jarlais et al., 14 and Flountzi et al. 45 and these give other options for modeling HIV interventions. To our knowledge our current article makes a research contribution in that it is the first time that disease awareness programs have been applied to mathematical models of HIV among PWIDs.

ACKNOWLEDGMENTS
Maha Alsharari would like to thank the Ministry of Education in the Kingdom of Saudi Arabia and the University of Tabuk for their funding as scholarship support (TBU258).

CONFLICT OF INTEREST
The authors declare no potential conflict of interest.

DATA AVAILABILITY STATEMENT
Data sharing is not applicable to this article as no new data were created or analyzed in this study.  (4) where , , , and are defined by Equation (14). It is straightforward to show that the function f is Lipschitz continuous in [0,1]. We shall split the proof into three different cases. The first one is (0) ∈ (0, 1), the second one is (0) = 1 and the third one is (0) = 0. First, suppose that (0) ∈ (0, 1). By applying the Picard-Lindelöf theorem there exists a unique local solution.

ORCID
Let us define [0, e ) to be the maximum interval where a solution exists and ∈ (0, 1) for all in [0, e ). We shall show that e = ∞ by using an argument by contradiction.
We have now finished the proof of Theorem 1. ▪ Proof of Theorem 4. We have shown that there is always a DFE possible which is locally asymptotically stable if R 0 < 1 and unstable if R 0 > 1. In the latter case there is a unique EE which is locally asymptotically stable. It is clear that (0) = 0 implies that (t) = 0 for all time.
If < pick > 0 such that ≤ min