Mathematical models of SIR disease spread with combined non-sexual and sexual transmission routes

The emergence of diseases such as Zika and Ebola has highlighted the need to understand the role of sexual transmission in the spread of diseases with a primarily non-sexual transmission route. In this paper we develop a number of low-dimensional models which are appropriate for a range of assumptions for how a disease will spread if it has sexual transmission through a sexual contact network combined with some other transmission mechanism, such as direct contact or vectors. The equations derived provide exact predictions for the dynamics of the corresponding simulations in the large population limit.


Introduction
The recent emergence of Ebola Mate et al., 2015) and Zika (Foy et al., 2011;Musso et al., 2015) demonstrates that diseases which spread primarily through other means can also have a sexual component to their spread.
Zika is a mosquito-borne virus which can cause birth defects if a pregnant woman is infected (Oliveira Melo et al., 2016;Ventura et al., 2016). Although more is being learned, it appears that Zika causes self-limiting infections (Duffy et al., 2009) and it is likely that infected individuals recover with immunity (Lessler et al., 2016).
Ebola is a directly transmitted disease which causes extreme morbidity and mortality (Baron, McCormick, & Zubeir, 1983;Emond, Evans, Bowen, & Lloyd, 1977;Heymann et al., 1980). It is spread through direct contact with bodily fluids from an infected individual. Individuals who survive appear to gain immunity (McElroy et al., 2015).
For both Zika and Ebola, evidence suggests that sexual transmission is possible Foy et al., 2011;Mate et al., 2015;Musso et al., 2015). Further, it appears that individuals may be able to transmit the viruses sexually even after they appear to no longer be infectious through the usual mechanisms (Mackay, & Arden, 2015;Deen et al., 2015;Nicastri et al., 2016).
To better predict the spread of such diseases, we need to capture these different transmission modes into a model. Unfortunately, many existing mathematical models of disease spread through networks require a large number of equations (Kiss, Miller, & Simon, 2017;Miller & Kiss, 2014). However, in the specific case of SusceptibleeInfectedeRecovered (SIR) diseases such as Ebola and Zika, a low dimensional model, the Edge-based compartmental model, exists which can capture a diverse set of assumptions about the network structure (Miller, Slim, & Volz, 2012;Miller & Volz, 2013). However, it is structurally very different from the usual models used for other transmission mechanisms. Consequently it is not immediately obvious that we can combine the different transmission models into a single low-dimensional mathematical model.
Our goal in this paper is to develop simple mathematical models which capture a range of different transmission mechanisms combined with a sexual mode of transmission, and to provide enough examples to show how various assumptions can be combined into simple, low-dimensional models.
We begin by revisiting existing models for an SIR disease spreading through mass-action mixing and for an SIR disease spreading through a sexual contact network. We explore a number of models of the spread of an SIR disease through a sexual contact network combined with another transmission mechanism. In all of the models we build, we assume that the sexually infectious period lasts longer than the period of infectiousness through the other mechanism. It is straightforward to modify this assumption. We divide the models presented into two broad classes: We consider a random static sexual network (a "configuration model" network (Newman, 2003)) combined with some other transmission model, in particular e mass action mixing, e vector-borne transmission, e or social contact network. We then consider a simple mass-action transmission model combined with a more complex sexual transmission network including: e dynamic partnership changes e preferential mixing.

Models with configuration model networks
Throughout we assume that S, I, and R (and any subdivisions of these classes) represent the proportion of the population in the susceptible, infected, or recovered state. We assume the outbreak is initialized with a fraction r of the population chosen uniformly at random and infected at time t ¼ 0.

The standard models
We briefly review the mass-action SIR model (Kermack & McKendrick, 1927) and the Edge-based compartmental model Miller, 2011). Although these models appear structurally different, a simple change of variables shows that they are closely related. This close relation will allow us to combine mixtures of the models.

