Within-Host Phenotypic Evolution and the Population-Level Control of Chronic Viral Infections by Treatment and Prophylaxis

Chronic viral infections can persist for decades spanning thousands of viral generations, leading to a highly diverse population of viruses with its own complex evolutionary history. We propose an expandable mathematical framework for understanding how the emergence of genetic and phenotypic diversity affects the population-level control of those infections by both non-curative treatment and chemo-prophylactic measures. Our frameworks allows both neutral and phenotypic evolution, and we consider the specific evolution of contagiousness, resistance to therapy, and efficacy of prophylaxis. We compute both the controlled and uncontrolled, population-level basic reproduction number accounting for the within-host evolutionary process where new phenotypes emerge and are lost in infected persons, which we also extend to include both treatment and prophylactic control efforts. We used these results to discuss the conditions under which the relative efficacy of prophylactic versus therapeutic methods of control are superior. Finally, we give expressions for the endemic equilibrium of these models for certain constrained versions of the within-host evolutionary model providing a potential method for estimating within-host evolutionary parameters from population-level genetic sequence data.


Introduction
Pathogens that lead to persistent chronic infection in people must mitigate both the innate and adaptive immune systems. Strategies for evading the innate immune system are complex including direct subversion of host signaling pathways [12]. Pathogens such as HIV avoid the adaptive immune system by simply evolving new phenotypes faster than the host immune system can adapt leading to a rapid co-evolutionary race. Because HIV has a short generation time and generates a massive number of new viral particles each day [16], this evolutionary race creates a large potential for the emergence of new viral strains and phenotypes. This rapid evolutionary process is one of the many reasons that HIV is exceedingly difficult to treat. Mutations that evolve in a single host are also known to be transmitted. In 2014-2016, 6 out of 11 countries looking for the presence of pre-treatment drug resistance (i.e. presence of a drug resistant phenotype in persons unexposed to the drug) reported greater than 10% of new infections were resistant to one or more non-nucleoside reverse-transcriptase inhibitor, which is related to both treatment failure and death [15]. The rapid emergence of new viral phenotypes within infected persons is not just a clinical problem, it is an epidemiological problem. Chemo-prophylactic measures that focus on protecting uninfected persons using similar drugs to those used for treating infected persons are not immune to evolutionary derailment. In King Country, Washington, 0.5% people living with HIV were found to have resistance mutations to the drugs used for prophylaxis [15]. However, with the emergence of chemo-and bio-prophylactic agents (i.e. anti-HIV antibodies for prevention), we must consider the possibility that population-level administration of these agents can shift the ever evolving landscape of chronic viral infections. This paper is motivated by a need for mathematical models that integrate within-host genetic diversity and phenotypic evolution with epidemiological dynamics and consider the effects of joint therapeutic and prophylactic controls. We also attempted to balance the complexity of the model to be usable as a data analysis tool with the desire to understand the mathematical and statistical properties of the model using analytical methods. Our model accounts for within-host evolution among multiple phenotypes characterized by variable contagiousness, resistance to prophylactic measures, and resistance to therapeutic measures. Our framework allows for new phenotypes to emerge in chronic infection that can be both transmitted and possibly lost in later hosts. We consider both the epidemiological and evolutionary effects of both therapy for infected persons and chemo-prophylaxis-type measures for uninfected persons.
There has been a number of results devoted to the analysis of different aspects of the evolutionary and epidemiological dynamics of a multi-strain pathogen. While there is a wide spectrum of different models covering different aspects of virus/immune system evolution and their interaction, most developed models are too complex to be analyzed analytically thus, their analysis restricts to carrying out and analyzing the results of numerical simulations. Our model is related to the approach to Lythgoe et al. [11] that considers the possibility of a person infected with virus of type i can transmit a virus of type j at a time-dependent rate β ij (t). While this approach presents a detailed model of the within-host viral evolution, it requires a substantial amount of data which is not readily available: virus reproduction, mutation and death rates. Furthermore, since we need to take into account the duration of the infection at the time when transmission occurs, the system dynamics is governed by integro-differential equations which are difficult to deal with. On the other hand, such a detailed approach turns out to be an overkill as the total pool of infected contain the individuals at different stages of disease and hence, the transmission rates undergo a sort of averaging over the whole set of infected. Therefore, we employ a simpler formalism in which we treat the virus evolutionary dynamics in a more coarse grained fashion. This allows us to balance our mutual goals of a sufficiently complex model that can still be aproached analytically.
Complex multi-strain models have been proposed for influenza [9,14] and dengue [1] that focus on both cross-reactivity among circulating strains and coinfection [18] rather than the emergence of new strains within infected hosts. Much of this work of this work is often based on complex models that are intended to explain specific biological phenomena that are too complex to be understood by applied analysis methods. On the other hand, there has been a number of papers devoted to the analytic analysis of certain aspects of multi-strain virus dynamics. However, most of the papers either deal with rather restricted setups or study only certain aspects of the system dynamics. We mention stability analysis of within-host multi-strain virus dynamics with mutations [3]; analysis of a multi-strain (actually two-strain) disease with environmental transmission, no mutations [2]; bifurcation analysis of a number of (rather simple) multi-strain epidemiological models without mutation [8].
Further information about different approaches to modeling the evolutionary and population-based dy-namics of multi-strain pathogens as well as the description of the problems that arise in this connection can be found in [10,21] It should be noted that most research effort aimed at studying the dynamics of multi-strain viruses does not take into account the possibility of within-host mutations and concentrate on modeling different immune system responses in reaction to re-infection or co-infection. In contrast to that, we are more concerned with the effect of mutations on both the within-host and population-level distribution of viral strains and on how both the emergence and loss of phenotypes within infected persons alters the population-level control of chronic viral infections.

