Stability Analysis of Delayed Tumor-Antigen-Activated Immune Response in combined BCG and IL-2 Immunotherapy of Bladder Cancer

We use a system biology approach to translate the interaction of Bacillus Calmette-Gurin (BCG) + interleukin 2 (IL-2) for the treatment of bladder cancer into a mathematical model. The model is presented as a system of differential equations with the following variables: number of tumor cells, bacterial cells, immune cells, and cytokines involved in the tumor-immune response. This work investigates the delay effect induced by the proliferation of tumor antigen-specific effector cells after the immune system destroys BCG-infected urothelium cells following BCG and IL-2 immunotherapy in the treatment of bladder cancer. For the proposed model, three equilibrium states are found analytically. The stability of all equilibria is analyzed using the method of Lyapunov functionals construction and the method of linear matrix inequalities (LMIs).


Introduction
Bladder cancer (BC) is the fourth most common cancer in males [1] and the 11th most common cancer in women [2]. The global prevalence of BC is estimated at more than one million and is steadily increasing [1]. The risks of BC appear to vary across world regions, correlating with smoking and occupational exposures to carcinogens [3].
The long-term nature of BC is similar to that of many chronic diseases and often requires invasive and careful long-term follow-up, and in some cases, adjuvant treatment. This translates into high patient costs, making BC one of the most expensive cancers to treat. The preferred treatment for BC depends on its grade at diagnosis. Transurethral resection (TURBT) is the standard primary treatment for BC low-grade (stages Ta, T1 and carcinoma in situ), where cancer growth occurs superficially on the inner surface of the bladder in the form of a polyp but does not extend to the muscle. After TURBT, in malignant case, the treatment is either chemotherapy or immunotherapy for eradication of any residual cancer cells. Chemotherapeutic agents such as mitomycin C, doxorubicin, and epirubicin have long been used as intravesical therapies for BC [3,4]. On basis of these key variables we build the equations in the system (2.1), described 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 Ag 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 ; -transforming growth factor-beta (TGF-β) denotes as F β .
Mathematical and biological interpretation of every equation from the system (2.1) are examined below: BCG dynamics: whereḂ(t) is the dynamical rate of BCG level changes with time. It is comprised of a positive term corresponding to BCG instillations, and of negative terms corresponding to the elimination of BCG by antigen-presenting cells (APCs) according to the rate coefficient p 1 , BCG tumor cell infection at a rate coefficient p 2 , and bacteria cell death with rate coefficient µ B . A quantity b of BCG is instilled into the bladder via a catheter inserted through the urethra once in a week during 6-8 weeks. In this study, we have chosen to simplify the problem by assuming that BCG is introduced into the bladder at a constant rate b. Dynamics of non-activated APC: where dynamic of A(t) is governed by two positive terms and three negative terms. The first positive term describes the normal influx of APCs to the tumor at a constant rate γ. The second positive term describes the recruitment of APCs due to bacterial infection at a rate coefficient η. The first negative term describes the activation of APCs by BCG at the rate coefficient p 1 . The second negative term is natural cell death at the rate coefficient µ A . The last negative term accounts for the two-stage elimination of tumor cells, according to recent knowledge, first by effector CTL activity upon BCG-infected tumor cells, which leads to lysis of these cells and flooding of the tumor micro-environment with tumor antigens. Activation of APC cells with tumor-specific antigens occurs with a delay of τ(t) after the destruction of infected tumor cells. The localized inflammatory response then attracts APCs, such as macrophages, which in turn eliminate uninfected tumor cells, according to the rate θ p 3 . Dynamics of APC activated by BCG: where dynamic of A B (t) is described by one positive term and two negative terms. The positive term is proportional to the numbers of non-activated APCs as well as BCG bacteria, with rate coefficient p 1 . The first negative term is the migration of the infected, activated APCs to the draining lymphoid tissues, at a rate of coefficient β 1 . The second negative term is the death of activated APCs at a rate of coefficient µ A 1 . Dynamics of tumor-Ag-activated APC (TAA-APC): where the dynamic TAA-APC is comprised of one positive term and three negative terms. The positive term describes the APCs which were activated by tumor antigen after eradication of infected tumor cells with the same τ(t) delay function. The first negative term represents the tumor-Ag-activated APCs cells which destroy the uninfected tumor cells, with a rate coefficient λ after τ(t) delay. This term is multiplied by an IL-2-dependent parameter with a saturation constant g I , to propose that in the absence of IL-2, A T production ceases, while in the presence of external IL-2, the production term is close to 1. The second negative term describes the migration of TAA-APC to the draining lymphoid tissues at a rate of coefficient β 1 . The third negative term denotes the natural death of TAA-APC at a rate coefficient µ A 1 . Dynamics of effector CTLs that react to BCG infection: where dynamic of E B (t) is comprised of their migration rate, determined by their creation in the lymph node and subsequent migration to the bladder, inactivation rate, and their death rate. The migration element is proportional to A B and IL-2, with a maximal rate of coefficient β B . This rate is brought to saturation by large numbers of A B , using a Michaelis-Menten saturation function, with Michaelis parameter g. The first negative term is inactivation of effector CTLs via their encounter with infected tumor cells (T i ) at a success rate coefficient p 3 . The second negative term corresponds to the BCG-effector CTL (E B ) cells' natural death rate µ E . Dynamics of effector CTLs that react to tumor Ag: whereĖ T (t) is the dynamic of effector cells reacting to tumor Ag after delay τ(t) time due to the eradication of infected tumor cells. It is comprised of their migration rate, inactivation rate, and death rate. The migration element is proportional to A T and IL-2 with a maximal rate coefficient β T . This rate is brought to saturation by large numbers of A T using a Michaelis-Menten saturation function, with Michaelis parameter g. The first negative term describes the inactivation of effector CTLs via their encounter with uninfected tumor cells (T u ), at success rate coefficient p 3 . The second negative term describes the E T natural death rate, with a rate coefficient µ E . IL-2 dynamics: where IL-2 dynamic is driven by a natural source (with constant rate coefficient q 1 ), an external source i 2 , as well as sink and degradation courses. I 2 is consumed by APCs and CTLs with q 2 rate. The consumption depends on I 2 and is limited in a Michaelis-Menten fashion, with the Michaelis constant g I . I 2 degradates with rate of µ I 2 .
Infected tumor cells:Ṫ where dynamic of T i (t) depend on two mechanisms: the rate of bacterial infection of uninfected tumor cells, (T u ), according to rate coefficient p 2 ; and the elimination of infected tumor cells (T i ) by their interaction with BCG-CTL effector cells (represented by E B ), at rate coefficient p 4 . Uninfected tumor cells: where dynamic of T u (t) is comprised of three processes: natural tumor growth, tumor infection by bacteria and tumor elimination by immune cells. The natural tumor growth is characterized by a maximal growth rate coefficient, r, which is limited by the maximal tumor cell number, K. The first negative term, due to bacterial infection, is characterized by a coefficient rate of p 2 . The second negative term is attributed both to the capture and elimination of T u cells by APCs cells, which were activated by tumor-Ag at rate coefficient λ, and to the activity of TAA-CTL effectors, (E T ), which destroy uninfected tumor cells, (T u ), at a rate coefficient α. Two these processes take place after delay τ(t). The dependence in the equation of T u on F β is decreasing from 1to a T,β with Michaelis constant e T,β [23]. And then there is a multiplication of those terms by an I 2 -dependent Michaelis-Menten term, with Michaelis parameter g I , to propose that in the absence of I 2 , T u cellular death does not occur. Since the tumor produces a variety of mechanisms in the biological settings that curtail the success of effector cell activity, they multiply I 2 I 2 + g I by g T T u + g T , to denote the inversely proportional reduction in effector cell activity rate, such that when T u = 0 the term is equal to 1 and lim although this factor can, in principle, nullify the efficacy of CTLs, this is not observed in cases of interest because T u ≤ K [24]. Dynamic of a transforming growth factor-beta TGF-β is proportional to the tumor cell population, T u , with a β,T and destroyed at a rate of µ β The model is as follows: Here it is supposed that the delay τ(t) is given by the equality τ(t) = ν 0 + ν 1 e −ν 2 t , ν i ≥ 0, i = 0, 1, 2. So, the delay is decreasing and τ(0) = ν 0 + ν 1 , τ(∞) = ν 0 . τ(t) is a time-varying function, representing the delay in immune response following treatment, and expressing the number of effector cells in the cancer region. Delay τ(t) means that the system in the current time moment t depends on her state in the past, i.e., in the time t − τ(t). In particular, the dynamics of non-activated APC in the moment t depends on A(t − τ(t)) (see the first equation in system (2.1)), dynamics of tumor-Ag-activated APC (TAA-APC) depends on A(t − τ(t)) and A T (t − τ(t)) (see the fourth equation in system (2.1)) and so on. The delay is measured in reference to the beginning of BCG treatment (t = 0), with a maximum delay of approximately 10 days.
The influence of BCG tends towards zero over time.
To use the mathematical model (2.1), it is important to estimate ranges of parameters that are realistic and consistent with values from the biological and medical literature. In Table 1, we present a list of all the estimated model parameters and the literature for their estimation. Methods for estimating most of the parameters used here are described in [24][25][26].

