Introducing pre-exposure prophylaxis to prevent HIV acquisition among men who have sex with men in Sweden: insights from a mathematical pair formation model

Objectives Since 2017, the Public Health Agency of Sweden recommends that pre-exposure prophylaxis (PrEP) for HIV should be offered to high-risk individuals, in particular to men who have sex with men (MSM). The objective of this study is to develop a mathematical model investigating the effect of introducing PrEP to MSM in Sweden. Design A pair formation model, including steady and casual sex partners, is developed to study the impact of introducing PrEP. Two groups are included in the model: sexually high active MSM and sexually low active MSM. Three mixing assumptions between the groups are considered. Setting A gay-friendly MSM HIV/sexually transmitted infection testing clinic in Stockholm, Sweden. This clinic started offering PrEP to MSM in October 2018. Participants The model is calibrated according to detailed sexual behaviour data gathered in 2015 among 403 MSM. Results By targeting sexually high active MSM, a PrEP coverage of 3.5% of the MSM population (10% of all high actives) would result in the long-term HIV prevalence to drop considerably (close to 0%). While targeting only low actives would require a PrEP coverage of 35% for a similar reduction. The main effect of PrEP is the reduced susceptibility, whereas the increased HIV testing rate (every third month) among PrEP users plays a lesser role. Conclusions To create a multifaceted picture of the effects of interventions against HIV, we need models that include the different stages of HIV infection and real-world data on detailed sexual behaviour to calibrate the mathematical models. Our findings conclude that targeting HIV high-risk individuals, within HIV risk populations such as MSM, with PrEP programmes could greatly decrease the long-term HIV prevalence in Sweden. Therefore, risk stratification of individuals is of importance in PrEP implementation programmes, to ensure optimising the effect and cost-effectiveness of such programmes.

S1 Formulation of the model S1.

Pair-formation model
We consider a sexually active same-sex population, where new individuals enter the sexually active population according to an exponential distribution with rate µn and each individual leaves the sexually active population at rate µ. The size of the sexually active population will therefore fluctuate around the value n, which is assumed to be large. Individuals enter the sexually active population without a steady partner. The rate at which an individual who is single enters into a partnership is ρP 0 , where P 0 is the fraction of single individuals in the population. This means that the higher the fraction of single individuals, the higher the pair-formation rate. Individuals can have at most one steady partner at a time and the separation rate for each partnership is denoted σ. Therefore, a partnership lasts for an exponential time with mean duration 1/(σ + 2µ).
The partnership network is assumed to be stable, i.e. the proportion of singles remains at P 0 for all t. We can then express P 0 (and the proportion P 1 = 1 − P 0 of individuals with a partner) in terms of model parameters [S1]: The rate of sexual acts within a partnership is denoted λ. Beside steady partners, individuals may have casual sex partners during steady partnerships as well as during single periods; the rate at which this occurs depends on the partnership status of the individual under consideration. Up to this point, the pair-formation model described is the same as in [S2]. The first extension of the model from [S2] is to allow for individuals to be either low-active or high-active with regards to the number of casual sex partners. In our application, an individual is assumed to be sexually high-active if they have 15 or more sex partners per year. The fractions of sexually high-active and sexually low-active in the population are denoted π h and π l , respectively (π h + π l = 1).
Let α rq ij be the rate an individual with activity degree r ∈ {l = low, h = high} and i ∈ {0, 1} steady partners tries to find a casual sex partner with activity degree q ∈ {l = low, h = high} and j ∈ {0, 1} steady partners. For this attempt to succeed the individual must actually meet an individual with activity degree q and j steady partners, and therefore, the rate of actual casual sex is α rq ij P j π q . For example, a single who is low-active has casual sex with another low-active single at rate α ll 00 P 0 π l , and with a high-active individual in a steady partnership at rate α lh 01 P 1 π h .

