Towards Predictive Computational Models of Oncolytic Virus Therapy: Basis for Experimental Validation and Model Selection

Oncolytic viruses are viruses that specifically infect cancer cells and kill them, while leaving healthy cells largely intact. Their ability to spread through the tumor makes them an attractive therapy approach. While promising results have been observed in clinical trials, solid success remains elusive since we lack understanding of the basic principles that govern the dynamical interactions between the virus and the cancer. In this respect, computational models can help experimental research at optimizing treatment regimes. Although preliminary mathematical work has been performed, this suffers from the fact that individual models are largely arbitrary and based on biologically uncertain assumptions. Here, we present a general framework to study the dynamics of oncolytic viruses that is independent of uncertain and arbitrary mathematical formulations. We find two categories of dynamics, depending on the assumptions about spatial constraints that govern that spread of the virus from cell to cell. If infected cells are mixed among uninfected cells, there exists a viral replication rate threshold beyond which tumor control is the only outcome. On the other hand, if infected cells are clustered together (e.g. in a solid tumor), then we observe more complicated dynamics in which the outcome of therapy might go either way, depending on the initial number of cells and viruses. We fit our models to previously published experimental data and discuss aspects of model validation, selection, and experimental design. This framework can be used as a basis for model selection and validation in the context of future, more detailed experimental studies. It can further serve as the basis for future, more complex models that take into account other clinically relevant factors such as immune responses.

5. The growth has to be saturated in both x and y, such that lim x→∞ yG(x, y) = H x (y), 0 < H x (y) < ∞, where H x (y) is a function of y independent of x. Similarly, with y.
6. For large values of x, the growth term cannot be positive in the limit of small y, that is, lim yG(x, y) > 0.
Note that this expression could be infinite.
The analyses below rely on these properties.

The number of equilibria
The fixed points of the virus-cancer system (equations (1-2) of the main text) are given by (0, 0) and all the solutions of the equations G(x, y) = a β . (2) The trivial point (0, 0) has eigenvalues F (0) and −a and is thus a saddle. The number of solutions of equations (1-2) depends on the particular shapes of the functions F and G. In order to find the nontrivial equilibria, we solve equation (1) to find y(x), and then substitute it into equation (2). The equilibria are thus defined by the roots of equation βG(x, y(x)) = a.
From equation (1) we can see that y(0) = 0. We know from assumption (2) on the function G that G(0, 0) = 0. The next step is to study the limiting behavior of G(x, y(x)) for large values of x. For that we need to know the behavior of y(x) for large x. We have from equation (1): There are three cases. (i) For a linear type growth, we have lim x→∞ xF (x) = c 0 , a nonzero constant. In this case, lim x→∞ y(x) = c 0 /a, with 0 < c 0 < ∞.
(ii) For any growth F which is superlinear but slower than exponential, we have lim x→∞ y(x) = ∞, but lim x→∞ y/x = lim x→∞ F (x + y)/a = 0, that is, y increases slower than x. (iii) Finally, for exponential growth, F = 1 and y(x) = x/a, such that y(x) ∼ x for large values of x. From the biological assumptions on the function G(x, y) listed above, it follows that for any of the possible dependencies y(x), the function G(x, y(x)) approaches a finite limiting value as x → ∞, and this value can be zero or nonzero. We use the exponential case F = 1 and y exp (x) = x/a to separate all functions G into two classes. If lim x→∞ G(x, x/a) = 0, then we will consider the virus spread to be slow, and if lim x→∞ G(x, x/a) = G ∞ exp > 0, with G ∞ exp < ∞, we will regard this as fast spread. Note that for all laws of cancer growth slower than exponential, we have G(x, y(x)) ≥ G(x, y exp (x)). This is because y(x) ≤ y exp (x), and G is a decreasing function of y.

