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 [Formula: see text], 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 listingIn 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.

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

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, 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. 1998b, 2001). This means that we start from the notion of state at 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.

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

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, 2012, 2013) 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 p-level equations, but used the i-level version to derive an expression for . 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: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.

The disease dynamics are fast relative to any partnership- or demographic changes. The network is static and everyone is susceptible at time .

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 .

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.

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: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.

0 - free

1 - occupied by a susceptible partner

2 - occupied by an infectious partner

3 - occupied by a recovered partner.

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 (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 , . Here 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 binding sites (in other words, by conditioning on the individual not getting infected through its 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 toIn 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 . We call the binding site through which the change in the owner’s disease status occurred the ‘exceptional’ binding site. The other binding sites are i.i.d. and, at time , they are distributed according to . 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.

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 . 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

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 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, 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 . 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, 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

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 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 at which a susceptible partner of a susceptible binding site becomes infected at time t. Here 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’ conditionsandTo 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 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 time t is with(recall (1.1) and note that in case I we have ). The probability that a binding site is in state 2, given that the binding site is susceptible, is . Therefore,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: , and .

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).

Next, we introduce , denoting the fraction of the population with label . Here denotes the number of partners of an individual with each of the different disease statuses, i.e. susceptible, infectious, and recovered partners. Furthermore, 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 .

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:Here the sums are over all possible configurations of and with , . The second factor in each term of this sum denotes the probability that a susceptible partner of a susceptible individual is in state . The number of infectious partners is then given by , and we find the expected number of infectious partners by summing over all possibilities.

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

In the third step, we let denote the probability that an individual is in state 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 susceptible, infectious, and recovered partners is given by the multinomial expressionThe 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 .

Finally in the last step, by combining steps 2, 3, and 4, we can express in terms of the x-probabilities.

Step 5. Environmental variables in terms of x-probabilities (combining 2, 3, 4)

By combining (2.8), and (2.7) we find that and . Then definition (2.6) yields the same expression for as (2.4).

Finally, steps 1 to 5 together yield the closed system (2.5) of ODE for x. The dynamics of the i-level probabilities are fully determined by the system of three ODE for x. We can use this three-dimensional system of ODE to determine r, , 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 for . 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.

Remark 2

One obtains the p-level ODE system by differentiation of (2.7) and use of (2.5) and (2.8). In doing so, one obtains a system of ODE for the p-level fractions concerning individuals with a − disease status: , with defined by (2.6) (compare with Lindquist et al. 2011, Eq. (13)).

The beginning and end of an epidemic: , , and final epidemic size

In this section we consider the beginning and end of an epidemic. We first focus on and r, so on the start of an epidemic.

Note that we can very easily find an expression for from the interpretation: when infected individuals are rare, a newly infected individual has exactly 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 isWe now rederive and derive r from the binding site system (2.5).

Note that the p-level fractions can be fully expressed in terms of the binding site level probabilities (Eqs. (2.8), (2.7)). Furthermore, the 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 state , , yields a decoupled ODE for the linearization of the ODE for . To avoid any confusion, let denote the linearized variable. Then the linearization yieldswith (for ) ‘far past’ behaviour for . In particular, the right-hand side of the ODE for depends only on .

To illustrate the method used in case II and III in Sects. 3.3 and 4.3, we derive expressions for and r from a special form of the characteristic equation. Variation of constants for the ODE of yieldsSubstituting the ansatz yields the characteristic equationThen 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 . Likewise, we can derive the expression (2.9) for by evaluating the integral with .

Next, we consider the final epidemic size. We do so by considering the dynamics of defined in (2.3). Recall (1.2), i.e. the probability that an individual is susceptible at time t, is given by . We observe that, by (2.8), is also equal to the fraction of susceptible individuals in the population at time t (Alternatively, one can show that by combining (2.8) and (2.7)). In fact, it is possible to describe the dynamics of in terms of only 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 relationWe can use the interpretation to derive (2.10); is the probability that a susceptible binding site with owner u is occupied by a susceptible partner v, 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 for :The fraction of susceptible individuals at the end of the outbreak is determined by the probability . Since satisfies (2.11) and is a constant, we find that necessarily is the unique solution in (0, 1) ofif . The final epidemic size is given byIn Sect. 2.5 we show that one can actually describe the dynamics of the probability 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 for . The final size equation is then obtained by simply taking the limit . We highly recommend reading Sect. 2.5 to understand the derivation of the renewal equation for 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 (and likewise p-level fractions given by (2.8)) for .

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 and that it does not recover in the period under consideration. Let denote the probability for the exceptional binding site to be in state i at time t, . Similarly, denotes the probability for a non-exceptional binding site to be in state i at time t, . Here the probabilities are defined only for . Note that y and 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 . 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 with probability 1, while the distribution of the state of a non-exceptional binding site at time is given by , i.e. we have boundary conditionsThe 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 for susceptible individuals: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 ).

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)

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 are described by the same ODE systemwithand case specific boundary conditions (2.13). Observe that this means that for all .This also immediately follows from the interpretation: at time , 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 - 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

