Stability Analysis of Delayed Tumor-Antigen-ActivatedImmune Response in Combined BCG and IL-2Immunotherapy of Bladder Cancer

Abstract

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 main goal of this research is to predict the outcome of BCG + IL-2 treatment combinations. We examined whether the delay effect caused by the proliferation of tumor antigen-specific effector cells after the immune system destroys BCG-infected urothelium cells after BCG and IL-2 immunotherapy influences success in bladder cancer treatment. To do this, we introduce a system of differential equations where the variables are the main participants in the immune response after BCG installations to fight cancer: the number of tumor cells, BCG cells, immune cells, and cytokines involved in the tumor-immune response. The relevant parameters describing the dynamics of the system are taken from a variety of biological, clinical literature and estimated using the mathematical models. We examine the local stability analysis of non-negative equilibrium states of the model. In theory, treatment could improve system stability, and we analyze the stability of all equilibria using the method of Lyapunov functionals construction and the method of linear matrix inequalities (LMIs). Our results prove that the period for the proliferation of tumor antigen-specific effector cells does not influence to the success of the non-responsive patients after an intensified combined BCG + IL-2 treatment.

1. 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].

Immunotherapy, BCG, a live attenuated strain of Mycobacterium bovis, has been shown to reduce both recurrence and progression of BC, and therefore represents an important tool in the treatment of BC [5]. BCG was first introduced as an intravesical therapy for BC in 1976 year by Morales and associates [6]. However, despite BCG favorable effects, a significant proportion of patients do not respond to BCG or tolerate treatment. Therefore, research has been pursued and efforts made to improve BCG therapy. During the past decades, cytokine-based therapies have been developed [7]. To date, multiple cytokines, such as IFN-$\alpha$, IL-2 and IL-12 [8,9], have been evaluated, alone, or in combination with BCG for the treatment of BC.

This research is based on the model of BC immunotherapy [10], focusing on the clinical use of BCG and IL-2. During this treatment, the malignant cells are attacked by the patient’s own immune system rather than through external chemical or surgical intervention. BCG is introduced into the bladder where it adheres to the damaged regions of the bladder wall. The cascade of immune response induced by intravesical BCG + IL-2 instillation consists of 2 phases:

I. BCG infects tumor cells, and stimulates the patient’s immune response, which mainly targets BCG infected tumor cells and may potentially eliminate the remaining tumor. At the same time, BCG infection stimulates the antigen (Ag)-presenting cells (APCs) to produce cytotoxic T-lymphocytes (CTL) associated with bacteria Ag. This causes the presentation of bacterial Ag on the tumor surface, which attracts APCs that ingest these cells.

II. Once a tumor cell has been ingested, tumor Ag’s are presented by the APCs. Due to the aforementioned inflammatory environment, created by the bacterial infection, APCs cause the effector cells (the most CTL cells) to either mature and capture tumor cells according to their tumor-associated-antigen (TAA) [11].

In other words, two CTL populations can destroy tumor cells—either via the TAA mechanism (uninfected tumor cells) or via bacteria-associated Ag (infected tumor cells). The addition of IL-2 is expected to create an inflammatory environment that stimulates the maturation of either CTL population [12]. During the second phase of the immune response, there is a delay in the entry of CTL cells into the bladder, because the CTL maturation cells process can take until 7 days. The modeling of biological systems with delays is based on differential equations with delay (Delay Differential Equations—DDE), where the derivative of an unknown function depends on the past states of the system [13,14].

The purpose of the current research is to analyze the stability of BCG model in DDE. Stability of equilibria is analyzed using the classical Lyapunov-Krasovskii functionals method together with the Kolmanovskii-Shaikhet general method of Lyapunov functionals construction [15,16] and the method of linear matrix inequalities (LMIs) [17,18,19,20]. The considered BCG-model is described by a system of nonlinear differential equations with delays. A stability investigation of systems of this type can be reduced to stability investigation of the linear part of a nonlinear system. The obtained sufficient conditions for asymptotic stability of the zero solution of an auxiliary linear system, at the same time, are sufficient conditions for the local stability of the corresponding equilibrium of the initial nonlinear system. Here standard definitions of stability are used (see, for instance [21]).

2. Description of the Model

The model that we analyze in this work describes interactions between BCG, the immune system and tumor cells within the bladder during combined BCG + IL-2 treatment (see Figure 1). BCG induces tumor removal by attaching to the urothelium and causing an inflammatory response there. The inflamed area attracts innate immune cells, which in turn attract CTLs and natural killer cells that attack tumor cells [22]. There are key variables interest us in the BCG + IL-2 treatment’s model: BCG and IL-2, APCs not activate and activated, effector T lymphocytes, tumor cells and tumor cells infected with BCG, transforming growth factor-beta (TGF-$\beta$).

On basis of these key variables we build the equations in the system (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-$\beta$) denotes as $F_{\beta }$.

Mathematical and biological interpretation of every equation from the system (1) are examined below:

BCG dynamics:

$$\dot{B}\left(t\right)=b-p_{1}A\left(t\right)B\left(t\right)-p_{2}B\left(t\right)T_u\left(t\right)-{\mu }_{B}B\left(t\right),$$

where $\dot{B}\left(t\right)$ 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 ${\mu }_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:

$$\dot{A}\left(t\right)=\gamma -(p_1-\eta )A\left(t\right)B\left(t\right)-{\mu }_{A}A\left(t\right)-\theta p_3{E}_B\left(t\right)T_i\left(t\right)A(t-\tau \left(t\right)),$$

where dynamic of $A\left(t\right)$ 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 $\gamma$. The second positive term describes the recruitment of APCs due to bacterial infection at a rate coefficient $\eta$. 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 ${\mu }_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 $\tau \left(t\right)$ 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 $\theta p_3$.