Fast virus spread
In this case, the function G exp (x) ≡ G(x, y exp (x)) is either a monotonically increasing function (figure 1(a) of the main text), or it can attain one or more local extrema before converging to its nonzero horizontal asymptote.
For a monotonically increasing G exp , low values of β correspond to zero roots in equation (2), which means that the cancer growth will continue indefinitely (figure 1(a) of the main text). As β crosses a critical value defined by a/G ∞ exp , there is one root. The value of x at this root drops as β increases (this is due to the convergence of G exp to an asymptote, G ∞ exp ). For large values of β, the value of x at the intersection tends to zero.
If the function G exp has an absolute maximum at point x exp , then an initial increase of β above a/c max with c max = G exp (x exp , y(x exp )) results in the appearance of two roots. Additional local extrema will result in an appearance and disappearance of pairs of roots. However, as β increases through the second threshold given by a/c 2 , only one (the lowest) root remains. The value of c 2 is given by the lower of the values {G ∞ exp , c min }, where c min is the value at the lowest local minimum, if it exists.
In both cases, for sufficiently large values of β, there will be only one root in equation (2). Introducing other cancer growth laws can increase the limiting value of G thus decreasing the value of β c . In the case of a monotonically increasing G exp , there will be no qualitative change. If G exp is one-or multiplehumped, the hump(s) may disappear. Whether this qualitative change happens depends on the relative size of the two spacial scales involved. The first scale is defined by the location of the maxima of G exp and is related to the virus spread scale, s v . The second scale is given by the size, s t , at which cancer growth law starts to deviate from exponential. Once s t ∼ s v , the limiting value of G becomes sufficiently large such that the "hump" disappears.
It is useful to investigate the value of x at the equilibrium as a function of β, for different values of s t . Suppose that the graph of G(x, x/a) is a monotonically growing function of x which approaches a limiting value, G ∞ exp . Suppose that the cancer growth slows down around the scales near s t . So near x ∼ s t , the function G(x, y(x)) deviates from the horizontal asymptote, G ∞ exp , and starts growing toward a different, and higher horizontal asymptote, which we will call  figure 1 for a particular example). The phase diagram as β increases is as follows. For β < a/G 1 , there are no roots. As β crosses the first threshold, a/G 1 , one root appears. The value of x at this equilibrium decays rapidly from infinity to values around s t , as β grows (because of the fact that G 1 is a horizontal asymptote). Then as β grows through its second threshold, a/G ∞ exp , the value of x at the equilibrium drops from s t to values of order 1. The second transition is sharp if the following is satisfied: In other words, x 1 is the value of x where the function G(x, y(x)) comes near its first horizontal asymptote ("near" means that it is at least as close to G ∞ exp , as G ∞ exp is to G 1 ). If s t x 1 , then the function G has a significant interval in x where it comes near the value G ∞ exp , before it deviates from it to start growing toward G 1 . This guarantees a threshold effect.
We conclude that for all cancer growth laws and for all functions G corresponding to a fast virus spread, increasing β beyond a threshold leads to the existence of only one equilibrium, whose value correlates negatively with the viral replication rate, β. For large enough s t , there is a "threshold" effect, such that the size at equilibrium decreases very sharply as β approaches a defined value.

