Free boundary problems associated with cancer treatment by combination therapy

Many mathematical models of biological processes can be represented as free boundary problems for systems of PDEs. In the radially symmetric case, the free boundary is a function of \begin{document}$ r = R(t) $\end{document} , and one can generally prove the existence of global in-time solutions. However, the asymptotic behavior of the solution and, in particular, of \begin{document}$ R(t) $\end{document} , has not been explored except in very special cases. In the present paper we consider two such models which arise in cancer treatment by combination therapy with two drugs. We study the asymptotic behavior of the solution and its dependence on the dose levels of the two drugs.

Here the X i represent densities of different types of cells for i ≤ i ≤ m, the X i represent concentrations of proteins and other molecules for m + 1 ≤ i ≤ n; some of the molecules may be drugs used in the treatment of the cancer. The variable c is the concentration of nutrients, for example oxygen, which is consumed by the cells. The constants δ i , δ c and λ i are all positive. The functions F i (X, c) represent the interactions among cells and molecules, and are cancer specific. We denote by ∂Ω ∞ the lateral boundary of Ω ∞ , that is, the boundary of Ω ∞ which lies in 0 < t < ∞. This is a free boundary, moving with velocity u; the velocity u in Ω ∞ arises from the proliferation of cancer cells, which induces internal pressure within the cells population. This pressure has negligible effect on the relatively much smaller proteins and the other molecules.
Here we shall consider only the case of spherically symmetric variables, namely, functions of (r, t) where r = |x|, with u = u(r, t) x |x| and Ω(t) = {r < R(t)}.
Then the free boundary r = R(t) moves with the velocity We assume that the total density of the cells is constant, taking By adding the equations for cells, we obtain the following equation for u(r, t): The system (1)-(5) is complemented by initial and boundary conditions, consistent with Eq. (4). It is not difficult to prove that such a system has a unique solution in Ω ∞ with 0 < R(t) < ∞ for all t > 0; see, for instance, [7]. Some example of systems (1)- (5) with cancer treatment by combination therapy are given in [12,10].
In the present paper we are interested in the asymptotic behavior of solutions, and in the existence of steady state, for cancer treatment by combination therapy.
The asymptotic behavior of solutions of the system (1)- (5) and the existence of steady state were studied extensively in the case n = m = 1, also in the non-radially symmetric case; see [14,15] and the references therein. The uniqueness of a steady state and its stability was also established for n = m = 2 in the hyperbolic case where δ = 0 [5,4].
Wound healing is a process that can also be modeled by the system (1), (2) but the velocity u is defined differently, taking into account the viscoelastic properties of the partially healed region surrounding the wound. In this case, for the axially symmetric case in two dimensions (i.e., for "flat" wounds) it was proved that dR dt ≤ 0 for all t > 0, and, if the oxygen supply is restricted, then R(t) = cosntant > 0 for all t sufficiently large [8].
In the present paper we consider two different cancer models (1)-(5) with δ > 0 in the case n = m = 2, c ≡ constant. We include in the models treatment of the cancer with two drugs. We address the question how the two drugs affect the asymptotic behavior and/or the existence of steady states (i.e., benign tumors). In both models we begin with two PDEs: in Ω ∞ , and boundary conditions ∂P ∂n on the free boundary ∂Ω ∞ , for some β > 0. We assume that so that div u = f 1 + f 2 . The outward normal velocity V n of the free boundary is given by where n is the outward normal. We assume that P 0 ≥ 0, Q 0 ≥ 0, P 0 + Q 0 = 1, and that f 1 , f 2 are continuously differentiable, and are such that if P ≥ 0, Q ≥ 0 at t = 0 then P ≥ 0, Q ≥ 0 for all t > 0.
We focus on the radially symmetric case, where We can reformulate the above system as a free boundary problem for P with nonlocal coefficients, namely, ∂P ∂t with the boundary condition (8), and with free boundary r = R(t) satisfying the equation We note that and f 1 and f 2 are uniformly bounded for 0 ≤ P, Q ≤ 1. Hence for some ν > 0, so that the a priori estimates hold; hence the solution can be extended to all t > 0. The above system, with specific functions f 1 , f 2 also arises in a simple model of granulomas in tuberculosis. In this case some asymptotic properties of the solution were proved in [9], and the existence of steady states and their linear stability were established in [11].
The free boundary problem (13)-(15), (8) for P is just one example of free boundary problems for one PDE with non-local coefficients. There are several other such models, but with different free boundary conditions. In [1] the PDE has the form In [13], with free boundary this model arises in nonlocal branching of brownian motion [6]. Some free boundary problems were recast in [3] as non-local parabolic equations. Evolution of free boundary problems with the fraction of the Laplacian is yet another different type of a problem with non-local coefficients [2].
In the present paper we consider systems which arise in modeling cancer treatment by combination therapy. These systems are simplified versions of more comprehensive models studied in [12,10]. One of the two variables in the PDE system is the density of cancer cells, C. Another variable is the density of T cells. The immune system consists of several different classes of cells, and 'effective T cells' or, briefly, T cells form one of these classes. When T cells are active, they kill cancer cells, but in doing so, they may generate too much inflammation. For this reason they have built-in 'brakes' on their activity, called checkpoints. A checkpoint is a protein, called 'receptor', on the membrane of the T cells. Cancer cells produce complementary proteins called 'ligands'. When a ligand combines with a receptor of a T cell, it initiates a signaling cascade within the T cell which leads to blocking the activity of the T cell. In recent years several drugs have been developed which disable checkpoint receptors on T cells, and thus enable the T cells to remain active; they are called: checkpoint inhibitors. One of the two drugs in our models will be a checkpoint inhibitor.
We shall consider two models. In the first model, the second drug is a dose of oncolytic virus, a virus that is genetically programmed to invade only cancer cells, multiplies within them, and cause their death. We denote the density of virusinfected cancer cells by C i , and shall consider the system (13)- (15) and (8) with P = C, Q = C i , T = constant, and β = 0.
The oncolytic virus is injected at a rate γ V and the checkpoint inhibitor is injected at a rate γ A . As shown in [12], there is antagonism between the two drugs. If we increase γ V we make the T cells more active, so that they kill more effectively cancer cells. However, since they kill more effectively also the virus-infected cancer cells, this may decrease the effect of γ A . The result of increasing the drug γ V may then be an increase in the cancer volume. This was in fact demonstrated in [12], and it sends as important message for those that design clinical trials, to stay away from the regions of (γ A , γ V ) where antagonism exists.
The demonstration in [12] was done numerically. It would be very interesting to prove it rigorously. In this paper we do it for a much simpler model than in [12], but the model still captures the essence of the biology.
We introduce two parameters, C * and C * * , which depend on γ V and γ A and prove the following result: There is a region in the (γ V , γ A )-plane such that as γ A increases, C * − C * * changes sign from negative to positive. This means that increasing the amount of drug γ A has actually a negative result: it increases the tumor radius instead of decreasing it. This antagonism between γ A and γ V must be considered when designing clinical trials. In the second model, the second drug is cancer vaccine. It induces anti-cancer immune response and, in particular, it increases the introduction of T cells into the tumor. We then take P = C, Q = T , and boundary conditions with P 0 = 0, Q 0 = 1 and β > 0.
In this case γ A still represents the injection rate of checkpoint inhibitor, but γ V now represents the injection rate of cancer vaccine. Under some conditions we show that R(t) → 0 as t → ∞, and, under other conditions, that R(t) → ∞ as t → ∞. We also show under a different set of conditions that include assumptions on the initial data, that R(t) > R(0) for all t > 0. The question then arises whether there exist stationary solutions, i.e., begin tumors. We prove for some parameter regimes of (γ V , γ A ) that, for any small R, there exists a stationary solution with free boundary r = R and β = β(R).
In Section 2 we consider the model with oncolytic virus and checkpoint inhibitor. The model with cancer vaccine and checkpoint inhibitor is considered in Section 3. We conclude with some open problems in Section 4.
2. Model 1. Co-treatment with oncolytic virus. We denote by C and C i the densities of uninfected and virus-infected cancer cells, respectively, and consider the following system in the tumor Ω ∞ = {0 ≤ r < R(t), 0 < t < ∞}: where T is the density of T cells, and The parameter γ V represents the oncolytic virus therapy, which transforms cancer cells into infected cancer cells. The killing rate γ A of cancer cells consist of γ 0 + γ 1 where γ 0 is fixed, and γ 1 is the dose level of the checkpoint inhibitor. For simplicity, we shall refer to γ A as the checkpoint inhibitor drug. The parameter δ is positive and reflects the fact that infected cells are more easily detected and killed by T cells than uninfected cancer cells. The parameter K is the carrying capacity of the tissue: cancer cells' density cannot exceed K. We note that the boundary conditions in (20) mean that there is no flux of cancer cells across the tumor boundary. Let We then can state the following: Proof. Adding Eqs. (17), (18) and using Eq. (19), we get It will be convenient to rewrite Eq. (17) in the form where and Since has a unique zero in the interval (0, 1), given by C * , and