A baseline model of a chronic multi-strain virus infection
In the baseline version of the model we consider the within-host evolution and transmission of distinct strains that have the same phenotype. Most observed mutations at the nucleotide level will not alter a pathogen's phenotype, but they can be used in combination to discern unique lineages of transmission. In order to account for virus variability, the whole set of viruses is divided into n groups (strains), V i , i = 1, . . . , n. An individual gets infected by a virus of particular type i. During the acute stage the patient's viral population consists of viruses of the same type, while during the chronic stage the original virus mutates thus producing strains of all types. Therefore, we assume that each chronically infected individual's viral population contains viruses of all n types in the proportion that depends on the initial type of the virus.
To model the process of disease propagation we employ an SI model and assume that there are two stages of the disease: the acute and the chronic one. Furthermore, we extend the set of state variables to include the individuals enrolled into treatment. In doing so we assume that the treatment is completely efficient and the patients are fully compliant with the treatment.
When writing the differential equations of the model we assume the inflow to be equal to the outflow hence, the total population size remains constant. Therefore, we write the model equations for the fractions of the respective cohorts in the total population. This implies, in particular, that the sum of all the states is equal to 1. We have the following set of DEs (state variables and parameters defined in Table 1): where X = [I A1 , . . . , I An , I C1 , . . . , I Cn , T, S] is the (2n + 2)-dimensional state vector, i ranges from 1 to n, and the respective forces of infection are defined as In (2), α ij ∈ [0, 1] denotes the average fraction of type i viruses in the viral population of an individual initially infected by the type j virus. It should thus hold that n i=1 α ij = 1 for all j = 1, . . . , n. Furthermore, we assume that α ii = 0 for all i = 1, . . . , n. This means that the viral population of an individual infected with the type i virus always contains a non-zero amount of the corresponding strain. Parameters β A and β C are the transmissibility rates of acute and chronically infected individuals. In this simple setting we assume that the probability that a susceptible individual contracts a disease depends only on the disease stage of the infected contact, but not on the type of virus. That is to say, a susceptible can be equally well infected by any virus. In the following, we will assume that β A and β C differ by a proportionality coefficient ξ: β A = ξβ C 1 . With this, expression (2) turns into

A generalized model with differentially effective control, variable transmissibility and prophylaxis
We generalize the baseline model by allowing different strains to have different phenotypes by relaxing the model assumptions along the following lines: • The treatment program does not ensure complete suppression of viral replication. That is to say, the treatment program fails with certain probability, which varies depending on the virus strain.
• Virus strains have different transmissibility.
• In addition to the treatment, we consider the effect of chemical or biological prophylaxis. While on prophylaxis, an individual acquires protection against the virus, the extent of which depends on the virus strain.
To account for different failure rates of treatment we divide the group of people on treatment into n compartments T i , where i corresponds to the virus strain. Furthermore, we add a cohort of people receiving prophylaxis, denoted by P . While on prophylaxis, the individuals acquire variable protection against different virus strains denoted by ψ i ∈ [0, 1] with ψ i = 1 corresponding to full protection. Thus we have the following model:İ where ζ i ≥ 0 is the failure rate associated with the ith control, δ is the inverse duration of prophylaxis, and u T , resp., u P are the rates at which people are administered to either treatment or prophylaxis. To account for variable transmissibility of different virus strains we define a set of transmissibility rates β Ai and β Ci , i = 1, . . . , n. Similarly to the baseline case, the transmissibility rates for the corresponding acute and chronic stages are assumed to be proportional, i.e., β Ai = ξβ Ci . The forces of infection φ i (X) are defined as Note that setting either ζ i = 0 or ψ i = 0 or β Ci = β C for all i = 1, . . . , n, we obtain different variations of the baseline model.
Notation. We let 0, 1, and E denote the matrices of zeros, ones, and the identity matrix 2 . The sizes of the respective matrices are indicated as subscripts. A single subscript, e.g., as in E n , denotes a square [n × n] matrix of respective type. Furthermore, I A and I C denote the column vectors of respective variables and A denotes the matrix of α's: Note that A is a non-negative, column stochastic matrix, i.e., all its columns sum to 1. All parameters and variables used in the model are listed in Table 1. Note that all quantities used are assumed to take on non-negative values and the index i always runs from 1 to n.

