Abstract
Background
Human influenza is characterized by seasonal epidemics, caused by rapid viral adaptation to population immunity. Vaccination against influenza must be updated annually, following surveillance of newly appearing viral strains. During an influenza season, several strains may be cocirculating, which will influence their individual evolution; furthermore, selective forces acting on the strains will be mediated by the transmission dynamics in the population. Clearly, viral evolution and public health policy are strongly interconnected. Understanding populationlevel dynamics of coexisting viral influenza infections, would be of great benefit in designing vaccination strategies.
Methods
We use a Markov network to extend a previous homogeneous model of two cocirculating influenza viral strains by including vaccination (either prior to or during an outbreak), age structure, and heterogeneity of the contact network. We explore the effects of changes in vaccination rate, crossimmunity, and delay in appearance of the second strain, on the size and timing of infection peaks, attack rates, and diseaseinduced mortality rate; and compare the outcomes of the network and corresponding homogeneous models.
Results
Prevaccination is more effective than vaccination during an outbreak, resulting in lower attack rates for the first strain but higher attack rates for the second strain, until a “threshold” vaccination level of ~3040% is reached, after which attack rates due to both strains sharply dropped. A small increase in mortality was found for increasing prevaccination coverage below about 40%, due to increasing numbers of strain 2 infections. The amount of crossimmunity present determines whether a second wave of infection will occur. Some significant differences were found between the homogeneous and network models, including timing and height of peak infection(s).
Conclusions
Contact and age structure significantly influence the propagation of disease in the population. The present model explores only qualitative behaviour, based on parameters derived for homogeneous influenza models, but may be used for realistic populations through statistical estimates of interage contact patterns. This could have significant implications for vaccination strategies in realistic models of populations in which more than one strain is circulating.
Background
Human influenza infection is characterized by seasonal epidemics. This occurs because influenza A is able to maintain its presence in human populations by evolutionary adaptations to populationwide immunity, resulting in mutations that gradually change viral antigens allowing the virus to evade immune detection, a process known as “antigenic drift”. On account of these rapid mutations, vaccination for influenza must be updated annually on a global basis, following surveillance to monitor the appearance of new strains [1]. Antigenic drift also diminishes vaccine efficacy for mutant strains, but may still confer partial immunity to these strains. Therefore, understanding the shortterm evolution of influenza virus is crucial to developing seasonal vaccines. Conversely, vaccination of a population may influence the shortterm evolution of the virus, for example by decreasing the number of hosts in which the virus may replicate.
In general, during a single influenza season, more than one viral strain is circulating. It is known [2,3] that when suitable invasion conditions are satisfied, stable coexistence of two different strains is possible. The coexistence of two or more strains in a population will influence their individual evolution; and furthermore, the selective forces acting on the strains will be mediated by the transmission dynamics in the population. For example, the infection of hosts by one strain will reduce susceptibility to other strains, thereby limiting their spread in the population [4]. In addition, the time lag in emergence of a second strain following onset of an epidemic by a first strain will be influenced by the strategy and timing of vaccination [5]. It is clear that viral evolution and public health policy are strongly interconnected, and understanding the populationlevel dynamics of coexisting viral influenza infections, when vaccination of the population is to be undertaken, would be of great benefit in designing such vaccination strategies [6].
In [6], a homogeneous model of two viral strains was developed, incorporating crossimmunity and delay in emergence of the second strain. It was found that for small delay and large crossimmunity, infections with both strains appeared as a single epidemic wave; on the other hand, with sufficient delay, a second epidemic wave is possible. Further, for sufficient delay and high crossimmunity, the population of susceptible hosts may become so depleted as to prevent a second wave. These findings, together with possible impact of vaccination on antigenic drift, suggest that vaccination would be an important factor to include [6].
In large populations, contacts between individuals are not uniform, as assumed in the homogeneous model [6]. Typically, the number of contacts per day per individual is much smaller than the population size, and the structure of the corresponding ‘contact matrix’ plays an important role in the development of the pattern of the disease [7]. The effects of spatial correlations [8], such as occur when community structures are present [9], were illustrated in the spread of drug resistance in a network with mild clustering [10]: the spread of the resistant strain occurred more rapidly, and at significantly lower treatment levels, than was predicted by the homogeneous model.
The present paper extends the model in [6] in a number of ways. The model includes either prevaccination or vaccination during the epidemic, of a predetermined part of the population. The contact structure is modelled as a Markov network [11], in which the distribution of degrees of the nodes (i.e., number of contacts for individuals in the population) is specified. In addition, the model allows a distribution of ages in the population by incorporating a prescribed number of age classes. The Markov assumption for the contact network allows the specification of structural parameters such as assortativity [12] and clustering [1315] that are important characteristics of social groupings. These generalizations enable vaccination to be targeted according to age group and ‘contact number’ (degree of node), which in general respond to the vaccine in different ways. The model inevitably contains many parameters and allows a wide range of network structures to be specified; in addition, initial conditions can be specified in many different ways. Therefore, in this paper, only a simplified network model will be investigated. The structure of the network is comprised of uncorrelated nodes, with degree distribution specified as a truncated scalefree form [7]. Furthermore, for simplicity only one or two ageclasses are considered, where, for the latter, the median age is chosen to separate the two classes. While a detailed age distribution, characteristic of a real population, could be specified, the present results are intended to be illustrative only and to allow comparison with the corresponding homogeneous model. The network model can potentially be useful in describing specific populations, such as a small or large city, in which case the network structure and age distribution would need to be determined from statistical analysis of demographic and census data [16].
Section 2 describes the model in broad terms, and lists some of the parameter values used; technical details are given in the Appendix. Section 3 presents the results of simulations, in which the crossimmunity and delay in appearance of the second strain infection are varied. These results are also compared with those produced by the corresponding homogeneous model, to ascertain the importance of structure in the network for determining the timecourse and final extent (“total attack rate”) of the disease. Finally, Section 4 discusses these results, some possible extensions of the model, and implications for vaccination strategies in more realistic models based on specific demographic data.
Methods
The state flow diagram of the model is given in Figure 1, which represents either the population counts in various compartments in the homogeneous model or the state of any given node (labelled by infection state, degree and ageclass) in the network model. The model describes the evolution of two concurrent strains of influenza infection, over a duration short compared to the natural lifespan of an individual in the population, and for this reason birth and natural death processes are ignored, and furthermore the number of individuals in each age class remains constant. Since we consider a static contact network, the degree class of each individual is also fixed. Therefore, each individual in the population belongs to a unique class (k,a), where k denotes the number of contacts, and a the age class, and the total number of individuals in each (k,a) class is constant. S denotes the susceptibles and V denotes individuals receiving vaccination either prior to the onset of the first infection or after this onset.
Figure 1. State flow diagram for the twostrain influenza model. S denotes the susceptible state, without prior vaccination. Other susceptible individuals may receive vaccination prior to the onset of infection, or after infection has appeared, and are denoted by V. Despite vaccination, some individuals become infected and follow a similar sequence of infection states to that of the susceptibles. States (and parameters) originating from vaccination are denoted by a subscript ‘V’. States representing symptomatic infection by strain j are denoted by I_{j}, and correspondingly those infected asymptomatically by A_{j}. Doublesubscripted states indicate that the individual was previously infected with one strain, and is now progressing through infected states of the other strain. P_{j} denotes individuals partially recovered from infection by strain j, but still susceptible to infection by the other strain. R denotes the class recovered from successive infection by both strains. In addition, infected individuals may die (at rates d or d_{A}) and transfer (via the dashed arrows) to a diseaseinduced death class D (not shown). See text for details and explanation of the parameters.
Vaccination prior to the onset of infection is specified by the fraction of susceptibles in each age class receiving vaccination. For vaccination occurring during an outbreak, the following model is used: for individuals in any given (k,a) class, the rate of vaccination at any given time is (i) proportional to the current number of susceptibles in the class; (ii) an increasing function of the total current (symptomatic) infection in the population as a whole, saturating at a prescribed rate. This was done to attempt to model the social response to an outbreak in the population, in which the greater the number of infected individuals the more likely that susceptible individuals would avail themselves of existing vaccination opportunities. The precise mathematical specification of this response is given in the Appendix.
The baseline transmission rate of infection between a susceptibleinfected pair of individuals is denoted by τ. The actual rate will depend on the ageclasses that these individuals belong to, and whether the susceptible individual of the pair is seeing infection (by either strain) for the first or second time. These various possibilities are accounted for by expressing the actual transmission rate as τ times a factor, which depends on age classes involved, whether this is the first or second infection, and whether the individual has received prior vaccination. Details are given in the Appendix and in Table 1.
States labelled with I denote symptomatic infection, and those labelled with A denote asymptomatic infection. The P states describe immunity to one strain but not the other: P_{j} is the state with immunity to strain j (j = 1, 2), and R the state with immunity to both strains. In this model, we exclude coinfection: at any given time, an individual may be infected with at most one strain. State I_{j} denotes infection with strain j; and I_{jk} denotes previous infection with (and subsequent recovery from) strain j and current infection with strain k (where k ≠ j). A similar notation applies to the Aclasses. The efficacy of the vaccine against strain j is denoted by σ_{j}.
Subscript ‘V’ denotes states of infection (or partial recovery) arising from failure of the vaccine; and as before, labels states with infection due to, or partial recovery from, one of the strains. Following vaccination, infection due to strain j occurs with probability (1σ_{j}). In general, for seasonal influenza, the vaccine is targeted against the earlieroccurring strain 1 virus; its efficacy against the lateroccurring strain 2 (mutated) virus is expected to be less, i.e., σ_{2} < σ_{1}. As in [6], the delay T* in appearance of strain 2 in the population is a parameter of the model.
In Figure 1, the diverging pairs of directed edges are labelled with branching ratios for each strain of infection, with two pairs of such edges emanating from S and V classes. For example, if S is infected with one of the strains, it has a probability p of being symptomatically infected, and 1p of being asymptomatically infected. (We assume that p is the same for both strains). Since S may be infected with either strain, there are two pairs of branches emanating from S in Figure 1. Similarly, there are two branch pairs for V, representing infection due to failure of the vaccine.
After recovery from one strain of infection, an individual is still, in general, susceptible to infection by the other strain: individuals in state P_{j} (i.e., recovered from infection with strain j), can become infected with strain k (≠ j) but with diminished probability δ_{jk}. The probability of such infection being symptomatic is denoted by p_{jk}. Similarly, for individuals who have received prior vaccination but still become infected by strain j, the probability of strain k infection is denoted by p_{Vjk}. Finally, the model allows for the possibility of diseaseinduced death, denoted by the state D. The rates at which these occur are assumed to be d or d_{A} for symptomatic and asymptomatic infections, respectively, regardless of which of the disease states precede death; furthermore, the death rates  as with other parameters of the model – may depend on the age group in which the death occurs.
The converging directed edges in this Figure are labelled with the recovery rates from infection: either µ (symptomatic infection) or µ_{A} (asymptomatic infection), where we assume that these rates are the same for both strains, regardless of whether this is the first or second infection for that individual. The parameter values used in the simulations are given in Table 1.
For the homogeneous model, we may apply the technique of the nextgeneration matrix [17] to derive the basic reproductive number R_{0}. In general, the second strain appears after infection due to the first strain has begun, so that R_{0} can be calculated using a onestrain submodel. With this assumption, we find
where S_{0}, V_{0} denote the initial numbers of susceptible and vaccinated individuals, respectively, in the population, and β is the transmission coefficient. To establish a relationship between β and the baseline transmission rate τ between individuals in contact, we construct a (singleage class) network in which the ‘edge probability’ of randomly choosing an edge, one of whose vertices has degree k, is uniform. By relating this to the mean field model we derive (see Appendix)
where k_{1} = vertex degree of population subclass into which the Strain1 infection is introduced at time t = 0, and k_{max} = maximum vertex degree in the finite network (k_{max} = 20 in the simulations). If we choose for V_{0} = 0.2, a conservative value R_{0} = 1.9 for influenza [10], then using the above expressions for β and R_{0} we derive τ = 3.5 d^{1} for the transmission rate to be used in the simulations. The value of R_{0} corresponding to this τ in the absence of vaccination is R_{0} = 2.34.
In keeping with the definition of the two age class model (see Appendix), the estimates of death rates [18,19] arising from symptomatic or asymptomatic infection (d, d_{A}, respectively) for the two ageclass model correspond to the general population above and below the median of the age distribution P_{a} which, for the city of Vancouver, is about 38 years [20]. We assume that the death rates due to natural causes are negligible, and choose nominal values for the diseaseinduced rates: d(a_{1} ) = d(a_{2} ) = 0.002 d^{1} (Ref.[10]). These rates vary with the particular circulating influenza strains. Furthermore, we set d = d_{A} in this illustrative study.
In the model described above, the total number of individuals N_{k}_{,}_{a} in each (k,a) class is fixed, and hence the total population N (summed over all (k,a) classes) is constant. Therefore, by dividing the number of individuals in class (k,a) in state X at any given time by N, we may express the model in terms of the probability X_{k}_{,}_{a}(t) that a randomly chosen individual is in state X, and belongs to class (k,a), at time t. The resulting set of ordinary differential equations describing this deterministic model is given in the Appendix.
Results
The initial state was specified as follows. For prevaccination, a prescribed fraction V_{0} (a) of individuals in each age class a receive vaccination. Infection by strain 1 is introduced into fraction ε_{1} of the remaining susceptibles residing in a single class (k_{1},a_{1}). After the strain 1 infection has spread through the population for a time T*, a strain 2 infection is introduced into a fraction ε_{2} of class (k_{2},a_{2}) individuals. In the simulations, we use ε_{1} = ε_{2} = 0.5; k_{1} = 5, k_{2} = 10, and for the two age class model, a_{1} = 1, a_{2} = 2. As previously mentioned, it is assumed that no individual may be infected with both strains simultaneously. The simulations were performed using three models: (1) network model with two age classes; (2) network model with one age class; and (3) the homogeneouslymixing 'mean field' model. For (1) and (2), the structure of the network was chosen to have a scalefree form [7], with the number of individuals (nodes of the contact network) with k contacts being proportional to k^{2.5}[21], and 1 ≤ k ≤ k_{max} = 20. Furthermore, the degrees of the nodes of the network were assumed to be uncorrelated: although real networks show significant correlation structure – e.g., clustering and associativity [1215]  the purpose of the present simulations is to illustrate the general effects of departure from the homogeneous mixing assumption. A value R_{0} = 1.9 was fixed for the mean field model with V_{0} = 0.2.
Figure 2 compares the results from the network model using one age class (top row), two age classes (middle row), and the mean field model (bottom row), for delays of T* = 10 days (left column) and 60 days (right column) for introducing strain2 infection into a prescribed subpopulation. Solid curves correspond to a large (60%) crossimmunity (with δ ≡ δ_{12} =δ_{21} = 0.4) between strains, and dashed curves to a low (10%) crossimmunity (δ = 0.9). These results show that the level of crossimmunity has a significant effect on whether a second wave appears: if it is too high (60%, δ = 0.4), then the second wave does not appear in any of the three models. For the 2age class model (middle row), the (first) peak infection occurs at a slightly lower value for the larger crossimmunity: 0.03 vs. 0.033 fraction of the population; however, the timing of these peaks (~75 days after initial infection) is not sensitive to the level of crossimmunity. Similar conclusions apply to the one age class network model and the mean field model, though the timing of these peaks is different for different models.
Figure 2. Comparison of models: Varying crossimmunity and delay for introducing strain. 2 Comparison of results from network model for one age class (top row), two age classes (middle row), and mean field (bottom row) models. Solid curves are for δ = 0.4 (60% crossimmunity) and dashed curves for δ = 0.9 (10% crossimmunity). Plots in the left column are for the appearance of the second strain 10 days after strain1 infection starts; right column plots for the appearance of the second strain 60 days after initial strain1 infection. For both network models, strain 1 was introduced into class k = 5 and strain 2 into class k = 10. For the two age class model, infections were started in different age classes. Note that a second wave appears only when the crossimmunity is small (dashed curves), apart from a small, delayed outbreak in the oneage class model.
As expected, if the second strain is introduced after the strain 1 infection has been largely cleared from the population (as is the case when T* = 60 days: right column), then the first and second waves behave as distinct, noninteracting onestrain epidemics. However, when T* is only 10 days (left column), there is still a significant presence of strain 1 infection in the population: the infections in the two ageclass model merge into a single broad peak, whereas the other two models show two distinct peaks, with the second peak occurring in both models ~50 days after initial infection.
It is therefore apparent that the two age class network model exhibits a larger delay in peak infection – for both first and second waves – compared to the one age class and meanfield models. This can be accounted for by the reduced transmissibility between classes compared to within one class, as well as reduced transmissibility within the second age class (see the M matrix in the Appendix). (Recall that strain 1 and 2 infections are introduced into different age classes in the two age class network model). Such differences in delays between meanfield and structured models have been observed elsewhere [10], and underline the importance of spatial structure in determining the course of an epidemic.
Figure 3 shows the effect of different levels of prevaccination (V_{0} = 0.2 [dashed curves] or 0.4 [solid curves]) on the infection profiles, for each of the three models, assuming a delay T* = 10 days in appearance of strain 2 infection. For all models, the peak infection is reduced by about 50%, and occurs slightly later, when increasing from 20% to 40% coverage. The peak infections in the two age class model, however, occur at significantly lower values than in the other two models: by a factor of 5 for the mean field model, and a factor of 10 for the one age class network model. Again, this may be accounted for by the reduced transmissibility between different age classes. Notice also that, as in Figure 2, the two age class model exhibits a single peak of infection for both levels of vaccination.
Figure 3. Comparison of models: Priorvaccination coverage. Comparison of meanfield (a), one and two age class network models ((b) and (c), respectively), for different prior vaccination coverage of 40% (solid curves) or 20% (dashed curves). For the two age class model, strain1 and strain2 infections occur in different age classes. The parameter δ = 0.9 (10% crossimmunity), and second strain infection begins at T* = 10 days after firststrain infection.
In Figure 4, the two age class network model is used to explore the effects of vaccination during an epidemic outbreak, with no vaccination prior to the initial appearance of strain 1 infection, where vaccination rates are determined according to the “social response” to total infection in the population, as described in Section 2 and the Appendix. The resulting infection profiles are remarkably insensitive to (i) level of crossimmunity, (ii) delay T*, and (iii) the rate of vaccination ω_{0}. Apart from a temporary decrease in the rate of disappearance of the infection when T* = 10 days and δ = 0.4 (top left in Figure 4), there is only one peak of infection which occurs consistently at very similar times (50 days after strain 1 infection) and with peak magnitudes between 0.07 and 0.08 of the fraction of population. Comparing these results to the two age class model with prevaccination (middle row of Figure 2), it can be seen that the peak infection occurs about 2030 days earlier for vaccination during the epidemic than when prevaccination is carried out.
Figure 4. Comparison of models: Vaccination during outbreak. Comparison of total clinical infection resulting when vaccination is introduced during a disease outbreak, for different vaccination rates. Left column is for T* = 10 days, right column for T* = 60 days. Top row: δ = 0.4; bottom row: δ = 0.9. The underlying model is a two age class network model, and vaccination rates in each (k,a) class are proportional to the number of susceptible individuals in that class. These rates are determined by the total infection in the general population, and saturate at a value of ω_{0} per individual per day (see Appendix), where ω_{0} = 1 (solid curve), 0.5 (dashed curve) or 0 (dotted curve).
Tables 2, 3, 4 compare the attack rates when individuals are infected by each strain separately (represented by the final values of the P_{1} and P_{2} compartments in Figure 1) and by both strains in succession (represented by the R compartment), for each of the three models. It is seen that in general, for given T* and δ, the attack rates predicted by the network models are quite different between the model types. However for prevaccination, within each model, strain 1 attack rates increase, while strain 2 attack rates decrease, with the delay T*. For vaccination during the epidemic, these trends are also seen though for strain 1 infection are less pronounced. This is reasonable, because for longer T*, the strain 2 infection draws upon a smaller pool of susceptibles (and individuals with failed vaccination) which are depleted by strain 1 infection for a longer time interval before strain 2 appears. Also, as expected, strain 2 attack rates are sensitive to the level of crossimmunity of the strains, decreasing sharply as crossimmunity increases from 10% to 60%. Variations between models are manifestations of the importance of heterogeneity of contact structure. The dependence on age distribution between network models in Tables 2, 3 is a consequence of our assumption (see matrix M in the Appendix) that transmissibility within age class 1 is greater than that within age class 2 or between age classes.
Table 2. Total attack rates for 2age class network model
Table 3. Total attack rates for 1age class network model
Table 4. Total attack rates for mean field model
A question of importance to public health is the dependence of both the attack rate and death rate due to infection. The attack rates are shown in Table 5 for the two age class network model with T* = 60 days, and crossimmunity of 10% and 60%. For prevaccination, as expected, strain 1 attack rates decrease with increasing level of vaccination, at first slowly then dropping off sharply between V_{0} = 0.3 and 0.4 regardless of the level of crossimmunity. Interestingly, the maximum attack rate due to strain 2 infections is reached in this range, and drops off sharply thereafter. Thus, for the parameter values used, it appears that vaccination is most effective around these values, with diminishing returns for higher V_{0}. Tables 6, 7 show that death rates due to prevaccination are lower than predicted for the entire range of vaccination rates during an epidemic; and in both cases the death rates decrease with increasing levels of vaccination. Analogous to the attack rates, there is a small increase in mortality for vaccination coverage around 20%, due to increased numbers of strain 2 infections at the expense of reduced numbers of strain 1 infections; however, the mortality rate drops sharply once the vaccination coverage exceeds about 40%.
Conclusions
We have considered extensions of the two viral strain mean field (homogenous) model introduced in [6], to explore the effects of both local network structure and the division of the population into different age classes. The present study used model parameter values (in particular, R_{0}) originally estimated for mean field models; and in order to translate these to the network models derived in this paper, a correspondence was established between the mean field model and a limiting case of the network model (see Appendix). The two age class model assumed the age boundary was located at the median age (about 38 years for Vancouver), with vaccine efficacy of 80% in the lower age group and 40% in the upper age group.
Several notable features were observed when comparing the network models to the corresponding mean field case. Firstly, the amount of crossimmunity present is significant in determining whether a second wave of infection occurs. Due to a lower transmission rate between age classes and within the second age class, compared to within the first age class, infection levels were found to be significantly less for the two age class model than for either the one age class or mean field models. The infections occurred as either a single wave or as two successive waves. A second wave is more likely to occur the longer the delay in introduction of the second strain, since when this delay is short (~10 days) infections due to both strains merge into a single, broad peak. When a second wave does occur, the shapes of the two waves depend on when the second strain infection is introduced. If it occurs well after the first infection has run its course, then the two waves behave as distinct, noninteracting infections. The second infection peak is delayed, and its amplitude reduced, in the network model, compared to the mean field case. This behaviour reflects a longer propagation time in the network model, and has been qualitatively observed in other models, reinforcing the importance of including local network structure in realistic models.
As expected, the amount of crossimmunity between the two strains is important in determining the size of the secondstrain outbreak. It was found that its size decreased sharply with increasing crossimmunity. As the level of vaccination increases, strain 1 attack rates decrease, with a sharp drop occurring around 3040% prevaccination coverage; at the same time, strain 2 infections increase with increasing vaccination coverage, reaching their maximum somewhere in this range, and drop off sharply for higher coverage levels. This phenomenon is reminiscent of the development of drug resistance, where there is an optimal level of drug treatment (compare: vaccination coverage) that minimizes the overall infection [10]. This could have significant implications for vaccination strategies in realistic models of populations in which more than one strain is circulating.
It was found that increasing either prevaccination or vaccination during an outbreak, reduces the diseaseinduced mortality. Furthermore, prevaccination appears to be more effective than vaccination during an outbreak in reducing overall mortality, though this needs further investigation as it may depend critically on how the latter is implemented. This study considered only a simple model in which at any given time vaccination rates during an outbreak were governed by the total infection in the population at that time, and considers only vaccination of the susceptible class S, neglecting vaccination of other classes (e.g., P_{1} and P_{2} and asymptomatic cases).
As mentioned earlier, the particular form of the terms included in the model to incorporate local network structure and the effects of age classes was chosen for illustrative purposes. This approach, though, can be used on a specific population if sufficient data are available to determine realistic estimates of the age classes and network structure present and of the parameters of the model. The main difficulty is in determining the form of the twopoint correlations between vertices of the contact network for a realistic particular population, and this must be derived indirectly from estimates of network structure extracted from the data [16]. An intermediate approach is to explore the effects of a few network structure parameters – e.g., clustering, associativity, betweenness, and centrality [7,16], obtaining expressions for the twopoint probabilities defining the Markov network directly in terms of these parameters. This is currently under investigation.
Appendix: Effects of vaccination and population structure on influenza epidemic spread in the presence of two circulating strains
The various parameters in the model (Figure 1 of main text) are defined below:
τ = baseline transmission rate between a susceptibleinfected pair
p = probability of developing symptomatic infection with no prior exposure
p_{V}_{1}, p_{V}_{2} = probabilities of prevaccinated individuals developing symptomatic infection from strains 1 and 2, respectively, with no prior exposure
σ_{1}, σ_{2} = effectiveness of vaccine to strains 1 and 2, respectively
δ_{V}_{1}, δ_{V}_{2} = reduction in transmissibility of strains 1 and 2, respectively, for vaccinated individuals
p_{12} = probability of developing symptomatic infection with prior exposure to strain 1
p_{V}_{12} = probability of prevaccinated individuals developing symptomatic infection with prior exposure to strain 1
p_{21} = probability of developing symptomatic infection with prior exposure to strain 2
p_{V}_{21} = probability of prevaccinated individuals developing symptomatic infection with prior exposure to strain 2
δ_{A} = reduction in infectiousness due to asymptomatic infection
µ, µ_{A} = recovery rates from symptomatic and asymptomatic infections
δ_{12}, δ_{21} = level of crossimmunity induced by previous exposure to strain 1 and strain 2, respectively
d, d_{A} = diseaseinduced death rates, assumed to be agedependent but the same for each type of infection.
Define age classes 1 ≤ a ≤ a_{max}, and network degree classes 1 ≤ k ≤ k_{max}. Let
for U ∈ {A_{1}, A_{V}_{1},I_{1}, I_{V}_{1},A_{21}, A_{V}_{21},I_{21}, I_{V}_{21}, A_{2}, A_{V}_{2}, I_{2}, I_{V}_{2}, A_{12}, A_{V}_{12}, I_{12}, I_{V}_{12}}, denote the force of infection for ageclass a and degreeclass k. Here, M(a,a′) denotes the relative transmission coefficient between agegroups, so that τM(a,a′) = transmission coefficient between a susceptible individual of ageclass a in contact with an infected individual of ageclass a′. Also, P(k′,a′k,a) is the probability that an individual (node) of ageclass a and degreeclass k has a neighbour (adjacent node) of ageclass a′ and degreeclass k′.
In the special case that the contact network is the same for all ageclasses, P(k′,a′k,a) = P_{a}(a′a)P(k′k), where P_{a} ( a′a) denotes the probability that an ageclass a individual has an ageclass a′ neighbour. The two conditional distributions obey the conditions
In the general case, . If the nodedegrees are uncorrelated, then P(k′k) = P_{e}(k′), where P_{e}(k) is the edge distribution [22], defined as the probability of randomly drawing an edge connected to a vertex of degree k. It is related to P(k), the vertex distribution, by Similarly, if the age distributions are uncorrelated, then P_{a}(a′a) = P_{a}(a′). Thus, for uncorrelated agestructured networks, which are considered in this paper, P(k′,a′k,a) = P_{a}(a′)P_{e}(k′). In the present study, the degree distribution follows a scalefree form [7]P(k) ~ k^{}^{γ}, where we have chosen γ = 3.5. In this paper, we consider only one or two age classes (a_{max} = 1, 2); more realistic models would incorporate demographic data on several ages classes (typically, a_{max} ≥ 4), but as discussed in the main text, the model simulations are only intended to illustrate the effects of heterogeneous contact and agestructure, in comparison with homogeneous models. For the simulations reported in the main paper, we have for the oneage class model: P_{a}(a) = 1, and for the two age class model we chose P_{a}(a_{1}) = P_{a}(a_{2}) = 0.5, and
We extend the homogeneous (meanfield) model in [6] to a Markov network model with age structure, and include vaccination of strain 1 prior to onset of influenza outbreak. Furthermore, we extend this model to allow vaccination during an epidemic, by introducing a flux φ of susceptibles from the S classes to their corresponding V classes, according to the prescription:
(i) φ is a function of the total (symptomatic) infection in the population, I_{tot} (summed over all k and a), and φ = 0 when I_{tot} = 0;
(ii) φ is proportional to the population in class S(k,a,t);
(iii) φ eventually saturates at a maximum value as I_{tot} increases.
A functional form that satisfies these criteria is
where ω_{0} is the (ageclass dependent) saturation rate of vaccination, α is the value of I_{tot} at halfsaturation, and n > 0 governs the steepness of the response curve. In the simulations, α = 0.4 (i.e., halfsaturation occurs when 40% of population is infected), and n = 2.
In order to incorporate death due to infection, we add a set of classes {D(k,a,t)} to the model, and (similar to the recovery rates) assume that death rates are either d (for all symptomatic infections) or d_{A} (for all asymptomatic infections).
Based on Figure 1 (main text), this gives rise to the following set of 24×(a_{max}×k_{max}) equations:
where S ≡ S(k,a,t), etc., and the force of infection Θ_{U} may be expressed as
where U ∈ {A_{1}, A_{V}_{1},I_{1}, I_{V}_{1}, A_{21}, A_{V}_{21},I_{21}, I_{V}_{21},A_{2}, A_{V}_{2},I_{2}, I_{V}_{2},A_{12}, A_{V}_{12},I_{12},I_{V}_{12}}, and is the connectivity matrix, defined by the contact structure of the population. Also,
where
This shows that the various Θ(k,a,t)’s describe the connectivity of a vertex of degree k and ageclass a to all the infected adjacent vertices.
Comparison with mean field model
In order to make this comparison, we need to obtain a “limiting” form of the network model that approximates the mean field model. This will enable us to obtain a relationship between the mean field model transmission rate β and the transmission rate τ between a susceptibleinfected pair in the network. We consider a simplified network model consisting of a single age class, and an uncorrelated network (so that P(k’k)=P_{e}(k’)), and further that P_{e} is a uniform distribution: P_{e}(k) = 1/k_{max}, where k_{max} is the maximum degree in the network. With these assumptions, (A3) simplifies to
Therefore, the equation for S(k,t) becomes
where the term in parentheses is independent of k. Assume that S(k,t) = P(k)S(t), etc.; then this becomes
Summing over k from 1 to k_{max}, we obtain
whereis the mean degree of nodes in the network.
Comparing this approximation with the Mean Field expression for dS/dt, suggests we make the following correspondence:
More realistically, since we introduce Strain1 infections into the subpopulation defined by (k,a) = (k_{1},a_{1}), and taking account of the fact that R_{0} is defined in terms of the first generation of infection, it would be more accurate to replace by k_{1}, so that
Using this approximate relationship enables us to compare the simulation of the behaviours of the network and meanfield models, by relating numerical values of the parameters β and τ through the simplified limiting case of a network in which the probability of drawing an edge at random from the network is uniform.
Initial conditions for the mean field model
In what follows, it is assumed that the total population (including deaths) is normalized to unity, which is permissible since for this model it is constant. The initial conditions for the mean field model must be consistent with those of the network model. The analysis that follows applies to an arbitrary number of age classes and degree distributions.
In the network model, at t = 0 a fraction ε_{1} of Strain1 infection is introduced into the subpopulation in class (k_{1},a_{1}). Therefore, we specify the initial conditions at t = 0 according to
where P_{ka}_{1} ≡ P(k_{1} )P_{a}(a_{1} ) represents the fraction of the total population in class (k_{1},a_{1}). Similarly, at t = T* when strain 2 infections are introduced to a fraction ε_{2} of the subpopulation in class (k_{2}, a_{2}), the corresponding mean field conditions are modified as follows:
whereare the changes in the susceptible and vaccinated subpopulations, respectively, and P_{ka2} ≡P(k_{2} )P_{a}(a_{2} ) is the fraction of the total population in class (k_{2},a_{2}).
In order to allow comparison between Mean Field and network models, all agedependent parameters δ_{A}, δ_{V}_{1}, δ_{V2}, σ_{1}, σ_{2}, µ_{A}, µ, d_{A}, d, etc., in the network model are replaced by their agedistributed averages:, etc., where without risk of ambiguity we may drop the ‘MF’ superscript.
For the network model, for all age classes we set ε_{1} = ε_{2} = 0.5, p = 0.6, V_{0} = 0.2, σ_{1} = 0.8, σ_{2} = 0.4. For the two ageclass model, we chose (k_{1},a_{1}) = (5,1), (k_{2},a_{2} ) = (10,2); and for the oneage class model k_{1} = 5, k_{2} = 10. The (truncated) scalefree distribution P(k) ~ k^{}^{3.5} with k_{max} = 20 yields P_{ka}_{1} = 0.0067, P_{ka}_{2} = 0.0012 (so that, in a population of N = 10,000, the number of infections is 67 and 12, respectively), where we are assuming P_{a} to be uniformly distributed in the 2ageclass model: P_{a} (a_{1}) = P_{a}(a_{2} ) = 0.5.
Substituting these values into the expression for R_{0} (Section 2 in main paper), and using R_{0} = 1.9, V_{0} = 0.2, k_{max} = 20, and k_{1} = 5, yields the values β = 0.8765, and τ = 3.5 d^{1}. For V_{0} = 0.4, using the same value τ = 3.5 d^{1}, the corresponding value of R_{0} is 2.34.
Authors’ contributions
MEA wrote the bulk of the manuscript, and designed, implemented and simulated the network model from a homogeneous one constructed with assistance from Dr S. Moghadas. RK implemented and simulated the homogeneous model, and wrote the Conclusions. MEA wrote the Appendix.
Competing interests
The authors declare that they have no competing interests
Acknowledgements
It is with great sadness that I report the sudden passing of my dear friend and coauthor, Dr Randy Kobes, whilst this paper was in revision. He will be sorely missed by his many friends and colleagues. I would like to thank Dr Seyed Moghadas for his helpful comments and assistance in designing the mean field model, and the referee for incisive and useful comments that have greatly improved the paper. This work was supported by a Natural Sciences and Engineering Research Council of Canada Discovery Grant.
This article has been published as part of BMC Public Health Volume 11 Supplement 1, 2011: Mathematical Modelling of Influenza. The full contents of the supplement are available online at http://www.biomedcentral.com/14712458/11?issue=S1.
References

