Correction: Complex Dynamics of Virus Spread from Low Infection Multiplicities: Implications for the Spread of Oncolytic Viruses

[This corrects the article DOI: 10.1371/journal.pcbi.1005241.].


Part B. Mathematical details Background models
We will first consider two widely used models of viral infection (Eqs. [1] and [2]) [Perelson andRibeiro, 2013, Nowak andMay, 2000]. These two models distinguish between two types of cell populations: uninfected cells, x 1 , and infected cells, y 1 . The infection rate is , the death rate of uninfected cells is d and the death rate of infected cells is a (where a d).
Model (1): The first model (Eq. [1]) assumes that uninfected cells are produced with a constant rate . This assumption is used to describe in vivo infection dynamics where target cells (x 1 ) are produced by other types of cells not directly modeled in the equations, or when the number of uninfected cells is large and the quantity x 1 represents a smaller subpopulation of these cells.
⇢ẋ 1 = dx 1 x 1 y 1 y 1 = x 1 y 1 ay 1 [1] Model (S1): The second model (Eq. [2]) describes a closed system that assumes density-dependent target cell division. This type of model is appropriate to describe the dynamics of in vitro colonies of cells. Note that adenoviruses lock cells into the S-phase of the cell cycle preventing them from undergoing mitosis, hence in the model only uninfected cells are able to divide. If r is the division rate and K is the carrying capacity of the system. Model S1 is described by the following equations: ⇢ẋ 1 = rx 1 1 x 1 +y 1 K dx 1 x 1 y 1 y 1 = x 1 y 1 ay 1 [2] Stability and equilibria of the background models The experimental data discussed in the main manuscript describes a biological system where there is more than one stable outcome: weak self-limiting infection vs. robust and persistent viral spread.
It is then relevant to analyze the stability of the classical models of viral spread (1) and (S1). We find that both models have at most one non-negative stable steady state.
The steady states of system [1] are P 1 = ( /d, 0) and P 2 = (a/ , ( ad + )/(a )). Note that P 2 is non-negative if and only if 1 = ad + 0. However, if 1 = 0, then P 1 = P 2 . Hence, if P 1 and P 2 are di↵erent and non-negative we must have 1 > 0. On the other hand the eigenvalues of the Jacobian of [1] evaluated at P 1 are d and a + /d. The second eigenvalue is positive when 1 > 0. Hence, P 1 and P 2 cannot be both di↵erent, non-negative and stable.
Similarly, the steady states of system [2] are Q 1 = (a/ , ( ar + K( d + r))/( 2 K + r)), Q 2 = (0, 0) and Q 3 = (K dK/r, 0). First, note that the steady state number of cells when there is no infection is K(1 d/r) and it follows that for this model to be biologically meaningful we must have d < r; however, one of the eigenvalues of Q 2 is r d and thus the trivial steady state, Q 2 , is never stable. Next we turn our attention to Q 1 and Q 3 . The point Q 1 non-negative if and only if 2 = ar + K( d + r) 0. However, if 2 = 0, then Q 1 = Q 3 . Hence, if Q 1 and Q 3 are di↵erent and non-negative we must have 2 > 0. Now, one of the eigenvalues of Q 3 is 2 /r and thus it is impossible to have Q 1 6 = Q 3 , both points non-negative, and Q 3 stable. It follows that system (2) has at most one positive stable steady state 1 .

Modeling the antiviral state
To incorporate the antiviral state in the models we begin by assuming the simplest scenario, i.e., cells in antiviral cannot be infected. Let us call cells in antiviral state x 0 , the rate of antiviral induction , and the rate of of reversal from the antiviral state g. Then by modifying models (1) and (S1), we find the models (2) and (S2).

