Abstract
Oncolytic viruses are genetically altered replicationcompetent viruses that infect, and reproduce in, cancer cells but do not harm normal cells. On lysis of the infected cells, the newly formed viruses burst out and infect other tumor cells. Experiments with injecting mutant herpes simplex virus 1 (hrR3) into glioma implanted in brains of rats show lack of efficacy in eradicating the cancer. This failure is attributed to interference by the immune system. Initial pretreatment with immunosuppressive agent cyclophosphamide reduces the percentage of immune cells. We introduce a mathematical model and use it to determine how different protocols of cyclophosphamide treatment and how increased burst size of the mutated virus will affect the growth of the cancer. One of our conclusions is that the diameter of the cancer will decrease from 4 mm to eventually 1 mm if the burst size of the virus is triple that which is currently available. The effect of repeated cyclophosphamide treatment is to maintain a low density of uninfected cells in the tumor, thus reducing the probability of migration of tumor cells to other locations in the brain. (Cancer Res 2006; 66(4): 23149)
 virotherapy
 glioma
 cyclophosphamide
 mathematical model
Introduction
Oncolytic viruses are genetically altered replicationcompetent viruses that can infect, and reproduce in, cancer cells but leave healthy normal cells unharmed. On lysis of an infected cell, a swarm of new viruses burst out of the dead cell and infect neighboring tumor cells. Tumor therapy by oncolytic viruses has been and continues to be actively tested in clinical trials for variety of malignancies ( 1– 8). Agents include ONYX015 ( 9, 10), the herpes virus G207 ( 11), the prostatespecific adenovirus CN706 and CN708 ( 12), and the reovirus that attacks tumors with activated Ras pathways ( 13). Additional oncolytic virus species are in development.
At present, despite producing many thousands of infectious virions per infected tumor cell, most virus species are unable to eradicate the majority of tumor modules. There is increasing evidence that the host response to an active viral infection plays a critical role in determining the overall efficacy of viral therapy. Indeed, it has been manifested that the innate immune system destroys infected cells as well as free virus particles, thus enabling the tumor to grow ( 14). Until very recently, the course of the immune response has been defined in the wildtype viral infection of normal organs in humans and animals but not in recombinant, oncolytic viral infections of tumors ( 15). A drug that seems effective in suppressing the innate immune response and thereby enhancing oncolysis is cyclophosphamide ( 16, 17).
In this article, we formulate and analyze a mathematical model of spherical glioma that has been injected at its center with oncolytic virus hrR3, which is a mutant of herpes simplex virus (HSV) 1. The model includes uninfected and infected tumor cells, necrotic cells, free virus particles, innate immune cells, and cyclophosphamide. Our model is similar to a model introduced by Wu et al. ( 18, 19) and Friedman and Tao ( 20); however, there are several important differences. First, in our model, we have included the effects of the immunosuppression drug cyclophosphamide. Indeed, one of the main aims of our article is to determine, by analysis of the mathematical model, the effect of administering cyclophosphamide under different protocols on the progress of glioma.
A second difference between our model and the model of Wu et al. ( 18) is that we include the presence of innate immune cells in the tumor tissue, whereas Wu et al. ( 18) includes only the immune response (tumor necrotic factor), which consists of molecules with negligible volume. This difference is important, because the immune cells make up to 50% of the total number of cells at some stages of the tumor progression.
Because the free virus particles are very small, they disperse in the fluid tissue like Brownian particles. We have therefore incorporated into our model a diffusion term for the free viruses. We have also added a term that accounts for the destruction of virus particles by the immune cells ( 21). Finally, in the mathematical model of Wu et al. ( 18), the parameters are estimated by using data from head and neck cancer ( 9, 22, 23). In our model, the parameters are estimated to conform with experimental results for glioma in article. ^{5} There is a substantial difference in the values of some of the parameters because glioma is a much more aggressive cancer.
Materials and Methods
The mathematical model describes a spherical tumor undergoing oncolytic virus therapy. The model parameters are based on experimental results of Fulci et al. ^{5} In these experiments, D74/HveC rat glioma cell lines were implanted into the brain of rats. After 7 days, the tumor reached the size of 4 mm in diameter, and the oncolytic virus hrR3 (which is a mutant of HSV) was injected into the center of the tumor. This mutant attacks tumor cells but does not attack healthy normal cells ( 6). Six hours after injection, some rats were sacrificed, the tumor was stained, the JPEG pictures were taken, and the infected area of the tumor was measured. This procedure was repeated after 72 to 76 hours from the time of virus injection and then again 1 week after virus injection. Five different stains were used: one for each of four innate immune cells, and one for the infected cells. The immune cells are CD11b (dendritic cells), CD68^{+} (main monocytes; also called ED1 in rats), CD163^{+} (also called ED2 in rats), and CD169 (also called ED3 in rats). To study the effects of the immunosuppression drug cyclophosphamide, some rats were preadministered by cyclophosphamide on the fifth day after tumor cell implantation (i.e., 2 days before virus injection).
The mathematical model includes infected tumor cells (y), uninfected tumor cells (x), necrotic cells (n), immune cells (z), and free virus particles (v). The quantity y represents the number density of infected cells (i.e., the number of infected cells in mm^{3}); the same meaning is associated with the quantities x, n, z, and v. Different immune cells may participate in the inflammatory response at different time points and may be cleared at different time points. Rapid upregulation was observed for certain immune cells, and depletion of some macrophages alone has also been produced. However, the cyclophosphamide effects on the different immune cells were only partially reproduced. For this reason, we shall take in our model the average response of all the immune cells rather than the response of each population of immune cells. The partial differential equations (PDE) of our model are listed in Appendix A. The parameters are derived as follows.
The uninfected cells proliferate at a rate λ, where λ is computed from Fulci et al., ^{5} for untreated gliomas. The infection rate of tumor cells by virus, β, is based on the observation ^{5} that ∼70% of the virus particles successfully invade uninfected cells. According to Fulci et al., ^{5} the mean lifetime of infected tumor cells is ∼18 hours. Hence, the infected tumor cell lysis rate, δ, is 1/18 h^{−1}. Necrotic cells are removed at the average rate of 2 to 3 days ( 24); we take the removal rate μ to be 1/48 h^{−1}. With regard to the diffusion coefficient of virus particles, D, according to Olmsted et al. ( 25) for HSV particles, D/D_{pbs} = 0.0089 ± 0.0052; according to Batchelor ( 26), D_{pbs} = 1.1 × 10^{−5} cm^{2}/s. We take D = 3.6 × 10^{−2} mm^{2}/h. This is nearly twice the diffusion coefficient in Chaplain et al. ( 27) as computed by the Einstein formula: D = (kT) / (6πRη), where η is the viscosity. However, this difference is to be expected because virus particles are not spherical; thus, their “effective” radius is smaller that the radius R used in the Einstein formula.
The burst size, b, of new viruses per one infected cell, is an important parameter. In wildtype HSV, it is in the thousands. However, for the oncolytic virus hrR3, it is much smaller, ranging from 10 to 100. Our default number is taken to be 50, but we shall also try to see how oncolytic virus therapy is effected if this number can be made larger.
The way immune cells are cleared is typically by membrane protein Fas and FasL, which, when combined on the cell surface, signal to caspase protein to execute the cell. The percentage of immune cells in the brain is typically 1% to 2%. When stimulated by infected cells in glioma, this percentage arises sharply. As the number of infected cells drops, the need for immune cells diminishes, so they undergo apoptosis either by killing themselves (using their own Fas and FasL to activate caspase) or by killing each other (Fas from one cell ligands to FasL from another cell; refs. 21, 28). The first process occurs when z is small and yields a linear clearance; the second process occurs when z is large and yields a quadratic clearance. Hence,
From ref. 28, one can infer the value of z_{0}, and combining this with the linear clearance rate in refs. 27, 29, we arrive at the number ω = c/z_{0} ≈ 20 × 10^{−8} mm^{3}/h cell. For simplicity, we shall assume only quadratic clearance with a rate ω as above.
It remains to estimate the following four parameters: a firstorder clearance rate of virus, γ; the killing rate of infected cells by the innate immune system, k; the immune takeup rate of virus, k_{0}; and the stimulation rate of innate immune cells by infected tumor cells, s. The experimental results ^{5} with and without cyclophosphamide treatments suggest that k = 2 × 10^{−8} mm^{3}/h immune cell. We assume that the immune cells take notice of, and are attracted to, the infected cells more than they are to the free virus. We therefore take k_{0} to be smaller than k: k_{0} = 10^{−8} mm^{3}/h immune cell.
Infected cells stimulate the innate immune system and cause a secondorder increase s · y · z of the number density of the immune cells. It is difficult to measure the stimulation rate s, but here again we derive an estimate from the experimental results, ^{5} which show that the percentage of immune cells increased from 1% before oncolytic virus injection to 30% to 50% after 6 hours. Accordingly, we take s = 5.6 × 10^{−7} mm^{3}/h immune cell. Finally, we take the firstorder clearance rate of virus to be γ = 2.5 × 10^{−2}/h. Table 1 summarizes the model parameters and their numerical values.
The parameters k, k_{0}, and s have been estimated rather crudely from Fulci et al., ^{5} and γ was determined by trying to fit the percentages of infected and immune cells derived by the mathematical model to the experimental percentages ^{5} 6 and 72 hours after injection of oncolytic viruses.
Results
Fitting of data. To compare our numerical simulation with the experimental results, ^{5} we take, as in Fulci et al., ^{5} the initial radius of the tumor to be 2 mm and the number of particleforming units of virus injected at the center to be 10^{8} to 10^{9}. The initial time is the seventh day after tumor implantation in the rat's brain.
Figure 1A compares the experimental measurements with our numerical simulation of the percentage of infected tumor cells (relative to the total number of all cells) without pretreatment of cyclophosphamide. Figure 1B compares the experimental measurement data with the numerical simulation of the percentage of the innate immune cells without pretreatment of cyclophosphamide.
After cyclophosphamide is given to the rats, the level of cyclophosphamide arises and reaches a plateau after 2 days. This level is maintained for ∼3 days and then begins to drop off to zero in the next 2 days. We simulate the cyclophosphamide level in the tumor by a piecewise linear functionand P(t) = 0 if t ≥ 120, where the unit of P(t) is 1/h.
Figure 1C compares the experimental measurements of immune cells, when cyclophosphamide was given, with our numerical simulation. The discrepancy between measurements and simulation develops only after a relatively long time both in Fig. 1B and C.
Figure 1 shows that our model fits quite well with experiments. We next proceed to compare the cancer progression with and without cyclophosphamide.
Comparison results. Figure 2 shows simulation results based on our model. Figure 2A shows the profiles of the averages over space of immune cell densities with and without cyclophosphamide pretreatment. Without cyclophosphamide, the immune cells reach the maximum 52% at the 26th hour after virus injection; with cyclophosphamide, the immune cells reach the maximum 34% at the 24th hour after virus injection. Thus, cyclophosphamide suppresses the maximum level of innate immune cells and shortens the time that the immune system reaches its peak. Because the effect of cyclophosphamide disappears after 120 hours, the percentage of the immune cells climbs up after 120 hours, thus forming a bimodal profile.
Figure 2B shows the profiles of the percentage of infected tumor cells with and without cyclophosphamide. Clearly, with cyclophosphamide, we expect to have more infected cells. As the simulation shows, without cyclophosphamide, the infected cells reach the maximum 46% at the 5th hour; with cyclophosphamide, the infected cells reach the maximum 50% at the 7th hour approximately. Figure 2C shows the profiles of the percentage of uninfected tumor cells with and without cyclophosphamide. The first thing to notice is that there is a time delay in the effect of cyclophosphamide; the immune suppression does not begin right away; in fact, there is a 3day delay. The effect of cyclophosphamide becomes negligible after ∼17 days. However, during the intermediate period, it is significant.
Discussion
We will use our model to make certain predictions on oncolytic therapy. We first examine what will happen within the glioma if we change the burst size without administering cyclophosphamide and then examine the effect of administering cyclophosphamide using different protocols.
Burst sizes. Suppose we inject into the tumor oncolytic virus that replicates at a faster rate than hrR3. A large burst size b of virus will increase the stimulation of the immune system, which will then attack the infected cells and the free viruses. As a result, the population of viruses will decrease and this will be followed by a decrease in the population of immune cells. The population of virus and infected cells will then be able to increase, and it will follow by a restimulation of the immune system, etc. Thus, we may expect a “feedback mechanism” that will cause an oscillatory behavior of the percentages of infected cells, immune cells, and uninfected cells, with slight time shift of the corresponding maxima. This is indeed shown in Fig. 3 for burst size b = 400 and 1,000. For smaller values of b, such as b = 200, oscillations occur only in the first 20 to 30 days. For b = 50 (data not shown), we do not see any oscillation.
Cyclophosphamide treatments. We shall compare two different protocols for administering cyclophosphamide. The first protocol is to administer a “normal” amount of cyclophosphamide weekly [as in P(t) above], and the second protocol is to administer twice the normal amount every 2 weeks. The simulation results of the percentage of uninfected tumor cells for burst sizes b = 100, 200, and 400 are shown in Fig. 4 .
These figures show that there is little difference between weekly and biweekly treatments, especially when we consider weekly averages.
It is instructive to compare these treatments with the “traditional” treatment, where we preadminister normal amount cyclophosphamide just once. The weekly treatment reduces the weekly averages of the uninfected cells, but due to oscillations (noted above), there are periods when this traditional treatment yields a small percentage of uninfected cells.
Tumor radius. All the preceding numerical results are based on solving the PDEs of the model given in Appendix A and then taking averages over the tumor region. We have obtained approximately the same numerical results using the much simpler ordinary differential equations obtained by neglecting spatial variations.
There is, however, one important quantity that we have not yet taken into account (i.e., the radius of the tumor), and to compute it, we need to work directly with the PDE system.
As the simulations presented above show, uninfected tumor cells will persist even with large burst sizes and with repeated cyclophosphamide treatments. However, if the radius of the tumor can be kept small enough, then longterm survival of the animals can be assured unless cell invasion and metastasis will occur as a result of cell shedding and migration.
Figure 5A simulates the growth of the tumor radius for b = 50 without cyclophosphamide and with one pretreatment of cyclophosphamide.
According to Fulci et al., ^{5} without cyclophosphamide, the rats die 8 to 10 days after the injection of viruses, and with one cyclophosphamide pretreatment, the rats die after 11 to 13 days after the injection of viruses; at death, the radius of the tumor is ∼6 mm. These experimental results roughly fit with the simulations in Fig. 5A.
Figure 5B simulates the growth of the tumor's radius without cyclophosphamide and with weekly cyclophosphamide treatment for different burst sizes.
We see that with weekly cyclophosphamide treatment the radius of the tumor will decrease as long as the burst size is bigger than 100; without any cyclophosphamide treatment and with burst size b already as high as 150, the radius of the tumor will decrease to 1 mm. If the burst size could be increased to 400, the tumor will shrink to an extremely small size. Finally, the weekly treatment by cyclophosphamide helps in decreasing the radius of the tumor as can be seen by comparing the profiles in Fig. 5B for burst size b = 130 and 150; for b = 130 with cyclophosphamide, we achieve a smaller radius than for b = 150 without cyclophosphamide.
Summary. We have shown that with the current burst size of oncolytic virus hrR3 the tumor cannot be eradicated with cyclophosphamide treatment; in fact, its radius will grow, and the rats will die within several weeks. If, however, the virus can be further altered to yield a burst size b ≥ 150, then the tumor will shrink to a very small size even with no cyclophosphamide treatment. It is well known that tumor cells in glioma may shed and migrate into other areas in the brain; this article does not consider this infiltration/invasion problem. Thus, even when the tumor size can be kept very small, there is still a chance of developing a secondary tumor. In this respect, a repeated treatment of the tumor by cyclophosphamide is important because it decreases the percentage of uninfected tumor cells and thus reduces the risk of secondary tumors.
As was shown by our model, there is little difference between weekly and biweekly cyclophosphamide treatments. The protocol of choice should therefore depend on the side effects to this chemotherapy.
Finally, our model considers only spherical gliomas with oncolytic virus injection at the center. However, we expect that the conclusions summarized above will hold also for nonspherical gliomas.
Appendix A
Consider a radially symmetrical tumor and denote by r the distance from a point to the center of the tumor. We denote the boundary of the tumor by r = R(t). Set