Boni MF: Vaccination and antigenic drift in influenza.
Vaccine 2008, 26(Suppl 3):C8C14. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

CastilloChavez C, Hethcote HW, Andreasen V, Levin SA, Liu WM: Epidemiological models with age structure, proportionate mixing, and crossimmunity.
J. Math. Biol. 1989, 27:233258. PubMed Abstract  Publisher Full Text

Dietz K: Epidemiologic interference of virus populations.
J. Math. Biol. 1979, 8:291300. PubMed Abstract  Publisher Full Text

Andreassen V, Lin J, Levin SA: The dynamics of cocirculating influenza strains conferring partial crossimmunity.
J. Math. Biol. 1997, 35:825842. PubMed Abstract  Publisher Full Text

Bowman CS, Arino J, Moghadas SM: Evaluation of vaccination strategies during pandemic outbreaks.
Math. Biosci. Eng. 2011, 8:113122. PubMed Abstract  Publisher Full Text

Moghadas SM, Bowman CS, Arino J: Competitive interference between influenza viral strains.

Newman MEJ: The structure and function of complex networks.
SIAM Rev. 2003, 45:167256. Publisher Full Text

Keeling M: The effects of local spatial structure on epidemiological invasions.
Proc Biol Sci. 1999, 266:859867. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Fortunato S: Community detection in graphs.
Physics Reports 2010, 486:75174. Publisher Full Text