The mass-action model
We begin with the mass-action SIR model. An infectious individual transmits infectious doses as a Poisson process with rate b. Each infectious dose is received by an individual randomly chosen from the population. If the recipient is susceptible, then she becomes infectious. Infected individuals recover as a one-step Poisson process with rate g. Once recovered they are immune to future infection. The diagram in the top-left of Fig. 1 leads us to the equations 2 shows that the solutions of system (1) accurately predicts the outcome of stochastic simulation in the largepopulation limit. The fit is excellent. Indeed, the equations correctly predict the large-population limit of the simulations. The convergence can be understood through the results of (Kurtz, 1971).

Edge-based compartmental model
For the simplest network-based model, we assume that a function PðkÞ is known which gives the probability a randomly chosen individual has k partners (its degree). If we assume that partners are randomly chosen, then the probability a random partner of a random individual has degree k is P n ðkÞ ¼ kPðkÞ=〈K〉 where 〈K〉 ¼ P kPðkÞ is the average degree.
This model differs from the mass-action model by assuming individuals transmit only through partnerships which do not change. Thus as an individual transmits to its partners, the number of available susceptible partners is reduced. The transmissions occur as Poisson processes with rate t. As before, if the recipient is susceptible he becomes infected, and eventually recovers to an immune state with rate g. The two diagrams on the right of Fig. 1 leads us to the Edge-based compartmental model (EBCM)  for a population with arbitrary initial fraction r infected Miller, 2011Miller, , 2014Kiss et al., 2017, chapter PðkÞx k is the probability generating function of the degree distribution. Here qðtÞ represents the probability that at time t a randomly chosen partner v of a randomly chosen individual u has not transmitted infection to u and f I represents the probability that v is infected and has not transmitted infection to u. These equations have been proven correct for random networks of given degree distribution in the infinite population limit (Janson, Luczak, & Windridge, 2014), and under stronger assumptions by (Decreusefond, Dhersin, Moyal, & Tran, 2012). See also (Ball & Neal, 2008) which derived a large system of equations that can be shown to be equivalent to system (2) (Miller & Kiss, 2014). We briefly outline a derivation of these equations, glossing over some more subtle details. We first note that the probability a degree k individual is susceptible is q k . Thus taking a weighted average of q k , we conclude that SðtÞ ¼ P k PðkÞqðtÞ k ¼ jðqÞ. That I ¼ 1 À S À R and _ R ¼ gI is straightforward. So it remains to find an equation for q. We divide q into three parts, q ¼ f S þ f I þ f R based on the status of v.
f S : the probability that the partner v of u is susceptible (and has not transmitted to u). f I : the probability that v is infected and has not transmitted to u. f R : the probability that v is recovered and did not transmit to u. We can find f S in terms of q starting with the observation v has degree k with probability P n ðkÞ ¼ kPðkÞ=〈K〉 which accounts for the fact that a random partner can be expected to have higher degree than a random individual (Feld, 1991). Following the derivation of S, we get f S ¼ ð1 À rÞ P kPðkÞq kÀ1 =〈K〉 where q kÀ1 is the probability that none of v's other partners have transmitted to it (we can assume u has not). Thus We can find f R in terms of q starting with the observation that _ f R ¼ gf I . Then _ f R and _ q are proportional. Using qð0Þ ¼ 1 and More details (including a subtle argument explaining that we can ignore transmissions from u to v) are found in  and (Kiss et al., 2017, chapter 6). The right hand side of Fig. 2 compares stochastic simulation with solutions of system (2). The fit is excellent.

