Coevolutionary Immune System Dynamics Driving Pathogen Speciation

We introduce and analyze a within-host dynamical model of the coevolution between rapidly mutating pathogens and the adaptive immune response. Pathogen mutation and a homeostatic constraint on lymphocytes both play a role in allowing the development of chronic infection, rather than quick pathogen clearance. The dynamics of these chronic infections display emergent structure, including branching patterns corresponding to asexual pathogen speciation, which is fundamentally driven by the coevolutionary interaction. Over time, continued branching creates an increasingly fragile immune system, and leads to the eventual catastrophic loss of immune control.

in the main text), drawing from previously measured and modeled characteristics of T-cells and well-studied viruses.
The parameter values shown in Table S1 are adopted or slightly modified from a previously published shape space model of T-cell dynamics [1] or a related model of thymic T-cell influx during chronic infection [2], in order to reproduce the basic acute infection immune system dynamics seen in these models. The rate of pathogen growth λ is consistent with [1], and similar to previously measured values of the initial expansion rate of HIV-1 [3,4]. The rate of pathogen killing by the immune response has been estimated to be between 0.1 and 0.5 day −1 near the peak of viremia in HIV and SIV infections [5,6]. The concentration of virus-specific T-cells responding to the primary infection near this peak can be roughly approximated as 10 4 µl −1 , based on measurements of the CD8 T-cell response to an infection of lymphocytic choriomeningitis virus (LCMV) in the mouse spleen [7,8]). Given these estimates, a value on the order of 10 −5 µl · day −1 , as used in [1] and [2], is reasonable for a pathogen killing rate per day, per concentration unit of T-cells in the adaptive response.
The T-cell replication and decay rates used here (denoted σ and δ, respectively) are consistent with those estimated from measurements of the immune response to LCMV in mice (replication rates between 1 and 3 day −1 , decay rates between 0.3 and 0.5 day −1 ) [8,9]. An approximate stimulation coecient κ is given by early-stage HIV data in [10], in which the significant T-cell response is observed to begin when viremia is on the order of 10 5 µl −1 .
From a previous shape space modeling paper [1], we adopt an initial naïve T-cell concentration L 0 of 3 µl −1 for each shape space site (T-cell clone). This is a biologically reasonable value with typical overall naïve T-cell concentrations. We can then estimate Γ, the naïve cell inux at a single shape space site (T-cell clone), from the initial equilibrium condition in the absence of pathogen, given in Eq. 10 in the main text: setting L0 = Γ/δ gives an approximate value of Γ = 1 (µl·site·day) −1 . We chose a value of φ = 10 10 µl −1 for the pathogen carrying capacity in order to give acute infection dynamics comparable to those in [10] and [1]. Table S1: Parameters adopted from previously published models [1,2], and consistent with several measurement and estimation studies on immune system and virus properties (see text).

Parameter
Sym. Value pathogen growth λ 3 day −1 pathogen killing We express mutation in terms of the overall mutation rate µ of the pathogens, in the commonly used units of substitutions per DNA base per replication cycle. Each epitope has on the order of 10 amino acids, for an estimate of about 30 base pairs; thus, under the assumption that all base pairs in the pathogen genome mutate independently at identical rates, we approximate the expected number of mutations per cycle for a given epitope x as Consequently, if every member of a population of P identical epitopes replicates once, we will have P × 30µ mutants. In terms of our model, then, µ ( x) can be interpreted as the fraction of pathogens at location x that differentiate to locations x = x, determined by the normalized mutation kernel Q and the simulation space size n: Here, the sum is performed over all n−1 sites x that are not equal to x. For our chosen kernel Q, as given in Eq. 7 in the main text of the paper, and n = 400, µ is nearly constant across the simulation space (i.e., < µ ( x) >≈ µ ( x), ∀ x), so we set its average value equal to our approximation in Eq. 1: Thus, having chosen Q and n, the value of the parameter χ is uniquely determined from µ. In this study, we use the standard value µ = 2.2 × 10 −5 substitutions per base pair per cycle, similar to commonly measured values for the rapidly mutating virus HIV-1 [11], and study the effects of variation between 1 × 10 −5 and 5 × 10 −5 . The binding specificity parameter b gives the width of the Gaussian affinity curve in shape space, setting the specificity of the naïve T-cell response. We choose b to correspond to typical values of the precursor frequency, or the probability of two different naïve T-cells responding to the same epitope (about 1/10 5 [12,13]). Starting from the uniformly distributed initial condition in which L( y) = Γ/δ ∀ y (Eq. 10 in the main text of the paper), the number of naïve T-cells under the affinity curve centered at some y 0 , and thus cross-reactive with the response to a pathogen maximally complementary to y 0 , is given in one dimension by Although we only simulate a small portion of the full shape space for computational efficiency, the total number of T-cells in our system is approximated by the initial naïve density times the total number of clones in the human immune system, with reasonable values between 10 6 − 10 7 [12].
To estimate b, we can then set the fraction of cross-reactive naïve T-cells equal to the precursor In this paper, we investigate values of b corresponding to the above range of total naïve T-cell clones in the system. We investigated the system under several values of ξ, which controls the strength of the damping homeostatic pressure. We chose values small enough that ξ << σR/|L tot − R| and ξ << δR/|L tot − R|, so that the homeostatic pressure would not significantly alter the expected chronic infection steady-state levels. In the main text of the paper we present results obtained at two different values of ξ. The higher value (Fig. 5 in main text) shows much stronger damping of total population density oscillations in acute infection, and converges to approximate chronic equilibrium quickly. Infections at the lower value of ξ (Fig. 4 in main text) take longer to converge to a chronic infection, branch, and eventually escape control. The branching and escape behaviors are observed consistently during chronic infections across values of ξ, and appear to be driven by the same mechanisms, although the timescale on which they occur changes, as well as the precise values of other parameters needed to generate a chronic infection.
In order to ensure that our choice to simulate only a limited region of the total shape space of naïve T-cell clones does not substantially affect our results, we investigated the behavior of several infections in a shape space with a substantially larger T-cell population (between 10 6 − 10 7 naïve clones). We found that extending the shape space in this way is effectively the same as lowering the parameter ξ: as discussed in the main text, the timescales of convergence to chronic infection, branching, and eventual immune escape are lengthened in a large shape space, while the overall effectiveness of the immune system is enhanced. However, the mechanisms leading to chronic infection, branching, and eventual pathogen escape are similar, though occurring on different timescales. In the main text and figures, we thus focus on presenting results from a small portion of shape space for computational efficiency and more effective visualization.