Mathematical analysis of tumor-free equilibrium in BCG treatment with effective IL-2 infusion for bladder cancer model

1.   Introduction

Mathematical models in biology and medicine are a powerful tool for physicians and biologists, as can be seen from the proliferation of both academic and clinical projects on this topic over the past few decades ^{[1,2,3,4,5,6,7,8]}.

Cancer is a complex clinical problem because each type of cancer has unique properties and biological dynamics. Therefore, the use of mathematical modeling in support of clinical research, as well as decision making in oncology, encompasses many problems and research questions ^[9].

Bladder cancer (BC) is a widespread disease: the 7th most common cancer (the 4th most common for men) with approximately 549,000 new cases each year and more than 200,000 deaths per year. The highest incidence occurs in industrialized and developed areas such as Europe, North America, and Australia ^[10]. The primary cause of about half of BC cases is occupational exposure to chemicals in industrial areas, processing paints, metals, dyes, and petroleum products. Tobacco smoking and environmental carcinogens are other risk factors for BC ^[11].

Historically, the treatment of BC has hardly advanced during the last 50 years with the use of Bacillus Calmette-Guerin (BCG) as a treatment suggested by Morales et al. (1976) ^[12]. BCG is an attenuated nonpathogenic strain of Mycobacterium bovis that was originally used as a vaccine against tuberculosis. In this treatment, bacterial instillations are introduced into the bladder with a lighted tube (catheter) that is inserted through the urethra ^[13]. BCG-immunotherapy has been the standard therapy and has been proven to reduce both recurrence and progression of BC following primary resection (TUR). However, a significant proportion of patients do not achieve remission after BCG. Therefore, the BCG treatment protocol has yet to be specifically optimized for those patients. The cytokine-based therapies, such as IL-12, IFN-$\alpha$, and IL-2, have been developed to improve BCG therapy in ^[14,15,16].

Over the past decades, there have been several attempts to develop models of BCG treatment using ODE and Partial Differential Equations (PDEs) to find the optimal treatment protocol for BC patients ^[17,18,19]. In silico models have been created and strategies have been developed to eliminate bladder cancer by distinguishing between patients based on their immunological parameters ^[20]. One of the possible improvements in treatment is to test patients for different immunological status to adjust the treatment protocol in order to increase efficiency in problem cases. This is the goal of our present work: identification of the dose of IL-2 and its effect on the outcome of treatment. Despite its benefits, IL-2 is notorious for its side effects. Hence, cautious dosing is required and only a low dosage is applicable. The first study of BCG+IL-2 combination therapy, ^[21], did not find a significant effect of IL-2 on tumor clearance. The next BCG+IL-2 combination therapy study, ^[22], found an effect of IL-2 on tumor clearance, but the dose of IL-2 was constant. In this work, presenting IL-2 treatment in integral form allows us to change the treatment dose and treatment interval for each patient.

One of the main problems in mathematical models described by differential equations is to find stability conditions for equilibria points. In previous works, we analyzed the stability using the Lyapunov method ^[17], the Kolmanovsky-Shaikhet method ^[24], and the linear matrix inequalities (LMI) method ^[23]. In this study, the effect of IL-2 on the immune response during BCG treatment is observed, and the stability condition of the tumor-free equilibrium is analyzed analytically via the Cauchy matrix ^[25,26,27].

This manuscript is organized as follows: In Section 2, we introduce a BCG treatment model with IL-2 distributed feedback control. In Section 3, the stability condition for the model with distributed control on the immune system in the form of an integral term is obtained. In Section 4, we build the Cauchy matrix of the linearized system for further research of our model. In Section 5, we examine the individual reaction of the patient to the IL-2 injections. For this aim we obtain the stability condition of the model with an uncertain coefficient, using the Cauchy matrix. In Section 6, a system with a delay in the upper and lower limits of distributed control in the form of an integral term is analyzed. In the last section, we discuss the results including the advantages and limitations of the proposed model, including the personalizing of treatment according to treatment protocol, in addition to potential future directions.

2.   BCG + IL-2 combined treatment model for bladder cancer

We based our work on a BCG treatment model ^[23] for BC. The model considers four main biological variables to describe the treatment process:

- $B$ - BCG bacteria introduced into the bladder,