S1.2 Model of infection
As explained in the main text, to model an infection on the network we use a so-called SIR compartmental model (for a survey on stochastic SIR models see [S3]). Individuals can either be susceptible (S), infectious in the acute stage (A), infectious in the chronic stage (C), or on ART-treatment (T ). The second extension of the model in [S2] is to allow for these two infectious stages. Once an individual becomes aware of their infection and starts ART-treatment they are interpreted as immune and can no longer transmit infection. The average duration of the acute infection stage is 2.9 months (= 0.24 years) [S4]. Hence, an individual goes from A to C at rate δ a = 1/0.24 years −1 . Given an unprotected sexual contact (in our case anal intercourse) between an infectious and susceptible individual, there is a probability of transmission depending on stage of infection: p A when in the acute stage and p C when in the chronic stage. Therefore, the transmission rate for an infectious individual in the acute stage in a steady partnership with a susceptible individual is p A λ, and the transmission rate in a casual sexual encounter is p A α rq ij . Note that the probabilities p A and p C of transmission are for the unprotected case, in reality some of the intercourses are with condom. Condom use may also differ with steady and casual sex partners.
The third extension to [S2] is that the time until diagnosis and the beginning of successful ART-treatment may depend on the degree of sexual activity. A sexually high-active individual is put on ART-treatment at rate γ h and a sexually low-active at rate γ l . Table S1: Summary of model parameters. The partnership formation model parameters are given in the first part of the table and the parameters connected to the epidemic in the second part.
Partnership parameters n average population size µ rate of leaving the sexually active population ρ partnership formation rate σ separation rate λ rate of sex acts within a steady partnership π h fraction of high-active individuals π l fraction of low-active individuals ξ rate for high-actives to start taking PrEP α rq ij rate for a r-active individual with i steady partners to try to have casual sex with a q-active with j steady partners Epidemic parameters p A probability of infection in one unprotected anal intercourse during the acute infectious stage p C probability of infection in one unprotected anal intercourse during the chronic infectious stage γ h ART-treatment rate for high-actives γ l ART-treatment rate for low-actives γ P ART-treatment rate for high-actives on PrEP δ a rate of going from acute infection stage to the chronic stage The fourth and main extension is to introduce the possibility for a high-active to take PrEP, which dramatically decreases the probability of getting infected with HIV. The rate a high-active initiate PrEP is denoted ξ, and in our model the use of PrEP reduces the per-act probability of infection by 86% [S5]. Note that the rate ξ does not include imperfect PrEP adherence; ξ is the rate for high-actives to initiate PrEP and to achieve its protective effect. A high-active on PrEP is tested and, if HIV-positive, put on ART-treatment at a rate γ P = 1/0.24 years −1 .
To summarise, the model is captured by 30 parameters (20 free parameters): n, µ, ρ, σ, λ, p A , p C , δ a , ξ; the fraction high-actives π h and the fraction and low-actives π l = 1 − π h ; the 16 parameters α rq ij (8 free), where r, q ∈ {l, h} and i, j ∈ {0, 1}; and the three γ P , γ h and γ l , where γ h = 2.349γ l as described later in Section S3.3. We provide an overview of the notation in Table S1 and an illustration of the model in Figure S1. Note that, we could instead allow for low-actives, instead of high-actives, to be offered PrEP. This possibility will be briefly explored to be compared to the effect of targeting HIV high-risk individuals for a PrEP intervention programme. Figure S1: Representation of possible states a low-active (a) and a high-active (b) can be in. Individuals enter the MSM population as singles into the S compartment. High-actives move to S P rEP at a rate ξ whereas low-actives can never start to use PrEP. Susceptible individuals who acquire infection move to the A compartment (acute infection). Individuals in the A compartment can move to the C compartment (chronic infection) at rate δ a or to the T compartment (ART-treatment). The rate an individual moves to the T compartment is γ P for a high-active on PrEP, γ h for a high-active not on PrEP, and γ l for a low-active. Individuals in the T compartment stay there until they leave the sexually active population. The rate of leaving the sexually active population is denoted µ.

S2 Rate of new casual sex partners
From our egocentric data on the sexual behaviour, we can determine if a participant is low-active or highactive and if he was single or in a steady partnership while having casual sex. It is therefore possible to obtain an estimate of the rate that an r−active individual with i steady partners finds new casual sex partners, let us denote this by α r· i· , where r ∈ {h = high, l = low} and i ∈ {0, 1}. These rates (α r· i· ) can be used to express the different α qr ij for different mixing assumptions. Before we give the equations of the observable α r· i· in terms of the sought α qr ij , we will present some symmetry arguments which reduces the number of free parameters. To begin with, we have 16 different parameters α rq ij . For symmetry reasons, when disregarding the activity degrees, the total rate in the population at which singles have casual sex with individuals in a partnership needs to equal the rate at which individuals in partnership have casual sex with singles. Assume that we have a population of size n = n 0 +n 1 , where n 0 = nP 0 is the number of individuals without a steady partner and n 1 = nP 1 is the number with a steady partner. A similar consistency criterion as for the rate of casual contacts disregarding heterogeneity in activity degree must also hold for the model with two activity degrees: the rate low-active singles (n 0 π l in the population) have casual sex with high-active singles must equal the rate high-active singles (n 0 π h in the population) have casual sex with low-active singles, i.e. n 0 π l α lh 00 P 0 π h = n 0 π h α hl 00 P 0 π l , which simplifies to α lh 00 = α hl 00 .

(S2.2)
This system of equations can later together with a proportionate mixing or complete assortativity assumption be solved, the solutions can be found in Section S2.1 and section S2.1, respectively. For the case of an assortativity between the proportionate mixing and complete assortativity, we need further information than the α r· i· provides. With the help of a proxy question on participants' partners sexual behaviour (explained in Section S3.5), we can additionally estimate more detailed rates than α r· i· : the rate that an r−active individual with i steady partners finds new casual sex partners that are q−active, α rq i· . The following hold for these rates This gives us 8 equations with 10 unknowns. The data does not provide information on whether a casual sex partner is single or in a steady partnership. Therefore, we need to make further assumptions concerning the relation (ratio) between the rate of finding a casual sex partner that are single and the rate of finding a casual sex partner that are in a steady partnership. Let us consider one high-active individual, we assume that the rate of casual contact with a high-active in a partnership compared to the rate with a high-active single, is the same regardless if the considered individual is in a partnership or not, i.e. However, it turns out that the two equations are linearly dependent, we therefore need one more equation to be able to solve the system of equations (S2.3). We make the same kind of assumption but for the case when a high-active meets a low-active: for a high-active individual, the rate of casual contact with a low-active in a partnership compared to the rate with a low-active being single is the same regardless if the high-active individual is in a partnership or not . From the consistency criteria given in Equation (S2.1), Equation (S2.4)-(S2.6) the system of Equations (S2.3) can therefore be solved. Let D qr = π q (α qr 0· P 0 + α qr 1· P 1 ), then we can express the rate for a q−active to try to find a r-active as: Note that, by consistency To conclude, the consistency criteria given in Equation (S2.1), that D hl = D lh , together with the assumption given in Equation (S2.4)-(S2.6) implies that we can write α qr ij as a product: where ω qr i = α qr i· D qr .

