Life cycle synchronization is a viral drug resistance mechanism

Viral infections are one of the major causes of death worldwide, with HIV infection alone resulting in over 1.2 million casualties per year. Antiviral drugs are now being administered for a variety of viral infections, including HIV, hepatitis B and C, and influenza. These therapies target a specific phase of the virus’s life cycle, yet their ultimate success depends on a variety of factors, such as adherence to a prescribed regimen and the emergence of viral drug resistance. The epidemiology and evolution of drug resistance have been extensively characterized, and it is generally assumed that drug resistance arises from mutations that alter the virus’s susceptibility to the direct action of the drug. In this paper, we consider the possibility that a virus population can evolve towards synchronizing its life cycle with the pattern of drug therapy. The periodicity of the drug treatment could then allow for a virus strain whose life cycle length is a multiple of the dosing interval to replicate only when the concentration of the drug is lowest. This process, referred to as “drug tolerance by synchronization”, could allow the virus population to maximize its overall fitness without having to alter drug binding or complete its life cycle in the drug’s presence. We use mathematical models and stochastic simulations to show that life cycle synchronization can indeed be a mechanism of viral drug tolerance. We show that this effect is more likely to occur when the variability in both viral life cycle and drug dose timing are low. More generally, we find that in the presence of periodic drug levels, time-averaged calculations of viral fitness do not accurately predict drug levels needed to eradicate infection, even if there is no synchronization. We derive an analytical expression for viral fitness that is sufficient to explain the drug-pattern-dependent survival of strains with any life cycle length. We discuss the implications of these findings for clinically relevant antiviral strategies.

1 Analytic results for single-strain model with constant drug level 1

.1 Maturation happens in n intermediate steps
Consider the single-strain viral dynamics model with n maturation steps (Eq. (1)) where m i = m n for all i. For a constant infection efficacy β, this system of equations has a unique steady state where infection does not occur, for which x * = λ/d x , y = w i = 0, ∀i ≤ n. This steady state is stable if R 0 < 1, where R 0 is the basic reproductive ratio for the multi-stage model as defined in Eq.(4).
The system also has a unique stable steady state where infection occurs (see S12 Fig) when R 0 > 1, such that:

Maturation happens after a fixed time delay
Consider the single-strain viral dynamics model with fixed maturation time τ = 1/m (Eq. (3)). For a constant infection efficacy β, this system of equations has a unique steady state where infection does not occur, for which x * = λ/d x , y = 0. This steady state is stable if R 0 < 1, where R 0 is the basic reproductive ratio for the fixed-time model as defined in Eq. (5).
The system also has a unique stable steady state where infection occurs when R 0 > 1, such that: 2 Analytic results for multi-strain model with constant drug levels 2

.1 Maturation happens in n intermediate steps
We used the following system to model competition between N viral strains with different maturation times that co-infect the same host, when the maturation process happens in n steps: Here i is the index of the maturation phase (1 ≤ i ≤ n) and k is the index of the strain (1 ≤ k ≤ N ). For a constant infectivity β, this model has a unique steady state where infection does not occur, for which x * = λ/d x , y * k = w * ki = 0, ∀ i ≤ n, ∀ k ≤ N . This steady state is stable if R k < 1, ∀ k = 1 . . . N .
When d w = 0, the system has an infinite number of stable steady-states where infection occurs (see S13 Fig), such that: Similarly, when d w > 0, the system has N steady states where infection occurs, for which: is the basic reproductive ratio of the strain with maturation rate m k .
R k is a monotonically increasing function of m k for d w = 0, and does not depend on m k when the death rate d w of immature infected cells is zero (d w = 0). The same thing can be said about the steady-state value of mature infected cells y k , while the steady-state value of immature infected cells is inversely related to the maturation rate m k for all d w .
The equilibrium described by Eq. (S.5) implies competitive exclusion: only one strain survives at equilibrium and all others go extinct. Only 1 out of N of the possible equilibria are stable, and this is the one in which the strain with the highest maturation rate and thus the highest R k survives. All other infected equilibria are unstable (see S14 Fig).