Slow virus spread
In this case, the function G exp (x) ≡ G(x, y exp (x)) is a one-or a multiple-humped function, which approaches zero for large x (figure 1(b) of the main text).
In the case of an exponential growth, the bifurcation diagram looks as fol-lows. As before, small values of β correspond to no equilibria (zero roots in equation (2)). As we increase β, a pair of roots appears after the threshold given by a/c max , where c max = G(x exp , y(x exp )), the maximum value of G exp . As β increases, other roots may appear and disappear in pairs. Since G exp (0) = 0 and the function G exp has zero as its horizontal asymptote, there will be two equilibria for all values of β larger than a/c min , where c min is the value of G at its lowest local minimum in the case that such a minimum exists; c min = c max otherwise. Let us next consider how non-exponential laws of cancer growth modify this picture. In some cases, a slower-than-exponential growth term, F , will lead to the horizontal asymptote of G(x, y(x)) becoming nonzero. In general whether this happens depends on the functional forms of both G and F . In the case of linear growth, y(x) ≡ y lin (x) converges to a nonzero constant, c 1 , and we have lim x→∞ G lin (x) = lim x→∞ G(x, c 1 ) = G x (c 1 ) = c 2 < ∞, which is a nonzero constant. Depending on the value of s t , G lin can be a one-or a multiple-humped function, or (for s t similar or smaller than x exp ) it is a monotonically increasing function of x. In either of these cases, there exists a finite value of β such that for all values of β larger than this value, there is only one root in equation (2). For growths faster than linear but slower than exponential, we have y → ∞ as x grows, but y = o(x), i.e. it grows slower than x. In some cases the function G will retain a zero asymptote (e.g. in the case where G = x/(x + 1)/(y + 1) and a surface growth law for F ). In other cases it will acquire a nonzero limit (e.g. with G = x/(x + 1 + √ x(y + 1)) and a surface growth law for F ). In the latter case we can say that the surface cancer growth is sufficiently slow to warrant successful treatment given the particular mode of viral spread.

Finite tumor size
Here we consider a growth term which becomes zero for a finite value of x + y. The growth starts off exponential (F (0) = 1) and at some size, s t , it slows down (we do not exclude the possibility that s t ∼ 1, that is, the growth becomes slower than exponential right away). Then there exists another characteristic size, W s t such that the growth slows down further and stops. We define W such that F (W ) = 0. Note that if s t ∼ W then there is no need to introduce the two scales, s t and W . Therefore, the assumption s t W must hold. The previous analysis holds on the scales intermediate between s t and W , such that s t x W . In particular, for values x W , the shape of the curve G(x, y(x)) is similar to that obtained for the corresponding unlimited growth. As x grows far beyond s t and approaches W , the function G approaches G(W, 0). If, for the unbounded growth, the limiting value of the G function is c 2 , we have in general G(W, 0) ≥ c 2 , and the curve G takes an upward turn in the vicinity of x = W . This means that equation (2) acquires an additional root corresponding to the cancer growing to its carrying capacity, W . In the systems with unrestricted growth this was equivalent to an unlimited growth of the cell population.
It is useful to note the following: in systems with a limited size, the function G(x, y(x)) is always bounded away from zero. Therefore, strictly speaking, we can always find a threshold value β t such that for β > β t , only one root is present. However, if W s v , such values of β are very large compared to β c , and in most cases are probably not achievable.

Stability
Let us suppose that (x 0 , y 0 ) with x 0 ≥ 0 and y 0 ≥ 0 is a solution of system (1-2), and consider the stability of the corresponding equilibrium. The Jacobian of the system can be written as a 2 × 2 matrix, {m ij }, with where the functions F and G and their derivatives are evaluated at the point (x 0 , y 0 ): G x = ∂G/∂x| x=x0,y=y0 , and similarly with G y and F . The equilibrium is stable if the following two conditions hold:

Saddle points
Condition (5) is equivalent to the positivity of the derivative of G in the direction defined by the implicit relation ya = xF (x + y), equation (1). Differentiating it, we get: ady = F dx + xF (dx + dy). The directional derivative is equal to . The denominator is positive, so this expression has the same sign as the left hand side of condition (5). The equilibria are defined by the roots of equation (3). Since G(0, 0) = 0, all the odd roots of equation (3) will correspond to a positive, and the even ones to a negative slope of the left hand side of equation (3). This means that all even equilibria are saddles. This is because in such cases the directional derivative is negative, condition (5) is violated, and therefore there are two real eigenvalues of opposite signs. On the other hand, an odd root can be either a sink, a source or a spiral (stable or unstable). This is because for such a root, condition (5) is always satisfied, so that we could either have complex eigenvalues, or real roots of the same sign (positive or negative).
In the presence of a saddle, an infinite outcome (corresponding to an unchecked cancer growth) is possible. For large values of x, we havė where lim x→∞ G(x, y) = G ∞ (y) and lim x→∞ F (x, y) = F ∞ . The growth of y becomes negative as y increases if lim y→∞ G ∞ (y) = 0, which suggests that y settles to a finite value which makes the right hand side of equation (7) zero, such that the outcome (∞, const) is observed. If lim y→∞ G ∞ (y) = const > 0, then for large enough values of β we can have an outcome of the form (∞, ∞).