S2.1 Proportionate Mixing with respect to activity degree
Proportionate mixing with respect to activity degree means that an individual has no preference regarding which type, high or low-active, it has casual sex with. An individual chooses at random of the potential casual sex attempts in the population. The fraction of potential high-active and low-active casual sex partners will not only depend on the sizes of the two groups, but also the rates at which they try to find new casual sex partners. If the sizes of the groups would be equal, someone trying to find a new casual sex partner would by chance meet a high-active more often since the high-active try to find a new casual sex partner more often than the low-active. Proportionate mixing then implies that the rate a low-active single has casual sex with a high-active single will be i.e. the rate a low-active single has casual sex, times the proportion of all casual sex partners that are from high-active singles. In terms of α qr ij we have that The expression for α qr ij can also be found by using that proportionate mixing implies that ω lh i = ω ll i and ω hh i = ω hl i . By dropping the second superscript and simply write ω l i and ω h i , yields that This together with the system of equations (S2.2) gives the above solution for α rq ij .

⇐⇒
For example, the rate a high-active single finds new casual sex partners that also are high-active, but in a partnership, becomes

S2.3 Mixing determined by the proxy question
To obtain the different α qr ij for the case of an assortativity between the proportionate mixing case and complete assortativity, we use the 8 different α qr i· estimated from data via the proxy question (see Table S3) and Equation S2.9 to get the 16 α qr ij , i.e.

S3 Data and parameter estimates
We will here describe the data gathering, calibration of the model and the parameter estimates obtained from the STI-clinic.

S3.1 Data description
The data used in this study was gathered at a gay-friendly STI-testing clinic in Stockholm, Sweden. Collection of data took place between February 2 and December 15, 2015. Participants first reported demographic information and the total number of sex partners during the last 12 months, then the participants were asked to fill in an app-based timeline follow-back (TLFB) questionnaire.
In the TLFB questionnaire participants were asked to mark up to 10 of their most recent sex partners on a 12-month timeline. Participants did themselves label their partners into one of four partnership types: 1) casual unknown sex partner, 2) casual known sex partner, 3) regular sex partner (regular sex partner but not a 'love' relationship), and 4) main sex partner (a loving relationship, e.g. boyfriend/husband). For casual sex partners, a partner was represented by a single point on the interactive timeline, and a steady sex partner was represented by marking the start and end date of the relationship. For each sex partner on the timeline the participants could report: the partnership type 1) to 4); age of partner; frequency of each sex type (oral/anal; receptive/insertive); frequency of condom use; if the sex took place in Sweden or abroad; drug use and transactions in connection to sex with each partner; and if the participant believed the sex partner had other sex partners concurrently. This last question on concurrency is by us here referred to as the proxy question (for activity-degree assortativity).
In total 403 participants completed the TLFB questionnaire, giving detailed information on 2112 different sex partners. However, for a participant to be included in this study the total number of sex partners and the proxy question need to be answered, as explained in the main manuscript, yielding the data-set of this study consisting of 368 participants and 1903 partners.

S3.2 Scaling of the rate of finding a new casual sex partner
Participants reported their total number of sex partners during a year. The maximum number of sex partners of a participant was 250. When dividing the population into a category of high-active (≥ 15 sex partners a year) and one category low-active (< 15 sex partners a year), 124 (33.7%) participants are defined as high-active and 244 (66.3%) are defined as low-active. The mean number of sex partners of high-actives is 33.21 (median 25, sd 32), and the mean number of low-actives is 5.96 (median 5, sd 3.2). The assumption that only the number of casual sex partners is affected by activity-degree is supported by our data: the mean number of steady partners for high-actives and low-actives is 1.37 and 1.39 respectively.
Additionally, the participants gave detailed information on their (up to) 10 most recent of these sex partners. Participants reported what type of partner these 10 most recent sex partners were, either casual or steady, and the timings of these partners on a timeline. As mentioned, a casual sex partner was reported as a cross on the timeline representing the date of sex and a steady sex partner was given by the start date and end date of the relationship. This timeline data is used to estimate, for example, the time until a new casual sex partner. However, the timeline data only considers data on up to 10 casual sex partners, when we in fact know that many participants have more than this. We therefore need to scale the rates of finding a new casual sex partner according to the total number of partners reported by the participants. The distribution of the number of steady sex partners of high-actives and low-actives are shown in Figure  S2, as seen the distributions are similar. The mean number of steady sex partners per year of a high-active is 1.37 (sd 1.38), from which we estimate that the mean number of casual sex partners is µ h = 31.84 (total number -number of steady). The mean number of steady sex partners of a low-active is 1.39 (sd 1.15) and the mean number of casual sex partners is therefore µ l = 4.57. We use µ h and µ l to scale the rates of finding casual sex partners, since the detailed data from where timings of casual sex partners was given only include (up to) the 10 most recent sexual partners. In the detailed data, the mean number of casual sex partners is 5.35 for high-active individuals and 3.12 for low-active individuals. Hence the scaling factor will be µ h /5.35 for high-active individuals and µ l /3.12 for low-active individuals.
Finally, participants reported to have (a mean value of) 1.4 sex acts with each casual sex partner.

S3.3 Time since last HIV-test
In Table S2 the time since the last HIV-test are shown, indicating that sexually high-active participants test themselves more often than sexually low-active participants. The rate to successful ART-treatment is denoted γ h and γ l for high-actives and low-actives, respectively. Calculating ML-estimates of the time since the last HIV-test, assuming an exponential distribution, we find that the testing rate of a high-active is 2.349 higher than the testing rate of a low-active. We will further assume that the same relationship holds for the time to successful ART-treatment, i.e. the rate of initiating successful ART-treatment for an infected highactive is 2.349 times higher than the rate of initiating ART-treatment for an infected low-active. Therefore, we fix γ h = 2.349γ l such that we only need to vary one of the parameters.