We let denote the probability that an infectious individual, infected at time , is in state 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)Note that for , i.e. for all at least one partner is not susceptible.

A susceptible individual becomes infected at time 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 , i.e. the rate at which a susceptible individual becomes infected at time , equals and consequently the incidence at time , i.e. the fraction of the population that becomes, per unit of time, infected at time , equals(recall that is the fraction of the population that is susceptible).

Furthermore, an infectious individual that is infected at time is still infectious at time t if it does not recover in the period . Since the infectious period of an individual is assumed to be exponentially distributed with rate , the probability that this happens isWe then find an expression for the unconditional i-level probabilities that a randomly chosen individual is in state at time t in terms of infectious binding site probabilities and the history of susceptible binding site probabilities:where is given by (2.17). The i-level probabilities 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 explicitly depend on the dynamics of x through the boundary conditions (2.13) and the environmental variable  (2.15). In addition, 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 ODE for the p-level fractions concerning individuals with a disease status: , with defined by (2.6) (compare with Lindquist et al. 2011, Eq. (13)).

In case of recovered individuals, one considers their binding sites and first conditions on time of infection and time of recovery . 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 and for non-exceptional and exceptional binding sites, i.e.The dynamics for z and can be described by a system of ODE identical to the ODE systems for y and , but with replaced by . The environmental variable is given byBy combining (2.21) with (2.8) and (2.7) we findWe find an expression for the probability that a recovered individual, infected at time and recovered at time , is in state at time , in terms of z and probabilities for recovered binding sites with the same reasoning as for (one can simply replace by , by , and by in (2.17)). Then, to arrive at an expression for the unconditional probability , we again need to take into account the incidence at . The probability that recovery does not occur in the time interval is given by and the rate at which an infectious individual recovers is , thereforewhere 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 . 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 towith 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 introduceIn the SIR example we haveIt 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, equalsFor the Markovian SIR example (2.24) boils down toa 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 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 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 ) 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).

Fig. 5

Volz focused on the variable 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

The variable corresponds to introduced in Sect. 2.1 and therefore we use that symbol also in this section. We reformulate (2.3) as(see also Fig. 6). There is an underlying stochastic process in the definition for 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 . 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.

Fig. 6

Schematic representation of . 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

To derive the consistency relation for 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 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 other binding sites are susceptible. Conditioning on its other binding sites not transmitting to v, a binding site of v is susceptible at time t with probability . Therefore, given susceptibility of partner u, v is susceptible at time t with probabilityThis just repeats the consistency relation (2.10) stating that the probability that a susceptible binding site is occupied by a susceptible partner is equal to the probability that a partner of a susceptible individual is susceptible.

