Optimal Control of Heterogeneous Mutating Viruses

Different strains of influenza viruses spread in human populations during every epidemic season. As the size of an infected population increases, the virus can mutate itself and grow in strength. The traditional epidemic SIR model does not capture virus mutations and, hence, the model is not sufficient to study epidemics where the virus mutates at the same time as it spreads. In this work, we establish a novel framework to study the epidemic process with mutations of influenza viruses, which couples the SIR model with replicator dynamics used for describing virus mutations. We formulated an optimal control problem to study the optimal strategies for medical treatment and quarantine decisions. We obtained structural results for the optimal strategies and used numerical examples to corroborate our results.


Introduction
An epidemic of infectious disease occurs when a virus population undergoes genetic mutations or new species of viruses are introduced into the host population, and host immunity to that change in the virus population is suddenly reduced below a certain threshold.Hence, epidemic modeling should take into account not only the population dynamics of the host population, but also those of the viruses.In traditional epidemiological models, differential equations are used to capture the dynamic evolution of different classes of host populations.In particular, susceptible (S) is the class of people who are not infected; infected (I) is the class of people who have the disease; and removed or recovered (R) represents the quarantined or immune population.The commonly used SIR model [1,2] is used to describe the population migrations between these three classes of models.In order to capture the interdependencies between a virus and host populations, we established a system framework that combines the SIR model with evolutionary models that describe virus mutations.
In this work, we studied an influenza epidemic in urban populations.We analyzed the evolutionary model for virus mutations together with SIR models for the evolution of susceptible, infected, and recovered subpopulations.Over the time, individuals from these subpopulations randomly interact with each other and change their state.We consider the epidemic process as a dynamic process of changing states from susceptible individuals to the infected and finally to the recovered.The influenza epidemic is a fast-spreading process, involving a large part of the total population.Hence, one key concern is protecting a population during the annual epidemic season.There exist methods of prevention that reduce the sickness rate to protect populations, and medical measures (pharmacological products, quarantine policies, etc.) that reduce the number of the infected in the population.
Another aspect of the influenza epidemic is that different strains of influenza viruses can spread in the population during each epidemic season.As it was shown in previous studies, influenza viruses have many modifications and can quickly transform every season.Thus, by including a simulation of the mutation process into the SIR model, we can preliminarily forecast treatment strategies.In this work, we used evolutionary-game tools to describe the evolution of the influenza virus because this approach allows us to estimate the survivability of virus modifications and to determine the strongest modification.Thereby, we can use treatment strategies more effectively.
More specifically, in our work, we assume that the virus has two types with different strains and fitness functions.Both types of viruses spread in urban populations, and, hence, during an epidemic, different parts of population are infected.In our model, we split infected subpopulations into two subgroups and considered a modified SIR model.Therefore, the epidemic process in urban population depends on changes in the virus population.
In our work, we formulated the SIR model under the mechanism of a virus mutation that affects a human population and considered the minimization of treatment costs and the number of infected in both subpopulations to reduce the speed of epidemics.This complex problem is formulated as an optimal-control problem, and the virus mutation is described by replicator dynamics.
The paper is organized as follows.In Section 2, we discuss related work to our model.Section 3 presents the evolutionary model of viruses.In Section 4, we establish the epidemic model for urban populations.In Section 5, we use Pontryagin's maximum principle to find the optimal control and present structural results of the optimal-control problem.In Section 6, we used a numerical simulation to illustrate our results.The paper is concluded in Section 7.

