{'fraction_non_alphanumeric': 0.02551496670608003, 'fraction_numerical': 0.013099757296658162, 'mean_word_length': 4.757613758509494, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The role of behavioral changes and prompt treatment in the control of STIs

In this paper, we study general recovery functions and treatment in the dynamics of an SIS model for sexually transmitted infections with nonzero partnership length. It is shown how partnership dynamics influences the predicted prevalence at the steady state and the basic reproduction number. Sobol's indices are used to evaluate the contribution of model parameters to the overall variance of R0. The recovery functions studied here take into account that society's capacity to provide treatment is limited when the number of infected individuals is large. Bifurcation analysis is used to establish a relationship between an alert level of prevalence and the minimum recovery time that guarantees the eradication of the disease. We also show that a backward bifurcation can occur when there are delays in the treatment of infected individuals.


Introduction
Mathematical models are an important tool for understanding the spread of diseases in populations. A better sense of the transmission characteristics of infectious diseases can lead to an improvement of our capacity to prevent and control them. Therefore, in recent years, there has been a growing interest in the field of mathematical epidemiology. In particular, many models have been proposed to examine the transmission of sexually transmitted infections (STIs), and the impact of various control policies (Chen & Ghani, 2010;Dietz & Hadeler, 1988;Heijne et al., 2011Heijne et al., , 2013Kretzschmar & Dietz, 1998;Kretzschmar & Heijne, 2017;Muller & Bauch, 2010).
Traditional models for STIs usually overlooked the existence of sexual partnerships. This is partly because they assume homogeneous mixing among individuals. In other words, these models suppose that each sexual contact occurs with a random person; therefore, the entire population is all the time at risk of contracting the infection. This assumption may be comprehensible when modeling a group of highly promiscuous individuals but is not reasonable for the average population. In fact, in real life scenarios, sexual partnerships usually have nonzero length; and there is also a positive time gap between partnerships (Foxman, Newman, Percha, Holmes, & Aral, 2006;Muller & Bauch, 2010).
The explicit inclusion of sexual partnerships is important to address the fact that transmission only takes place when a susceptible individual and an infected one form a sexual partnership. Therefore, as remarked in (Dietz & Hadeler, 1988;Muller & Bauch, 2010), two susceptible individuals that form a pair can be considered temporarily immune as long as they do not separate. For a pair of two infected individuals, if the partnership lasts enough they can clear the infection before the partnership ends and avoid transmission to future partners. These examples show that transmission in a model that includes sexual partnerships (pair model) could be slower than in a classical homogeneous mixing model.
Another interesting subject of research that has been shown to bring rich dynamics in epidemiological models (Villavicencio-Pulido, Barradas, & Luna, 2017;Wang, 2006;Wang, Liu, Zheng, & Takeuchi, 2012;Zhang & Liu, 2008) is the role played by the recovery function on the spread of infections. Assuming a recovery rate proportional to the size of the infective class is not satisfactory due to the fact that the health-care system is limited when the number of infected individuals is large. Incorporating a general recovery function is necessary to evaluate situations like these. Moreover, a general recovery function can also describe different control policies and behavioral changes in the population related to the number of infectious individuals.
The present work includes pair formation considering nonzero length partnerships as an explicit variable of the model and capturing the dynamics of partnership formation and dissolution. The pair model proposed below generalizes the work presented in Kretzschmar and Dietz (1998) to an SIS (susceptible-infectious-susceptible) structure. The new model incorporates a general recovery function that allows studying situations in which, for instance, behavioral changes occur or medical resources are limited. We shall analyze how the recovery function affects the equilibrium level of prevalence and other possible effects in the context of the persistence of the disease in the population.
The structure of this paper is as follows. In the next section, we describe the basic assumptions of the model and calculate an expression for the basic reproduction number. In section 3, a global sensibility analysis is made to identify the parameters or combinations of parameters that contribute the most to the variance of R 0 . In section 4, we propose and study two different recovery functions. Bifurcation analysis is used as a tool to explore the structure of the model's solution for a plausible range of parameter values. The role of the pair formation process is analyzed in section 5. Finally, the conclusions and a discussion are presented in section 6.

The model
The compartmental SIS model proposed below describes the pair formation process and the spread of an infection within partnerships. Single individuals form monogamous sexual partnerships at a constant rate r per unit of time and the partnerships break up at a rate s per unit of time. Transmission of the infection can only occur within partnerships, with f being the number of sex acts per unit of time and h the transmission probability per contact.
As pointed out in (Dietz & Hadeler, 1988), if the pair formation process is considered only as a social act and sexual contacts occur within a pair with a certain rate, for a high separation rate there could be partnerships without a sexual contact during their existence. These pairs would be irrelevant for the infection process. Hence, we shall assume that every partnership starts with a sexual contact.
Individuals enter the sexually active phase of their lives at a constant rate n as singles and leave the sexually active population at a rate m. Moreover, for simplicity, we are going to omit the relation between the infectious disease and the pair formation process. Therefore, being infected does not bias individual's tendency neither to form or break partnerships, nor to have sexual contacts (Kretzschmar & Dietz, 1998).
The model equations resulting from these assumptions are: where IðtÞ ¼ X 1 ðtÞ þ P 01 ðtÞ þ 2P 11 ðtÞ is the total prevalence at time t. Table 1 summarizes the model variables. System (1) is an extension of the pair model proposed by Kretzschmar and Heijne (Kretzschmar & Heijne, 2017). The main difference is that we consider a density-dependent recovery function FðIÞ instead of a constant recovery rate. We also assumed that every partnership starts with a sexual contact which is not the case in Kretzschmar and Heijne (2017). Since the recovery function is directly involved in the transmission of the disease, the inclusion of the recovery function FðIÞ has a major impact on the dynamics of the model. For instance, as it is shown in section 4.2, allowing FðIÞ to be density-dependent may induce the appearance of a backward bifurcation which is an important phenomenon in the context of disease control in epidemiology (Villavicencio-Pulido et al., 2017;Wang, 2006;Wang et al., 2012).
For convenience, we shall assume that the total population size N ¼ X þ 2P is constant with a value of N ¼ n=m.
Assuming that the total resources are fixed and all infected individuals have the same chance to be treated; the likelihood of each infected individual to be treated depends on the density of infected individuals. Therefore, the recovery rate is represented by a function FðIÞ, such that Fð0Þ > 0, and FðIÞ ! 0 for all I > 0. This last condition implies that there is always a nonzero recovery rate; for example, because of the natural immune response or a permanent health care program (Villavicencio-Pulido et al., 2017). Notice that although Fð0Þ > 0, the number of individuals recovering per unit of time is FðIÞI; therefore, for I ¼ 0 there are no recoveries.
Different forms for FðIÞ have been proposed in the literature to model different scenarios. For example, in Villavicencio-Pulido et al. (2017) studied a model where the recovery function is a Michaelis-Menten equation which corresponds to a nonconvex saturation function. They also considered an exponential recovery function whose effect caused the appearance of a backward bifurcation.
Researchers have also explored treatment functions in epidemiological models (Eckalbar & Eckalbar, 2011;Wang, 2006;Wang et al., 2012;Zhang & Liu, 2008). The idea is that public health authorities will mobilize their resources to fight perceived infections depending on the level of prevalence. For instance, in (Eckalbar & Eckalbar, 2011) a quadratic treatment function was proposed to model the fact that society's capacity for providing treatment is limited and can decline after critical equipment and supplies are exhausted or health-care workers fall victim to the disease. In this work, we shall analyze two different forms for the recovery function related to particular control policies to understand how recovery influences the equilibrium level of prevalence.
Under the assumption that the number of pairs is at equilibrium; model (1) can be reduced to the following three dimensional system in terms of the proportions are the proportions of singles and pairs at the equilibrium level of partnership dynamics. To avoid clumsy notation, from now on we are going to omit the bars in the variables of the model (2).
A summary of baseline values, units and interpretation of the parameters is given in Table 2. Behavioral parameters were estimated under the assumption that individuals have on average 1.5 new partners per year and that 70% of the population is in a partnership at the steady state, based on published data (Heijne et al., 2011;Johnson et al., 2001). In addition, we are going to consider human papillomavirus (HPV) infection as an example of a curable STI and parameters consistent with that.
The basic reproduction number (Diekmann, Dietz, & Heesterbeek, 1991;Kretzschmar & Dietz, 1998) for model (2) is given by where Fð0Þ ¼ g represents the infection clearance rate under normal conditions. These conditions might include the natural immune response and any other permanent health policy that influences the recovery time of a single infected individual.
It is worth mentioning that when the recovery function FðIÞ is assumed constant, then the infection clearance rate, and consequently the mean recovery time, depends only on the biology of the pair pathogen-human and permanent health policies. In particular, no additional control interventions depending on the number of infected individuals are included.
Allowing the recovery function FðIÞ to change with I takes into consideration behavioral changes in the population and modifications in control policies due to a perceived hazard: the level of prevalence of the infection. 1=FðIÞ can be interpreted as the mean recovery time when the level of prevalence is I. In particular, 1=g is the recovery time when no additional densitydependent treatment or actions are taken by the infectious individual or society's health sectors. The basic reproduction number R 0 for pair models is defined as the expected number of secondary infections one typical infectious individual will produce during his/her infectious period starting in a P 11 partnership in a completely susceptible population (Heijne et al., 2013). A point to consider when computing R 0 for pair models with nonzero recovery rate is that infected individuals can clear their infection before the partnership ends, but they can get reinfected if their partner is still infectious. These reinfections should be included in the computation of R 0 .
To derive the expression (4) for the basic reproduction number we need the following components: (i) The probability that the initial infectious individual (initial case) is still infectious when separating from a partner.
(ii) The probability that the initial case is still infectious when he/she forms a new partnership.
(iii) The transmission probability per partnership.
(iv) The number of new partners during the infectious period and the number of reinfections for each partner.
First, we compute the probability that the initial case is still infected when separating from a P 11 partnership. The initial case can reach the state X 1 by direct separation from the P 11 with a probability (5) or by first passing through P 01 with a probability The probability p t of the initial case still being infected after separation when there are no reinfections is the sum of equations (5) and (6), In addition, the initial case can reach the state X 1 after m loops of clearance-reinfection of his/her partner. For that reason, we need to examine the case that one or more reinfections take place before separation. A reinfection occurs with a probability Thus, the probability that clearance and reinfection happen exactly m times before separation of the partnership is p m r , m2ℕ. Therefore, the probability that at least one reinfection occurs is As a consequence, the probability p s that an individual who started in a P 11 partnership is still infectious after separation is given by Seeing that is the probability that the initial case is still infectious and sexually active when he/she forms a new partnership, and is the transmission probability per partnership, we conclude that (13) is the average number of individuals that the initial case will infect during his/her infectious period in a completely susceptible population. R c is known as the case reproduction number (Heijne et al., 2013). If we add the expected number of reinfections in the starting and the subsequent partnerships of the initial case ðp r =ð1 À p r ÞÞðR c þ 1Þ to the case reproduction number R c , we obtain expression (4) for the basic reproduction number R 0 .
Another approach to get expression (4) for R 0 is to compute the Jacobian matrix of system (2) at the disease-free equilibrium and pose a condition that guarantees that all the eigenvalues have negative real parts. Under this approach, it is easy to see that R 0 < 1 guarantees the local stability of the disease-free equilibrium. However, R 0 < 1 does not implies neither uniqueness nor global stability of the disease-free equilibrium.

Sobol's indices for R 0
In this section, we will perform a Sobol sensibility analysis (Saltelli et al., 1999(Saltelli et al., , 2004 to evaluate the relative contribution of each individual parameter, as well as the interactions among parameters to the overall variance of the basic reproduction number (4). This will allow us to identify the parameters or combinations of parameters that influence the most the value of the basic reproduction number. This is important in order to plan appropriate control policies.
Due to the uncertainty in parameters' values that appear in the definition of R 0 , we are going to explore plausible ranges for them. The assumptions made for the sensitivity analysis and the corresponding ranges for the parameters are listed below: The number of new partners varies from 1 up to 4 per year. The percentage of people in a partnership at the equilibrium level ranges from 50% to a maximum of 80%.
The pair formation rate is between 2=year and 20=year. The separation rate, s, lies between 1:027=year and 7:77=year. (Equation (3) i.e the fraction of singles and pairs at the steady state was used to estimate the ranges for r and s). The results shown in Fig. 1 indicate that the dominant parameter contributing with about 35% of the variability of R 0 is g, the recovery rate under normal conditions. The transmission probability per contact, h, is also a very influential parameter with a first-order index of 0.2449. The influence of the behavioral parameters r and s on R 0 is smaller but not negligible. This is because although their first-order indices are not very high, the sum of their total-order sensibility indices is above 0.30. This suggests that both r, and s have strong compound interactions with the remaining parameters. On the other hand, the contact frequency within partnerships, f, has only a weak influence on R 0 . The dark marks on top of the bars in Fig. 1 represent 95% confidence intervals for the sensibility indices. Notice that they are very small.
We also computed the second-order sensibility indices to measure the contribution to the variance of R 0 caused by the interaction of two model inputs. Within the ranges explored, the combination of g and s influences the variance of R 0 more than any other combination of two parameters. Our sensibility analysis also indicates that there is a significant interaction between the parameters g and h, and also between r and s. Given the relevance of g and h, we plotted the basic reproduction number (4) as a function of them. As expected, R 0 is a non-increasing function of the recovery rate g. The opposite pattern is seen for R 0 as a function of the transmission probability per contact, i.e. R 0 is a non-decreasing function of h (see Fig. 2).
These results suggest that public health efforts should focus primarily on increasing g, the recovery rate under normal conditions. This can be achieved, for example, through a permanent program of screening, diagnosis, and treatment of cases.
An improvement of clinical services and training of health personnel can also help to increase g. In addition, a reduction of h, the transmission probability per contact, is likewise essential to eradicate the disease. Among the different ways of doing that, barrier methods stand out because they are relatively low-cost, accessible and effective in reducing h.
Finally, the control of the parameters r and s can also be significant to reduce the transmission of the infection. For this reason, efforts to maintain public awareness of STIs and health education are an essential component to control them. These behavioral interventions are useful for reducing individuals' risk of contracting and transmitting STIs. Although in real life situations these strategies are often difficult to implement, their benefits can be considerable. Among the expected effects are a reduced number of sexual partners, delayed sexual debut and mutual monogamy.

Equilibrium prevalence versus recovery function
In this section, we will perform numerical bifurcation analysis of the pair model (2). We shall study through numerical examples the relationship between the prevalence at the equilibrium level and different forms of the recovery function corresponding to various scenarios.

A sigmoid recovery function
Here, the nonlinear recovery function which depends on treatment FðIÞ takes the form of a sigmoid curve, The quotient 1=M is the minimum time needed to recover from the disease when undergoing treatment (therefore g < M). I 0 represents an alert level of prevalence after which the response of the public health authorities begins to grow faster until saturation begins. k models the speed of resource allocation per new infected case. For example, when k is big, the reaction of the system is almost negligible for I < I 0 and, nearest to the maximum capacity for I > I 0 . On the other hand, when k is small, the system takes more into account the gradual increase or decrease in the prevalence. Note that despite the fact that there are no estimates of k, given that Fð0Þ ¼ g, the parameter k can be used to calibrate the value of R 0 . From the definition of FðIÞ, it is not difficult to understand the effects of the parameter I 0 . For low values of I 0 , the system tries to control the epidemic since the first cases with a strong response whose growth is weakened over time by the consumption of resources. In contrast, for high values of I 0 , the initial stage of the response is slow but becomes faster when the alert level of prevalence is reached and grows to the carrying capacity of the system. We are interested in finding conditions on the parameters I 0 and M that guarantee the eradication of the disease. This is essential because it will allow us to know how fast we should act to control an outbreak, given the recovery time from the infection. To find these conditions, we shall explore by means of a bifurcation diagram how the structure of the solution of model (2) depends on the alert level of prevalence, I 0 , for a fixed value of M. The parameter I 0 will vary between 0 and 1, while the other parameters are fixed at their baseline values and k ¼ 6:74. In particular, the value of M will be 52 so that the minimum recovery time is close to one week.
The resulting bifurcation diagram is plotted in part (a) of Fig. 3. The bifurcation parameter I 0 is shown on the horizontal axis of the plot and the vertical axis shows the prevalence of system (2) at the steady state. As usual, stable solutions are represented by a solid line and unstable ones with a dotted line. White squares symbolize static bifurcation points. The results indicate that the disease-free equilibrium point begins being stable for small values of I 0 , but loses its stability when I 0 reaches the value 0.6006. At this value, the system has a transcritical bifurcation with bifurcation point 2 and consequently a stable endemic equilibrium appears. This indicates that the alert level of prevalence should be less than 60% of the population to successfully control the disease when the minimum time to recover is roughly one week. Note that in these conditions, the maximum percentage of the population that is infected at the steady state is 35:051%.
The bifurcation point 2 is of paramount importance because it gives the exact relationship between I 0 and M for which the value of the basic reproduction number R 0 is equal to 1. The bifurcation diagram (a) in Fig. 3 shows this relationship for M ¼ 52. But, this can be obtained for any value of M. We proceed to show that, through the continuation in two parameters of the bifurcation point 2. The result is plotted in Fig. 3 part (b), where the horizontal axis corresponds to M and the vertical axis to I 0 . From this figure, we can deduce the alert level of prevalence needed to control the disease according to the minimum recovery time. For example, when the minimum time needed to recover is nearly one month, that is, when M ¼ 12, the alert level of prevalence should be less than 38% to successfully control the disease. This can be found drawing a vertical line at M ¼ 12 in Fig. 3 part (b). Fig. 3. Part (a) shows how the prevalence at the steady state depends on the alert level of prevalence. Part (b) shows a relationship between I 0 and M that guarantees R 0 ¼ 1. Other parameters were fixed with their baseline value (see Table 2), additionally, we fixed k ¼ 6:74.

A saturated treatment function
Treatment is one of the most efficient ways to control the spread of a variety of infectious diseases, see (Wang, 2006) and the references therein. Nevertheless, due to logistic and economic constraints, any community or city has a maximum capacity for treatment of a disease. We are going to considered this phenomenon in our model introducing the saturated treatment function proposed in Zhang and Liu (2008). The treatment function can be interpreted as the product of a recovery function and the prevalence, that is, TðIÞ ¼ FðIÞI. Thus, to include (15) in our model it is enough to assume that FðIÞ ¼ g=ð1 þ aIÞ.
The function TðIÞ describes the effect of delayed treatment when the population of infected individuals is large, and medical resources are limited. This is reflected through the parameter a ! 0, which measure the extent of the effect of there being a delay in the treatment of infected individuals (Wang et al., 2012;Zhang & Liu, 2008). Note that TðIÞ $ gI for small I; therefore, the treatment rate is proportional to the number of infected individuals when the prevalence is low. But, when the fraction of infected individuals is close to 1, it tends to a saturation level, since TðIÞ $ g=ð1 þ aÞ. This appears to be more acceptable than the conventional constant rate.
In order to get a better understanding of how delayed treatment can affect the dynamics of our model, we will perform numerical bifurcation analysis with respect to the parameter g and two different values of a. The first value proposed is a ¼ 0:1; hence, in this case, treatment delays are negligible. For the second value, a ¼ 2, the effect of delays in treatment is more pronounced. We will keep other parameters fixed at their baseline values, see Table 2.
The results can be observed in Fig. 4. In part (a), that is, when there are no treatment delays, the system has a transcritical bifurcation at g ¼ 0:8917 designated by the bifurcation point 2. Hence, in this case, the model presents classical behavior in the sense that the basic reproduction number aptly determines the threshold for disease's eradication. However, the model exhibits a backward bifurcation in the presence of treatment delays, see part (b) in Fig. 4. Therefore, public health authorities should guarantee a minimum level of efficiency for the treatment of infected individuals to avoid the danger that a backward bifurcation represents.

The role of the pair formation process
The explicit inclusion of sexual partnerships in epidemiological models for STIs is necessary because most of these infections are transmitted within a partnership of two individuals who engage in sexual intercourse and have repeated sexual contacts with each other (Heijne et al., 2013). The pair formation process can impact the transmission dynamics in different ways. For example, if the mean time of partnership duration is short, the number of sexual acts within the partnership could not be enough to transmit the infection. On the other hand, if the partnership duration is long, infected individuals can clear the infection before the partnership ends, and the number of partners during the infectious period is low.
In this section, we shall study how partnership duration affects the spread of the disease. In order to do that, we shall analyze the dynamics of the system (2) for different values of the separation rate s. Clearly, the duration of partnerships is inversely related to the value of the separation rate. Consequently, we can explore the formulation of pairs in comparison to what would occur in a non-pair model taking a large value for s. It should be noted that changes in the separation rate when the rest of the parameters are fixed, not only change the duration of partnerships. It also modifies the expected number of partners per unit of time and the total number of pairs. Hence, the rate of pair formation should also be changing to ensure that the percentage of people in a partnership remains constant. If r is fixed, it can be difficult to know if dynamical changes when s increases are due to the duration of partnerships or solely to the decline in the total number of pairs (Heijne et al., 2011;Muller & Bauch, 2010).
Let us now analyze how the structure of the solution depends on the separation rate. We check numerically the stability of the equilibria and the dependence of the basic reproduction number on the parameter s, see Fig. 5. The recovery function used was FðIÞ ¼ g=ð1 þ aIÞ with g ¼ 1, a ¼ 0:1. In the bifurcation diagram (a) in Fig. 5, we varied the parameter s in the interval ½0; 5. In this diagram the value of r changes with s in such a way that the percentage of people in a partnership is constant consistent with published data (Johnson et al., 2001). In part (b), the basic reproduction number is plotted as a function of the mean duration of partnerships which in turn depends on s. In fact, partnership's duration is long for small values of s and short when s is big.
In Fig. 5 part (b) it can be observed that R 0 is a non-increasing function of partnership duration, this implies that non-pair models overestimate the value of the basic reproduction number which is logical since they implicitly assume a highly promiscuous behavior in the population. The bifurcation diagram (a) confirms that something similar happens for the equilibrium prevalence. Since for small values of s (i.e. large partnership duration) the disease-free equilibrium is stable because, among other things, the number of partners in the infectious period is not enough to maintain the infection.
However, as s increases and reaches the value of 2.3836 a forward bifurcation occurs causing the loss of stability of the disease-free equilibrium and giving rise to an endemic equilibrium whose value increases with s.

Discussion
In this work, we have analyzed the relationship between the recovery function and the prevalence in an epidemic SIS model with nonzero partnership length. The main purpose of incorporating general recovery functions was to study the different scenarios that can occur when treating infected individuals. This contemplates, for instance, logistic limitations to treat infected individuals when their number is large or behavioral changes related to the prevalence of the infection.
Due to the uncertainty in parameters' values and the complexity of the expression for the basic reproduction number R 0 , we have performed a sensibility analysis in order to get a better insight of how input parameters influence the variance of R 0 . The results of this analysis suggest that control strategies should center principally on increasing the recovery rate under normal conditions, g, which is natural. However, a non-obvious conclusion of the sensibility analysis is that the parameters r and s have a significant influence on the basic reproduction number, R 0 . Therefore, behavioral interventions that aim to reduce individuals' risk to transmit and contract the infection are relevant to control STIs. Moreover, in section 5, we showed that R 0 is a non-increasing function of partnership duration and that the endemic equilibrium is stable when s is large. This coincides with previous studies that showed that non-pair models overestimate the value of R 0 compared to models where partnerships are explicitly included (Chen & Ghani, 2010;Heijne et al., 2011;Muller & Bauch, 2010).
The inclusion of a general recovery function provided interesting insight into how saturation and delays in treatment affect the dynamics of epidemiological models. First, from the sigmoid recovery function (14), we established a relationship between the inverse of the minimum recovery time, M, and the alert level of prevalence, I 0 , that guarantees R 0 ¼ 1. In practical terms, this relationship determines when public health authorities should act to successfully control disease's transmission given the minimum recovery time 1=M. Second, when the saturated treatment function TðIÞ ¼ gI=ð1 þ aIÞ is used to describe delays in treatment, we showed numerically that a backward bifurcation can occur when such delays are considerable. This is undesirable in terms of control strategies because driving R 0 below 1 is no longer enough to eradicate the disease. Therefore, timely treatment is of paramount importance to avoid the risk that represents a backward bifurcation.
Our model focuses on capturing the pair formation process and the effects of the recovery function. Thus, many natural extensions are possible to improve the realism of the model. For example, we did not take into consideration that sometimes individuals have contacts outside their partnerships. The existence of sexual risk group is another important aspect that should be addressed when modeling STIs. In addition, omitting the relationship between the infectious disease and the pair formation process is not realistic and deserves further studies. Currently, we are working on these extensions and the inclusion of optimal control theory to study the balance between cost and effectiveness of public health interventions.