Structural analysis
In this subsection we consider only the baseline model (1) since the extended model has the same properties and can be readily analyzed along the same lines.
Non-negativity of the solutions. The equations (1) can be written as The vector-valued functions Ψ(X) and Ψ u (X) are essentially non-negative, i.e., for all j = 1, . . . , m, m = 2n + 2, it holds that Ψ j (X) ≥ 0 (resp., Ψ u j (X) ≥ 0) for anyX ∈ R m ≥0 such thatX j = 0 (see [4] for details). This implies that solutions of (1) are non-negative. That is to say, for any non-negative initial condition X(0) = X 0 ∈ R m ≥0 and any non-negative control u T the solution of (1) belongs to R m ≥0 for all t ≥ 0.

Boundedness of solutions.
Observe that the m-simplex ∆ m , formed as the convex hull of m unit vectors e j , j = 1, . . . , m, is invariant with respect to (1): . This result follows immediately from the fact that the states X i represent the fractions of the respective groups within the total population and hence sum to 1.

Local analysis at a disease-free equilibrium
Below, we will compute the basic reproduction number for two considered models and present a number of related results. To distinguish between the basic reproduction numbers related to different models we will add a superscript denoting the particular model: α for the baseline model (1) and β for the extended model (4).

Basic reproduction number for the baseline model
The system (1) has a unique disease-free equilibrium (DFE) X DF E = [0, . . . , 0, 1]. To analyze the stability property of the system (1) at the DFE we compute the controlled basic reproduction number R 0 using the classical next-generation method [17] (see [5] for an extension of the method that takes into account the action of a control).
Theorem 3.1. For any choice of parameters α ij ≥ 0 such that i α ij = 1 and α ii = 0 for all i, j = 1, . . . , n, the controlled basic reproduction number of the system (1) is given by Proof. See Appendix B.
Note that the α ij values do not affect the basic reproduction number, which makes sense in this context as mutation from one strain into another does not imply any change in a relevant phenotype such as contagiousness or resistance to therapy. In this context, a different strain simply carries a distinct mutation (or pattern of mutations) that makes it identifiable from other strains. However, understanding the distribution of strains with the same phenotype is an important aspect of molecular epidemiology, which is dependent on the specific α ij values. This relationship between within-host mutations and endemic equilibrium of infection types is discussed below.
Sensitivity analysis. When devising an intervention strategy, the main question to be answered is whether the proposed treatment or prophylaxis scheme is capable of eliminating the infection, i.e., driving the basic reproduction number below 1. To address this issue we introduce the sensitivity parameter(s) R 1 that quantify the efficiency of sufficiently small controls in reducing the value of R 0 , [5]. In particular, the controlled basic reproduction number R α 0 (u T ) is expanded as where . Before proceeding with the further analysis, we define the notion of efficiency of a control.
Definition 3.1. Let the uncontrolled basic reproduction number be larger than 1, i.e., R 0 (0) > 1. A control u is said to be 1. locally efficient if the respective sensitivity parameter is negative, i.e., R u 1 < 0; 2. (globally) efficient if there exists a non-negative value u * such that R 0 (u * ) = 1.
Furthermore, we say that a control is unconditionally locally (globally) efficient if 1. (2.) holds for all admissible values of parameters. Otherwise the control is said to be conditionally efficient.
We can immediately observe that u T is unconditionally locally efficient. However, an unconditionally locally efficient control may fail to reach the stated goal of eliminating the infection, i.e., reducing R 0 below 1. The following result illustrates that.
Proof. This result can be easily checked by observing the expression for R α (u T ) in (7) and noting that the second summand vanishes as u T tends to infinity.
Remark 3.1. Note that the condition (9) can be alternatively rewritten as β A θ A < 1, where θ A = 1/(γ +µ) denotes the average duration of the acute stage.
The result of Lemma 3.2 implies that the control u T is only conditionally globally efficient. That is, it can be used to completely eliminate the infection only if the transmissibility β C satisfies (9).

