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

The emergence of Zika and Ebola demonstrates the importance of understanding the role of sexual transmission in the spread of diseases with a primarily non-sexual transmission route. In this paper, we develop low-dimensional models for how an SIR disease will spread if it transmits through a sexual contact network and some other transmission mechanism, such as direct contact or vectors. We show that the models derived accurately predict the dynamics of simulations in the large population limit, and investigate ℛ0 and final size relations.

Introduction

The recent emergence of Ebola (Christie et al., 2015, 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 et al., 1983, Emond et al., 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 (Christie et al., 2015, 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 and 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 et al., 2017, Miller and Kiss, 2014). However, in the specific case of Susceptible–Infected–Recovered (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 et al., 2012, Miller and 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

mass action mixing,

vector-borne transmission,

or social contact network.

We then consider a simple mass-action transmission model combined with a more complex sexual transmission network including:

dynamic partnership changes

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 ρ of the population chosen uniformly at random and infected at time .

The standard models

We briefly review the mass-action SIR model (Kermack & McKendrick, 1927) and the Edge-based compartmental model (Miller et al., 2012, 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 β. 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 γ. Once recovered they are immune to future infection. The diagram in the top-left of Fig. 1 leads us to the equationswith initial conditions

Fig. 1

Flow diagrams leading to systems (1), (2), and (3). (top left) A flow diagram demonstrating the standard model of mass action SIR dynamics. (right) Two flow diagrams showing the EBCM model of network SIR dynamics. (bottom left) A flow diagram, that when combined with the diagram above it allows for a system of equations similar to the EBCM model.

Fig. 2 shows that the solutions of system (1) accurately predicts the outcome of stochastic simulation in the large-population 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).

Fig. 2

(left) A comparison of stochastic agent-based simulation with ODE solution for the mass action model (1) or (3). (right) A comparison of stochastic agent-based simulation with ODE solution for the network-based model (2). In both cases, dashed curve denotes the solution to the ODE. A cloud made up of 200 simulations in populations of individuals is shown lightly in color, with 3 of these simulations highlighted. An additional simulation in a population of individuals is shown in darker color, almost exactly matching the ODE solution. In the mass action model, the parameters are and . The network has , so . The parameters are and .

Edge-based compartmental model

For the simplest network-based model, we assume that a function 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 where 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 τ. As before, if the recipient is susceptible he becomes infected, and eventually recovers to an immune state with rate γ. The two diagrams on the right of Fig. 1 leads us to the Edge-based compartmental model (EBCM) (Miller et al., 2012) for a population with arbitrary initial fraction ρ infected (Miller et al., 2012, Miller, 2011, Miller, 2014, Kiss et al., 2017, chapter 6)with initial conditionswhere is the probability generating function of the degree distribution. Here 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 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 . Thus taking a weighted average of , we conclude that . That and is straightforward. So it remains to find an equation for θ. We divide θ into three parts, based on the status of v.

: the probability that the partner v of u is susceptible (and has not transmitted to u).

: the probability that v is infected and has not transmitted to u.

: the probability that v is recovered and did not transmit to u.

We can find in terms of θ starting with the observation v has degree k with probability 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 where is the probability that none of v's other partners have transmitted to it (we can assume u has not). Thus

We can find in terms of θ starting with the observation that . Then and are proportional. Using and , we can show thatThusMore details (including a subtle argument explaining that we can ignore transmissions from u to v) are found in (Miller et al., 2012) 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.

We define to be the expected number of transmissions a random individual has received by time t. Using an integrating factor we seebecomes

From this it follows thatFrom the definition of ξ we conclude . Because , and , we can show that . This leads us towith the initial conditionThe variable ξ plays a similar role to θ [more precisely, plays the same role as ].

We provide a direct derivation of system (3), for which we consider both diagrams on the left of Fig. 1. We interpret ξ 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 ξ, and each is randomly assigned, the probability of escaping transmission is . Thus . To find ξ, we see that the rate of receiving transmissions is , so . Then I and R follow as before, and we have a direct derivation of system (3) without going through system (1) first.

An important quantity of mathematical epidemiology is the basic reproductive number . 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 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 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 et al., 1990, Diekmann et al., 2009, Diekmann et al., 2012).

Mass-action model

In the mass-action model, we assume infected individuals recover with rate γ, thus having an average infection duration of . All infected individuals transmit at rate β 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 . The probability of transmitting along a partnership prior to recovery is . Thus the expected number of transmissions is

Final size relation

There is a well-known final size relation for the mass-action model (1) linking the number of infections to . The usual derivation involves dividing by , finding I as a function of S and then finding the values of S for which . We find that system (3) makes this more direct. We simply note that at equilibrium , and solve to get Substituting , we concludewhich is the final size relation. The easiest way to find a solution is to iterate, starting with a guess .

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 (2a), (2b), (2c), (2d), (2e). Taking , we findWe also solve this through iteration, starting with .

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 θ to be the probability a random sexual partner of an individual u has not transmitted to u, while we define 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 () is infectious through both mass action and sexual contact. The second stage () is infectious only through sexual contact. They finally reach an immune recovered state (R). The symbols S, , , and R denote both the stage and the fraction of individuals in that stage, so .

The per-individual transition rate from to is , and the per-individual transition rate from to R is . The mass action transmission rate is β, while the per-partnership transmission rate during is and during is . Fig. 3 shows the corresponding flow diagrams for the combined model.

Fig. 3

Flow diagrams leading to system (5). The top diagram shows transitions of individuals between the states. The middle diagram shows transitions of the status of a partner v of a randomly chosen individual u as well as whether v has transmitted to u (all ignoring transmissions to u to v). The bottom diagram shows the change in ξ and how it relates to . These diagrams are related together by the dependence of the flux into ξ on and on the fact that we can express S and in terms of ξ and θ.

We can find , 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 and . Conservation of probability gives . We are only missing equations for θ and ξ.

As before the equation for ξ is simply . The equation for θ is very similar to before. We have . We will use . We can again find explicitly, the only change is that there is an additional factor accounting for mass-action transmissions: . We cannot explicitly find and , but we have differential equations for them: and . The final system is

If we take or , 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).