S3.4 Distribution of parameters
The time durations in the model, such as the time until finding a new steady sex partner, are assumed to be exponential and are hence specified by their rates. With this assumption we calculate a point estimate and its standard error. In general, if we are looking at something occurring at an exponential rate, α say, then  (αt)). Observing n events during a time t leads to the ML-estimateα = n/t. The variance of the estimate equal V ar (α) = V ar (N/t) = V ar(N )/t 2 = α/t, which leads to a standard error of the estimate of α/t. For example, in calculating the estimate of the rate for a high-active in a steady partnership to find a high-active casual sex partner, we do as follows: find the total time high-active individuals are in a steady partnership (T h 1 ), then find the number of casual sex partners that occur during that time that is with someone that also is high-active (N hh 1 ) and multiply it with the scaling factor from Section S3.2 (N hh 1 × µh 5.35 ). The point estimate is given byα hh 1· = N hh 1 × µh 5.35 /T h 1 and its standard error by s.e.(α hh 1· ) = α hh 1· /T h 1 . For the number of occurrences of anal intercourses (AIs) in a steady partnership, the participants reported the number of acts during a 1-month period. Let m be the number of steady partners among all participants and let a i denote the number of occurrences of AI with partner i = 1, . . . , m. The estimated rate of AI, in units of months, isλ and the standard error is s.e(λ) = λ m .
In the data for casual sex partners, it is recorded if a condom was used (1) or not used (0) during receptive anal intercourse (RAI) and during insertive anal intercourse (IAI). To estimate the condom use in casual contacts we use a Bernoulli assumption and calculate the mean condom use during RAI and IAI, respectively. With a Bernoulli assumption we mean that, in each new casual sex act a condom is used with a probability, p c say, independently of previous sex acts. Then the estimatep c is given by the mean number times a condom was used. The standard error is given by where n here is the number of observations, i.e. the number of casual sex partners where a binary response on condom use was given. For condom use with a steady sex partner, participants could choose from a five-degree scale on how often a condom was used during RAI and during IAI: always (100%), often (75%), half of the times (50%), seldom (25%), and never (0%). Here, the participants did themselves, in a sense, give the mean number of times they used condom with a partner. Assume the data consist of m such steady partners with corresponding responses (y 1 , . . . , y m ). The estimated condom use in steady partnerships,p s , is then the mean of the m reported values on the five-degree scale, and the distribution ofp s is approximated by the normal distribution Where p s is the true expected value of the condom use in steady partnerships and τ is the standard deviation of the condom use. The estimate of τ 2 is 1 The standard error of the estimated condom use with steady sex partners is then given by With the estimated condom use during RAI and during IAI, we calculated the overall mean condom use by taking the weighted average of these two estimates. See Table S3 for condom use estimates and the proportion (weights) of sex acts that are RAI and IAI, from these values we can calculate the overall condom use for steady partners and for casual sex partners

S3.5 Proxy of partners' activity degree
For the 10 most recent sex partners the participants responded to detailed questions. One of the detailed questions was what we referred to as the proxy question: 'Do you think that your sex partner had other sex partners than you during the same time frame that he/she met you?', with the possible answers 1. I know this person had sex with others 2. I think this person had sex with others 3. No, this person only had sex with me 4. I don't know The percentages of the answers to the proxy question are given in Table S4. Table S4: Distribution of proxy variable on partners' activity degree.
One of the questions participants answered concerning their 10 most recent sex partners was whether they believed their partners had other concurrent sex partners. We use this as a proxy for partners activity-degree, as shown in the fourth and last column. With the help of a consistency criterion, Equation (S2.8), we assigned the partners of answer 4 as low-active.

Answer
Of If a participant answered either 1 or 2 for a partner, we will take this as a proxy for that the partner is high-active. If a participant answered 3, we will use this as a proxy for that the partner is low-active. We now need to decide what to do with the 25% of the partners that participants labelled as No. 4 on the proxy question, the partners that participants did not know whether they had other sex partners.
In the total population, consistency requires that the the number of casual sex acts high-actives have with low-actives needs to equal the number of casual sex acts low-actives have with high-actives. This criterion can be written, in terms of rates, as (Equation (S2.8)) With this consistency we get help in determining how to assign the partners labelled as No. 4 on the proxy question. If we simply remove these partners entirely, the left-hand side of Equation (S2.8) becomes 0.18 and the right-hand side 1.93, very far from each other. Hence, to remove the partners of which the participants do not know (No. 4) yields a too big inconsistency. If we instead assign all these partners as low-active, we get that the left-hand side of Equation (S2.8) becomes 2.19 the right-hand side 1.93. This suggests that many of the partners participants labelled as No. 4 should be categorised as low-active.
In Table S5 we show (for the two choices of actions of answer No. 4 "I don't know"), the proportion of high-active individuals' casual sex partners that will be with low-actives and with high-actives, respectively; and the proportion of low-active individuals' casual sex partners that will be with low-actives and with high-actives, respectively. For example, if we assign all partners that participants labelled as No. 4 on the proxy question as low-active, we find that 34.3% of low-active individuals' casual sex partners will be with low-actives, and 65.7% will be with high-actives. Table S5: Consequence of the two assignments of the partners of which participants do not know if the partner did have other partners. The column named Removed means that the partners of participants of which we do not know (answer 4 on the proxy question) were removed, and the column As low-active means we assigned those partners as low-active. First, for the two assignments, we show the proportion of low-active and high-active partners of participants who are low-active, then the same kind of proportion but for participants who are high-active. Then we show the values of the right-hand side (D lh ) and left-hand side D hl of the consistency criterion Equation (