Model (2):
8 < :ẋ Once again we find that models (2) and (S2) have each at most one non-negative stable steady state. The analysis proceeds in a similar way to the discussion found in the previous section. Model (S2) has three steady states: Q 1 = (0, 0, 0), Q 2 = (K(1 d/r), 0, 0) and Q 3 = (x 1 ,x 0 ,ỹ 1 ), wherẽ x 0 = a r( a + K(1 d/r)) (a r + 2 K(d + g) + (d( K + r) + gr)) [6] y 1 = r(d + g)( a + K(1 d/r)) a r + 2 K(d + g) + (d( K + r) + gr) [7] Since systems [2] and [4] are identical when there are no infected cells, it follows again that we must have r > d and thus, since one of the eigenvalues of Q 1 is r d, this point cannot be stable. From the expression forx 0 in [6] and the definition of Q 2 , it follows that Q 3 is non-negative and di↵erent from Q 2 if and only if a + K(1 d/r) > 0. However, a + K(1 d/r) is also an eigenvalue of Q 2 and hence, Q 2 and Q 3 cannot both di↵erent, non-negative and stable. A similar argument proves that model (2) has at most one non-negative stable steady state.

1
In fact, it is easy to prove that each systems has exactly one non-negative stable steady state when 1 and 2 are di↵erent than zero.

Saturation of anti-viral defenses by multiple infection
The previous analysis suggests that multiple stability in the model is di cult to obtain if antiviral cells are completely immune to infection. Another possibility is that cells in the antiviral state can still become infected, but the antiviral state inhibits virus production. In this case multiple infection by more than one virus particle can increase the amount of virus produced and eventually saturate the intracellular defenses, causing multiply infected cells in the antiviral state to become infectious. The number of viruses particles that need to enter the cell before the antiviral defenses are overridden is not known. Hence, the first and simplest assumption to investigate is that infection by two or more virus particles overrides the antiviral state. If we adapt models (2) and (S2) to incorporate these set of assumptions we arrive at equations [8] and [9], where y 0 is the number of cells in antiviral state infected by just one virus particle. Note that y 0 cells are not infectious, but once they are infected a second time they turn infectious and behave like y 1 cells.
We first focus on model (3). The system's trivial steady state is: P 1 = ( /d, 0, 0, 0). This point is asymptotically stable if and only if the basic reproductive ratio of the virus R 0 = ( )/(da) is less than one. Let us know look for the other steady states. First, if for a steady state P we have y 1 = 0, it is straightforward to see that P = P 1 . Hence, we can assume that y 1 6 = 0 and settinġ y 1 = 0 in Eq. [8] we have: Using the previous equation for y 0 and settingẏ 0 = 0 we have: Substituting these expressions for y 0 and y 1 into the equation forẋ 1 = 0, we find: Note that the previous equations are not valid when x 0 y 0 = a/ + (x 0 + x 1 ) = 0. However, if this equality holds, it is straightforward to verify that there is a steady state only if y 1 = 0, in which case we simply recover the trivial steady state P 1 .
Expression (13) is a quadratic equation on x 0 and x 1 , which can be rewritten as: where, A 1 = 2 (a d)+ a , B 1 = 2 (g d), C 1 = 2 g, D 1 = 2 + a(d a) a 2 , E 1 = 2 ag and F 1 = a . Similarly substituting (10) and (11) into the equation forẋ 1 = 0 we find the condition: where, A 0 = a , B 0 = 2 ( a + d + g), C 0 = 2 (d + g), D 0 = a 2 , E 0 = ( a + d + g)a and F 0 = 0. Hence, the non-trivial steady states lies in the intersection of the two conic sections (14) and (15). Two conic sections intersect in 0,1,2 or 4 points, which gives the possible number of non-trivial steady states of the model. An explicit solution for the intersection of two conic sections can be written down. However, the formulas are too unwieldy to be of practical use. Instead we use numerical methods to find and classify the non-trivial steady states according to their stability (see the Numerical Method section in this Supplementary Materials).
We find that depending on the values of the parameters, model (3) can have multiple nonnegative stable steady states. That is, some sets of parameter values result in system [8] having more than one non-negative stable steady state; while other parameter sets result in system [8] having just one non-negative stable steady state. See Figure 4 in the main manuscript for an instance where the model is bistable.
Next we turn our attention to model (S3). In this case the trivial steady states are P 0 = (0, 0, 0, 0) and P 1 = (K(1 d/r), 0, 0, 0), where again we have d < r. Since one of the eigenvalues of P 0 is r d it follows that this point is not stable. It is also easy to see that P 1 is asymptotically stable if and only if K(1 d/r) < a/ .
To analyze the non-trivial steady states, we first note that systems [8] and [9] di↵er only in their equations forẋ 1 . It then follows that equations [10], [11] and [13] hold for both models. Thus, the non-trivial steady states of model (S3) lie in the intersection of the conic sections defined by [15] and the equation: where, H 1 = 2 r, H 2 = K(a( + ) + ( d + r)), H 3 = 2 r, H 4 = 2 K( d + g + r), H 5 = aK(a( + ) + ( d + r)), H 6 = 2 gK and H 7 = a gK. Once again it is possible to write explicit solutions for the steady states, but the expressions are to unwieldy to be useful. Instead we use numerical methods to analyze the stability and non-negativity of the non-trivial steady states (see methods section).
We find that depending on the values of the parameters, model (S3) can have multiple nonnegative stable steady states, or just one such steady state. Figure B(i) depicts an instance where the system is bistable. With the same set of parameters, but with di↵erent initial conditions, the trajectories converge to two di↵erent equilibria. In this figure, when the initial number of infected cells is four (plot labeled with "1") the infection first takes o↵, but eventually stalls and regresses to a very low level; when the initial number of infected cells is five (plots labeled with "2"), the infection is significantly more robust and it persists a↵ecting a large fraction of the cell population.

