Modeling optimal treatment strategies in a heterogeneous mixing model
 Seoyun Choe^{1} and
 Sunmi Lee^{2}Email author
https://doi.org/10.1186/s129760150026x
© Choe and Lee. 2015
Received: 20 September 2015
Accepted: 16 November 2015
Published: 25 November 2015
Abstract
Background
Many mathematical models assume random or homogeneous mixing for various infectious diseases. Homogeneous mixing can be generalized to mathematical models with multipatches or age structure by incorporating contact matrices to capture the dynamics of the heterogeneously mixing populations. Contact or mixing patterns are difficult to measure in many infectious diseases including influenza. Mixing patterns are considered to be one of the critical factors for infectious disease modeling.
Methods
A twogroup influenza model is considered to evaluate the impact of heterogeneous mixing on the influenza transmission dynamics. Heterogeneous mixing between two groups with two different activity levels includes proportionate mixing, preferred mixing and likewithlike mixing. Furthermore, the optimal control problem is formulated in this twogroup influenza model to identify the groupspecific optimal treatment strategies at a minimal cost. We investigate groupspecific optimal treatment strategies under various mixing scenarios.
Results
The characteristics of the twogroup influenza dynamics have been investigated in terms of the basic reproduction number and the final epidemic size under various mixing scenarios. As the mixing patterns become proportionate mixing, the basic reproduction number becomes smaller; however, the final epidemic size becomes larger. This is due to the fact that the number of infected people increases only slightly in the higher activity level group, while the number of infected people increases more significantly in the lower activity level group. Our results indicate that more intensive treatment of both groups at the early stage is the most effective treatment regardless of the mixing scenario. However, proportionate mixing requires more treated cases for all combinations of different group activity levels and group population sizes.
Conclusions
Mixing patterns can play a critical role in the effectiveness of optimal treatments. As the mixing becomes more likewithlike mixing, treating the higher activity group in the population is almost as effective as treating the entire populations since it reduces the number of disease cases effectively but only requires similar treatments. The gain becomes more pronounced as the basic reproduction number increases. This can be a critical issue which must be considered for future pandemic influenza interventions, especially when there are limited resources available.
Keywords
A twogroup influenza model Heterogeneous mixing Preferred mixing Optimal control theory Targeted treatment strategiesBackground
The timely and effective countermeasures of influenza challenge global health experts around the world, especially when limited resources are available. Mathematical modeling has made significant contributions to understanding the spread of influenza, and also providing useful insights to control or decrease the disease burden [1–4]. A number of mathematical models assume random or homogeneous mixing for the influenza dynamics, which can provide a good approximation to real epidemiological phenomenon [5, 6]. This simple assumption of homogeneous mixing can be extended to more general mathematical models with multipatches or age structure by incorporating contact matrices to capture the dynamics of the heterogeneously mixing populations [1, 7]. In general, contact or mixing patterns are difficult to measure in many infectious diseases including influenza. There is no doubt that mixing patterns are considered to be one of the critical factors for infectious disease modeling.
There are many different approaches that allow us to investigate the impact of contact patterns on the transmission dynamics of infectious diseases. The agedependent contact matrices based on empirical social data have been estimated [8, 9]. Agedependent transmission matrices that describe the mixing and the probability of infection are studied using synthetic data [10]. In these studies, it has been noted that contact patterns are strongly dependent on distinct age groups, and therefore, the heterogeneity of contact patterns should be recognized as an important feature for the realistic modeling of many infectious diseases. Moreover, individual based models or network based models can provide more details on the disease dynamics by studying the effects of heterogeneous and clustered contact patterns. Contact patterns and their underlying network structures have shown to be one of the critical factors for determining the characteristics of infectious disease transmission [11, 12]. Also, several different heterogeneity types in infectious disease models have been incorporated, such as susceptibility, infectivity and mixing patterns [13, 14]. Agent or network based models have been developed to study effective controls in the influenza pandemic [2, 15, 16].
Deterministic models which are much simpler than network based models have been successfully employed to study the transmission dynamics of various infectious diseases and continued to produce valuable insights. Preferred mixing has been used to highlight the role of contact patterns in the HIV transmission dynamics [17, 18]. The impact of selective mixing is studied in the transmission of STDs [19]. In these studies, the contact rate matrices are formulated in terms of activity levels and subpopulation sizes by using a proportionate mixing assumption. The relation between the basic reproduction number and the initial exponential growth rate of an epidemic to models with heterogeneous mixing has been studied [20–22]. The authors show that an epidemic with heterogeneous mixing may have a quite different epidemic size than an epidemic with homogeneous mixing, even though they may have the same reproduction number and initial exponential growth rate. Determination of the final size of an epidemic under the assumption of heterogeneous mixing requires additional data from the initial exponential growth stage of the epidemic [21]. More recently, a twogroup influenza model has been used to study the impact of heterogeneous mixing on the probability of the extinction of influenza [23, 24]. It has been pointed out that heterogeneous mixing between two subgroups would play a key role to explain the delays in the geographic spread of the 2009 H1N1 pandemic observed in Mexico and Japan.
In this manuscript, a deterministic twogroup model is used to study the influenza transmission dynamics in heterogeneous environments. A twogroup model allows for different activity levels and heterogeneous mixing between subgroups. In particular, two groups are coupled by a mixing matrix whose entries p _{ ij }, i, j=1,2 represent the proportion of individuals in group i that contact individuals in the other group j. This twogroup model of influenza is an extension of the prototype model by Brauer [20] in which a twogroup model is used to investigate the impact of proportionate mixing on the basic reproduction number and the final epidemic size. Now, our model involves more extensive heterogeneous mixing scenarios between the two groups. Specifically, several mixing scenarios are considered, including proportionate mixing, preferred mixing and likewithlike mixing, by varying the group mixing fractions.
We explore how this mixing pattern can affect the basic reproduction number and the final epidemic size. Heterogeneous mixing certainly changes the reproduction number and the final epidemic size, but it is not trivial to determine whether different mixing assumptions can change them substantially or not. The level of transmissibility measured by the basic reproduction number \(\mathcal {R}_{0}\) is varied to highlight the differences and similarities for the results under several heterogeneous mixing scenarios. Moreover, we formulate an optimal control framework to investigate how these mixing patterns will influence the effectiveness of groupspecific treatment strategies in the twogroup model. Under various mixing scenarios, optimal groupspecific treatment strategies and the corresponding influenza outcomes are compared. This can help us address some of the important issues such as allocating optimal treatments for future pandemic preparedness plans.
Methods
A heterogeneous mixing model
For each group i, κ _{ i } is the rate of passage from the latent to the symptomatic infective or asymptomatic infective classes; p is the fraction of latent members who become symptomatic infectious, and the fraction (1−p) progress to the asymptomatic stage; δ is the infectivity reduction factor for the asymptomatic class and σ is the infectivity reduction factor for treated members; α _{ i }(α _{ T,i }) is the natural recovery rate from the infected (treated) to the removed stage and η _{ i } is the rate of passage from the asymptomatic to the removed stage. Also, u _{ i } is a constant treatment rate, which will be modified as a timedependent treatment rate in the next section. d _{ i } (d _{ T,i }) is the diseaseinduced death rate from the infected class (the treated class).
More details on the preferred mixing formulation can be found in previous studies [18, 20].
The impact of mixing patterns on the contact matrix
When the group mixing fraction is π _{1}=π _{2}=0, we have proportionate mixing which is a special case of preferred mixing (C _{1}). It is also possible to have likewithlike mixing when π _{1}=π _{2}=1, in which members of each group mixes only with members of the same group. That is, for likewithlike mixing, p _{11}=p _{22}=1 and p _{12}=p _{21}=0 (C _{5}). For likewithlike mixing, the contact matrix is a diagonal matrix.
The basic reproduction number
One of the most important factors in mathematical epidemiology is the basic reproduction number, which is the average number of secondary infectious cases when one infectious individual is introduced to a whole susceptible population. The basic reproduction number can be calculated by using the next generation matrix approach, outlined in [25, 26]. Since the model includes treatments, we also compute the controlled reproduction number in the presence of constant treatment rates (u _{ i }).
Let x=(L _{1},I _{1},A _{1},T _{1},L _{2},I _{2},A _{2},T _{2})^{ T } and F(x) represent all the new infection rates. The net transition rates out of the corresponding compartment are represented by V(x). Then, we find the Jacobian matrix of F(x) and V(x) evaluated at the diseasefree equilibrium point x ^{∗}, which consists of S _{1}=N _{1},S _{2}=N _{2} and the rest of the components zero. The spectral radius of the matrix F V ^{−1} yields the basic reproduction number (\(\mathcal {R}_{0}\)) and the controlled reproduction number (\(\mathcal {R}_{c}\)) in the presence of treatments (more details are given in Appendix A).
The expressions for the basic reproduction number \(\mathcal {R}_{0}\) and the controlled reproduction number \(\mathcal {R}_{c}\) have been generalized from the ones with proportional mixing [20] to the ones with preferred mixing. The activity levels and the mixing fractions play a critical role in the basic reproduction number.
For instance, taking partial derivatives of p _{ ij } with respect to π _{ i }, we can show that p _{11} and p _{22} increase and p _{12} and p _{21} decrease as either π _{1} or π _{2} increase to 1. This results in the basic reproduction number increasing as preferred mixing becomes likewithlike mixing. Numerical sensitivity analysis of \({\mathcal {R}}_{0}\) and \({\mathcal {R}}_{c}\) is carried out in the next section.
The final size relation
When these values are substituted into the final sized system, S _{1}(∞) and S _{2}(∞) can be expressed in terms of the model parameters. As seen in the analytic expression above, the group specific final sizes are coupled with each other in a complex way. In the previous study [21], it can be simplified under proportionate mixing and shown that the final epidemic size in group 1 is larger than in group 2 when a _{1}>a _{2}. Also, it has been pointed out that \(\mathcal {R}_{0}\) alone is not enough to determine the final epidemic size due to this complex coupling. It is difficult to observe how different mixing patterns affect the final size relation. Therefore, we carry out sensitivity analysis numerically as mixing patterns are varied in the following section. The details on the computation of the final size relation are given in Appendix A and the references [20, 21].
Modeling optimal treatment strategy
Optimal control theory has been used frequently in a number of biological and epidemiological models (see [29] and the references therein). For influenza transmission models, optimal interventions are identified and the impact of optimal interventions on the influenza dynamics are investigated [30–32]. Various intervention strategies such as vaccination, antiviral treatment, and isolation controls are studied; optimal strategies for the 1918 influenza pandemic with limited resources [33] and agedependent optimal vaccination strategies are investigated in context of the transmission dynamics of the 2009 influenza pandemic [34, 35].
We choose to model the control efforts via a linear combination of quadratic terms, u _{ i } ^{2}(t)(i=1,2). The constants C _{1},C _{2} are the weight constants for infected individuals and W _{1},W _{2} are the relative costs of the interventions. We might include the cost of deaths in the objective functional so that we would emphasize the cost of diseaseinduced deaths. However, it turns out that the results including the cost of deaths and the ones without the cost of deaths are almost indistinguishable (results are not shown here).
where Ω={(u _{1}(t),u _{2}(t))∈(L ^{1}(0,T))^{2}∥0≤u _{1}(t),u _{2}(t)≤b, t∈ [ 0,T]} subject to the state equations given by (1) with initial conditions. The existence of optimal controls is guaranteed from standard results on optimal control theory [36]. Pontryagin’s Maximum Principle is used to establish necessary conditions that must be satisfied by an optimal solution [37]. Derivations of the necessary conditions are shown in Appendix B.
Parameter definitions and baseline values used in the numerical simulations
Parameter  Description  Values 

