Stability Analysis of Delayed Immune Response BCG Infection in Bladder Cancer Treatment Model by Stochastic Perturbations

We present a revised mathematical model of the immune response to Bacillus Calmette-Guérin (BCG) treatment of bladder cancer, optimized according to biological and clinical data accumulated during the last decade. The improved model accounts for cytotoxic T lymphocyte differentiation as an integral element of the delayed immune response, as well as the logistic growth terms for cancer cell proliferation. Three equilibria are demonstrated for the proposed model, which is assumed to be influenced by white noise stochastic perturbations that are directly proportional to the system state deviation from an equilibrium. Stability conditions for all equilibria are analyzed using the Kolmanovskii-Shaikhet general method of Lyapunov functionals construction.


Introduction
Bladder cancer (BC) is 7th most common cancer (the 4th most common for men) with approximately 356,000 new cases each year and more than 145,000 deaths per year. The highest incidence occurs in industrialized and developed areas such as Europe, North America, and Australia (Jemal et al., [1]). Tobacco smoking is the main BC risk factor, accounting for at least 50% of BC cases. Roughly 10% of all BC cases have been related to occupational exposure to chemicals and dye, mostly in industrial areas processing paint, metal, and petroleum products (Bunimovich-Mendrazitsky et al. [2]).
The treatment of the BC has improved during last 40 years due to development of high definition of cystoscopy, newly technology in the bladder drugs instillation. However, the prognosis of advanced bladder cancer has not improved during the last years (Alexandroff et al. [3]). The high rates of recurrence, invasive surveillance strategies, and high treatment costs combine to make bladder cancer the single most expensive cancer in both England and the United States (Eylert et al. [4]).
BC is most frequently treated with intravesical instillations of an adjuvant immunotherapy with the Bacillus Calmette-Guérin (BCG) bacteria. BCG immunotherapy, originally established by Morales et al. [5], is administered after surgical removal of the tumor at a point where no visual lesions or morphologically evident malignant cells present in random biopsies (Brandau and Suttman, [6]). BCG immunotherapy has proven its superiority over chemotherapy in reducing tumor recurrence rates for patients with high grade or high risk nonmuscle invasive BC. While Lamm et al. [7] found BCG to even reduce disease progression, there is a need to understand why the standard BCG treatment protocol is not effective for nonresponding or relapsing patients. The BCG treatment protocol remains to be optimized specifically for those patients who do not achieve remission from treatment with the standard regimen.
In last three decades, it has been commonly accepted that a qualitative understanding of dynamic cancer treatment requires a mathematical framework in which the essential features of this complex process are represented (Byrne,[8]). The model proposed in [9] was the first mathematical model to describe tumor-immune system interactions in the bladder as a result of continuous BCG therapy. Bunimovich-Mendrazitsky et al. [9][10][11][12] have modeled the use of BCG in noninvasive bladder cancer, identifying stability points of the mathematical system to ensure durability of the simulated results, and found a tolerable bacteria threshold and effective treatment regimens to minimize undesirable side effects. A system of ordinary differential equations (ODE) was used for effective description of BCG treatment dynamics. In this manuscript, we present an improved BCG model based upon that of B-M et al. [9] which describes the tumor-immune system interactions in the bladder in response to BCG therapy, updated according to newly published biological and clinical data. Three equilibria for the optimized model are demonstrated.
One of the main problems encountered in mathematical models described by differential equations is that of their stability. In this work, equilibria stability is analyzed using the Kolmanovskii-Shaikhet general method of Lyapunov functionals construction [13][14][15][16][17] and the method of linear matrix inequalities (LMIs). Equilibria stability in probability of a system of nonlinear stochastic differential equations possessing the order of nonlinearity higher than one can be reduced to investigation of asymptotic mean square stability of the linear part of the considered nonlinear system. The proposed method of stability investigation is based upon a preliminary approximation in which the nonlinear system is centralized and linearized around an equilibrium point, and the zero solution of the obtained linear system is investigated for asymptotic mean square stability. Obtained asymptotic mean square stability conditions of the linear system zero solution at the same time are sufficient conditions for stability in probability of a corresponding equilibrium of the initial nonlinear system. This method can be applied to a system of arbitrary nonlinear differential equations with the order of nonlinearity higher than one.
The manuscript is organized as follows: in Section 2, we give general information about the BC and the biological background on which our model is based. The mathematical derivation of the improved model is described in Section 3. In Section 4 we discuss the existence of three equilibria. Section 5 is devoted to stochastic perturbations, centering, and linearization of the considered system around equilibria including description of an appropriate linear system for each from three equilibria. The conditions of stability for all three equilibria and some numerical examples are discussed in Section 6. Some auxiliary information is presented in Appendix. Conclusions from a number of clinical scenarios leading to successful treatment or to progression of cancer and some future projects are outlined in Section 7.

