Dynamical properties and tumor clearance conditions for a nine-dimensional model of bladder cancer immunotherapy.

Understanding the global interaction dynamics between tumor and the immune system plays a key role in the advancement of cancer therapy. Bunimovich-Mendrazitsky et al. (2015) developed a mathematical model for the study of the immune system response to combined therapy for bladder cancer with Bacillus Calmette-Guérin (BCG) and interleukin-2 (IL-2) . We utilized a mathematical approach for bladder cancer treatment model for derivation of ultimate upper and lower bounds and proving dissipativity property in the sense of Levinson. Furthermore, tumor clearance conditions for BCG treatment of bladder cancer are presented. Our method is based on localization of compact invariant sets and may be exploited for a prediction of the cells populations dynamics involved into the model.


1.
Introduction. Most mathematical models for population dynamics are defined by ordinary differential equations (ODEs). For systems with large number of equations, it is highly important to determine their dynamics in detail, although it is frequently extremely difficult or impossible to obtain accurate results. For example, rigorous studies of global interactions dynamics between tumors and the immune system were performed by Starkov and coauthors in [15] for the Kirschner-Panetta model [6] and in [16] for the Owen-Sherratt model [11]. In-depth analysis of both local and global dynamics in Bunimovich-Mendrazitsky's model [1] of bladder cancer were recently accomplished in [4,17].
The current work is dedicated to dynamical analysis of the tumor cell population model elaborated by Bunimovich-Mendrazitsky et al. in [3]. This model uses nine non-linear ODEs investigated numerically [3], in order to describe the effects of combined BCG and IL-2 immunotherapy for bladder cancer cells and upon other cells populations in order to evaluate comparative therapeutic scenarios. For the class of models defined by ODEs we propose to examine the following properties as indicators of the model's validity: 1) positive invariance property of the nonnegative orthant R n +,0 ; 2) dissipativity property in the sense of Levinson, see e.g. [5]; 3) ultimate upper and lower bounds for cells populations involved into the model. Our approach is based on the localization method of compact invariant sets.
Here localization refers to description of the location of all compact invariant sets within a chosen domain by utilizing equations and inequalities which depend upon the system parameters. The localization method based on using the first order extremal conditions was proposed in [7] and was subsequently developed and applied for analysis of the Lorenz system [8], and later for study of the global dynamics of various systems within mathematical biology and physics, see e.g. [9,12,13,14,15,16,17].
Intravesical BCG administration is a type of immunotherapy in use to treat superficial bladder cancer for more than 40 years [10]. BCG is an attenuated, nonpathogenic strain of Mycobacterium bovis that was originally used as a vaccine against tuberculosis. In this treatment, bacterial instillations are introduced into the bladder via catheter inserted through the urethra. BCG is internalized and processed by both antigen-presenting cells (APC) and uninfected tumor cells, leading to presentation of the bacterial antigen (Ag) on the tumor surface, attracting APCs which in turn ingest the entire tumor host cell. Once a tumor cell has been ingested, tumor Ag's are presented by the APCs, inducing further recruitment of immune effector cells to the tumor. BCG antigens stimulate an immune response characterized by a surge in pro-inflammatory cytokine levels (such as IL-2) in the infected areas, which is measurable in the urine. Due to inflammatory environment created by the bacterial infection, APCs cause the cytotoxic T-lymphocytes (CTLs) to track bacterial antigen and induce apoptosis of tumor cells according to their affinity to tumor-associated-antigen (TAA). Hence, two types of CTL populations can destroy tumor cells -either via the TAA mechanism, targeting uninfected tumor cells, or via bacteria-associated Ag expressed by infected tumor cells. The addition of exogenous IL-2 is used to exacerbate the local inflammatory environment, which in turn stimulates the maturation of CTLs and prolongs their life span.
The model (1) describes the rates of change in concentrations of molecules or cell populations. The mathematical equations are as follows: Equations (1) describe rates of change in concentrations of molecules or cell populations using the following notations: BCG bacteria within the bladder as B; APCs (dendritic cells (DCs) and macrophages) as A; activated/matured APC's after BCG internalization and processing as A B ; activated/matured APC's specific to tumor Ags as A T ; effector T lymphocytes consisting mostly of CTLs that react to BCG as E B ; effector T lymphocytes consisting mostly of CTLs that react to tumor Ags as E T ; IL-2 units injected inside the bladder as I 2 ; tumor cells infected with BCG as T i ; tumor cells not infected by BCG as T u .
A full explanation of the model is given in [3]; here we present a modified version described below. In the 1st and the 7th equations piecewise constant functions D 1 (t) and D 2 (t) possess a finite number of points of discontinuity. Also, we assume that D 1 (t) > 0, with t ≥ 0, and D 2 (t) ≥ 0, with t ≥ 0. Functions D 1 (t) and D 2 (t) describe processes of instillations of BCG and IL-2 respectively. After BCG instillation bacteria accumulate adjacent to the bladder wall. Upon binding to wall cells, BCG is internalized into the bladder and is processed by APCs at a rate p 1 . The free BCG binds to tumor cells infecting them at a rate p 2 , see the Eq.1 in (1).
When bacteria are present, additional resting APCs (from the closest organs) are recruited at a rate η in response to activated and infected APCs. This parameter also includes the rate at which dead tumor cells infected by bacteria are ingested by APCs. BCG is readily internalized by APCs at a rate p 1 and the APCs become activated by BCG, see the Eq. 2 in (1). The activation of APCs by BCG (A B ) mainly depends on the production rate p 1 , and migration rate β of recruitment of APCs to the draining lymphoid tissues, see the Eq. 3 in (1). The activation of TAA-APC (A T ) is proportional to the number of non-activated APCs (A) and uninfected tumor, with a rate coefficient λ. This term is multiplied by the IL-2 dependent term I 2 I 2 + g I , such that in the absence of IL-2 the production of A T halts, while in the presence of external IL-2 the production term is ∼ 1 representing the migration of activated APCs to the draining lymphoid tissues. Tumor cells are divided into two sub populations: those that have been infected (T i ) with BCG (B) and those that are still uninfected (T u ), see Eq. 8-9 in (1). The growth of the tumor is logistic with a rate r and the carrying capacity K. Effector cells as E B and E T target and destroy infected tumor cells (T i ) at a rate p 4 and uninfected tumor cells (T u ) at a rate α accordingly.
In addition to the model from [3], we add to Eq. 8 the rate of infected tumor cells death rate µ Ti based on empirical findings [18]. Activation of the immune response resulting from encounter between activated APC cells and BCG (A B ) is controlled by parameters β B and β T accordingly, see Eq. 5-6 in (1) describing the rates of change in concentrations of molecules or cell populations.
Model parameters estimated from biological data are summarized in Appendix 1 of [3]. We added to this table the rate µ Ti = 1 2 d −1 . This set of parameters is also used for checking the conditions guaranteeing global tumor clearance conditions and the ultimate bounds. In order to simplify calculations of the model dynamics the variable I 2 is replaced by the constant value obtained from the simulation in Eq. 5-6 of (1). We recall that model equation trajectories of the BCG-IL-2 dynamic system must remain in the positive orthant for all positive times [3]. The current work is dedicated to studies of ultimate dynamics of the model from [3], especially for finding the upper and lower bounds for ultimate densities of interacting cells populations, proving the dissipativity property in the sense of Levinson and finding global asymptotic tumor clearance conditions.
The structure of this article is as follows. The model from [3] is recalled in Section 2. In Section 3, we find ultimate upper bounds for all cells populations in this model. The lower bounds of BCG, APC, activated APC cells and IL-2 concentrations are computed as well. These bounds depend on treatment administration (number of BCG and IL-2 instillations). In Section 4, we obtain convergent sequences of ultimate upper and lower bounds in order to improve bounds presented in Section 3. These bound sequences are obtained with help of localizing functions used in the cyclic way. In Section 5, we introduce an additional convergent sequence of ultimate upper bounds for the infected tumor cells population, further improving bounds. This sequence is derived by increasing ultimate lower bounds for the population. Section 6 describes the conditions under which the system [3] has dissipativity property in the sense of Levinson. In Section 7, we derive conditions under which the ω-limit set of any trajectory in R 9 + is located within the tumor-free plane. In Section 8, we present results of numerical simulation for verification of our computations. Obtained results are then discussed in Section 9.
2. Some preliminaries. In what follows we apply the following notation for variables and parameters in equations (1): and a 1 = µ B , a 12 = a 13 = p 1 , a 19 = a 91 = a 80 = p 2 , Below we assume everywhere that Thus the final form of the model equations is: x 7 x 7 + g 1 , In this section we present some helpful propositions concerning the localization of compact invariant sets, see [7,8,9]. Let us consider a nonlinear systeṁ be a function such that h is not the first integral of the system (4). The function h is exploited in the solution of the localization problem of compact invariant sets and is called a localizing function. By h| U we denote the restriction of h on a set Lie derivative with respect to F. Assume that we are interested in the localization of all compact invariant sets contained in the set U . Further, we define Assertion 1. For any h(x) ∈ C 1 (R n ) all compact invariant sets of the system (4) located in U are contained in the set K(h; U ) defined by the formula as well. If U ∩ S(h) = ∅ then there are no compact invariant sets located in U.
Any of sets K(h; U ) is called a localization set. We notice that the intersection of localization sets is a localization set as well. This observation is refined in contain all compact invariant sets of the system (4) located in U and