Equilibria
The concept of equilibrium refers to the theory of dynamical systems, that is, systems developing in time. Equilibrium means a state of rest in which, in the absence of external influences, the system can be indefinitely. In our case, the investigated dynamic system is the human organism and different equilibria means different state of human health. Equilibria of the model (2.1) are found by setting all derivatives to zero and solving for B,A,A B ,A T ,E B ,E T ,I 2 ,T i ,T u and F β (where the asterisk indicates that the variables are at their steady state values). The system (2.1) have several fixed points, we only need to focus on the non-negative equilibria. Equilibria of the model (2.1) are defined by the system of the algebraic equations that follows from (2.1) by the assumption that Note that the solution of the system (3.1) can be not unique. Let us get some solutions of the system (3.1) in two different situations: b > 0 and b = 0.
Consider the following way to get a solution of the system (3.1), i.e., an equilibrium of the system (2.1) for the "tumor-free" case: (1) From (9') it follows that one of the possible T u is T u = 0. (2) (10') it follows F β = 0 (via T u = 0).
with the solution (see Appendix A.1) From (5') and (7') the system for E * B , I * 2 it follows (via A T = E T = T i = 0) with the solution (see Appendix A.2) .
As a result we obtain a tumor-free equilibrium Consider another way to get equilibria of the system (2.1): From (10') it follows F β = 0 or F β = α β,T µ β K (via T u = 0 or T u = K). As a result we obtain two following equilibria: (2) not tumor-free (T * u = 0)

Remark 1.
Suppose that E B = E T = 0. Then from the Equations (5') and (6') of the system (3.1) it follows that A B = A T = 0. From (7') it follows that I 2 = i 2 /µ I 2 . From (3') it follows that AB = 0. From (1') it follows that A cannot be zero by γ > 0 and from (2') it follows that B is zero by b = 0. So, we obtain again the equilibria E 2 , E 3 .

Centralization and Linearization
Calculating the Jacobian matrices J 1 and J 2 of the initial system (2.1), we obtain (see Appendix A.3) the linear approximation of this system in the forṁ z(t) = Hz(t) + Dz(t − τ(t)), (4.1) where z = {z 1 , ..., z 10 } , H = J 1 and D = J 2 are the matrices of dimension 10 × 10, such that   (4.2) and the nonzero elements a ij and d ij of these matrices respectively are , and where It is clear that stability of the zero solution of the system (4.1), (4.2) is equivalent to a local stability of the equilibrium of the initial system (2.1).

Stability
The stability of the equilibrium means that small deviations from the equilibrium do not greatly affect the system and allow it to remain at state of rest. Instability of equilibrium means that the slightest deviation from equilibrium does not allow the system to return to this state. In this case, the system can go into one of the other equilibria (if such states there are) or it can be completely destroyed. In our case, the investigated dynamic system is the human organism. The task of physicians is to determine the equilibrium favorable for human health and to organize such treatment in which this equilibrium becomes stable. In [16] stability conditions for the Equation (4.1) are obtained in the form of nonlinear matrix Riccati equations. Via Schur complement (see Appendix A.4) similarly to [19,20] these conditions can be reformulated in the form of LMIs: Lemma 1. Put Φ 0 (P) = H P + PH. Ifτ(t) ≤ 0 and for some positive definite matrices P and R at least one of the LMIs holds then the zero solution of the Equation (4.1) is asymptotically stable.

Corollary 1.
If at least one of the LMIs (5.1) holds then the appropriate equilibrium of the system (2.1) is locally asymptotically stable.

Remark 2.
For LMIs (5.1) the matrix H has to be the Hurwitz matrix.

Example 2.
Let be again r = 0.0048 but i 2 = 0, µ E = 0.19 and all other parameters are given in Table 1. Similarly to Example 1 it is shown that the equilibrium E 1 is locally asymptotically stable in the same interval b ∈ B 1 = [2.2 × 10 4 , 58.9 × 10 8 ]. In particular, the equilibria a) b =2.2 × 10 4 ,  are locally asymptotically stable. Example 3. Let be r = 0.0085 and all other parameters as in Example 1. In this case, the equilibrium E 1 is locally asymptotically stable for b ∈ B 2 = [3.6 × 10 4 , 53.9 × 10 8 ]. In particular, the equilibria b =3.6 × 10 4 ,  are locally asymptotically stable. For b ≤ 3.5 × 10 4 and b ≥ 54 × 10 9 the equilibrium E 1 is unstable.