Maturation happens after a fixed time delay
We used the following system to model competition between N viral strains with different fixed maturation times τ k = 1/m k , k = 1 . . . N , that co-infect the same host: For a constant infectivity β, the multi-strain competition model with fixed maturation times has a unique steady state where infection does not occur, for which x * = λ/d x , y * k = 0, ∀k ≤ N . This steady state is stable if R k < 1, ∀k = 1 . . . N , where R k is the basic reproductive ratio for the model with fixed maturation time as defined in Eq. (5).
When d w = 0, the system has an infinite number of stable steady-states where infection occurs such that: When d w > 0, the system has N steady states where infection occurs, for which: , for a unique strain k y * j = 0, ∀τ j = τ k , j = 1 . . . N (S. 9) Similarly to Eq. (S.5), the only stable equilibrium resulting from Eq. (S.9) is the one where only the strain with the shortest maturation time (highest R k ) survives. Therefore at steady state, y * j = 0, where τ k < τ j , ∀j = k, j = 1 . . . N . All other infected equilibria are unstable.

Initial conditions
The initial conditions for our analysis correspond to the levels of uninfected cells and infected cells of different life cycle stages present before drug treatment is administered. For the single-strain models, we use as initial conditions the infected equilibrium of the model in the absence of drug (i.e. when the drug efficacy f = 0). Explicitly, the initial condition for the single-strain model with n maturation steps is calculated using Eq. (S.1). For the single-strain model with fixed maturation time, the initial condition (x(t ≤ 0) and y(t ≤ 0)) is calculated using Eq. (S.2). For the multi-strain competition models (with N strains), the situation is more complicated, because either there exists an infinite number of possible equilibria containing all strains (d w = 0) or no equilibrium including all strains (d w > 0). Therefore, we had to either choose or create an initial condition which contained a balanced representation of each strain.
For the multi-strain competition model with n maturation steps, we set the following initial conditions: For the multi-strain competition model with fixed maturation times, we set the following initial conditions (for all pre-treatment time points): These initial conditions were set for both the deterministic and stochastic multistrain models; for the stochastic implementation, the cell concentrations were rounded to the nearest integer.
When d w = 0, these expressions represent one of the (infinite) possible no-drug equilibria (Eqs. (S.4) and(S.8)). When d w > 0, the only stable no-drug equilibria are ones where only a single strain persists (Eqs. (S.5) and(S.9)), but here, we construct a distribution of all N strains such that each one exists at 1/N of the level it would exist at if it were the only strain in the population.

Equilibrium conditions
Both the single-strain and multi-strain deterministic simulations were terminated when all strains had a relative change of less than 5% over consecutive 10 day periods in each cell type (healthy, infected-mature, infected-immature).
The stochastic competition simulation was terminated when a single strain survived in the population, and the relative change in each cell type (healthy, infected-mature, infected-immature) was less that 0.5% over consecutive 5 day periods. These conditions also allow for no strains to fix in the population.
Plots of the cell concentrations as a function of time were used to verify that the above thresholds were necessary and sufficient in identifying the equilibrium state of the simulated systems.

