Sequential mutations in exponentially growing populations

Abstract

Stochastic models of sequential mutation acquisition are widely used to quantify cancer and bacterial evolution. Across manifold scenarios, recurrent research questions are: how many cells are there with n alterations, and how long will it take for these cells to appear. For exponentially growing populations, these questions have been tackled only in special cases so far. Here, within a multitype branching process framework, we consider a general mutational path where mutations may be advantageous, neutral or deleterious. In the biologically relevant limiting regimes of large times and small mutation rates, we derive probability distributions for the number, and arrival time, of cells with n mutations. Surprisingly, the two quantities respectively follow Mittag-Leffler and logistic distributions regardless of n or the mutations' selective effects. Our results provide a rapid method to assess how altering the fundamental division, death, and mutation rates impacts the arrival time, and number, of mutant cells. We highlight consequences for mutation rate inference in fluctuation assays.

Fig 1. Model schematic.

A: We consider a multitype branching process in which cells can divide, die, or mutate to a new type. B: We study the waiting time until a cell of the nth type exists, τ_n, starting with a single cell of type 1. C: Stochastic simulation of the number of cells over time, with dashed lines indicating the large-time trajectories given by Eq (1). Grey horizontal line occurs at the inverse of the mutation rate, while the grey vertical lines indicate the time at which the type n population size reaches the inverse of the mutation rate, which gives the arrival time of the type n + 1 cells to leading order. Parameters: α₁ = α₃ = 1.1, α₂ = 1, β₁ = 0.8, β₂ = 0.9, β₃ = 0.5, ν₁ = ν₂ = 0.01. Thus, the net growth rates are λ₁ = 0.3, λ₂ = 0.1, λ₃ = 0.6 and the running-max fitness follows δ₁ = δ₂ = λ₁, δ₃ = λ₃.

Fig 2. Comparison with prior work and motivating examples.

A. Previous work has considered special cases of growth rate sequences, here we consider general sequences as long as λ₁ > 0. B. Two biological scenarios in which the growth rate sequences covered in this paper are relevant: the acquisition of driver mutations in the canonical carcinogenesis pathway of colorectal cancer, and the accumulation of neoantigens by cancer cells which results in increased cell death due to immune system surveillance.

Fig 3. Comparison of limiting Mittag-Leffler distribution for the number of type n cells with stochastic simulations.

Eq (1), states that for large times and small mutation rates, the scaled number of type n cells, $e^{-{\delta }_{n}t}{t}^{-\left(r_n-1\right)}{Z}_n(t)\approx V_n$, is approximately Mittag-Leffler distributed with scale ω_n and tail λ₁/δ_n. Here, we compare simulations of the scaled number of type n divided by ω_n, to the density of V_n/ω_n which is Mittag-Leffler with scale parameter 1, and tail parameter λ₁/δ_n ∈ (0, 1]. We chose three tail parameter values λ₁/δ_n = 0.25, 0.5, 1.0, and these curves are depicted with solid lines. The simulation parameter were always α₁ = 1.2, β₁ = 0.2, ν₁ = 0.01, β₂ = 0.3 and for n = 2 types sim 1: α₂ = 4.3, t = 5; sim 2: α₂ = 2.3, t = 7; sim 3: α₂ = 1.0, t = 12. Then for n = 3 types sim 4: as in sim 3 plus α₃ = 2.4, β₃ = 0.4, ν₃ = 0.001, t = 12. Density lines were created in Mathematica using x^γ−1MittagLefflerE[γ, γ, −x^γ].

Fig 4. Comparison of limiting logistic distribution for arrival times with stochastic simulations.

Normalized histogram for the arrival times of types 1–3 obtained from 1000 simulations of the exact model versus the probability density corresponding to the logistic distribution of Eq (4). Note the shape of the distribution remains unchanged. Parameters: α₁ = α₃ = 1, α₂ = 1.4, ν₁ = ν₂ = ν₃ = 0.01, β₁ = β₂ = 0.3, β₃ = 1.5.

Fig 5. Statistical inference for an n-mutation fluctuation assay.

A. Schematic of a fluctuation assay for the measurement of mutation rates when n mutations are required for resistance. Drug sensitive cells are initially cultured, and after growth for a given time t, the cells are exposed to a selective medium. Non-resistant cells are killed, revealing the number of mutants. This experiment is conducted over replicates, and the number of replicates without resistance and the mutant numbers are recorded. B. Likelihood inference on a simulated fluctuation assay assuming: 2 mutations are required for resistance, 100 replicates, no death, α_i = 1 for each i, t = 10, and the mutation rate ν stated on the x-axis. Wide error bars are expected when using the p₀ method for $t\gg t_{1/2}^{\left(n\right)}$ as only a small number of replicates have no resistant cells; in such a setting using the mutant counts (right panel) provides superior inference. Likewise, if $t\approx t_{1/2}^{\left(n\right)}$ the approximation of Eq (1) is not appropriate, which explains the inaccurate inference for log₁₀(ν) = −3 when using the mutant counts; the p₀ method provides improved inference in this scenario.