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.

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

2. 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.

3. 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 one-site 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:

BCG bacteria within the bladder as B.

Effector T-lymphocytes, principally CTLs that react to BCG and tumor antigens as E.

Tumor cells infected with BCG as T .

Tumor cells not infected with BCG as T .

To summarize briefly, we assume BCG to be introduced into the bladder at a constant rate b. The free BCG binds to tumor cells, infecting them at rate p 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 (E) target and destroy infected and noninfected tumor cells (T and T ) at a rate p 3 and take up BCG at a rate p 1. Activation of the immune response is dependent upon:

the encounter between immune cells and BCG, controlled by parameter p 4

the encounter between immune cells and tumor cells, controlled by parameter α.

The rate of inactivation of E cells via encounter with T is given by p 5. Natural logistic tumor growth is characterized by a maximal growth rate coefficient r and is limited by the maximum number of tumor cells K. 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:In system (1) in the second equation (effector cells ) BCG induces the recruitment of the effectors with rate p 4 B(t)E(t − τ(t)), while the tumor induces the recruitment of the effectors at a linear rate αT (t)E(t − τ(t)), with delay function τ(t) = τ 0 e −, λ > 0. Correspondingly, eradication of infected T (in the third equation) and uninfected T (in the fourth equation) tumor cells by effector cells will produce with the same delay τ(t) function.

τ(t) is a time-varying function, representing the immune response delay in response to treatment, and expressing the number of effector cells in the cancer area. The delay is measured with reference to beginning of BCG treatment (t = 0), with a maximum delay of approximately 10 days. The influence of BCG tends towards zero over time, such that the delay function is similar to fading phenomenon: τ(t) = τ 0 e −.

4. Three Equilibria

The equilibria of system (1) are defined by the conditions , , , and that is equivalent to the system of four algebraic equationswith three solutions: whereand T 3 is defined by the algebraic equationwith B 3 , R , and F defined in (5) and (6) (see Appendix).

Remark 1 .

Putting x = T 3 /K, β = μ 1/p 2, and γ = (1 + p 1 p 5/p 3 p 4)K, via (6), and a positivity of R for the third equilibrium one can obtain the estimation for b: b < F H = r(1 − x)(β + γx) ≤ r(β + γ)2/4γ.

5. 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 (B(t), E(t), T (t), T (t)) from the point of equilibrium (B , E , T , T ) and influence immediately. In this way, system (1) takes the formwhere σ are constants and w (t), i = 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 (B , E , T , T ) be one from the equilibriums of system (8). Using (2) substitutes new variables y 1(t) = B(t) − B , y 2(t) = E(t) − E , y 3(t) = T (t) − T , and y 4(t) = T (t) − T , into system (8) for each equilibrium.

For the first equilibrium we obtain the nonlinear systemSimilarly for the second equilibriumand for the third oneIt 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 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.

Following Remark 2, represent the linear parts of systems (9) - (11) in the formwhere z = (z 1, z 2, z 3, z 4)′, w = (w 1, w 2, w 3, w 4)′, C(z) = diag{σ 1 z 1, σ 2 z 2, σ 3 z 3, σ 4 z 4}, and the matrices A, D are, respectively, for the first equilibriumfor the second equilibriumand for the third equilibrium

6. 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 A = ‖a ‖ to be the Hurwitz matrix.

Definition 3 .

The trace of the k−th order of a n × n-matrix A is defined as follows: Here, in particular, S 1 = Tr(A), S = det(A), and S = ∑ A , where A is the algebraic complement of the diagonal element a of the matrix A.

Lemma 4 (see [17]).

A 3 × 3-matrix A is the Hurwitz matrix if and only ifA 4 × 4-matrix A is the Hurwitz matrix if and only if

Lemma 5 (see [17]).

Put Φ0(P) = A′P + PA + diag{p 11 σ 1 2, p 22 σ 2 2, p 33 σ 3 2, p 44 σ 4 2}. If and for some positive definite matrices P and R at least one of the inequalities,holds, then the zero solution of (12) is asymptotically mean square stable.

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):

Remark 7 .

For conditions (19) or (20) the matrix A has be the Hurwitz matrix; i.e., the conditions of Lemma 4 must be hold.

6.1. Stability of the First Equilibrium