S3.6 Final estimates of the rates of acquiring new casual sex partners
We will now show the final estimates of the rates of finding new casual sex partners, using the three different mixing assumptions with respect to activity degree. We found that the estimates for α qr ij that utilises the proxy question could be written as Equation S2.9 where consistency requires that D hl = D lh (Equation (S2.8)) needs to be fulfilled. Utilising Table S5, we found that assigning all partners that participants answered "I don't know" (answer No. 4) as low-active on the proxy question yielded a value of 2.19 for the left-hand side and a value of 1.93 of the right-hand side of Equation (S2.8), i.e. similar but not equal. We choose to work with the left-hand side, setting D hl = D lh = 2.19 when estimating the different α qr ij , which can be seen in Table S6. To make the comparison between the different mixing assumptions, we also use that D hl = D lh = 2.19 for all mixing assumptions. In Table S6 we show the estimates for proportionate and complete mixing under the assumption that D hl = D lh , where the values in parenthesis are the ones not requiring that D hl = D lh .

13
Supplementary material BMJ Open doi: 10.1136/bmjopen-2019-033852 :e033852. The value of θ that corresponds to these values is 0.141. Doing the same kind of calculations but for α lh and α ll yields the same θ. Note that, there could exist disassortative with respect to activity-degree, however, we disregard from this since it seems unlikely in our application.

S4 Finding the endemic prevalence, a deterministic approximation
Both the basic reproduction number and the endemic level can be obtained by using a deterministic approximation of the stochastic model. Assuming that the population is large, it is then enough to consider expected values of the fraction of individuals that are susceptible, infectious, or recovered to obtain these quantities. In the following section, we explain how we construct the compartments and give the differential equations governed by the possible transitions within the network model. To find the endemic (or equilibrium) prevalence one needs to find the non-trivial steady state of the system of differential equations. The trivial solution is that everyone is susceptible.
In Figure S1 we showed the different infectious states an individual can be in; the four different infectious states are susceptible S, acute infectious A, chronic infectious C, and treated (on ART-treatment) T . We further divide the population into different types, these types specify: if an individual is single or in a partnership, if the individual is low-active or high-active in having casual sex partners, the infectious state of the individual, and the partner's infectious state. We will study the fraction of the population belonging to each type, each individual will therefore contribute with 1/n to the type it belongs to.
The fraction of all individuals that are susceptible, single and r-active is denoted by S r 0 ; the fraction of all individuals that are single, r-active, and infectious in the acute stage is denoted by A r 0 ; the fraction of all individuals that are single, r-active, and infectious in the chronic stage is denoted by C r 0 ; and recovered singles (on ART-treatment) that are r-active is denoted by T r 0 . Let X = {S, A, C, T } be the set of possible states not including PrEP, let X P = {SP, AP, CP } be the possible states when being on PrEP, and let D = {l, h} be the set of possible activity degrees with regards to the casual contacts. Furthermore, let X rq Y denote the fraction of individuals that are r-active of type X ∈ X ∪ X P , with a q-active partner of type Y ∈ X ∪ X P . Note that, this counts each individual in the fraction X rq Y , not each pair. E.g. S rq S is the fraction of all individuals that are susceptible r-active and in a partnership with a q-active susceptible. The reason for taking this individual-based perspective is that the data is individual based. Moreover, the individual-based perspective makes it simpler to extend the model by allowing for more than one steady partner at a time in future work.
Disregarding the use of PrEP, there are in total 8 different single types and 64 types of partnerships (hence 72 equations). Some of these types' fraction must by consistency be equal, namely • the 4 equations X hl X = X lh X where X ∈ X • the 12 equations X qr S = S rq X where X ∈ {A, C, T } and q, r ∈ D, • the 8 equations X qr A = A rq X where X ∈ {C, T } and q, r ∈ D, • the 4 equations T qr C = C rq T where q, r ∈ D, which reduces the number of equations to 44.
Including the use of PrEP for high-actives creates: 3 more single types (denoted SP h 0 , AP h 0 , and CP h 0 ); 48 partnership types between one participant on PrEP and one not on PrEP; and 9 partnership types where both participants in the steady partnership are on PrEP. Hence, introducing PrEP increases the number of types by 60; however, some of the PrEP types' fraction must also by consistency be equal • the 8 equations X qh SP = SP hq X where X ∈ X and q ∈ D, • the 8 equations X qh AP = AP hq X where X ∈ X and q ∈ D, • the 8 equations X qh CP = CP hq X where X ∈ X and q ∈ D, Introducing PrEP increases the number of equations needed to be specified by 33, from 44 to 77. There are some additional facts that will reduce the number of equations. Recall that the fraction without a steady partner is denoted P 0 , the fraction with a steady partner is denoted P 1 , and the fraction high-active individuals and low-active individuals in the population are denoted π h and π l , respectively. The following two constrains concerning singles must hold: The following three constrains for individuals in steady partnerships must hold. (I) The fraction of the population that is low-active in a steady partnership with a low-active is That is, we sum over all possible states X ∈ X , the first low-active individual in the relationship can have, and over all possible states that Y ∈ X , the second low-active individual in the relationship can have. (II) The fraction of the population that is low-active in a steady partnership with a high-active is here we sum over all possible states X ∈ X , the first low-active individual in the relationship can have, and over all possible states that Y ∈ X ∪ X P , the second high-active individual in the relationship can have. The difference from the previous sum is that a high-active individual can be on PrEP. The fraction of the population that is high-active in a steady partnership with a low-active (which is the same as fraction low with high above) is (III) The fraction of the population that is high-active in a steady partnership with a high-active is This leads to that we can reduce the number of equations further, from 77 to 72. Before we show the system of differential equations that describe our model, we define some quantities that will improve readability. Let us write the fraction of the population that is r-active susceptible in a partnership as S r 1 (where 1 refer to having a steady partner), Similarly, we write the fraction r-active acute infectious individuals in a partnership as A r 1 , the fraction r-active chronic infectious individuals in a partnership as C r 1 , and the fraction r-active diagnosed and on ART-treatment in a partnership as T r 1 , For the PrEP states we have Remember that the transmission probability in one sex act in the acute stage is denoted by p A and in the chronic stage by p C . An r-active individual who is single, susceptible, and not on PrEP acquire infection at rate . Similarly, an r-active susceptible not on PrEP with a steady partner acquire infection via casual sex at rate . Introducing PrEP will make some susceptible less likely to acquire infection, let us denote the reduction by ζ. Assuming a susceptible is protected by 86% by PrEP yields a value of ζ = 1 − 0.86 = 0.14. A high-active susceptible that is single and on PrEP will acquire infection via casual sex at rate ζβ h 0 . The model can now be described by a set of 72 differential equations. We will begin by specifying all single state equations, followed by the different partnership state equations. We remind the reader of the parameter definitions that can be found in Table S1.