Basic reproduction number for the extended model
In contrast to the baseline case, the disease free equilibrium for the modified model (4) is shifted due to the action of the control u P . So, we have where S DF E (u P ) = δ+µ δ+µ+uP and P DF E (u P ) = 1 − S DF E (u P ) = uP δ+µ+uP . Local stability of the DFE (10) is determined by R β 0 (u T , u P ). Before we proceed with the analysis, we note that the results to follow will be formulated using matrix notation. In particular, we will write B C = diag(β C1 , . . . , β Cn ), Ψ = diag(ψ 1 , . . . , ψ n ), and Z = diag(ζ 1 , . . . , ζ n ) to denote the diagonal matrices of transmissibility rates, protection factors and treatment failure rates.
Theorem 3.3. The controlled basic reproduction number of the system (4) is given by Note that the basic reproduction number of the extended system is a product of two terms: the first one closely resembles R α 0 as in (8), while the second term is a spectral radius of the product of two matrices, where the first one depends only on u P and the second one depends only on u T .
Before we proceed with the analysis, we formulate an important result on stochastic matrices.
The above result implies that N (0) = 1 γ+ξµ [ξµ E n + γ A] is a column stochastic matrix, whose left and right dominant vectors coincide with those of A. We will thus write N (0) =Ā.
Sensitivity analysis. We begin this paragraph by writing down an expansion of R β 0 (u T , u P ).
Theorem 3.5. Let A be irreducible and let w 0 and v 0 be the right and the left dominant eigenvectors of Q(0)N (0) =B CĀ , corresponding to ρ B CĀ and normalized such that w 0 v 0 = 1. The controlled basic reproduction number R β 0 (u T , u P ) can be written as where This result has a number of important consequences as formulated below. We first consider a slightly simplified setup. Let there be no variability in transmission rates, i.e., B C = β C E n . Then the vectors w 0 and v 0 coincide with those of A. Furthermore, we have w 0 = 1 [n×1] . The respective coefficients turn into Obviously, both controls are unconditionally locally efficient. We can also observe that the control u T is locally more efficient than u P if it holds that Obviously, we have that u P is locally more efficient if the opposite holds true. The inequality (13) can be interpreted as follows. Note that τ i = 1/(ζ i + µ) and π = 1/(δ + µ) are the average duration of being either on treatment or on prophylaxis and recall that w 0 = [1, . . . , 1]. Then we can write (13) as Here, the factor γ/(γ + ξµ) is interpreted as the degree of protection given by the treatment. Note that this number decreases with increasing ξ, i.e., when the acute stage is much more contagious compared to the chronic stage. If ξ = 1, the fraction γ/(γ + µ) merely corresponds to the fraction of people that survive to the chronic stage. Note that this interpretation has to do with the fact that we assume the acute stage is short enough that people will not start treatment conditional on being in the chronic stage of infection and therefore the acute to chronic stage contagiousness of infection is a major determinant of the efficacy of therapy as a population-level control. This assumption is reasonable for diseases like HIV, but may need to be revisited for application to other diseases. Next, we note that the components of the vector v 0 are proportional to the stationary distribution of different strains of the virus in the baseline model (see Sec. 4 for more details on that). Thus, we can interpret the sensitivity parameters R β 1,T and R β 1,P as a sum of products average duration of the medical intervention × protection conferred by the intervention taken with the weights corresponding to the stationary distribution of the virus strains.
Following the same line, one can attempt to compare the efficiency of two controls in the general case. To start with, we write (12) as As above, we say that u T is locally more efficient than u P if Similarly to the previous case, we interpret the expression in front of τ i as the degree of protection given by the treatment to those infected with the i-type virus. Note that a sufficient condition for this expression to be positive is β Ai θ A < 1. In contrast to the previous case, the components w 0i v 0i do not have that clear interpretation. However, their behavior is pretty close to that of v 0i as the numerical simulation presented in Fig. 2 illustrates. Finally, we observe that for sufficiently large controls u T and u P we have which yields the result that agrees with the result of Lemma 3.2.
Lemma 3.6. The controls u T and u P are jointly globally efficient if

Endemic equilibrium
In contrast to the unique disease-free equilibrium, there can be one, many (a continuum), and no endemic equilibria at all. Which case realizes in our system depends on the value of the basic reproduction number and on the structure of the matrix A as will be shown below. For the general case, the endemic equilibrium can be computed using a rather involved semi-analytic procedure and offers a little insight into the structure of the respective equilibrium. Therefore, we restrict ourselves to the baseline model and a particular extension thereof. The general model is considered in Section 5 that is devoted to the numerical simulations.