Discussion and Conclusions
In this work we present the improved model of combined therapy BCG + IL-2 immunotherapy for BC. There are two outcomes in the presented model analysis: (1) analytical analysis of stationary system's points and (2) examining the delay influence to stability of these stationary system's points.
A detailed description of the results is presented below: (1). The current manuscript describes the outcome of analytical methods used to derive the equilibria points and especially the tumor-free equilibrium point, at which cancer cells are effectively eliminated. The model demonstrates several equilibria which depend on biologically related parameters and initial conditions.
In our previous works [10,24,42], we provided a local analysis of the stability of the equilibrium states of the model only in simulation ways, because the model consists of 10 equations, which made the analytical analysis of stationary points very difficult. Using the general method of Lyapunov functionals construction [15,16] and the method of linear matrix inequalities (LMIs) [17][18][19][20] gives a possibility for research of the "big" system's stability.
The biological system (bladder is in our case) in the presence of cancer cells without treatment is unstable [25]. The clinical importance of steady state stability stems from the uncertainty typically encountered by the clinician in assessing the exact number of tumor and immune cells present at the start of treatment. In theory, treatment could improve system stability. The physician will be interested in the local stable system with great interest. This means that if a tumor can be localized or eradicated, minor violations of its components will not cancel it.
It is shown that the considered system has three equilibria describing the different states of the patient. Only in the E 1 equilibrium do get cancer cell eradication (T u = 0), meaning successful treatment. Stability analysis of the system (2.1) shows the equilibrium E 1 is stable if the BCG dose is reflected in the condition depend on the growth of cancer cells (as indicated in Examples 1-4).
In equilibrium E 1 we obtain a strong immune response because A B = 10 5 and I 2 = 10 5 that help to arrive at a tumor-free fixed point. By the basic parameters of BCG, maximum tumor size, tumor growth rate, and immune response parameters, we found the BCG dose where E 1 will be stable (Example 1-4). Equilibria E 2 and E 3 were obtained only for IL-2 therapy. E 2 equilibrium do get cancer cell eradication (T u = 0), but the system is not stable in the equilibrium E 2 (as shown in Example 5 and Remark 4). E 3 equilibrium do not get cancer cell eradication (T u = 10 11 ) and system is not stable in the equilibrium too. Hence, we found that administering only IL-2 did not result in the elimination of tumor cells. The same result we received via simulations in the [24].
(2). In this work, for the first time, we consider the real picture of the immune response in the treatment of BCG + IL-2, taking into account the time of formation of effector cells due to BCG infection, which leads to a delay in the response to the destruction of cancer cells. We investigate this influence using the method of Lyapunov functionals and the method of linear matrix inequalities (LMIs). The delay effect during cancer cell eradication does not influence the stability in the equilibrium E 1 . As it is shown in the Remark 3, each element of the matrix D for E 1 before the term with delay or equals zero or close enough to zero. It means that dependence on delay is low enough.
We would like to raise awareness in the community of urological-oncological doctors about the possibilities of mathematical modeling and receive quantitative data to improve this model. The ability to plan and predict by calculating a modulated dose of treatment can benefit patients who are unable to take routine treatment because of its serious side effects, as well as to patients who were previously not considered treatable. It is necessary to note also that three equilibria that are investigated in this work are equilibria obtained from the system (3.1) in an analytical way. So, there is a possibility of continuing stability investigation of the considered model via getting additional equilibria by numerical methods and using additional results of stability theory [16]. Thus, it will be an interest of experts in this direction to the obtained here results and this research will be continued .
where lim |y|→0 |o(y)| |y| = 0, |y| is the Euclidean norm in R n , and the equality F(x * , x * ) = 0, we obtain the linear approximation of the Equation (A2). So, a condition for asymptotic stability of the zero solution of the Equation (A3) is also a condition for local stability of the equilibrium x * of the initial Equation (A1).

Appendix A.4. Schur Complement
Schur complement [43]. The symmetric matrix A B B C is negative definite if and only if C and A − BC −1 B are both negative definite.