Next, suppose that v gets infected at some time , then u is not infected by v before time t if no transmission occurs in the time interval of length . The expression (2.27) has as a corollary that the probability per unit of time that v becomes infected at time equalsNoting that the probability of no transmission to u in the time interval is we conclude that necessarily,Finally, by integration by parts, we obtain the renewal equationFor a configuration network with general degree distribution 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))withThe solution , , of (2.30) corresponds to the disease free situation. If we put and assume h is small, we easily deduce that the linearized equation is given byThe corresponding Euler-Lotka characteristic equation readsIf we evaluate the right hand side of (2.32) at , we obtaincf. 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)).

To derive an equation for the final size is even simpler, one takes the limit in (2.30) to deduceand next observes that the escape probability is given by(to compare to Diekmann et al. 2013, Eqs. (12.36)–(12.38), p. 295) identify , , and rewrite (2.33) as ).

In the case that is given by (2.25), the REcan be transformed into an ODE for by differentiation:In the special case of Sects. 2.1–2.4, we have and for all so 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 is of the formwhere, for some , and V are non-negative vectors in while is a positive-off-diagonal (POD) 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 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, , etcetera.) If is given by (2.34), the variablesatisfies the ODEand, since (2.30) can be rewritten asthe Eq. (2.35) is a closed system once we replace 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 about . 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 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 . 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 , while an occupied binding site becomes free at rate . Similar to (Leung et al. 2012) (set ), we find that F satisfies the ODESo we find that F converges to a constant for . Therefore, we assume that the fraction of free binding sites is constant, and this constant is again denoted by the symbol F. Then F satisfiesAlthough 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 . 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 ).

Finally, later on in Sect. 3.2.2, we also need the probability 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 is the unique solution of the initial value problem:we havewhere we used (3.1) in the second equality.

Susceptibles

We describe the dynamics of susceptible binding sites in terms of x-probabilities. Consider a susceptible binding site and suppose its owner does not become infected through one of its other binding sites for the period under consideration. An occupied binding site (in states 1, 2, or 3) becomes free if it loses its partner (corresponding to a transition to state 0). This occurs at rate . A binding site that is free, i.e. a binding site in state 0, can acquire a partner. The rate at which this occurs is where F is the fraction of free binding sites defined by (3.1). Free binding sites either have a susceptible, infectious, or recovered owner. So there are three additional environmental variables, viz. the fraction of binding sites that are free and have disease status d (i.e. having owners with disease status d), we denote these by , . Then . Finally, there are infection and recovery events that can cause state transitions (as in the case of a static network in Sect. 2).

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 , i.e. we have ‘far past’ conditionsLet . The environmental variables 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), andNext, 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 with .

In step 3 we define the i-level probabilities in terms of the x probabilities by using the conditional independence of binding sites: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 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 thatexactly as the interpretations of and x would suggest.

Before we can specify and in terms of (the history of) x we need to define p-level fractions and . 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 , the binding site through which infection is transmitted is from that point on ‘exceptional’. Then, given that its owner became infected at time and that it does not recover for the time under consideration, we consider an infectious binding site. Let denote the probability that the exceptional binding site is in state i at time t and 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 . The exceptional binding site is with probability one in state 2 at time . The states of all other binding sites are distributed according to . Therefore, we put boundary conditionsSince the individual does not recover in the period under consideration, the infectious binding sites behave independently of one another.

The dynamics of and y are both governed by the systemwithNote 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 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 for .

Next, similarly to Sect. 2.4, by combinatorics (but now probabilities and are not equal to zero for ), we find that the probability that an individual, infected at time , is in state at time is given byThe probability that a randomly chosen individual is in state at time t is obtained by taking into account the time of infection and the probability (2.19) that an individual has not recovered time after infection. The definition (2.18) for the incidence carries over (but now with the x defined by the ODE system (3.4)). SoBy combining (2.8) and the expression (3.11) for in terms of y and we can redefine 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 in terms of y and :where and . Since and y are probability vectors,Next, we consider the probabilities . Note thatwith given by (3.2). The dynamics of are described in terms of and the history of (by means of the boundary condition). We can solve for . This yields(note that time of infection matters in this probability and not only the length of the time interval). We can further simplify (3.12) to obtainwhich only depends on the model parameters and past probabilities for susceptible binding sites.