Model similarity
The two models (1) and (2) appear quit dissimilar structurally, which would make it difficult to derive a hybrid model. However, we can modify the mass action system (1) making it easier to combined the models.
Ið b tÞ d b t to be the expected number of transmissions a random individual has received by time t.
Using an integrating factor we see From this it follows that individuals is shown in darker color, almost exactly matching the ODE solution. In the mass action model, the parameters are b ¼ 10 and g ¼ 1. The network has Pð2Þ ¼ Pð4Þ ¼ 0:5, so jðxÞ ¼ ðx 2 þ x 4 Þ=2. The parameters are t ¼ 2 and g ¼ 1.
From the definition of x we conclude _ x ¼ bI. Because _ R ¼ gI, and Rð0Þ ¼ xð0Þ ¼ 0, we can show that R ¼ gx=b. This leads us to S ¼ ð1 À rÞe Àx (3a) with the initial condition xð0Þ ¼ 0 : The variable x plays a similar role to q [more precisely, e Àx plays the same role as jðqÞ].
We provide a direct derivation of system (3), for which we consider both diagrams on the left of Fig. 1. We interpret x as the expected number of transmissions an individual has received. An individual is susceptible if she has received no transmissions and was not initially infected. If the expected number of transmissions is x, and each is randomly assigned, the probability of escaping transmission is e Àx . Thus S ¼ ð1 À rÞe Àx . To find x, we see that the rate of receiving transmissions is bI, so _ x ¼ bI. Then I and R follow as before, and we have a direct derivation of system (3) without going through system (1) first.
2.1.4. ℛ 0 An important quantity of mathematical epidemiology is the basic reproductive number ℛ 0 . This is often defined to be the expected number of infections caused when a single infected individual is introduced into a fully susceptible population. When the population exhibits heterogeneity, or partnerships have non-negligible duration, this definition breaks down. If we want ℛ 0 to accurately capture the transmission in the early stages of an epidemic, then its definition must account for the fact that the typical infected individual early in an epidemic may not be a typical individual in the population, and we must account for the fact that an individual cannot reinfect his infector. In defining ℛ 0 we must account for what the typical infected individuals looks like once the system settles down rather than what the typical introduced infected individual looks like (Diekmann, Heesterbeek, & Metz, 1990, 2012.
2.1.4.1. Mass-action model. In the mass-action model, we assume infected individuals recover with rate g, thus having an average infection duration of 1=g. All infected individuals transmit at rate b and the recipient is chosen uniformly at random from the population. Early in the epidemic, all transmissions reach a susceptible individual, thus causing infection. The expected number of infections caused by an individual infected early in the epidemic is thus Network-based model. In the network model, transmissions go only to direct partners of the infected individuals, and a non-negligible proportion of these transmissions may go to an individual's infector, thus having no effect. Early in the epidemic, an individual is chosen by the disease with probability proportional to their degree. So the early infections have degree k with probability P n ðkÞ ¼ kPðkÞ=〈K〉. The probability of transmitting along a partnership prior to recovery is t=ðt þ gÞ.
Thus the expected number of transmissions is There is a well-known final size relation for the mass-action model (1) linking the number of infections to ℛ 0 . The usual derivation involves dividing _ I by _ S, finding I as a function of S and then finding the values of S for which I ¼ 0. We find that system (3) makes this more direct. We simply note that at equilibrium _ x ¼ 0, and solve to get xð∞Þ=ℛ 0 ¼ 1 À ð1 À rÞe Àxð∞Þ Substituting R ¼ x=ℛ 0 , we conclude Rð∞Þ ¼ 1 À ð1 À rÞe Àℛ0Rð∞Þ which is the final size relation. The easiest way to find a solution is to iterate, starting with a guess Rð∞Þ ¼ 0.
There is a similar final size relation (Newman, 2002). The usual derivation is based on directly deriving a consistency relation without having a dynamic system of equations. We derive it from the EBCM equation (2). Taking f I ¼ 0, we find We also solve this through iteration, starting with qð0Þ ¼ 1.
Both the network and mass action final size relations can be intepreted as calculating the number of transmissions occurring given the number of infected individuals, and then in turn calculating the number of infected individuals assuming that those transmissions are randomly distributed. In general, this method can be used to calculate all final size relations of which we are aware.
It is instructional to recognize that with appropriate initial conditions the iteration corresponds directly to solving the dynamics of a particular discrete-time (Reed-Frost) disease model. 1 Each iteration gives the successive generation's size. Thus iterating until convergence gives the result as the number of generations tends to ∞. This is discussed in more detail in (Pellis, Ferguson, & Fraser, 2008) (see also (Ludwig, 1975) and some cautionary discussion in (Fennell, Melnik, & Gleeson, 2016)).
If, however, the timing of v's infection affects whether it would transmit to u (as may happen if there are seasonal effects or if the number of other infections in the same timestep somehow influences v's likelihood of transmitting), then we cannot derive a final size relation in the implicit way we have done it here. Instead, we must solve the dynamical equations. This represents the fact that we cannot predict the number of transmissions v has caused until we know the time at which v becomes infected.

