The role of TNF-α inhibitor in glioma virotherapy: A mathematical model.

Virotherapy, using herpes simplex virus, represents a promising therapy of glioma. But the innate immune response, which includes TNF-α produced by macrophages, reduces the effectiveness of the treatment. Hence treatment with TNF-α inhibitor may increase the effectiveness of the virotherapy. In the present paper we develop a mathematical model that includes continuous infusion of the virus in combination with TNF-α inhibitor. We study the efficacy of the treatment under different combinations of the two drugs for different scenarios of the burst size of newly formed virus emerging from dying infected cancer cells. The model may serve as a first step toward developing an optimal strategy for the treatment of glioma by the combination of TNF-α inhibitor and oncolytic virus injection.


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-α. It was demonstrated there that inhibition of TNF-α 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-α 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-α 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-α 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-bloodbarriers (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. CED has recently reached phase I/II 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-α inhibitor, have negative sideeffects. 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-α 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-α 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-α. 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-α 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 ≤ b ≤ 150. Using unit of g cm 3 , the parameter b is the burst size defined as b × mass of virus mass of cells = b × 10 −6 . The dynamics of the model is expressed mathematically by the following system of ODEs: All the densities and concentrations are in unit of g cm 3 . In Eq. (1) α represents the proliferation rate of uninfected cancer cells and δ x is the death rate; β is the infection rate of x by viruses v. In Eq. (2) the term ky T K+T represents the necrotic  4), the first term on the right-hand side is the production of TNF-α 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. ρxv) and natural degradation/clearance (δ 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-α inhibitor (in Eq. (4 ).
In the present paper, we take u 1 ≡ const = C and u 2 ≡ 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 g cm 3 , ρ β = virus mass cell mass is a small number, and similarly κ 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 where µ 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 dm dt . Adding up the equations (1)- (3) and (6) we get that The total mass of the tumor then increases at rate where z denotes the average of z. Let θ 0 is be the sum of averages, θ 0 = x+ y+ M + n.
We assume that an increase in total mass causes the tumor volume to increase proportionally, that is by θ 0 dV dt and that yM = y M , so that Finally, in order to simplify the calculations, we assume that the averages x, y, M and v satisfy the same equations that x, y, M , v. Hence we get the following formula for the tumor radius: 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: 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) = 43cm. For b = 90, R(50) is still large, namely R(50) = 6cm. 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-α, 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.

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) ≡ C = 5 · 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 ≤ b ≤ 110 in Fig. 4 and 110 ≤ b ≤ 150 in Fig. 5) for different values of C (above and below C = 5 · 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 · 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-α 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 ∈ [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) ≡ C and u 2 (t) ≡ D. For our simulations we will take C = 5 · 10 −7 and D = 15. Again we will study the response of the system for various levels of viral production by the infected cells.  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 · 10 −7 in comparison to 7 · 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.1cm. Now it becomes natural to look closer at the therapy itself and see how the two main agents, viral infusion and TNF-α 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 · 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-α 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 · 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 ≥ 123, with the low dose of D, D = 3 (i.e., the production of TNF-α  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) ≡ C and u 2 (t) ≡ D. The efficacy of the combined  therapy is defined by 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-α 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 ∈ [10 −7 , 5 · 10 −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 macrophagesproduced TNF-α can enhance virotherapy treatment in glioma. In the present paper we developed a mathematical model which includes both drugs, oncolytic viruses and TNF-α 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-α 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.
We assume that the infected cells die faster than x cells and take δ y = 0.2/day. The diameter of the virus is 0.13µm. The diameter of the cell is 10µm. Hence the mass of 5 · 10 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 · 10 −15 g. From [5,8] we get that 7 · 10 −10 < β < 1.7 · 10 −8 in mm 3 virus · day .
The value of K = 5 · 10 −7 g cm 3 is taken from [17] and the value of δ 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 δ y , which gives k = 0.4/day. The ratio κ k is the same as the ratio of mass of T mass of cell , which we take to be 10 −9 . Hence κ = 10 −9 k = 0.4 · 10 −9 /day = 4 · 10 −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 δ M = 0.015/day. The source of activated macrophages in the healthy brain is constitutively small; we take A = 9 · 10 −7 g cm 3 ·day . The parameter s is difficult to measure; we take s = 0.15 cm 3 g·day . According to [9] macrophages infected by M. tuberculosis produce TNF-α at the rate λ = 1.07 · 10 −3 /day. We assume that the macrophages activated by the infected tumor cells are capable of producing TNF-α at a larger rate, taking λ = 2.86 · 10 −3 /day.