 Research
 Open Access
 Published:
Modeling optimal treatment strategies in a heterogeneous mixing model
Theoretical Biology and Medical Modelling volume 12, Article number: 28 (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.
Background
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
We consider a twogroup influenza model based on a standard compartmental SITR model. Two additional compartments, a latent class and an asymptomatic class, are included due to the characteristics of influenza. Each class is divided into two subpopulations of sizes N _{1} and N _{2}. For each group i=1,2, we have a susceptible class S _{ i }, a latent class L _{ i }, an infected class with symptoms I _{ i }, an asymptomatic infected class without symptoms A _{ i } and a treated class T _{ i }. A twogroup influenza model involves two different age groups, which are connected by a mixing matrix (p _{ ij }) for i,j=1,2, by allowing for the possibility of subgroups with different activity levels and heterogeneous mixing between these subgroups. A twogroup influenza model with two subgroups can be written as
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).
In this manuscript, we generalize the proportional mixing assumption to the preferred mixing one and we carry out mathematical analysis under different mixing patterns. Suppose that the members of group i make a _{ i } contacts per unit time and that the fraction of contacts made by the members of group i with the members of group j is p _{ ij }, for i, j=1,2; then we have the following:
For our twogroup influenza model, we consider preferred mixing, in which a fraction π _{ i } of each group mixes randomly with its own group and the remaining members mix proportionately. Thus, preferred mixing is given by
where
More details on the preferred mixing formulation can be found in previous studies [18, 20].
The impact of mixing patterns on the contact matrix
Let us investigate the impact of different mixing patterns on the contact matrix. The contact matrix is defined as the product of group activity level, a _{ i } and the group mixing proportions p _{ ij } (i,j=1,2) given in (2).
To illustrate the impact of different mixing patterns, several group mixing fractions are chosen; C _{1} (π _{1}=π _{2}=0), C _{2} (π _{1}=0.25,π _{2}=0.75), C _{3} (π _{1}=π _{2}=0.5), C _{4} (π _{1}=0.75,π _{2}=0.25) and C _{5} (π _{1}=π _{2}=1) using a _{1}=0.5260,a _{2}=0.2670 and N _{1}=1750 and N _{2}=250. Then, we get the following contact matrices:
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).
As a result, the basic reproduction number \(\mathcal {R}_{0}\) with u _{1}=u _{2}=0 is
where \(\Phi _{1}=\left (\frac {\delta (1p)}{\eta _{1}}+\frac {p}{\alpha _{1}} \right), \Phi _{2}=\left (\frac {\delta (1p)}{\eta _{2}}+\frac {p}{\alpha _{2}} \right)\).
The controlled reproduction number \(\mathcal {R}_{c}\) is
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)\).
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
For a onegroup epidemic model, there is a final size relation that makes it possible to calculate the size of the epidemic from the reproduction number [5, 22, 27, 28]. In this section, we establish a final size relation for the twogroup model (1) with u _{1}=u _{2}=0. This relation does not involve the basic reproduction number explicitly but still makes it possible to calculate the size of the epidemic from the model parameters. The final size relation of the model (1) can be obtained as
For the relation between the final size relation and the basic reproduction number, we use the eigenvector v of \(\mathcal {R}_{0}\) as in [21], then
where
The eigenvalue and the eigenvector can be written as
Also, the activity levels can be found in terms of \(\mathcal {R}_{0}\) using (4) and (5),
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 employ optimal control theory to explore the impact of antiviral treatment in situations that mimic 1918like influenza pandemic scenarios. We modify model (1) by incorporating timedependent control functions to measure the effectiveness of groupspecific treatment strategies. Intervention strategies (policies) are modeled by the functions u _{ i }(t)(i=1,2) that externally control the number of treated cases. The objective functional \(\mathcal {F}\) over a finite time interval [ 0,T] is given by the expression:
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).
The optimal control problem is that of finding optimal functions (u _{1} ^{∗}(t),u _{2} ^{∗}(t)) such that
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.
A two point boundary method [29] is employed to find numerical solutions to (7). First, the state system (1) is solved forward with initial conditions. Then, the adjoint system with transversality conditions is solved backward in time. Finally, the optimality condition is updated and whole steps are iterated until convergence is achieved. The baseline parameter values are given in Table 1, which has been taken [20].
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
First, we illustrate the influenza dynamics of (1) in the absence of treatments (u _{1}=u _{2}=0). Figure 1 compares the group specific incidence under three mixing scenarios: proportionate mixing (dotted curve), half mixing (solid curve) and likewithlike mixing (dashed curve). Furthermore, the results are shown under two different values of \(\mathcal {R}_{0}\) (using a moderate value and a higher value on the left and right, respectively). The group 2 incidence is smaller under likewithlike mixing than the one using proportionate mixing, while the group 1 incidence is larger. 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. Hence, this leads to the total incidence or the final epidemic size getting smaller as the mixing becomes likewithlike mixing. Also, it clearly shows more significant differences in the final epidemic size as \({\mathcal {R}}_{0}\) gets larger in the right panels.
Next, the basic reproduction number \({\mathcal {R}}_{0}\) is displayed as a function of group mixing fractions in Fig. 2. The left panel shows the basic reproduction number \(\mathcal {R}_{0}\) using a moderate value of activity levels (\(\mathcal {R}_{0} \in \, [\!1.35, 1.45]\)) while the right one using a higher value of activity levels (\(\mathcal {R}_{0} \in \, [\!2.45, 2.55]\)). Both panels show that the basic reproduction number gets slightly larger as preferred mixing becomes likewithlike mixing (either π _{1} or π _{2} becomes 1). This is consistent with the analytic expression for \(\mathcal {R}_{0}\). Since p _{11} and p _{22} increase and p _{12} and p _{21} decrease as either π _{1} or π _{2} become 1, the basic reproduction number increases as preferred mixing becomes likewithlike mixing. Using the parameter values given here, it is worth mentioning that the effect of the lower activity group mixing fraction (π _{2}) on the values of \({\mathcal {R}}_{0}\) is slightly more significant than the higher activity group mixing fraction (π _{1}). The slope for the axis of π _{2} increases more than the slope for the axis of π _{1} in both panels.
We compute the final epidemic size, which is the number of members of the population who are infected over the course of the epidemic, N−S _{ ∞ } with \(S_{\infty }={\lim }_{t\rightarrow \infty } S(t)\). This can be described in terms of the final attack ratio, (1−S _{ ∞ }/N). In Fig. 3, the final attack ratio is displayed under various mixing scenarios (six different combinations of π _{1} and π _{2}). The left and middle panels show the final attack ratios for group 1 and group 2, respectively, while the right panel shows the total final attack ratio. Note that the range of \({\mathcal {R}}_{0}\) is between 1.72 and 1.82 (xaxis) using a _{1}=0.526 and a _{2}=0.267 as the mixing fractions are varied. The final attack ratio for group 1 gets larger as the mixing becomes likewithlike mixing (π _{1}=π _{2}=1), while it becomes significantly smaller in group 2. Consequently, the total final attack ratio becomes smaller as the mixing becomes likewithlike mixing. Proportionate mixing makes the individuals in group 2 more likely to get infected than likewithlike mixing resulting in a significantly increased the final attack ratio in group 2. This leads to the result that the total final attack ratio follows exactly the same order (i.e. the total final attack ratio becomes smaller as π _{2} becomes to 1 in the right panel). Moreover, all results follow the order of a mixing fraction for the group 2 (π _{2}) whether decreasing or increasing in the final attack ratio (all panels). Therefore, the basic reproduction number and the final attack ratio are not consistent and this reconfirms that the basic reproduction number alone is not sufficient to determine whether preferred mixing increases the final epidemic size or not [21].
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
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. Also, the impact of different levels of transmissibility is investigated by varying the basic reproduction number. Figures 4 and 5 show the results under three different mixing patterns using a moderate value of R _{0}∈ [ 1.73,1.79]. Likewise in the previous section, three mixing patterns are chosen as proportionate mixing, half mixing and likewithlike mixing. In Fig. 4, the proportion of incidence and cumulative incidence in the presence of optimal treatments (red curves) are compared with the results in the absence of treatments (black curves). The results show that there are no outbreak in the presence of treatments and this indicates that groupspecific optimal treatments are effective enough to prevent outbreaks regardless of mixing patterns.
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).
Figures 6 and 7 show the results under three different mixing patterns using a higher value of \(\mathcal {R}_{0} \in \, [\!2.45,2.53]\) and higher activity levels a _{1}=0.742,a _{2}=0.377. Again, Fig. 6 shows the proportion of incidence and cumulative in the presence of optimal treatments (red curves) are compared with the results in the absence of treatments (black curves). Since the basic reproduction number becomes higher, optimal treatments can not stop the outbreaks under all mixing patterns. Figure 7 displays the groupspecific optimal treatment controls and the proportion of incidence. We observe that the time period of treatment gets smaller for group 1 but larger for group 2, than the ones using a moderate \(\mathcal {R}_{0}\). Hence, this results in the cumulative treated cases increasing significantly in both groups. As \(\mathcal {R}_{0}\) becomes higher, the number of infected individuals in group 2 increases dramatically as the mixing becomes proportionate. Note that group 2 (the lower activity group) is more sensitive to mixing patterns.
The impact of mixing patterns on the groupspecific final epidemic size
Figure 8 displays the groupspecific final attack ratio and the total final attack ratio as a function of \({\mathcal {R}}_{0}\) in the absence of treatment under three distinct mixing patterns. The basic reproduction number is increased as the activity levels are increased in the ranges of a _{1}∈ [ 0.1,0.827] and a _{2}∈ [ 0.05,0.42]. The left panel shows that the final attack ratio for group 1 becomes almost indistinguishable regardless of mixing. The final attack ratio for group 2 has the largest value under proportionate mixing (circled), while it becomes smaller as the mixing becomes likewithlike (triangle). This results in the fact that the total (both groups) final attack ratio follows the order of group 2.
Figure 9 presents the comparisons of the final attack ratio and the proportion of cumulative treated as a function of \({\mathcal {R}}_{0}\) in the presence of treatment. It is clear how optimal treatment strategies and the mixing fractions affect the final attack ratio. Similar to the results in the absence of treatment, the final attack ratio for group 1 is almost indistinguishable under all mixing patterns. Again, the final attack ratio for group 2 becomes smaller as the mixing becomes likewithlike. However, these results show that 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}\)). Particularly, for the lower activity group, the reduction is dramatic (triangle in the top middle panel). The cumulative treated results are consistent with the final attack ratio results (more infected, more treatment needed in the bottom panels).
The impact of control parameters
There are two critical control parameters that change the corresponding dynamics greatly. One of them is the control upper bound, b _{ i }, which represents the maximum level of effectiveness for implementing the treatment. The influenza outcomes are dependent on the control upper bound in a straightforward fashion. As the control upper bound is decreased, the magnitude of control decreases for all mixing patterns as shown in Fig. 10. This leads to a longer time of treatment in both groups and interestingly, it results in larger costs (more treated cases) and larger infected cases. This indicates that the higher treatment rate is more effective (less treated cases and less infected as seen in Fig. 5).
The other control parameters are the weight constants or the relative costs of treatments which can play a critical role in the influenza dynamics. As we increase these parameters, the relative cost of control increases and therefore, the magnitude of the optimal controls decreases, resulting in the increase of infected individuals and the final epidemic size in both groups regardless of the mixing pattern in Fig. 11.
The impact of different activity levels and subpopulation sizes
All simulation results so far have been based on the case where group 1 is a higher activity group and a larger population size than the ones for group 2 (a _{1}>a _{2}, N _{1}>N _{2}). Now we investigate the impact of different group activity levels and subpopulation sizes on the optimal treatments and the resulting two group influenza dynamics. There are a total of nine scenarios as we vary group activity levels (a _{1},a _{2}) and subpopulation sizes (N _{1},N _{2}). Only some selected results are presented because the rest of the cases are identical.

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
The basic reproductive number is calculated by using the methodology (the next generation matrix approach) outlined in [26]. Now, let F(x) represent the rate of appearance of new infections. 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) and denote them \({\mathbf {F}}=\left [\frac {\partial F}{\partial \mathbf {x}_{j}}\right ]\) and \(\textbf {V}=\left [\frac {\partial V}{\partial \mathbf {x}_{j}}\right ]\), evaluated at the disease free equilibrium point x ^{∗}, which consists of S _{1}=N _{1}, S _{2}=N _{2} with the rest of them zero.
In F, \(a_{1}p_{12}\frac {S_{1}}{N_{2}}\) is replaced by the balance relation \(\frac {a_{2} p_{1}}{N_{1}}=\frac {a_{1} p_{2}}{N_{2}}\). Then, we get
Also,
The matrix F V ^{−1} has six zero eigenvalues and the remaining two eigenvalues are the roots of the following quadratic equation:
The controlled reproduction \(\mathcal {R}_{c}\) is the largest of these two eigenvalues, which is
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)\).
The basic reproduction number \(\mathcal {R}_{0}\) is \(\mathcal {R}_{c}\) with u _{1}=u _{2}=0:
where \(\Phi _{1}=\left (\frac {\delta (1p)}{\eta _{1}}+\frac {p}{\alpha _{1}} \right), \Phi _{2}=\left (\frac {\delta (1p)}{\eta _{2}}+\frac {p}{\alpha _{2}} \right)\).
Next, we compute the final size relation by introducing the notation g(∞) for \({\lim }_{\textit {t}\rightarrow \infty }g(t)\) and \(\hat {g}\) for \(\int _{0}^{\infty }g(t)dt\) assuming g is a nonnegative integrable function defined for 0≤t<∞.
Also, we have the following:
Using the first equation in (1), we have
Integrating Eq. (9),
Similarly,
Also, integrating Eq. (11),
The final size relation in the absence of treatment (u _{1}=u _{2}=0) can be written as:
Appendix B
The optimal control problem for the twogroup influenza model is formulated to minimize the number of infected individuals for a finite time interval at a minimal cost. We define our objective functional as follows:
Then, we seek an optimal pair (U ^{∗},X ^{∗}) such that
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 in optimal control theory [36]. The necessity conditions of optimal solutions are derived from Pontryagin’s Maximum Principle [37]. This principle converts the systems (1), (6), (7) into the problem of minimizing the Hamiltonian H given by
From this Hamiltonian and Pontryagin’s Maximum Principle [37], we obtain the following theorem:
Theorem.
There exist optimal controls \(u_{1}^{*}(t), u_{2}^{*}(t)\) and corresponding solutions, \(S^{*}_{i}, L^{*}_{i}\), \(I^{*}_{i}\), \(A^{*}_{i}\) and \(T^{*}_{i}\) that minimizes \(\mathcal {F}(u_{1}(t), u_{2}(t))\) over the domain Ω. In order for the above statement to be true, it is necessary that there exist continuous functions \(\lambda _{S_{i}(t)}\), \(\lambda _{L_{i}(t)}\), \(\lambda _{I_{i}(t)}\), \(\lambda _{A_{i}(t)}\) and \(\lambda _{T_{i}(t)}\) for i=1,2 such that
with the transversality conditions,
Furthermore,
Proof.
The existence of optimal controls follows from Corollary 4.1 of [36] since the integrand of J is a convex function of U(t) and the state system satisfies the Lipschitz property with respect to the state variables. The following can be derived from the Pontryagin’s Maximum Principle [37]:
with \(\lambda _{S_{i}}\), \(\lambda _{L_{i}}\), \(\lambda _{I_{i}}\), \(\lambda _{A_{i}}\), \(\lambda _{T_{i}} (i=1,2)\) and evaluated at the optimal controls and corresponding states, which results in the adjoint system (17). The Hamiltonian H is minimized with respect to the controls at the optimal controls, so we differentiate H with respect to u _{1} and u _{2} on the set Ω, giving the following optimality conditions:
Solving for u _{1} ^{∗},u _{2} ^{∗}, then
By using the standard argument for bounds a≤u _{ i }(t)≤b for i=1,2, we have the optimality conditions (15). □
References
 1