Related Works
The recent literature has seen a surge of interest in using optimal control and game-theoretic methods to study disease control of influenza for public health [3][4][5].This research problem can be traced back to Reference [6], where an SIR type of mathematical framework was proposed to study epidemic spread in a homogeneous population.It provides a deterministic dynamical system model as the mean field approximation of the underlying stochastic evolution of the host subpopulations.In Reference [7], a control problem was formulated for a model of the carrier-borne epidemic model, and it was shown that the optimal-control effort switches at maximum once between the maximum and the minimum control effort.In Reference [8], many variants of optimal-control models of SIR epidemics were investigated for the application of medical vaccination and health-promotion campaigns.In References [9,10], a dynamic SIR epidemic model was used to identify the optimal vaccination-policy mixes for flu season.
Epidemic models have also been used in computer science and engineering to describe the temporal evolution of worm propagation in computer networks [11][12][13].Methods such as stochastic system analysis and optimal control have been applied to provide insights on the spread of epidemics, as well as disease-control policies for protecting a population with quarantine and removal.In Reference [14], sequential hypothesis testing was adopted to detect worm-epidemic propagation over the Internet under a stochastic density-dependent Markov jump process propagation model.In References [15,16], optimal-control methods were used to study the class of epidemic models in mobile wireless networks, and Pontryagin's maximum principle was used to quantify the damage that malware could inflict on the network by deploying optimum decision rules.
Game-theoretic approaches were also used to analyze the strategic interactions between malicious worms and the system under attack.In Reference [17], a virus-protection game was proposed based on two-state epidemic models for N nodes and the characterization of the equilibrium focus on the steady state of the response.In Reference [18], static-and dynamic-game frameworks were used to design equilibrium-revocation strategies for defending sensor networks from node-capturing and cloning attacks.It was shown that the N + 1 non-zero-sum differential game framework was equivalent to a zero-sum differential game between a team of N attackers and the system.In References [19,20] the SIR model was combined with a game-theoretical approach to define the optimal medical approach: vaccination or treatment of seasonal influenza.
In many modern studies, different modifications of SIR models were used to estimate the behavior of infectious diseases, such as Ebola and Severe acute respiratory syndrome (SARS) [21].
Different from past works, this paper considers a coupled system framework composed of the SIR epidemic model and the evolutionary dynamic model for virus mutations.This framework is motivated by the fact that the epidemic spread of a virus can facilitate virus mutations, strengthening its virulence, which, in turn, expedites the spread and worsens epidemics.The SIR epidemic model with virus mutations can capture this complex interactions between virus and host, and allows us to explain more complex phenomena through analysis.

