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 [1]. 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 [2], 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 [3]. 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% of people living with HIV were found to have resistance mutations to the drugs used for prophylaxis [3]. 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 to more resistant variants.

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 of Lythgoe et al. [4] that considers the possibility of a person infected with a virus of type i can transmit a virus of type j at a time-dependent rate β(t). While this approach presents a detailed model of the within-host viral evolution, it requires a substantial amount of data that 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 approached analytically.

Complex multi-strain models have been proposed for influenza [5,6] and dengue [7] that focus on both cross-reactivity among circulating strains and coinfection [8] rather than the emergence of new strains within infected hosts. Much 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 [9]; analysis of a multi-strain (actually two-strain) disease with environmental transmission, no mutations [10]; bifurcation analysis of a number of (rather simple) multi-strain epidemiological models without mutation [11].

Further information about different approaches to modeling the evolutionary and population-based dynamics of multi-strain pathogens as well as the description of the problems that arise in this connection can be found in [12,13]

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 coinfection. 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.

Description of the Models and Their Structural Properties

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. This type of neutral genotypic evolution is the basis of both the molecular epidemiology and phylodynamics of viral infections. In order to account for virus variability, the whole set of viruses is divided into n types (strains), V, i = 1, … , n (in the following, we will use the words virus strain and virus type interchangeably). This space defines the evolutionary domain of the model (i.e., only these pre-specified strains can exists); however, the analytical results we present are valid for any number of strains. We assume that all infections are founded by a single virus of type i. During the acute infection stage, we assume that the patient’s viral population consists of only the founding-type virus, while during the chronic stage, the original virus is allowed to mutate thus producing new strains according to the within-host evolutionary model. Therefore, we assume that each chronically infected individual’s viral population contains a distribution of viral strains dependent on both the founding viral type and the time from infection.

To model the population-level process of disease propagation we assume an SI model with two stages of disease: an initial acute stage followed by a longer chronic stage. 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, … , I, I, … , I, T, S] is the (2n + 2)-dimensional state vector, i ranges from 1 to n, and the respective forces of infection are defined as

Table 1.

Model parameters. Parameters indicated with an asterisk are used only in the extended model Equation (4).

State Variable	Range	Description
IAi	[0, 1]	Fraction of acutely infected individuals infected by the virus of type i.
ICi	[0, 1]	Fraction of chronically infected individuals infected by the virus of type i.
S	[0, 1]	Fraction of susceptible individuals
T	[0, 1]	Fraction of patients involved in treatment
Ti*	[0, 1]	Fraction of patients infected by the virus of type i that are involved in treatment
P*	[0, 1]	Fraction of patients involved in prophylaxis
Parameter	Range	Description
uT		Rate at which chronically infected are enrolled into treatment (controlled parameter)
uP*		Rate at which susceptible individuals are enrolled into prophylaxis (controlled parameter)
γ		Inverse duration of the acute phase
μ		Mortality rate
αij	[0, 1]	Fraction of type i viruses in the viral population of an individual initially infected by the type j virus.
βA, βC		Transmissibility rates of acute and chronically infected individuals.
ξ		Proportionality coefficient of the transmissibility in acute and chronic stages
ζi*		Failure rate of treatment for individuals infected by the virus of type i
δ*		Failure rate of prophylaxis
ψi*	[0, 1]	The level of protection against the virus strain i, which is conferred by prophylaxis; ψi = 1 corresponds to full protection