A baseline model
We begin by stating a general result on the endemic equilibrium.
Proof. See Appendix B.
Note that the only additional property of the matrix A that is involved in this theorem is that A is irreducible. For the definition of irreducibility and further details see Appendix A.
The obtained result can be used to compute a number of derived quantities. For instance, we have that the total prevalence at the endemic equilibrium is equal to and the the ratio of transmissions through acutely infected to the transmission through chronically infected is given by Using the statistical estimations of these two parameters one can recover ξ and β C . Before we proceed to the next result we recall that α ij can be interpreted as the probability of catching a virus of type i through the contact with an individual initially infected by the virus of type j. So, we can make the following observation.

Lemma 4.2.
At the endemic equilibrium, the probability of encountering a chronically infected in the ith category is equal to the probability of catching the type i virus through the contact with a randomly chosen chronically infected individual: Proof. Using the expression for I * Ci , we can write However, since v is the stationary eigenvector of A, it holds that n j=1 α ij v j = v i , whence the result follows.
If the matrix A is reducible, the results of Theorem 4.1 do not apply any longer. However, we can formulate a weaker version of the theorem. First, we note that a reducible matrix can be transformed to the normal form by means of a properly chosen permutation matrix: where A i , i = 1, . . . , k are irreducible matrices and asterisks denote arbitrary non-negative matrices.
Theorem 4.3. Let A be a reducible non-negative matrix with non-zero diagonal elements such that it can be transformed into the normal form (17) by an appropriate simultaneous permutation of rows and columns. Then A has at most k unit eigenvalues. Furthermore, let v = {v 1 , . . . , v q } be the set of normalized eigenvectors corresponding to the unit eigenvalues, q ≤ k. Then the set of endemic equilibria is defined as follows: where the vectorv belongs to the linear hull of vectors from v:v ∈ Span(v).
Theorem 4.3 implies that the set of endemic equilibria can form a linear subspace of the system's state space. The case when the matrix A is reducible corresponds to the situation when there are some particular groups of virus strains, say, two groups G 1 and G 2 . Reducibility implies that the mutations between these groups are either not possible at all, G 1 G 2 or possible only in one direction, G 1 ← G 2 , but G 1 G 2 (or vice versa). Such a setup allows for considering directed patterns of viral evolution. However, this question is beyond the scope of this paper and will be addressed in our future work.
Structure of the matrix A: uniform within host mutations. An important observation that follows from the preceding analysis is that one cannot unambiguously determine all n 2 parameters α ij from the observations made at the endemic equilibrium. The reason for that is that the equilibrium values depend on the n components of the vector v (see (23)) of which only n − 1 values are independent. We thus restrict ourselves to considering one particular structure of the matrix A that can be formulated in terms of only n parameters. More complex structures are possible and can be treated using the same results. In particular, Theorem 5.1 offers a convenient tool for computing the respective dominant eigenvector.
Assume that during the chronic infection stage the viral population of the individual, initially infected with the type i virus contains the fraction α i of the original virus while the remaining strains of the virus are distributed uniformly. This means that the matrix A has the following form: The matrix (19) is positive hence, the Perron-Frobenius theorem applies. There is a unique dominant eigenvalue that is equal to 1, and the components of the dominant eigenvector have the following form: .
The respective expressions for the system states at endemic equilibrium are pretty bulky. However, we can compute the relative ratios of infected in different groups, which turn out to have a simple form: Note that the condition i v i = 1 implies that there are only (n − 1) independent equations. Thus, one parameter α i can be set to an arbitrary value within the range [0, 1]. Let, for instance, α n be used as a free parameter. In this case, all remaining probabilities can be expressed in terms of α n and v i : In the following we will consider a slightly more realistic scenario in which all viruses are ordered according to their genetic similarity and any virus can mutate only to its "neighbors". The respective matrix A has the following form: Quite remarkably, the respective expressions do not change that much compared to the previous case. Setting v i and assuming that α n can be freely chosen we get