Evolutionary Model of Virus Mutations
Infectious diseases, such as influenza and SARS, are urgent public health problems in modern urban environments.Influenza spreads faster, especially in large urban populations, and affects people's lifestyle and working facilities.The occurrence of epidemics depends on many factors, such as the size of the human population, and the virus strain and virulence, and it has become important to use effective tools to reduce their influence on human populations [22][23][24].A mathematical model of virus infections in a population can be used to study the factors that influence epidemic growth to improve existing treatments and evaluate new effective prevention measures and treatment.In earlier research in the literature, it was shown that, during epidemic season, the influenza virus can mutate, and several types of the influenza virus circulate in the human population.Different mutations of the influenza virus affect human beings with different intensities, and epidemics evolve depending on virus type and intensity.Hence, the evolution of virus mutations should be taken into account when the SIR model is used to model influenza epidemics.
In this work, we coupled together two dynamic processes, i.e., the evolution of virus mutations and the epidemic process in human populations, as one dynamical system.The corresponding scheme of the system is illustrated in Figure 1.At the top level, two different types of the influenza virus compete to infect the host to continue their life cycles, and thereby lead to the spread of epidemics in the human population.The total population would contain several infected subpopulations that correspond to different virus types.On the bottom level, the human population is divided into subpopulations: the susceptible (S), the infected (I), and the recovered (R).
Spreading of viruses can be controlled using prevention measures such as medical treatment or isolation of infected individuals of population.Thus, at this level, we considered the SIR model with those control parameters.
At the top level of coupled dynamical system, we used evolutionary dynamics to describe the mutations within the virus population.We first describe the interactions between two virus types using an evolutionary-game model, for which we defined pure strategies, fitness, and rule of changes in a population.In the game, two types of viruses compete for human organisms, and, depending on the strength of the virus, one type can survive or vanish from the virus population.
We assume that the virus has two types or strains, denoted by V 1 and V 2 , and without loss of generality, we assume that V 2 dominates virus V 1 .The fitness of virus type V i in the population is J i (V i , V j ), i, j = 1, 2, which depends on the survivability of the virus among its infected population (e.g., human beings).The life cycle of viruses requires a host organism, and occupation of such an organism leads to energy costs.Hence, virus payoff J i is composed of two components: one is the utility of occupying the host organism, and the other is the cost, i.e., energy costs, Utility of occupation b i = b(I i ), as well as energy costs C I = C(I i ), are the functions of population state I i and, hence, mixed strategies x 1 , x 2 ∈ [0, 1] over set (V 1 , V 2 ) are also dependent on the population states.Here, a mixed strategy is defined as the fraction of corresponding virus types circulating in a population.Thereby, the state of the virus population is defined as value x(t) = (x 1 (t), x 2 (t)), where x 1 (t) + x 2 (t) = 1, ∀t ∈ [0, T].Transition rule: Scheme describes the reaction to two heterogeneous viruses circulating in a population.We assumed that the epidemic process could be controlled using treatment or quarantine methods.These measures can be considered as the control parameters in the system, which are used to reduce the size of the infected population and terminate the epidemic process.
Depending on virus strength, the number of people infected by different virus types is different.We used two flows of epidemic processes in a human population to describe our model and the population.We used I 1 to denote the population state for the subpopulation infected by the Type 1 virus, and I 2 is the state for the subpopulation infected by the Type 2 virus.Both viruses spread over the entire human population, and the interactions between two viruses when attacking the same human organism are described with the following four scenarios: The virus incurs energy costs C 1 with probability 1/2 if it cannot occupy a host organism, and achieves utility of b 1 with probability 1/2 if it succeeds in occupation.
The above four competition cases between the two types of viruses are summarized in the following matrix representation: ) .
According to the theory of evolutionary games [25], we compared the payoff of the i-th pure strategy with the average payoff of the total population.If the difference was positive, then the number of individuals using this pure strategy increases, or otherwise decreases.The average payoff of population u : R 2 → R is defined as where e i ∈ R 2 , i = 1, 2, is a vector with i-th element being one and 0 otherwise, indicating the i-th pure strategy; u is a continuous function.The payoff of the i-th pure strategy is defined by u(e i , x) = e i Ax, where A is the payoff matrix of the current symmetric game, and x i (t) is a fraction of virus V i [26].
To describe the evolution of the virus, we need to use a system of differential equations.In the current work, we focused on replicator dynamics [27] to describe changes of states in virus populations.
where ε ∈ R + is the time-scaling factor.Since the mutation process in a virus population and the epidemic process in a human population may develop with different speeds (e.g., virus can mutate faster than spread in human populations), ε can be used to describe such difference in the time scale of the two dynamics.
The stationary state of System of differential Equation ( 1) leads to a symmetric Nash equilibrium [26,27].Therefore, depending on parameters b i and C i , the game has two asymmetric Nash equilibria (1, 0), (0, 1) corresponding to the strategies in which all populations are types V 1 and V 2 , respectively; and one symmetric Nash equilibrium (x, x), where x = (x 1 , x 2 ), The symmetric case is more interesting since both virus types have influence on the human population.By using the evolutionary scenario of virus behaviors, we can predict which type of virus will prevail in a human population.Hence, this knowledge can correct treatment therapy and, at the same time, decrease the total cost of epidemics.

Epidemic Process for an Urban Population
Consider a total urban population of size N with two types of viruses circulating in the population during epidemic season.The human population is divided into four groups: the susceptible, the infected by virus V 1 , the infected by virus V 2 , and the recovered.The susceptible are a subpopulation of human beings that are not infected by viruses but could be infected by one or both types of viruses, and they do not have immunity to them.We assume that, in human populations, two types of viruses coexist at the same time.Human organisms can be occupied by both types of viruses, and, hence, this leads to competition between viruses for the host.Depending on virus strength, we observed that the number of people infected by virus i or by virus j can be different, and people infected by virus V 1 or V 2 belong to the infected subpopulation.The recovered subpopulation consists of people recovered from being infected.The mixing of urban populations allows viruses to spread quickly, and each person in the population is assumed to be in contact with others with equal probability.Hence, when an infected individual interacts with a susceptible one, the virus spreading is made possible.A virus with higher virulence, by our assumption, has a higher probability of success in spreading when interaction occurs between an infected individual and a susceptible individual.