- $E$ - effector $T$ cells - lymphocytes that respond to both BCG antigens and cancer antigens,

- $T_{i}$ - BC cells infected with BCG,

- $T_{u}$ - BC cells not infected with BCG.

During treatment, BCG with rate $b(t)$ was entered into the bladder by a catheter and invade the antigen-presenting cells (APCs) (macrophages or dendritic cells (DCs)) $E$ at a rate $p_{1}$. BCG is internalized with APC via endocytosis and immune response activation starts. At the same time, BCG infect cancer cells $T_{u}$ that remain after surgery cells at rate $p_{2}$. The tumor cells $T_{u}$ infected by BCG are referred to as $T_{i}$.

Infected tumor cells $T_{i}$ stimulate recruitment of cytotoxic effector cells (CTL) $E$ from the bone marrow ^[28] at a rate $\alpha$. APCs induce CTL according to BCG antigen. Hence the encounter between $E$ and $B$ is controlled by the parameter $p_{4}$. After the tumor infected cell $T_{i}$ has been ingested, the tumor antigens are presented to the APC. APCs induce CTL according to their tumor antigen. In other words, there is a two-stage elimination of tumor cells: first, destruction of infected tumor cells by the action of effector CTLs infected with BCG, and second, destruction of uninfected tumor cells by CTLs with tumor antigens. Hence, $E$ cells target and destroy infected and non-infected tumor cells $T_{i}$ and $T_{u}$ at a rate $p_{3}$. We define the death rate of infected tumor cells by $\mu_{3}$. The natural death rate of BCG is denoted by $\mu _{1}$. All parameter values described in this section are in Appendix B, Table 1.

The tumor growth rate in the bladder is characterized by a coefficient $r$. $T_u$ is limited by the tumor cell carrying capacity in the bladder, $K$, where $T_u \leq K$. The scheme of "immune system$<$—BCG—$>$tumor" response is presented in Figure 1.

The ODE system describing biological interactions between $B$, $E$, $T_{i}$, $T_{u}$ is presented in (2.1):

where the model begins with non-negative initial conditions:

We introduce the cytokine IL-2 into the system (2.1), into the effector cell equation inducing the proliferation of cytotoxic immune cells, as a distributed feedback control in an integral form and get the following system of IDE:

where $G(t) = \int_{0}^{t}e^{-\gamma \left(t-s\right) }E\left(s\right) ds$, $\gamma$ and $A$ are constant values, characterizing the IL-2 dose as indicated by the dotted line in Figure 1.

3.   Stability analysis of system (2.2)

The system of IDE (2.2) can be reduced to the following ODE system:

We assume that all coefficients and parameters of the system (3.1) are positive and

The equilibria of the model are found by setting all derivatives to zero and solving for $B^{\ast }, E^{\ast }, T_{i}^{\ast }, T_{u}^{\ast }$ and $G^{\ast }$ their equilibrium values. Equations have multiple equilibria, but we need only focus on the nonnegative equilibria assuming all initial conditions are positive.

In this work, we want to focus on the tumor-free equilibrium: $B^{\ast } = \frac{b}{\mu _{1}},$ $\ \ E^{\ast } = T_{i}^{\ast } = T_{u}^{\ast } = G^{\ast } = 0.$

Linearization of system (3.1) in the neighborhood of the tumor free equilibrium, where $x_{1} = B-B^{\ast },$ $x_{2} = E,$ $x_{3} = T_{i},$ $x_{4} = T_{u},$ $x_{5} = G$ is

In this case, the Jacobian is

where $H = \frac{bp_{4}}{\mu _{1}}-\mu _{2} < 0$. Therefore, we obtain the following eigenvalues of matrix $J$

Theorem 3.1. Let inequalities $r < \frac{bp_{2}}{\mu _{1}}, \ \ H < 0, \ \ -\gamma H > A > 0$ be fulfilled, then the tumor-free equilibrium is exponentially stable.