The basic combined model
Our first new model assumes two modes of transmission: one is mass-action, while the other is through a sexual contact network. We take q to be the probability a random sexual partner of an individual u has not transmitted to u, while we define xðtÞ to be the expected number of transmissions u has received by time t through mass action transmission. We assume there are four possible statuses. Individuals begin susceptible (S). Through either sexual or mass-action interactions, they become infectious. The infectious phase occurs in two stages. The initial stage (I 1 ) is infectious through both mass action and sexual contact. The second stage (I 2 ) is infectious only through sexual contact. They finally reach an immune recovered state (R). The symbols S, I 1 , I 2 , and R denote both the stage and the fraction of individuals in that stage, so S þ I 1 þ I 2 þ R ¼ 1.
The per-individual transition rate from I 1 to I 2 is g 1 , and the per-individual transition rate from I 2 to R is g 2 . The mass action transmission rate is b, while the per-partnership transmission rate during I 1 is t 1 and during I 2 is t 2 . Fig. 3 shows the corresponding flow diagrams for the combined model. We can find SðtÞ ¼ ð1 À rÞe ÀxðtÞ jðqðtÞÞ, representing the fact that S is the probability an individual was not initially infected, has not been infected through mass-action transmission, and has not been infected through sexual transmission. Using the top diagram in Fig. 3 then shows that _ I 2 ¼ g 1 I 1 À g 2 I 2 and _ R ¼ g 2 I 2 . Conservation of probability gives I 1 ¼ 1 À S À I 2 À R. We are only missing equations for q and x.
As before the equation for x is simply _ x ¼ bI 1 . The equation for q is very similar to before. We have _ q ¼ t 1 f I;1 À t 2 f I;2 . We will use f I;1 ¼ q À f S À f I;2 À f R . We can again find f S explicitly, the only change is that there is an additional factor e Àx accounting for mass-action transmissions: f S ¼ ð1 À rÞe Àx j 0 ðqÞ=〈K〉. We cannot explicitly find f I;2 and f R , but we have differential equations for them: _ f I;2 ¼ g 1 f I;1 À ðt 2 þ g 2 Þf I;2 and _ f R ¼ g 2 f I;2 . The final system is 1 For example, with (4), if we take q 0 ¼ 1 and f S ð0Þ ¼ 1 À r, then iteration directly gives us q 1 , the value after one generation from which we can infer the number infected in the first generation. Repeating, we can capture the entire dynamics.
If we take t 2 ¼ 0 or g 2 ¼ ∞, then assuming a randomly introduced infection, this model is equivalent to the model considered by (Ball & Neal, 2008), although it is much more compact. The equivalence follows from techniques of (Miller & Kiss, 2014). Fig. 4 compares simulated epidemics in large populations with predictions from system (5).
2.2.1. ℛ 0 In deriving ℛ 0 , we separate the individuals by their "generation", that is, the length of the transmission chain from the initial case. The relation between one generation and the next is found through a matrix multiplication. The result quickly settles down to a multiple of the dominant eigenvector, with ℛ 0 being the dominant eigenvalue. We demonstrate this calculation using two different approaches, because each will be useful in the following models. The approaches track different variables, resulting in different matrices.
Before giving the derivations, we note that the probability a sexual partnership will transmit during the first infectious stage if one individual is infected is t 1 =ðt 1 þ g 1 Þ, and if transmission does not occur in the first stage the probability of transmission in the second stage is t 2 =ðt 2 þ g 2 Þ. Combining this, the full probability of sexual transmission within a partnership is The expected number of mass-action transmissions from a single infected individual is R ma ¼ b=g 1 . The expected number of sexual transmissions from an individual infected through the mass-action route is R sejma ¼ T se 〈K〉 and the expected number of sexual transmissions from an individual infected through the sexual route is R sejse ¼ T se 〈K 2 À K〉=〈K〉.
2.2.1.1. First derivation. For our first derivation, we distinguish individuals based on whether they were infected through a mass-action transmission or through a sexual transmission. The distinction is needed because early in the epidemic individuals infected through a mass-action transmission have on average 〈K〉 susceptible sexual partners, while those infected through a sexual transmission have on average 〈K 2 À K〉=〈K〉 susceptible sexual partners. The number of mass-action transmissions caused by each class of individual is the same. We define N ma ðgÞ and N se ðgÞ to be the number of individuals in generation g which were infected through mass-action and sexual transmissions respectively. Then The dominant eigenvalue is 2.2.1.2. Second derivation. An alternate pair of variables that could be considered is the number of infected individuals and the number of discordant (SI) partnerships. We use NðgÞ to denote the number infected in generation g and yðgÞ to be the number of discordant partnerships in generation g. Then Nðg þ 1Þ is the sum of the number of infections caused by the mass action route, R ma NðgÞ and the number caused by sexual transmission T se yðgÞ ¼ R sejma yðgÞ=〈K〉. Similarly the number of new discordant partnerships is the sum of the number created through the mass-action route 〈K〉R ma NðgÞ and the sexual route ½T se 〈K 2 À K〉=〈K〉yðgÞ ¼ R sejse yðgÞ. Thus The eigenvalues of this matrix are the same as before.
2.2.1.3. Approximations of ℛ 0 . We highlight the fact that ℛ 0 is not simply a sum of R ma and R sejse . We will discuss this more in the vector-borne disease model. It is useful to recognize that by expanding the square root as a binomial series, our expression for ℛ 0 can be approximated as If either R ma or R sejse is very small, then this approximates as R sejse or R ma respectively If R ma R sejma s0, then ℛ 0 is larger than the maximum of R ma and R se .
However, if R ma zR sejse , this approximation breaks down. Then writing R ma ¼ R sejse þ ε, we find