Variable transmissibility
In this subsection we consider an extension of the baseline model in which different strains are assumed to have different transmissibility rates. This implies that we consider the model (4) with u P = 0 and ζ i = 0. The respective set of DEs is written (in vector notation) aṡ For this case, the results of Theorems 3.3 and 4.1 can be restated as follows: Theorem 4.4. The basic reproduction number R β 0 (u T ) for the system (21) is given by Let, furthermore, A be an irreducible non-negative column stochastic matrix such that all diagonal elements are non-zero. Let, furthermore, v = [v 1 , . . . , v n ] be the right normalized dominant eigenvector ofB CĀ satisfying n i=1 v i = 1. Then the endemic equilibrium exists and unique if R 0 (u T ) > 1 and the respective components of the endemic equilibrium state are given by Proof. The expression for R 0 (u T ) is obtained from (11) after some algebraic manipulations. The rest of the proof closely follows the proof of Thm. 4.1.
We note thatĀ is a stochastic matrix and has the same right and left dominant eigenvectors as A due to Lemma 3.4. The productB CĀ corresponds to multiplying the rows ofĀ with the respectiveβ Ci . For this expression one cannot find the dominant eigenvector analytically, so we consider some numerical examples to illustrate the influence of variable transmissibility on the distribution of different virus strains. To see how varying transmissibility influences the endemic distribution we fix the transmissibilities of the first 3 strains to be equal to β C , while the transmissibility of the fourth strain is β C4 = aβ C , where a changes from 0.7 to 2. The resulting endemic frequencies are shown in Fig. 1a. Note that at a = 1 the endemic distribution coincides with the baseline one. To make a comparison, we consider a different set of parameters α i . In this case, we assume that these parameters are chosen such that the endemic distribution in the baseline case satisfies v i /v i+1 = 7. The respective values are α 1 = 0.9985, α 2 = 0.9796, α 3 = 0.8571 and α 4 = 0.25 (the last parameter was chosen to be equal to α 4 in the previous case). The trajectories of the endemic frequencies are shown in Fig. 1b. It is interesting to observe that variation in transmissibility of one strain leads to substantial variation in the frequencies of the other strains as facilitated by within-host mutation.
For both cases, Fig. 2 illustrates the dependence of the products w 0i v 0i on the parameter α. One can observe that the behavior of the respective terms (that enter the expression for the the sensitivity terms R β 1,T and R β 1,P ) closely follows that of v 0i . a) b) Figure 2: The change of the weights w 0i v 0i that enter the sensitivity terms R β 1,T and R β 1,P with the parameter a.
Propylaxis. First, we simulate the baseline case with treatment and see that the endemic distribution agrees with the assumed values (Fig. 3a). Next, we turn on prophylaxis at t = 50 and observe how it changes the process (Fig. 3b). The respective control u P was set to 0.1 and the vector of protections against different virus strains is assumed to be ψ = [0.9, 0.6, 0.3, 0.1]. We see that the endemic frequencies change in the way that the frequency of the first strain decreases (since prophylaxis confers almost total protection against this virus) and the frequency of the remaining strains increases. The (normalized) proportion of different strains in the population of acutely infected is [0.616, 0.255, 0.096, 0.033]. For u P = 0.2 this effect becomes even more pronounced. The respective frequencies are [0.547, 0.288, 0.121, 0.043]. On the other hand, as u P grows, the total fraction of infected individuals decreases and asymptotically approaches zero.
Finally, we consider quite a radical setup in which prophylaxis confers total protection against two first virus strains and no protection against the remaining two strains. For u P = 0.2, the respective frequencies are [0.527, 0.277, 0.144, 0.053]. Comparing with the previous result we see that the difference is rather minute. This implies that while imperfect prophylaxis leads to some increase in the frequencies of the viruses that evade it, this increase is rather restricted. The main reason is that when prophylaxis covers a small fraction of the population it does not create sufficient evolutionary pressure, while when it increases it eventually contributes to the complete eradication of the disease. This result is potentially very encouraging as new prevention methods for HIV based on administration of broadly neutralizing antibodies are predicted to have highly differential levels of protection to diverse viral panels [19]. Although further work is needed to explore the potential of selective prophylactic agents to cause strain-level selection in populations in the context of within-host mutation.
Variable transmissibility. The next experiment is aimed at estimating the effect of variable transmissibility of the final distribution of the virus strains. We assumed that the transmissibility of the fourth strain is either 1.25 or 1.5 times larger than the transmissibility of other strains. The results turn out to be quite impressive (see Fig. 4). While for the first case (Fig. 4a)  a complete reshuffle of the strains (Fig. 4b). It should also be noted that an increase in the transmissibility of one strain not only led to an increase of its endemic frequency, but also resulted in an increase in the endemic frequencies of other strains! While the total fraction of acutely infected at the endemic equilibrium is equal to i I Ai = 7.72 · 10 −4 in the first case, it increases to 13.0 · 10 −4 in the second case -almost a doubling! The respective endemic fractions are [0.103, 0.259, 0.35, 0.289].
Imperfect treatment. Now we consider the case when there is a chance of treatment failure. We set the respective rates to be rather small: ζ = [0.025, 0.025, 0.05, 0.1] while keeping all remaining parameters as in the baseline case. We also assume that there is no treatment. The result turns out to be quite surprising: not only the endemic frequencies reshuffle, but also the total proportion of acutely infected increases dramatically, see