Proof. It is clear that $\lambda_1, \lambda_2, \lambda_3$ are real and negative eigenvalues. $\lambda _{4}, \lambda _{5}$ are real and negative eigenvalues and $\lambda _{4}\neq \lambda _{5}$, because for $\lambda_4$, $\left(\gamma -H\right) > \sqrt{\left(\gamma +H\right) ^{2}+4A}$ $\Rightarrow$ $\left(\gamma -H\right) ^{2} > \left(\gamma +H\right) ^{2}+4A$ $\Rightarrow$ $-\gamma H > A$ and for $\lambda_5$, $\left(\gamma +H\right) ^{2}+4A > 0,$ $-\left(\gamma -H\right) < 0$.

4.   The differences between convergent speed of systems (2.1) and current model (3.1)

Linearizing system (3.1) in the neighborhood of the tumor-free equilibrium: $B^{\ast } = \frac{b}{\mu _{1}},$ $\ \ E^{\ast } = T_{i}^{\ast } = T_{u}^{\ast } = G^{\ast } = 0$, we obtain (3.2), and system (2.1) we obtain

where $x_{1} = B-B^{\ast },$ $x_{2} = E,$ $x_{3} = T_{i},$ $x_{4} = T_{u}$.

The corresponding matrix of coefficients of system (4.1) are

Eigenvalues of matrix (4.2) are

Theorem 4.1. If all coefficients of system (2.1) are positive and inequalities $r < \frac{b p_2}{\mu_1}$, $\frac{bp_{4}}{\mu _{1}} < \mu _{2}$ are fulfilled, then system (2.1) is exponentially stable in tumor-free equilibrium.

Proof. The proof of this theorem follows from the negativity of eigenvalues.

Let us denote spectral radius of base system (3.2) and of moderated system (4.1) correspondingly by

Theorem 4.2. If all coefficients of system (3.1) are positive, and the inequalities

are fulfilled, then system (3.1) is exponentially stable in tumor-free equilibrium and $\tilde{\rho} < \rho$.

Proof. It is clear that $\lambda_1 = \tilde{\lambda_1} = -\mu_1, \ \ \lambda_2 = \tilde{\lambda_2} = -\mu_3, \ \ \lambda_3 = \tilde{\lambda_3} = r-\frac{bp_{2}}{\mu _{1}}$. We have to show that $|\lambda_5| > |\tilde{\lambda_4}|$, because $|\lambda_5| > |\lambda_4|$.

But $|\lambda_5| = \frac{\left(\gamma -H\right) +\sqrt{\left(\gamma +H\right) ^{2}+4A}}{2}$ and $|\tilde{\lambda_4}| = -H$ and we have to show that

$\frac{\left(\gamma -H\right) +\sqrt{\left(\gamma +H\right)^{2}+4A}}{2} > -H,$ i.e. $\sqrt{\left(\gamma +H\right)^{2}+4A} > -H-\gamma$.

If $-H-\gamma \leq 0$ it is true.

If $-H-\gamma > 0$, then $\left(\gamma +H\right)^{2}+4A > (H+\gamma)^2$ is also true. Hence

5.   Constructing the Cauchy matrix of system (3.2)

Let us denote $R = \sqrt{(\gamma+H)^2+4A}$. Solving the homogeneous system (3.2) (see Appendix A), we obtain

where

Let us denote $X(t) = \left\{ x_{1}(t), x_{2}(t), x_{3}(t), x_{4}(t), x_{5}(t)\right\}$

Taking into account that $C (s, s) = I$ we obtain the cases below.

1) If $X(0) = \left\{ 1, 0, 0, 0, 0\right\}$, we obtain

The first column of the Cauchy matrix is the following

2) If $X(0) = \left\{ 0, 1, 0, 0, 0\right\}$, we obtain

The second column of the Cauchy matrix is the following

3) If $X(0) = \left\{ 0, 0, 1, 0, 0\right\}$, we obtain

The third column of the Cauchy matrix is the following

4) If $X(0) = \left\{ 0, 0, 0, 1, 0\right\}$, we obtain

The fourth column of the Cauchy matrix is the following

5) If $X(0) = \left\{ 0, 0, 0, 0, 1\right\}$, we obtain

The fifth column of the Cauchy matrix is the following

6.   System with uncertain coefficient

Consider the following system of IDE

where $a\left(t\right)$ is the individual addition to the constant dose of IL-2 injections for BC patient. Using the Cauchy matrix we obtain the conditions under which the system (6.1) will be exponentially stable in tumor-free equilibrium. The system of IDE (6.1) can be reduced to the following ODE system:

where

Linearizing system (6.2) we obtain

Consider the system

where

Denote

The general solution of the system

can be represented in the following form (see, for example, ^[26])

Without loss of generality, $X(0) = col\left\{ 0, 0, 0, 0, 0\right\}.$ Substituting (6.6) into (6.4) we obtain

which can be rewritten in the operator form

where the operator $\Omega :L_{\infty }^{5}\rightarrow L_{\infty }^{5}$, where $L_{\infty }^{5}$ is the space of 5 vector-functions with essentially bounded components, is defined by

where $C\left(t, s\right)$ is the Cauchy matrix of system (3.2).

Let $\|\Omega\| = sup_{1 \leq i \leq 5} \sum_{j = 1}^5 |w_{ij}|$ be the norm of the operator $\Omega$.

Denote $P_j = ess \sup_{t \geq 0} \int_0^t \sum_{i = 1}^5 |(L_a(t)C(t, s))_{ij}|ds$ and $a^{\ast } = ess\sup_{t \geq 0}\left\vert a\left(t\right) \right\vert$, we obtain the estimates:

It is clear that $P_{j}\leq P_{j}^{\ast }$ for $1\leq j\leq 5.$

Theorem 6.1. Let inequalities $r < \frac{bp_{2}}{\mu _{1}}, \ \ H < 0, \ \ -\gamma H > A$ be fulfilled, and the inequality $\max_{1\leq j\leq 5}\left\{ P_{j}^{\ast }\right\} < 1$ be true, then the tumor-free equilibrium is exponentially stable.

Proof. The proof follows from estimates of $P_j, \ \ 1 \leq j \leq 5.$

Let us denote $Q = \frac{1}{a^*}\max_{1\leq j\leq 5}\left\{ P_{j}^{\ast}\right\}$, then a system with an uncertain coefficient is still exponentially stable in tumor-free equilibrium, if for a maximum value of function $a(t)$ the following inequality is true: $a^* < \frac{1}{Q}$.

7.   Exponential stability of system (2.2) with delay in upper and lower limit of control function

Consider system (2.2), where

where $\tau_2(t) < \tau_1(t), \ \tau_1(t) < t, \ \tau_2(t) < t$.

Let us denote

We can write (7.1) in the following form

Reducing IDE system (2.2), with $G(t)$ defined by (7.1) to the ODE system, we obtain

Linearizing system (7.4) in the neighborhood of the tumor free equilibrium, we obtain

where $x_{1} = B-B^{\ast },$ $x_{2} = E,$ $x_{3} = T_{i},$ $x_{4} = T_{u},$ $x_{5} = \tilde{G}$.

Let us denote $\tau_1^* = {\rm{ess}}\; \inf_{t \geq 0} |\tau_1(t)|$ and $\tau_2^* = {\rm{ess}}\; \inf_{t \geq 0} |\tau_2(t)|$ and

Theorem 7.1. Let inequalities $r < \frac{bp_{2}}{\mu _{1}}, \ \ H < 0, \ \ -\gamma H > A$ be fulfilled and $Q < 1$, then the tumor-free equilibrium is exponentially stable.

Proof. We can write system (3.2) in the following form

where

It is known that the general solution of the system

can be written in the following form

where $C(t, s)$ is a Cauchy matrix of system (3.2). We can rewrite (7.5) in the following form

Without loss of generality, we can assume that $X(0) = 0$. Substituting $X(t) = \int_{0}^{t} C(t, s)W(s) ds$ into system (7.9), we obtain

and for $W(t)$ the following system (see (7.7))

where $W(t) = col\{w_1(t), w_2(t), w_3(t), w_4(t), w_5(t)\}$, and $\Omega:L_{\infty }^{5}\rightarrow L_{\infty }^{5}$ is the operator that is defined by

where

Estimating the norm of operator $\Omega$ we obtain the assertion of this theorem.

Example 1. We use the parameter values from Table 1:

From the condition $-\gamma H > A$, we take $A = 10^5$, and we obtain the following condition:

This inequality are fulfilled, for example, in the point $\tau_1^* = 10^{-9}, \ \tau_2^* = 0.5$.

8.   Discussion and Conclusions