Anderson RM, May RM, Anderson B. Infectious Diseases of Humans: Dynamics and Control. Oxford: Oxford university press; 1970.
 2
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.
 3
Heesterbeek JAP, Vol. 5. Mathematical Epidemiology of Infectious Diseases: Model Building, Analysis and Interpretation. New York: John Wiley & Sons; 2000.
 4
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.
 5
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.
 6
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.
 7
Hethcote HW, Van Ark JW. Epidemiological models for heterogeneous populations: proportionate mixing, parameter estimation, and immunization programs. Math Biosci. 1987; 84(1):85–118.
 8
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.
 9
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.
 10
Del Valle SY, Hyman JM, Hethcote HW, Eubank SG. Mixing patterns between age groups in social networks. Soc Netw. 2007; 29(4):539–54.
 11
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.
 12
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.
 13
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.
 14
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.
 15
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.
 16
Longini IM, Halloran ME, Nizam A, Yang Y. Containing pandemic influenza with antiviral agents. Am J Epidemiol. 2004; 159(7):623–33.
 17
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.
 18
Nold A. Heterogeneity in diseasetransmission modeling. Math Biosci. 1980; 52(3):227–40.
 19
Hyman JM, Li J. Behavior changes in sis std models with selective mixing. SIAM J Appl Math. 1997; 57(4):1082–94.
 20
