The paper deals with mathematical modeling of a novel mode of cancer therapy whereby a tumor is treated with oncolytic viruses that have a selective affinity to cancer cells. The model is couched as a system of coupled differential equations for the sizes of populations of uninfected an infected cancer cells. The main thrust of the work is on introducing parametric heterogeneity and analytic/numerical study of its effects on the dynamics of tumor cell populations. The problems in questions are of significant theoretical and practical interest. The authors suggested a few compelling general ideas, demonstrated a good knowledge of both biomedical and mathematical aspects of the general field, and provided ample references.
However, the paper concerns me on several levels.
1. The model was derived not from a detailed theoretical or experimental study of the underlying biological phenomena but rather from "general considerations." Also, it was not fitted to any set of real data on the sizes of tumor cell populations interacting with oncolytic viruses in vivo or in vitro. Thus, the model lacks biological specificity, and its adequacy has not been established. In particular, is there any strong evidence of the assumed heterogeneity of tumor cell populations with regard to killing action of oncolytic viruses, and does this assumption improve the quality of the model and its predictive power?
: The model was developed to describe qualitative behaviors of the system "tumor cells – oncolytic virus", and we by no means claim that this model should be used for quantitative predictions. The current empirical data is insufficient to obtain robust estimates of the model parameters. These points are addressed in greater detail in our previous paper in this series (Ref
The hypothesis that tumor cell populations are heterogeneous with regard to the killing action of oncolytic viruses is a general one, and our model strongly suggests that the assumption of homogeneity(that is, that all the tumor cells behave in the same way under oncolytic virus infection) is an oversimplification. Indeed, there is no strong evidence that tumors are heterogeneous in terms of their response to oncolytic virus infection but this seems highly likely given the well-demonstrated heterogeneity of tumors in other respects, in particular, chemotherapy. We investigate the implications of heterogeneity for the outcome of oncolytic virus treatment and hope that our model stimulates experimental studies in this direction.
2. The way parametric heterogeneity was introduced has a fundamental flaw. Speaking about the distribution xβ of the susceptibility parameter β (that is, the rate of virus transmission from infected to uninfected tumor cells) in the population of uninfected tumor cells of size X, the authors assumed that the distribution of β is a function of time t alone. In reality, this distribution depends on the population size X:xβ = xβ (·, X) or more generally, xβ = xβ (·, X, t). For example, even in the simplest case of a uniformly distributed β, its density is given by X/|B|1B and thus depends on X. Omission of the dependence of the distribution of β on X is tantamount to assuming that every first order ODE is of the form y = f(t), where t is independent variable and y is unknown function, on the grounds that for each solution the right-hand side of such an ODE is a function of t. As a result, the model of distributed susceptibility, which in reality is non-linear, is trivialized to a linear model which, naturally, boils down to the time course of the expected value of parameter β. Treatment of the heterogeneity of other parameters (cytotoxicity and virulence) suffers from the same defect.
Author response: This is, in our opinion, the main objection of the reviewer to the subject and conclusions of our work, and the only one with which we completely disagree.
i) We never stated or assumed
"that the distribution of β
is a function of time t alone"; as a matter of fact, we do not make any explicit assumptions on the distribution dynamics except for the given initial distribution and the form of the model equations. Under the proposed mathematical formalism (Theorem 1 in MA [see Additional file 1]) and equations used to describe the system, we can unambiguously determine the evolution of the parameter distribution given the initial distribution. In particular, we never assumed that the parameter (β) is gamma-distributed at any time moment with a time-dependent mean and variance; on the contrary, we suppose only that β is gamma-distributed at the initial moment and prove that then it must be gamma-distributed at any time moment due to the system dynamics, and compute its mean and variance depending on time (Theorem 1 and Example 1, MA [see Additional file 1])
ii) We emphasize that xβ is not a probability density function, that is, its integral equals the total population size and not 1. Theorem 1 and its corollary show how one can calculate the density xβ , the total population size X
), the probability density function xβ
/X(t) and its moment generation function. To calculate all these quantities we use nonlinear systems ((A.4)–(A.6) in MA [see Additional file 1]) which evidently depend on the current population sizes in a complex way. The same is true for all other parameters and versions of the model
3. The model formulation contains arbitrary assumptions that are not justified and at times not even explicitly stated:
(a) The authors assume logistic growth of the sizes X and Y of infected and uninfected tumor cell populations, respectively. Why logistic and not Gompertz, which is generally believed to be more adequate?
Author response: This issue was considered in our previous work. In short, it is quite disputable that "Gompertz's growth is generally believed to be more adequate". We chose the logistic law because it is the simplest form whose predictions agree with the empirical data. Moreover, it seems most unlikely that qualitative results change if we use the Gompertz growth law instead of the generalized logistic law.
(b) The interaction between the two populations is assumed to be governed by the term
βXY/(X + Y). Why is the infection rate assumed proportional to the relative size Y/(X + Y) of the infected population rather than to its absolute size Y? Is there any theoretical or experimental rationale for such an assumption?
: This issue was discussed in detail in Ref
(c) The assumption that "transmission coefficient β is the product of the susceptibility of uninfected cells β1 and the virus replication rate β2" (Section Distributed susceptibility ...) is quite problematic because it contradicts the structure of the model. Imagine the virus replication rate is doubled. However, the term βXY/(X + Y) would not necessarily double! Therefore, model (11) is also incorrect. The above assumption would be plausible only if the interaction term in the model were βXY.
Author response: We would like to clarify that the system (11) is correct if the parameter β2 is attributed to any trait that describes the ability of the virus to infect tumor cells. To justify the use of model (11), we slightly changed the text, in particular, by replacing the term 'virus replication rate' with the more abstract term 'virus virulence' that reflects the general ability of a virus to infect tumor cells. In this setting, the model (11) is correct.
(d) The formula E(t) = Eβ1 (t) Eβ2 (t) on p. 17 means that the authors tacitly assumed that parameters β1 and β2, viewed as random variables, are independent. Are they? It seems likely that "virulence" of the virus may affect both virus replication rate β2 and susceptibility of uninfected cells β1.
: Yes, it was assumed that the two parameters are independent. We do not know anything about the dependence between these parameters although it is a natural extension of the model to assume some form of dependence. The independence assumption was chosen to simplify the math. More general situations are also amenable to mathematical analysis given the initial joint distribution (similar to Ref
4. Moment generating function MP(λ) of a probability distribution P is generally not defined for all real λ. This brings about artificial difficulties that are not properly accounted for in the paper. Also, why not to deal with the characteristic function instead and avoid these difficulties altogether?
: Importantly, the current parameter distribution is determined through the mgf of the initial distribution, and not through the characteristic function. Indeed, a mgf cannot be defined for all real λ, but this does not amount to "artificial difficulties". We would like to emphasize that, in some cases, the non-existence of the mgf reflects intrinsic and important problems. In brief, in the MA, it is shown that heterogeneous model (A.2) is equivalent to the Cauchy problem (A.4)–(A.5) if the latter has a unique global solution in
). Obviously, this solution does not exist if the corresponding mgf M
) does not exist for some λ = q
), and, consequently, solution of the heterogeneous model does not exist at this t (e.g., there can be a blow-up; for a simple example, see ref
]in the MA). This is now mentioned in the text (section on
5. The authors should draw a careful distinction between mathematical results (in which case they should be proved) and the results of their simulation studies.
Author response: All new mathematical results that are used in the text are now proved in the MA.
6. Bifurcation analysis and phase portraits of the homogeneous model are too sketchy for their validity to be evaluated. The authors should describe the types I-VIII of system behavior in more detail. Are all of them observed in reality? Also, what happens in the case γ = 1?
: Detailed bifurcation analysis of the homogeneous model is presented in our recent publication (Ref
]), so we did not perceive it necessary to reiterate the descriptions here
7. How were parameters for the numerical simulations selected? Are they representative of the outcomes?
Author response: The parameters for the numerical simulations were selected such as to illustrate new dynamical regimes of heterogeneous models which are unobservable in homogeneous settings. The parameters used are representative in the sense that a set of close parameter values yields qualitatively similar results.
8. The content of the appendix provides some mathematical background for the main text but does not seem to support it in more specific ways. In particular, asymptotic analysis that was discussed in the main text was not carried out in the appendix.
Author response: The corresponding changes were made in the main text (section on Distributed susceptibility) and in the Mathematical Appendix)
Reviewer's report 2
Natalia Komarova, Department of Mathematics, University of California-Irvine (nominated by Orly Alter)
The paper by Karev et al "Mathematical modeling of tumor therapy with oncolytic viruses: Effects of parametric heterogeneity on cell dynamics" creates a mathematical framework for studying dynamical systems with distributed parameters. This is used to model treatment regimes of cancer with oncolytic viruses. There, tumor heterogeneity is shown to play an important role. It is demonstrated that, depending on the distribution width of the parameters, very complex behaviors can be observed, including tumor dormancy, tumor recurrence and a reversal of treatment success.
I believe that developing a systematic modeling framework for systems with continuously distributed parameters is very important. I therefore recommend the paper for publication. However, I would like the authors to clarify the following points.
(1) The authors emphasized the importance of the variance of the initial distribution of the heterogeneous parameter. I believe that in this model, there is a characteristic of the initial distribution which is even more important than that. It is the support of the initial distribution of the parameter. Simply speaking, since there are no mutations in the model, knowing the types present in the system initially reveals with certainty what will happen in the end.
In fact, the behavior of the system with a distributed parameter can be predicted qualitatively from the following two pieces of information: (i) the direction of selection and (ii) the support of the initial distribution function. For instance, if the transmission coefficient, beta, exhibits heterogeneity, then we know that (i) cells with a lower beta values will be selected for and (ii) the lowest beta from the initial distribution is given by η. Therefore, we know that eventually the cells with β = η will outcompete the rest, and the dynamics of the system can thus be predicted.
Author response: We agree that the support of initial distribution is an equally important characteristic of the system. While the parameters of the initial distribution together with system dynamics determine the speed with which the final state is reached, the support determines, in the models considered in the text, the final outcome. We added some text to clarify this point.
In general, the assertion that the behavior of the system with a distributed parameter can be predicted given that the direction of selection and the support are known is not valid (although it is true for our particular models). The simplifying fact in our presentation was that the mean parameter values are monotone functions (i.e., we can unambiguously identify the direction of selection). In general, these functions depend on the system dynamics and can be arbitrary (e.g., it is possible to construct a system where they are periodic (Ref
(2) Intuitively speaking, the extreme values of the distributed parameter (e.g. η in the gamma-distribution used by the authors) must be responsible for the final outcome of the dynamics, and the width of the distribution should correlate with the speed with which the system attains its final homogeneous state. Is this true? Is this possible to demonstrate?
Author response: The final outcome in our models depends, as correctly stated in Komarova's remark (1), on the direction of selection and the extreme values of the distributed parameters. E.g., if we have a gamma-distribution on [η, ∞), then the value of η determines the final outcome if individuals with smaller parameter values are more fit, but is not essential for the asymptotic system dynamics in the opposite case.
We do not have an analytical solution to the question on correlation between the speed of movement in the parameter space and the width of the distribution. It is generally true, however, that the speed with which the system attains its final state depends on the width of the distribution and on the current parameter values, and the sizes of the cell populations.
(3) The elegant model presented by the authors deals with continuously distributed parameters. In reality, many parameters can only attain a discrete (and small) set of values. It would be very nice to see (and easy to show) how the mathematical analysis should be modified to deal with discrete distributions. Some of the consequences of the analysis may look more intuitive in this case.
Author response: We explicitly indicate in the revised text that the mathematical framework can equally be applied to continuous or discrete distributions. The mathematical analysis does not have be modified (discrete distributions were used, e.g., in (Karev, 2003, Ecol. Model. 160, 23–37). For example, if we assume that the initial distribution of the parameter β that describes susceptibility of uninfected tumor cells is Poissonian with initial mean β0 then, using mgf for the Poisson distribution, we obtain Eβ (t) = β0 exp(q(t)). The latter expression can be used in system (5)–(6).
(4) I found the mathematical derivation presented in the appendix a little hard to read. It would help if the authors showed that the distribution x(t, β) (normalized) is equal to the distribution exp(β q) p (0, β) (normalized) if the variables x and q satisfy the given equations.
: The proof of Theorem 1 contained some typos which might have led to confusion; this was corrected [see Additional file 1]
Reviewer's report 3
David Krakauer, Santa Fe Institute
In this paper Karev et al extend their earlier work on the dynamical properties of an implicit oncolytic dynamic (virus is not treated and assumed through a separation of time scales to be stationary) to include continuous variation in the viral traits transmission coefficients and cytotoxicity. The model comprises a system of ODEs tracking mean densities of uninfected and infected cells including time varying mean transmissibility and cytotoxicity.
In this more complex model with heterogeneity, the range of dynamical behaviors is expanded, manifesting patterns of tumor recurrence, and quasi-periodic orbits in cell densities. The paper raises interesting questions about the possibility of systematic, positive intervention into such a complex system.
While the distributional approach of this paper is timely and expands the range of model behaviors, it remains a continuous, deterministic treatment and the stochastic implication of rare events stemming from the tails of distributions can not be understood. And I would hypothesize that it is these improbable events that are more typically associated with recurrence expanding from small populations of concealed cells.
Author response: Indeed, stochastic implications of rare random events cannot be understood within the framework of this paper because the models are deterministic. We use a probabilistic distribution to describe the parametric heterogeneity of tumor cells. The tails of distributions considered in the paper show only that there are subpopulations of cells with relatively high or low parameter values, and the proportion of these subpopulations is small. We agree that "improbable events" implicated by Krakauer are important in carcinogenesis, and in cancer recurrence in particular. However, we would like to stress that our model shows that not only such events can lead to, e.g., tumor recurrence, but also existence of small populations of cells whose probability to be infected by an oncolytic virus is relatively small but still essentially non-zero.
A general problem I have with the paper, that will probably diminish its impact, is that no clearly important insights are generated by the model, but rather a range of complex behaviors which probably need further analysis. Since this paper offers qualitative results, I view this as a failing, as the purpose of such models is largely to sharpen intuition. One problem throughout is identifying clearly the source of heterogeneity – is this of viral or cellular origin? In the model this remains ambiguous.
Author response: Obviously, most biological systems are heterogeneous including the population of tumor cells. The whole point of this paper is whether a heterogeneous model allows to describe qualitatively new phenomena in comparison to the corresponding homogeneous model. The models described here possess, e.g., regimes of i) tumor recurrence, ii) transiently constant tumor size, iii) quasi-chaotic behavior. Thus, it seems that this work, actually, does sharpen intuition and even yields results that might not be intuitively obvious. Furthermore, we would like to emphasize that we propose a new modeling approach to deal with parametric heterogeneity, which is computationally and, in some cases, analytically feasible, and which can hopefully be applied to other existing models of ODEs for cancer progression and treatment.
The source of heterogeneity in our models can be of viral or cellar origin. When describing each specific model, we tried to explain possible sources of heterogeneity. For instance, when we deal with distributed susceptibility, the source is cellular. Viral heterogeneity is included in the models through the infected cell population because we do not treat the virus population explicitly. We agree that explicit virus dynamics would improve and clarify the models, but it also would complicate considerably the full parametric analysis of the homogeneous model.
I wonder if the authors could not shorten the paper and really focus on its most important insights. I think that these would include: (1) Increased resilience of tumors stemming from heterogeneity, (2) The increased success of oncolytic treatment with heterogeneous virus (if this is a real result of the model), (3) the impact of variance of heterogeneity on dynamics, (4) the impact of the distribution, (5) the clearly stated implications for therapy.
Author response: We believe that all these issues are addressed in the paper; we do not see how the article could be shortened significantly without losing information.
I also have some specific comments which overlap with those I made in a review of the author's previous paper in this series.
1. I feel that not explicitly treating the virus misses important problems and would make the interpretation of model parameters much clearer.
Author response: As indicated above, we agree but this would make the model substantially less tractable.
2. I am worried about the requirement that complete information about the distributions needs to be known in order to track the means in the ODEs. I am guessing that not only the shape but the form of the distributions could change over the course of infection.
Author response: Formally, to obtain a solution of the system at any time moment, including t → ∞, we have to know the exact form of initial distributions. It should be noted, however, that, to predict the behavior of the system on relatively short time spans, several moments of the distribution might be enough.
The question on the form of distributions can be treated mathematically (for a closely related model, see Ref
]. The fact is that a number of important distributions maintain their form over time. For example, if the initial distribution is gamma-distribution, it can be rigorously proved that it will always remain a gamma-distribution with parameters changing with time. The same statement is valid for normal distribution and others. However, there are distributions that change their form. For example, a uniform distribution becomes exponential
3. The model is strictly speaking ecological as all variants are assumed to be present at t0. My understanding is that cancer progresses through mutational events and these are not treated in this model.
Author response: We agree with this remark. Our model does not consider mutations. In our model, all variants are present at the initial time moment, perhaps, with extremely low densities.
4. Typically heterogeneity is spatial across a tissue.
Author response: The question on relation of spatial and parametric heterogeneity is important and interesting. However, in this work, we do not consider the spatial structure.
5. I still feel that a very important aspect of oncolytic therapy is the problem of target-degeneracy leading to infection of non-cancer cells. As these are not treated in this paper I remain suspicious of the efficacy of the complete clearance equilibrium. I think it would be interesting to explicitly include cell heterogeneity in order to treat both cell populations – healthy and cancerous.
: This issue was, actually, addressed in our response to Krakauer's comments on our previous paper (Ref
]). There are, indeed, many ways to expand these models. We stuck to a minimal version in order to keep the model within analytic solvability