Fig. 4

A comparison of simulation with solutions to system (5) for network-based transmission combined with mass-action transmission. The colored cloud consists of 200 simulations in a population of individuals, with three highlighted. A darker solid colored curve shows a simulation in a population of . This is almost exactly overlain by the solution to the ODE system (5) shown in a black dashed curve. The parameters are , , , , and . The function ψ takes the form .

In deriving , 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 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 , and if transmission does not occur in the first stage the probability of transmission in the second stage is . 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 . The expected number of sexual transmissions from an individual infected through the mass-action route is and the expected number of sexual transmissions from an individual infected through the sexual route is .

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 susceptible sexual partners, while those infected through a sexual transmission have on average susceptible sexual partners. The number of mass-action transmissions caused by each class of individual is the same.

We define and to be the number of individuals in generation g which were infected through mass-action and sexual transmissions respectively. Then

The dominant eigenvalue is

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 to denote the number infected in generation g and to be the number of discordant partnerships in generation g. Then is the sum of the number of infections caused by the mass action route, and the number caused by sexual transmission . Similarly the number of new discordant partnerships is the sum of the number created through the mass-action route and the sexual route . ThusThe eigenvalues of this matrix are the same as before.

Approximations of

We highlight the fact that is not simply a sum of and . 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 can be approximated as

If either or is very small, then this approximates as or respectively If , then is larger than the maximum of and .

However, if , this approximation breaks down. Then writing , we find

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.

Indirect derivation

We can derive a final size relation. We take and note that all the variables corresponding to active infectious states must go to zero. ThenWe use this to express and in terms of and . This yields an implicit relation to solve.

We first find . We have and . It is straightforward to see that and , so integrating the equation for yields . We conclude

We similarly find . We have , and is given above. Thus . Integrating the equation for yields

We attack the integrals sequentially. We have . Integrating from 0 to and noting , we concludeSo we haveSimilarly since and , we have . ThusRearranging this final expression, we concludeWe can solve this system iteratively and substitute the result into our expression for to find the final size of the epidemic.

Direct derivation

It turns out we can derive this directly: We interpret as the probability that a sexual partner either would never transmit even if infected (), or would transmit, but is never infected (). Adding these yields the expression above. We can similarly derive . This is interpreted as the expected number of mass action transmissions an individual receives. The expected number caused by each infected individual is . So . Substituting for S, we haveand 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 , corresponding to a death-rate of δ. 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 , and the mosquito to human transmission rate is . We use and to represent the number of susceptible and infected mosquitos per person. We haveGuided by the mixed network-mass-action model, we set ξ 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

Fig. 5

Flow diagrams leading to system (6) for the combined sexual network and vector transmission routes. These are very similar to Fig. 3, except for the introduction of the vector component at the bottom.

Fig. 6 compares simulated epidemics in large populations with solutions to the equations. We again see excellent fit.

Fig. 6

A comparison of simulation with solutions to system (6). Because resolving individual mosquitos significantly increases calculation time, we run fewer calculations with smaller populations. The cloud of simulations consist of 100 simulations in a population of individuals with typically 20 mosquitos per individual. Again 3 are highlighted. An additional solid colored curve is given for a population of individuals. The dashed black curve represents the solution to the ODEs of system (6). The parameters are , , , , , , and , again with .