Final size relation
For this model, we can derive a consistency relation which leads to the final size. We can do this in two ways. We first show an indirect derivation by integrating the full dynamical equations. However, these equations can be derived directly through a mechanistic argument, and so we show this argument as well.
We first find xð∞Þ. We have Rð∞Þ ¼ g 2 Z ∞ 0 I 2 dt and xð∞Þ ¼ b It is straightforward to see that I 2 ð0Þ ¼ 0 and I 2 ð∞Þ ¼ 0, so integrating the equation for _ We attack the integrals sequentially. We have _ f R þ _ f I;2 ¼ g 1 f I;1 À t 2 f I;2 . Integrating from 0 to ∞ and noting f I;2 ð0Þ ¼ f I;2 ð∞Þ ¼ 0, we conclude Rearranging this final expression, we conclude We can solve this system iteratively and substitute the result into our expression for Sð∞Þ to find the final size of the epidemic.

Direct derivation.
It turns out we can derive this directly: We interpret qð∞Þ as the probability that a sexual partner either would never transmit even if infected (1 À T se ), or would transmit, but is never infected (T se ð1ÀrÞe Àxð∞Þ j 0 ðqð∞ÞÞ 〈K〉 ). Adding these yields the expression above. We can similarly derive xð∞Þ. This is interpreted as the expected number of mass action transmissions an individual receives. The expected number caused by each infected individual is b=g 1 . So xð∞Þ ¼ bð1 À Sð∞ÞÞ=g 1 . Substituting for S, we have h 1 À ð1 À rÞe Àxð∞Þ jðqð∞ÞÞ i and thus we have a direct derivation of the consistency relation for the final size.