Properties of the internal equilibrium
Let us first show that for large values of β, there will be an equilibrium, (x 0 , y 0 ), such that lim β→∞ x 0 = 0 and lim β→∞ y 0 = 0. We call this equilibrium the "internal equilibrium". Its existence follows from equation (3) and the properties of the function G. We know that y(0) = 0, and also that G(0, 0) = 0. It is also clear that there is an interval of x, [0, ξ], where G is a growing function. Therefore, by continuity, for all β ≥ a/G(ξ, y(ξ)), there will be a solution of equation (3). From monotonicity of the function G, the value of x at the intersection with a/β decays with β. From equation (1) it follows that there is an interval of x, [0, ξ 1 ], where y is a growing function of x. Therefore, we conclude that for large enough β, there is an equilibrium whose x and y values decay with β and approach 0 in the limit β → ∞. Let us evaluate the left hand sides of inequalities (4) and (5) for small values of x 0 and y 0 . First, we express β from equation (2): β = a/G(x 0 , y 0 ). Then we approximate the curve y(x) by its Taylor series for small values of x 0 : where the function F and its derivatives are evaluated at 0. This expression follows from expanding both sides of equation (1) in Taylor series in terms of x 0 and y 0 , solving for y 0 and using a Taylor expansion of this expression. Now, let us multiply the left hand side of inequality (4) by G(x 0 , y 0 ), and use expression (8). Expanding in terms of small x 0 , we obtain: Here the functions F and G and their derivatives 2 are evaluated at zero. To derive the above expression we also used the fact that the function G and its y-derivatives are equal to zero if x = 0, and F (0) = 1. Next, we evaluate the left hand side of inequality (5) in the same manner: We can see that the expression above is always positive, so condition (5)   as follows from expression (9). The expansion can be positive or negative, depending on the particular properties of the functions F and G. Next, we would like to investigate whether the eigenvalues are real or complex. For the eigenvalues to have an imaginary part, the following condition has to be satisfied: (m 11 − m 22 ) + 4m 12 m 21 < 0.
Performing a Taylor expansion of the above expression for small values of x 0 and y 0 at internal equilibrium, we obtain: We can see that this quantity is always negative. Therefore, we conclude that the internal equilibrium has complex eigenvalues for sufficiently large values of β.

Fitting experimental data
The experimental data of [1], figures 1 and 5, was used to fit with the models. The fitting was performed by using Mathematica, which allows to perform nonlinear least squares fitting procedure on solutions of ODEs defined on a set of unknown parameters. The root mean square (RMS) error was calculated as Here n is the number of experimental data points in the time-series, and (t i , x i ) are the coordinates of the experimental points, where t i values correspond to time and x i correspond to the cancer size, in mm 4 . The fitted function obtained as a solution of ODEs is denoted as f (.). First we used various growth laws to fit with the data of [1] (figure 1), where the growth of uninfected cancer cells was measured. The results are presented in figure 6(a) of the main text; the best-fitting parameters are listed in Table 1. We then used the Logistic growth model to fit the data pertaining to cancer-virus interactions.
The virus growth terms used were G(x, y) = x/(x + y + ) for the fast spread model and G = x/(xy 1/3 + ) for the slow spread model (both contain only one parameter). Thus, there are 4 unknown parameters to fit: , β, a, and f 0 , the initial fraction of infected cells at the time when the first measurement took place. Some examples of well-fitting parameter sets are presented in table 2.