Conditions (18) for the matrix A given in (13) do not hold since S 4 = 0. But putting T 1 = 0 and rejecting the equation for we obtain system (12) for the equilibrium (B 1 , E 1 , T 1 ) = (bμ 1 −1, 0,0) with the 3 × 3 matricesHere S 1 = −(μ 1 + μ 2 + p 2 B 1 − r), S 2 = μ 1 μ 2 + (μ 1 + μ 2)(p 2 B 1 − r), and S 3 = −μ 1 μ 2(p 2 B 1 − r). It is easy to see that if p 2 B 1 > r, i.e., bp 2 > rμ 1, then conditions (17) hold; i.e., the matrix A given in (21) is the Hurwitz matrix.

Using [17] we obtain that the inequalitiesare 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.

Lemma 8 .

Ifthen the first equilibrium (B 1 , E 1 , T 1 ) = (bμ 1 −1, 0,0) of system (12) is stable in probability.

Example 9 .

Let μ 1 = 1, μ 2 = 0.41, p 2 = 0.28, p 4 = 0.12, r = 0.078. Then via (23) possible interval for b ∈ (0.279,3, 417). If b = 1 and σ 1 2 < 2, σ 2 2 < 0.58, σ 4 2 < 0.404, then the equilibrium (1,0, 0,0) is stable in probability.

6.2. Stability of the Second Equilibrium

For positivity of the second equilibrium it is necessary to suppose that p 4 b > μ 1 μ 2. Calculating for system (12), (14) S , i = 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 orthen 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) withFor this aim one can solve the nonlinear inequalities (19) or to check feasibility of the LMIs (20) via MATLAB.

Example 10 .

Let μ 1 = 1, μ 2 = 0.41, p 1 = 1.25, p 2 = 0.28, p 3 = 1.1, p 4 = 0.12, p 5 = 0.003, r = 0.078, and α = 3.7. The second equilibrium exists if b > μ 1 μ 2 p 4 −1 = 3.4167. Checking the LMIs (20) for system (12), (27) we obtain that in deterministic case (σ 1 = σ 2 = 0) the zero solution of the system is asymptotically stable for b ≥ 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 b ≥ 3.896. It means, for example, that for b = 3.896 the equilibrium (B 2 , E 2 , T 2 , T 2 ) = (3.4167,0.1122,0, 0) is stable in probability.

Checking the LMIs (20) for the initial system (12), (14) we obtain the close result: in the deterministic case the zero solution of system (12), (14) is asymptotically stable for b ≥ 3.417. In the stochastic case for σ 1 = 0.45, σ 2 = 0.31, σ 3 = 0.5, and σ 4 = 1.4 the zero solution of system (12), (14) is asymptotically mean square stable for b ≥ 3.902. Conditions (18) and (26) in this case naturally hold. For example, for b = 3.902 the equilibrium (B 2 , E 2 , T 2 , T 2 ) = (3.4167,0.1136,0, 0) is stable in probability.

Remark 11 .

Note that via (23), (26), for stability of the first equilibrium must be rμ 1/p 2 < b < μ 1 μ 2/p 4 and for stability of the second equilibrium must be b > μ 1 μ 2/p 4.

6.3. Stability of the Third Equilibrium

Example 12 .

Put μ 1 = 1, μ 2 = 0.41, p 1 = 1.25, p 2 = 0.28, p 3 = 1.1, p 4 = 0.12, p 5 = 0.003, α = 3.7, K = 80, r = 0.078, and b = 0.99. By given values of the parameters via (6), (7) we obtain T 3 = 11.19, F = 0.0671, G = 40.993, H = 15.0793, and R = 1247.6 and via (5) the equilibrium is (B , E , T 3 , T 3 ) = (7.4126,0.0015,13961,11.19). Checking LMIs (20) we can conclude instability of this equilibrium even in the deterministic case. Similar calculations for different b < 1.7466 (via Remark 1) show that the third equilibrium is unstable.

7. 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 E 1 = 0. In addition BCG-infected tumor cells were remained (T 1 > 0); that indicates cancer will not obligatorily be eradicated, although this equilibrium is stable. In the second equilibrium we receive T 2 = T 2 = 0; that means the cancer could be eradicated after BCG immunotherapy with dose b > (μ 2/p 4)(μ 1 + (p 1/p 3)δ). This equilibrium has very strong immune response E 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:

This model takes into account 3-biological processes (tumor, immune system, and BCG interactions) only, which are captured by 4 differential equations.

We examine continued BCG treatment with a constant rate b, although in a real life pulsed therapy and variable rate can be used.

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.