Alexander ME, Dietrich SM, Hua Y, Moghadas SM: A comparative evaluation of modelling strategies for the effect of treatment and host interactions on the spread of drug resistance.
J. Theor. Biol. 2009, 259:253263. PubMed Abstract  Publisher Full Text

Boguñá M, PastorSatorras R, Vespignani A: Epidemic spreading in complex networks with degree correlations. In Lecture Notes in Physics. Volume 625. Edited by PastorSatorras R, Rubi M, DiazGuilera A. New York: Springer; 2003::127147. Publisher Full Text

Newman MEJ: Assortative mixing in networks.
Phys. Rev. Letts 2002, 89:208701. Publisher Full Text

Park J, Newman MEJ: Solution for the properties of a clustered network.
Phys. Rev. E 2005, 72:026136. Publisher Full Text

Serrano MÁ, Boguñá M: Clustering in complex networks. I. General formalism.
Phys. Rev. E 2006, 74:056114. Publisher Full Text

Serrano MA, Boguñá M: Clustering in complex networks. II. Percolation properties.
Phys. Rev. E 2006, 74:056115. Publisher Full Text

Kolaczyk ED: Statistical Analysis of Network Data. New York: Springer; 2009. PubMed Abstract  Publisher Full Text

van den Driessche P, Watmough J: Reproduction numbers and subthreshold endemic equilibria for compartmental models of disease transmission.
Math. Biosci 2002, 180:2948. PubMed Abstract  Publisher Full Text