S4.1 Single states
For the other single states, we have Due to constrain (S4.1), the fraction of the population that is high-active on ART-treatment is equal to In a similar way we get the equations for a low-active single. Note that a low-active never starts to use PrEP; however, we could switch which activity-group is targeted for the intervention.
And the fraction low-active on ART-treatment is equal to T l 0 = π l P 0 − S l 0 + A l 0 + C l 0 .

S4.2 Partnership states
We will specify the partnership state equations in the following order: first all possible combinations of a susceptible not on PrEP, S S , S A , S C , S T , S SP , S AP , S CP ; then all possible combinations of an acute infectious individual not on PrEP that has not previously been specified, A A , A C , A T , A SP , A AP , A CP ; then all possible combinations of a chronic infectious individual not on PrEP that has not previously been specified, C C , C T , C SP , C AP , C CP ; then the treated individuals not previously specified, T T , T SP , T AP , T CP . Then we specify the partnership type equations of individuals on PrEP: for a susceptible on PrEP, SP SP , SP AP , SP CP ; for an acute infectious on PrEP, AP AP , AP CP ; and finally for a chronic infectious on PrEP, CP CP .

S4.2.1 Susceptible not on PrEP
For susceptible individuals not on PrEP in a partnership with another susceptible not on PrEP we have that

For susceptible individuals not on PrEP in a partnership with an acute infectious individual not on PrEP
Note that A rq S = S qr A . For susceptible individuals not on PrEP in a partnership with a chronic infectious individual not on PrEP And additionally, C rq S = S qr C . For susceptible individuals not on PrEP in a partnership with an individual on ART-treatment We also have that T rq S = S qr T . For susceptible individuals not on PrEP with a steady susceptible partner on PrEP we have Note that SP hh S = S hh SP and SP hl S = S lh SP . For susceptible individuals not on PrEP in a partnership with an acute infectious individual on PrEP Note that AP hh S = S hh AP and AP hl S = S lh AP . For susceptible individuals not on PrEP in a partnership with a chronic infectious individual on PrEP Note that CP hh S = S hh CP and CP hl S = S lh CP .

S4.2.2 Acute infectious individuals not on PrEP
For acute infectious individuals not on PrEP in a partnership with another acute infectious not on PrEP we have where additionally A lh A = A hl A . For acute infectious individuals in a partnership with a chronic infectious individual Note that C rq A = A qr C . For acute infectious individuals in a partnership with an individual on ART-treatment Also, T rq A = A qr T . For acute infectious individuals not on PrEP in a partnership with a susceptible individual on PrEP Note that SP hh A = A hh SP and SP hl A = A lh SP . For acute infectious individuals not on PrEP in a partnership with an acute infectious individual on PrEP Note that AP hh A = A hh AP and AP hl A = A lh AP .

For acute infectious individuals not on PrEP in a partnership with a chronic infectious individual on PrEP
Note that CP hh A = A hh CP and CP hl A = A lh CP .

S4.2.3 Chronic infectious not on PrEP
For an individual with a chronic infection not on PrEP in a steady partnership with another individual with a chronic infection where C lh C = C hl C . For an individual with a chronic infection not on PrEP in a steady partnership with an individual under ART-treatment where T rq C = C qr T . For an individual with a chronic infection not on PrEP in a partnership with a susceptible individual on PrEP And additionally SP hh C = C hh SP and SP hl C = C lh SP . For an individual with a chronic infection not on PrEP in a partnership with an acute infectious individual on PrEP Note that AP hh C = C hh AP and AP hl C = C lh AP .

For an individual with a chronic infection not on PrEP in a partnership with a chronic infectious individual on PrEP
where CP hh C = C hh CP and CP hl C = C lh CP .