Stochastic formulation
In the experiments the di↵erent infection outcomes (limited vs robust infections) occur under identical conditions. Hence, we investigate if a stochastic formulation of the previous models can produce the di↵erent outcomes under identical initial conditions. The idea being that early stochastic events can push the system into either one of the domains of attraction of each of the two infection outcomes. We implement the stochastic formulations in the standard way using Gillespie's algorithm [Gillespie, 1977]. Figure B(ii) depicts two stochastic simulations of model (S3), both of which use the same set of parameters and identical initial conditions. In the first simulation (plots labeled with "1") the infection first takes and at remains at a low level for a long time until it eventually goes extinct. In the second simulation (labeled with "2") the infection is never extinguished, instead the total number of cells significantly diminishes and infection persists at high levels. It is important to note that the infection extinction outcomes are not in general the result of stochastic fluctuations around the initial very low numbers of infected cells. This observation is verified by Figure  B(iii). This figure present the distribution of the maximum number of infected cells for simulations where the ultimate outcome was viral extinction. Note that although the initial number of infected cells is very small (four cells) at their maximum extension limited infections average approximately 80 cells. Analogous results for the stochastic version of model (3) are found in the main manuscript.
As explained in the main manuscript, the inhibition of IFN signaling in our experiments resulted in a shift in outcome towards a prevalence of robust infections vs. limited infections. The same qualitative behavior is reproduced by the stochastic version of model (S3). Figure B(iv) plots the probability of the emergence of a robust infection as a function of the relative strength of the antiviral induction rate, . As this figure indicates, the reduction of the strength of the antiviral induction rate, which would occur as a consequence of the down modulation of IFN signaling, results in an increased probability of establishing a robust infection. The same qualitative behavior occurs for model (3) (discussed in the main manuscript).