Alternative model with costs to viral life cycle length
Throughout this paper we have assumed that time spent by infected cells in immature phases is either cost-less or costly, depending on the death rate of immature infected cells (d w ). However, in reality there may be a benefit to longer time spent in the immature phases, because it is during these phases that the precursor compounds used to construct new viral particles are produced, and hence maturation time may be positively correlated with the eventual viral burst size and infectivity. Here we present a simple model to take these effects into account.
Consider the viral dynamics model with fixed maturation time (Eqs. (3)). Assume that immature infected cells produce a precursor molecule (concentration defined by state variable P ), that is needed for infectious virions to be produced when the cell advances to a mature phase. Let's assume that P is produced in a single immature infected cell at rate α, and that it decays at a per capita rate δ. Therefore, in a given immature infected cell, P followṡ Since an infected cell only stays in the immature phase for time τ , this equation only holds for t ∈ (0, τ ). We assume that there is initially no precursor molecule, so P (0) = 0. Therefore, the concentration of P over time follows We want to be able to vary the production rate of P without varying the eventual steady-state level as t → ∞, and we are not worried about the units of P , so we set α = δ and have the final concentration of P when the cell becomes mature as Now we assume that the amount of this precursor molecule directly translates into the amount of virions that can be produced by mature infected cells, and therefore into the viral burst size k and net infectivity β. We then get that viral infectivity depends on the maturation time τ , even in the absence of fluctuating drug levels: When drug levels are not constant, β would also depend on the time-dependent drug efficacy, (t), which is given by D(t) in the on-off drug model (Eq.(6)), and 1/(1 + (D(t)/IC 50 ) M ) in the pharmacological model (Eq. (7)).
S10 Fig A shows the dependence of viral infectivity on maturation time under this model with constant drug levels. When maturation time τ is less than 1/δ, infectivity increases approximately linearly with maturation time, and so there is a clear cost to maturing faster. However when maturation time is much longer than 1/δ, infectivity is roughly constant with maturation time, so there is little cost to changing maturation time. Therefore this model can be used to go between situations where there is and where there is not a cost to speeding up maturation.
We can also examine how this cost to maturation time influences net fitness by considering the basic reproductive ratio (R 0 ). Including β(τ ) in the expression for R 0 for the fixed maturation time model, we get While allowing death of immature infected cells (d w > 0) increases fitness of strains who mature faster, including production of the precursor decreases fitness of quickly maturing strains. In the presence of constant drug levels there is trade-off between these two effects. The optimal maturation time is We additionally want to ensure that R 0 (τ c ) is independent of δ, so that δ only controls what the optimal maturation time would be in the absence of drug, but not what the viral fitness would be at this optimum. To do this, we need to choose a different value for the production rate α, and we instead use the following τ -dependent infectivity β With this scaling, the basic reproductive ratio in the absence of drug becomes:  Table 1, we consider three different regimes based on the rate of precursor production and decay, δ