S4.2.4 Treated individual
For an individual under ART-treatment in a steady partnership with another individual under ART-treatment we have where T lh T = T hl T . For treated individuals in a partnership with a susceptible individual on PrEP We also have that SP hh T = T hh SP and SP hl T = T lh SP . For treated individuals in a partnership with an acute infectious individual on PrEP Also, AP hh T = T hh AP and AP hl T = T lh AP . For treated individuals in a partnership with a chronic infectious individual on PrEP where CP hh T = T hh CP and CP hl T = T lh CP .

S4.2.5 Individuals on PrEP
For susceptible individuals on PrEP with a steady partner on PrEP we have Note that AP hh SP = SP hh AP and that CP hh SP = SP hh CP . For acute infectious individuals on PrEP in a steady partnership with an individual on PrEP Note that CP hh AP = AP hh CP . And very much finally, for an individual with a chronic infection on PrEP in a partnership with a chronic infectious individual on PrEP we have

S5 Additional results
We will in this Section go through some additional results mentioned in the main manuscript.

S5.1 Not distinguishing individuals according to activity degree
Here we examine what happens if we do not divide the population according to active-degree but assume that everyone behaves in the same way regarding the number of casual sex partners and regarding the rate to ART-treatment. For the case when no one yet is on PrEP, we find that R 0 = 1 when the mean time to successful ART-treatment is 3.28 years. A prevalence of 5% is obtained for a mean time to ART-treatment of 3.57 years. In Figure S3 we show the effect of introducing PrEP in this model without high-actives and low-actives. We see that the PrEP coverage in such a population would need to exceed 5% to reach a prevalence close to 0 (in contrast to 3.5% as in the model where we have two activity degrees).  Figure S3: Effect of introducing PrEP in a population that is not separated according to activity degree, but where everyone is assumed to behave the same with regards of finding new casual sex partners.

S5.2 Not including the different infectious stages acute and chronic
In Figure 2b) in the main text, we saw that the reduction in susceptibility due to PrEP had a larger effect in reducing the prevalence than the increased testing rate of those on PrEP. The transmission probability of HIV is much higher in the acute infectious phase, the first 3 months following infection, than in the chronic phase. The reason for the lesser effect of an increased testing rate could be that it misses a large proportion of the acute stage.
To help verify that this is the case, we modified the model to not make a distinction between the acute and chronic stage; to only include one transmission probability during the whole infectious lifetime of an infected individual. We calibrate this transmission probability so that when no one is on PrEP, and the mean time to successful ART-treatment is 1.77 years for high-actives, the prevalence is equal to 5%. This is done to match the set-up of the analysis in Figure 2b). The transmission probability is then 0.0208 for the whole infectious time, instead of 0.1301 for the acute stage and 0.0098 for the chronic stage.
In Figure S4 it is seen that when only one infectious stage is included, the increased testing and diagnosis rate has as equally big impact on the reduction of the prevalence as the reduced susceptibility of PrEP. This implies that the lesser effect of the increased testing rate, that we found in Figure 2b) in the main manuscript, can be assigned to it missing the 3 month long acute stage. Because, when we in this analysis distributed the increased transmission probability of the acute stage over an infected individual's lifetime,

S5.6 PrEP coverage for alternative transmission probabilities
The transmission probabilities for unprotected receptive anal intercourse (URAI) during the acute and chronic stage are taken from the literature [S6] (0.1835 and 0.0138, respectively). To get the transmission probabilities during unprotected insertive anal intercourse (UIAI), we use an estimate of the relationship between the transmission probability between URAI and UIAI. The transmission probability for URAI is 2.39 times larger than the transmission probability for UIAI [S7]. Assuming equally many insertive as receptive acts, the transmission probability during the acute stage was set to p A = 0.1835 × 0.5 + 0.1835 × 0.5/2.39 = 0.1301, and during the chronic stage p C = 0.0138 × 0.5 + 0.0138 × 0.5/2.39 = 0.0098.
We now want to study how robust our conclusion concerning PrEP coverage, to achieve an endemic prevalence close to 0%, is. We do this by altering the two transmission probabilities by setting them to 50% − 150% of their estimated values. For example: with 50% of the transmission probabilities, we have that p A = 0.0651 and p C = 0.0049; with 150% of the transmission probabilities, we have that p A = 0.1952 and p C = 0.0147. With given transmission probabilities p A and p C , we find the mean time to ART-treatment corresponding to a prevalence of 5%. The given mean time to ART-treatment are for high-actives, for lowactives it is 2.35 times larger. With the different set-ups generating a prevalence of 5%, PrEP is introduced to high-active individuals. In Table S7 we show the different set-ups and the PrEP coverage needed to get an endemic prevalence close to 0. Additionally, we alter the two transmission probabilities one at a time in Table S8 and Table S9. In Table S8 we vary p A but let p C stay fixed at 0.0098. In Table S9 we vary p C but let p A stay fixed at 0.1301. We conclude by noting that the results are almost invariant to which set-up is used. Table S7: Robustness of the PrEP coverage to obtain a long-term prevalence close to 0: alteration of p A and p C . Assuming no one is on PrEP, we first find other combinations than the one used in the main manuscript of the transmission probabilities and mean time to ART-treatment for high-actives that generates a prevalence of 5%. With these different scenarios that generates a prevalence of 5%, we then study the needed PrEP coverage to obtain a long-term prevalence close to 0. PrEP coverage for 0 prevalence 3.42% 3.43% 3.44% 3.47% 3.49% 3.52% 3.55% 3.58% 3.61% 3.65% 3.68% Table S8: Robustness of the PrEP coverage to obtain a long-term prevalence close to 0: alteration of p A . Same procedure as in Table S7, but only the transmission probability in the acute stage is altered. The transmission probability in the chronic stage is remained fixed at p C = 0.0098.  Table S9: Robustness of the PrEP coverage to obtain a long-term prevalence close to 0: alteration of p C . Same procedure as in Table S7, but only the transmission probability in the chronic stage is altered. The transmission probability in the acute stage is remained fixed at p A = 0.1301. PrEP coverage for 0 prevalence 3.5% 3.5% 3.51% 3.51% 3.52% 3.52% 3.53% 3.53% 3.53% 3.54% 3.54%