x(r,t) = number density of uninfected tumor cells,

y(r,t) = number density of infected tumor cells,

n(r,t) = number density of dead cells,

z(r,t) = number density of immune cells (the different types of immune cells are lumped together),

v(r,t) = number density of free virus particles, which are not contained in tumor cells,

P(t) = concentration of cyclophosphamide.
The proliferation and removal of cells cause a movement of cells within the tumor, with a convection term, for tumor cells x, in the form , where u(r,t) is the radial velocity; u(0,t) = 0. The other cells undergo the same convection. The HSV particles have a small diameter, typically 0.13 μm, compared with the typical diameter of 10 μm of a cell. Hence, free virus particles undergo diffusion rather than convection. By mass conservation law and convection of cells and by dispersion of virus particles, we deduce the following equations in r < R(t):
We assume that all the cells have the same size and that they are uniformly distributed in the tumor. Then, x(r,t) + y(r,t) + n(r,t) + z(r,t) = constant = ϑ, and by ref. 30, ϑ ≈ 10^{6} cells/mm^{3}.
Combining the above equations, we obtain
We also have ∂v/∂r(0,t) = 0, and ∂v/∂r[R(t),t] = 0 for t > 0 because the free viruses remain in the tumor. Finally, the free boundary is subject to the kinematic condition dR(t)/dt = u[R(t),t].
In our simulation, we have taken R(0) = 2 mm and the following initial values: x(r,0) = 0.84 × 10^{6}, y(r,0) = 0.1 × 10^{6}, z(r,0) = 0.06 × 10^{6}, v(r,0) = , 0 ≤ r ≤ 2, 5.2π × 10^{8} ≤ 4π∫_{0}^{2} rdr ≤ 11.2π × 10^{8}.
Acknowledgments
Grant support: National Science Foundation agreement 0112050.
The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.
Footnotes

↵5 G. Fulci et al. Experimental glioma virotherapy is enhanced by selective depletion of the innate immune response mediated by tumorinfiltrating CD163^{+} macrophages, submitted for publication, June 2005.
 Received July 27, 2005.
 Revision received October 7, 2005.
 Accepted December 16, 2005.
 ©2006 American Association for Cancer Research.