The maturation process happens with no intermediate maturation steps
It is instructive to first consider the case where β(t) = β(t + T ) is time-dependent and periodic, but there are no intermediate maturation steps, so that newly infected cells begin producing new virions immediately. This simpler case helps to foster intuition for the more intricate calculations that follow. The viral dynamics are specified by the following system: We assume that the system starts from the uninfected equilibrium, and a small amount of the virus is introduced. If we focus on the early-time dynamics, when the density of infected cells is low, then we haveẏ This equation can be rewritten as Our goal is to use Equation (S.14) to determine how y(L) is related to y(L+T ). (Here, L is any early time in the dynamics such that the density of target cells is approximately at its uninfected value, so that Equation (S.14) holds.) There is a simple intuition for why we are interested in how y(L) is related to y(L + T ): Depending on the values of the model parameters, when tracked over a longenough time period, y(t) will tend to either grow or decay, since the particular viral strain under consideration will either establish a persistent infection or go extinct. Moreover, β(L) = β(L+T ) is periodic, and this periodicity of β(t) influences the infection level, y(t). Therefore, we expect that y(t) responds to β(t) in an irregularly periodic pattern. Indeed, y(t) shows quasiperiodic behavior in response to β(t). Specifically, y(t) has the following form [1,2]: So, if we determine how y(L) is related to y(L + T ), then the periodic prefactor, Y (L) = Y (L + T ), cancels out of the analysis (for example, when computing y(L + T )/y(L)), and we are able to solve forr 0 , which is key: Ifr 0 > 0, then, when sampled at times that are integer multiples of the drugdosing period, the density of infected cells increases in time, and infection is established. Ifr 0 < 0, then, when sampled at times that are integer multiples of the drug-dosing period, the density of infected cells decreases in time, and infection is eliminated.
To find a connection between y(L) and y(L + T ), we begin by integrating Equation (S.14) between time L and time L + t k ; we obtain (In Equation (S.16), similarly to L, t k is any early time in the dynamics such that the density of infected cells is low, so that Equation (S.14) holds. t k−1 is just an integration variable, but the notation t k−1 for the integration variable is convenient, as it implies the manipulations that are to follow. The intuition behind the notation for the integration variables will become clearer in the following two sections. Also, for simplicity of notation, we assume that L is equal to an integer multiple of T , so that β(L + t k−1 ) = β(t k−1 ).) We proceed by noting that Equation (S.16) can be repeatedly substituted into itself. To see how this works, we set k = 0 and t 0 = T in Equation (S.16); we obtain Next, we use Equation (S.16) to substitute for e dyt −1 y(L + t −1 ) in Equation (S.17); we obtain Continuing, we use Equation (S.16) to substitute for e dyt −2 y(L + t −2 ) in Equation (S.18); we obtain Continuing further, we use Equation (S.16) to substitute for e dyt −3 y(L + t −3 ) in Equation (S.19), then we use Equation (S.16) to substitute for e dyt −4 y(L + t −4 ) in the resulting equation, etc.; continuing this process ad infinitum, we arrive at where M is given by (Throughout this Supplement, we define the empty product, as for j = 1 in the product in Equation (S.21), as being equal to 1.) In Equation (S.21), t −1 , t −2 , t −3 , etc. are all integration variables. The lowest-order term on the right-hand side features no integration. The next-lowest-order term on the right-hand side features a single integration over t −1 . The next-next-lowest-order term on the right-hand side features a double integration over t −2 and t −1 . This continues, so that higher-order terms on the right-hand side have higher-dimensional integrals that must be computed. These integrals over β(t) in Equation (S.21) can be done analytically, if possible, or can be performed using, e.g., Monte Carlo integration [7,8]. Monte Carlo integration is a numerical technique in which random or pseudorandom numbers are used to approximately determine the value of a definite integral. To see how it works, consider β(t) given by Equation (6) with β 0 = 1. For this choice of β(t), if we want to numerically evaluate T 0 dt −1 β(t −1 ), then we can pick a random number, u −1 , between 0 and T . If β(u −1 ) = 0, then u −1 is recorded as "failure"; otherwise, u −1 is recorded as "success". This Bernoulli trial is repeated for many random values of u −1 . Our estimate for the value of the definite integral is then given simply by T multiplied by the fraction of trials that were successes. By increasing the number of such Bernoulli trials, we can achieve arbitrarily high precision in our numerical estimate of the definite integral. Also for this choice of β(t), if we want to numerically evaluate , then we can pick two random numbers, u −1 and u −2 , each between 0 and T . If u −2 > u −1 , β(u −1 ) = 0, or β(u −2 ) = 0, then the trial pair of values u −1 and u −2 are recorded as "failure"; otherwise, the trial pair of values u −1 and u −2 are recorded as "success". This Bernoulli trial is repeated for many random pairs of values of u −1 and u −2 . Our estimate for the value of the definite integral is then given simply by T 2 multiplied by the fraction of trials that were successes. Monte Carlo integration is especially useful for evaluating high-dimensional integrals, which often do not feature nice analytical solutions; indeed, this approach is useful for numerically evaluating the high-dimensional integrals in the high-order terms resulting from the series expansion of Equation (S.21). (Although we assume the simple form of β(t) given by Equation (6) for these examples of how to implement Monte Carlo integration, the analytical results in this section and the following two sections hold for arbitrary functional forms of β(t), and for any case of β(t), a variation of the Monte Carlo integration technique described above can be applied to evaluate the definite integrals that appear in our formulae.) The series in Equation (S.21) converges to the quantity, M , that we are seeking: y(L), when multiplied by M , yields y(L + T ). Thus, if M > 1, then, when sampled at times that are integer multiples of the drug-dosing period, the density of infected cells increases in time, and infection is established. Or, if M < 1, then, when sampled at times that are integer multiples of the drug-dosing period, the density of infected cells decreases in time, and infection is eliminated.
We can proceed further by substituting Equations (S.15) into Equation (S.20), which allows us to solve for the growth rate,r 0 , of the viral strain when sampled at integer multiples of the drug-dosing period. Performing these substitutions, we obtain where M is given by For a nontrivial solution to Equation (S.22) (i.e., a solution for which y(L) is nonzero), we can solve explicitly forr 0 : It is also helpful to define a quantity,R 0 , which has the following properties: •R 0 is nonnegative.
Thus, the conditionR 0 > 1 indicates the establishment of infection, while the conditionR 0 < 1 indicates the elimination of infection.
• If β(t) = β is constant, thenR 0 simply reduces to the basic reproductive ratio, denoted here as R 0 .
Thus, the quantityR 0 for the case of a time-varying drug effectiveness, β(t), is analogous to the basic reproductive ratio that is widely used for analyzing the more-familiar case of a constant value of β [9]. We therefore define the parameterR 0 as The exercise above aids in understanding the derivations that follow for more complicated cases, but as a consistency check,r 0 in Equation (S.24) andR 0 in Equation (S.25) must also be given much more simply by r 0 and R 0 for the case of constant drug levels if we replace β in those simple formulae with the time-average of β(t) over a single drug-dosing period, denoted by β(t) .