Metapopulation models
To study the dynamics of the multiple infection and antiviral response in a spatial setting we consider a meta population approach. The model consists of a collection of n local patches. Within individual patches, local dynamics occur according to mass-action rules. The patches are coupled to each other by populations migrating between them. Here, we consider one-dimensional meta population models where populations in a given patch can only migrate to the nearest patches. In the equations k is the carrying capacity of each patch and cell populations with subscript i refer to cells in the ith patch. We assumed that the patches are arranged in a 1-dimensional linear array, and both target cells and infected cells can migrate to the neighboring patches to the left and to the right of a given patch with migration rate equal to m. The stochastic implementation of the metapopulation models follows directly from the reformulation of models (4) and (S4) using Gillespie's method.
The metapopulation model (4), which assumes a constant production of uninfected cells is described in the main text. Here we present the metapopulation model, which assumes densitydependent target cell division, which we call model (S4). First, let us define the operator f (z [17] Figure

Two-dimensional stochastic agent-based model
The rationale and algorithmic description of the agent-based model is discussed in the main manuscript. In particular, panels E-G in Figure 4 Figure D(iii) depicts the distribution of the maximum number of infected cells for simulations where the ultimate outcome was viral extinction (limited infections). The distinction between limited and robust infection outcomes is more di cult to make in the agent-based model because we cannot use ODEs as a guide to inform us about the stability properties of the the system. Hence, we chose the the threshold of 1000 infected cells to distinguish between robust and limited infections. If the cell population reached 1000 infected cells it was scored as robust, if it did not reach 1000 infected cells before going extinct it was scored as limited. In the experiments limited infections peaked at less than 1000 cells before they started to regress. With these criteria we found that out of 10,000 simulations 4094 resulted in limited infections. The e↵ect of inhibiting IFN signaling was modeled through a reduction of the anti-viral induction rate ( ). In agreement with experiments the down modulation of IFN signaling results in an increased proportion of the robust infection outcomes (Figure D(iv)).

Implementation of the agent-based model
Next, we provide details on the implementation of the two-dimensional stochastic model. First, in the simulation we choose a lattice distance of two cells to define local neighborhoods (i.e., if a cell has a lattice coordinates (i, j) its neighbors have coordinates (k, l), where i 2  k  i + 2 and j 2  l  j + 2). Let us refer to all cells that can be infected as susceptible cells (types x 0 , x 1 , y 0 ), and to cells that can transmit the virus as infectious cells (type y 1 ). For each lattice location (i, j) we define four quantities: R ij , I ij , D ij and A ij . The algorithm proceeds as follows: I) If there is a type x 0 cell at location (i, j), make R ij equal to the number of unoccupied lattice point in its neighborhood, otherwise make R ij equal to zero. II) If there is an infectious cell at location (i, j), make I ij equal to the number of susceptible cell in its neighborhood, otherwise make I ij equal to zero. III) If there is an infected cell at location (i, j), make D ij = 1, otherwise make D ij = 0. IV) If there is an infectious cell at location (i, j), make A ij equal to the number of type x 0 cells in the neighborhood, otherwise make A ij equal to zero. Now, if we define a = rR ij + I ij + dD ij + A ij , then the probability that at location (i, j) the next reaction will be cell death is dD ij /a, and the probability that it will be reproduction is rR ij /a. Similarly, the probability that the next event is infection by a cell in location (i, j) is I ij /a, and the probability that it is the conferral of the antiviral state is A i,j /a. In case of reproduction, the placement of one of the two daughters cells is chosen randomly from the unoccupied locations in the neighborhood of the dividing cell. Finally, if a cell at location (i, j) either infects or confers the antiviral state to another cell, the target cell is chosen randomly amongst the suitable cells in the neighborhood. To iteratively select the time and identity of the next reaction, we build and indexed priority queue and a reaction dependency graph as outlined in the so called "Next reaction method" [Gibson and Bruck, 2000]. This method has a O(log(m)) time complexity, where m is the number of lattice sites. Thus the method significantly improves the computational performance of the algorithm used in our previous study, which had a O(m) time complexity. The algorithm was coded in C using subroutines from GNU Scientific Library.