We can use the consistency condition for the total fraction of free binding sites:to express in terms of the history of x (use (3.8) for and (3.13) for ). 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 , recovers at time . After recovery, we still distinguish between the exceptional binding site and the other binding sites. We introduce probabilities and for recovered binding sites. The y and probabilities yield the conditions for z and at time , i.e.The dynamics of z and are described by the system of ODE for y and , 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 denote the probability that a recovered individual is in state given that it was infected at time and recovered at time . Then can be expressed in terms of z and by replacing in (3.10) by , by , and by .

The unconditional probability is then obtained by taking into account time of infection and recovery time :which, by relation (2.8), is equal to the p-level fraction . Note that we can also use this definition of to define in terms of x similar to the way we did for 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 in terms of x (Eqs. (2.4), (3.8)). Subsequently, in Sect. 3.2.2, by considering infectious binding site probabilities y, and , we also defined and 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 given by (3.13) and with ‘far past’ conditionThe ODE (3.16) for x together with the expression (3.13) for 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 i-level probabilities (hence p-level fractions) 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 . Interpreting , , and 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 .

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 . In order to do so, we first consider the fraction 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 . Next, we consider . 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 at which free binding sites become occupied. The fraction of infectious binding sites that are occupied is given by and the rate at which these binding sites become free is . Then, finally, a susceptible individual with infectious partners becomes infected at rate , taking into account all we find probability per unit of time at which a susceptible individual becomes infected. The probability that a non-exceptional binding site is free and susceptible upon infection is , so the expected fraction of free binding sites created upon infection of a susceptible individual is . Hence there is a flow into . We have the following ODE for :with ‘far past’ condition .

Alternatively, we can derive the ODE (3.19) for 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 (similar to in Sect. 3.2.2). This yields .

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 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 that prepares for the characterization of r and 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 and final size