The maturation process happens in n intermediate maturation steps
We now apply a variation of the procedure used for the case of no intermediate maturation steps to the case of a finite number, n, of intermediate maturation steps. Again, we have a time-dependent and periodic β(t) = β(t + T ). The viral dynamics are given bẏ We assume that the system starts from the uninfected equilibrium, and a small amount of the virus is introduced. If we focus on the early-time dynamics, when the density of infected cells is low, then we haveẏ (Here, we choose to write the equation forẏ(t) first simply for notational convenience, as will become apparent.) Similarly to the analysis in the previous section, we can use an integrating factor on each equation and then integrate each resulting equation between time L and time L + t k to obtain (Here, just as in the previous section, L and t k are any early times in the dynamics such that the density of infected cells is low, so that Equations (S.28) hold. The notation t k−1 for the integration variable is convenient, as it implies the manipulations that are to follow. Also, for simplicity of notation, we assume that L is equal to an integer multiple of T , so that β(L + t k−1 ) = β(t k−1 ).) Our approach is similar to that of the previous section: We would like to find a linear operator that tells us how y(L) and w i (L) for all 1 ≤ i ≤ n are related to y(L + T ) and w i (L + T ) for all 1 ≤ i ≤ n. It is helpful to start by exploring Equations (S.29) for some small values of n.
For the case n = 1, there is a single intermediate maturation step. In this case, Equations (S.29) reduce to the following two equations: If we set k = 0 and t 0 = T in Equations (S.30), then we obtain Next, we substitute the first of Equations (S.30) into the second of Equations (S.31), and we substitute the second of Equations (S.30) into the first of Equations (S.31); we obtain e dyT y(L + T ) = y(L) Next, we substitute the first of Equations (S.30) into the first of Equations (S.32), and we substitute the second of Equations (S.30) into the second of Equations (S.32); we obtain e dyT y(L + T ) = y(L) Next, we substitute the first of Equations (S.30) into the second of Equations (S.33), and we substitute the second of Equations (S.30) into the first of Equations (S.33). Notice that if we continue this procedure ad infinitum, then y(L + T ) and w 1 (L + T ) can be expressed as where the M ij are infinite series that are functions of the model parameters and functionals of the infectivity, β(t). reduce to the following three equations: If we set k = 0 and t 0 = T in Equations (S.35), then we obtain e dyT y(L + T ) = y(L) Next e (m 2 +d 2 )T w 2 (L + T ) = w 2 (L) Next, we substitute the first of Equations (S.35) into the first of Equations (S.38), we substitute the second of Equations (S.35) into the second of Equations (S.38), and we substitute the third of Equations (S.35) into the third of Equations (S.38). Notice that if we continue this procedure ad infinitum, then y(L + T ), w 1 (L + T ), and w 2 (L + T ) can be expressed as where the M ij are infinite series that are functions of the model parameters and functionals of the infectivity, β(t).
Case of any n ≥ 1 Following the procedure described above, for any n ≥ 1, we obtain Recall that our goal is to solve for the linear operator, M ij . Although the procedure outlined above shows us how to do this, it is evident that the mathematical expressions are symbolically quite cumbersome. Therefore, we simplify notation substantially by making several definitions.
First, the notation in Equations (S.40) is more natural if we define Then, Equations (S.40) simply feature w i (L) and w i (L + T ) for 0 ≤ i ≤ n, which is notationally much more convenient. Next, since the dynamical quantities, w i (t), are indexed by 0 ≤ i ≤ n, and since each M ij is an infinite series in which factors of w i (L) appear in terms in cyclical fashion, it is natural to define Λ(i, n) ≡ i mod (n + 1) Then, for any integer value of i, the function Λ(i, n) simply outputs an integer between 0 inclusive and n inclusive.
Again, noticing that each M ij is an infinite series in which factors appear in terms in cyclical fashion, we strive to simplify notation in the series expansion of each M ij to the maximum extent possible. We therefore make the following definitions, whose notational utility will become apparent: Now, with the preceding definitions, notice that Equations (S.29) can be written very compactly: In Equations (S.41), k ≤ 0, and q can be any integer. Finally, notice that in Equations (S.33) for n = 1 and in Equations (S.38) for n = 2, each integral over t −1 has an upper limit of t −1 = T . Therefore, it simplifies notation further if we define Let us now revisit our calculations for the simple cases of n = 1 and n = 2 to see how the notation with the above substitutions is dramatically simplified.