Numerical Method
As we previously mentioned, we can write closed analytical expressions for the steady states of models (3) and (S3), however they are too unwieldy to be of practical use. Hence, we need a robust and e cient numerical method to establish whether a parameter set yields a system with multiple stable equilibria or not. This is particularly important to be able to conduct thorough parameter explorations as the ones conducted in the next sections.
We now describe a method to find and score the equilibria of models (3) and (S3). We will explain the method for system [8]. We have shown the steady states must lie in the intersection of the two conic sections described by equations [14] and [15]. Note that for a fixed x 0 equation [15] is a quadratic equation on x 1 , which will have real solutions if and only if its discriminant It is easy to see that for non-negative steady states we must have 0  x 0  /d. Since is a quadratic function of x 0 , on the interval [0, /d], the condition 0 defines a set of points D, which is equal to either: i) a closed interval, ii) the union of two disjoint closed intervals, iii) a single point, or iv) the empty set. On the set D equation [15] defines two branches of x 1 as a real-valued functions of x 0 , corresponding to the two quadratic solutions for x 1 of [15]. On each branch of x 1 we can then define a function g : D ! R where g(x 0 , x 1 (x 0 )) is equal to the righthand side of equation [14]. If they exist, the non-trivial steady states must be zeros of the functions g. Since, for the non-degenerate cases D is either a closed interval or the union of closed intervals and, the functions g are continuous and smooth in the interior of D, it is straightforward to adapt a conventional one-dimensional search algorithm (e.g., bisection, interpolation, newton's method, etc) to find all the roots. Finally, to evaluate if the roots are non-negative exponentially stable steady states of [8], we evaluate x 1 , y 1 and y 0 for non-negativity and compute the eigenvalues numerically.

Non-negativity of solutions
To establish that the models are well posed from a biological perspective it is important to demonstrate the non-negativity of the solutions. We do so here for models (3) and (S3). The non-negativity of all other ODE based models follows directly from the same type of arguments presented here.
• Case 4 (y 0 (t) = 0). The analysis is the same as in Case 3: If y 1 (t) = 0, there is nothing to prove; if y 1 (t) > 0, the definition oft as the supremum of S is violated.
• Case 2 (y 1 (t) = 0). Following the proof of Proposition 1 (Case 2), if y 1 (t) = 0, then y 1 (t), y 0 (t) and x 0 (t) are all non-negative for t t . Next, if x 1 (t) 0 8t t we have z(t) 0 for t t and there is nothing left to prove. Suppose there is a t ⇤ >t such that x 1 (t ⇤ ) < 0, then the set R = {t t : x 1 (s) 0 for s 2 [t, t]} is bounded and non-empty and has a least upper boundt. Since x 1 is continuous we must have , which contradicts the definition oft as the supremum of R. We must have then x 1 (t) 0 8t t and it follows that z(t) 0 for t 0.
• Case 3 (x 0 (t) = 0). If either x 1 (t) or y 1 (t) are equal to zero, then from Cases 1 and 2 it follows that z(t) 0 for t 0. Furthermore, from the proof of Proposition 1 (Case 3), x 1 (t) and y 1 (t) cannot both be positive when x 0 (t) = 0. It follows that z(t) 0 for t 0.
• Case 4 (y 0 (t) = 0). The analysis follows the previous arguments. If y 1 (t) or x 0 (t) are zero, then from Cases 3 and 4, z(t) 0 for t t . If x 0 (t)y 1 (t) > 0, then the definition oft as the supremum of S is violated. Hence, z(t) 0 for t 0.

Parameter units
It is valuable here to discuss the units of each of the parameters. First, r, a and d have units of 1/time; K has units of cells; has units of cells/time; has units of 1/(time⇥cells); and has also units of .1/(time⇥cells) For comparison purposes we can rescale , , and so that every parameter has either units of cells or units of 1/time. Let L = /d be the equilibrium number of cells prior to infection in models (1),(2), (3), and (4), then = Ld and if we make¯ = d, the rescaled parameter (¯ ) has units of 1/time. Similarly, if we set¯ = L and¯ = L, then the rescaled parameters¯ and¯ have also units of 1/time. In the same way for model (S1), (S2), (S3), and (S4) we can set¯ = K and ¯ = K. We can then substitute in (1) the parameters and in terms of¯ ,¯ and L. Clearly the models can be rewritten in terms of these rescaled parameters. For example, model (3) becomes: 8 > > > > > > < > > > > > > :ẋ

Parameter exploration
We have demonstrated that models (3) and (S3) can have multiply stable steady states depending on the choice of parameters. It is then important to investigate the parameter regions where multiple stability might or might not occur.

Model (3)
First note that the units of time are arbitrary, therefore we can fix the value of one of the parameters and vary the others in relation to that fix parameter. In the simulations we chose to make = 0.001. It is also trivial to verify that the stability and number equilibria in system [18] does not depend on the value of L. For the simulations we choose a value of L = 1000 (note that we have then¯ = 1). Next we generated 400,000 random sets of parameters (a, d, , g) and color coded them according the number of stable equilibria. 2 Figure E depicts the projections in the six planes determine by the four free parameters. The shape of projection in the da plane is caused by the restriction d < a imposed on the model. Points in yellow represent parameters for which virus extinction was the only stable steady state; points in blue correspond to parameter where the single stable equilibrium included virus persistence; red points correspond to parameters that produce two stable equilibria (bistability); and finally, green points, of which they are very few, correspond to parameters where there was no stable steady state and the system instead displayed oscillatory behavior (see Figure F(i)).
The basic reproductive ratio of the virus R 0 =¯ /a is drawn for the da, a, and ga planes ( Figure E). Interestingly R 0 provides some information on the behavior of the system. For example, we did not find a single parameter in the region defined by R 0 < 1 where virus persistence occurred. However, a reciprocal strong statement cannot be said about the region defined by R 0 > 1. Indeed as this figure shows, all possible outcomes with regards to the persistence or extinction of the virus can occur in the region R 0 > 1. (Note that these results correspond to the ODE model and hence the outcomes are not stochastic.)

Model (S3)
To conduct a parameter exploration of model (S3) we again fixed = 0.001 and made K = 1000. The methods and color coding used was the same as for model (3). However, in this figure we distinguished an additional outcome related to the extinction of the target cell population. In the ODEs the number of infected cells cannot reach zero, but depending on the parameters, pronounced oscillations in the cell population can drive the number of target cells below one, which in a biological scenario would imply their extinction. For this reason we included a new color code, grey, to represent those parameters (with a single stable equilibria) that resulted in oscillations that pushed the target cell population below one before the infected cell population reached the same threshold, starting from the initial conditions x 1 (0) = K and y 1 (0) = 3 (consistent with our paper's focus on low infection multiplicities). We did not included this option (target cell extinction) in model (3) because the constant rate of production of uninfected cells in [8] does not allow for the extinction of target cells. Figure G presents the results for the parameter exploration, where to limit the number of panels we focused our attention on the case d = 0 (no uninfected cell death).
Finally, an interesting observation, is that even though in both models ( (3) and (S3)) the number of steady states can number up to five, we did not find any parameter set that produced more than two positive stable equilibria.

Agent-based model
For the stochastic agent-based model we do not have the type of analytical tools that allowed us to classify the parameter sets for models (3) and (S3). For this reason in Figure H we only color code two outcomes: (i) limited infections (yellow), and (ii) robust infections (blue). A point in the parameter space was recorded as a robust infection if in the simulation the infected cell population reached 1000 infected cells. Conversely, points in the parameter space were recorded as limited if in the simulation the infection died out before reaching an infected cell population of 1000 cells. As before we fix the value of (in this case to = 1). Simulations were conducted on large grid of size 1000 ⇥ 1000 cells. In this case the reproductive ratio was defined as R 0 = k/a, where k = 25 is the size of the local neighborhood of each cell. With the same set of parameters, but di↵erent initial conditions the trajectories converge to di↵erent equilibria. In this panel, if the initial number of infected cells is equal to three (y 1 (0) = 3), the infection first takes o↵, but eventually regresses and stabilizes at a very low level (plots labeled with "1" in the figure). If the initial number of infected cells is equal to four the system converges to a persistent and robust infection (plots labeled with "2"). (ii) Bistability in the stochastic mass action model. This panel shows two simulations with identical initial conditions. The simulations results in either a weak, limited infection that is eventually extinguished (labeled with "1"), or in a persistent and robust infection (labeled with "2"). (iii) Distribution of the maximum number of infected cells when there is a limited infection. Although the initial number of infected cells in this panel is very small (y 1 (0) = 3), at their maximum extension limited infections average approximately 80 infected cells. (iv) Probability of the emergence of a robust infection as a function of the normalized level of the the antiviral induction rate . The reduction of the strength of the antiviral induction rate results in an increased probability of establishing a robust infection (1000 simulations per antiviral induction level, error bars represent 95% confidence interval). In panel (i): r = 10, d = 0, = 8.94 ⇥ 10 4 , = 1.92 ⇥ 10 2 , g = 4.76 ⇥ 10 3 , a = 0.136. In panels (ii), (iii) and (iv): r = 40, d = 0, = 1.45 ⇥ 10 4 , g = 2 ⇥ 10 2 , a = 2 ⇥ 10 2 . In panel (ii) the value of = 2 ⇥ 10 2 . In panel (iii) the value of = 2.34 ⇥ 10 2 . In panel (iv) the 100% level of antiviral induction rate corresponds to = 2.88 ⇥ 10 2 . In all panels: x 1 (0) = 1000, and x 0 (0) = y 0 (0) = 0. In panels (ii),(iii), and (iv): y 1 (0) = 3.             Figure H