Brauer F. Epidemic models with heterogeneous mixing and treatment. Bull Math Biol. 2008; 70(7):1869–85.
 21
Brauer F. Heterogeneous mixing in epidemic models. Can Appl Math Q. 2012; 20(1):1–13.
 22
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.
 23
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.
 24
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.
 25
Diekmann O HJ, MG R. The construction of nextgeneration matrices for compartmental epidemic models. J R Soc Interface. 2010; 7(47):873–85.
 26
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.
 27
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.
 28
Brauer F. The kermackmckendrick epidemic model revisited. Math Biosci. 2005; 198(2):119–31.
 29
Lenhart S, Workman JT. Optimal control applied to biological models. New York: Chapman & Hall/CRC; 2007.
 30
Rowthorn RE, Laxminarayan R, Gilligan CA. Optimal control of epidemics in metapopulations. J R Soc Interface. 2009; 6(41):1135–44.
 31
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.
 32
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.
 33
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.
 34
Lee S, Golinski M, Chowell G. Modeling optimal agespecific vaccination strategies against pandemic influenza. Bull Math Biol. 2012; 74(4):958–80.
 35
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.
 36
Fleming WH, Rishel RW. Deterministic and stochastic optimal control. New York: Springer Verlag; 1975.
 37
Pontryagin LS, Boltyanskii VG, Gamkrelidze RV. The mathematical theory of optimal processes. New Jersey: Wiley; 1962.
Acknowledgements
This work was supported a grant from Kyung Hee University in 2013 (KHU20130683).
Author information
Affiliations
Corresponding author
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
Both authors have equally contributed to the development of the model, the analysis of the results and the writing of the paper. Both authors read and approved the final manuscript.
Rights and permissions
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.
About this article
Cite this article
Choe, S., Lee, S. Modeling optimal treatment strategies in a heterogeneous mixing model. Theor Biol Med Model 12, 28 (2015). https://doi.org/10.1186/s129760150026x
Received:
Accepted:
Published:
Keywords
 A twogroup influenza model
 Heterogeneous mixing
 Preferred mixing
 Optimal control theory
 Targeted treatment strategies