Because of the similarity to the mass action model we see that , the expected number of sexual transmissions from an individual infected through sexual transmission early in the epidemic, is as before. We use to be the expected number of sexual transmissions caused by an individual infected by mosquitos early in the epidemic. We define 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 . Then is the leading eigenvalue ofand 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 . See, for example, (Gao et al., 2016). However, the resulting threshold would still be the same.

The fact that is not simply the sum has important implications. Consider a population in which , but , 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 . If we conflate this measured value with , then we will underestimate the possibility of a purely sexual epidemic in a population without mosquitos (Allard, Althouse, Hébert-Dufresne, & Scarpino, 2016).

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 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 to sexual partners and β to social partners. The second stage transmits with rate to sexual partners and not at all to social partners.

We assume a joint degree distribution , allowing the number of sexual partners to be correlated with the number of social partners, and we define . We define to be the probability a sexual partner has not transmitted to an individual u and 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 , , , and , , and 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 as simply .

The diagrams in Fig. 7 give

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.

Fig. 8 compares simulated epidemics in large populations, with two networks along which disease transmits. Again the fit is excellent.

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 individuals, with three simulations highlighted. The other solid colored curves represent a single simulation in a population of 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, . The parameters are , , , , and .

To find we define to be the expected number of social transmissions caused given that an individual is infected through a social transmission, to be the expected number of social transmissions caused given that an individual is infected through a sexual transmission, and similarly and . For notational simplicity, we define to be the probability a social contact will transmit before the infected individual moves into the second infectious stage. We again define to be the probability a sexual contact will transmit in (at least) one of the two stages. We haveNote that , , and Taking and to be the number infected by each method in generation g we haveThe dominant eigenvalue of this matrix is .

The final size relation of this model is straightforward to derive directly. Given and , the probability a sexual partner of u is never infected is . Similarly the probability a social partner of u is never infected is . Translated into equations, these becomeWe can solve this iteratively. From the solution, we have .

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.

Dynamic network

We begin by assuming that partnerships change in time. We assume the partnership duration is exponentially distributed with mean . 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 et al., 2015)) which connect to stubs of other individuals to form temporary partnerships. Partnerships terminate with rate η.

We use a refined definition of θ: 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 , but determining , , and 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 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 ζ is : 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 and use a differential equation for ζ. However, we choose to use a differential equation for and eliminate ζ. We write

We require new variables , , and , representing the probability a random stub belongs to an individual of each status. We can directly show that . Fig. 9 leads to a new system of equationsOur initial conditions are

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 explicitly and we must calculate all the fluxes in and out of . 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.

Fig. 10 compares simulated epidemics in a large population. Transmissions occur through a mass action mechanism and through transmission in a dynamic sexual network.

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 individuals, with three simulations highlighted. The additional solid colored curve is a single simulation in a population of individuals. The dashed black curve is the solution to system (8). The parameters are , , , , and . We take .

Although 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 calculated for the static network case is a lower bound. We can find an upper bound by assuming the limit such the partnerships always change prior to the next transmission. Then when an individual becomes infected along a sexual contact we expect partnerships available to transmit at all times. However when an individual becomes infected through a mass action transmission we expect partnerships to be available. The number of transmissions occurring for each is .

Following the second approach of section 2.2.1 we set to be the number of individuals infected in generation g and to be the number of partnerships those infected individuals are in. We findand the number of partnerships they are in isThus we have defined a matrix problem, and the dominant eigenvalue ofis an upper bound for .

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ómez, & 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.

The flow diagrams of Fig. 11 lead toThe initial conditions are

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 so there is a corresponding equation for each degree k. The bottom diagram is as before.

Fig. 12 compares simulated epidemics in a large population, with a mass action transmission mechanism and a static sexual network exhibiting preferential mixing. The sexual network has , with the partnerships distributed so that

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 individuals, with three simulations highlighted. The additional solid colored curve comes from a simulation of individuals. The dashed black curve is the solution to system (9). The parameters are given in the text.

As before, the infectious period has two stages. The first infectious stage transmits through both mechanisms with , , and . The second stage has and .

To find for mass action transmission combined with sexual transmission on a static preferentially mixing network, we define to be the number of individuals having k partners that are infected in generation g. We again use to be the probability that an infected individual transmits to a susceptible partner before recovering. The expected number of infections of individuals with degree caused through sexual transmission by an infected individual with degree k is . The expected number of infections of individuals with degree caused through mass action transmission by any infected individual is . Thus we conclude thatWriting this as a matrix equation, we haveThen 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.

We note thatand

This provides the relation to solve. Once we have this, then