1. Introduction

Oncolytic viruses are genetically altered replication-competent viruses which infect and reproduce in cancer cells but do not harm normal cells. When an infected cell dies many newly formed viruses are released and spread out, infecting neighbouring tumor cells. Treatment by oncolytic viruses (OV) has been and continues to be actively tested in clinical trials for various types of cancer with the use of variety of viruses, [4], [10], [12], [13], [14], including ONYX-15 [7], [11], herpes simplex virus (HSV) [16] and prostate-specific adenovirus CN706 and CN708 [18].

This therapy, although based on quite promising assumptions, encounters one major obstacle: the innate immune system recognizes the infected cells and destroys them before the viruses within them get a chance to multiply [3].

It was reported in [6] that CD 163+ macrophages, in rat experiments for glioma, inhibited OV therapy making it unsuccessful. The solution suggested in [6] was to use cyclophosphamide (CPA) as a suppressant of the immune response through the inhibition of CD 163+ and thus enhance the effectiveness of the OV therapy. This approach has been studied from the mathematical point of view by Friedman et al. in [5]. The authors formulated a mathematical model of virotherapy following the earlier work by Wu et al.[21,22], but focusing on the data from glioma rather than from head and neck cancer. The model in [5] was described by a system of PDEs and the effects of the therapy with and without CPA was analysed.

In the present paper we intend to pick up on this work, but pursue a different avenue based on a very recent paper by Auffinger et al. [2]. In that paper it was suggested that in order to enable the effective action of the virotherapy one should try to block the main "weapon" used by macrophages, namely the TNF-$\alpha$. It was demonstrated there that inhibition of TNF-$\alpha$ could significantly enhance virus replication and the efficacy of the overall treatment.

Thus our goal here is to construct a model which captures the interactions between healthy tumor cells, infected tumor cells, viruses, and macrophages and the TNF-$\alpha$ they produce. The model is based on the work of Friedman et al. [5] and the PDE model presented there. However, here we will formulate a reduction of this model from the spatial PDE model to a population type ODE model. This reduction will enable us to pursue a detailed dynamical system analysis of the model as well as analysis of it as an optimal control problem. Although the spatial element has to be compromised for that, it is not entirely removed from the analysis. Indeed, we will be able to estimate the tumor radius in terms of the cells population. The efficacy of both treatments by injection of virus and TNF-$\alpha$ inhibitor will be analysed in the context of the radius of the tumor bringing the spatial aspect back into the picture.

The approach pursued in our paper will be to target the tumor by combining the two therapies: the viral injection and the TNF-$\alpha$ inhibitor. We will analyze the response of the system to various doses, particularly the efficacy of the therapy, having as a goal the minimization of the tumor radius. The administration of the virus will be pursued through a continuous injection of a fixed dose. This is motivated by the fact that the therapeutic efficacy of systemic chemotherapeutic agents is significantly limited due to the presence of highly selective brain-blood-barriers (BBB). In order to overcome this limitation, researchers have developed a new delivery method, Convention Enhanced Delivery (CED), which relies on intracranial delivery of viral vectors through continuous infusion via catheters. CEDhas recently reached phase Ⅰ/Ⅱ in clinical trials [2], although so far this procedure has been shown to be only moderately effective in recurrent glioma patients [2].

Both therapeutic agents, the virus and the TNF-$\alpha$ inhibitor, have negative side-effects. Side-effects of virus such as HSV used in tumor virotherapy are significant in rodent tumor models and it is not known if these will occur also in humans [20]. In the case of TNF-$\alpha$ inhibitors the side effects include opportunistic infections like tuberculosis, listerosis and pneumocystis. There are also adverse events including systemic lupus, congestive heart failure and general infections [1]. Because of these side effects, the amount of dosing should be optimized to maximize the effect and minimize the negative side effects. In this paper we will pursue study of the efficacy which is a natural first step towards analysis of this model as an optimal control problem.

---

2. A mathematical model for the treatment of glioma with virotherapy and TNF-$\alpha$ inhibitors

Let $x(t)$ denote the density of cancer cells (uninfected), $y(t)$ the density of infected cancer cells, $v(t)$ the density of the virus, $M(t)$ the density of the macrophages, and $T(t)$ the concentration of TNF-$\alpha$. The model includes two controls, $u_{1}$ and $u_{2}$. The control $u_{1}(t)$ represents the amount of the virus that is injected into the tumor, and the control $u_{2}(t)$ stands for the dosage of the TNF-$\alpha$ inhibitor. The "in" and "out" flows between the compartments and the actions of the two drugs are schematically illustrated in Figure 1.