Ultimate upper and lower bounds.
Considering the direction of the vector field of the system (3) on the boundary of R 9 +,0 it is evident that this system possesses positive invariance property of R 9 +,0 . Finding the ultimate upper bounds for trajectories of cancer tumor growth models has been already achieved for some other models, see e.g. [15,16,17] to verify the constructed model and to estimate its parameters.
Let max D j (t) = d j max ; min D j (t) = d j min , j = 1, 2. Assertions given in Section 2 allow us to derive upper and lower bounds for the ultimate dynamics. These bounds are calculated using several localizing functions as follows. Firstly, we notice that if D 1 (t) ≡ 0 then all trajectories in R 9 + are attracted to the maximal invariant set contained in the plane x 1 = 0. a) Applying h 1 = x 9 we derive the set S (h 1 ) ∩ {x 9 > 0} given by Thus we have Hence on the set S (h 2 ) we may write and we can take Hence in virtue of (2) we may write on the set S (h 3 ) that and we can take d) Applying h 4 = x 3 we derive the set S (h 4 ) to be given by −a 3 x 3 +a 13 x 1 x 2 = 0. Thus we have the following inequality within e) Applying h 5 = x 4 we derive the set S (h 5 ) to be given by Thus we have the following inequality within Hence we have the following inequality within S (h 3 ) ∩ {x 9 > 0} x 2 ≥ x 2 min := γ a 2 + a 22 x 1 max + a 29 x 9 max .
g) Applying h 6 = x 8 we derive the set S (h 6 ) to be given by Hence we have the following inequality within S (h 6 )∩{x 9 max ≥ x 9 }∩{x 1 ≤ x 1 max } h) Here we apply h 7 = x 5 . Then we derive the set S (h 7 ) to be given by So we have the following inequality within S (h 7 ) i) Here we apply h 8 = x 6 . Then we derive the set S (h 8 ) to be given by So we have the following inequality within S (h 8 ) j) By applying h 9 = x 7 we derive that the inequality holds on S (h 9 ) . So we come to Suppose now that a 70 ≥ a 77 . Then we obtain that the inequality a 7 x 7 ≥ d 2 min is fulfilled on S (h 9 ) . So we come to k) Applying h 2 = x 1 we derive that the inequality Other lower bounds are taken as zero. So we come to Theorem 3.1. All compact invariant sets are located in the polytope P := ∩ 9 j=1 {x j min ≤ x j ≤ x j max } . Remark 1. If d 1 max = 0 hence x 1 max = 0, and x 3 max = 0. Therefore using the function h 7 again we obtain that x 5 max = 0. Despite x 1 max = 0 we derive that x 4 max > 0 and x 6 max > 0. In the case of no BCG treatment, no immune response was associated with BCG instillations (x 3 max = 0 and x 5 max = 0). However, an immune response associated with the tumor is observed (x 4 max > 0 and x 6 max > 0). 4. Iterative construction of ultimate bounds. We start from using x 9 max and x 2 min . a) Applying h 2 = x 1 we have on the set S (h 2 ) and x (1) b) Applying h 1 = x 9 we have on the set S (h 1 ) c) Applying h 3 = x 2 we have on the set S (h 3 ) d) Applying h 6 = x 8 we have on the set S (h 6 ) e) Applying h 3 = x 2 we have on the set S (h 3 ) 1 min x 2 min .
1 min x 3 max := g) Applying h 5 = x 4 we have on the set S (h 5 ) h) Applying h 7 = x 5 we have on the set S (h 7 ) i) Applying h 8 = x 6 we have on the set S (h 8 ) j) Applying h 9 = x 7 we have on the set S (h 9 ) Now using x 2 min instead of x 2 min in (5) we obtain x 1 max . Then using x (1) 9 max instead of x 9 max in (6) we obtain x 8 max . With this goal let us apply the function h 7 = x 5 and derive that the set S (h 7 ) is given by .
As a result, we get the lower bound Next, we exploit the function h 6 = x 8 and derive that the set S (h 8 ) is given by This gives the bound Further, using the expression for x Further, using the expression for x