S5.7 Short-term effect of different PrEP strategies
In our model, we assume that a sexually high-active start using PrEP at rate ξ. We then determine the lowest possible rate ξ that yields an equilibrium prevalence of 0% and calculate which PrEP coverage this corresponds to. The lowest PrEP coverage that eventually results in a 0% HIV prevalence is 3.5% of the population. This 'eventually' is a very long time in the future-it would take centuries. If no new HIV cases would occur, it would still take many years before the prevalence reaches 0%; the HIV prevalence would not reach 0% until the last person with HIV dies. However, HIV will effectively disappear when no new infections occur. Remember that, in our model we assume that diagnosis and the beginning of ART-treatment is the same as being uninfectious, and consequently, only individuals with undiagnosed HIV can transmit the infection. Hence, we will here study not only the prevalence but also the percentage undiagnosed HIV cases for different PrEP initiation rates, ξ.
In what follows, we will look at different rates ξ, where all rates result in an HIV prevalence of 0% in the equilibrium steady state. The lowest ξ we look at will therefore corresponds to an equilibrium PrEP coverage of 3.5% of the population (≈ 10% of high-actives). As a starting point, before any high-active accepts PrEP, the prevalence is set to 5% and the percentage undiagnosed HIV cases to 0.21% (the model with ξ = 0 calibrated to data). In the left panel of Figure S8, we show the HIV prevalence (%) at different PrEP initiation rates, ξ, for 50 years after the beginning of a PrEP implementation programme. In the middle panel of Figure S8, we show the percentage of individuals that are infectious and undiagnosed. In the right panel of S8, we show the corresponding percentages of the population that are on PrEP for 50 years after the beginning of the PrEP programme. For the lowest rate ξ that results in an equilibrium HIV prevalence of 0%, we see that after 50 years the PrEP coverage has only had time to reach 2%, but the percentage undiagnosed has more than halved, and the prevalence has dropped from 5% to almost 4%. This can be compared to the scenario with a ξ that results in an equilibrium PrEP coverage of 11%. Then the PrEP coverage has reached a bit over 7% after 50 years, and the percentage undiagnosed HIV cases is only 1/20 of its value before the initiation of a PrEP programme (from 0.21% to 0.01%).
In Figure S9 we study, in more detail, the effects of different PrEP scenarios 10 and 20 years after their initiation. We look at both the prevalence and the percentage undiagnosed HIV cases. After 10 years, if the PrEP coverage has reached 5% (15% of all high-actives), the percentage undiagnosed HIV cases is reduced from 0.21% to 0.14%. Looking at 20 years after a PrEP programmes initiation and where the PrEP coverage has reached 5%, the percentage undiagnosed HIV cases is reduced to the low level of 0.04%. If the PrEP coverage on the other hand has reached almost all high-actives after 10 years, that is 30% of the population, then the percentage undiagnosed HIV cases is reduced to 0.03%. The same percentage of coverage after 20 years yields a percentage of 0.004% infectious and undiagnosed HIV cases; that is, almost no new HIV infections occur.

S5.8 Allowing individuals to change activity-group
The conclusions of the PrEP coverage needed to eventually eliminate HIV from the community (an equilibrium HIV prevalence of 0%) does not vary much if participants are allowed to switch activity-group.
We further extended our main model to ascertain the effect of allowing high-active individuals to become low-active and allowing low-active individuals to become high-active. The fraction high-active and the fraction low-active in the population were held constant during this analysis. This was achieved by introducing a parameter ν governing the switching, then letting π h ν be the rate for one low-active to change to being high-active and letting π l ν be the rate for one high-active to change to being low-active. The number of lowactive MSM in the population is nπ l and the number of high-active MSM is nπ h . Therefore, the total rate for low-actives to switch will be nπ l π h ν and the total rate for high-actives to switch will be nπ h π l ν. Hence, the fractions being high-active and low-active will fluctuate around the values π h and π l . If a high-active on PrEP switch to low-active we additionally assume that this individual stops taking PrEP.
We tested different switching rates ν = 0, 0.01, 0.1, 0.2, 0.5, 0.75, 1. If, for example, ν = 0.1 then one randomly chosen individual will switch on average every tenth year and if ν = 0.5 then one randomly chosen individual will switch on average every second year. The case when ν = 0 corresponds to no switching. For each scenario we calibrate the mean time until ART-treatment so that without anyone on PrEP the prevalence is 5%. As before, the given mean time to ART-treatment is for high-actives, for low-actives it is 2.35 times larger. The results can be seen in Table S10, from there we see that the PrEP coverage needed to eventually eliminate HIV from the community is not that different between the scenarios, it varies between ≈ 2.5% to ≈ 3.5% of the total population. However, when people on PrEP are allowed to become low-active and thereby stop taking PrEP a higher PrEP-initiation rate is needed to obtain a certain PrEP coverage.
In Table S11 we instead show the PrEP coverage needed when low-actives are targeted. We see the same tendency in which scenarios that easiest eliminate HIV as in Table S10. We also find the same kind of conclusion as in the main text-targeting high-actives for PrEP is much more effective than targeting low-actives.