Particular case: n = 1
For the case n = 1, recall that, depending on the value of i, the function Λ(i, n) reduces to 0 or 1. Therefore, with the definitions above, Equations (S.41) reduce to the following two equations: If we set k = 0 and t 0 = T in Equations (S.42), then we obtain F Λ(0,1) (T ) = F Λ(0,1) (0) Next, we substitute the first of Equations (S.42) into the second of Equations (S.43), and we substitute the second of Equations (S.42) into the first of Equations (S.43); we obtain Next, we substitute the first of Equations (S.42) into the first of Equations (S.44), and we substitute the second of Equations (S.42) into the second of Equations (S.44); we obtain F Λ(0,1) (T ) = F Λ(0,1) (0) Next, we substitute the first of Equations (S.42) into the second of Equations (S.45), and we substitute the second of Equations (S.42) into the first of Equations (S.45). We continue this procedure ad infinitum. Using the definitions above, it can be verified by direct substitution that Equations (S.42), (S.43), (S.44), and (S.45) reduce to Equations (S.30), (S.31), (S.32), and (S.33), respectively. However, the notational and conceptual utility of the substitutions described above for simplifying the writing is evident.
We see that the M ij in Equations (S.34) are given by Particular case: n = 2 For the case n = 2, recall that, depending on the value of i, the function Λ(i, n) reduces to 0, 1, or 2. Therefore, with the definitions above, Equations (S.41) reduce to the following three equations: If we set k = 0 and t 0 = T in Equations (S.46), then we obtain Next, we substitute the first of Equations (S.46) into the second of Equations (S.47), we substitute the second of Equations (S.46) into the third of Equations (S.47), and we substitute the third of Equations (S.46) into the first of Equations (S.47); we obtain Next, we substitute the first of Equations (S.46) into the third of Equations (S.48), we substitute the second of Equations (S.46) into the first of Equations (S.48), and we substitute the third of Equations (S.46) into the second of Equations (S.48); we obtain F Λ(0,2) (T ) = F Λ(0,2) (0) Next, we substitute the first of Equations (S.46) into the first of Equations (S.49), we substitute the second of Equations (S.46) into the second of Equations (S.49), and we substitute the third of Equations (S.46) into the third of Equations (S.49). We continue this procedure ad infinitum. Using the definitions above, it can be verified by direct substitution that Equations (S.46), (S.47), (S.48), and (S.49) reduce to Equations (S.35), (S.36), (S.37), and (S.38), respectively. However, the notational and conceptual utility of the substitutions described above for simplifying the writing is again evident.
We see that the M ij in Equations (S.39) are given by Case of any n ≥ 1 If we set k = 0 and t 0 = T in Equations (S.41), then we obtain n + 1 equations: By repeatedly substituting for F Λ(q−1,n) (t −1 ), F Λ(q−2,n) (t −2 ), etc. in Equations (S.50), the n + 1 equations become Equivalently, From Equations (S.51), we see that the M ij in Equations (S.40) are given by The integrals over β(t) in Equation (S.52) can be done analytically, if possible, or can be performed using, e.g., Monte Carlo integration. (Please see the section on the derivation ofr 0 for the case of no intermediate maturation steps for a description of how Monte Carlo integration works.) The series in Equation (S.52) converges to the linear operator, M ij , that we are seeking. (In the language of Floquet theory, M ij is the monodromy matrix.) Thus, if the largest eigenvalue of M ij is greater than 1, then, when sampled at times that are integer multiples of the drug-dosing period, the density of infected cells increases in time, and infection is established. Or, if the largest eigenvalue of M ij is less than 1, then, when sampled at times that are integer multiples of the drug-dosing period, the density of infected cells decreases in time, and infection is eliminated.
Similarly to the case of no intermediate maturation steps, each w i (t) shows quasiperiodic behavior in response to β(t). Specifically, w i (t) for 0 ≤ i ≤ n have the following form [1,2]: We can proceed further by substituting Equations (S.53) into Equation (S.40), which allows us to solve for the growth rate,r 0 , of the viral strain when sampled at integer multiples of the drug-dosing period. Performing these substitutions, we obtain The maximum value ofr 0 for which Equation (S.55) is satisfied specifies the growth rate of the virus when each w i (t) is sampled at times that are integer multiples of the drug-dosing period, T . Ifr 0 > 0, then, when sampled at times that are integer multiples of the drug-dosing period, the density of infected cells increases in time, and infection is established. Ifr 0 < 0, then, when sampled at times that are integer multiples of the drug-dosing period, the density of infected cells decreases in time, and infection is eliminated.
Using what we have already derived, we can write M ij for the simple cases of n = 1 and n = 2.
Particular case: n = 1 If n = 1, then the M ij are given by Particular case: n = 2 If n = 2, then the M ij are given by Case of any n ≥ 1 For any n ≥ 1, the M ij are given by