Conclusions
In this paper, we described two models of joint evolutionary and epidemiological dynamics of a viral pathogen. While the first, baseline model did not take into account the phenotypic variability of the virus, the extended model addressed the within-host evolution among multiple phenotypes characterized by variable contagiousness, resistance to prophylactic measures, and resistance to therapeutic measures. We presented an analytic expression for the controlled basic reproduction number for both cases and carried out sensitivity analysis of the derived expression to the changes of the control actions. It turned out that the sensitivity coefficients R T 1 and R P 1 have a straightforward interpretation that can be used when assessing the relative efficacy of the controls. Further, we characterized the endemic equilibria for the baseline model and an extension thereof and shown that a sole assumption of variable transmissibility of different virus strains can lead to wide variations in the endemic distribution of the respective strains. Finally, we carried out a numerical analysis aimed at analyzing the effects of phenotypic diversity of virus strains on the population level dynamics and distribution of different virus strains wihtin the population.

Appendix A. Necessary ingredients from matrix algebra
In this appendix, we present some facts about non-negative and stochastic matrices that will be used in the sequel. The interested reader can find a thorough treatment of non-negative matrices, in particular the Perron-Frobenius theory in [13,Ch. 8]. The theory of stochastic matrices within the context of Markov chains is detailed in [7].
Non-negative matrices. A matrix M is said to be non-negative (positive), denoted by M 0 (M 0), if it is element-wise non-negative (positive). The matrix M is said to be reducible if there exists a permutation matrix P such that the conjugated matrix P M P has a block upper-triangular form. Otherwise the matrix M is said to be irreducible. The matrix M is primitive if it is non-negative and there exists k ∈ N such that M k 0. A non-negative irreducible matrix is primitive if at least one diagonal element is non-zero.
For irreducible non-negative matrices, there exists a real eigenvalue, called dominant that is equal to the spectral radius of the matrix. The corresponding left and right dominant eigenvectors are positive. This result follows from the celebrated Perron-Frobenius theorem [13,Ch. 8]. In a reducible case, the above results should be substantially weakened to remain true. In particular, there can be multiple eigenvalues corresponding to the spectral radius r and the respective eigenvectors are merely non-negative, rather than positive. If a non-negative matrix M is reducible, it can be transformed to the normal form, which corresponds to a block upper diagonal matrix such that the diagonal blocks are irreducible.
Stochastic matrices. A non-negative matrix Q is said to be column stochastic (row stochastic) is all its columns (rows) sum to 1. The notions of (ir)reducibility, primitivity and the Perron-Frobenius theoren can be extended to stochastic matrices in a straightforward manner. Below we mention several properties that are specific for stochastic matrices.
A stochastic matrix is typically used to describe the transition structure of a Markov chain. The spectral radius of a stochastic matrix is equal to 1. The respective normalized right eigenvector v is called the stationary distribution of the respective Markov chain, i.e., Qv = v. Here, normalization means that the components of v must sum to 1. If the stochastic matrix is irreducible, then due to the Perron-Frobenius theorem the stationary distribution is unique and component-wise positive. Finally, we present the result on computing the stationary distribution of an irreducible stochastic matrix. This is a version of the Markov chain tree theorem [6], formulated using the results from matrix theory (cf. [20]).
Theorem 5.1. Given an [n × n] irreducible column stochastic matrix Q , the ith element of the right dominant eigenvector of Q is defined as the ith principal minor of the corresponding Laplacian Λ = Q − E n : Proof. We have Qw = w, which is equivalent to (Q − E n )w = Λw = 0. That is, w is the eigenvector corresponding to the zero eigenvalue of Λ or, alternatively, w ∈ ker(Λ). By the Perron-Frobenius theorem, the eigenspace associated with the dominant eigenvalue of Q, and hence, the kernel of Λ is one-dimensional.
By the definition of the adjugate, we have adj(Λ)Λ = Λ adj(Λ) = det(Λ)E n = 0 n . This means, in particular, that the columns of adj(Λ) are linearly dependent and proportional to the stationary distribution v.
By transposing the first expression we get Λ adj(Λ) = 0. The columns of Λ and hence, the rows of Λ sum to 0, which implies that the kernel of Λ is a one-dimensional space spanned by 1 [n×1] . This means that each column of adj(Λ) has the form col i adj( Respectively, each column of adj(Λ) has the form col i (adj(Λ)) = [Λ] 1,1 , . . . , [Λ] n,n . This concludes the proof.

Appendix B. Proofs
Proof (Theorem 3.1). The Jacobian matrix of (1) evaluated at the DFE has the form We observe that the stability of (24) is determined by the eigenvalues of its [2n × 2n] leading submatrix. As a side remark, we mention that this implies that in our case the computation of R 0 requires considering both I A and I C as infected states (cf. the discussion at the end of Sec. 2 in [5]). Following the standard procedure, we split the respective submatrix in two thus obtaining The basic reproductive number is defined as the spectral radius of the product F V −1 , i.e., R 0 = ρ(F V −1 ). Using the standard result on the block matrix inversion we get The product F V −1 is equal to and hence, R 0 is found as the spectral radius of the [n×n] matrix (γ +µ We use the well known facts that if λ is an eigenvalue of the matrix M , then aλ is an eigenvalue of the matrix aM and k + λ is an eigenvalue of the matrix [k E + M ] for any α ∈ R, k ∈ R. This implies that . Finally, since A is a column stochastic matrix, it holds that ρ(A) = 1 and hence, we obtain (7).
Proof (Theorem 3.3). The Jacobian matrix of (4) evaluated at the DFE has the form whereΨ = E n − P DF E Ψ. The Jacobian (26) is a block lower-triangular matrix, whose bottom-right [2 × 2] block is a negated Mmatrix and hence Hurwitz. Thus stability of the DFE is determined by the eigenvalues of the leading [3n×3n] submatrix. As a side remark, we mention that this implies that the T i compartments must be considered as infected. That differs from what we observed in the baseline case and emphasizes the importance of the right choice of the infected compartments (see the discussion at the end of Section 2 in [5]).
Following the standard procedure, we split the respective submatrix in two: A complete expression for the inverse of V is rather bulky. However, we note that since R 0 is computed as the spectral radius of F V −1 , we need to compute only those blocks of the inverse that enter the leading [n × n] submatrix of the product F V −1 . So, we write where we used asterisks to denote the blocks that are not relevant to our problem. The diagonal matrix ∆(u T ) is defined as ∆(u T ) = ( Z + (µ + u T ) E n ) −1 ( Z + µ E n ).
Proof (Theorem 3.5). The proof is similar to the proof of Theorem 3.4 3 in [5] and hence will only be sketched. First, we note that R β 0 = R β 0 (0, 0) =β C (γ+ξµ) (γ+µ)µ ρ B CĀ . Computation of R β 1,T and R β 1,P coefficients boils down to computing partial derivatives of R β 0 (u T , u P ) w.r.t. either u T or u P at u T = u P = 0. Using the same approach as in [5], we get Noting that Q(0) =B C , R(0) =Ā, Q (0) = d duP Q(u P ) The expressions (27a) and (27b) can be further transformed using the fact that w 0 and v 0 are the left and the right eigenvectors ofB CĀ corresponding to the spectral radius of this matrix and expressing ρ B CĀ through the basic reproduction number R β 0 as shown below.
Proof (Theorem 4.1). When computing the endemic equilibrium we let the equilibrium value of S be equal to some (not yet known) value S * . The respective equilibrium values of I Ai and I Ci are found as the solution of the following system of 2n algebraic equations: From the first equation we get Expressing I * C from the second equation of (28) and substituting into the above equation we obtain the expression for I * A : that is I * A belongs to the kernel of the matrix Γ = ((γ+µ)−ξβCS * )(uT+µ) βCγS * E n − A. The matrix A is nonsingular, thus Γ has a non-trivial kernel only if the factor in front of E n is equal to one of the eigenvalues of A. The solution I * A would be then equal (up to a positive factor) to the corresponding eigenvector. However, as follows from the Perron-Frobenius theorem, the only positive eigenvector corresponds to its dominant eigenvalue, which is equal to 1 for a stochastic matrix. Any other eigenvector would contain negative components which contradict the assumption that all system's states are non-negative. This implies that the equilibrium value of S must satisfy the equation ((γ+µ)−ξβCS * )(uT+µ) βCγS * − 1 = 0, whence we obtain The equilibrium solution I * A corresponds to the right dominant eigenvector of the column stochastic matrix A. Since this eigenvector is defined up to the multiplication by a positive scalar, we can further specify it using the following argument. Let i AΣ be the sum of all components of I A , i.e., i AΣ = n j=1 I * Ai . Similarly, we write i CΣ = γ (uT+µ) i AΣ and subsequently, we get T * = γu (uT+µ)µ i AΣ . Since all states sum to 1, we have the following equation Let v be the (normalized) dominant eigenvector of A such that n i=1 v i = 1. Multiplying the components of v with the respective factors we obtain the expressions for I Ai and I Ci .
Finally, we observe that R 0 < 1 implies that the respective components of the endemic equilibrium state turn negative which implies that there is no admissible endemic equilibrium. This concludes the proof.