Biological Framework for Optimization of the Model
BCG is thought to encourage tumor elimination by attachment of the BCG to the urothelium and initiation of localized inflammation, which attracts innate immune cells that in turn draw cytotoxic T-lymphocytes (CTLs) and natural killer (NK) cells which attack the tumor cells [6,18]. Bacteria can be engulfed by macrophages and dendritic cells (DCs), two types of antigen-presenting cells (APCs). Bacteria may also infect occasional residual cancer cells, which present bacterial protein fragments, attracting APCs that will ingest the infected host. Once a BCG-infected tumor cell has been ingested by an APC, its tumor antigens are now presented by the APCs, which further stimulate CTL proliferation, together with the bacterial infection itself. These effects upon CTL proliferation were not taken into account in the original BCG model [9]. Recent studies have described the BCG-immune system interactions: Biot et al. [18], Kikamura and Tsukamoto [19], and Bunimovich-Mendrazitsky et al. [12], which describe two distinct CTL populations capable of destroying tumor cells, either via tumor associated Ag mechanism or via bacterial associated Ag. Another important factor not accounted for in [9] is the rate of maturation and activation of CTL cells [20,21]. A limiting factor of CTL development is the time required to identify bacterial and tumor antigens and to mobilize an effective immune response [22]. While the modeling of each relevant stage involved in this process is not possible within the current ODE system, a simplified model representing tumor-immune system interplay will provide insight to this fundamental yet complex interaction.

Description of the Model
Our model describes the parameters governing the efficacy of BCG treatment for bladder cancer. Based in part upon previous study [9], we have further optimized the model to account for the delay inherent to the immune response to BCG injection, taking into account the latest experimental findings regarding BCG immunotherapy of bladder cancer.
We model the stage in pathogenesis wherein no metastases have yet occurred, such that the entire dynamics of the system take place within the bladder lumen. Therefore, a onesite mathematical model is sufficient.
The BCG treatment model is composed of four nonlinear ODEs to characterize the interactions between the four different biological components, with the local quantity of each noted as follows: (i) BCG bacteria within the bladder as .
(ii) Effector T-lymphocytes, principally CTLs that react to BCG and tumor antigens as .
(iii) Tumor cells infected with BCG as .
(iv) Tumor cells not infected with BCG as .
To summarize briefly, we assume BCG to be introduced into the bladder at a constant rate . The free BCG binds to tumor cells, infecting them at rate 2 , with 1 denoting the natural death rate of BCG. Tumor cells are tracked by a continuous variable, due to their large number, and the homogeneous nature of the tumor population. Effector cells ( ) target and destroy infected and noninfected tumor cells ( and ) at a rate 3 and take up BCG at a rate 1 . Activation of the immune response is dependent upon: (1) the encounter between immune cells and BCG, controlled by parameter 4 (2) the encounter between immune cells and tumor cells, controlled by parameter .
The rate of inactivation of cells via encounter with is given by 5 . Natural logistic tumor growth is characterized by a maximal growth rate coefficient and is limited by 3 the maximum number of tumor cells . Finally, we denote by 2 the natural death rate of effector cells. The equations describing the interactions between these four variables are given in the following system: (1) In system (1)

Stochastic Perturbations, Centering, and Linearization
Let us assume that system (1) is exposed to stochastic perturbations that are of the white noise type and are directly proportional to the deviation of the system state ( ( ), ( ), ( ), ( )) from the point of equilibrium ( * , * , * , * ) and influence (( ),( ),̇( ),̇( )) immediately. In this way, system (1) takes the forṁ where are constants and ( ), = 1, 2, 3, 4, are mutually independent standard Wiener processes. System (8) is a system of stochastic differential equations [17,23]. Stochastic perturbations of a such type were first proposed in [24] and successfully used later in other researches (see [16,17] and references therein). An important feature of this type of perturbations is that the equilibrium of the deterministic system is also a solution of a system with stochastic perturbations. Let ( * , * , * , * ) be one from the equilibriums of system (8). Using (2) Similarly for the second equilibriuṁ and for the third onė It is clear that stability of the zero solution of system (9) ((10), (11)) is equivalent to stability of the first (the second, the third) equilibrium of system (8).

Remark 2.
To obtain sufficient conditions of stability in probability for nonlinear system with the order of nonlinearity Computational and Mathematical Methods in Medicine 5 higher than one it is enough [17] to get sufficient conditions of asymptotic mean square stability for the linear part of the considered nonlinear system.