The maturation process happens after a fixed time delay
If the number of intermediate maturation steps, n, is large, then the distribution of maturation times for newly infected cells becomes sharply peaked about the mean maturation time. Therefore, we suppose that an infected cell only begins producing new virions after a fixed maturation time, τ , has elapsed since its infection. We thus consider the viral dynamics specified by the following system of equations:ẋ We assume that the system starts from the uninfected equilibrium, and a small amount of the virus is introduced. If we focus on the early-time dynamics, when the density of infected cells is low, then the equation forẏ(t) becomeṡ Similarly to the analysis in the previous section, we can use an integrating factor on Equation (S.58) and then integrate this equation between time t and time t to obtain e dyt y(t ) = e dyt y(t) ds β(s)e dys y(s) (S.59) (Here, t and t are any early times in the dynamics such that the density of infected cells is low, so that Equation (S.58) holds.) Before proceeding, it is helpful to explain our intuition for how to determine the evolutionary success of a viral strain with a fixed maturation time.
In the previous section, for n intermediate maturation steps, our aim was to solve for a linear operator, M ij , that, when acting on the vector w i (L), produces the vector w i (L + T ) (where recall that w 0 (t) is simply a convenient notation for y(t)). The largest eigenvalue of M ij for a particular patient, viral strain, and drug-dosing regimen then determines if that viral strain establishes persistent infection or goes extinct. Such an approach is not directly applicable here because, for a fixed delay, τ , we are only considering the evolution of y(t). However, notice that we also cannot use the simple procedure for the case with no intermediate maturation steps because Equation (S.59) contains an integration in which the upper and lower integration limits are both offset by an amount τ . We thus need a slightly different approach.
Physical intuition is helpful. Suppose that the system begins evolving according to Equation (S.58) at time t = 0. The initial data for this problem is given by specifying β(t) and y(t) between times t = −τ and t = 0. At times shortly after the evolution begins, the initial data will certainly affect the evolution of the density of infected cells. As the evolution progresses, the initial data will become less important, and eventually, the density of infected cells will settle into a quasiperiodic temporal pattern that is determined by the form of the periodic infectivity, β(t). But at times that are too long after the evolution begins, the density of target cells can no longer be assumed to be at its uninfected level.
For performing a meaningful and tractable analytical analysis of this problem, we therefore want to focus on times that are sufficiently long after the evolution begins that the influence of the initial data on the temporal pattern of y(t) has decayed away, yet we must still assume that times are sufficiently short that the density of target cells remains approximately at its uninfected level. Within this time regime, y(t) has a similar form as was described in the previous two sections: To see how to simplify the analysis of Equations (S.59) and (S.60), it is helpful to begin by example.