First, just as in case I, the final size is given byBut while in case I we derived a simple scalar equation for ((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 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 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 by , , . Linearization in this state yields a decoupled system of equations for the linearized and equations. We let and denote the variables in the linearization in the disease free steady state. Note that, in the disease free steady state , i.e. 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 is equal to the probability F that a randomly chosen binding site is free. Thenwhich can be viewed as a linear delay differential equation for . In order to obtain an informative version of the corresponding characteristic equation, we rewrite it as a renewal equation for .

Variation of constants for the ODE for yields:Substituting (3.20) into this expression yields the renewal equation for :with(where the rearrangement of the terms in the integrals is in preparation for the interpretation). Next, we substitute the ansatz , and obtain the characteristic equationThere is a unique real root to (3.21) and this root is by definition the Malthusian parameter r. We define . Then sign(r)=sign(), and we find that is a threshold parameter with threshold value one for the stability of the disease free steady state, with equal toWe can evaluate the integrals and find an explicit expression for . 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 . 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 . 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 (note that, in addition to the susceptible partners, u has 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 . The exceptional binding site of u is free at time with probability (see (3.2)). Each of the 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 (at the beginning of the epidemic). Integrating over all possible lengths of the infectious period, we find that 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 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).

In a world with demographic turnover, next to calendar time, also age matters. We keep track of both age a and time of birth of an individual (calendar time is then given by ). 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 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 denote the probability that a binding site is free at age , given that it was free at age , and denotes the probability that a binding site is free at age , 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 as a function of a is described bywith initial conditions, respectively,The explicit expressions for the are given by See also (Leung et al. 2012, Eq. (10)) (where can be identified with ) and (Leung et al. 2015, Eq. (67)) (where and can be identified with and , 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 bySimilarly, 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 in Sect. 4.3 below, we remark that and 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 and concern the age; the probability distributions, however, provide the same information).

Demography does not give rise to any additional environmental variables, we still deal with the mean field at distance one variable , and the fractions of free binding sites with disease status 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 , and suppose that its owner, for the period under consideration, does not die and does not become infected through one of its other binding sites. Let . The dynamics of x as a function of age are described by the following system of equations:with given byAn 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 conditionsGiven the environmental variables and , we can formally view as a function of the environmental variables:i.e. is completely determined byWe 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 in terms of x. Given that an individual was born at time and does not die in the period under consideration, the probability that an individual is at age a in state 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 to i-level probabilities . In order to do so, we use the stationary age distribution with density . The fraction of the population that is in state 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 in terms of x. This can be done by combining (4.10) and (4.11) with (2.6) (for ) or (3.6) (for ). We find thatandIn order to complete step 5 (expressing the environmental variables and in terms of x) we need to consider infectious and recovered binding sites.

Consider a binding site that was born at time and infected at age and its owner remains alive and infectious for the period under consideration. Note that age for this individual corresponds to calendar time . Let denote the probability that the exceptional binding site is in state i at age a and the same probability for a non-exceptional binding site.

Then, at age , the exceptional binding site is for certain in state 2, while the other binding site states are distributed according to :The dynamics of infectious binding sites are described by:withAgain, there is no rate in 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 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 . The probability that an individual, born at time and infected at age , is in state at age is given byThe contribution to the incidence of individuals of age , born at time and alive for the period under consideration, is given bywhere the reasoning is similar to cases I and II. Then, taking into account all possible ages of infection , and the probability that as yet recovery did not occur, the probability that an individual, born at time , is in state at age a is given byThe p-level fractions 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 is defined in terms of infectious (and susceptible) binding sites as follows:Since and y are probability vectors, they sum to one: . Moreover, with given by (4.3), since we can express in terms of the history of x:We can use the consistency condition for the total fraction of free binding sites:to express 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 , infected at age , and recovered at age (and as usual, suppose its owner does not die in the period under consideration).We consider probabilities and for recovered exceptional and non-exceptional binding sites in state i, respectively. The y and probabilities yield the conditions for z and at age , i.e.The dynamics for z and can be described by a system of ODE similar to the ODE systems (4.14) for y and . 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 denote the probability that a recovered individual is in state given that it was born at time , infected at age and recovered at age , and does not die in the period under consideration. Then can be expressed in terms of z and by replacing by , by , and by in (4.15).

The probability is then obtained by taking into account all possibilities for age of infection and age of recovery :Finally, by relation (4.11), we obtain

A system of three renewal equations

To summarize, by replacing by (4.19), we are left with three environmental variables , , and which are defined by Recall that is completely determined byvia (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 i-level probabilities and p-level fractions . (Note that there are in total states of the form , with , .)

One may not particularly like renewal equations to work with. However, the ODE system (4.7) has 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 and (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 and in terms of p-level fractions, the system is then closed. However, there are in total variables of the form .

As the system of three renewal Eqs. (4.20)–(4.22) has a clear interpretation, and , 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: and r

To describe the beginning of an epidemic, we are interested in characterizing and r. We have done so for the full p-level ODE system in Leung et al. (2015). In that paper, the characterization of 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 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 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 by , , , where , 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 variable 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 is independent of and given by (4.16). On the other hand, 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 with probability . But this is equal to the probability 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 is exactly ). 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 that satisfiesNext, linearization of yieldswhere we used relation (4.4) between F and .

We now derive two renewal equations for and . Variation of constants yields an expression for in terms of and :We substitute this in the expressions for and to find the system of two renewal equations:In preparation for defining and interpreting we write these integrals in convolution form: (the and appear by multiplying with F / F and ). This is a system of two renewal equations of the formwith non-negative kernel .

From these two renewal equations (4.27) and (4.28), we can obtain the characteristic equation and deduce threshold parameters r and . We defineNote that is a matrix that can be evaluated explicitly so we have an explicit expression for . We define r to be the real root (if it exists) of the characteristic equation such that the spectral radius of 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( sign(r) so the definition (4.30) of indeed has the right threshold property.

For , to see that sign( sign(r), one uses that each matrix element of is a strictly monotonically decreasing function of and therefore the dominant eigenvalue of is strictly monotonically decreasing as a function of  (Li and Schneider 2002; Diekmann et al. 2013, Sect. 8.2 the intrinsic growth rate). For , 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) 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):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.

Fig. 7

The two birth-types of links between individuals u and v

for type 0 this distribution has density since the − binding site was free until that moment

for type 1 this distribution has density since the − binding site was (and remains) occupied

In the nonlinear setting, the total rate in the population at which links of type 0 are formed is equal to (note the asymmetry here, which is in preparation for the linearization). The rate at which type 1 links are formed is equal to , respectively. Indeed, the expected number of free infectious binding sites in the population is , and the rate at which a free and infectious binding site acquires a susceptible partner is . The expected number configurations per ‘middle’ − individual is (see also Fig. 8) and the rate of transmission is .

Fig. 8

An example of a configuration: u is one of the − partners of the ‘middle’ − individual v in state and w is one of the partners of v

Linearization in the disease free steady state yield and , respectively (use p-level definition (2.6) for ). These observations motivate us to scale the two renewal Eqs. (4.27) and (4.28). LetThat such a rescaling does not affect the definition of r and follows from the following observation:

Observation

In general, if we have a system of renewal equations of the form (4.29), we may ‘scale’ , i.e. put and consider the renewal equationwith and C the diagonal matrix with non-zero entries . ThenMoreover, the matrices and are similar, so they have the same eigenvalues. In particular, they have the same dominant eigenvalue .

Rescaling (4.32) yields a system of renewal equationswith , and a matrix with matrix elements(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 asHence the probability per unit of time that u is infected at age is given by .

the partner is alive and infectious

separation did not occur

u is alive and not yet infected

When u is infected at age an expected number of offspring of type 1 is produced. A schematic representation is given in Fig. 9. This is how and in (4.34) can be interpreted.

Fig. 9

The production of type 1 offspring

Offspring of type 0 are (potentially) produced by both free and occupied (at the time of infection) binding sites of u. To calculate the mean number of offspring of type 0, suppose u is infected at age . Then u is alive and infectious at age with probability . The expected number of free binding sites it has at age is equal to . A free binding site becomes occupied at rate . Integrating over all possible , we find that the expected offspring of type 0 isThis is how and in (4.34) can be interpreted. A schematic representation is given in Fig. 10.

Fig. 10

The production of type 0 offspring

The infectious binding site perspective

In Leung et al. (2015) the focus was on infectious binding sites. As exhibited by the densities and 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 next-generation-matrix involving mean times spent with a susceptible partner with k partners, , (in the form of the inverse of an matrix). We were able to find an explicit expression for although it required quite a lot of work to deal with this 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 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 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 so in that sense it does not matter which perspective we take. However, while the the dominant eigenvalue 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 be the vector of environmental variables. Note that we use consistency relation (4.19) to substitute environmental variable for . This choice of environmental variables leads to the disease free steady state corresponding to . Then we have a system of three renewal equations for E.

Letwhere is completely determined by 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 ). 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 . Note that a solution E to (4.35) needs to have biological meaning. Therefore, we only consider solutions that satisfy , , and .

Conjecture: If , then the only solution to the fixed point problem is the trivial solution. If , 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 . 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, 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 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 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 (and that no such fixed point exists for ) 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 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 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 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, (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 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 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 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.