Vector-borne transmission
We now consider transmission with a vector component, which we will refer to as a "mosquito". We assume there is no latent phase: when an individual or a mosquito is infected it is immediately infectious. We again assume two infectious phases in humans, with human to mosquito transmission possible only in the first phase. We make the same assumptions as before about sexual transmission.
We assume mosquito lifetimes are shorter than epidemic duration. We assume the mosquito population is always close to equilibrium, but more realistic dynamics could be added. Once infected, a mosquito remains infectious until death. The average lifespan of a mosquito is taken to be 1=d, corresponding to a death-rate of d. The influx of new mosquitos is at a constant rate B, measured in mosquitos per unit time per person. The human to mosquito transmission rate is b 1 , and the mosquito to human transmission rate is b 2 . We use V S and V I to represent the number of susceptible and infected mosquitos per person. We have Guided by the mixed network-mass-action model, we set x to be the expected number of transmissions received by a human from mosquitos. Taking other variables as before, and following the flow diagram in Fig. 5, we have S ¼ ð1 À rÞe Àx jðqÞ (6a) (6l) Fig. 6 compares simulated epidemics in large populations with solutions to the equations. We again see excellent fit.
2.3.1. ℛ 0 Because of the similarity to the mass action model we see that R sejse , the expected number of sexual transmissions from an individual infected through sexual transmission early in the epidemic, is as before. We use R sejmo to be the expected number of sexual transmissions caused by an individual infected by mosquitos early in the epidemic. We define R mo to be the number of human infections resulting from a single infected human through mosquito-based transmission. We assume that the mosquito population is at equilibrium, so R mo ¼ ðB=dÞðb 1 =g 1 Þðb 2 =dÞ. Then ℛ 0 is the leading eigenvalue of R mo R mo R sejmo R sejse and so it takes the same form as before. Note that in this definition, we have considered human to mosquito to human to constitute a single generation. An alternate formulation might consider human to mosquito and mosquito to human each as a separate generation. This will lead to a different expression for ℛ 0 . See, for example, (Gao et al., 2016). However, the resulting threshold ℛ 0 ¼ 1 would still be the same.
The fact that ℛ 0 is not simply the sum R mo þ R sejse has important implications. Consider a population in which R sejse > 1, but R sejmo < 1, and assume that the vast majority of transmissions are through mosquitos. If we measure the expected number of sexual transmissions an infected individual causes, this will be close to R sejmo . If we conflate this measured value with R sejse , then we will underestimate the possibility of a purely sexual epidemic in a population without mosquitos (Allard, Althouse, H ebert-Dufresne, & Scarpino, 2016). Fig. 7. Flow diagrams leading to system (7) for a population in which interactions all lie within a sexual or a social network. We assume that the infectious stage can be subdivided into two stages, the first of which transmits through either contact, and the second transmits only through sexual contact.

Final size relation
We are not able to find a simple recurrence leading to a finite size relation, and do not believe such a recurrence exists. Generally final size relations can be derived by calculating the number of infections given the number of transmissions occurring, and then relating the number of transmissions occurring given the number of infections (Miller, 2012). In this model such a relation is not possible because the timing of infections matter.
If we assume that all infections happen across a brief window in which there is negligible mosquito turnover, the number of infected mosquitos is bounded by the equilibrium mosquito population size: many mosquitos might bite multiple infected individuals, and many bites have no effect. However, if the epidemic is much slower, then this blocking mechanism is reduced, and so the total number of mosquitos infected may be larger, thus causing a different number of total human infections. It would be possible to determine the number of humans infected if the total number of infected mosquitos is known.
By appropriately changing the parameters of the mosquitos, we could keep R mo unchanged, but still alter the total number of infected mosquitos, thus altering the total number of infected humans. The usual techniques to derive a final size relation will fail.

Multi-layer networks
We now consider a population in which the potentially transmitting contacts occur within one of two networks. The individuals in the networks are the same, but the partnerships have different interpretations. The first network is the network of sexual partnerships, while the second is a network of social interactions. As before we assume two infectious stages. The first infectious stage transmits with rate t 1 to sexual partners and b to social partners. The second stage transmits with rate t 2 to sexual partners and not at all to social partners.
We assume a joint degree distribution Pðk se ; k so Þ, allowing the number of sexual partners to be correlated with the number of social partners, and we define jðx; yÞ ¼ P kse;kso Pðk se ; k so Þx kse y kso . We define q se to be the probability a sexual partner has not transmitted to an individual u and q so to be the probability a social partner has not transmitted to u (ignoring the possibility of u transmitting to its partners). As before we use f se S , f se I;1 , f se I;2 , f se R and f so S , f so I;1 , and f so R with the superscript denoting the type of partnership and the subscript denoting the status of the partner. As the second infected state does not cause social infections, we lump it in with R for social infections, and write f so I;1 as simply f so I . The diagrams in Fig. 7 give S ¼ ð1 À rÞjðq se ; q so Þ (7a) (7e) Fig. 8. A comparison of simulation with solutions to system (7) for transmission through two overlapping static networks. The cloud consists of simulations in a population of 10 3 individuals, with three simulations highlighted. The other solid colored curves represent a single simulation in a population of 10 5 individuals. The black dashed curve represents the solution to the ODE of system (7). In the sexual network, a quarter of the population has 2 partners, half have 3, and a quarter have 4. In the social network, half have 10 partners and half have 20. Their values are correlated so that the least active individuals in both networks are the same. That is, jðx; yÞ ¼ ðx 2 y 10 þ x 3 y 10 þ x 3 y 20 þ x 4 y 20 Þ=4. The parameters are b ¼ 0:1, t 1 ¼ 0:5, g 1 ¼ 1, t 2 ¼ 0:2, and g 2 ¼ 2.
f se Fig. 8 compares simulated epidemics in large populations, with two networks along which disease transmits. Again the fit is excellent.

