Dangerous connections: on binding site models of infectious disease dynamics

We formulate models for the spread of infection on networks that are amenable to analysis in the large population limit. We distinguish three different levels: (1) binding sites, (2) individuals, and (3) the population. In the tradition of physiologically structured population models, the formulation starts on the individual level. Influences from the ‘outside world’ on an individual are captured by environmental variables. These environmental variables are population level quantities. A key characteristic of the network models is that individuals can be decomposed into a number of conditionally independent components: each individual has a fixed number of ‘binding sites’ for partners. The Markov chain dynamics of binding sites are described by only a few equations. In particular, individual-level probabilities are obtained from binding-site-level probabilities by combinatorics while population-level quantities are obtained by averaging over individuals in the population. Thus we are able to characterize population-level epidemiological quantities, such as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_0$$\end{document}R0, r, the final size, and the endemic equilibrium, in terms of the corresponding variables.


Introduction
Consider an empirical network consisting of individuals that form partnerships with other individuals. Suppose an infectious disease can be transmitted from an infectious individual to any of its susceptible partners and thus spread over the network. Consider an individual in the network at a particular point in time. We are interested in the disease status of the individual, but also in the presence of the infection in its immediate surroundings that are formed by the individual's partners. We may label this individual by listing -its disease status in terms of the S, I, R classification, where, as usual, S stands for susceptible, I for infectious and R for recovered (implying immunity) -how many partners this individual has -the disease status of these partners In this spirit, we may provide a statistical description of the network at a particular point in time by listing, for each such label, the fraction of the population carrying it.
Is it possible to predict the future spread of the disease on the basis of this statistical description? The answer is 'no', simply because the precise network structure is important for transmission and we cannot recover the structure from the description. But if we are willing to make assumptions about the structure (and to consider the limit of the number of individuals going to infinity), the answer might be 'yes'. And even if the true answer is still 'no', we may indulge in wishful thinking and answer 'to good approximation'.
When considering an outbreak of a rapidly spreading disease, we can consider the network as static. If we are willing to assume that the network is constructed by the configuration procedure (Durrett 2006;van der Hofstad 2015), the answer is indeed 'yes' (Decreusefond et al. 2012; Barbour and Reinert 2013;Janson et al. 2014). But if the disease spreads at the time scale of formation and dissolution of partnerships, we need to take these partnership dynamics into account and next indeed rely on wishful thinking (though the answer may very well be 'yes'). In case of HIV, the disease spreads on the time scale of demographic turnover and this motivated our earlier work (Leung et al. 2012(Leung et al. , 2015 that also takes birth and death into account [(here we know that the answer is 'no', see Leung et al. (2015, Appendix B)].
In the rest of this introduction we first discuss the model formulation used and the relation between our work and existing literature. Next, we consider three different settings based on the time scales of disease spread, partnership dynamics, and demographic turnover. Individuals are decomposed into conditionally independent components (the 'binding sites') and we discuss how the dynamics of these binding sites can be specified. We end the introduction with an outline of the structure of the rest of the paper.

Physiologically structured population models
As in our earlier paper (Leung et al. 2015), our model formulation is in the tradition of physiologically structured population models (PSPM Metz and Diekmann 1986;Diekmann et al. 1998bDiekmann et al. , 2001. This means that we start from the notion of state at u v w u v w Fig. 1 An illustration of binding sites with three individuals u, v, and w. In this example, u, v, and w have four, three, and two binding sites for partners, respectively. On the left, all binding sites are free. On the right, a partnership between u and w is formed and they both have one occupied binding site the individual level, called i-state (where i stands for individual). Model specification involves, first of all, a description of changes in time of the i-state as influenced by i-state itself and the relevant environmental variables that capture the influence of the outside world. Next the model specifies the impact of individuals on the environmental variables. Thus the feedback loop that creates density dependence, i.e. dependence among individuals, is described in a two step procedure. To lift the i-level model to the population level (p-level) is just a matter of bookkeeping, see Diekmann and Metz (2010) for a recent account.
In the setting considered here, i-state ranges over a finite set. As a consequence, the p-level equations are ordinary differential equations (ODE). These ODE describe, apart from death and birth of individuals, the dynamical changes of i-state, i.e. how individuals jump back and forth between the various states. In the spirit of the theory of Markov chains (Taylor and Karlin 1998), we describe an individual not by its actual state but by the probability distribution, i.e. the probability of being in the various states. Equating a p-level fraction to an i-level probability provides the link between the two levels.
The approach of both earlier work and this paper is to pretend that the label can be considered as the i-state, the information about the individual that is relevant for predicting its future. The i-state contains information about partners, but not about partners of partners. Implicitly this entails that we use a mean field description of partners of partners. We call this the 'mean field at distance one' assumption. The description of partners of partners is incorporated in an environmental variable, the information about the 'outside world' that is relevant for a prediction of the future of the individual.
A rather special feature of the models considered here is that i-state involves a number of conditionally independent components: the binding sites. An individual has binding sites for partners. Two free binding sites can be joined together to form a partnership between two individuals (see Fig. 1 for an illustration). In graph theory the words 'half-edge' or 'stub' are often used. We think that for static networks these terms capture the essence much better than the word 'binding site'. But the latter provides, in our opinion, a better description for dynamic networks. The fact that our research started with dynamic networks is responsible for our choice of terminology.
It is attractive to model the dynamics of one binding site and next use combinatorics to describe the full i-state. It is precisely this aspect that we did not yet elaborate in Leung et al. (2015) but highlight now. It is precisely this aspect that uncovers the link/relationship between the work of (Lindquist et al. 2011;Leung et al. 2015) on the one hand and the edge-based modelling approach of Volz and Miller et al. (2007, 2008, 2009 on the other hand. Volz and Miller focus on the binding site (=half-edge/stub) and individual level and draw p-level conclusions by a clever use of probabilistic arguments to determine the relevant environmental variables. Lindquist et al. (2011) systematically formulate and analyze the p-level equations. In Leung et al. (2015) we too emphasized the plevel equations, but used the i-level version to derive an expression for R 0 . The link between the two was established by somewhat contrived linear algebra arguments. In the present paper we build our way upwards from binding site-via individual-to population level. One of the secondary aims of this paper is to show that the systematic methodology of PSPM is also very useful when i-state space is discrete, rather than a continuum, and when i-state involves multiple identical components.
In the present paper we do not specify any stochastic model that in the large population limit may lead to the equations that we consider. But we do specify, by way of time-dependent transition rate matrices, various continuous-time Markov chains (where the time-dependence is through the environmental variables). These describe how the states of binding sites change in time. Most of the dependence among binding sites is heuristically captured by the time-dependent coefficients in the transition rate matrix (partly based on the 'mean-field-at-distance-one' assumption). Binding sites that belong to the same individual (the 'owner') experience additional dependence, most notably through transmission of infection: if the owner is infected along one of them, it affects all of them.
When we use the words 'probability', 'expectation', and 'conditioning', we refer to the Markov processes describing binding sites (and hence individuals) for given environmental conditions (i.e. given time-dependent coefficients in the transition rate matrices) and not to the population process as a whole.

Three network cases
Now, consider a network. An epidemic starts when, at some point in time, a small fraction of the population is infected from outside. Our idealized description shifts the 'point in time' towards −∞ while letting the fraction become smaller and smaller. In other words, our story starts 'far back' in time when all individuals are still susceptible (see Appendix 1 for elucidation). We consider three different situations, characterized by the relation between the time scales of, respectively, transmission, partnership dynamics and demographic turnover: I The disease dynamics are fast relative to any partnership-or demographic changes.
The network is static and everyone is susceptible at time t = −∞. II The disease dynamics are on the same time scale as the partnership dynamics, but fast relative to demographic turnover. In this network individuals can acquire and lose partners over time. Everyone is susceptible at time t = −∞. III The disease dynamics and partnership-and demographic changes are on the same time scale. Here the age of an individual matters and we assume that, at birth, an individual enters the population as a susceptible without any partners.
We assume that infection is transmitted from an infectious individual to a susceptible partner at rate β and infectious individuals recover at rate γ (but see Sect. 2.5 for a far more general setting). We also assume that infection does not influence the partnership dynamics or the probability per unit of time of dying in any way.
Each individual in the population is assumed to have a so-called partnership capacity n which denotes the number of binding sites it has (so n is the maximum number of simultaneous partners it may have). Throughout the life of the individual this partnership capacity does not change. An individual with partnership capacity n can be thought of as having n binding sites for partners (in Fig. 1, individuals u, v, and w have partnership capacities 4, 3, and 2, respectively). We call the individual to which a binding site belongs its owner. For the purpose of this paper, we will assume that all individuals have the same partnership capacity n. One can generalize this by allowing individuals to have different partnership capacities; in that case, one needs to average over n in the correct way (see Sect. 2.5 for the static case).

Binding sites
An individual with partnership capacity n is to some extent just a collection of n binding sites. These n binding sites are coupled through the disease status (or death) of their owner. We assume that this is the only manner in which the binding sites of an individual are coupled. As long as the disease status of the owner does not change (and the owner does not die), binding sites behave independently of one another and the 'rules' for changes in binding site states are the same for each binding site. Obviously the latter depends on the network dynamics under consideration (either case I, II, or III). As a port to the world, a binding site can be in one of four states: -0 -free -1 -occupied by a susceptible partner -2 -occupied by an infectious partner -3 -occupied by a recovered partner.
Here (and in the remainder of this introduction) our formulation is precise for case II while sometimes requiring minor adaptations to capture cases I and III.
A key component of the model is the description of the dynamics of a binding site. The state of an individual is specified by listing its disease status and the states of each of its n binding sites. So it makes sense to first consider a binding site as a separate and independent entity and to only take the dependence (by way of a change in the disease status of the owner) into account when we combine n binding sites into one individual.
The case of a susceptible binding site (i.e. a binding site with a susceptible owner) is, as will become clear, far more important than the other cases. This is partly due to our assumption that all individuals start out susceptible, i.e. are susceptible at time t = −∞ (I and II) or at birth (III).
It is tempting to think of the dynamics of susceptible binding sites in terms of jumps between the various states until, sooner or later, transmission occurs along the binding site under consideration (with 'transmission' having no effect at all if the owner was infected earlier along another of its binding sites). There is nothing wrong with this mental image, but it ignores that the binding site under consideration acts as a source of infection once the owner is infected. So one should focus on the first binding site that facilitates transmission and label all binding sites as infected once this happens.
The dynamics of a susceptible binding site are described by a differential equation for the variable x(t) = (x i (t)), i = 0, 1, 2, 3. Here x i can be interpreted as the probability that a binding site is susceptible and has state i at time t, given that its owner does not become infected through one of its other n − 1 binding sites (in other words, by conditioning on the individual not getting infected through its n − 1 other binding sites, the only way the individual could get infected is through the binding site under consideration). In particular, given that its owner does not become infected through one of its other binding sites, is the probability that the binding site is susceptible at time t (or, in other words, that the owner is not infected along this binding site before time t). Accordingly, the probability that an individual is susceptible at time t is equal tō In order to arrive at a closed system of equations for x, we need to go through several steps. The variable x contains information about a partner. Consequently the dynamics of x is partly determined by partners of partners, hence by one or more environmental variables. The 'mean field at distance one' assumption yields expressions for environmental variables in terms of subpopulation sizes (for a given label, the corresponding subpopulation size is the fraction of the population that carries this label). In turn, p-level fractions can be expressed in terms of i-level probabilities. And since a susceptible individual is in essence a collection of n conditionally i.i.d. binding sites, we can use combinatorics to express i-level probabilities in terms of binding-site-level probabilities as incorporated in x.
The exchangeability of the binding sites is broken by the infection event. There is exactly one binding site along which infection took place, viz. the binding site occupied by the individual's epidemiological parent, and for this binding site we know with certainty that it is in state 2 at time of infection t + . We call the binding site through which the change in the owner's disease status occurred the 'exceptional' binding site. The other n − 1 binding sites are i.i.d. and, at time t + , they are distributed according to x(t + )/x(t + ). Recovery (and death) is an event that occurs at a constant rate for an infectious individual so they are independent of binding site states. Therefore, also after recovery, there remains exactly one exceptional binding site, viz. the one through which transmission occurred. See also Fig. 2 for an illustration of the exceptional binding site.

Structure of the paper
In Sects. 2, 3, and 4 below, we will discuss the three network cases I, II, and III separately. For each of the three cases we will explain how the model can be formulated Fig. 2 An illustration of the exceptional binding site. Susceptible, infectious, and recovered individuals are displayed in black, red, and blue, respectively. Three time points in the life of individual u are displayed. Suppose u is susceptible and becomes infected by an infectious partner v at time t + . From that moment on, the binding site along which transmission occurred is the exceptional binding site. This binding site remains exceptional throughout u's life and no other binding site can become exceptional, regardless of whether or not v is still a partner or u is still infectious and described in terms of susceptible binding sites. By considering the susceptible binding site perspective we can write a closed system of only a few equations that fully determine the dynamics of i-level probabilities and p-level fractions. This system is then used to determine epidemiological quantities of interest: R 0 , r , the final size (in cases I and II), and the endemic steady state (in case III). In all three cases, an explicit expression can be given for R 0 . In case I, one can derive a simple scalar equation for the final size. In cases II and III, we could only implicitly characterize the final size and endemic equilibrium, respectively. In Sect. 2 case I of a static network is considered. This is the simplest case among the three. The relative simplicity allows for the derivation of an ODE system for susceptible binding sites directly from the interpretation. This will be the first way in which we formulate the model for this case. But case I will also serve to illustrate the systematic procedure for model formulation in the spirit of PSPM. This systematic procedure allows us to connect the three different levels, viz. (1) binding sites, (2) individuals, and (3) the population, to each other.
In network case I, since it is relatively simple, one can derive a one-dimensional renewal equation from which R 0 , r , and the final size almost immediately follow. This renewal equation will be treated in Sect. 2.5 for a much more general class of infectious disease models than only SIR.
Part of the systematic procedure in cases II and III focuses on infectious binding sites. We use case I to illustrate the model formulation concerning infectious (and recovered) binding sites, even though, for case I these are not needed to obtain a closed system for susceptible binding sites. However, depending on the network features of interest (e.g. fractions of infectious individuals) one may still want to consider infectious (and recovered) binding sites.
In network cases II and III, there are also network dynamics in absence of infection due to partnership changes (and demographic changes). We will only describe the essential characteristics of the network dynamics that we use in this paper. Certainly, much more can be said about the networks in absence of infection (Leung et al. 2012).
Finally, in Sect. 5 we discuss the issues that we have encountered in the three different network cases and pose some open problems. We end the discussion by considering a few generalizations that can be implemented using the systematic model formulation of Sect. 2.2.

Part I: static network 2.1 Model formulation
We derive a closed system of ODE for x purely on the basis of the interpretation of binding sites (without explicitly taking into account i-level probabilities or p-level fractions). The relatively simple setting of a static network allows us to do so. We are able to consider a binding site as a separate and independent entity all throughout its susceptible life. Implicitly, this uses (2.8) below. One can show that the system of ODE for x indeed captures the appropriate large population limit of a stochastic SIR epidemic on a configuration network. This requires quite some work; see (Decreusefond et al. 2012;Barbour and Reinert 2013;Janson et al. 2014).
Consider a susceptible binding site and assume its owner does not become infected through one of its other n − 1 binding sites for the period under consideration. If a susceptible binding site is in state 2, it can become infected by the corresponding infectious partner. This happens at rate β and when it happens, the binding site is no longer susceptible so it 'leaves' the x-system. It is also possible that the infectious partner recovers. This happens at rate γ . Finally, there is the possibility that a susceptible partner of a susceptible binding site becomes infectious (corresponding to a transition from state 1 to state 2). The rate at which this occurs depends on the number of infectious partners that this susceptible partner has. So here we use the mean field at distance one assumption: we average over all possibilities at the p-level to obtain one rate at which a susceptible partner of a susceptible binding site becomes infected. More specifically, we assume that there is a rate βΛ − (t) at which a susceptible partner of a susceptible binding site becomes infected at time t. Here Λ − (t) has the interpretation of the expected number of infectious partners of a susceptible partner of a susceptible individual. Inserting an expectation into a stochastic i-level model in order to lift it to the p-level is reminiscent of the work of Nåsell and Hirsch around 1972, see the book (Nåsell 1985) Then, putting together the various assumptions described above, the dynamics of x is governed by the following system (please note that the environmental variable Λ − is a p-level quantity that we have yet to specify): with 'far past' conditions To express Λ − in terms of x we use the interpretation. Consider a susceptible partner v of a susceptible individual u. Then, since u is susceptible, we know that v has at most n − 1 binding sites that are possibly in state 2 (i.e. occupied by infectious partners). Since v is known to be susceptible, also all its binding sites are susceptible (in the sense that their owner v is). The probability that a binding site is susceptible at (recall (1.1) and note that in case I we have x 0 (t) = 0). The probability that a binding site is in state 2, given that the binding site is susceptible, is x 2 (t)/x(t). Therefore, (2.4) By inserting (2.4) into (2.1) we find that the x-system is fully described by an ODE system in terms of the x-variables only: with 'far past' conditions Remark 1 In the pioneering paper (Volz 2008) an equivalent system of three coupled ODE was introduced to describe the binding-site level of the model. The variables of Volz are connected to our x-system as follows: θ =x, p S = x 1 /x and p I = x 2 /x.

Systematic procedure for closing the feedback loop
Before analyzing (2.5) in the next section, we describe a systematic procedure, consisting of five steps, for deriving the complete model formulation. A key aim is to rederive the crucial relationship (2.4) in a manner that can be extended to the dynamic networks. Thus the present section serves to prepare for a quick and streamlined presentation of the cases II and III in Sects. 3 and 4, respectively. The various steps reveal the relation between binding site probabilities, i-level probabilities and p-level fractions. In addition we introduce some notation.
Step 1. Susceptible binding sites: x-probabilities The first step is to describe the dynamics of x while specifying the environmental variable Λ − only conceptually, i.e. in terms of the interpretation. We then arrive at system (2.1)-(2.2).

Fig. 3
The susceptible partner v of a susceptible individual u has a certain number of infectious partners and the mean of this number equals Λ − V Λ_ U Next, we introduce P (d,k) (t), denoting the fraction of the population with label (d, k). Here k = (k 1 , k 2 , k 3 ) denotes the number of partners of an individual with each of the different disease statuses, i.e. k 1 susceptible, k 2 infectious, and k 3 recovered partners. Furthermore, d ∈ {−, +, * } denotes the disease status of the individual itself, with − corresponding to S, + to I, and * to R.
In the second step, the environmental variable Λ − is, on the basis of its interpretation, redefined in terms of p-level fractions P (−,k) (t).
Step 2. Environmental variables: definition in terms of p-level fractions The mean field at distance one assumption concerns the environmental variable Λ − . This variable is interpreted as the mean number of infectious partners of a susceptible partner v of a susceptible individual u (in terms of Fig. 3: we envisage the uv-connection as a random choice among all such connections). We define it in terms of p-level fractions as follows: . (2.6) Here the sums are over all possible configurations of m and k with m 1 + m 2 + m 3 = n, k 1 + k 2 + k 3 = n. The second factor in each term of this sum denotes the probability that a susceptible partner of a susceptible individual is in state (−, m). The number of infectious partners is then given by m 2 , and we find the expected number of infectious partners Λ − by summing over all possibilities.
In the third step, we let p (−,k) (t) denote the probability that an individual is in state (−, k) at time t. This i-level probability can be expressed in terms of x-probabilities.
Step 3. i-level probabilities in terms of x-probabilities The probability that an individual is, at time t, susceptible with k 1 susceptible, k 2 infectious, and k 3 recovered partners is given by the multinomial expression The solution of the x-system then gives us a complete Markovian description of the i-state dynamics of susceptible individuals.
In this setting of a static network age does not play a role. Therefore, i-level probabilities can immediately be linked to p-level fractions in step 4 below.

Step 4 . p-level fractions in terms of i-level probabilities
The i-level probabilities and p-level fractions coincide, i.e.
In a way, individuals are interchangeable as they all start off in the same state at t = −∞. Finally in the last step, by combining steps 2, 3, and 4, we can express Λ − in terms of the x-probabilities.
Finally, steps 1 to 5 together yield the closed system (2.5) of ODE for x. The dynamics of the 1/2(n + 1)(n + 2) i-level probabilities p (−,k) (t) are fully determined by the system of three ODE for x. We can use this three-dimensional system of ODE to determine r , R 0 , and the final size as we will show in Sect. 2.3. In this particular case of a static network, we can do even better by considering one renewal equation forx. This one equation then allows us to determine the epidemiological quantities as well. This is the topic of Sect. 2.5 where we consider epidemic spread on a static configuration network in greater generality.

The beginning and end of an epidemic: R 0 , r, and final epidemic size
In this section we consider the beginning and end of an epidemic. We first focus on R 0 and r , so on the start of an epidemic. Note that we can very easily find an expression for R 0 from the interpretation: when infected individuals are rare, a newly infected individual has exactly n − 1 susceptible partners. It infects one such partner before recovering from infection with probability β/(β + γ ). Therefore, the expected number of secondary infections caused by one newly infected individual is (2.9) We now rederive R 0 and derive r from the binding site system (2.5).
Note that the p-level fractions P (−,k) (t) can be fully expressed in terms of the binding site level probabilities x i (Eqs. (2.8), (2.7)). Furthermore, the P (−,k) (t) fractions, i.e. the fractions concerning individuals with a − disease status, form a closed system. Therefore, a threshold parameter for the disease free steady state of the binding-site system x is also a threshold parameter for the disease free steady state of the p-level system. (This argument extends to the dynamic network cases II and III in Sects. 3 and 4) Linearization of system (2.5) in the disease free steady statex 1 = 1,x 2 = 0 =x 3 , yields a decoupled ODE for the linearization of the ODE for x 2 . To avoid any confusion, letx 2 denote the linearized x 2 variable. Then the linearization yieldŝ with (for R 0 > 1) 'far past' behaviourx 2 (t) ∼ e (β(n−1)−(β+γ ))t for t → −∞. In particular, the right-hand side of the ODE forx 2 depends only onx 2 .
To illustrate the method used in case II and III in Sects. 3.3 and 4.3, we derive expressions for R 0 and r from a special form of the characteristic equation. Variation of constants for the ODE ofx 2 yieldŝ Substituting the ansatzx 2 (t) = e λt yields the characteristic equation Then there is a unique real root to this equation for λ that we denote by r and call the Malthusian parameter. Evaluating the integral we find that r = β(n −1)−(β +γ ). Likewise, we can derive the expression (2.9) for R 0 by evaluating the integral with λ = 0.
Next, we consider the final epidemic size. We do so by considering the dynamics of x defined in (2.3). Recall (1.2), i.e. the probability that an individual is susceptible at time t, is given byx(t) n . We observe that, by (2.8),x(t) n is also equal to the fraction of susceptible individuals in the population at time t (Alternatively, one can show that k P (−,k) (t) =x(t) n by combining (2.8) and (2.7)). In fact, it is possible to describe the dynamics ofx in terms of onlyx itself. This was first observed in Miller (2011), where the Volz equations of (Volz 2008) were taken as a starting point. The most important observation is the consistency relation (2.10) We can use the interpretation to derive (2.10); x 1 is the probability that a susceptible binding site with owner u is occupied by a susceptible partner v,x n−1 is the probability that v is susceptible given that it is a partner of a susceptible individual u (see also (2.27) below).
Then, using (2.10) together with algebraic manipulation of the ODE system (2.5) (see Miller 2011 for details), one is able to find a decoupled equation forx: (2.11) The fraction of susceptible individuals at the end of the outbreak is determined by the probabilityx(∞). Sincex satisfies (2.11) andx(∞) is a constant, we find that necessarilyx(∞) is the unique solution in (0, 1) of The final epidemic size is given by In Sect. 2.5 we show that one can actually describe the dynamics of the probabilitȳ x for deterministic epidemics on configuration networks for a much larger class of submodels for infectiousness. The SIR infection that we consider here is a very special case of the situation considered in Sect. 2.5. There we show that it is possible to derive a renewal equation forx. The final size equation is then obtained by simply taking the limit t → ∞. We highly recommend reading Sect. 2.5 to understand the derivation of the renewal equation forx based on the interpretation of the model (with a minimum of calculations).

After susceptibility is lost
In the preceding section we have seen that the x-system (2.5) for susceptible binding sites is all that is needed to determine several epidemiological quantities of immediate interest. On the other hand, we might not only be interested in the fraction (1.2) of susceptibles in the population, but also in the dynamics of i-level probabilities p (d,k) (t) (and likewise p-level fractions P (d,k) (t) given by (2.8)) for d = +, * .
So what happens after an individual becomes infected? We work out the details for infectious individuals and only briefly describe recovered individuals. Again, we are able to formulate the model following steps 1-5 of Sect. 2.2 (where the word 'susceptible' should be replaced by 'infectious' or 'recovered' whenever appropriate and step 3 should be replaced by a slightly different step 3', but we will come back to this later on in this section). But now we need to take into account the exceptional binding site, i.e. the binding site through which infection was transmitted to the owner (see also Fig. 2).
In step 1 one considers the dynamics of infectious binding sites, i.e. binding sites having infectious owners. Suppose that the owner became infected at time t + and that it does not recover in the period under consideration. Let y e i (t | t + ) denote the probability for the exceptional binding site to be in state i at time t, i = 1, 2, 3. Similarly, y i (t | t + ) denotes the probability for a non-exceptional binding site to be in state i at time t, i = 1, 2, 3. Here the probabilities are defined only for t ≥ t + . Note V Λ + U Fig. 4 The susceptible partner v of an infectious individual u has a certain number of infectious partners and the mean of this number equals Λ + (note that this number is always larger or equal to 1 since u is a partner) that y and y e are probability vectors, i.e. the components are nonnegative and sum to one.
Instead of 'far past' conditions we now have to take into account the distribution of binding site states at time of infection t + . Whether or not an infectious binding site is exceptional has an influence on the state it has at epidemiological birth. Indeed, the exceptional binding site is in state 2 at time t + with probability 1, while the distribution of the state of a non-exceptional binding site at time t + is given by we have boundary conditions (2.13) The mean field at distance one assumption again plays a role. Here, we need to deal with the environmental variable Λ + that is defined as the expected number of infectious partners of a susceptible partner of an infectious individual (in terms of Fig. 4: we envisage the uv connection as a random choice among all such connections, compare with Fig. 3). We can redefine Λ + in terms of p-level fractions P (−,k) for susceptible individuals: . (2.14) In particular, once again, Λ + can be expressed in terms of x by combining steps 2, 3, and 4. Using (2.14), (2.8), and (2.7) we find that (alternatively, one can find the same expression for Λ + in terms of x-probabilities directly from the interpretation, exactly as before in the case of Λ − ).
The rates at which changes in the states (1, 2, 3) of infectious binding sites occur is the same for each binding site, including the exceptional one. There is a rate γ at which an infectious partner of an infectious binding site recovers (this corresponds to a change in state from 2 to 3). And there is a rate at which a susceptible partner of an infectious binding site becomes infected (either along the binding site under consideration or by one of its other infectious partners) corresponding to a change in state from 1 to 2. The rate at which this occurs is βΛ + where Λ + is defined by (2.14) and hence (2.15).
Recall that we condition on the infectious binding site under consideration not recovering, therefore, these are all state changes that can occur. So we find that the dynamics of y and y e are described by the same ODE system and case specific boundary conditions (2.13). Observe that this means that y e 1 (t | t + ) = 0 for all t ≥ t + .This also immediately follows from the interpretation: at time t + , the binding site is occupied by an infectious partner, the network is static, and an infectious individual can not become susceptible again. Next, we turn to infectious individuals. Compared to susceptible i-level probabilities, it is more involved to express infectious i-level probabilities in terms of y e -and y-probabilities. Therefore, we first consider conditional i-level probabilities before finding an expression for the unconditional probabilities. We replace step 3 by step 3'.

Step 3' Infectious i-level probabilities in terms of y and y e
We let φ (+,k) (t | t + ) denote the probability that an infectious individual, infected at time t + , is in state (+, k) at time t, given no recovery. As in the case of a susceptible individual, this is given by a multinomial expression (note that there is one exceptional binding site which is either in state 2 or 3) (2.17) Note that φ (+,k) (t | t + ) = 0 for k = (n, 0, 0), i.e. for all t ≥ t + at least one partner is not susceptible.
A susceptible individual becomes infected at time t + if infection is transmitted to this individual through one of its n binding sites. Infection is transmitted at rate β. Therefore, the force of infection at time t + , i.e. the rate at which a susceptible individual becomes infected at time t + , equals βn x 2 x (t + ) and consequently the incidence at time t + , i.e. the fraction of the population that becomes, per unit of time, infected at time t + , equals βn (recall thatx n is the fraction of the population that is susceptible). Furthermore, an infectious individual that is infected at time t + is still infectious at time t if it does not recover in the period (t + , t). Since the infectious period of an individual is assumed to be exponentially distributed with rate γ , the probability that We then find an expression for the unconditional i-level probabilities p (+,k) (t) that a randomly chosen individual is in state (+, k) at time t in terms of infectious binding site probabilities and the history of susceptible binding site probabilities: is given by (2.17). The i-level probabilities p (+,k) (t) are lifted to the p-level by (2.8).
In this way we can use infectious binding sites as building blocks for infectious individuals. We see that y and y e explicitly depend on the dynamics of x through the boundary conditions (2.13) and the environmental variable Λ + (2.15). In addition, x 2 plays a role in determining the time of infection of an individual.
Remark 3 Similar to the ODE system for − individuals considered in Remark 2, one obtains the p-level ODE system by differentiation of (2.20) and use of (2.16), (2.7) and (2.8). In doing so, one obtains a system of 1/2(n + 1)(n + 2) ODE for the p-level fractions concerning individuals with a + disease status: Lindquist et al. 2011, Eq. (13)).
In case of recovered individuals, one considers their binding sites and first conditions on time of infection t + and time of recovery t * . Again one needs to distinguish between the exceptional and the non-exceptional binding sites. The dynamics of recovered binding sites are described by taking into account the mean field at distance one assumption for the mean number Λ * of infectious partners of a susceptible partner of a recovered individual. Boundary conditions are given by the y(t * | t + ) and y e (t * | t + ) for non-exceptional and exceptional binding sites, i.e.
The dynamics for z and z e can be described by a system of ODE identical to the ODE systems for y and y e , but with Λ + replaced by Λ * . The environmental variable Λ * is given by . (2.21) By combining (2.21) with (2.8) and (2.7) we find ( 2.22) We find an expression for the probability ψ ( * ,k) (t | t + , t * ) that a recovered individual, infected at time t + and recovered at time t * , is in state ( * , k) at time t ≥ t * , in terms of z and z e probabilities for recovered binding sites with the same reasoning as for φ (+,k) (t | t + ) (one can simply replace φ by ψ, y i by z i , and y e i by z e i in (2.17)). Then, to arrive at an expression for the unconditional probability p ( * ,k) (t), we again need to take into account the incidence βnx 2x n−1 (t + ) at t + . The probability that recovery does not occur in the time interval (t + , t * ) is given by e −γ (t * −t + ) and the rate at which an infectious individual recovers is γ , therefore where the first equality in (2.23) follows from (2.8).

The renewal equation for the Volz variable
So far we dealt with the SIR situation, where an individual becomes infectious immediately upon becoming infected and stays infectious for an exponentially distributed amount of time, with rate parameter γ , hence mean γ −1 . During the infectious period any susceptible partner is infected with rate (=probability per unit of time) β.
Here we incorporate randomness in infectiousness via a variable ξ taking values in a set Ω according to a distribution specified by a measure m on Ω. This sounds abstract at first, but hopefully less so if we mention that the SIR situation corresponds to with ξ corresponding to the length of the infectious period. In this section we only consider the setting where the 'R' characteristic holds, i.e. after becoming infected, individuals can not become susceptible for infection any more.
In order to describe how the probability of transmission to a susceptible partner depends on ξ , we need the auxiliary variable τ corresponding to the 'age of infection', i.e. the time on a clock that starts when an individual becomes infected. As a key model ingredient we introduce π(τ, ξ ) = the probability that transmission to a susceptible partner happens before τ, given ξ.
In the SIR example we have It is important to note a certain asymmetry. On the one hand, there is dependence in the risk of infection of partners of an infectious individual u. Their risk of getting infected by u depends on the length of the infectious period of u (and, possibly, other aspects of infectiousness encoded in ξ ). On the other hand, if u is susceptible, the risk that u itself becomes infected depends on the length of the infectious periods of its various infectious partners. But these partners are independent of one another when it comes to the length of their infectious period (see also Diekmann et al. 2013, Sect. 2.3 'The pitfall of overlooking dependence'). In particular, the probability that an individual escapes infection from its partner, up to at least τ units of time after the partner became infected, equals ( 2.24) For the Markovian SIR example (2.24) boils down to a formula that can also be understood in terms of two competing events (transmission versus ending of the infectious period) that occur at respective rates β and γ . As in Diekmann et al. (1998a) and earlier subsections, we consider a static configuration network with uniform degree distribution: every individual is connected to exactly n other individuals. At the end of this section we shall formulate the renewal equation for arbitrary degree distribution. In Diekmann et al. (1998a) an expression for R 0 and equations for both final size and the probability of a minor outbreak were derived. In addition, it was sketched how to formulate a nonlinear renewal equation for a scalar quantity, but the procedure is actually that complicated that the resulting equation was not written down.
The brilliant idea of Volz (2008) is to focus on the variable θ(t) corresponding to the probability that along a randomly chosen partnership between individuals u and v no transmission occurred from v to u before time t, given that no transmission occurred from u to v (see also Fig. 5 for a schematic representation). Here one should think of 'probability of transmission' as being defined by π (and hence F) and not require that the individual at the receiving end of the link is indeed susceptible (though, if it actually is, or has been, infectious, the condition of no transmission in the opposite direction is indeed a nontrivial condition).
The variable θ corresponds tox introduced in Sect. 2.1 and therefore we use that symbol also in this section. We reformulate (2.3) as x(t) = prob{a binding site is susceptible at time t | its owner does not become infected through one of its other binding sites before time t} (2.26) 5 Volz focused on the variable θ(t) corresponding to the probability that along a randomly chosen partnership between individuals u and v no transmission occurred from v to u before time t, given that no transmission occurred from u to v Fig. 6 Schematic representation ofx. In this figure, the binding site under consideration is indicated in green. Its owner has three binding sites in total. It is given that no transmission occurs through its other two binding sites (see also Fig. 6). There is an underlying stochastic process in the definition forx that we have not carefully defined here. Yet we shall use the words from the definition to derive a consistency relation that takes the form of a nonlinear renewal equation for x(t). The renewal equation describes the stochastic process starting 'far back' in time when all individuals were still susceptible. A precise mathematical definition and an in-depth analysis of a more general stochastic process (that allows the contact intensity between two connected individuals to depend on the number of binding sites of both of them) can be found in Barbour and Reinert (2013). See (Karrer and Newman 2010, Sect. V) for a different way of specifying initial conditions.
To derive the consistency relation forx(t) we shift our focus to the partner that occupies the binding site under consideration. For convenience we call the owner of the binding site under consideration u and the partner that occupies this binding site v. Then, given that u does not become infected through one of its n − 1 other binding sites, u is susceptible at time t if (1) v is susceptible at time t or (2) v is not susceptible at time t but has not transmitted infection to u up to time t.
We begin by determining (1). Given its susceptible partner u, individual v is susceptible if its n − 1 other binding sites are susceptible. Conditioning on its n − 1 other binding sites not transmitting to v, a binding site of v is susceptible at time t with probabilityx(t). Therefore, given susceptibility of partner u, v is susceptible at time t with probabilityx (t) n−1 .
( 2.27) This just repeats the consistency relation (2.10) x 1 =x n−1 stating that the probability x 1 that a susceptible binding site is occupied by a susceptible partner is equal to the probabilityx n−1 that a partner of a susceptible individual is susceptible. Next, suppose that v gets infected at some time η < t, then u is not infected by v before time t if no transmission occurs in the time interval of length t − η. The expression (2.27) has as a corollary that the probability per unit of time that v becomes infected at time η equals Noting that the probability of no transmission to u in the time interval (η, t) is F(t −η) we conclude that necessarily, (2.28) Finally, by integration by parts, we obtain the renewal equation For a configuration network with general degree distribution ( p n ) for the number of binding sites n of an individual, exactly the same arguments hold. But now there is randomness of the partnership capacity n of susceptible partners. This leads to the renewal equation (compare with (2.29)) The solutionx(t) = 1, −∞ < t < ∞, of (2.30) corresponds to the disease free situation. If we putx(t) = 1 − h(t) and assume h is small, we easily deduce that the linearized equation is given by The corresponding Euler-Lotka characteristic equation reads If we evaluate the right hand side of (2.32) at λ = 0, we obtain cf. Diekmann et al. (2013, Eq. (12.32), p. 294). In short, the relevant characteristics of the initial phase of an epidemic outbreak are easily obtained from the linearized renewal equation (2.31) (see Pellis et al. 2015 for a study of the Malthusian parameter, i.e. the real root of (2.31)).
In the case that F is given by (2.25), the RĒ can be transformed into an ODE forx by differentiation: In the special case of Sects. 2.1-2.4, we have p n = 1 and p k = 0 for all k = n so g(x) = x n−1 and we recover (2.11).
As explained in (Diekmann et al. Finite dimensional state representation of linear and nonlinear delay systems. Submitted), the natural generalization of (2.25) assumes that F is of the form where, for some m ∈ N, β and V are non-negative vectors in R m while Σ is a positive-off-diagonal (POD) m × m matrix. (In Diekmann et al. Finite dimensional state representation of linear and nonlinear delay systems. Submitted) it is explained how this relates to the method of stages from queueing theory, see (Asmussen 1987), but with the kernel F not only incorporating the progress of infectiousness of the focus individual but also the fact that any partner can be infected at most once. Concerning infectiousness, one might think of SEIR, SEI 1 I 2 , etcetera.) If F is given by (2.34), the variable and, since (2.30) can be rewritten as the Eq. (2.35) is a closed system once we replacex at the right hand side of (2.35) by the right hand side of (2.36) So one can solve/analyze (2.35) and next use the identity (2.36) to draw conclusions aboutx. We conclude that various ODE systems as derived in Miller et al. (2012) are subsumed in (2.30) and can be deduced from (2.30) by a special choice of F and differentiation.

Part II: dynamic network without demographic turnover
In Sect. 2, only one environmental variable Λ − is involved in the specification of the dynamics of the susceptible binding sites. In dynamic networks, additional environmental variables play a role. Notably, we have to specify the (probability distribution of the) disease status of a new partner. Before formulating the model for susceptible binding sites, we first consider the network itself in Sect. 3.1. This is needed in order to determine the appropriate 'far past' conditions of the susceptible binding site system. In Sect. 3.2, the model formulation is divided into three subsections. First, we formulate the model in terms of susceptible binding site probabilities x by following the scheme of five steps presented in Sect. 2.2. This allows us to express in terms of x those environmental variables that are defined in terms of susceptible p-level fractions P (−,k) . We then consider infectious and recovered binding site systems and these allow us to express the other environmental variables in terms of (the history of) x as well.

Network dynamics
Binding sites are either free or occupied. We denote the fraction of free binding sites in the population by F. We assume that a binding site that is free becomes occupied at rate ρ F, while an occupied binding site becomes free at rate σ . Similar to (Leung et al. 2012) (set μ = 0), we find that F satisfies the ODE So we find that F converges to a constant for t → ∞. Therefore, we assume that the fraction of free binding sites is constant, and this constant is again denoted by the symbol F. Then F satisfies Although we could give an explicit expression in terms of σ and ρ for F, we prefer to state the more useful identity (3.1) that, viewed as an equation, has F as its unique positive root. The network structure, although dynamic, is stable. A randomly chosen binding site (in the pool of all binding sites) is free with probability F and occupied by a partner with probability 1 − F. Later on we shall use that, given that a binding site is free with probability F at time τ , the probability that a binding site is free at time ξ + τ is F (and the probability that it is occupied at time ξ + τ is 1 − F).
Finally, later on in Sect. 3.2.2, we also need the probability ϕ 1 (ξ ) that a binding site is free at time ξ + τ if it is occupied at time τ . Note that, by the Markov property, this probability only depends on the length ξ of the time interval. Since ϕ 1 (ξ ) is the unique solution of the initial value problem: where we used (3.1) in the second equality.

Susceptibles
We describe the dynamics of susceptible binding sites in terms of x-probabilities. Long ago in time, by assumption, all individuals (and therefore binding sites) are susceptible. In accordance with Sect. 3.1 the fraction of free and susceptible binding sites is equal to F and the fraction of susceptible binding sites occupied by susceptible partners is equal to 1 − F, i.e. we have 'far past' conditions The environmental variables F and Λ − are p-level quantities that we have yet to specify. Putting together the various assumptions described above, the dynamics of x is governed by the system: with 'far past' conditions (3.3), and Next, in step 2, we define the environmental variables in terms of p-level fractions. The definition (2.6) of Λ − in terms of p-level fractions carries over. We define the fractions of free binding sites in terms of p-level fractions as follows: where the sum is over all possible configurations of k with 0 ≤ k 1 + k 2 + k 3 ≤ n.
In step 3 we define the i-level probabilities p (−,k) (t) in terms of the x probabilities by using the conditional independence of binding sites: (3.7) As in the static network case I, the i-level probabilities coincide with the p-level fractions, i.e. (2.8) holds. This is step 4 in our model formulation.
Then, in step 5, we can express the environmental variables Λ − and F − in terms of x-probabilities. By combining (3.7) with (2.8) and (2.6), we again find (2.4) to hold (only now the x are defined by the system of ODE (3.4)). By combining (3.7) with (2.8) and (3.6) we find that exactly as the interpretations of F − and x would suggest. Before we can specify F + and F * in terms of (the history of) x we need to define p-level fractions P (+,k) (t) and P ( * ,k) (t). We do so in the next section where we turn to infectious and recovered binding site systems.

After susceptibility is lost
If an individual becomes infected at time t + , the binding site through which infection is transmitted is from that point on 'exceptional'. Then, given that its owner became infected at time t + and that it does not recover for the time under consideration, we consider an infectious binding site. Let y e i (t | t + ) denote the probability that the exceptional binding site is in state i at time t and y i (t | t + ) this same probability for a non-exceptional binding site.
As in Sect. 2.4, the exceptionalness plays a role only in the states at epidemiological birth, i.e. at time t + . The exceptional binding site is with probability one in state 2 at time t + . The states of all other binding sites are distributed according to x(t + )/x(t + ). Therefore, we put boundary conditions Since the individual does not recover in the period under consideration, the infectious binding sites behave independently of one another.
The dynamics of y e and y are both governed by the system Note that there is no rate γ of leaving the infectious state by the owner as we condition on the owner remaining infectious in the period under consideration (this is compensated by a factor e −γ (t−t + ) in various integrals below). Furthermore, note that, contrary to the network case I of Sect. 2.4, the exceptional binding site can lose its epidemiological parent by separation. Therefore y e 1 (t | t + ) > 0 for t > t + . Next, similarly to Sect. 2.4, by combinatorics (but now probabilities y e 0 and y e 1 are not equal to zero for t ≥ t + ), we find that the probability φ (+,k) (t | t + ) that an individual, infected at time t + , is in state (+, k) at time t ≥ t + is given by y e 0 y k 0 −1 0 y k 1 1 y k 2 2 y k 3 3 + k 1 n y e 1 y k 0 0 y k 1 −1 1 y k 2 2 y k 3 3 + k 2 n y e 2 y k 0 0 y k 1 1 y k 2 −1 2 y k 3 3 + k 3 n y e 3 y k 0 0 y k 1 1 y k 2 2 y k 3 −1 3 (t | t + ). (3.10) The probability P (+,k) (t) that a randomly chosen individual is in state (+, k) at time t is obtained by taking into account the time of infection t + and the probability (2.19) that an individual has not recovered time t − t + after infection. The definition (2.18) for the incidence carries over (but now with the x defined by the ODE system (3.4)). So By combining (2.8) and the expression (3.11) for P (+, f O f f k) in terms of y and y e we can redefine F + in terms of the history of x as we will show now. First of all, combining (3.10), (3.11) and (3.6) we express F + in terms of y and y e : n−1 + (n − 1)ȳ e y 0ȳ n−2 (t | t + )dt + .
Next, we consider the probabilities y 0 , y e 0 . Note that with ϕ 1 given by (3.2). The dynamics of y 0 are described in terms of y 0 and the history of x 0 /x (by means of the boundary condition). We can solve for y 0 . This yields (note that time of infection t + matters in this probability and not only the length t − t + of the time interval). We can further simplify (3.12) to obtain which only depends on the model parameters and past probabilities x i for susceptible binding sites. We can use the consistency condition for the total fraction of free binding sites: (3.14) to express F * in terms of the history of x (use (3.8) for F − and (3.13) for F + ). So this specifies all environmental variables for the susceptible binding site system x in terms of (the history of x).
Next, similar to case I of Sect. 2.4, we consider recovered individuals and their binding sites. Suppose that the infectious individual, that was infected at time t + , recovers at time t * . After recovery, we still distinguish between the exceptional binding site and the n − 1 other binding sites. We introduce probabilities z e i (t | t + , t * ) and z i (t | t + , t * ) for recovered binding sites. The y and y e probabilities yield the conditions for z and z e at time t = t * , i.e.
The dynamics of z and z e are described by the system of ODE for y and y e , with the mean field at distance one quantity Λ + replaced by Λ * where Λ * is defined in terms of p-level fractions by (2.21) and hence is given by (2.22) in terms of x-probabilities.
Let ψ ( * ,k) (t | t + , t * ) denote the probability that a recovered individual is in state ( * , k) given that it was infected at time t + and recovered at time t * . Then ψ ( * ,k) (t | t + , t * ) can be expressed in terms of z and z e by replacing φ in (3.10) by ψ, y i by z i , and y e i by z e i . The unconditional probability p ( * ,k) (t) is then obtained by taking into account time of infection t + and recovery time t * : which, by relation (2.8), is equal to the p-level fraction P ( * ,k) (t). Note that we can also use this definition of P ( * ,k) (t) to define F * in terms of x similar to the way we did for F + in (3.13).

One renewal equation or a system of six ODE, whatever you like
We ended the model formulation in Sect. 3.2.1 by defining the environmental variables Λ − and F − in terms of x (Eqs. (2.4), (3.8)). Subsequently, in Sect. 3.2.2, by considering infectious binding site probabilities y, and y e , we also defined F + and F * in terms of x (Eqs. (3.13), (3.14)). Combining these formulas, we find that the system describing the dynamics of susceptible binding sites is given by: with F + given by (3.13) and with 'far past' condition The ODE (3.16) for x together with the expression (3.13) for F + yields a closed system of five equations. By substituting expression (3.13) in system (3.16), one can view (3.16) as a system of four delay differential equations. The dynamics of the 1/6(n + 1)(n + 2)(n + 3) i-level probabilities (hence p-level fractions) p (−,k) for susceptible individuals are fully determined by this set of four delay differential equations (regardless of n).
Alternatively, we can view the solution x(t) of (3.16)-(3.17) as fully determined by F + | (−∞,t] . Interpreting x 2 ,x, and x 0 at the right hand side of (3.13) in this manner, we arrive at the conclusion that the dynamics are fully determined by a single renewal equation for F + .
One may prefer a system consisting only of ODE rather than a delay system. We can in fact reason directly in terms of the interpretation to derive an ODE for F + . In order to do so, we first consider the fraction I (t) = k P (+,k) of infectious binding sites in the population. This fraction decreases when infecteds recover. Infecteds recover at a constant rate γ . The fraction I increases when a susceptible individual becomes infected so there is the positive term (2.18) in the ODE for I (combine (3.7) with (2.8) and (2.18), the x are defined by the ODE system (3.4))). We find that the dynamics of I are described by the following ODE: with 'far past' condition I (−∞) = 0. Next, we consider F + . Any infectious owner recovers at constant rate γ . In addition, partnership formation and separation affect the fraction of free infectious binding sites. There is a rate ρ F at which free binding sites become occupied. The fraction of infectious binding sites that are occupied is given by I − F + and the rate at which these binding sites become free is σ . Then, finally, a susceptible individual with k 2 infectious partners becomes infected at rate βk 2 , taking into account all 0 ≤ k 2 ≤ n we find probability per unit of time βnx 2x n−1 at which a susceptible individual becomes infected. The probability that a non-exceptional binding site is free and susceptible upon infection is x 0 /x, so the expected fraction of free binding sites created upon infection of a susceptible individual is 1 n (n − 1)x 0 /x. Hence there is a flow β(n − 1)x 0 x 2x n−2 into F + . We have the following ODE for F + : with 'far past' condition F + (−∞) = 0. Alternatively, we can derive the ODE (3.19) for F + by differentiating (3.13) with respect to t. Note that we can express I in terms of x by first expressing it in terms of y and y e (similar to F + in Sect. 3.2.2). This yields n−1 (t + )dt + . The combination of (3.16) with (3.18) and (3.19) yields a six-dimensional closed system of ODE. (Compare with the slightly different but related network model called the 'dormant contacts' model of (Miller et al. 2012). Presumably (3.16), (3.18), (3.19) is a transformed but equivalent version of their system (3.11)-(3.19) for constant n, i.e. for a degree distribution concentrated in one point).
Both (3.13) and (3.16)-(3.19) can be used to represent the system. In terms of the number of equations, it does not matter too much which system one considers. In the first case, one renewal equation is needed compared to six ODE in the second case. In both formulations one can determine r and R 0 with not too much effort. In Sect. 3.3 below, we will use a pragmatic mixture. This gives us a way of determining r and R 0 that prepares for the characterization of r and R 0 in case III in Sect. 4.3 (where a model formulation in terms of only ODE becomes troublesome).

The beginning and end of an epidemic: R 0 , r and final size
First, just as in case I, the final size is given by But while in case I we derived a simple scalar equation forx(∞) ((2.12) or (2.33)), depending explicitly on the parameters, we did not, despite fanatical efforts, manage to derive such an equation from the implicit characterization by (3.16), (3.18), (3.19); see also Appendix 1.
Next, in the rest of this section, we use the binding site level system (3.16)-(3.19) to consider the beginning of an epidemic and determine R 0 and r . The point here is not only to use (3.16)-(3.19) to find threshold parameters but to find threshold parameters with their usual interpretation of R 0 and r .
Using the same arguments as in network case I of Sect. 2, we find that a threshold parameter for the disease free steady state of system (3.16)-(3.19) on the binding site level is also a threshold parameter for the disease free steady state of the p-level system.
The disease free steady state of (3.16) is given byx 0 = F,x 1 = 1− F,x 2 = 0 =x 3 . Linearization in this state yields a decoupled system of equations for the linearized x 2 and F + equations. We letx 2 andF + denote the variables in the linearization in the disease free steady state. Note that, in the disease free steady stateỹ 0 in the disease free steady state, the probability that an infectious binding site is free at time t given that it is free at time t + is equal to the probability F that a randomly chosen binding site is free. Then which can be viewed as a linear delay differential equation forx 2 . In order to obtain an informative version of the corresponding characteristic equation, we rewrite it as a renewal equation forx 2 . Variation of constants for the ODE forx 2 yields: Substituting (3.20) into this expression yields the renewal equation forx 2 : where the rearrangement of the terms in the integrals is in preparation for the interpretation). Next, we substitute the ansatzx 2 (t) = e λt , and obtain the characteristic equation There is a unique real root to (3.21) and this root is by definition the Malthusian parameter r . We define R 0 = ∞ 0 k(ξ )dξ . Then sign(r )=sign(R 0 − 1), and we find that R 0 is a threshold parameter with threshold value one for the stability of the disease free steady state, We can evaluate the integrals and find an explicit expression for R 0 . However, the interpretation is easier in the form it is written now. First of all, consider a newly infected individual u. Individual u transmits infection to a susceptible partner with probability ∞ 0 βe −(σ +β+γ )ξ dξ = β/(β + σ + γ ). By multiplying this probability with the expected number of susceptible partners u has at epidemiological birth plus the expected number of susceptible partners u acquires during its infectious period after epidemiological birth, we obtain R 0 . As we will explain now, these are exactly the two terms in {. . .} of (3.22).
The mean number of susceptible partners of u at epidemiological birth is (n − 1)(1 − F) (note that, in addition to the susceptible partners, u has (n − 1)F free and 1 exceptional binding site). This is the first term in {. . .} of (3.22). We are left with determining the expected number of susceptible partners u acquires after epidemiological birth. This goes as follows. At time τ after u became infected, u has not recovered yet with probability e −γ τ . The exceptional binding site of u is free at time τ with probability ϕ 1 (τ ) (see (3.2)). Each of the n − 1 non-exceptional binding sites of u are free with probability F regardless of whether they were free or occupied at epidemiological birth (recall Sect. 3.1). Note that a free binding site becomes occupied by a susceptible partner at rate ρ F (at the beginning of the epidemic). Integrating over all possible lengths τ > 0 of the infectious period, we find that ∞ 0 e −γ τ ρ F ϕ 1 (τ ) + (n − 1)F dτ is the expected number of additional susceptible partners of u in its infectious period after epidemiological birth.
Note that we made the distinction of the susceptible partners at and after epidemiological birth of u but what really matters is the total number of susceptible partners in the infectious period of u. So really, we did not need to make any distinction between at and after epidemiological birth. But this distinction is essential in Sect. 4.3 of case III. The distinction here serves both to illustrate this difference with case III and as a preparation for case III.
Finally, in the same spirit, we would like to mention that rather than taking the perspective of an infectious individual/binding site, we can also take the perspective of a susceptible binding sites 'at risk' of infection, i.e. susceptible binding sites occupied by infectious partners, and interpret R 0 in that way. In the present context this does not change much. Therefore we refrain from elaborating. We leave this for Sect. 4.3 of case III where this different perspective leads to a major simplification compared to the 'standard' perspective of infectious binding sites that we took here.

Part III: dynamic network with demography
In this part, the network is not only dynamic due to partnership formation and separation but also due to demographic turnover. We assume that there is a constant per capita death rate μ and a constant population birth rate so that the population size is in equilibrium and the age of individuals is exponentially distributed with parameter μ. At birth, an individual does not have any partners. Details are presented in Leung et al. (2012).

Network dynamics
In a world with demographic turnover, next to calendar time, also age matters. We keep track of both age a and time of birth t b of an individual (calendar time is then given by t = a + t b ). When we speak about the age and time of birth of a binding site, we mean the age and time of birth of its owner. Often, we assume that the owner of a binding site does not die in the period under consideration. By assumption, at age zero, a binding site is free. A free binding site becomes occupied at rate ρ F where F denotes the total fraction of free binding sites in the population. This F is assumed to be constant (see Leung et al. 2012 for the justification) and satisfies (compare with (3.1)). If the binding site is occupied, then it becomes free at rate σ + μ where σ and μ represent separation and death of partner, respectively.
In this section we will also make use of the following binding site probabilities (where, as usual, we condition on the owner not dying in the period under consideration). We let ϕ 0 (a) denote the probability that a binding site is free at age a + α, given that it was free at age α, and ϕ 1 (a) denotes the probability that a binding site is free at age a + α, given that it is occupied at age α. Note that, by the Markov property, these probabilities only depend on the time interval a (recall that F is constant). The dynamics of ϕ i as a function of a is described by with initial conditions, respectively, The explicit expressions for the ϕ i are given by See also (Leung et al. 2012, Eq. (10)) (where (a) can be identified with 1 − ϕ 0 (a)) and (Leung et al. 2015, Eq. (67)) (where 0 (t) and 1 (t) can be identified with 1−ϕ 0 (t) and 1 − ϕ 1 (t), respectively). Furthermore, we have the identity (use (4.1)), expressing that a randomly chosen binding site is free with probability F. So, according to Bayes' Theorem, the probability density function of the age of (the owner of) a free binding site is given by Similarly, the probability density function of the age of (the owner of) a randomly chosen occupied binding site is (in view of the derivation of a formula for R 0 in Sect. 4.3 below, we remark that π 0 and π 1 should be compared to probability distributions q and Q, respectively, in Leung et al. (2015); the difference is that q and Q concern the number of partners while π 0 and π 1 concern the age; the probability distributions, however, provide the same information).

Susceptibles
Demography does not give rise to any additional environmental variables, we still deal with the mean field at distance one variable Λ − , and the fractions F d of free binding sites with disease status d, d ∈ {−, +, * }. We follow the steps 1-5 of Sect. 2.2. In step 1 we consider x-probabilities. Consider a susceptible binding site, born at time t b , and suppose that its owner, for the period under consideration, does not die and does not become infected through one of its other n − 1 binding sites. Let F = (F − , F + , F * ). The dynamics of x as a function of age are described by the following system of equations: An individual is assumed to be susceptible without any partners at birth (and therefore the same applies to all its binding sites). So we have the birth conditions Given the environmental variables F and Λ − , we can formally view x(a | t b ) as a function of the environmental variables: We now define the environmental variables in terms of p-level fractions. Note that Λ − has the exact same interpretation as in network cases I and II. It should therefore come as no surprise that the definition of Λ − in terms of p-level fractions is again (2.6). The fractions of free binding sites with disease status d are again defined by (3.6). This is step 2.
Next, in step 3, we define the i-level probabilities p (−,k) (a | t b ) in terms of x. Given that an individual was born at time t b and does not die in the period under consideration, the probability that an individual is at age a in state (−, k) is given by the multinomial expression (compare with Eq. (3.7) and note that we condition on the survival of the individual).
In step 4 we relate p-level fractions P (d,k) . In order to do so, we use the stationary age distribution with density a → μe −μa . The fraction of the population that is in state (d, k) at time t is obtained by adding all individuals in that state that are born before time t and are still alive at time t. We find that In step 5, we express the environmental variables Λ − and F − in terms of x. This can be done by combining (4.10) and (4.11) with (2.6) (for Λ − ) or (3.6) (for F − ). We find that and (4.13) In order to complete step 5 (expressing the environmental variables F + and Λ − in terms of x) we need to consider infectious and recovered binding sites.

After susceptibility is lost
Consider a binding site that was born at time t b and infected at age a + and its owner remains alive and infectious for the period under consideration. Note that age a + for this individual corresponds to calendar time t + = t b + a + . Let y e i (a | t b , a + ) denote the probability that the exceptional binding site is in state i at age a and y i (a | t b , a + ) the same probability for a non-exceptional binding site.
Then, at age a + , the exceptional binding site is for certain in state 2, while the other n − 1 binding site states are distributed according to x(a + | t b )/x(a + | t b ): , The dynamics of infectious binding sites are described by: with Again, there is no rate γ in M + (F, Λ + ) of the owner leaving the system of infectious binding sites as we assume that infectious binding sites remain infectious in the period under consideration (which is compensated by a factor e −γ (a−a + ) in various integrals below). In (4.14) we can consider Λ + as 'known'. Indeed, by combining (2.14) with (4.11) and (4.10), we can express Λ + in terms of x as follows: We now set out to derive an expression for F + . The probability φ (+,k) (t b , a | a + ) that an individual, born at time t b and infected at age a + , is in state (+, k) at age a ≥ a + is given by y e 0 y k 0 −1 0 y k 1 1 y k 2 2 y k 3 3 + k 1 n y e 1 y k 0 0 y k 1 −1 1 y k 2 2 y k 3 3 + k 2 n y e 2 y k 0 0 y k 1 1 y k 2 −1 2 y k 3 3 + k 3 n y e 3 y k 0 0 y k 1 1 y k 2 2 y k 3 −1 3 (a | t b , a + ). (4.15) The contribution to the incidence of individuals of age a + , born at time t b and alive for the period under consideration, is given by where the reasoning is similar to cases I and II. Then, taking into account all possible ages of infection 0 ≤ a + ≤ a, and the probability that as yet recovery did not occur, the probability that an individual, born at time t b , is in state (+, k) at age a is given by The p-level fractions P (+,k) (t) at time t are obtained through relation (4.11). In this way, the dynamics of infectious binding sites describe the dynamics of infectious individuals and the population of such individuals.
In particular, we find that F + is defined in terms of infectious (and susceptible) binding sites as follows: y e 0ȳ n−1 + (n − 1)ȳ e y 0ȳ n−2 (a | t − a, a + )da + da.
Since y e and y are probability vectors, they sum to one:ȳ e (a|t b , a + ) = 1 = y(a|t b , a + ). Moreover, with ϕ 1 given by (4.3), since we can express F + in terms of the history of x: (4.18) We can use the consistency condition for the total fraction of free binding sites: (4.19) to express F * in terms of the history of x by using (4.13) and (4.18). Thus we have specified all environmental variables for (4.7) in terms of (the history of) x. For completeness we briefly consider recovered binding sites.
Suppose a recovered binding site was born at time t b , infected at age a + , and recovered at age a * (and as usual, suppose its owner does not die in the period under consideration).We consider probabilities z e i (a | t b , a + , a * ) and z i (a | t b , a + , a * ) for recovered exceptional and non-exceptional binding sites in state i, respectively. The y and y e probabilities yield the conditions for z and z e at age a = a * , i.e.
The dynamics for z and z e can be described by a system of ODE similar to the ODE systems (4.14) for y and y e . Only now the mean field at distance one quantity Λ + needs to be replaced by Λ * where Λ * is defined in terms of p-level fractions by (2.21). By combining (2.21) with (4.11) and (4.10) we can express Λ * in terms of x-probabilities: Let ψ ( * ,k) (a | t b , a + , a * ) denote the probability that a recovered individual is in state ( * , k) given that it was born at time t b , infected at age a + and recovered at age a * , and does not die in the period under consideration. Then ψ ( * ,k) (a | t b , a + , a * ) can be expressed in terms of z and z e by replacing φ by ψ, y i by z i , and y e i by z e i in (4.15). The probability p ( * ,k) (a | t b ) is then obtained by taking into account all possibilities for age of infection a + and age of recovery a * : Finally, by relation (4.11), we obtain

A system of three renewal equations
To summarize, by replacing F * by (4.19), we are left with three environmental variables Λ − , F − , and F + which are defined by (4.22) Recall that x(a | t − a) is completely determined by ,t] , and Λ −| [t−a,t] , via (4.7)-(4.9). Therefore (4.20)-(4.22) is a closed system of three renewal equations.
Together, the three renewal Eqs. (4.20)-(4.22) fully determine the dynamics of ilevel probabilities p (−,k) (a | t b ) and p-level fractions P (−,k) (t). (Note that there are in total 1/6(n + 1)(n + 2)(n + 3) states of the form (−, k) One may not particularly like renewal equations to work with. However, the ODE system (4.7) has t b as a parameter, so is not finite dimensional. Therefore, contrary to Sect. 3, in order to describe the model with a closed finite system of ODE one needs to turn to p-level fractions P (−,k) and P (+,k) (the p-level system of ODE can be written down directly from the interpretation; see also (Leung et al. 2015) and Remarks 2 and 3). Together with the definition of the environmental variables F ± and Λ ± in terms of p-level fractions, the system is then closed. However, there are in total 1/3(n + 1)(n + 2)(n + 3) variables of the form P (±,k) .
As the system of three renewal Eqs. (4.20)-(4.22) has a clear interpretation, and R 0 , r , and the endemic steady state can very nicely be characterized from this system (see Sect. 4.3 below), we strongly advocate this formulation of the model rather than a (very high-dimensional) system with only ODE.

The beginning of an epidemic: R 0 and r
To describe the beginning of an epidemic, we are interested in characterizing R 0 and r . We have done so for the full p-level ODE system in Leung et al. (2015). In that paper, the characterization of R 0 involved the dynamics of infectious binding sites in the beginning of the epidemic. This infectious binding site system was then, via a linear map, coupled to the linearized p-system to show that the definition of R 0 via the interpretation indeed yields a threshold parameter with threshold value one for the p-level system.
In this section, we use the system of three renewal Eqs. (4.20)-(4.22) to characterize R 0 and r . Using the same arguments as in Sects. 2.3 and 3.3 of network cases I and II, we deduce that, in order to find a threshold parameter for the disease free steady state of the p-level system, we can focus on a threshold parameter for the stability of the disease free steady state of the binding site level system (4.7). Hence we can focus on (4.20)-(4.22).
The linearization of (4.20)-(4.22) involves the linearization of (4.7). The disease free steady state of (4.7) is given byx 0 where ϕ 0 (a), the probability that a binding site is free at age a given that it was born free (i.e. free at age 0), is given by (4.2).
We again put a ∧ on the symbols to denote the variables in the linearized system. The ODE for the linearized variablex 2 is straightforward: In the following we condition (as usual) on the owner of the binding site staying alive in the period under consideration. The probability y e 0 (a | t b , a + ) is independent of t b and given by (4.16). On the other hand, y 0 (a | t b , a + ) in the disease free steady state can be interpreted as the probability that a binding site is free at age a given that it is free at age a + with probability ϕ 0 (a + ). But this is equal to the probability ϕ 0 (a) that a binding site is free at age a given that it was born free at age 0 (since then, the probability that it is free at age a + is exactly ϕ 0 (a + )). So we find that, in the disease free steady state, where the first equality follows from simply evaluating (4.17) in the disease free steady state and the second can be deduced (as above) from the interpretation (or by algebraic manipulation). So we find thatF + satisfieŝ where we used relation (4.4) between F and ϕ 0 . We now derive two renewal equations forF + andΛ − . Variation of constants yields an expression forx 2 in terms ofF + andΛ − : (4.26) We substitute this in the expressions forF + andΛ − to find the system of two renewal equations: In preparation for defining and interpreting R 0 we write these integrals in convolution form:F (4.28) (the π 0 and π 1 appear by multiplying with F/F and (1− F)/(1− F)). This is a system of two renewal equations of the form (4.29) with non-negative kernelK . From these two renewal equations (4.27) and (4.28), we can obtain the characteristic equation and deduce threshold parameters r and R 0 . We define (4.30) Note that ∞ 0K (τ ) is a 2 × 2 matrix that can be evaluated explicitly so we have an explicit expression for R 0 . We define r to be the real root (if it exists) of the characteristic equation such that the spectral radius of ∞ 0 e −λτK (τ )dτ equals 1. Note that r is necessarily the rightmost solution of the characteristic Eq. (4.31).
Then r is a threshold parameter with threshold value zero for the stability of the disease free steady state of the system of renewal Eqs. (4.20)-(4.22). Furthermore sign(R 0 − 1) = sign(r ) so the definition (4.30) of R 0 indeed has the right threshold property. For R 0 > 1, to see that sign(R 0 − 1) = sign(r ), one uses that each matrix element of ∞ 0 e −λτ K (τ )dτ is a strictly monotonically decreasing function of λ and therefore the dominant eigenvalue of ∞ 0 e −λτ K (τ )dτ is strictly monotonically decreasing as a function of λ (Li and Schneider 2002;Diekmann et al. 2013, Sect. 8.2 the intrinsic growth rate). For R 0 < 1, one uses that the rightmost real solution r of (4.31) (if it exists) is strictly less than zero and this establishes the stability of the disease free steady state (Heijmans 1986;Inaba 1990;Thieme 2009).
In the epidemic context, 'reproduction' corresponds to transmission of the infectious agent to another host. The definition of (and the derivation of an expression for) R 0 in Leung et al. (2015) is in this spirit: it follows infectious binding sites in time and counts how many new infectious binding sites are formed when transmission occurs. A slight modification of the derivation in Leung et al. (2015) is required to generalize from SI to SIR. We did check that (4.30) is identical to the appropriately modified version of the dominant eigenvalue of (59) in Appendix C of (Leung et al. 2015).
Yet we would like to have a direct interpretation of the would-be reproduction number (4.30). To achieve this, it is helpful to think in terms of reproduction 'opportunities'. In the present context, these consist of +− links. In Leung et al. (2015) the spotlight is on the + side of the link. The present bookkeeping scheme focuses on x, so on − binding sites. So now the spotlight is on the − side of the link. The difference is just a matter perspective. A key point, however, is that after transmission the link disappears from the x stage. This forces us to formulate the interpretation in terms of reproduction opportunities rather than reproductions. (Note that, in traditional epidemiological models involving the random mixing assumption, contacts between individuals are instantaneous so there are no −+ links or 'reproduction opportunities' in the above sense.) We distinguish two birth-types of −+ links, according to the way they originate: Type 0: the −+ link was formed when a − binding site and a + binding site linked up Type 1: the −+ link is a transformed − link (one of the two owners got infected by one of its other partners) The relevant difference is the age distribution of the − binding site at the 'birth' of the −+ link (see Fig. 7): -for type 0 this distribution has density π 0 since the − binding site was free until that moment -for type 1 this distribution has density π 1 since the − binding site was (and remains) occupied v u w (-, k , k , k ) 1 2 3 Fig. 8 An example of a − − + configuration: u is one of the k 1 − partners of the 'middle' − individual v in state (−, k 1 , k 2 , k 3 ) and w is one of the k 2 + partners of v So the density of the age distribution of the − binding site at birth depends on the birth-type, making it necessary to distinguish between the two birth-types 0 and 1, while in case II there is no need to make this distinction.
In the nonlinear setting, the total rate in the population at which −+ links of type 0 are formed is equal to ρ F − k 0 P (+,k) = ρ F − n F + (note the −+ asymmetry here, which is in preparation for the linearization). The rate at which type 1 −+ links are formed is equal to β k 1 k 2 P (−,k) , respectively. Indeed, the expected number of free infectious binding sites in the population is k 0 P (+,k) , and the rate at which a free and infectious binding site acquires a susceptible partner is ρ F − . The expected number − − + configurations per 'middle' − individual is k 1 k 2 P (−,k) (see also Fig. 8) and the rate of transmission is β.
Rescaling (4.32) yields a system of renewal equations (4.34) (Again we note that each of these four integrals can be evaluated explicitly.) We now explain how (4.34) can be interpreted in terms of reproduction opportunities of types 0 and 1. A −+ link has no 'descendants' when transmission does not occur. When transmission occurs, it has at that very moment descendants of type 1, because the 'other' partners of the owner u of the − link then all of a sudden are connected to a + individual. In addition, it has descendants of type 0 when empty binding sites of u get occupied (necessarily by a − partner, since we consider the initial phase when + individuals are rare). Note that we should follow all binding sites of u until u either dies or becomes removed, since occupied binding sites may become free, occupied again, etcetera.
We now compute the expected number of descendants of either type for a −+ link given that the owner u of the − binding site has age a at the birth of the −+ link. The force of infection on u along the link equals β as long as -the + partner is alive and infectious -separation did not occur u is alive and not yet infected Hence the probability per unit of time that u is infected at age α + τ is given by βe −(μ+γ +σ +μ+β)τ .

The infectious binding site perspective
In Leung et al. (2015) the focus was on infectious binding sites. As exhibited by the densities π 0 and π 1 of the age-distribution of − individuals in the birth-types of −+ links, it matters whether a newly created −+ links is type 0 or type 1. To take this into account, in Leung et al. (2015), we kept track of the number of partners of susceptible partners of infectious binding sites. This led to the reduction to a 2 × 2 next-generation-matrix involving mean times spent with a susceptible partner with k partners, k = 1, . . . , n, (in the form of the inverse of an n × n matrix). We were able to find an explicit expression for R 0 although it required quite a lot of work to deal with this n × n inverse matrix.
The results in this paper teach us that, to take into account the birth-types of −+ links, we can also keep track of the age of susceptible partners rather than partners of partners. While age can be anything from zero to infinity, it can only move forward in time, i.e. individuals can only grow older. The same 2 × 2 NGM is obtained in a much more straightforward manner.
Whether doing the bookkeeping of partners of partners or of the age of partners, a big downside of taking the infectious binding site perspective is that it takes quite some work to prove that so-obtained R 0 is actually a threshold parameter for the stability of the disease free steady state of the p-level system (see (Leung et al. 2015)). This comes almost for free when taking the − perspective as we did in this paper.
Finally note that, whether we consider actual 'reproductions' (taking the + perspective) or 'reproduction opportunities' (taking the − perspective), both yield the exact same threshold parameter R 0 so in that sense it does not matter which perspective we take. However, while the the dominant eigenvalue R 0 of the next-generation-matrix is the same with both perspectives, the matrices themselves are different. And so are the underlying interpretations.

Endemic steady state
Let E = (F + , F * , Λ − ) be the vector of environmental variables. Note that we use consistency relation (4.19) to substitute environmental variable F * for F − . This choice of environmental variables leads to the disease free steady state corresponding to E = (0, 0, 0). Then we have a system of three renewal equations for E.
Let ,a] via (4.7)-(4.9). Therefore, is a closed system of three renewal equations. In endemic equilibrium, the environmental variable E is constant (note that, if E is constant, then also p-level fractions are constant and binding-site-and i-level probabilities are constant as functions of time of birth t b ). So the endemic steady state is characterized as a solution to the fixed point problem (4.35) where now the symbols denote the values of constant functions. The fixed point problem always has a trivial solution given by the disease free steady state E = (0, 0, 0). Note that a solution E to (4.35) needs to have biological meaning. Therefore, we only consider solutions that satisfy F + , F * ≥ 0, 0 ≤ F + + F * ≤ F, and 0 ≤ Λ − ≤ n − 1.
Conjecture: If R 0 < 1, then the only solution to the fixed point problem is the trivial solution. If R 0 > 1, then there is a unique nontrivial solution.
Open problem: Prove (or disprove) the conjecture.
In Appendix 2 we elaborate on an unsuccessful attempt at a proof of the conjecture for the simpler case of an SI infection, rather than an SIR infection, obtained by setting γ = 0. This attempt tried to use Krasnoselskii's method Krasnoselskii (1964) (see also Hethcote and Thieme 1985).
Note that the three-dimensional fixed point problem (4.35) provides a way to find the endemic steady state numerically. Furthermore, even though we did not manage to prove the conjecture, numerical investigations strongly suggest that all conditions for Krasnoselskii's method are satisfied and that the conjecture holds true.

Conclusions and discussion
In this paper we formulated binding site models for the spread of infection on networks. The binding sites serve as building blocks for individuals. In fact we considered three different levels: (1) binding sites, (2) individuals, and (3) the population. On both the binding site and individual level, we have a Markov chain description of the dynamics, where feedback from the population is captured by environmental variables. These environmental variables are population-level quantities. By lifting the individual level to the population level (where the model is deterministic), the feedback loop can be closed. In the end, this leads to a model description in terms of susceptible binding sites in case I and in terms of just environmental variables in cases II and III.
The systematic model formulation leads, in all three cases, to only a few equations that determine the binding site, individual, and the population dynamics. Moreover, from these equations we derive the epidemiological quantities of interest, i.e. R 0 , r , the final size (in cases I and II) and the endemic steady state (in case III).
Quite in general, understanding is enhanced by an elaboration of the interpretation of R 0 in a specific context. In cases I and II we have taken the obvious perspective of a + binding site to do so. But in case III, cf. Sect. 4.3, we reasoned in terms of 'reproduction opportunities'. These consist of +− links. From these links we took the − perspective. Somewhat surprisingly, this turned out to lead quickly and efficiently to a simple interpretation. Moreover, the derivation of R 0 follows from the system of equations in a natural manner. One can adopt the − perspective in cases I and II too, but there it does not change much. Yet we wouldn't be surprised if the − perspective turns out to be powerful in other dynamic network models of infectious disease transmission.
Several open problems remain. Although we are able to implicitly characterize the final size in case II, we have not been able to make it more explicit. We would like a characterization in the same spirit as (2.33) for case I, but we have not succeeded and our optimism subsided. A more useful characterization of the endemic steady state was given for case III as a three-dimensional fixed point problem. Unfortunately, we have not (yet) been able to prove the existence and uniqueness of a nontrivial fixed point for R 0 > 1 (and that no such fixed point exists for R 0 < 1) and therefore we posed this as a conjecture in Sect. 4.4.
Of another nature are open problems related to the mean field at distance one assumption. While, in case I, the mean field at distance one assumption is proven to be exact in the appropriate large population limit of a stochastic SIR epidemic on a configuration network (Decreusefond et al. 2012;Barbour and Reinert 2013;Janson et al. 2014), it remains an open problem whether or not this also holds for the dynamic network case II (we conjecture it does). In the dynamic network case III, we know that the mean field at distance one assumption is really an approximation of the true dynamics as we pointed out in the introduction of this paper (see also Leung et al. 2015). What we have not discussed is how good or bad of an approximation it is. In particular, are there conditions for which the approximation works nicely and can we understand intuitively the extent to which this assumption violates the truth?
In both cases II and III, we ended the model formulation with renewal equations. In case II one can just as easily consider a system of ODE, and we represented this view also in Sect. 3.2.3. In case III, a system of ODE clearly becomes inconvenient. An ODE formulation in that case would require at least 1/3(n + 1)(n + 2)(n + 3) variables, while, by considering a system of renewal equations, only three equations are needed. More importantly, the system of renewal equations has the huge advantage that R 0 and r more or less immediately follow from the linearization of the system in the disease free steady state. The calculations are straightforward, the expressions are interpretable biologically, and the proof that R 0 and r are threshold parameters for the disease free steady state of the p-level system comes more or less for free.
By distinguishing the three different levels, and formulating the model on the binding site level, one can consider several generalizations (see also Leung et al. 2012, 2015 for a discussion). In principle, any generalization that maintains the (conditional) independence assumption for binding sites of an individual fits within this framework. One can think of generalizations concerning the network or generalizations concerning the infectious disease. For the infectious disease, one can easily take any compartmental model such as SIR, SEIR, SI, SI 1 I 2 (as long as infected individuals can not return to the susceptible class within their lifetime). The main difference is in the different states that a binding site can be in. Generalizations of the network that one can think of are (i) a heterosexual population rather than a homosexual population, (ii) allowing for different n in the population, i.e. letting n be a random variable (which we already considered in the static network case in Sect. 2.5) (iii) allowing for multiple types of binding sites, e.g. binding sites for casual and steady partnerships, and combinations of the three. One can formulate models incorporating these generalizations by following the five steps described in Sect. 2.2. The main added difficulty is in the bookkeeping that becomes more involved. But in terms of the characterization of R 0 and the endemic steady state, mathematically speaking the situation does not become more complex.
Generalization (ii) is in a sense different from the other generalizations. At first sight this generalization seems rather straightforward: one only needs to average over n in the right way. But things are more subtle than they seem. One needs to also take into account the partnership capacity of partners in the bookkeeping. In particular, the mean field at distance one assumption is expressed in environmental variables Λ − (n) that can be interpreted as the expected number of infectious partners of a susceptible partner with partnership capacity n of a susceptible individual. However tempting it is, we cannot simply assume that we can average over all possible partnership capacities of susceptible partners, i.e. we cannot average over these terms Λ − (n) and consider only one termΛ − (at least in the static network case such an assumption would yield a system different from the Volz equations Volz 2008). Our modelling framework is powerful in the way that such subtleties are exposed. In particular, the systematic approach that connects the binding site level, i-level, and p-level yields better insight into the ramifications that generalizations have at the different levels.
Finally, in the current framework, and as usual in literature, demographic turnover as considered in case III takes the individual's age to be exponentially distributed. This assumption is mainly for mathematical convenience and is not realistic for many populations. We believe that it is possible to relax the assumption on the age distribution to consider more general survival functions. In that case, lifting the i-level to the p-level changes, and one needs to take into account the age of partners (but hopefully this may be done by simply averaging in the right way). Moreover, in the current framework, disease does not impact mortality. In the context of HIV, disease-related mortality is certainly very relevant. We believe that the framework presented in this paper provides a way to incorporate this by means of the infectious y binding sites. While these generalizations relating to the demographic process are less straightforward to implement than the ones described in the previous paragraph, the current framework provides an excellent starting point.
Linearization of an epidemic system in the disease free steady state leads to a linear system that leaves a cone, characterized by positivity, invariant. Perron-Frobenius theory, or its infinite dimensional Krein-Rutman variant, yields the existence of a simple eigenvalue r such that (i) the corresponding eigenvector is positive (ii) Re λ < r for all eigenvalues λ = r The theory of stable and unstable manifolds yields a nonlinear analogue: the nonlinear system has exactly one orbit that is tangent to the eigenvector corresponding to eigenvalue r . If r > 0 then this orbit belongs to the unstable manifold and tends to the disease free steady state for t → −∞. If r < 0 then the orbit belongs to the stable manifold and tends to the disease free steady state for t → +∞. Our interest is in the case r > 0.
Note that one orbit of an autonomous dynamical system corresponds to a family of solutions that are translates of each other. See (Diekmann 1977) for an early example of this type of result (but note that the proof in that paper has a flaw; see (Diekmann and van Gils 1984, Sect. 7) for a flawless proof).
These ideas apply directly to the three-dimensional ODE system (2.5) in case I. For the scalar renewal Eq. (2.30) we can refer to Sect. 1 of (Diekmann and van Gils 1984) provided that we are willing to assume that F has compact support. For the ODE system of case II there exists an eigenvalue zero (corresponding to conservation of binding sites). This eigenvalue zero creates havoc. Presumably, the difficulties can be overcome by the introduction of a tailor-made cone, but we did not elaborate this in all required detail. The alternative is to consider the scalar renewal Eq. (3.13) for F + and to combine ideas from Diekmann et al. (2007) with theory developed in Diekmann and Gyllenberg (2012). This combination should, we think, also cover the system of renewal Eqs. (4.20)-(4.22) for case III.

Appendix 2: Endemic steady state: unsuccessful attempt at a proof
We explain our attempt to prove the conjecture of Sect. 4.4 about the existence and uniqueness of solutions to the fixed point problem (4.35) for the simpler case of an SI infection rather than an SIR infection (set γ = 0). We only need to consider two environmental variables, rather than three, as we will explain. This attempt to prove the conjecture uses the sublinearity method of Krasnoselskii (1964) (see also Hethcote and Thieme 1985), the idea of which for one dimension is represented in Fig. 11.
(5.1) with boundary condition (5.2) As before in Sect. 4, x(a | t b ) is completely determined by via (5.1)-(5.2). In particular, there are now only two environmental variables F + and Λ − . These environmental variables satisfy renewal equations. Let then we obtain a fixed point problem for the environmental variables F + and Λ − : Note that in endemic equilibrium the environment is constant, i.e. F + (t) =F + , Λ − (t) =Λ − . Therefore x no longer depends on time of birth. In what follows we write x = x(a).
The fixed point problem (5.3) can be related to R 0 by considering the linearizaton of the right hand side of (5.3) in the disease free steady state (F + , Λ − ) = (0, 0). Indeed, the linearization DG(0, 0) has dominant eigenvalue R 0 .
Monotonicity and sublinearity of G 1 in both variables F + and Λ − is easily proven. One can show that the derivatives of x 0 , x 1 , and x 1 + x 2 with respect to F + and Λ − are nonpositive while the mixed second order derivatives are all nonnegative. Then one can easily prove that the derivatives D i G 1 (F + , Λ − ) ≥ 0 showing that G 1 is a monotonically increasing function of both F + and Λ − . Sublinearity can be proven by showing that the function f (t) = G 1 (t (F + , Λ − ))−tG 1 (F + , Λ − ) satisfies f (t) < 0 for 0 < t < 1.
We work out only the proof to show that D 1 G 1 (F + , Λ − ) ≥ 0. The derivative of G 1 with respect to F + is equal to n−2 (a)da.

Remark 4
The variable x 2 is not necessarily monotone in F + or Λ − . One can find parameter values for which we find that ∂ x 2 /∂ E(a), E = F + , Λ − , is neither nonpositive nor nonnegative as a function of a.
Note that the feedback function G 2 for Λ − involves x 2 . The arguments to prove monotonicity and sublinearity do not seem to work for G 2 . Numerical investigation strongly suggest that G 2 is indeed monotonically increasing as a function of both F + and Λ − as well as sublinear. So far, we have not been able to provide a proof.
Nevertheless, once we show that both G 1 and G 2 are monotonically increasing functions of environmental variables F + and Λ − and sublinear, Krasnoselskii's method then provides a proof that for R 0 < 1 only the trivial solution exists and for R 0 > 1 there exists a unique nontrivial solution to (5.3).