6.
Dissipativity in the sense of Levinson of the system (3). Firstly, we notice that since the vector field of the system (3) is directed inside the positive orthant at points of its boundary dynamics of (3) is positively invariant. Thus, for each point x ∈ R 9 +,0 = {x j ≥ 0, j = 1, ..., 9} the corresponding positive half trajectory remains in R 9 +,0 . Further, as a population dynamics model the system (3) cannot possess positive half trajectories in R 9 +,0 which escape to infinity because of the boundedness of the available food. The formalization of this property is dissipativity in the sense of Levinson, see e.g. in [5]. We recall that the system (4) is called dissipative in the sense of Levinson if there exists r > 0 such that for any x ∈ R n we have that lim here |x| is the Euclidean norm of x ∈ R n . In this case there exists a bounded set which attracts any point in R n , see [5]. Now our goal is to prove x j x 7 x 7 + g 1 + ξ 8 a 85 x 5 x 8 + (a 96 x 6 x 9 + a 92 x 2 x 9 ) x 7 g 3 (x 7 + g 1 ) (x 9 + g 3 ) .
Then we compute that x j − ξ 8 a 8 x 8 + a 9 x 9 − a 99 x 2 9 − Q, Let us examine conditions ξ 3 a 13 − ξ 1 a 12 − ξ 2 a 22 ≤0, (20) This system of inequalities is consistent. In order to see this, we take any positive ξ 1 , ξ 2 and then we choose positive ξ 8 , ξ 4 , ξ 3 satisfying the 2nd inequality/ the 3rd inequality/ the 1st inequality respectively: Then we define the domain U in R Let C {U } be the complement to U in the orthant R 9 +,0 . By the definition, C {U } is a bounded domain and L f h 10 | U < 0. This implies that any trajectory in R 9 +,0 is attracted to the domain U .
Finally, we derive from Theorem 6.1 that for each point in R 9 +,0 its ω-limit set is a compact invariant set that is not empty. Hence, it is located in the polytope P ⊂ U .
This theorem serves as an important ingredient in the proof of the global asymptotic tumor clearance conditions, see the beginning of the proof of Theorem 7.1.
Theorem 7.1. Supposing that Then the ω− limit set of any positive half trajectory in R 9 + is contained in the tumor-free plane x 9 = 0.
Proof. Since each positive half trajectory in R 9 + is bounded in virtue of Theorem 6.1 we notice that its ω− limit set is a compact invariant set located in the domain P. Therefore each positive half trajectory in R 9 + eventually enters into the domain P and remains there or approaches to its boundary asymptotically. Now we can define the function x 9 , in the domain {x 1 > 0} where positive parameter ζ satisfies the condition ζ ≤ a 99 a −1 19 . Next, we introduce the function and we notice that Then we get that the following inequalities are fulfilled on the domain x 1 − a 91 x 1 + a 1 ζ 2 + a 9 + ζ 2 a 12 x 2 − x 9 a 99 − a 19 ζ 2 − G min ≤h 7 −2ζ d 1 min a 91 + a 1 ζ 2 + a 9 + ζ 2 a 12 x 2 − x 9 a 99 − a 19 ζ 2 − G min ≤h 7 −2ζ d 1 min a 91 + a 1 ζ 2 + a 9 + ζ 2 a 12 x 2 − G min .
Since each positive half trajectory in R 9 + eventually enters into the domain {x 1 > 0; x 2 min ≤ x 2 ≤ x 2 max ; x 7 min ≤ x 7 } as well and remain there or tends to its boundary asymptotically we continue computations in the inequality (24) as follows: Now let us examine the inequality P (ζ) := ζ 2 a 1 + a 12 γ a 2 − 2ζ d 1 min a 91 + a 9 − G min < 0.
In virtue of (21) both of roots of the quadratic polynomial P are real and distinct. By ζ − (ζ + ) we denote the lower (the greater) root of P correspondingly; ζ − > 0. We notice that if then one can find ζ * > 0 such that Taking h 7 with this ζ * we derive from (25) that ;x2≤x2 max }∩S(h7) = 0. Therefore, all compact invariant sets in {x 1 > 0} ∩ {x 2 ≤ x 2 max } are contained in the plane x 9 = 0 as well and we get the desirable conclusion. Now let us consider the condition (26) which entails θ := a 2 √ d 1 min a 91 a 1 a 2 + a 12 γ − a 2 2 d 1 min a 91 (a 1 a 2 + a 12 γ) 2 − a 2 (a 9 − G min ) a 1 a 2 + a 12 γ − a 99 a 19 < 0.
It is easy to check that since (22)-(23) are fulfilled the inequality (27) is satisfied. Thus the proof is completed.
8. Simulation results. In our work, complex biological interactions between the tumor, immune system, and BCG were represented in a nine-dimensional system. As the dimension of the system is relatively high, all calculations were performed with the aid of Matlab subroutines. In order to verify the mathematical model of this biological system, it is widely accepted to estimate the conditions defined in the presented theorems based on parameter sets adopted from the biological and medical literature (these parameters are borrowed from Tables 1 and 2 in [3]).
We examine equations (1) with the following initial conditions: The latter is the initial number of uninfected tumor cells in the urothelium before BCG treatment.
The calculations performed for instillation BCG at 2.2×10 6 c.f.u (colony-forming units, i.e., number of viable bacterial cells) per week, together with 200000 IU (international units) of IL-2. Analysis of global tumor clearance conditions was performed, assuming the growth of cancer cells in the bladder to be logistic.
Upper and lower bounds for ultimate densities of BCG, immune and tumor cell populations derived as in Sections 3 and 4 are presented in Table 1. These bounds are given in a form of a polytope in R 9 +,0 in order to state Theorem 3.1. In Table 1  Upper bounds for ultimate densities of effector cells populations x 5 ; x 6 and cytokine x 7 concentration derived as in Sections 3 and 4 are presented in Table 2. (iii) twenty BCG instillations. 3. BCG+ IL-2 combined therapy: (iv) six BCG and three IL-2 instillations; (v) six BCG and six IL-2 instillations.
In Table 2, effector cells that react to BCG (x 5 ) and effector cells reacting to tumor Ags (x 6 ) in protocol (i) reach the maximum value 6.8 × 10 15 and 7.1 × 10 13 respectively. However in protocol (iv) which includes IL-2 treatment, maximum value of effector cells reaches 6.6 × 10 16 and 5.2 × 10 14 . Treatment efficacy is further improved when the number of IL-2 instillations increases to six: the effector cells population grows significantly to 6.6×10 16 and 5.2×10 14 for x 5 and x 6 respectively. Hence, combined treatment with BCG and IL-2 induces a greater number of effector cells than does the standard BCG treatment alone.
In Section 5 the conditions of global asymptotic tumor clearance were obtained, assuming that only BCG treatment is applied. Here we notice that the inequality (27) is slightly violated, i.e. the value of the expression standing in the left side of (27) is a relatively small positive number under model parameter values taken from the Table 1 in [3].
Computing the left side of (27) we receive values of θ = 0.00044, θ = 0.00015, θ = 0.00009 for protocols (i), (ii) and (iii), respectively. Hence, by increasing the number of BCG instillations during treatment, the probability of cancer cells eradication is increased. Thus, for patients with clinical parameters sufficiently close values to those of model values, BCG treatment may be applied in accordance with (21)-(22)-(23), such that condition (27) is satisfied, achieving asymptotic global tumor clearance. 9. Conclusions. The model was first developed by Bunimovich-Mendrazitsky et al. [3] as an analytical tool to predict the treatment outcome before it is performed, to assist in patient screening and in treatment protocol optimization, especially in problematic cases. Previous mathematically-based simulations [3] individually tailored the BCG+ IL-2 combined therapy to augment low doses of BCG with IL-2 maintenance, resulting in predicted therapy success approaching 95%.
While the results presented in [3] were obtained on the basis of a simulation platform, the current manuscript describes the outcome of analytical methods used to achieve mathematical validation of the model, and to derive the tumor-free equilibrium point, at which cancer cells are effectively eliminated.
The main contribution of the present work lies in the mathematical validation of the nine-dimensional model of bladder cancer immunotherapy, and in obtaining global tumor clearance conditions via the localization method of compact invariant sets [7,8].
In summary: (a) we derived ultimate upper bounds for all components of the state vector of the model (3); (b) we found explicit formulae for conservative ultimate bounds x j min ; j = 1; 2; 3; 5; 7, corresponding to BCG levels, APC population, activated APC population, EB population and complex IL -2 concentration, respectively; (c) we demonstrated the existence of the bounded positively invariant domain. Immune cells populations' bounds are presented in Table 1 The analysis of the model (3) may be used to monitor efficacy of individual patient treatment, based on the conditions required for global tumor clearance. However, the question remains as to whether a biologically significant set of model parameters meets these conditions. The present findings revealed that extremely high doses of BCG+ IL-2 would be required to reach the tumor-free equilibrium point, beyond that which is medically feasible. While the equilibrium point is not currently clinically achievable, these findings nonetheless represent a significant advancement towards the development of a clinical tool for monitoring and improving the efficacy of individualized therapy for bladder cancer with BCG and IL-2.
The derived clearance conditions are based on the assumed treatment schedule and model parameters for an individual patient. Importantly, these clearance conditions are global, i.e. they do not depend upon initial tumor parameters. Thus, a most promising direction for further study is to derive attractivity conditions to the tumor-free equilibrium point under the variation of supply parameters d 1 min and d 1 max within a biologically feasible range and to calculate an attraction domain for such conditions.
Based upon results of the simulation, we conclude that: 1. increasing the number of BCG instillations does not significantly increase the number of effector cells; 2. introduction of IL-2 therapy to the treatment protocol increases the number of effector cells by an order of magnitude.
Importantly, the global tumor clearance conditions are formulated in terms of the parameter d 1 min which is bounded from below, as in formulae (21)-(22)-(23) depending upon the BCG instillation and the biological parameters of the individual patient. These findings support the validity of the described mathematical model of BCG therapy and its optimization for treatment of non-invasive bladder cancer. Therefore the present study may lead to the further development of clinical tools for monitoring and improving the efficacy of individualized therapy for bladder cancer with BCG and IL-2.