ℛ 0
To find ℛ 0 we define R sojso to be the expected number of social transmissions caused given that an individual is infected through a social transmission, R sojse to be the expected number of social transmissions caused given that an individual is infected through a sexual transmission, and similarly R sejse and R sejso . For notational simplicity, we define T so ¼ b=ðb þ g 1 Þ to be the probability a social contact will transmit before the infected individual moves into the second infectious stage. We again define T se ¼ 1 À g 1 t1þg 1 g 2 t2þg 2 to be the probability a sexual contact will transmit in (at least) one of the two stages. We have Fig. 9. Flow diagrams leading to system (8) for mass action transmission combined with transmission through a dynamic sexual network. The top diagram is as before. The middle diagram is much as before for the network component, but because partnerships can change there are additional paths to take. Because of this we cannot solve for f S explicitly and we must calculate all the fluxes in and out of f S . The fluxes due to finding new partners depend on the probability that the new partner has a given status. These are tracked using the diagram in the bottom left. The bottom right is as before.
Note that j xx ð1; 1Þ ¼ 〈K 2 se À K se 〉, j xy ð1; 1Þ ¼ 〈K so K se 〉, and j yy ð1; 1Þ ¼ 〈K 2 so À K so 〉 Taking N se ðgÞ and N so ðgÞ to be the number infected by each method in generation g we have The dominant eigenvalue of this matrix is ℛ 0 .

Models with more complex network structure
We have considered disease spread with a fairly simple network structure combined with various other transmission mechanisms. Now we take the other transmission mechanism to be mass action mixing and consider more complex network structure. Fig. 10. A comparison of simulation with solutions to system (8) for mass action transmission combined with sexual transmission through a dynamic network. The cloud consists of simulations in a population of 10 3 individuals, with three simulations highlighted. The additional solid colored curve is a single simulation in a population of 10 5 individuals. The dashed black curve is the solution to system (8). The parameters are b ¼ 5, t 1 ¼ 1, g 1 ¼ 5, t 2 ¼ 0:3, and g 2 ¼ 1. We take jðxÞ ¼ ðx 2 þ 2x 3 þ x 4 Þ=4.