Consider the ODE system
From the above properties of f (C), we deduce that then, by comparing C(r, t) withĈ(t) forĈ(0) = c 1 and forĈ(0) = c 2 where 0 < c 1 < C(r, 0) < c 2 < 1, we see that From Eqs. (21)(25) we then get, We next observe that and the two zeros of h(C) are , then C * * ± are complex roots and, hence, h(C) < 0 for all C ≥ 0. On the other hand, if (λ + γ A δ) 2 > 4 λ K γ A (1 + δ), then C * * ± are two positive roots; the smaller positive root, C * * , may be smaller or larger than 1, but the larger positive root is definitely larger than 1 (since K ≥ 2). Fig. 1 shows the profiles of f (C) and h(C), and the proof of Theorem 2.1 then follows from (29).
We can interpret the results of Theorem 2.1 in terms of the drugs γ A and γ V . This is illustrated in Fig. 2 for the case K = 2, λ = 2, δ = 1. We see that if γ V is increased, we go from the region {R → ∞} to the region {R → 0}; thus by increasing the oncolytic drug γ V we get improved results for the anti-cancer treatment. However, by increasing γ A (the checkpoint inhibitor) we may go from the region {R → 0} to {R → ∞}; this means that increasing γ A does not decrease the cancer; on the contrary, it may increase the cancer. This conclusion was also established by simulations for the more comprehensive model in [10]; what it says is that "more is not always better," and in combination therapy one should be aware of the possibility that the interaction of the two agents may be antagonistic.
3. Model 2. Co-treatment with cancer vaccine. We consider the model with cancer cells (C) and T cells (T ), in Ω ∞ = {0 < r < R(t), t > 0}, and assume that We take boundary conditions ∂C ∂n and Here β > 0, and the boundary condition for T accounts for T cells migration from the lymph nodes into the tumor. The drug γ V is a cancer vaccine, which increases the population of T cells, while γ A accounts for the checkpoint inhibitor, as in the previous model. We deleted the carrying capacity K in Eq. (30) in order to somewhat simplify the subsequent computations. We assume that As before we have and we rewrite Eq. (30) in the form where and We are interested in the asymptotic behavior of the solutions of the system (30)-(35) and in the existence of stationary solutions. The results to be proved in what follows depend on the profiles of the functions h(C) and g(C). Figures 3(a), 3(b) show two scenarios of these profiles, depending on parameters C * , C * * in Fig. 3(a) and on parameters C * ± , C * * in Fig. 3(b); these parameters depend on doses γ A , γ V . We proceed to define these parameters.
We now show that under some conditions, R(t) → 0 as t → ∞.
We next show that R(t) remains uniformly positive under the following conditions,: Theorem 3.3. If the conditions (45)-(47) hold, then for all t > 0.
Proof. The proof is similar to the proof of Theorem 5 in [9]. We differentiate Eq. (30) in r to obtain an equation for G(r, t) ≡ ∂C(r,t) ∂r . Noting that G(r, 0) ≤ 0 and G(R(t), t) = −βC < 0, and applying the maximum principle to G, the assertion (48) follows.
We next integrate Eq. (30) over 0 ≤ r ≤ R(t), 0 ≤ t ≤ T , and use the boundary condition for C. We get Since ∂C ∂r ≤ 0, the last integral is larger than By assumption (46), the inequality (51) holds for small T . We claim that it holds for all T > 0. Indeed, otherwise, let T * be the smallest T such that We claim that R(t) > R(0) for all t ≤ T * . Indeed, otherwise let t 1 (t 1 ≤ T * ) be the smallest t for which R(t 1 ) = R(0). Then dR(t1) dt ≤ 0, so that and since t 1 ≤ T * we can use (52) for T = t 1 . We conclude that where R(t 1 ) = R(0), which is a contradiction to the assumption (47). Hence, R(t) > R(0) for all t ≤ T * . We then also have which is a contradiction to (53). We have thus proved that (51) holds for all T > 0, and hence (52) also holds for all T > 0, and R(t) > R(0) for all t > 0.

4.
Conclusion. The results of Sections 2,3 suggest several open problems. In the case of Section 2, one may extend the model by taking in equation (17) a variable T satisfying Eq. (31) (with its boundary condition in (33)), and C + C i + T ≡ 1. Then, after substituting T = 1 − C − C i we get, instead of Eq. (37), two equations, with C + C i < 1. The question whether the corresponding ODE system for C(t),Ĉ i (t) can be used to resolve the asymptotic behavior of R(t) remains to be addressed.
In the model of Section 3, we have shown, under very stringent conditions, that R(t) → 0, R(t) remains uniformly positive, or R(t) → ∞ as t → ∞. It would be interesting to establish these different behaviors under more general assumptions on the parameters of the system.