a _{ i }  Groupspecific activity level for age group i  0.1−0.84 
π _{ i }  Groupspecific mixing fraction for age group i  0−1 
α _{ i }  Recovery rate for infected class for age group i (days ^{−1})  0.244 
η _{ i }  Rate of progression from asymptomatic to recovered class for age group i (days ^{−1})  0.244 
α _{ T,i }  Recovery rate for treated class for age group i (days ^{−1})  0.323 
κ _{ i }  Rate of progression from latent to infective or asymptomatic class (days ^{−1})  0.526 
p  Fraction of latent individuals who become infected  0.667 
σ  Infectivity reduction for treated class  0.2 
δ  Infectivity reduction for asymptomatic class  0.5 
d _{ i }  Mortality rate for infected class for age group i (days ^{−1})  0.01 (i=1) 
0.08 (i=2)  
d _{ T,i }  Mortality rate for recovered class for age group i (days ^{−1})  0.005 (i=1) 
0.04 (i=2)  
b  The upper bound of control  0.2, 0.5, 1 
C _{ i }  Weight constants on I _{ i }(i=1,2)  1 
W _{ i }  Weight constants on controls (i=1,2)  1,100 
N _{1}  Population size for group 1  1750 
N _{2}  Population size for group 2  250 
Results and discussion
We present the numerical simulations associated with implementing optimal treatment control functions as well as their effect on twogroup influenza dynamics under different mixing patterns. In order to investigate the impact of mixing patterns, the group mixing fractions π _{ i } are varied from 0 to 1 for each i=1,2, including proportionate mixing (π _{ i }=0), half mixing (π _{ i }=0.5) and likewithlike mixing (π _{ i }=1).
The results in the absence of treatments
These results are dependent on the group activity level, the group mixing fraction, and the group population size. We present a summary of the impact of these parameters on the final attack ratio as group activity levels and group population sizes are varied from the baseline scenarios in the next section.
The impact of mixing patterns on the groupspecific optimal treatment
Figure 5 illustrates the impact of mixing patterns on the groupspecific optimal treatment controls and the proportion of incidence under the three mixing patterns. Note that the time window of treatment is wider for group 1 under all mixing patterns while the time period of treatment becomes smaller for group 2 as the mixing becomes likewithlike (top panels). This is due to the fact that group 1 has a higher activity level and a larger population size than group 2. Also, the cumulative treated proportion becomes smaller as mixing becomes likewithlike for group 2, while it is the opposite for group 1. As noted in the previous section, the incidence in the lower activity level group 2 gets significantly larger as the mixing becomes more proportionate than the one in the higher activity group 1. Hence, this leads to the final epidemic size getting smaller and it requires less treatment as mixing becomes likewithlike mixing (bottom panels).
The impact of mixing patterns on the groupspecific final epidemic size
The impact of control parameters
The impact of different activity levels and subpopulation sizes