Case of any τ
Notice from above that regardless of the particular value of τ , we obtain a solution forr 0 by solving det M ij = 0 Also notice that the dimensionality of the matrix M ij depends on the value of the maturation time, τ , in relation to the drug-dosing period, T . We therefore denote the dimensionality of the matrix M ij as U (τ, T ). Let us again consider the simple cases described above. If τ = T , then the dimensionality of M ij is equal to 1. If 2τ = T , then the dimensionality of M ij is equal to 2. If 3τ = T , then the dimensionality of M ij is equal to 3. In each of these cases, we have that U (τ, T ) is given by the lowest positive integer such that τ U (τ, T ) is equal to an integer multiple of T . This holds generally. For example, if 2τ = 3T , then the dimensionality of M ij is equal to 2. As another example, if 3τ = 2T , then the dimensionality of M ij is equal to 3.
Following the procedure described above, for any value of τ , we obtain Although the procedure outlined above shows us how to solve for M ij , it is evident that the mathematical expressions for M ij are symbolically quite cumbersome. Therefore, we simplify notation substantially by making several definitions.
First, rather than writing Y ((L + 1)τ ), Y ((L + 2)τ ), Y ((L + 3)τ ), etc., we can more simply write Y L+0 , Y L+1 , Y L+2 , etc. Therefore, we define Y L+q ≡ Y ((L + q + 1)τ ) Next, since only the Y values Y L+0 , Y L+1 , . . ., Y U (τ,T )−2 , Y U (τ,T )−1 appear in Equations (S.85), it is notationally useful to have a function, Λ (q, τ, T ), that takes an integer, q, along with τ and T as inputs and that outputs an integer between 0 and U (τ, T ) − 1. We therefore define Λ (q, τ, T ) ≡ (q − 1) mod U (τ, T ) We would also like to simplify the notation for the integrations appearing in the infinite series in the solution for M ij . To this end, consider Equation (S.59) with the following substitutions: t = (L + q)τ + t k , t = (L + q)τ , and s = (L + q − 1)τ + t k−1 : e dyt k y((L + q)τ + t k ) = y((L + q)τ ) Here, Lτ and t k are any times within the time regime of interest. Since we will be using the function Λ (q, τ, T ) to simplify the writing, q can be any integer in Equation (S.86). The notation t k−1 for the integration variable is convenient, as it implies the manipulations that are to follow.
Since each M ij is an infinite series in which factors appear in terms in cyclical fashion, we strive to simplify notation in the series expansion of each M ij to the maximum extent possible. We therefore make the following definitions, whose notational utility will become apparent: Now, with the preceding definitions, notice that Equations (S.86) can be written very compactly: In Equations (S.87), k ≤ 0, and q can be any integer. Also, notice that in Equation (S.65) for U (τ, T ) = 1, in Equations (S.73) for U (τ, T ) = 2, and in Equations (S.82) for U (τ, T ) = 3, each integral over t −1 has an upper limit of t −1 = τ . Therefore, it simplifies notation further if we define ξ i (q) ≡ ξ i (q, τ ) Setting k = 0 and t 0 = τ in Equation (S.87), we have Let us now revisit our calculations for the simple cases of τ = T , 2τ = T , and 3τ = T . More generally, since the respective mathematical steps are identical, let us consider the calculations for the simple cases of U (τ, T ) = 1, U (τ, T ) = 2, and U (τ, T ) = 3 to see how the notation with the above substitutions is dramatically simplified. For the case U (τ, T ) = 1, recall that the function Λ (q, τ, T ) reduces to 0. Therefore, with the definitions above, Equation (S.90) reduces to the following equation: M 00 Y L+0 = 0 (S.91) It can be verified by direct substitution that Equations (S.95) and (S.96) reduce to Equations (S.83) and (S.84), respectively, for the case 3τ = T . However, the notational and conceptual utility of the substitutions described above for simplifying the writing is evident.
Solution for any U (τ, T ) ≥ 1 For any U (τ, T ) ≥ 1, we follow the same steps to obtain U (τ, T ) equations: In S15 Fig, for two different values of τ , we plot the infection level versus time obtained from simulations of the dynamics. In both cases, we also show the growth rate,r 0 , of the infection level as predicted from Equation (S.99) when the infection level is sampled at integer multiples of the drug-dosing period, T . We find excellent agreement between theory and simulation.
With the same reasoning described in the previous two sections, we define a parameter,R 0 , inspired by the basic reproductive ratio R 0 that can be calculated for the case of constant β (Eqs.