The burst number is the number of virus that emerge from a dying cancer cell. We shall take it in the range of $50 \leq b \leq 150$. Using unit of $\frac{g}{cm^3}$, the parameter $b$ is the burst size defined as $b \times \frac{\text{mass of virus}}{\text{mass of cells}}= b\times 10^{-6}$.

The dynamics of the model is expressed mathematically by the following system of ODEs:

$\frac{dx}{dt} =\alpha x-\beta xv-\delta_{x}x,$

(1)

$\frac{dy}{dt} =\beta xv-ky\frac{T}{K+T}-\delta_{y}y,$

(2)

$\frac{dM}{dt}=A+syM-\delta_{M}M,$

(3)

$\frac{dT}{dt} =\frac{\lambda}{1+u_{2}}M-\kappa y\frac{T}{K+T}-\delta_{T}T,$

(4)

$\frac{dv}{dt} =b_{1}ky\frac{T}{K+T}+b\delta_{y}y-\rho xv-\delta_{v}v+u_{1},$

(5)

All the densities and concentrations are in unit of $\frac{g}{cm^3}$. In Eq. (1) $\alpha$ represents the proliferation rate of uninfected cancer cells and $\delta_x$ is the death rate; $\beta$ is the infection rate of $x$ by viruses $v$. In Eq. (2) the term $ky\frac{T}{K+T}$ represents the necrotic death of infected cells caused by TNF-$\alpha$, while $\delta_y$ is the death by apoptosis. When a cell $y$ dies by apoptosis, $b$ virus particles are released, while if it dies by necrosis a very small number, $b_1$, of viruses emerge. These are accounted by the first two terms on the right hand side of Eq. (5). In Eq. (3) the terms $A$ and $\delta_M M$ represent the source and death of macrophages under healthy normal conditions, while $syM$ accounts for the tumorigenic response of the immune system invoked by the infected cells $y$. In Eq. (4), the first term on the right-hand side is the production of TNF-$\alpha$ by macrophages, while the remaining two terms are loss by absorption within y cells and loss by natural degradation. The virus equation (5) includes virus particles from dead $y$ cells and loss from absorption by $x$ (i.e. $\rho xv$) and natural degradation/clearance $(\delta_v v)$. We also included in the model a continuous injection $u_1$ of virus, as virotherapeutic drug, and a continuous injection $u_2$ as TNF-$\alpha$ inhibitor (in Eq. (4).

In the present paper, we take $u_1 \equiv const = C$ and $u_2 \equiv const = D$, but in future work we shall use $u_1$ and $u_2$ as time-dependent adaptive control variables in order to determine optimal strategies for glioma treatment. Since all the density and concentration in (1)-(5) are in the same unit of $\frac{g}{cm^3}$, $\frac{\rho}{\beta} = \frac{virus\;mass}{cell\;mass}$ is a small number, and similarly $\frac{\kappa}{k}$ is a small number. We expect $b_1$ to be very small (viruses are damaged during necrosis) so for simplicity we shall take $b_1 = 0$. As in [5] the viruses $burst$ (or $replication$) $number$ $b$ will play a major role in the progression of the disease and its treatment.

We denote by $n(t)$ the density of all the dead cells. Then, in addition to the dynamics given by Eqs. (1)-(5), we have the equation

$\frac{dn}{dt}=ky\frac{T}{K+T}+\delta_{y}y+\delta_{x}x+\delta_{M}M-\mu n \label{Dyn6},$

(6)

where $\mu$ is the rate by which dead cells are cleared out of the tumor.

Although we do not intend to directly focus on this equation, it will be taken into account in the calculation of the tumor radius.

---

2.1. Calculation of the tumor radius

We assume that the tumor is spherical with variable radius $R(t)$, volume $V(t)$ and mass $m(t)=x(t)+y(t)+M(t)+n(t)$. At each point of the sphere, the density of the cells increases at the rate $\frac{dm}{dt}$. Adding up the equations (1)-(3) and (6) we get that

$\frac{dm}{dt}=\alpha x-\delta_{x}x+A+syM-\mu n.$

The total mass of the tumor then increases at rate

$V(t)(\alpha\widetilde{x}-\delta_{x}\widetilde{x}+A+s\widetilde{yM} -\mu\widetilde{n})$

where $\widetilde{z}$ denotes the average of $z$. Let $\theta_{0}$ is be the sum of averages, $\theta_{0}=\widetilde{x}+\widetilde{y}+\widetilde{M}+\widetilde{n}$. We assume that an increase in total mass causes the tumor volume to increase proportionally, that is by $\theta_{0}\frac{dV}{dt}$ and that $\widetilde{yM}=\widetilde{y}\widetilde{M}$, so that

$\theta_{0}\frac{dV}{dt}=V(t)(\alpha\widetilde{x}-\delta_{x}\widetilde {x}+A+s\widetilde{y}\widetilde{M}-\mu\widetilde{n}).$

Since

$\frac{1}{V}\frac{dV}{dt}=\frac{3}{R}\frac{dR}{dt},$

we get

$\theta_{0}\frac{3}{R}\frac{dR}{dt}=(\alpha\widetilde{x}-\delta_{x} \widetilde{x}+A+s\widetilde{y}\widetilde{M})-\mu\widetilde{n}.$

Hence

$\theta_{0}\frac{3}{R}\frac{dR}{dt}=(\alpha\widetilde{x}-\delta_{x} \widetilde{x}+A+s\widetilde{y}\widetilde{M})-\mu(\theta_{0}-\widetilde {x}-\widetilde{y}-\widetilde{M}).$

Finally, in order to simplify the calculations, we assume that the averages $\widetilde{x}$, $\widetilde{y}$, $\widetilde{M}$ and $\widetilde{v}$ satisfy the same equations that $x$, $y$, $M$, $v$. Hence we get the following formula for the tumor radius:

$\theta_{0}\frac{3}{R}\frac{dR}{dt}=(\alpha x-\delta_{x}x+A+syM)-\mu(\theta_{0}-x-y-M). \label{Dyn7}$

(7)

---

2.2. Determination of the parameters

Table 1 gives the values of the parameters which will be used in our analysis. Explanations concerning the calculations of the parameter values are given in the Appendix.

---

3. Results

In this section we simulate the model (1)-(7) in order to determine how the state of the system responds to a combined therapy. Values of the parameters for the simulations will be taken from Table 1, unless specified otherwise. We assume that the process starts with the following initial conditions:

$x(0)=0.7,\;\;y(0)=0.1,\;\;n(0)=0,\;\;M(0)=0.1,\;\;T(0)=10^{-7}, \frac{4\pi}{3}R(0)^3=0.9,$

so that approximately $R(0) = 0.6$.

---

3.1. Behavior of the model without therapy (uncontrolled system)

We assume that initially a dose of viruses given by $v(0) = 10^{-6}$ injected and we start by studying the behavior of the system when no additional therapy is given.

Figure 2 shows the profile of all the variables of the model for the first 50 days for different values of the burst number $b$. If $b$ is small, namely $b=70$, the tumor radius increases from its initial value $R(0)=0.6$ to $43$ at $T=50$, that is, $R(50) = 43 cm$. For $b=90$, $R(50)$ is still large, namely $R(50)=6 cm$. It is only when $b=150$ that $R(50)$ does no longer increase relative to $R(0)$. These results are in agreement with the mouse experiments in [5] in the sense that if $b=50$ the tumor radius quickly increases while if $b=150$ the tumor radius begins to decrease. As in [5], we see also here the oscillations that occur in each of the variables. These oscillations are natural. For example, as $y$ oscillates, so will the viruses which are released from dying infected cells, and so will the rate of macrophages whose activation depends on $y$.

The simulation of $T(t)$ with $b=70$ suggests that the initial load $v(0)$ results in a massive increase of TNF-$\alpha$, which later seems to normalize. We also see that after the initial injection of the virus, the virus is multiplying to achieve its maximum by day $28$; after that time the virus density keeps decreasing -the drug is 'too weak'. As a result, the growth of $R$ after day $30$ is limited. However, the increasing amount of virus in the first $30$ days causes $M$ to grow -which leads to growth of $T$ mentioned above. And, as earlier, large amounts of $T$ restrain the replication of virus residing in infected cells $y$ -which is bad for the therapy.

Moving on to the case with $b = 90$, shown in Figure 2, we observe that the behaviour of $T$ is almost the same as in the previous case. The radius $R$ is still increasing in time, but not as sharply as before. Finally, the case of $b = 150$ presented in Figure 2 shows that the growth of $v$ is faster, but similarly it decays later to finally restart growing again. We can also see that $R(t)$ is first slightly growing, then it decreases, and finally it slightly increases again, mimicking the oscillation of $v$.

The observed behavior shows that if we would want to achieve the tumor size reduction without continuous therapy, but just with initial viral injection at time $t = 0$, we would need to inject a very high (probably unrealistic) amount of the virus. This suggests that continuous therapy of virus injection is needed.

---

3.2. Behavior of the model with continuous constant viral infusion

We will now study the behavior of the system assuming that we apply constant viral injection $u_1(t) \equiv C = 5\cdot 10^{-7}$ and the initial load $v(0) = C$, again for three viral replication numbers, $b = 70$, $b = 90$, $b = 150$. Figure 3 shows the evolution of the system under these assumptions. We see the same oscillating behaviour as before. But, more importantly, we can see the effect of the drug on the tumor radius: for $b=70$ and $b=90$ the radius $R=R(t)$ is still increasing, although much less than in the case of $C=0$ (no continuous infusion). But for the burst number $b=150$ the tumor radius, after a small initial increase, is strictly decreasing, with $R(50)$ approximately half the initial tumor radius.

In case $b=70$ the radius $R(t$) is increasing, but the size after $50$ days is almost $5$ times smaller then in Figure 2, though we gave smaller initial dose than in the case without therapy. Again $v(t)$ reaches its maximum value (now at day $25$), but it does not die out later on. Furthermore, the growth of $M$, and consequently of $T$, is smaller, and it has less impact on killing $y(t)$. Heading to the case of $b = 90$, we are still not in the most desirable situation -the radius is still growing, although $3$ times less than in Fig. 2 at day $50$. As might have been expected, for $b=150$ the tumor radius is decreasing, but only after a brief increase. This initial increase is probably caused by the fact that the initial virus load is small. The graphs of $v(t)$ suggest that the bigger $b$ is, the earlier the maximum of $v$ is reached, and the larger is its value.

Figures 4 and 5 show the profile of $R(50)$ as a function of $b$ (for $70 \leq b \leq 110$ in Fig. 4 and $110 \leq b \leq 150$ in Fig. 5) for different values of $C$ (above and below $C= 5 \cdot 10^{-7}$). We see that $R(50)$ is a monotone decreasing function of $b$; furthermore, for smaller $b$'s the decline in $R(50)$, as $b$ increases, is more steep. For each value of $b$, we can determine the exact value of $C$ for which the drug will decrease $R(50)$ below its initial size. For example, if $C = 5 \cdot 10^{-7}$, then $R(50)$ will be smaller than $R(0)$ only if $b>125$.

Figure 5 suggests that, for $b<120$, the doses around $10^{-7}$ do not decrease $R$ by day 50 -but for $b>120$, the smaller doses will make $R(50)$ smaller than $R(0)$ (without using the TNF-$\alpha$ inhibitor $D$).

Additionally, in that case we do not have to use $D$ to cure the patient. Similar behavior is shown in Figure 4, for $b\in [70,110]$. We can also conclude that the greater $C$ is, the less the effect of additional units of viruses, especially in the second case. Hence, administrating more $C$ and counting that it will solve the problem is not a good idea. $R(50)$ as a function of $b$, for every fixed $C$, is of course decreasing, and also it is a convex function, in both cases. However, our goal, which is the tumor size reduction, cannot be fulfilled by applying only constant viral injection, unless we assume very high, probably unrealistic rate of production of the virus.

---

3.3. Behavior of the model with combined therapy

Next we will assume that the combined therapy is given at the constant rates, i.e., $u_1(t) \equiv C$ and $u_2(t) \equiv D$. For our simulations we will take $C = 5 ·\cdot 10^{-7}$ and $D = 15$. Again we will study the response of the system for various levels of viral production by the infected cells. The responses for $b = 70$, $b = 90$, $b = 150$, are shown in Figure 6, which is analogous to Fig. 3.

Firstly, for $b = 70$, the tumor size increases, but not as much as in the case $D=0$. The graph of $T(t)$ looks similar to the case without $D$ therapy, the values of $T(t)$ are of an order of magnitude lower than earlier -now the maximum is close to $3.5 \cdot 10^{-7}$ in comparison to $7 \cdot 10^{-6}$ where $D=0$. Applying $u_2$ also supports growth of $v$ -now it is about $5$ times higher. What is interesting, is that although $T$ is decreased, $M$ is larger than before.

For $b=90$ the tumor radius, after initial increase, goes down to reach its original size. Finally for $b=150$, we obtain our desired goal for tumor reduction; the radius decreases significantly to about $R(50)= 0.1 cm$.

Now it becomes natural to look closer at the therapy itself and see how the two main agents, viral infusion and TNF-$\alpha$ inhibitor, contribute to the success of the therapy. In other words, what is the contribution of each of them into the goal of reducing the tumor radius. Figure 7 shows how the system evolves in terms of $R(50)$, for variable replication number $b$, while applying different amounts $D$. We take $C = 5 \cdot 10^{-7}$ and study the range of $b$ of $[70,90]$. We see that even with $b=80$ we can decrease $R(50)$ from $R(0)$ if we take $D = 30$, that is, if we restrict the production of TNF-$\alpha$ by 97$\%$.

In Figure 8, we take the range of $b$ to be in the interval $[110,130]$ and use a lower $C$, namely $C = 3 \cdot 10^{-7}$. Even with such high burst number, with $D=1$, $R(50)>R(0)$. This case shows the significant effect of small doses of D may have, namely, for $b\geq 123$, with the low dose of $D$, $D=3$ (i.e., the production of TNF-$\alpha$ was restricted by only 75%) we get the desired shrinkage of the tumor radius i.e. $R(50)<R(0)$.

Summarizing, our goal of shrinking the tumor size can be achieved under different scenarios, so in choosing the best combination, one should take into consideration potential serious side-effects of the treatment. This brings us to the main topic of this work, which is the evaluation of the efficacy of both treatments.

---

3.4. Efficacy of the combined therapy

In order to capture more clearly the benefits of the dual therapy by $C$ and $D$, we introduce the concept of efficacy. Let $T$ denote the duration of the therapy. For our analysis we will choose the window of $T=50$. We denote by $R(C,D)$ the radius of the tumor at the day $50$ under the combined treatment with $u_{1}(t)\equiv C$ and $u_{2}(t)\equiv D$. The efficacy of the combined therapy is defined by

$E(C,D)=\frac{R(0,0)-R(C,D)}{R(0,0)}$

Figure 9 is an efficacy map for the case $b=90$, showing the efficacy of the combined treatment with $C$ varying along the horizontal axis and $D$ along the vertical one. We can see that the efficacy of the combined treatment increases with either $C$ or $D$. For small $C$, the efficacy goes up sharply with D. For $D>4$ (TNF-$\alpha$ production is inhibited by more that 80%) the efficacy increases slowly with $C$. An efficacy map for $b=70$ has the same features as in Fig. 9 (not shown here). Thus, for small doses of $D$ the efficacy goes sharply up with $C$, but for larger doses of $D$ the efficacy increases slowly with $C$. This confirms the observation made in section 3.3.

Figure 10 is the efficacy map for the case $b=150$. Here we take a smaller range of C, $C \in [10^{-7},5\cdot10^{-7} ]$. The reason is, that for higher values of $C$ so the effect of $D$ is not that visible; $R(50)$ will become very small even if $D=0$.

---

4. Discussion and conclusion

Virotherapy represents a promising treatment of glioma. However the innate immune response, in particular by macrophages, reduces significantly the effectiveness of the treatment, by killing virus infected cancer cells prematurely. Recent experiments [15] show that blocking macrophages-produced TNF-$\alpha$ can enhance virotherapy treatment in glioma. In the present paper we developed a mathematical model which includes both drugs, oncolytic viruses and TNF-$\alpha$ inhibitor. We observed that the burst number b, that is, the number of new viruses released from the dead cancer cell plays a critical role in the model. Without any treatment, if $b$ is too 'small' the tumor radius will continue to grow, but for larger $b$, say $b>150$ the tumor radius will decrease. More generally, given any combined therapy $(C,D)$ and the terminal treatment time $T$, the tumor radius $R(t)$ at time $t=T$ will be smaller then the initial radius $R(0)$ if and only if $b$ exceeds certain threshold number, $b(C,D)$. Furthermore, we developed an efficacy map of treatment that depends on $b$. We illustrated it for small $b=90$ and large $b=150$. For $b=90$ the efficacy increases more sharply with the increase in $C$, whereas for $b=150$ the efficacy is almost independent on $C$. This means that for smaller burst number the increase in $C$ has a much larger effect on reducing the tumor radius than if the burst number is larger. Generally the maps show the importance of the TNF-$\alpha$ inhibitor application in the combined treatment. Even a small dose applied we can obtain better efficacy than raising the level of viral infusion very significantly. However, side effects of both drugs in combined treatment will require to limit the amount of doses. Under such limitations the efficacy map could be used as a prognosis tool. In this paper, we use constant amounts of each of the drugs but this does not have to be the case, because the patient situation may require to modify these levels as the treatment goes on. Future work should include developing schemes for evolving treatment of glioma with $u_1=u_1(t)=C(t)$ and $u_2=u_2(t)=D(t)$ as functions of time keeping as a goal minimizing the tumor and the side effects. This will require formulating the model as an optimal control problem. We plan to pursue analysis using methodologies as in [19] of this model with the goal to compute optimal time varying regimens for this problem.

---

5. Appendix: Parameter estimation

The model in [5] was based on mice experiments. Prorating some of the parameters to humans, we take

$\alpha=0.2/day,\;\;\;\;\delta_{x}=0.1/day,\;\;\;\; \mu=0.2/day.$

We assume that the infected cells die faster than $x$ cells and take $\delta_{y}=0.2/day.$ The diameter of the virus is $0.13\mu m$. The diameter of the cell is $10\mu m$. Hence the mass of $5\cdot10^{5}$ viruses is approximately equal to the mass of one cell. The mass of the cell is $10^{-9}g$. Hence the mass of one virus is $2\cdot10^{-15}g.$

From [5,8] we get that

$7\cdot10^{-10}<\beta<1.7\cdot10^{-8}\;\;\;\;\text{in}\frac{mm^{3} }{virus\cdot day}.$

We take

$\beta=10^{-9}\frac{mm^{3}}{virus\cdot day}$

which in units of $\frac{cm^{3}}{g\cdot day}$ gives

$\beta=10^{-8}\frac{10^{-3}}{2\cdot10^{-15}}\frac{cm^{3}}{g\cdot day} =2\cdot10^{4}\frac{cm^{3}}{g\cdot day}.$

The parameter $\rho$ is in units of $\frac{cm^{3}}{cell\cdot day}$. Hence

$\rho=\frac{\beta}{5\cdot10^{5}}=\frac{2\cdot10^{4}}{5\cdot10^{5}} =4\cdot10^{-2}\frac{cm^{3}}{g\cdot day}.$

The value of $K=5\cdot10^{-7}\frac{g}{cm^{3}}$ is taken from [17] and the value of $\delta_{T}=55.45/day$ is taken from [9]. We assume that if $T=K$ then the death rate of $y$ by $T$ is twice the death rate $\delta_{y}$, which gives $k=0.4/day.$ The ratio $\frac{\kappa}{k}$ is the same as the ratio of $\frac{mass\text{ of }T}{mass\text{ of cell}}$, which we take to be $10^{-9}.$ Hence

$\kappa=10^{-9}k=0.4\cdot10^{-9}/day=4\cdot10^{-10}/day$

In [9] the death rate of macrophages depends of their phenotype: for $M_{1}$ phenotype the death rate is $0.02/day$ and for $M_{2}$ phenotype, it is $0.008/day.$ Accordingly we take $\delta_{M}=0.015/day$. The source of activated macrophages in the healthy brain is constitutively small; we take $A=9\cdot10^{-7}\frac{g}{cm^{3}\cdot day}.$ The parameter $s$ is difficult to measure; we take $s=0.15\frac{cm^{3}}{g\cdot day}.$

According to [9] macrophages infected by M. tuberculosis produce TNF-$\alpha$ at the rate $\lambda=1.07\cdot10^{-3}/day$. We assume that the macrophages activated by the infected tumor cells are capable of producing TNF-$\alpha$ at a larger rate, taking $\lambda=2.86\cdot10^{-3}/day$.

---

Acknowledgments

The research of U. Ledzewicz is partially supported by the National Science Foundation under research Grant No. DMS 1311729. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation. A. Friedman has been partially supported by the Mathematical Biosciences Institute (MBI) of The Ohio State University.

---