Alexander ME, Bowman C, Moghadas SM, Summers R, Gumel AB, Sahai BM: A vaccination model for transmission dynamics for influenza.
SIAM J. Appl. Dyn. Sys. 2004, 3:503524. Publisher Full Text

Dushoff J: Mortality due to Influenza in the United States–An Annualized Regression Approach Using MultipleCause Mortality Data.
American Journal of Epidemiology 2005, 163:181187. PubMed Abstract  Publisher Full Text

Data from the Canada Census, Population by Age In Community Services Group. Vancouver; 1996.
City of Vancouver. Vancouver Local Areas 1996

PastorSatorras R, Vespignani A: Epidemic spreading in scalefree networks.
Phys. Rev. Letts. 2001, 86:32003203. Publisher Full Text

Weber S, Porto M: Generation of arbitrarily twopointcorrelated random networks.
Phys. Rev. E 2003, 76:046111. Publisher Full Text

Arino J, Bowman C, Moghadas SM: Antiviral resistance during pandemic influenza: implications for stockpiling and drug use.
BMC Infect. Dis. 2009.
DOI:10.1186/1471233498
PubMed Abstract  Publisher Full Text  PubMed Central Full Text 
Cauchemez S, Carrat F, Viboud C, Valleron AJ, Boelle PY: A Bayesian MCMC approach to study transmission of influenza : application to household longitudinal data.
Stat.Med. 2004, 23:34693487. PubMed Abstract  Publisher Full Text

Gojovic MZ, Sander B, Fisman D, Krahn MD, Bauch CT: Modelling mitigation strategies for pandemic (H1N1).
CMAJ 2009.
DOI:10.1503/cmaj.091641
PubMed Abstract  Publisher Full Text  PubMed Central Full Text