Dynamics of APC activated by BCG:

$${\dot{A}}_B\left(t\right)=p_{1}A\left(t\right)B\left(t\right)-(\beta +{\mu }_{A_1})A_B\left(t\right),$$

where dynamic of $A_B\left(t\right)$ 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 ${\beta }_1$. The second negative term is the death of activated APCs at a rate of coefficient ${\mu }_{A_1}$.

Dynamics of tumor-Ag-activated APC (TAA-APC):

$${\dot{A}}_T\left(t\right)=\theta p_3{E}_B\left(t\right)T_i\left(t\right)A(t-\tau \left(t\right))-\lambda A_T(t-\tau \left(t\right))T_u\left(t\right)\left({\displaystyle \frac{I_2\left(t\right)}{I_2\left(t\right)+g_I}}\right)-(\beta +{\mu }_{A_1})A_T\left(t\right),$$

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 $\tau \left(t\right)$ delay function. The first negative term represents the tumor-Ag-activated APCs cells which destroy the uninfected tumor cells, with a rate coefficient $\lambda$ after $\tau \left(t\right)$ 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 ${\beta }_1$. The third negative term denotes the natural death of TAA-APC at a rate coefficient ${\mu }_{A_1}$.

Dynamics of effector CTLs that react to BCG infection:

$${\dot{E}}_B\left(t\right)={\displaystyle \frac{{\beta }_B{A}_B\left(t\right)I_2\left(t\right)}{A_B\left(t\right)+g}}-p_3{T}_i\left(t\right)E_B\left(t\right)-{\mu }_E{E}_B\left(t\right),$$

where dynamic of $E_B\left(t\right)$ 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 ${\beta }_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 ${\mu }_E$.

Dynamics of effector CTLs that react to tumor Ag:

$${\dot{E}}_T\left(t\right)={\displaystyle \frac{{\beta }_T{A}_T(t-\tau \left(t\right))I_2\left(t\right)}{A_T(t-\tau \left(t\right))+g}}-p_3{T}_u\left(t\right)E_T(t-\tau \left(t\right))-{\mu }_E{E}_T\left(t\right),$$

where ${\dot{E}}_T\left(t\right)$ is the dynamic of effector cells reacting to tumor Ag after delay $\tau \left(t\right)$ 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 ${\beta }_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 ${\mu }_E$.

IL-2 dynamics:

$${\dot{I}}_2\left(t\right)=(A_B\left(t\right)+A_T\left(t\right)+E_B\left(t\right)+E_T\left(t\right))\left(q_1-{\displaystyle \frac{q_2{I}_2\left(t\right)}{I_2\left(t\right)+g_I}}\right)+i_2-{\mu }_{I_2}{I}_2\left(t\right),$$

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 ${\mu }_{I_2}$.

Infected tumor cells:

$${\dot{T}}_i\left(t\right)=p_{2}B\left(t\right)T_u\left(t\right)-p_4{E}_B\left(t\right)T_i\left(t\right),$$

where dynamic of $T_i\left(t\right)$ 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 $\left(T_i\right)$ by their interaction with BCG-CTL effector cells (represented by $E_B$), at rate coefficient $p_4$.

Uninfected tumor cells:

$$\begin{array}{cc}\hfill {\dot{T}}_u\left(t\right)=& r{T}_u\left(t\right)\left(1-{\displaystyle \frac{T_u\left(t\right)}{K}}\right)-p_{2}B\left(t\right)T_u\left(t\right)-(\lambda A_T(t-\tau \left(t\right))T_u\left(t\right)\hfill \\ \hfill & +\alpha E_T(t-\tau \left(t\right))T_u\left(t\right){\displaystyle \frac{{\alpha }_{T,\beta }{F}_{\beta }\left(t\right)+e_{T,\beta }}{F_{\beta }\left(t\right)+e_{T,\beta }}})\left({\displaystyle \frac{I_2\left(t\right)}{I_2\left(t\right)+g_I}}\right)\left({\displaystyle \frac{g_T}{T_u\left(t\right)+g_T}}\right),\hfill \end{array}$$