Epidemic Dynamics
We modeled a virus spreading in an urban population using the epidemiological SIR model, where a system of differential equations is used to describe the fraction of the urban population as a function of time.Then, at time t, n S (t), n I 1 (t), n I 2 (t), n R (t) correspond to fractions of the population who are susceptible, infected by virus V 1 , infected by virus V 2 , and recovered, respectively; for all t, condition as portions of the susceptible, the infected and the recovered in the population.At the beginning of epidemic process t = 0, most people in the population belong to the Susceptible subpopulation, a small group in the total population is infected, and the other people are in the recovered subpopulation.Hence, the initial states are: 0 < S(0 We extended the simple SIR model introduced by References [6,28] to describe the situation with two virus types: where δ i are the infection rates for virus V i , i = 1, 2, We can interpret self-recovery rate σ 1 for virus V 1 or σ 2 for virus V 2 , which shows the probability that infected nodes from subgroups I 1 or I 2 are recovered from the infection without causing any costs to our system.Infection rate is defined as the product of infection transmissibility, i.e., the probability of infection being transmitted during contact: where δ i (0) = δ i 0 I i (0).Here δ i 0 , i = 1, 2 determine the virulence of the particular virus or the ability to infect a susceptible host.
In this work, changes in the virus population influence the parameters of the SIR model; therefore, the number of infected is a function of corresponding virus subpopulation x i , I i (t) = I i (x i , t), i = 1, 2. We let I i (x i , t) be linear and δ i (t) take the following form: Then, the SIR model can be rewritten as follows: In the model above, the infection rate is captured in the evolution of the mutation process in the urban population.Medical treatment or quarantine isolation reduces the number of infected individuals in the urban population.These prevention measures can be interpreted as control parameters in the system defined as u = (u 1 , u 2 ); here, u i are fractions of the infected that are quarantined or under intensive medical treatment: 0

Objective Function
In this work, we minimize overall cost in time interval [0, T].At any given t, the following costs exist in the system: f 1 (I 1 (t)), f 2 (I 2 (t)), infected costs; g(R(t)), benefit rate; h 1 (u 1 (t)), h 2 (u 2 (t)), costs for medical treatments (i.e., quarantine or removal) that help reduce the epidemic spread; k I 1 , k I 1 , k R , costs and benefit for invective and recovered at the end of the epidemic.Here, functions f i (I i ) are non-decreasing and twice-differentiable, convex functions, i.e., f i (0) = 0, f i (I i ) > 0 for I i > 0, i = 1, 2., g(R) is non-decreasing and differentiable function and g(0) = 0, h i (u i (t)) is twice-differentiable and increasing function in u i (t) such as h i (0) = 0, h i (u i (t)) > 0, i = 1, 2, when u i (t) > 0. We show the structure of the optimal control strategies for the general case of costs functions.Hence, in particular cases this conditions will be satisfied for any functions with the same properties.
The cost for the aggregated system is given by and the optimal-control problem is to minimize the cost, i.e., min {u 1 ,u 2 } J.To simplify the analysis, we consider the case where

Optimal Control of Epidemics
We used Pontryagin's maximum principle [29] to find optimal control u(t) = (u 1 (t), u 2 (t)) to the problem described in Section 4 above.Define associated Hamiltonian H and adjoint functions λ S , λ I 1 , λ I 2 , λ R as follows: Here, we used condition We constructed an adjoint system as follows: λS with transversality conditions given by: According to Pontryagin's maximum principle, there exist continuous and piecewise continuously differentiable costate functions λ i that, at every point t ∈ [0, T], where u 1 and u 2 are continuous, satisfy System (7) and Equation (8).In addition, we have:

Structure of Optimal Control
Based on previous research [12,16,29], in this subsection we show that optimal control u(t) = (u 1 (t), u 2 (t)) has the following structural results: Proposition 1.The following statements hold for the optimal-control problem described in Section 4:
• When h i (•) are strictly convex functions, then there exist time moments t 0 , t 1 , 0 < t 0 < t 1 < T such that: The proof of Proposition 1 requires auxiliary Lemma 1 and it is discussed in detail in Section 5.2.Before stating Lemma 1, we define functions φ i as follows: Rewrite the Hamiltonian in terms of function φ and we obtain: For any admissible control u 1 , u 2 and according to Equation ( 9), for all t ∈ [0, T]: then, we obtain We observe that min Since u 1 = u 2 = 0 are admissible control, using Equation (11), we obtain To prove Proposition 1, we consider the following auxiliary lemma: Lemma 1. Functions φ i , i = 1, 2 are decreasing functions of t, for t ∈ [t 0 , T], t 0 ≥ t ≥ 0, while Proof of Lemma 1.The state and costate functions are differentiable functions; then, φ i are also differentiable functions at each time t, t ∈ [0, T] at which functions u 1 , u 2 are continuous.We have to show that φi < 0 at each time t ∈ [t 0 , T], t 0 ≥ t ≥ 0. Consider function φ 1 given by and, likewise, φ 2 as follows: Here, f 1 (I 1 ) ≥ 0, f 2 (I 2 ) ≥ 0, g (R) ≥ 0, δ i ≥ 0, I 1 , I 2 , S, R ≥ 0, then, the right-hand side of Expressions ( 15) and ( 16) are negative if conditions δ 1 SI 1 − σ 1 I 1 ≥ u 1 and δ 2 SI 2 − σ 2 I 2 ≥ u 2 are satisfied; otherwise, functions φ i are increasing.Term δ i SI i − σ i I i ≥ u i ≥ 0, i = 1, 2 can be interpreted as a condition for the beginning of the epidemic (see Reference [6]).The proof of Lemma 1 is completed.

Proof of Proposition 1
In this subsection, we prove Proposition 1 under two cases of cost functions h i (u i ), i = 1, 2.
From Equation ( 10), we have that optimal u i satisfies , where u i are any admissible control, u i ∈ [0, 1].If u i = 1, then switching functions are satisfied φ i ≥ h i (1), and, if Lemma 1 suggests that φ i are decreasing functions; then, there can be, at most, one time moment t ∈ [0, T] at which φ i (t) = h i (u i (t)).Moreover, if such moment exists, for example, t 1 , then φ i (t) < h i (1) on 0 ≤ t < t 1 and φ i (t) ≥ h i (1) on t 1 ≤ t < T.Then, Control (17) is satisfied.

h i (•) Are
Function φ i , h i , u i is continuous at all t ∈ [0, T].In this case, h i is strictly convex and h i is a strictly increasing function, so h (0) < h (1).Thus, there exist points t 0 , t 1 , 0 < t 0 < t 1 < T, so that Conditions ( 14) and ( 18) are satisfied; according to that, φ i is a decreasing function.In a time interval where functions φ i are increasing and Conditions (18) are rewritten.There may exist such time interval [0, t 0 ) that u i = 0 and φ i ≤ h i (0), i = 1, 2; then, for time interval [t 0 , T], Conditions (18) continue to be satisfied.
Using auxiliary Lemma 2, we completed the proof of Proposition 1. From Lemma 1, we need to check that multipliers (λ 15) and ( 16) are non-negative.Lemma 2. For all 0 ≤ t ≤ T, we have (λ Property 1.Let w(t) be a continuous and piecewise differential function of t.Let w(t 1 ) = L and w(t) > L for all t ∈ (t 1 , . . ., t 0 ].Then ẇ(t Property 2. For any convex and differentiable function y(x), which is 0 at x = 0, y (x)x − y(x) ≥ 0 for all x ≥ 0.
Proof: We first prove the case for t = T and then for t < T.
If δ 1 > 0 and δ 2 > 0, then the system of O.D.E. is autonomous, and hence the Hamiltonian and the control do not have the dependence of variable independent t.
Case III.In this case, we prove that (λ R (t) − λ I 1 (t)) > 0 in a similar way.
Together with Lemma 1, proof of Lemma 2 completes the proof of Proposition 1.

Quadratic Cost Functions
In this subsection, we consider a particular case of Proposition 1, where cost functions h i (u), i = 1, 2 are quadratic, i.e., The quadratic function is strictly convex if coefficient a 0 > 0, and we can apply the same arguments as in Case II of Lemma 1.Consider ∂ ∂x (h i (x) − φ i x) | x=x 1 = 0 from Proposition 1, where h i (u) is defined as in Function (28); then, we obtain: Hence, we arrive at the following form of optimal control: Functions φ i , h i , u i are continuous at all t ∈ [0, T].In this case, h i is strictly convex and h i is a strictly increasing function, so h (0) < h (1).Thus, there exist points t 0 , t 1 , 0 < t 0 < t 1 < T, such that conditions ( 29) and ( 14) are satisfied, and φ i is a decreasing function.If Conditions ( 14) are broken, then there may exist such time interval [0, t 0 ) that u i = 0 and φ i ≤ h i (0), i = 1, 2; then, for time interval [t 0 , T], Conditions (29) continue to be satisfied.

Numerical Simulation
In this section, we present numerical simulations to corroborate our results.Consider a city with population N = 100,000 people, where two viruses, of different strengths, spread (δ 1 0 = 0.4 and δ 2 0 = 0.5).At time moment t = 0, half of the population are susceptible to the infection, i.e., S(0) = 0.5.The initial Recovered subpopulation is R(0) = 0. We set 18% of population as infected by virus V 1 , and 32% of the population are infected by virus V 2 , i.e., I 1 (0) = 0.18 and I 2 (0) = 0.32.The epidemic lasted for 45 days.We assumed that, in this experiment, the self-recovery rates were equal to σ 1 = 0.001 and σ 2 = 0.002.During the epidemic, people from the infected population incur infected costs and, hence, we defined cost functions as f 1 (I 1 ) = 5I 1 , f 2 (I 2 ) = 6I 2 ; the benefit rate from the recovered subpopulation was g(R) = 0.1R.In Section 5, we have shown that medical-treatment cost functions to describe the value of treatment can be chosen as concave or strictly convex.In our simulations, or concave cost functions, we used h 1 (u 1 ) = 7u 1 , and h 2 (u 2 ) = 9u 2 ; for convex cost functions, we used h 1 (u 1 ) = 15u 2  1 and h 2 (u 2 ) = 10u 2 2 .The costs here were measured in the same monetary units (m.u.), which could be US dollars, Chinese RMB, or Euros, depending on the context.
The first experiment shows the behavior of the system in the absence of treatment (control).As a control strategy, we used the medical treatment of the infected host, and the convex form of cost functions, i.e., h 1 (u 1 ) = 15u 2  1 and h 2 (u 2 ) = 10u 2 2 .After simulations, we obtained that the maximum amount of replicas in the case without applying treatment were I 1 (t 1 ) = 0.3031 at t 1 = 11.25, and I 2 (t 2 ) = 0.6742 at t 2 = 10.25 (Figure 2).From Experiment 1, we can see that, in the absence of treatment in the end of epidemic period, our population had the following distribution of Infected hosts: 29% are infected by first type of virus and 64% are infected by the second type of virus.As we have self-recovery rate σ 1 = 0.001 and σ 2 = 0.002, a fraction of the Recovered is R(T) = 0.07.In the case when we applied the control, the fraction of infected hosts was zero at the end of interval T = 45, and the fraction of the recovered was 92%.There were also some susceptible nodes (S(T) = 0.07) that were not affected by the epidemic.Comparing the controlled and uncontrolled case, the aggregated system costs were: J = 1004.43 in the uncontrolled model and J = 268.35 in the controlled model (Figure 3).Experiment 2 represents the SIR model with virus mutation.Infection rates δ i (t), i = 1, 2 indicate the speed of viruses spreading in the population.In our work, it means that we take the competition between viruses for the host into account.Here, we supposed that the stronger virus captures more hosts.The strength of the virus depends on the infection rates, which change under a mutation process.In this experiment, infection rates changed by the formula δ Evolutionary dynamics allow us to estimate how viruses evolve in an urban population, and to calculate the Nash equilibrium.The Nash equilibrium shows the expected proportions of infected populations after a certain amount of time.This estimation can be useful when building forecasts and tactics to prevent an epidemic.We can allocate available resources in such way that the aggregated costs of the system are minimized.Using the SIR model, we illustrate how the system develops under various types of control.
Here, we have utility of occupation functions as b 1 (I 1 ) = I 1 and b 2 (I 2 ) = I 2 , and energy-cost functions are C 1 (I 1 ) = 5I 1 and C 2 (I 2 ) = 2I 2 .According to these data, we observed that the Nash equilibrium was equal to (0.24; 0.76) in the case when we did not apply any treatment.Figure 4 (left) shows the behavior of the system in the uncontrolled case, where the maximum population of V 1 was I 1 (t 1 ) = 0.2165 at t 1 = 11.25 and the maximum population of V 2 was I 2 (t 2 ) = 0.7746 at t 2 = 8.5.By applying optimal treatment strategies, we observed that the maximal values for both types of viruses were equal to their initial state (I 1 (0) = 0.12 and I 2 (0) = 0.38).The structure of the optimal control in the case with the virus mutation is shown in Figure 5 (left).Comparison of aggregated costs is presented on Figure 5 (right).The aggregated system cost was J = 1048.5 in the uncontrolled model and J = 294.59 in the controlled model.Experiment 3 shows the SIR model with the virus mutation (Figure 6).Utility-of-occupation functions were b 1 (I 1 ) = 2I 1 and b 2 (I 2 ) = 2I 2 , and energy-cost functions were C 1 (I 1 ) = 5I 1 and C 2 (I 2 ) = 3I 2 .The expected proportions of infected subpopulations were equal to the Nash equilibrium (0.17; 0.83), according to simulations we received that the difference between parameters b i and C i influences the equilibrium state.From Equation (1), we have that the fraction of the stronger virus is higher if the difference is smaller.Figure 7 demonstrates the change of the infection rates of Experiments 2 and 3.Under the application of optimal control, we have that the maximum values for both types of viruses are equal to their initial state (I 1 (0) = 0.2 and I 2 (0) = 0.3).
The fourth experiment describes the case when the costs for medical treatments are concave functions.In this case, we have proven in Section 5 that the optimal control has a "bang-bang" structure, i.e., Figure 8 (right).The control for the first type of virus is turned off on the seventh day, and, for the second type, on the ninth day.At the end of interval T = 45, the proportion of the recovered hosts was 60%, while the remaining population was still susceptible to infection.The comparisons of the aggregated costs are presented on Figure 9 (left).Aggregated system cost was J = 1035.4 in the uncontrolled model and J = 289.42 in the controlled model.
Figure 9 (right) demonstrates the evolution of viruses over time in human populations.Here, we can see that there are three stationary states corresponding to three Nash equilibria, and that the convergence of solution trajectories of ODE (1) depends on the initial states.).Vertical axis shows the aggregated costs at time moment t in m.u.Right: Simplex of mixed strategies of the symmetric bimatrix game for modeling virus mutations.In this numerical example, the set of Nash equilibrium was found to be {(1, 0), (0, 1), (0.5, 0.5)}.

Conclusions and Discussion
In this paper, we studied an epidemic model that takes into account the evolutionary dynamics of virus mutations.Different from other studies, we did not consider a group of latent individuals who are already infected but they do not have any clinical symptoms.Many studies have shown that this subgroup can intensively influence epidemic dynamics, but, in the current study, we concentrated on the extended SIR model, including the virus-mutation process.However, we believe that this additional group could enrich our future research.Classical SIR epidemic dynamics are strongly coupled with the replicator dynamics of the virus.We formulated an optimal-control problem in which we sought to minimize the costs of preventing outbreaks by describing the structure of the optimal treatment and quarantine strategies against infection from two different types of viruses.Using Pontryagin's maximum principle, we showed that, depending on the structure of the cost functions, optimal control has a threshold structure.We corroborated our results with numerical examples, observing different switching times for the control strategies under models with and without virus mutations.As future work, we would extend this study to multiple types of viruses and apply different evolutionary dynamics to model the process of virus mutations, including imitative dynamics and best-response dynamics.

Figure 1 .
Figure 1.Transition rule: Scheme describes the reaction to two heterogeneous viruses circulating in a population.We assumed that the epidemic process could be controlled using treatment or quarantine methods.These measures can be considered as the control parameters in the system, which are used to reduce the size of the infected population and terminate the epidemic process.

Figure 2 .
Figure 2. Experiment 1. Left: SIR model without virus mutation and without applying control.Initial states are I 1 (0) = 0.18, I 2 (0) = 0.32, maximum values are I 1max = 0.3031 and I 2max = 0.6742.Epidemic peaks a reached at 11th and 10th days.Right: SIR model without virus mutation with application of control.Vertical axes show the fractions of the subpopulations.

Figure 3 .
Figure 3. Experiment 1. Left: Optimal control in SIR model without virus mutation; cost functions are convex h i (u i ).Vertical axis shows the amount of applied control.Right: Comparison of aggregated costs of SIR model without virus mutation (Controlled model: J = 268.35;uncontrolled model: J = 1004.43).Vertical axis shows the aggregated costs at time moment t in m.u.

Figure 4 .
Figure 4. Experiment 2. Left: SIR model with virus mutation and without applying control.Maximum values are I 1max = 0.2165, I 2max = 0.7746.Epidemic peaks reached at the 11th and 8th days.Right: SIR model with virus mutation and application of control.Functions h i (u i ) were convex.Vertical axes show the fractions of the subpopulations.

Figure 5 .
Figure 5. Experiment 2. Left: Optimal control in the SIR model with virus mutation and convex cost functions h i (u i ).Vertical axis shows the amount of applied control.Right: Comparison of aggregated costs of the SIR model with virus mutation (Controlled model J = 294.59,uncontrolled model J = 1048.5).Vertical axis shows the aggregated costs at time moment t in m.u.

Figure 6 .
Figure 6.Experiment 3. Left: SIR model with virus mutation and application of control.Functions h i (u i ) are convex.Vertical axis shows the fractions of the subpopulations.Right: Optimal control in SIR model with virus mutation and convex-cost functions h i (u i ).Vertical axis shows the amount of applied control at time moment t.

Figure 8 .
Figure 8. Experiment 4. Left: SIR model with virus mutation and application of control.Functions h i (u i ) are concave.Vertical axis shows the fractions of the subpopulations.Right: Optimal control in SIR model with virus mutation and concave-cost functions h i (u i ).Control was switched off at the seventh day for V 1 , and at the ninth day for V 2 .Vertical axis shows the amount of applied control at time moment t.

Figure 9 .
Figure 9. Left: Experiment 4. Comparison of aggregated costs of the SIR model with virus mutation (Controlled model J = 289.42,uncontrolled model J = 1035.4).Vertical axis shows the aggregated costs at time moment t in m.u.Right: Simplex of mixed strategies of the symmetric bimatrix game for modeling virus mutations.In this numerical example, the set of Nash equilibrium was found to be {(1, 0), (0, 1), (0.5, 0.5)}.