Stability Investigation
For equations with delay there are two types of stability conditions: delay independent conditions and delay dependent conditions. Following the assumption that in the considered system (1) the delay is decreasing in infinity to zero, below we consider stability conditions that are delay independent. Nevertheless, the research method used here can be used also to obtain delay dependent stability conditions [17].
To investigate stability of the linear stochastic delay differential equations (12) consider the conditions for a matrix = ‖ ‖ to be the Hurwitz matrix.
Definition 3. The trace of the −th order of a × -matrix is defined as follows: Here, in particular, 1 = Tr( ), = det( ), and where is the algebraic complement of the diagonal element of the matrix .
Remark 6. Via Schur complement (see Appendix) instead of the nonlinear Riccati type inequalities (19) one can use, respectively, the conditions with the following linear matrix inequalities (LMIs): 6

Computational and Mathematical Methods in Medicine
Remark 7. For conditions (19) or (20) the matrix has be the Hurwitz matrix; i.e., the conditions of Lemma 4 must be hold. (18) for the matrix given in (13)

Stability of the First Equilibrium. Conditions
It is easy to see that if 2 * 1 > , i.e., 2 > 1 , then conditions (17) hold; i.e., the matrix given in (21) is the Hurwitz matrix.
Using [17] we obtain that the inequalities are sufficient conditions for asymptotic mean square stability of the zero solution of system (12), (21). From conditions (22) via Remark 2 the following statement follows.

Stability of the Second Equilibrium.
For positivity of the second equilibrium it is necessary to suppose that 4 > 1 2 . Calculating for system (12), (14) , = 1, 2, 3, 4, via Definition 3, we obtain that conditions (18) hold: Note however that two last equations of this system can be considered separately. So, if or then the zero solution of two last equations of system (12), (14) is asymptotically mean square stable [17]. From this it follows that for stability investigation of the second equilibrium it is enough to investigate asymptotic mean square stability of system (12) (12), (27) we obtain that in deterministic case ( 1 = 2 = 0) the zero solution of the system is asymptotically stable for ≥ 3.417. In the stochastic case for 1 = 0.45, 2 = 0.31 the zero solution of system (12), (27) is asymptotically mean square stable for ≥ 3.896. It means, for example, that for = 3.896 the equilibrium ( * 2 , * 2 , * 2 , * 2 ) = (3.4167, 0.1122, 0, 0) is stable in probability.
Remark 11. Note that via (23), (26), for stability of the first equilibrium must be 1 / 2 < < 1 2 / 4 and for stability of the second equilibrium must be > 1 2 / 4 . . Checking LMIs (20) we can conclude instability of this equilibrium even in the deterministic case. Similar calculations for different < 1.7466 (via Remark 1) show that the third equilibrium is unstable.

Conclusion
In this work we present the improved model of BCG immunotherapy in superficial bladder cancer. This study investigates stability of new treatment model with constant instillations of BCG under stochastic perturbations. The model demonstrates several stable states which depend on biologically related parameters and initial conditions.
The innovation of our study is the involvement of CTL cells, specific for the tumor antigen due to BCG infection. Adding BCG to the tumor-immune interaction may increase the immune response, which most benefit from the addition of the antigen. These effector cells capture tumor cells after time delay which has been used for maturation of these effector cells. The entire reaction can take place only with the presence of BCG. We added new terms in two equations (the second and the forth) in system (1). The new terms stand for the killing of uninfected tumor cells by effector cells.
It is shown that the considered system has three equilibria describing the different states of the patient. Stability of these states under stochastic perturbations which are directly proportional to the deviation of the patient's current state from the equilibrium state is investigated. New sufficient conditions for stability in probability of two equilibria and instability of the third equilibrium were obtained using the theoretical method of Lyapunov functionals and the numerical method of linear matrix inequalities (LMIs).
In the first equilibrium we obtain weak immune response because * 1 = 0. In addition BCG-infected tumor cells were remained ( * 1 > 0); that indicates cancer will not obligatorily be eradicated, although this equilibrium is stable. In the second equilibrium we receive * 2 = * 2 = 0; that means the cancer could be eradicated after BCG immunotherapy with dose > ( 2 / 4 )( 1 + ( 1 / 3 ) ). This equilibrium has very strong immune response * 2 > 0 and it is stable. The third equilibrium is unstable. The dynamic behavior of the system investigated from the point of view of local stability and a detailed analysis on stability of equilibrium was examined. As proposed in the paper research methods allow continuing and specifying investigation of the considered model and getting its new useful properties.
However, we have to note that this study has following limitations: (1) This model takes into account 3-biological processes (tumor, immune system, and BCG interactions) only, which are captured by 4 differential equations. (2) We examine continued BCG treatment with a constant rate , although in a real life pulsed therapy and variable rate can be used. (3) The duration of BCG treatment is not time limited in our model.
By capturing the main parameters of the BCG, maximum tumor size, tumor growth rate, and immune response parameters, we create a silicone model and develop an algorithm for eliminating cancer. 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 considered not to respond.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Disclosure
Note that the information about the results described above was presented only in the form of a short abstract on the conference [25].