In Equation (2), α ∈ [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 for all j = 1, … , n. Furthermore, we assume that α ≠ 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 β and β are the transmission rates of acute and chronically infected individuals. In the baseline model, the viruses are phenotypically homogeneous, therefore 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.

Furthermore, we assume a homogeneous contact structure, i.e., we assume that a susceptible individual can have contact with any infected person with the same probability. The only distinction is made between acutely and chronically infected individuals that are assumed to have different transmissibility coefficients β and β. Specifically, we assume that β and β differ by a proportionality coefficient ξ: β = ξβ (we shall still occasionally write β if it makes the notation more straightforward). With this, expression Equation (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 efficacy of the treatment program depends on the viral strain. That is, the treatment program fails with certain probability, which varies depending on the virus strain, causing the treated individuals to thus revert back to active chronic infection.

Virus strains have different levels of contagiousness.

The efficacy of prophylactic measures depends on the viral strain. While on prophylaxis, an individual acquires protection against the virus depending on the specific viral strain.

To account for different failure rates of treatment we divide the group of people on treatment into n compartments T, 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 ψ ∈ [0, 1] with ψ = 1 corresponding to full protection. Thus we have the following model: where ζ ≥ 0 is the failure rate associated with the ith control, δ is the inverse duration of prophylaxis, and u, resp., u 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 β and β, 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., β = ξβ. The proportionality coefficient ξ does not depend on the virus type i and is determined by the number of virions in the blood, which is assumed to be the same for all virus types. The forces of infection ϕ(X) are defined as

Note that setting either ζ = 0 or ψ = 0 or β = β 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 (the use of notation E instead of I for the identity matrix is common in German and Russian mathematical texts (Germ., Einheitsmatrix); here we use it to avoid confusing notation I with the letter I used for infected compartments). The sizes of the respective matrices are indicated as subscripts. A single subscript, for example, as in E, denotes a square [n × n] matrix of respective type. Furthermore, I and I 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. Necessary facts about special classes of matrices that will be used throughout the text are presented in Appendix A.

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 Equation (1) since the extended model has the same properties and can be readily analyzed along the same lines.

Non-negativity of the solutions. The Equation (1) can be written as

The vector-valued functions Ψ(X) and Ψ (X) are essentially non-negative, i.e., for all j = 1, … , m, m = 2n + 2, it holds that (resp., ) for any such that (see [14] for details). This implies that solutions of Equation (1) are non-negative. That is to say, for any non-negative initial condition and any non-negative control u the solution of Equation (1) belongs to for all t ≥ 0.

Boundedness of solutions. Observe that the m-simplex Δ, formed as the convex hull of m unit vectors e, j = 1, … , m, is invariant with respect to Equation (1): where . This result follows immediately from the fact that the states X 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 both the baseline and extended 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 Equation (1) and β for the extended model Equation (4).

Basic Reproduction Number for the Baseline Model

The system Equation (1) has a unique disease-free equilibrium (DFE) X = [0, … , 0, 1]. To analyze the stability property of the system Equation (1) at the DFE we compute the controlled basic reproduction number R0 using the classical next-generation method [15] (see [16] for an extension of the method that takes into account the action of a control).

Theorem 1. For any choice of parameters α ≥ 0 such that Σ α = 1 and α ≠ 0 for all i, j = 1, … , n, the controlled basic reproduction number of the system Equation (1) is given by

Proof. See Appendix B. □

Note that the α 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 α 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) R1 that quantify the efficiency of sufficiently small controls in reducing the value of R0, [16]. In particular, the controlled basic reproduction number is expanded as where , , and is a high-order term, which is proportional in magnitude to the square of the control u. Before proceeding with further analysis, we define the notion of efficiency of a control.

Definition 1. Let the uncontrolled basic reproduction number be larger than 1, i.e., R0 (0) > 1. A control u is said to be

Locally efficient if the respective sensitivity parameter is negative, i.e., R1 < 0;

(Globally) efficient if there exists a non-negative value u* such that R0 (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 is unconditionally locally efficient. However, an unconditionally locally efficient control may fail to reach the stated goal of eliminating the infection, i.e., reducing R0 below 1. The following result illustrates that.

Lemma 1. The control u is globally efficient if β

Proof. This result can be easily checked by observing the expression for R (u) in Equation (8) and noting that the second summand vanishes as u tends to infinity. □

Remark 1. Note that the condition θ < 1, where θ = 1/(γ + μ) denotes the average duration of the acute stage.

The result of Lemma 1 implies that the control u is only conditionally globally efficient. That is, it can be used to completely eliminate the infection only if the transmissibility β satisfies Equation (10).

Basic Reproduction Number for the Extended Model

In contrast to the baseline case, the disease free equilibrium for the modified model Equation (4) is shifted due to the action of the control u. So, we have where and . Local stability of the DFE Equation (11) is determined by . 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 = diag(β, … , β), Ψ = diag(ψ1, … , ψ), and Z = diag(ζ1, … , ζ) to denote the diagonal matrices of transmissibility rates, protection factors and treatment failure rates.

Theorem 2. The controlled basic reproduction number of the system Equation (4) is given by where , , , , and Δ(u) = (Z + (μ + u) E)−1 (Z + μ E).

Proof. See Appendix B. □

Note that the basic reproduction number of the extended system is a product of two terms: the first one closely resembles as in Equation (9), while the second term is the spectral radius of the product of two matrices, where the first one depends only on u and the second one depends only on u.

Before we proceed with the analysis, we formulate an important result on stochastic matrices that we need to obtain further results.

Lemma 2. Let Σ be a non-negative, column stochastic matrix. Then for any α ∈ [0, 1], the convex combination Σ = α E + (1 – α)Σ is a column stochastic matrix as well. Furthermore, the left and right dominant eigenvectors of Σ coincide with those of Σ.

Proof. Consider the ith column of the matrix Σ, i = 1, … , n. Summing its components and using the fact that Σ is column stochastic we get α1 + (1 – α)1 = 1. This implies that Σ is column stochastic as well. Further, let v be the right dominant eigenvector of Σ, i.e., Σv = v. We have hence, it is the right dominant eigenvector of Σ as well. The case of a left dominant eigenvector is shown analogously. □

Sensitivity analysis. We begin this paragraph by writing down an expansion of .

Theorem 3. Let A be irreducible and let w0 and v0 be the left and the right dominant eigenvectors of , corresponding to and normalized such that . The controlled basic reproduction number can be written as where , , and .

Proof. See Appendix B. □

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 = β E and . Then, according to Lemma 2, the vectors w0 and v0 coincide with those of A. In particular, we have w0 = 1[ since the matrix A is column stochastic. The respective coefficients turn into , , and . That is, we can write

Obviously, both controls are unconditionally locally efficient. We can also observe that the control u is locally more efficient than u if it holds that

The inequality Equation (14) implies that the control u decreases the basic reproduction number to a larger extent, as the respective sensitivity coefficient is larger in absolute value than . Obviously, we have that u is locally more efficient if the opposite holds true. The inequality Equation (14) can be interpreted as follows. Note that τ = 1/(ζ + μ) and π = 1/(δ + μ) are the average duration of being either on treatment or on prophylaxis and recall that . Then we can write Equation (14) 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 while they are in the acute stage of infection. Therefore, the duration and contagiousness of the acute stage of infection are potentially strong determinants 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 v0 are proportional to the stationary distribution of different strains of the virus in the baseline model (see Section 4 for more details on that). Thus, we can interpret the sensitivity parameters and 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 Equation (13) as

As above, we say that u is locally more efficient than u if where β is the contagiousness of the ith strain during the acute stage and θ is the average duration of the acute stage (cf. Remark 1).

Similarly to the previous case, we interpret the expression in front of τ 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 βθ < 1. In contrast to the previous case, the components w0v0 do not have that clear interpretation. However, their behavior is pretty close to that of v0.

Finally, we observe that for sufficiently large controls u and u we have which yields the result that agrees with the result of Lemma 1.

Lemma 3. The controls u and u 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. The general model is considered in Section 5 that is devoted to the numerical simulations.

We begin by stating a general result on the endemic equilibrium.

Theorem 4. Let A be an irreducible non-negative column stochastic matrix such that all diagonal elements are non-zero. Then the endemic equilibrium for the system Equation (1) exists and is unique if R0 > 1. Let, furthermore, v⊤ = [v1, … , v] be the right normalized eigenvector of A corresponding to the dominant eigenvalue of A and satisfying . The components of the endemic equilibrium state are given by

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 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 β.

Before we proceed to the next result we recall that α 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. 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 from Equation (16), we can write

However, since v is the dominant eigenvector of A, it holds that , whence the result follows. □

If the matrix A is reducible, the results of Theorem 4 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 = 1, … , k are irreducible matrices and asterisks denote arbitrary non-negative matrices.

Theorem 5. Let A be a reducible non-negative matrix with non-zero diagonal elements such that it can be transformed into the normal form Equation (18) by an appropriate simultaneous permutation of rows and columns. Then A has at most k unit eigenvalues. Furthermore, let v = {v1, … , v} 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 vector belongs to the linear hull of vectors from .

Theorem 5 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 G1 and G2. Reducibility implies that the mutations between these groups are either not possible at all, G1 ↮ G2 or possible only in one direction, G1 ← G2, but G1 ↛ G2 (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 n2 parameters α from the observations made at the endemic equilibrium. The reason for this is that the equilibrium values depend on the n components of the vector v (see Equation (16)) 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 A1 in Appendix A 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 π of the original virus while the remaining strains of the virus are distributed uniformly. This means that the matrix A has the following form:

To make a connection to the standard form of the matrix A Equation (6) we note that with such parametrization we have α = π and for all j ≠ i.

The matrix Equation (20) 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 ratios of infected in different groups, which turn out to have a simple form:

Note that the condition Σ v = 1 implies that there are only (n – 1) independent equations. Thus, one parameter π can be set to an arbitrary value within the range [0, 1]. Let, for instance, π be used as a free parameter. In this case, all remaining probabilities can be expressed in terms of π and v:

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:

Matrix A in Equation (21) is non-negative, irreducible and acyclic. Hence, the Perron–Frobenius theorem applies as well. Quite remarkably, the respective expressions do not change that much compared to the previous case. Setting v and assuming that π can be freely chosen we get

Numerical Simulation for Different Scenarios and Illustration of Results

To illustrate our results we will present several scenarios that are aimed at illustrating different aspects of the considered problem. As a testbed for our numerical analysis we consider a model with four virus strains. On the one hand, this model is complex enough to illustrate different interesting features of the studied model, but on the other hand, it can be easily visualized and analyzed. We believe that such setup might be a reasonable use case of the framework where we have a limited number of putative strains that are characterized by well-defined phenotypes. A particularly important fact is that our results are valid for any number of strains.

First consider the baseline model (no prophylaxis, single treatment, and uniform contagiousness). Later on we will extend this model along several directions. We take the following values of the parameters: μ = 0.025, γ = 3 (i.e., the acute phase takes about four months); ξ = 5 (during the acute phase an individual is five times more contagious as in the chronic one); u = 0.4 (it takes 2.5 years on average until the treatment begins). The baseline transmissibility rate β was chosen such that R0(u = 0.4, u = 0) ≈ 1.2: β = 0.25. The matrix A is assumed to have the form Equation (21). When choosing the values of the probabilities π, we imposed the following conditions on the endemic distribution of the different virus strains:

Case 1. v = 3, j = 1, 2, 3. Assuming that π4 = 0.25 one can compute the remaining probabilities using Equation (22): π1 ≈ 0.97, π2 ≈ 0.83, and π3 = 0.5. Finally, the endemic frequencies are [v1, v2, v3, v4] = [0.675, 0.225, 0.075, 0.025].

Case 2. v = 7, j = 1, 2, 3. Similarly to the previous case, we fix π4 = 0.25 and compute the remaining probabilities π1 = 0.9985, π2 = 0.9796, and π3 = 0.8571. The respective endemic frequencies are [v1, v2, v3, v4] = [0.8575, 0.1225, 0.0175, 0.0025].

Now we proceed with a detailed qualitative analysis of the two described models.

Controlled Basic Reproduction Number

Equation (12) that gives the controlled basic reproduction number for the extended model does not have an immediately apparent intuitive interpretation, however, it allows us to measure the influence of within-host evolution on the controllability of pathogens by therapy and prophylaxis. Figure 1 illustrates how Equation (12) can be used to measure the effects of different levels of resistance to prophylactic interventions on the extent of controls needed to bring a pathogen to sub-critical levels. In the “No resistance” case, we assume that resistance to prophylaxis cannot be evolved (ψ = 1, i = 1, … , 4), and we see a clear synergistic interaction between the therapy and prophylactic controls; the level of therapeutic control can be reduced by more than half by even a modest investment in the prophylactic control. In the “Resistance, high cost” we assume that ψ4 = 0, that is, a very rare strain is resistant to the prophylactic control. We refer to this scenario as ‘high cost’ because the strain quickly reverts to a wild-type variant suggesting that the mutations involved have a high evolutionary cost and, in the absence of the prophylactic control, the mutations defining this strain are slightly deleterious. In this scenario, the benefit of increasing the prophylactic control saturates when the selection for the resistant phenotype balances out the deleterious effects of the resistance mutations and there is no further population-level benefit to prophylaxis. Minor variants that confer resistance, even when that resistance is associated with a non-trivial evolutionary cost, can have a major effect on the control properties of the system. In the “Resistance, low cost” scenario we assume that the evolutionary cost of mutations that confer resistance is lower (i.e., once the mutation(s) occur, they tend to be lost at a slower rate), which we parameterize by setting ψ1 = 0. In this scenario, even a very small level of prophylaxis leads to the resistant strain becoming the dominant strain and has only a small effect on R0. The results in Case 1 and Case 2 were nearly identical suggesting that the control properties of the system are robust to some level of variation in the underlay within-host evolutionary dynamics. Note that we consider the uniform transmissibility/uniformly efficient treatment case, hence there is a little phenotypic variability in the model.

Figure 1.

The panel shows the values of R0 (u, u) as a function of two controls for two cases described above. The red color corresponds to the case R0 ≤ 1. We assume a uniform rate of transmission, i.e., β = β = 0.3 for all i = 1, … , 4 and fully efficient treatment, i.e., ζ = 0, i = 1, … , 4. Remaining parameters are: ξ = 5; γ = 3; μ = 0.025; and δ = 0. The subfigures forming the panel correspond to the following values of prophylaxis efficiency coefficients: left, ψ = [1, 1, 1, 1]; central, ψ = [1, 1, 1, 0]; right, ψ = [, 1, 1, 1].

Endemic Distribution with Variable Transmissibility

To see how varying transmissibility influences the endemic distribution we fix the transmissibilities of the first three strains to be equal to β, while the transmissibility of the fourth strain is β = aβ, where a changes from 0.7 to 2. This particular choice is dictated by the wish to have a well expressed example of the variation of endemic frequencies. Obviously, if we varied the transmissibility of a different strain, the result would be the same, but less expressed.

The resulting relative endemic frequencies for both cases are shown in Figure 2. The relative frequencies were computed as . Note that at a = 1 the endemic distribution coincides with the baseline one. 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.

Figure 2.

The relative endemic distribution of infected individuals for different values of the transmissibility rate of the 4th strain, parametrized with a: β = aβ. The values at a = 1 (marked by a red dashed line) correspond to the baseline case, where all transmissibility rates are equal. Subfigures (a) and (b) correspond to different values of mutation probabilities π.

Endemic Distribution with Variable Prophylaxis Effects

To study the effect of prophylaxis on the endemic distribution of different strains we assume that prophylaxis confers a full protection against the three first strains, while providing no protection against the last strain. That is to say, the matrix Ψ has the form

Furthermore, we assume that the prophylaxis program takes three months, i.e., δ = 4. We change the rate at which people are recruited to prophylaxis and study its effect on the relative frequencies of respective strains. The results are presented in Figure 3.

Figure 3.

The relative endemic distribution of infected individuals for different values of u. Subfigures (a) and (b) correspond to different values of mutation probabilities π.

We see that as u grows, the fraction of the last, resistant strain grows as well. In general, the frequencies of strains tends to a more uniform distribution. On the other hand, as u grows, the total fraction of infected individuals decreases and approaches zero for u = 0.88 (Case 1) or u = 0.8 (Case 2). It is worth noting that even for a sufficiently large value of u = 0.75, the fraction of the total population being on prophylaxis does not exceed 15%. This can be explained by a relatively fast turnover: one cycle of prophylaxis lasts three months, after which the individual returns to the group of susceptible.

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 [17]. 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.

Endemic Distribution with Imperfect Treatment

To study the situation with imperfect treatment we follow the same route and assume that the treatment is absolutely efficient for all strains except the last one. We vary the rate at which the treatment fails and compute the endemic distribution of strains as shown in Figure 4. We also assume that there is no prophylaxis. The result turns out to be quite surprising: not only the endemic frequencies reshuffle, but also the total proportion of infected individuals increases dramatically, see Figure 5.

Figure 4.

The relative endemic distribution of infected individuals for different values of ζ4. Subfigures (a) and (b) correspond to different values of mutation probabilities π.

Figure 5.

The total fraction of infected as a function of ζ4.

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 and 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 showed 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 within the population. Our presented framework can be used as a basic foundation for studying the complex interventions such as imperfect vaccines, antibody-based prophylaxis, and new small-molecule therapeutics for a variety of chronic infections such as HIV, herpes, and HPV. It is possible that this framework may be useful for studying even reality short lived infections such as COVID-19. Early data suggests that the virus accumulates mutations over the course of a single infection [18] and that some of those mutations may affect the viruses’ contagiousness [19]. Likewise, this model could be used to study the implications of strain specific vaccine effects or how different vaccines could change the genetic landscape of SARS-CoV-2.