Within-Host Stochastic Emergence Dynamics of Immune-Escape Mutants

Predicting the emergence of new pathogenic strains is a key goal of evolutionary epidemiology. However, the majority of existing studies have focussed on emergence at the population level, and not within a host. In particular, the coexistence of pre-existing and mutated strains triggers a heightened immune response due to the larger total pathogen population; this feedback can smother mutated strains before they reach an ample size and establish. Here, we extend previous work for measuring emergence probabilities in non-equilibrium populations, to within-host models of acute infections. We create a mathematical model to investigate the emergence probability of a fitter strain if it mutates from a self-limiting strain that is guaranteed to go extinct in the long-term. We show that ongoing immune cell proliferation during the initial stages of infection causes a drastic reduction in the probability of emergence of mutated strains; we further outline how this effect can be accurately measured. Further analysis of the model shows that, in the short-term, mutant strains that enlarge their replication rate due to evolving an increased growth rate are more favoured than strains that suffer a lower immune-mediated death rate (‘immune tolerance’), as the latter does not completely evade ongoing immune proliferation due to inter-parasitic competition. We end by discussing the model in relation to within-host evolution of human pathogens (including HIV, hepatitis C virus, and cancer), and how ongoing immune growth can affect their evolutionary dynamics.

Where is the growth rate, the death rate due to immunity, and y the density/population of immune cells. Immune growth is represnted with a logistic-growth model given the presence of infection: Here, r is the logistic growth rate of immunity, and K is the maximum/carrying capacity of the immune response.
In this model, if / < K (see below for a rational for this), then the first strain will increase in frequency until the immune-cells reach a critical size. After this point, the infection cells will decrease and go extinct, while the immune response will be maintained at a non-zero size. Below shows a numerical evaluation of this ODE set, with = 1, = r = 0.1, and K = 100. It is clear that with this example, the infection line goes extinct very quickly, after only 6-7 generations. This is because the maximum size that the infection line can achieve is equal to / unless K is very small, which we will show below.
In this model, if / < K (see below for a rational for this), then the first strain will increase in frequency until the immune-cells reach a critical size. After this point, the infection cells will decrease and go extinct, while the immune response will be maintained at a non-zero size. Below shows a numerical evaluation of this ODE set, with = 1, = r = 0.1, and K = 100. It is clear that with this example, the infection line goes extinct very quickly, after only 6-7 generations. This is because the maximum size that the infection line can achieve is equal to / unless K is very small, which we will show below.

sol = NDSolvex '[t] ⩵ (1 -− 0.1 y[t]) x[t], y '[t] ⩵ 0.1 x[t] y[t] 1 -− y[t]
100 , With / > K, then the parasite can escape immunity and rapildly increase in size. This is exemplified below with = 1, = r = 0.009 (so / ≃ 111). Solving differential equations, for use in potential analytical solution To make sure we know the general behaviour of the system, as well as work out when realistic behaviour arises, we will now analyse the basic properties of the system to determine the maximum growth rates, how to create analytical solutions for the infected cell-line, and so on. We start with the initial system of equations.
To make sure we know the general behaviour of the system, as well as work out when realistic behaviour arises, we will now analyse the basic properties of the system to determine the maximum growth rates, how to create analytical solutions for the infected cell-line, and so on. We start with the initial system of equations.
On the emergence probability, for a static and changing population From the form of x' [y], it is clear that the first infection will increase in frequency if R-y > 0, where R is the reproductive ratio in the absence of immunity, . For this model to make sense, R needs to be high initially, since the pathogen dies out if the number of immune cells y > R. For example, if R = 10, then the initial extinction probability equals 1 R = 0.1, quite low. But once y exceeds 10, extinction is certain. Therefore the same logic must apply with regards to emergence of the new pathogen. Say it has growth rate Φ and death rate (in absence of immunity) of Σ, then its baseline reproductive rate R2 = Φ Σ must exceed K in order for it to have any chance of completely escaping immune proliferation.
The total emergence probability is given by the 'evolutionary rescue' equation, 1-Exp[-∫x[y]Π[y]dy], where x is the population size of first cell-line and Π is the emergence probability. Preliminary results showed that using Π = 1-y/R2 greatly overestimated emergence probability compared to simulations, so we needed to account for the proliferation of the first parasite. We will first show how it is not possible to derive exact analytical solutions for Π, then how one can create an approximation based on previous solutions.
This term was then multiplied by the diffusion solution for emergence, 1 -− Exp-− 2 (R2-−y) R2+y , to create a scaled equation for the emergence probability. We can numerically find the y value causing this discontinuity, which represents the maximum possible y before emergence becomes inpermissable. Caution: If re-using this code, it is imperative to test whether the above function has a 'good' initial guess for finding the root of the denominator. Otherwise a different zero might be found, leading to inaccurate computation. For example, if the first guess at a root was set at 200 instead: Which is not correct.

Plots of Equations 5 and 10 (in the main text) for different parameter values
Equation 5 (dynamics of first strain as a function of immune size) We first re-define the function x1(y) and the maximum value of y, for ease of plotting: Below are plots for K = 1000, R1 = 100, = 0.5 (blue) and 9 (red). Below are plots for K = 10,000, R1 = 1000, = 0.5 (blue) and 9 (red). We also slightly redefine the function to find the maximum permissible value of y to change the starting value for finding the root, otherwise inaccurate values would be obtained. Plots for K = 100, R1 = 60, = 0.5, and R2 equals either 110, 130, or 150. Redefining the function to find the maximum permissible value of y to change the starting value for finding the root. Redefining the function to find the maximum permissible value of y to change the starting value for finding the root.

Pfix
Text S1.nb K = 1,000 simulations dataF = Import["results_6thFeb.dat", " Table"]; dataF2 = Import["results_10thFeb.dat", " Comparing models with and without feedbacks 1) Comparing this model to one not assuming feedbacks affecting emergence probability First, we redefine the system of numerical equations to be solved, to create a solution using a 'naive' estimate of emergence probability. The ratio function below compares the current model estimate to the naive solution that does not consider epidemilogical feedbacks. Irrespective of the parameters used, one sees that the naive estimate (no feedbacks) greatly underestimates actual emergence probability.