Baseline scenario: a _{1}>a _{2}, N _{1}>N _{2}

Scenario 1: a _{1}>a _{2}, N _{1}=N _{2}

Scenario 2: a _{1}>a _{2}, N _{1}<N _{2}

Scenario 3: a _{1}=a _{2}, N _{1}=N _{2}

Scenario 4: a _{1}<a _{2}, N _{1}<N _{2}
1. Let us consider the first scenario for a _{1}>a _{2} and N _{1}=N _{2}: the number of infected and treated cases for group 1 under the baseline scenario is larger than the ones under the first scenario, while it is the opposite for group 2 under all mixing patterns. However, the total number of infected and treated cases is larger under the baseline scenario than the first scenario. When mixing is proportionate (π _{ i }=0), the total cumulative treated (infected) number of the baseline scenario is 403 (894) and the one for the first scenario is 376 (837), respectively. Now, when the mixing is likewithlike (π _{ i }=1), the total cumulative treated number of the baseline scenario is 399 (885) and the one for the first scenario is 236 (522), respectively. Interestingly, the baseline scenario requires more treatment, but the number of infected people is larger than the first scenario, and the difference becomes more significant as mixing becomes more likewithlike.
2. The second scenario is a _{1}>a _{2} and N _{1}<N _{2}: this scenario is for the case where the higher activity group has a smaller population size than the lower activity group. The total cumulative treated (infected) number is 360 (799) for proportionated mixing and 73 (159) for likewithlike mixing, respectively. This result demonstrates that optimal treatment has more significant impact to reduce the epidemic size as mixing become likewithlike when a higher activity group has a smaller subpopulation size. The time using treatment control for group 2 is longer than the one of group 1 under proportionate mixing. However, the time using treatment of group 2 decreases by increasing π _{1}, so that the time of treatment for group 1 is longer than the one for group 2 as mixing becomes likewithlike. This suggests that the time duration of implementing treatment depends on the group activity levels (the higher activity requires a longer time to implement).
3. The third scenario is that the activity level and the population size for group 1 and group 2 are the same. The basic reproduction number, the final attack ratios and optimal groupspecific treatments are almost the same under all mixing patterns. It shows that the impact of mixing patterns is not significant when groups have the same activity level and the same population size. Lastly, the scenario 4 for a _{1}<a _{2} and N _{1}<N _{2} is exactly identical as the baseline scenario, hence, the results are omitted.
Conclusions
We have studied the dynamics of influenza transmission in a twogroup model in which two groups are connected via a mixing matrix. The model proposed here represents two age groups with different activity levels and distinct mixing patterns. Several mixing patterns are considered such as proportionate and preferred mixing by varying the group mixing fractions π _{ i } for i=1,2. The impact of these mixing patterns is illustrated on the basic reproduction number and the group specific final epidemic size. Also, the intensity of \({\mathcal {R}}_{0}\) is varied by using different values of the group specific activity levels.
The basic reproduction number \({\mathcal {R}}_{0}\) increases as the mixing becomes likewithlike mixing. Interestingly, in the absence of treatments, the opposite is true for the final epidemic size, which gets smaller as mixing becomes likewithlike as reported [23]. This is consistent with the observations that the basic reproduction number alone is not enough to determine the final epidemic size in a heterogeneous model [21]. However, the basic reproduction number and the final epidemic size depend on the group activity level and the group population size as well. Under our baseline scenarios (a _{1}>a _{2} and N _{1}>N _{2}), the final attack ratios decreases as π _{2} increases to 1, which implies that the group 2 mixing fraction determines the order of the final attack ratio. Using different sets of parameters, the order might changes depending on either π _{1} or π _{2}.
Furthermore, we formulated an optimal framework to investigate groupspecific optimal treatment strategies under various mixing scenarios. For a moderate value of \(\mathcal {R}_{0}\), the optimal treatment can prevent the outbreaks in both groups under all mixing patterns. The treatment time is longer for group 1 and the treated cases are larger since it has a higher activity level and a large population size, regardless of mixing patterns. For a higher value of \(\mathcal {R}_{0}\), the optimal treatment can not stop the outbreaks in both groups under all mixing patterns. Compared with the moderate \(\mathcal {R}_{0}\) results, the treatment time for group 1 decreases, but for group 2 increases. This is due to the fact that the number of those infected in the lower activity group gets significantly larger as the mixing becomes more proportionate. Therefore, treating more people is necessary with emphasis on the lower activity group, when two groups mix proportionately.
Optimal treatment strategies can significantly limit the severity of outbreaks when \({\mathcal {R}}_{0}\) is brought below a certain threshold (the controlled reproduction number, \({\mathcal {R}}_{c}\)). Under optimal treatments in both groups, the controlled reproduction number \({\mathcal {R}}_{c}\) and the final attack ratio decrease slightly as mixing becomes likewithlike. Again, the final attack ratio is in the order of either π _{1} or π _{2}. Preferred mixing changes the basic reproduction number, the controlled reproduction number and the final epidemic size in a rather complex way and the effect gets more substantial as the epidemic gets more severe.
Further, the impact of different group activity levels and subpopulation sizes are explored for the optimal treatments and the resulting two group influenza dynamics. First, the group mixing fractions can play a key role in the final attack ratio. For the case of a _{1}>a _{2} and N _{1}=N _{2}, π _{2} determines the order of the final attack ratio, i.e., the total final attack ratio becomes smaller as π _{2} becomes 1. For the case of a _{1}>a _{2} and N _{1}<N _{2}, the total final attack ratio becomes smaller as π _{1} becomes 1. Based on these results, the effectiveness of optimal treatments is dependent on the group specific parameters. In general, the optimal treatment becomes more efficient as the mixing becomes more likewithlike. The efficiency of optimal treatments becomes more substantial when a higher activity group has a smaller population size (a _{1}>a _{2} and N _{1}<N _{2}). Also, the time duration of implementing treatment depends on the group activity levels while the final attack ratios are more sensitive to the group population size.
Sensitivity analysis for the control upper bound and the weight constant has been carried out. The results using a lower upper bound (b=0.2) and a higher weight constant (W=100) show that the magnitude of the treatment controls decreases and therefore, the total amount of treated and infected cases are increased in both groups under all mixing patterns. Clearly, this indicates that a more intensive treatment or a higher treatment rate is able to more efficiently reduce the total number of infected individuals with less treatment.
For the parameters used here, our results indicate that treatment of both groups with a higher rate is the most effective, regardless of mixing scenarios. However, proportionate mixing requires more treated cases for all combinations of different group activity levels and group population sizes. In other word, as the mixing becomes more likewithlike mixing, treatment of the more active group in the population is almost as effective as treating the entire population, since it reduces the number of disease cases effectively but requires the similar treatments. The gain is more pronounced as the basic reproduction number increases. This can be a critical issue which has to be considered for future epidemic interventions, especially when there are limited resources.
This study focuses on the twogroup influenza model to explore the effect of heterogeneous mixing on the groupspecific optimal treatment. This simple twogroup influenza model can be used for any general disease which consists of two different activity levels and different mixing patterns. Furthermore, this work can be generalized to a multigroup influenza model (with more age groups) so that it can capture more interesting and realistic epidemiological scenarios. This will be carried out in our future study.
Appendix A
where \(\Gamma _{1}=\left (\frac {\delta (1p)}{\eta _{1}}+\frac {p(\alpha _{T,1}+\sigma u_{1})}{\alpha _{T,1}(\alpha _{1}+u_{1})} \right), \Gamma _{2}=\left (\frac {\delta (1p)}{\eta _{2}}+\frac {p(\alpha _{T,2}+\sigma u_{2})}{\alpha _{T,2}(\alpha _{2}+u_{2})} \right)\).
Appendix B
From this Hamiltonian and Pontryagin’s Maximum Principle [37], we obtain the following theorem:
Theorem.
Proof.
By using the standard argument for bounds a≤u _{ i }(t)≤b for i=1,2, we have the optimality conditions (15). □
Declarations
Acknowledgements
This work was supported a grant from Kyung Hee University in 2013 (KHU20130683).
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License(http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
Authors’ Affiliations
References
 Anderson RM, May RM, Anderson B. Infectious Diseases of Humans: Dynamics and Control. Oxford: Oxford university press; 1970.Google Scholar
 Ferguson NM, Cummings DAT, Cauchemez S, Fraser C, Riley S, Meeyai A, et al. Strategies for containing an emerging influenza pandemic in southeast asia. Nature. 2005; 437(7056):209–14.View ArticlePubMedGoogle Scholar
 Heesterbeek JAP, Vol. 5. Mathematical Epidemiology of Infectious Diseases: Model Building, Analysis and Interpretation. New York: John Wiley & Sons; 2000.Google Scholar
 Longini IM, Nizam A, Xu S, Ungchusak K, Hanshaoworakul W, Cummings DAT, et al. Containing pandemic influenza at the source. Science. 2005; 309(5737):1083–7.View ArticlePubMedGoogle Scholar
 Arino J, Brauer F, Van den Driessche P, Watmough J, Wu J. Simple models for containment of a pandemic. J R Soc Interface. 2006; 3(8):453–7.PubMed CentralView ArticlePubMedGoogle Scholar
 Arino J, Brauer F, Van Den Driessche P, Watmough J, Wu J. A model for influenza with vaccination and antiviral treatment. J Theor Biol. 2008; 253(1):118–30.View ArticlePubMedGoogle Scholar
 Hethcote HW, Van Ark JW. Epidemiological models for heterogeneous populations: proportionate mixing, parameter estimation, and immunization programs. Math Biosci. 1987; 84(1):85–118.View ArticleGoogle Scholar
 Mossong J, Hens N, Jit M, Beutels P, Auranen K, Mikolajczyk R, et al. Social contacts and mixing patterns relevant to the spread of infectious diseases. PLoS Med. 2008; 5(3):74.View ArticleGoogle Scholar
 Wallinga TPJ, Kretzschmar M. Using data on social contacts to estimate agespecific transmission parameters for respiratoryspread infectious agents. Am J Epidemiol. 2006; 164(10):936–44.View ArticlePubMedGoogle Scholar
 Del Valle SY, Hyman JM, Hethcote HW, Eubank SG. Mixing patterns between age groups in social networks. Soc Netw. 2007; 29(4):539–54.View ArticleGoogle Scholar
 Bansal S, Grenfell BT, Meyers LA. When individual behaviour matters: homogeneous and network models in epidemiology. J R Soc Interface. 2007; 4(16):879–91.PubMed CentralView ArticlePubMedGoogle Scholar
 Volz EM, Miller JC, Galvani A, Meyers LA. Effects of heterogeneous and clustered contact patterns on infectious disease dynamics. PLoS Comput Biol. 2011; 7(6):1002042.View ArticleGoogle Scholar
 Yates A, Antia R, Regoes RR. How do pathogen evolution and host heterogeneity interact in disease emergence?Proc R Soc London B: Biol Sci. 2006; 273(1605):3075–83.View ArticleGoogle Scholar
 Apolloni A, Poletto C, Ramasco J, Jensen P, Colizza V. Metapopulation epidemic models with heterogeneous mixing and travel behaviour. Theor Biol Med Modell. 2014; 11(1):3.View ArticleGoogle Scholar
 Gani R, Hughes H, Fleming D, Griffin T, Medlock J, Leach S. Potential impact of antiviral drug use during influenza pandemic. Emerg Infect Dis. 2005; 11(9):1355–62.PubMed CentralView ArticlePubMedGoogle Scholar
 Longini IM, Halloran ME, Nizam A, Yang Y. Containing pandemic influenza with antiviral agents. Am J Epidemiol. 2004; 159(7):623–33.View ArticlePubMedGoogle Scholar
 Jacquez JA, Simon CP, Koopman J, Sattenspiel L, Perry T. Modeling and analyzing hiv transmission: the effect of contact patterns. Math Biosci. 1988; 92(2):119–99.View ArticleGoogle Scholar
 Nold A. Heterogeneity in diseasetransmission modeling. Math Biosci. 1980; 52(3):227–40.View ArticleGoogle Scholar
 Hyman JM, Li J. Behavior changes in sis std models with selective mixing. SIAM J Appl Math. 1997; 57(4):1082–94.View ArticleGoogle Scholar
 Brauer F. Epidemic models with heterogeneous mixing and treatment. Bull Math Biol. 2008; 70(7):1869–85.View ArticlePubMedGoogle Scholar
 Brauer F. Heterogeneous mixing in epidemic models. Can Appl Math Q. 2012; 20(1):1–13.Google Scholar
 Ma J, Earn DJD. Generality of the final size formula for an epidemic of a newly invading infectious disease. Bull Math Biol. 2006; 68(3):679–702.View ArticlePubMedGoogle Scholar
 Nishiura H, Chowell G, Safan M, CastilloChavez C. Pros and cons of estimating the reproduction number from early epidemic growth rate of influenza a (h1n1) 2009. Theor Biol Med Model. 2010; 7(1):1.PubMed CentralView ArticlePubMedGoogle Scholar
 Nishiura H, Cook AR, Cowling BJ, Ramasco J, Jensen P, Colizza V. Assortativity and the probability of epidemic extinction: A case study of pandemic influenza a (h1n12009). Interdiscip Perspect Infect Dis. 2011:194507. doi:10.1155/2011/194507.
 Diekmann O HJ, MG R. The construction of nextgeneration matrices for compartmental epidemic models. J R Soc Interface. 2010; 7(47):873–85.PubMed CentralView ArticlePubMedGoogle Scholar
 Van den Driessche P, Watmough J. Reproduction numbers and subthreshold endemic equilibria for compartmental models of disease transmission. Math Biosci. 2002; 180(1):29–48.View ArticlePubMedGoogle Scholar
 Arino J, Brauer F, Van Den Driessche P, Watmough J, Wu J. A final size relation for epidemic models. Math Biosci Eng. 2007; 4(2):159.View ArticlePubMedGoogle Scholar
 Brauer F. The kermackmckendrick epidemic model revisited. Math Biosci. 2005; 198(2):119–31.View ArticlePubMedGoogle Scholar
 Lenhart S, Workman JT. Optimal control applied to biological models. New York: Chapman & Hall/CRC; 2007.Google Scholar
 Rowthorn RE, Laxminarayan R, Gilligan CA. Optimal control of epidemics in metapopulations. J R Soc Interface. 2009; 6(41):1135–44.PubMed CentralView ArticlePubMedGoogle Scholar
 GonzálezParra PA, Lee S, Velazquez L, CastilloChavez C. A note on the use of optimal control on a discrete time model of influenza dynamics. Math Biosci Eng. 2011; 8:183–97.View ArticlePubMedGoogle Scholar
 Lee S, Chowell G, CastilloChavez C. Optimal control for pandemic influenza: the role of limited antiviral treatment and isolation. J Theor Biol. 2010; 265(2):136–50.View ArticlePubMedGoogle Scholar
 Lee S, Morales R, CastilloChavez C. A note on the use of influenza vaccination strategies when supply is limited. Math Biosci Eng. 2011; 8(1):171–82.View ArticlePubMedGoogle Scholar
 Lee S, Golinski M, Chowell G. Modeling optimal agespecific vaccination strategies against pandemic influenza. Bull Math Biol. 2012; 74(4):958–80.View ArticlePubMedGoogle Scholar
 Lee J, Kim J, Kwon HD. Optimal control of an influenza model with seasonal forcing and agedependent transmission rates. J Theor Biol. 2013; 317:310–20.View ArticlePubMedGoogle Scholar
 Fleming WH, Rishel RW. Deterministic and stochastic optimal control. New York: Springer Verlag; 1975.View ArticleGoogle Scholar
 Pontryagin LS, Boltyanskii VG, Gamkrelidze RV. The mathematical theory of optimal processes. New Jersey: Wiley; 1962.Google Scholar