BCG therapy has become the standard BC treatment since 1976, but a significant proportion of patients still experience cancer progression and/or recurrence despite this therapy. One of the attempts to improve BCG efficacy is a combined therapy of BCG and IL-2. IL-2 is a good candidate, despite its 'toxic' side effects. Hence the dose of IL-2 used for therapy should take into account the individual effect of every patient in order to keep the proper balance between the activities of immune-killing cells and the eradication of BC cells.

Mathematical modeling is a useful tool in oncology, allowing us to explore possible treatments to find personalized treatment protocols for each patient. In this research, we added IL-2 therapy to an effector cell equation for better treatment outcomes as a distributed feedback control in an integral form. We obtained stability conditions in the neighborhood of a tumor-free point. The model (2.2) considers changes in the population of BC cells and immune cells under the influence of external factors, such as BCG and IL-2 immunotherapy. The ability to plan and predict by calculating the modulated dose of IL-2 treatment may benefit patients who cannot take standard treatment due to its severe side effects, as well as patients who were previously considered non-responsive to treatment.

A clinically significant feature of the model is the non-trivial dependence of treatment success on the dose of BCG and IL-2. In particular, the following requirements can be distinguished, clearly dependent on the appropriate treatment regimen and tumor growth obtained for the exponential stability condition of the BCG + IL-2 model (2.2) formulated in Theorem 3.1:

1. Treatment rate b of BCG: $\frac{r\mu _{1}}{p_{2}} < b < \frac{\mu _{1}\mu _{2}}{p_{4}}$;

2. Treatment rate A of IL-2: $\gamma \cdot (\mu _{2}-\frac{bp_{4}}{\mu _{1}}) > A$.

If one of these conditions is not met, then this corresponds to a clinical failure of the BCG + IL-2 treatment.

We compared the effectiveness of the treatment presented in the model (2.1) from the work ^[23] with the current model (3.1) using the spectral radius theory. We proved in Theorem 4.2 that therapy based on the current model (3.1) is much more effective because the spectral radius of the system (3.1) is greater (4.5) than the spectral radius of the system (2.1).

To control IL-2 in each injection, we defined the function a(t), which shows the individual change in the constant dose A for each patient with breast cancer in the system (6.1). To obtain the maximum variations in IL-2 doses, we constructed the Cauchy matrix of system (3.2) and, using the representation of the solution of system (5.3) in operator form, obtained the conditions under which system (6.1) will be exponentially stable in a tumor-free equilibrium.

To get a more flexible definition of the IL-2 treatment period, we introduce in the BCG + IL-2 model (2.2) the delay functions in the upper and lower limit of control function, which presents IL-2 immunotherapy defined in (7.1). These delays indicate the changes in the start and in the end of a patient's IL-2 treatment period as shown in Example 1. However, these results were obtained on a simplified version of the biological system with the participation of the immune system and personalized combination therapy for BC patients. In particular, there is a delay between the introduction of IL-2 into the bladder and an increase in the number of effector cells, which is neglected in order to avoid the delayed ODE system. This and other assumptions may lead to slightly different results.

The proposed model provides a mathematical framework for investigating the impact of different BCG+IL-2 treatment protocols, but its ability to predict outcomes at the patient level is limited. However, our model can serve as a baseline that can be extended to study more detailed dynamics in a broader context. Complementary therapies, such as various immunotherapies, including vaccines or modified cells, may improve both targeted and symptomatic outcomes in BC patients. New and effective therapies are required to increase the specificity and strength of the immune system against cancer^[31]. Our model is just such a new contribution to this area.

Conflict of interest

The authors declare no conflict of interest.

Appendix A

Solution of the fourth equation of (3.2), is

Solving the third equation of (3.2)

we obtain solution

Combining second and fifth equation of (3.2)

we obtain solutions

where

and

Solving the first equation of (3.2)

we obtain solution

where

Let us denote Jacobian of the system (3.1) in the tumor-free equilibrium by $J$. Linearization of system (3.1) in the neighborhood of the tumor free equilibrium we obtain

where where $x_{1} = B-B^{\ast }, \ \ x_{2} = E, \ \ x_{3} = T_{i}, \ \ x_{4} = T_{u}, \ \ x_{5} = G$.

Appendix B