where dynamic of $T_u\left(t\right)$ 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 $\lambda$, and to the activity of TAA-CTL effectors, ($E_T$), which destroy uninfected tumor cells, ($T_u$), at a rate coefficient $\alpha$. Two these processes take place after delay $\tau \left(t\right)$. The dependence in the equation of $T_u$ on $F_{\beta }$ is decreasing from 1to $a_{T,\beta }$ with Michaelis constant $e_{T,\beta }$[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 ${\displaystyle \frac{I_2}{I_2+g_I}}$ by ${\displaystyle \frac{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 $\underset{T_u\to \infty }{lim}{\displaystyle \frac{g_T}{T_u+g_T}}=0$. Note that although this factor can, in principle, nullify the efficacy of CTLs, this is not observed in cases of interest because $T_u\le K$ [24].

Dynamic of a transforming growth factor-beta TGF-$\beta$ is proportional to the tumor cell population, $T_u$, with $a_{\beta ,T}$ and destroyed at a rate of ${\mu }_{\beta }$

$${\dot{F}}_{\beta }\left(t\right)={\alpha }_{\beta ,T}{T}_u\left(t\right)-{\mu }_{\beta }{F}_{\beta }\left(t\right).$$

The model is as follows:

$$\begin{array}{cc}\hfill \dot{A}\left(t\right)=& \gamma -(p_1-\eta )A\left(t\right)B\left(t\right)-{\mu }_{A}A\left(t\right)-\theta p_3{E}_B\left(t\right)T_i\left(t\right)A(t-\tau \left(t\right)),\hfill \\ \hfill \dot{B}\left(t\right)=& b-p_{1}A\left(t\right)B\left(t\right)-p_{2}B\left(t\right)T_u\left(t\right)-{\mu }_{B}B\left(t\right),\hfill \\ \hfill {\dot{A}}_B\left(t\right)=& p_{1}A\left(t\right)B\left(t\right)-(\beta +{\mu }_{A_1})A_B\left(t\right),\hfill \\ \hfill {\dot{A}}_T\left(t\right)=& \theta p_3{E}_B\left(t\right)T_i\left(t\right)A(t-\tau \left(t\right))-\lambda A_T(t-\tau \left(t\right))T_u\left(t\right)\left({\displaystyle \frac{I_2\left(t\right)}{I_2\left(t\right)+g_I}}\right)-(\beta +{\mu }_{A_1})A_T\left(t\right),\hfill \\ \hfill {\dot{E}}_B\left(t\right)=& {\displaystyle \frac{{\beta }_B{A}_B\left(t\right)I_2\left(t\right)}{A_B\left(t\right)+g}}-p_3{T}_i\left(t\right)E_B\left(t\right)-{\mu }_E{E}_B\left(t\right),\hfill \\ \hfill {\dot{E}}_T\left(t\right)=& {\displaystyle \frac{{\beta }_T{A}_T(t-\tau \left(t\right))I_2\left(t\right)}{A_T(t-\tau \left(t\right))+g}}-p_3{T}_u\left(t\right)E_T(t-\tau \left(t\right))-{\mu }_E{E}_T\left(t\right),\hfill \\ \hfill {\dot{I}}_2\left(t\right)=& (A_B\left(t\right)+A_T\left(t\right)+E_B\left(t\right)+E_T\left(t\right))\left(q_1-{\displaystyle \frac{q_2{I}_2\left(t\right)}{I_2\left(t\right)+g_I}}\right)+i_2-{\mu }_{I_2}{I}_2\left(t\right),\hfill \\ \hfill {\dot{T}}_i\left(t\right)=& p_{2}B\left(t\right)T_u\left(t\right)-p_4{E}_B\left(t\right)T_i\left(t\right),\hfill \\ \hfill {\dot{T}}_u\left(t\right)=& r{T}_u\left(t\right)\left(1-{\displaystyle \frac{T_u\left(t\right)}{K}}\right)-p_{2}B\left(t\right)T_u\left(t\right)-(\lambda A_T(t-\tau \left(t\right))T_u\left(t\right)\hfill \\ \hfill & +\alpha E_T(t-\tau \left(t\right))T_u\left(t\right){\displaystyle \frac{{\alpha }_{T,\beta }{F}_{\beta }\left(t\right)+e_{T,\beta }}{F_{\beta }\left(t\right)+e_{T,\beta }}})\left({\displaystyle \frac{I_2\left(t\right)}{I_2\left(t\right)+g_I}}\right)\left({\displaystyle \frac{g_T}{T_u\left(t\right)+g_T}}\right),\hfill \\ \hfill {\dot{F}}_{\beta }\left(t\right)=& {\alpha }_{\beta ,T}{T}_u\left(t\right)-{\mu }_{\beta }{F}_{\beta }\left(t\right).\hfill \end{array}$$

(1)

Here it is supposed that the delay $\tau \left(t\right)$ is given by the equality $\tau \left(t\right)={\nu }_0+{\nu }_1{e}^{-{\nu }_{2}t}$, ${\nu }_i\ge 0$, $i=0,1,2$. So, the delay is decreasing and $\tau \left(0\right)={\nu }_0+{\nu }_1$, $\tau (\infty )={\nu }_0$. $\tau \left(t\right)$ 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 $\tau \left(t\right)$ means that the system in the current time moment t depends on her state in the past, i.e., in the time $t-\tau \left(t\right)$. In particular, the dynamics of non-activated APC in the moment t depends on $A(t-\tau (t\left)\right)$ (see the first equation in system (1)), dynamics of tumor-Ag-activated APC (TAA-APC) depends on $A(t-\tau (t\left)\right)$ and $A_T(t-\tau \left(t\right))$ (see the fourth equation in system (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 (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. List of all parameters.

Parameters	Physical Interpretation (Units)	Estimated Value	Reference
${\mu }_A$	APC half life $\left[day{s}^{-1}\right]$	0.038	[27]
${\mu }_{A_1}$	Activated APC half life $\left[day{s}^{-1}\right]$	0.138	[28]
${\mu }_{E_1}$	Effector cells mortality rate w/o IL-2 $\left[day{s}^{-1}\right]$	0.19	[23] and calculated
${\mu }_{E_2}$	Effector cells mortality rate with IL-2 $\left[day{s}^{-1}\right]$	0.034	[29]
${\mu }_B$	BCG half life $\left[day{s}^{-1}\right]$	0.1	[30]
$p_1$	The rate of BCG binding with APC $\left[cell{s}^{-1}\right]\left[day{s}^{-1}\right]$	$1.25\times {10}^{-4}$	[31] adjusted for liters
$p_2$	Infection rate of tumor cells by BCG $\left[cell{s}^{-1}\right]\left[day{s}^{-1}\right]$	$0.028\times {10}^{-6}$	From model simulation
$p_3$	Rate of E deactivation after binding with infected tumor cells $\left[cell{s}^{-1}\right]\left[day{s}^{-1}\right]$	$1.03\times {10}^{-10}$	[32]
$p_4$	Rate of destruction of infected tumor cells by effector cells $\left[cell{s}^{-1}\right]\left[day{s}^{-1}\right]$	$1.1\times {10}^{-6}$	[32]
$\lambda$	Production rate of TAA-APC $\left[day{s}^{-1}\right]$	${10}^{-8}$	[33]
${\beta }_B$	Recruitment rate of effector cells in response to signals released by BCG-infected and activated APC $\left[cell{s}^{-1}\right]\left[day{s}^{-1}\right]\left[I_2^{-1}\right]$	$1.45\times {10}^8$	[34]
${\beta }_T$	Recruitment rate of effector cells in response to signals released by TAA-infected and activated APC $\left[cell{s}^{-1}\right]\left[day{s}^{-1}\right]\left[I_2^{-1}\right]$	$1.514\times {10}^6$	[35]
$\gamma$	Initial APC cell numbers $\left[cell{s}^{-1}\right]\left[day{s}^{-1}\right]$	4700	[28]
$\eta$	Rate of recruited additional resting APCs $\left[cell{s}^{-1}\right]\left[day{s}^{-1}\right]$	$2.8\times {10}^{-6}$	[27]
r	Tumor growth rate $\left[day{s}^{-1}\right]$	0.0048–0.0085	[36]
b	Bio-effective dose of BCG [c.f.u./week]	$2.2\times {10}^8$	From clinical data provided by Dr. Sarel Halachmi
$\beta$	Migration rate of TAA-APC and bacteria activated APC to the lymph node $\left[cell{s}^{-1}\right]\left[day{s}^{-1}\right]$	0.034	[27]
$\alpha$	Efficacy of an effector cell on tumor cell $\left[cell{s}^{-1}\right]\left[day{s}^{-1}\right]$	$3.7\times {10}^{-6}$	[37]
g	Michaelis-Menten constant for BCG activated CTLs and for TAA-CTLs[cells]	${10}^{13}$	From model simulation
$g_T$	Michaelis-Menten constant for tumor cells[cells]	5200	[24]
K	Maximal tumor cell population [cells]	${10}^{11}$	[38]
$q_1$	Rate of IL-2 production IU $\left[cell{s}^{-1}\right]\left[day{s}^{-1}\right]$	0.007	[39] and simulations
$q_2$	The proportion of IL-2 used for differentiation of effector cells IU $\left[cell{s}^{-1}\right]\left[day{s}^{-1}\right]$	$1.2\times {10}^{-3}$	[35]
${\mu }_{I_2}$	Degradation rate $\left[day{s}^{-1}\right]$	11.5	[35,40]
$\theta$	Recruitment rate of Tumor-Ag-activated APC cells in response to signals released after binding effector cells, that react to BCG infection, with infected tumor cells [$1/cel{l}^{-1}$]	0.01	From model simulation
${\alpha }_{\beta ,T}$	The release term per tumor cell [$pg/cel{l}^{-1}\ast d^{-1}$]	$1.38\times {10}^{-4}$	[23]
${\alpha }_{T,\beta }$	Michaelis-Menten saturation dynamics. The dependence on $F_{\beta }$ is decreasing from 0 to ${\alpha }_{T,\beta }$ [none]	0.69	[23]
$e_{T,\beta }$	Michaelis constant [$pg$]	10000	[23]
${\mu }_{\beta }$	The constant rate, accounts for degradation of $F_{\beta }$ [$d^{-1}$]	166.32	[23]
$g_I$	Michaelis-Menten constant for IL-2 [cells]	10000	From model calculations
$i_2$	Rate of external source [units per treatment]	$8\times {10}^5$–$7.7\times {10}^6$	[9]

, 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].

3. 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 (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_{\beta }$ (where the asterisk indicates that the variables are at their steady state values). The system (1) have several fixed points, we only need to focus on the non-negative equilibria.

Equilibria of the model (1) are defined by the system of the algebraic equations

$$\begin{array}{c}\left(1^{\prime }\right)\gamma =[(p_1-\eta )B+{\mu }_A+\theta p_3{E}_B{T}_i]A,\hfill \\ \left(2^{\prime }\right)b=(p_{1}A+p_2{T}_u+{\mu }_B)B,\hfill \\ \left(3^{\prime }\right)p_{1}AB=(\beta +{\mu }_{A_1})A_B,\hfill \\ \left(4^{\prime }\right)\theta p_3{E}_B{T}_{i}A=\lambda A_T{T}_u\left({\displaystyle \frac{I_2}{I_2+g_I}}\right)+(\beta +{\mu }_{A_1})A_T,\hfill \\ \left(5^{\prime }\right){\displaystyle \frac{{\beta }_B{A}_B{I}_2}{A_B+g}}=(p_3{T}_i+{\mu }_E)E_B,\hfill \\ \left(6^{\prime }\right){\displaystyle \frac{{\beta }_T{A}_T{I}_2}{A_T+g}}=(p_3{T}_u+{\mu }_E)E_T,\hfill \\ \left(7^{\prime }\right)(A_B+A_T+E_B+E_T)\left(q_1-{\displaystyle \frac{q_2{I}_2}{I_2+g_I}}\right)+i_2={\mu }_{I_2}{I}_2,\hfill \\ \left(8^{\prime }\right)p_{2}B{T}_u=p_4{E}_B{T}_i,\hfill \\ \left(9^{\prime }\right)r{T}_u\left(1-{\displaystyle \frac{T_u}{K}}\right)=p_{2}B{T}_u+T_u\left(\lambda A_T+\alpha E_T{\displaystyle \frac{{\alpha }_{T,\beta }{F}_{\beta }+e_{T,\beta }}{F_{\beta }+e_{T,\beta }}}\right)\left({\displaystyle \frac{I_2}{I_2+g_I}}\right)\left({\displaystyle \frac{g_T}{T_u+g_T}}\right),\hfill \\ \left({10}^{\prime }\right){\alpha }_{\beta ,T}{T}_u={\mu }_{\beta }{F}_{\beta },\hfill \end{array}$$

(2)

that follows from (1) by the assumption that $A\left(t\right),B\left(t\right),A_B\left(t\right),A_T\left(t\right),E_B\left(t\right),E_T\left(t\right),I_2\left(t\right),T_i\left(t\right),T_u\left(t\right),F_{\beta }\left(t\right)$ are constants.

Note that the solution of the system (2) can be not unique. Let us get some solutions of the system (2) in two different situations: $b>0$ and $b=0$.

3.1. Equilibrium with $b>0$, $i_2\ge 0$

Consider the following way to get a solution of the system (2), i.e., an equilibrium of the system (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_{\beta }=0$ (via $T_u=0$).

(3)

From (8’) it follows $E_B{T}_i=0$ (via $T_u=0$).

(4)

From (4’) it follows $A_T=0$ (via $T_u=0$ and $E_B{T}_i=0$).

(5)

From (6’) it follows $E_T=0$ (via $T_u=0$ and $A_T=0$).

(6)

From (1’), (2’) the system for A, B it follows (via $E_B{T}_i=0$ and $T_u=0$)

$$\begin{array}{cc}\hfill & (p_{1}A+{\mu }_B)B=b,\\ \hfill & [(p_1-\eta )B+{\mu }_A]A=\gamma ,\end{array}$$

(3)

with the solution (see Appendix A.1)

$$\begin{array}{c}{A}^{*}={\displaystyle \frac{\sqrt{a_1^2+4a_0{a}_2}-a_1}{2a_0}},B^{*}={\displaystyle \frac{b}{p_1{A}^{*}+{\mu }_B}},\hfill \\ a_0=p_1{\mu }_A,a_1=b(p_1-\eta )+{\mu }_A{\mu }_B-\gamma p_1,a_2=\gamma {\mu }_B.\hfill \end{array}$$

(7)

From (3’) it follows $A_B^{*}={\displaystyle \frac{p_1}{\beta +{\mu }_{A_1}}}{A}^{*}{B}^{*}$ (via $A^{*}$, $B^{*}$).

(8)

From (5’) it follows that if $E_B=0$ then $I_2=0$ but via (7’) it is impossible. So, from $E_B{T}_i=0$ it follows $T_i=0$.

(9)

From (5’) and (7’) the system for $E_B^{*}$, $I_2^{*}$ it follows (via $A_T=E_T=T_i=0$)

$$\begin{array}{c}{\displaystyle \frac{{\beta }_B{A}_B^{*}{I}_2}{A_B^{*}+g}}={\mu }_E{E}_B,(A_B^{*}+E_B)\left(q_1-{\displaystyle \frac{q_2{I}_2}{I_2+g_I}}\right)+i_2={\mu }_{I_2}{I}_2,\hfill \end{array}$$

(4)

with the solution (see Appendix A.2)

$$\begin{array}{c}{I}_2^{*}={\displaystyle \frac{\sqrt{c_1^2+4c_0{c}_2}-c_1}{2c_0}},E_B^{*}=\nu I_2^{*},\hfill \\ c_0={\mu }_{I_2}-\nu (q_1-q_2),c_1=({\mu }_{I_2}-\nu q_1)g_I-i_2-A_B^{*}(q_1-q_2),\hfill \\ c_2=(i_2+A_B^{*}{q}_1)g_I,\nu ={\displaystyle \frac{{\beta }_B{A}^{*}}{{\mu }_E(A_B^{*}+g)}}.\hfill \end{array}$$

As a result we obtain a tumor-free equilibrium

$$\begin{array}{cc}\hfill E_1=& (A^{*},B^{*},A_B^{*},A_T^{*},E_B^{*},E_T^{*},I_2^{*},T_i^{*},T_u^{*},F_{\beta }^{*})\hfill \\ \hfill =& (A^{*},B^{*},A_B^{*},0,E_B^{*},0,I_2^{*},0,0,0).\hfill \end{array}$$

(5)

3.2. Equilibria with $b=0$, $i_2\ge 0$

Consider another way to get equilibria of the system (1):

(1)

From (2’) it follows $B=0$.

(2)

From (3’) it follows $A_B=0$ (via $B=0$).

(3)

From (5’) it follows $E_B=0$ (via $A_B=0$).

(4)

From (4’) it follows $A_T=0$ (via $E_B=0$).

(5)

From (6’) it follows $E_T=0$ (via $A_T=0$).

(6)

From (1’) it follows $A={\displaystyle \frac{\gamma }{{\mu }_A}}$ (via $B=E_B=0$).

(7)

From (7’) it follows $I_2={\displaystyle \frac{i_2}{{\mu }_{I_2}}}$ (via $A_B=A_T=E_B=E_T=0$).

(8)

From (9’) it follows $T_u=0$ or $T_u=K$ (via $B=A_T=E_T=0$).

(9)

From (10’) it follows $F_{\beta }=0$ or $F_{\beta }={\displaystyle \frac{{\alpha }_{\beta ,T}}{{\mu }_{\beta }}}K$ (via $T_u=0$ or $T_u=K$).

(10)

From (8’) it follows $T_i=C=const$ (via $B=E_B=0$).

As a result we obtain two following equilibria:

(1)

tumor-free ($T_u^{*}=0$)

$$\begin{array}{cc}\hfill E_2=& (A^{*},B^{*},A_B^{*},A_T^{*},E_B^{*},E_T^{*},I_2^{*},T_i^{*},T_u^{*},F_{\beta }^{*})\hfill \\ \hfill =& \left({\displaystyle \frac{\gamma }{{\mu }_A}},0,0,0,0,0,{\displaystyle \frac{i_2}{{\mu }_{I_2}}},C,0,0\right),\hfill \end{array}$$

(6)

(2)

not tumor-free ($T_u^{*}\ne 0$)

$$\begin{array}{cc}\hfill E_3=& (A^{*},B^{*},A_B^{*},A_T^{*},E_B^{*},E_T^{*},I_2^{*},T_i^{*},T_u^{*},F_{\beta }^{*})\hfill \\ \hfill =& \left({\displaystyle \frac{\gamma }{{\mu }_A}},0,0,0,0,0,{\displaystyle \frac{i_2}{{\mu }_{I_2}}},C,K,{\displaystyle \frac{{\alpha }_{\beta ,T}}{{\mu }_{\beta }}}K\right).\hfill \end{array}$$

(7)

Remark 1.

Suppose that $E_B=E_T=0$. Then from the Equations (5’) and (6’) of the system (2) it follows that $A_B=A_T=0$. From (7’) it follows that $I_2=i_2/{\mu }_{I_2}$. From (3’) it follows that $AB=0$. From (1’) it follows that A cannot be zero by $\gamma >0$ and from (2’) it follows that B is zero by $b=0$. So, we obtain again the equilibria $E_2$, $E_3$.

4. Centralization and Linearization

Calculating the Jacobian matrices $J_1$ and $J_2$ of the initial system (1), we obtain (see Appendix A.3) the linear approximation of this system in the form

$$\begin{array}{c}\dot{z}\left(t\right)=Hz\left(t\right)+Dz(t-\tau \left(t\right)),\hfill \end{array}$$

(8)

where $z={\{z_1,\dots ,z_{10}\}}^{\prime }$, $H=J_1$ and $D=J_2$ are the matrices of dimension $10\times 10$, such that

$$\begin{array}{c}H=\left[\begin{array}{cccccccccc}{a}_{11}& a_{12}& 0& 0& a_{15}& 0& 0& a_{18}& 0& 0\\ a_{21}& a_{22}& 0& 0& 0& 0& 0& 0& a_{29}& 0\\ a_{31}& a_{32}& a_{33}& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& a_{44}& a_{45}& 0& a_{47}& a_{48}& a_{49}& 0\\ 0& 0& a_{53}& 0& a_{55}& 0& a_{57}& a_{58}& 0& 0\\ 0& 0& 0& 0& 0& a_{66}& a_{67}& 0& a_{69}& 0\\ 0& 0& a_{73}& a_{74}& a_{75}& a_{76}& a_{77}& 0& 0& 0\\ 0& a_{82}& 0& 0& a_{85}& 0& 0& a_{88}& a_{89}& 0\\ 0& a_{92}& 0& 0& 0& 0& a_{97}& 0& a_{99}& a_{9,10}\\ 0& 0& 0& 0& 0& 0& 0& 0& a_{10,9}& a_{10,10}\end{array}\right],D=\left[\begin{array}{cccccccccc}{d}_{11}& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ d_{41}& 0& 0& d_{44}& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& d_{64}& 0& d_{66}& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& d_{94}& 0& d_{96}& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right],\hfill \end{array}$$

(9)

and the nonzero elements $a_{ij}$ and $d_{ij}$ of these matrices respectively are

$$\begin{array}{cc}\hfill a_{11}=& -[(p_1-\eta )B^{*}+{\mu }_A],a_{12}=-(p_1-\eta )A^{*},a_{15}=-\theta p_3{A}^{*}{T}_i^{*},a_{18}=-\theta p_3{A}^{*}{E}_B^{*},\hfill \\ \hfill a_{21}=& -p_1{B}^{*},a_{22}=-(p_1{A}^{*}+p_2{T}_u^{*}+{\mu }_B),a_{29}=-p_2{B}^{*},\hfill \\ \hfill a_{31}=& p_1{B}^{*},a_{32}=p_1{A}^{*},a_{33}=-(\beta +{\mu }_{A_1}),\hfill \\ \hfill a_{44}=& -(\beta +{\mu }_{A_1}),a_{45}=\theta p_3{A}^{*}{T}_i^{*},a_{47}=-{\displaystyle \frac{\lambda g_I{A}_T^{*}{T}_u^{*}}{{(I_2^{*}+g_I)}^2}},a_{48}=\theta p_3{A}^{*}{E}_B^{*},a_{49}=-{\displaystyle \frac{\lambda A_T^{*}{I}_2^{*}}{I_2^{*}+g_I}},\hfill \\ \hfill a_{53}=& {\displaystyle \frac{{\beta }_{B}g{I}_2^{*}}{{(A_B^{*}+g)}^2}},a_{55}=-(p_3{T}_i^{*}+{\mu }_E),a_{57}={\displaystyle \frac{{\beta }_B{A}_B^{*}}{A_B^{*}+g}},a_{58}=-p_3{E}_B^{*},\hfill \\ \hfill a_{66}=& -{\mu }_E,a_{67}={\displaystyle \frac{{\beta }_T{A}_T^{*}}{A_T^{*}+g}},a_{69}=-p_3{E}_T^{*},\hfill \\ \hfill a_{73}=& a_{74}=a_{75}=a_{76}=q_1-{\displaystyle \frac{q_2{I}_2^{*}}{I_2^{*}+g_I}},a_{77}=-\left({\mu }_{I_2}+q_2{g}_I{\displaystyle \frac{A_B^{*}+A_T^{*}+E_B^{*}+E_T^{*}}{{(I_2^{*}+g_I)}^2}}\right),\hfill \\ \hfill a_{82}=& p_2{T}_u^{*},a_{85}=-p_4{T}_i^{*},a_{88}=-p_4{E}_B^{*},a_{89}=p_2{B}^{*},\hfill \\ \hfill a_{92}=& -p_2{T}_u^{*},a_{97}=-{\displaystyle \frac{Q{T}_u^{*}{g}_I}{I_2^{*}+g_I}},a_{99}=-\left(p_2{B}^{*}+r\left({\displaystyle \frac{2T_u^{*}}{K}}-1\right)+{\displaystyle \frac{Q{I}_2^{*}{g}_T}{T_u^{*}+g_T}}\right),\hfill \\ \hfill a_{9,10}=& \alpha E_T^{*}{T}_u^{*}{I}_2^{*}R{\displaystyle \frac{(1-{\alpha }_{T,\beta })e_{T,\beta }}{{(F_{\beta }^{*}+e_{T,\beta })}^2}},a_{10,9}={\alpha }_{\beta ,T},a_{10,10}=-{\mu }_{\beta },\hfill \end{array}$$

(10)

and

$$\begin{array}{cc}{d}_{11}=\hfill & -\theta p_3{E}_B^{*}{T}_i^{*},d_{41}=\theta p_3{E}_B^{*}{T}_i^{*},d_{44}=-{\displaystyle \frac{\lambda g{T}_u^{*}{I}_2^{*}}{I_2^{*}+g_I}},\hfill \\ d_{64}=\hfill & {\displaystyle \frac{{\beta }_{T}g{I}_2^{*}}{{(A_T^{*}+g)}^2}},d_{66}=-p_3{T}_u^{*},d_{94}=-\lambda T_u^{*}{I}_2^{*}R,d_{96}=-\alpha T_u^{*}{I}_2^{*}R{\displaystyle \frac{{\alpha }_{T,\beta }{F}_{\beta }^{*}+e_{T,\beta }}{F_{\beta }^{*}+e_{T,\beta }}}\hfill \end{array}$$

(11)

where

$$\begin{array}{cc}Q=\hfill & \left(\lambda A_T^{*}+\alpha E_T^{*}{\displaystyle \frac{{\alpha }_{T,\beta }{F}_{\beta }^{*}+e_{T,\beta }}{F_{\beta }^{*}+e_{T,\beta }}}\right)R,R=\left({\displaystyle \frac{1}{I_2^{*}+g_I}}\right)\left({\displaystyle \frac{g_T}{T_u^{*}+g_T}}\right).\end{array}$$

(12)

It is clear that stability of the zero solution of the system (8), (9) is equivalent to a local stability of the equilibrium of the initial system (1).

5. 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 (8) 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 ${\Phi }_0\left(P\right)=H^{\prime }P+PH$. If $\dot{\tau }\left(t\right)\le 0$ and for some positive definite matrices P and R at least one of the LMIs

$$\left[\begin{array}{cc}\hfill {\Phi }_0\left(P\right)+R\hfill & \hfill PD\hfill \\ \hfill D^{\prime }P\hfill & \hfill -R\hfill \end{array}\right]<0,\left[\begin{array}{cc}\hfill {\Phi }_0\left(P\right)+R\hfill & \hfill D^{\prime }P\hfill \\ \hfill PD\hfill & \hfill -R\hfill \end{array}\right]<0,\left[\begin{array}{cc}\hfill {\Phi }_0\left(P\right)+D^{\prime }RD\hfill & \hfill P\hfill \\ \hfill P\hfill & \hfill -R\hfill \end{array}\right]<0,$$

(13)

holds then the zero solution of the Equation (8) is asymptotically stable.

Corollary 1.

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

Remark 2.

For LMIs (13) the matrix H has to be the Hurwitz matrix.

In the following examples, we check the stability of E1-E3 using the standard (from Table 1) and reduced (up to one third of the usual dose) BCG dose, taking into account Lamm’s [41] experience in dose reduction to avoid side effects. We distinguished between patients in the examples below using a clinically relevant range of tumor growth rate r.

Example 1.

Let be $r=0.0048$, $i_2=5\times {10}^6$, ${\mu }_E=0.034$ and all other parameters are given in Table 1. Via the LMIs (13) and MATLAB it is shown that the equilibrium $E_1$ is locally asymptotically stable for $b\in B_1=[2.2\times {10}^4,58.9\times {10}^8]$. In particular, the equilibria

$$\begin{array}{cc}\hfill \left(a\right)b=& 2.2\times {10}^4,\hfill \\ \hfill E_1=& (A^{*},B^{*},A_B^{*},A_T^{*},E_B^{*},E_T^{*},I_2^{*},T_i^{*},T_u^{*},F_{\beta }^{*})\hfill \\ \hfill =& (223.2,1.72\times {10}^5,2.790\times {10}^4,0,4.14\times {10}^4,0,4.348\times {10}^5,0,0,0)\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill \left(b\right)b=& 6\times {10}^6,\hfill \\ \hfill E_1=& (A^{*},B^{*},A_B^{*},A_T^{*},E_B^{*},E_T^{*},I_2^{*},T_i^{*},T_u^{*},F_{\beta }^{*})\hfill \\ \hfill =& (0.6415,5.995\times {10}^7,2.795\times {10}^4,0,118.96,0,4.348\times {10}^5,0,0,0)\hfill \end{array}$$

are locally asymptotically stable. For $b\le 2.1\times {10}^4$ and $b\ge 59\times {10}^8$ the equilibrium $E_1$ is unstable.

Example 2.

Let be again $r=0.0048$ but $i_2=0$, ${\mu }_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\in B_1=[2.2\times {10}^4,58.9\times {10}^8]$. In particular, the equilibria

$$\begin{array}{cc}\hfill \left(a\right)b=& 2.2\times {10}^4,\hfill \\ \hfill E_1=& (A^{*},B^{*},A_B^{*},A_T^{*},E_B^{*},E_T^{*},I_2^{*},T_i^{*},T_u^{*},F_{\beta }^{*})\hfill \\ \hfill =& (223.2,1.72\times {10}^5,2.790\times {10}^4,0,0.289,0,16.98,0,0,0)\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill \left(b\right)b=& 6\times {10}^6,\hfill \\ \hfill E_1=& (A^{*},B^{*},A_B^{*},A_T^{*},E_B^{*},E_T^{*},I_2^{*},T_i^{*},T_u^{*},F_{\beta }^{*})\hfill \\ \hfill =& (0.6415,5.995\times {10}^7,2.795\times {10}^4,0,8.328\times {10}^{-4},0,17.01,0,0,0)\hfill \end{array}$$

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\in B_2=[3.6\times {10}^4,53.9\times {10}^8]$. In particular, the equilibria

$$\begin{array}{cc}\hfill b=& 3.6\times {10}^4,\hfill \\ \hfill E_1=& (A^{*},B^{*},A_B^{*},A_T^{*},E_B^{*},E_T^{*},I_2^{*},T_i^{*},T_u^{*},F_{\beta }^{*})\hfill \\ \hfill =& (123.16,3.12\times {10}^5,2.792\times {10}^4,0,2.284\times {10}^4,0,4.348\times {10}^5,0,0,0)\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill b=& 6\times {10}^6,\hfill \\ \hfill E_1=& (A^{*},B^{*},A_B^{*},A_T^{*},E_B^{*},E_T^{*},I_2^{*},T_i^{*},T_u^{*},F_{\beta }^{*})\hfill \\ \hfill =& (0.6415,5.995\times {10}^7,2.795\times {10}^4,0,118.96,0,4.348\times {10}^5,0,0,0).\hfill \end{array}$$

are locally asymptotically stable. For $b\le 3.5\times {10}^4$ and $b\ge 54\times {10}^9$ the equilibrium $E_1$ is unstable.

Example 4.

Let be again $r=0.0085$ but $i_2=0$, ${\mu }_E=0.19$ and all other parameters as in Example 3. In this case the equilibrium $E_1$ is locally asymptotically stable for $b\in B_2=[3.6\times {10}^4,53.9\times {10}^8]$. In particular, the equilibria

$$\begin{array}{cc}\hfill b=& 3.6\times {10}^4,\hfill \\ \hfill E_1=& (A^{*},B^{*},A_B^{*},A_T^{*},E_B^{*},E_T^{*},I_2^{*},T_i^{*},T_u^{*},F_{\beta }^{*})\hfill \\ \hfill =& (123.16,3.12\times {10}^5,2.792\times {10}^4,0,0.1597,0,16.99,0,0,0)\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill b=& 6\times {10}^6,\hfill \\ \hfill E_1=& (A^{*},B^{*},A_B^{*},A_T^{*},E_B^{*},E_T^{*},I_2^{*},T_i^{*},T_u^{*},F_{\beta }^{*})\hfill \\ \hfill =& (0.6415,5.995\times {10}^7,2.795\times {10}^4,0,8.328\times {10}^{-4},0,17.01,0,0,0)\hfill \end{array}$$

are locally asymptotically stable.

Remark 3.

Note that in the equilibrium $E_1$ we have $A_T^{*}=E_T^{*}=T_i^{*}=T_u^{*}=F_{\beta }^{*}=0$. From (11) and (12) it follows that $d_{11}=d_{14}=d_{44}=d_{66}=d_{94}=d_{96}=0$. If $i_2=5\times {10}^6$ then $d_{64}=0.0658$ for both $r=0.0048$ and $r=0.0085$. If $i_2=0$ then $d_{64}=2.575\times {10}^{-6}$ for $r=0.0085$ and $d_{64}=0$ for $r=0.0048$. So, each element of the matrix D or equals zero or close enough to zero. It means that dependence on delay is low enough.

To examine the net effect of IL-2 on tumor dynamics we consider the equilibrium without the presence of BCG. We use physiologically acceptable doses of 200,000–250,000 IU per instillation ($8\times {10}^5-7.7\times {10}^6$ units per treatment) [9,24].

Example 5.

By $b=0$, $i_2=5\times {10}^6$, ${\mu }_T=0.034$ and the same values of all other parameters the following two equilibria are unstable:

$$\begin{array}{cc}\hfill E_2=& (A^{*},B^{*},A_B^{*},A_T^{*},E_B^{*},E_T^{*},I_2^{*},T_i^{*},T_u^{*},F_{\beta }^{*})\hfill \\ \hfill =& (123680,0,0,0,0,0,434783,0,0,0),\hfill \\ \hfill E_3=& (A^{*},B^{*},A_B^{*},A_T^{*},E_B^{*},E_T^{*},I_2^{*},T_i^{*},T_u^{*},F_{\beta }^{*})\hfill \\ \hfill =& (123680,0,0,0,0,0,434783,0,{10}^{11},82973).\hfill \end{array}$$

(14)

Remark 4.

If the matrix H is the Hurwitz matrix then $det\left(H\right)>0$. For the equilibria $E_2$ and $E_3$ we have $E_B=0$. From (10) it follows that $a_{18}=a_{48}=a_{58}=a_{88}=0$, i.e., all elements of the eighth column of the matrix H are zeros and therefore $det\left(H\right)=0$. It means that the matrix H is not the Hurwitz matrix and the equilibrium $E_2$ and $E_3$ are unstable for all values of the parameters.

6. 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 (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 (Examples 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 (2) 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.