Dynamic network
We begin by assuming that partnerships change in time. We assume the partnership duration is exponentially distributed with mean 1=h. When a partnership ends, the two individuals immediately form new partnerships. Conceptually, we can think of an individual with k partners as having k stubs (binding sites in the language of (Leung, 2016;Leung, Kretzschmar, & Diekmann, 2015)) which connect to stubs of other individuals to form temporary partnerships. Partnerships terminate with rate h.
We use a refined definition of q: It is the probability that for a given stub of u no partner connected to that stub has ever transmitted to u. As before the probability u is susceptible is SðtÞ ¼ ð1 À rÞe ÀxðtÞ jðqðtÞÞ, but determining f S , f I , and f R is different. For this we follow (Kiss et al., 2017, chapter 8).
Let v be a random partner of u at time t. We define zðtÞ to be the probability that the stubs currently forming their partnership have not brought infection to either u or v prior to joining. Then The derivative of z is hq 2 À hz: the rate at which such partnerships are created minus the rate at which such partnerships are destroyed by partnership turnover.
We have an option of how to build the model. We could write f S ¼ ð1 À rÞe Àx z j 0 ðqÞ 〈K〉 and use a differential equation for z. However, we choose to use a differential equation for f S and eliminate z. We write We require new variables p S , p I;1 , p I;2 and p R , representing the probability a random stub belongs to an individual of each status. We can directly show that p S ¼ ð1 À rÞ P k kPðkÞq k ¼ ð1 À rÞqj 0 ðqÞ. Fig. 9 leads to a new system of equations Fig. 11. Flow diagrams leading to system (9) for mass action transmission combined with sexual transmission in a preferentially mixing sexual network. The top diagram is as before. The middle diagram tracks q k so there is a corresponding equation for each degree k. The bottom diagram is as before.
S ¼ ð1 À rÞe Àx jðqÞ (8a) Our initial conditions are 10 compares simulated epidemics in a large population. Transmissions occur through a mass action mechanism and through transmission in a dynamic sexual network. Fig. 12. A comparison of simulation with solutions to system (9) for mass action transmission combined with a sexual network exhibiting preferential partnership formation. The cloud represents 200 simulations in a population of 10 3 individuals, with three simulations highlighted. The additional solid colored curve comes from a simulation of 10 5 individuals. The dashed black curve is the solution to system (9). The parameters are given in the text.
3.1.1. ℛ 0 Although ℛ 0 can be calculated for this model, the calculation is quite tedious. We must account for the fact that when an SI partnership forms, it can transmit, the partnership can dissolve and be replaced, and retransmit again. It become particularly tedious as we must consider the possible state of the partnership when the first transition happens and calculate different outcomes for whether it is II or SI, and take an appropriate weighted average.
For now, we note that the ℛ 0 calculated for the static network case is a lower bound. We can find an upper bound by assuming the h/∞ limit such the partnerships always change prior to the next transmission. Then when an individual becomes infected along a sexual contact we expect 〈K 2 〉=〈K〉 partnerships available to transmit at all times. However when an individual becomes infected through a mass action transmission we expect 〈K〉 partnerships to be available. The number of transmissions occurring for each is t1 Following the second approach of section 2.2.1 we set NðgÞ to be the number of individuals infected in generation g and yðgÞ to be the number of partnerships those infected individuals are in. We find NðgÞ þ t 1 g 1 þ t 2 g 2 yðgÞ and the number of partnerships they are in is Thus we have defined a matrix problem, and the dominant eigenvalue of is an upper bound for ℛ 0 .

Final size relation
Again we are unable to find a final size relation for much the same reason as in the vector-borne case. The relative timing of the epidemic and the partnership turnover mean that the dynamics of the process play an important role. We cannot express how many transmissions some number of infected individuals will cause because.

Sexual network with preferential mixing
Finally we consider a model in which individuals select their partners preferentially according to degree. This has been studied using pair-approximations (Moreno, G omez, & Pacheco, 2003), but to our knowledge, the corresponding EBCM model has not been explicitly considered. It can be inferred from model 2.2.3 of (Miller & Volz, 2013), and appears as an exercise in (Kiss et al., 2017, chapter 6). We will study the slightly more general case in which there are two infectious stages as before, and allow for the first infectious stage to also exhibit mass action transmission.

ℛ 0
To find ℛ 0 for mass action transmission combined with sexual transmission on a static preferentially mixing network, we define N k ðgÞ to be the number of individuals having k partners that are infected in generation g. We again use T se ¼ 1 À g 1 g 2 ðt1þg 1 Þðt2þg 2 Þ to be the probability that an infected individual transmits to a susceptible partner before recovering. The expected number of infections of individuals with degree b k caused through sexual transmission by an infected individual with degree k is T se kPð b k kÞ. The expected number of infections of individuals with degree b k caused through mass action transmission by any infected individual is Pð b kÞb=g 1 . Thus we conclude that Then ℛ 0 is the dominant eigenvalue of the matrix given by the sum within the square brackets.

Final size
Unlike the dynamic case, we can derive a final size relation for the preferential mixing model. This is very similar to our basic combined mass action and network model.