BCG and IL − 2 model for bladder cancer treatment with fast and slow dynamics based on SPVF method—stability analysis

1.   Introduction

Biological, chemical, physical and other phenomena with applications to various engineering science, are described by mathematical models that are described by a large set of complex linear or non-linear ordinary or partial differential equations. These systems of equations can be solved numerically. But one of the main drawbacks of a numerical solution of a system of equations is the loss of singular points of the system which can be critical for engineering applications. Another drawback in solving a complete set of equations is the time of computer running. Hence, there are different methods to reduce the system of the equations without losing information of the original one. These mathematical models have a number of essentially different time scales (i.e. rates of change) which correspond to sub-processes. These different time scales, in general, are not explicit for a given models, i.e., the hierarchy of the model is hidden. Once we expose the hierarchy, we can decompose the system of equations into fast and slow subsystems, and these decomposition allows one to apply different asymptotic methods. There are several asymptotic methods and numerical tools that can be applied to multi-scale systems. For example, the method of integral invariant manifold ($MIM$) ^[1,2,3,4], the iteration method of Fraser and Roussel ^[5,6,7,8], the computational singular perturbation ($CSP$) method ^[9,10], geometric singular perturbation theory ^[11,12,13], and the intrinsic low dimensional method ($ILDM$) which is a numerical method ^{[14,15,16,17,18]}. Each method has its advantages and disadvantages. The main disadvantage of all of these methods is that they are based on the fact that the hierarchy of the set of equations is given, i.e., all these model have the form of singularly perturbed system ($SPS$). According to what we have discussed above, we first need to expose the hierarchy of the given system. In this paper we introduced and applied a method called Singular Perturbed Vector Field ($SPVF$) ^[19,20] which is a new version of $ILDM$ method. Since a large number of differential equations (a mathematical model) is presented in a hidden hierarchy form, no explicit time scale of the system is known. And hence, our aim is first to expose the time-scale of a given system, and then apply a reduction method. This is the main result of the $SPVF$ method. The $SPVF$ method transfers the original system to a form of singularly perturbed system ($SPS$), i.e., a system with an explicit hierarchy of the dynamical variables of the model. Once we transfer the system to $SPS$, the new system can be treated by the very powerful machinery of the standard $SPS$ theory for model reduction and decomposition, as we have mentioned above without losing the essential dynamics of the original system.

In addition, exposing the hierarchy of the model and presenting the system of equations in $SPS$, allow one to decompose the model into fast and slow subsystems and hence enables one to investigate the stability of the reduced model, which in itself is very important in models that describe cancer in general ^[21,22,23].

In this research, we investigate the model of bladder cancer treatment that takes into account $BCG$ and $IL-2$ treatment presented in ^[24,25]. During bladder cancer therapy in the superficial phase the tumor is amenable to local excision, where small regions of cancerous tissue are surgically resected by direct inspection through the urethra, in a procedure called transurethral resection ($TUR$) ^[26]. To complete the procedure in the superficial stage, an adjuvant treatment is administered into the cavity is generally recommended to destroy any malignant cells that remain following the resection. The two complementary approaches to adjuvant treatment of superficial bladder cancer are the intravesical chemotherapy and the intravesical application of bacillus Calmette-Guerin ($BCG$) (immunotherapy).

The mathematical model ^[24,25] describes the dynamics of $BCG$, $APC$, effector and tumor cells as a result of $BCG$ and $IL-2$ instillations. During this therapy, the dose of $BCG$ and $IL-2$ instills into the bladder via a catheter at every scheduled time. These instillations are described in the model via a set of Dirac delta functions. The disadvantage of using this function is that it is impossible to apply different asymptotic methods, except for solving the model numerically. Hence, we wrote an explicit function that describes the amount of $BCG$ that is instilled into the bladder and which is depended on time. An explicit expression of function allowed us to overcome the disadvantages mentioned above as well as that we could control the dosages at different periods of $BCG$ and $IL-2$ pulsing.

We investigated the mathematical model ^[24,25] that contains ten nonlinear ordinary differential equations by applying the $SPVF$ method and exposing the hierarchy of the system explicitly. We used the eigenvectors of the system to change the coordinates and to present the model at $SPS$ form with two fast equations and eighth slow equations. Revealing the hierarchy of the system allowed us to investigate the fast sub-system of the model, which reduces the complexity of the calculations and decreases the run-time of the computer. In addition, we presented the stability analysis of the reduced model. In general, knowledge about global fast manifolds is very important in the stability analysis of a reduced model described by a slow manifold. Now, depending on the investigated dynamical regime one can use the fast manifold as a manifold equation for the reduced space that can be represented by a low-dimensional manifold in the detailed linear vector space, which is given in an explicit form and proceed with the reduction procedure as projected to the manifold.

2.   Mathematical model of bladder cancer

In this section, we present the system of the governing equations (ODE) for bladder cancer treatment from ^[24,25], with $BCG$ and $IL-2$ combination therapy. The main assumptions of the model are as follows: A dose $b_{k}$ of $BCG$ is instilled into the bladder every $t$ days. It is given in $6$-weekly intravesical instillations of $2\cdot 10^{8}$ c.f.u (colony-forming units, i.e., number of viable bacterial cells) ^[28]. After instillation, $BCG$ accumulates close to the bladder wall. Upon binding to wall cells, $BCG$ is internalized into the bladder and is processed by $APC$s ^[29] as well as the $BCG$ binds and enters into malignant tumor cells ^[26]. The number of $APC$s increase as a result of recruitment in the response to bacterial infection. The bacterial infection stimulates the $APC$s to produce inflammatory cytokines such as $IL-2$. Simultaneously, bacteria infect the occasional non-visualized residual cancer cell that remains after ^[30]. This causes the presentation of bacterial $Ag$ on the tumor cell surface, which attracts $APC$s that ingest the entire host. Once a tumor cell has been ingested, tumor $Ag$s are presented by the $APC$s. Due to the aforementioned inflammatory environment, created by the bacterial infection, $APC$s cause the $CTL$s to either mature and track bacteria $Ag$ or mature and capture tumor cells according to their $TAA$ ^[31]. This means that 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 exogenous $IL-2$ is expected to create an inflammatory environment that stimulates the maturation of either $CTL$ population and prolongs the life span of $CTL$s ^[32].

We describe the instillation of $BCG$ and $IL-2$ into the bladder by a periodic function that at first time rises very rapidly (as exponentially). This actually describes the rapid instillation of the $BCG$ and $IL-2$ into the body, and then the function decreases in the form of decay exponential, which describes the removal of the $BCG$ from the body. In order to define the function as described above in a periodic form, we used a variation of the Gamma probability density function in different intervals, then we summed all these functions.

In the papers ^[24,25], they used the well known Dirac Delta function $\delta(\cdot)$ to describe the instillation of a dose of $BCG$ and $IL-2$ into the bladder. Using this function, makes calculations very difficult and in some cases impossible to implement different asymptotic methods as we mentioned in the introduction. Hence, we proposed a new function (see Equation (2.11) bellow), instead of the Delta function, in order to overcome all the drawbacks listed above.

Under the above assumptions and descriptions of the mechanism of the $BCG$-$IL2$ immune effect on the bladder the overall dynamics model and the interactions among ten different biological elements were presented the form:

where $\mathcal{J}(t) = \mathcal{J}(t; \alpha, \beta)$ is the sum of variation of Gamma Distributions functions with two free parameters $\alpha$ and $\beta$ of the form

For more details about the Gamma Distributions functions please refer to the Appendix.

$Comment$: At equation (2.1) we take into account a combination of different dose of $BCG$ at a different period of BCG pulsing using the function $\mathcal{J}$. At equation (2.7) we used with the same function $\mathcal{J}$ but only with one form of protocol: same dose of $BCG$ at the same period of BCG pulsing that the $BCG$ is instilled. From biological point of view, this means that we take into account different dosages of $BGC$ at different times but $IL-2$ is given to the patient at the same time when the $BCG$ is given but during all the treatment the same dose of $IL-2$ is given.

The parameter $b_{k}$ indicates the dose $BCG$ that instilled $m$ times into the bladder inserted through the urethra every $t$ days-period of BCG pulsing. In our analysis, we take into account different combinations of the amount of $b_{k}$ instilled into the bladder at different periods of BCG pulsing. Figure (1.1) present the function $\mathcal{J}$ that indicates given the same dose of $BCG$ at the permanent period of BCG pulsing. Figure (1.2) present the function $\mathcal{J}$ that indicate a different dose of $BCG$ at the same period of BCG pulsing. Figure (1.3) present the function $\mathcal{J}$ that indicates given a different dose of $BCG$ at different periods of BCG pulsing. Figure (1.4) present the function $\mathcal{J}$ that indicates given the same dose of $BCG$ at different periods of BCG pulsing.

The initial conditions of the problem are

$Comment:$ We assume that the tumor cells is an $ellipsoid$ form with one radius of $1$ $mm$ ^[25], hence,

Each $ml$ contains $10^{9}$ cells, and $1$ $ml$ is $1000$ $mm^{3}$. For example, lesion of tumor $13\times 15$ $mm$ in diameters, i.e. $r_{1} = 6.5$, $r_{2} = 7.5$, and hence $T_{u}(0) = 204$ $mm^{3}$ or contains about $0.204\times 10^{9}$ cells.

2.1.   Mathematical and biological interpretation of Equations (2.1)–(2.10), ^[24]

Equation (2.1) describes the dynamical rate of $BCG$ level changes with time. This equation includes the function $\mathcal{J}$ that describes the combination between the different dose of $BCG$ that is instilled into the bladder via a catheter inserted through the urethra at different periods of BCG pulsing. In addition this equation comprised of a positive term corresponding to $BCG$ instillations, and of negative terms corresponding to the elimination of $BCG$ by antigen-presenting cells ($APC$s) 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$.

Equation (2.2) is the dynamics of non-activated $APC$s, as described in B-M et al. ^[24], it is governed by two positive terms and three negative terms. The first positive term describes the normal influx of $APC$s to the tumor at a constant rate $\gamma$. The second positive term describes the recruitment of $APC$s due to bacterial infection at a rate coefficient $\eta$. The first negative term describes the activation of $APC$s 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-stages 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. The localized inflammatory response then attracts $APC$s, such as macrophages, which in turn eliminate uninfected tumor cells, according to the rate $\theta p_{3}$.

Equation (2.3) describes the dynamics of $BCG$-activated $APC$s. It is described by one positive term and two negative terms. The positive term is proportional to the number of non-activated APCs as well as BCG bacteria, with rate coefficient $p_{1}$ (as in Equation (2.1)). The first negative term is the migration of the infected, activated $APC$s to the draining lymphoid tissues, at a rate of coefficient $\beta_{1}$. The second negative term is the death of activated $APC$s at a rate of coefficient $\mu_{A_{1}}$.

Equation (2.4) describes the tumor-Ag-activated $APC$ ($TAA-APC$) dynamics. It is comprised of one positive term and three negative terms. The positive term describes the $APC$s which were activated by tumor antigen. The first negative term represents the tumor-Ag-activated $APC$s cells which destroy the uninfected tumor cells, with a rate coefficient $\lambda$. This term is multiplied by an $IL-2$ dependent parameter with a saturation constant $g_L$, 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}}$.

Equation (2.5) describes the dynamics of effector $CTL$s that react to $BCG$ infection. It is comprised of their migration rate, determined by their creation in the lymph nodes 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 the inactivation of effector $CTL$s 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_{1}}$.

Equation (2.6) describes the dynamics of effector cells reacting to tumor $Ag$. 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$ (as in Equation (2.5)). The first negative term describes the inactivation of effector $CTL$s 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_{2}}$.

Equation (2.7) describes the $IL-2$ dynamics. It is driven by a natural source, an external source, as well as sinking and degradation sources. The first two processes are positive and the last two are negative. They assume equal expression at the constant rate coefficient $q_1$. This equation includes the function $\mathcal{J}$ as we described in Equation (2.1). $I_2$ is consumed by $APC$s and $CTL$s. They assume that the rate of consumption is similar for both types of cells and denote its coefficient by $q_2$. The consumption depends on $I_2$ and is limited in a Michaelis-Menten fashion, with the Michaelis constant $g_L$. They also introduce $\mu_{I_2}$, the $I_2$ degradation rate coefficient.

Equation (2.8) describes the dynamics of infected tumor cells depended on two mechanisms. The first corresponds only to the rate of bacterial infection of uninfected tumor cells, ($T_{u}$), according to rate coefficient $p_{2}$. The second mechanism is the elimination of infected tumor cells $(T_{i})$ by their interaction with $BCG-CTL$ effector cells (represented by $E_{B}$), at rate coefficient $p_{4}$.

Equation (2.9) describes the dynamics of uninfected tumor cells. It is comprised of three processes: one positive term, corresponding to natural tumor growth and two negative terms, corresponding to 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 $APC$s 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$. The dependence in the equation of $T_u$ on $F_\beta$ is decreasing from $1$ to $a_{T, \beta}$ with Michaelis constant $b_{T, \beta}$ ^[33]. And then, there is a multiplication of those terms by an $I_2$-dependent Michaelis Menten term, with Michaelis parameter $g_T$, to propose that in the absence of $I_2$, $T_u$ cellular death does not occur. Since the tumor produces a variety of mechanisms in the biological settings that curtail the success of effector cell activity, they multiply $I_2 / (I_2 + g_L)$ by $g_T / (T_u + g_T)$, to denote the inversely proportional reduction in effector cell activity rate, such that when $T_u = 0$ the term is equal to $1$ and when $\lim_{T_u\rightarrow \infty} g_T / (T_u + g_T) = 0$. Note that although this factor can, in principle, nullify the efficacy of $CTL$s, this is not observed in cases of interest because $T_u \leq k$ ^[24].

Equation (2.10) describes the dynamics of a transforming growth factor-beta ($TGF$-$\beta$, $F_\beta$), as proportional to the tumor cell population, $T_u$, with $a_{\beta, T}$ as a proportion coefficient and is destroyed at a rate of $\mu_{\beta}$ proportional to $F_{\beta}$.

2.2.   The effect of the function $\mathcal{J}$ on the bladder cancer model

In this section, we present the effect of the function $\mathcal{J}$ on the model presented by ^[24].

Figures (1) present the different combinations of $BCG$ at different times.

Figures (2)-(11) present the solution profiles of the model for different combinations of the $BCG$ doses and periods of BCG pulsing (different/equal doses of $BCG$, different/ equal periods of BCG pulsing that the $BCG$ introduced into the bladder).

Figure (2) presents the dynamical change of $BCG$ in time. As we can see from the solution profiles of $BCG$ for different combinations of $BGC$ and time, the graphs behave in a wavy manner. This wavy behavior is due to the periodic function $\mathcal{J}(t)$ that takes into account different/ equal dose of $BGC$ at different/ equal period of BCG pulsing. At the first stage, before the treatment, a patient is uninfected by $BCG$ then $B(0) = 0$. Then after the $BCG$ instillation into the bladder and via a catheter inserted through the urethra. The $BCG$ stimulates the body's immune system to destroy cancer cells, and then it leaves the body (by antigen-presenting cells ($APC$s), in addition $BCG$ concentration decreases as a result of natural decay) in an exponential manner and repeat the process several times. Figure (2.1) refers to the solution profile of $BCG$ with respect to the same amount at the same time. The $BCG$ decreases to zero at the same rate of amount and time. Figure (2.2) refers to the solution profile of $BCG$ with respect to different $BCG$ but at the same time. At the beginning of treatment, the range of change of $BCG$ is very large, and after that this range decreases to zero with respect to the same time. Figure (2.3) refers to the solution profile of $BCG$ with different amounts at different time. According to this figure one can see that the quantities of $BCG$ are change with time and reduced to zero according to the varying of the amount of $BCG$. Figure (2.4) refers to the solution profile of $BCG$ with the same quantities of $BCG$at different time. It can be seen from this graph that indeed the $BCG$ is reduced to zero but the amounts of the $BCG$ are approximately the same with respect to time.

Figure (3) presents the dynamics of non-activated $APC$s denoted by $A$. In all possible combinations of amount of $BCG$ and time, we see that $A$, which starts from the initial condition $10^{4}$, decreases very fast, but after that $A$ increases and decreases periodically throughout the treatment and finally stabilized. The biological explanation for the reduction of the variable $A$ from the initial condition is as follow: the activation of $APC$s by $BCG$ that described by the term $p_{1}AB$, and the capture of uninfected tumour cells ($Tu$) by $APC$s as well as the natural cell death of $A$, all these terms are very dominant compared to the sum of the term $\eta AB$, which describes the recruitment of $APC$ due to bacterial infection, and the parameter $\gamma$, we obtain that the dynamical variable is decreasing to its stable value.

Figure (4) presents the dynamics of the variable $A_{B}$ the ($APC$s) cells that are activated by $BCG$. At the beginning of the process, there is a very large activity and interaction that behave exponentially and then a decrease of the activity as a result of the decline of the $APC$ as seen from the graph $(3)$. In addition, there is a decrease in activity as a result of migration of the infected and activated $APC$s to the draining lymphoid tissues, as well as the death of activated $APC$s, but this activity does not cause the variable to change. On the contrary, $A_{B}$ stabilizes at some constant value.

Figure (5) presents the dynamics and interaction of Tumor Associated Antigen ($TAA$) with $APC$s cells. In general, $T$ cells are activated by the link between $TCR$ ($T$-Cell Receptor) and the tumor associated tumor ($TAA$) associated with the $MHC$ (Major Histocompatibility Complex) which is presented by the Antigen ($APC$ Presenting Cell), which is the initial "signal" of the $T$ cell. The biological dynamics of this variable include a "competition" between the interaction of non-activated $APC$s, uninfected tumor cells and $IL-2$ to the migration of the $TAA-APC$ to the draining lymphoid tissues, and the natural death of $TAA-APC$. At first stage, the interaction $TAA-APC-IL-2$ overcomes to migration of the $TAA-APC$ and the natural death of $TAA-APC$ and after that, the process is reversed and $A_{T}$ decrease.

Figure (6) presents the dynamical interaction of $BCG$ and $CTL$ effector cells. The biological dynamics of this variable $E_{B}$ is not like the normal dynamics of increase and decrease smoothly and continuously, but during the treatment, the interaction between $BCG$ and $CTL$ increases and decreases alternately until this interaction decreases. Since the main interaction depends on the migration rate that is determined by the creation in the lymph node and migration to the bladder, inactivation rate and their death rate. The term that causes to the variable $E_{B}$ to grow is an expression obeying Michaelis$-$Menten saturated function's. This expression depends directly on $A_{B}$, $IL-2$ and Michaelis$-$Menten constant $g$. The parameter $g$ appears in the denominator of that expression. The numerator of this expression containing the rate coefficient $\beta_{B}\approx 10^{8}$. This ratio between these two parameters, makes this expression very small relative to the (two) other terms of the equation that appear in it with a negative sign, what causes the derivative of $E_{B}$ to be negative and positive in certain points. As we will see later of the analysis of the results, the biological dynamics of $IL-2$ rise and fall in a wave manner with a relatively small wavelength and very big amplitude. This amplitude is the main cause of the oscillations of the variable $E_{B}$. The same behavior can be seen in $A_{T}$'s graph but in a much less prominent way. When we reduced the resolution of this graph we did see wavy behavior but with a very small wavelength and a very small amplitude. This is because of the interaction between $BCG$ and $CTL$, where $BCG$ have also a wavy behavior.

The wavelength and amplitude of graphs $(6.3)$ and $(6.4)$ which refer to different doses of $BCG$ at different periods of $BCG$ pulsing and the same dose of $BCG$ at permanent periods of BCG pulsing respectively, are more prominent compared to the other combinations of $BCG$ and time, where the amplitude there and wavelength are smaller. In addition, the amplitude and wavelength of these graph are stabilized approximately to a certain value. In Figure (6.1) which refers to same dose of $BCG$ at different periods of BCG pulsing, $E_{B}$ increases periodically very fast, then it is stabilized to a constant value, finally it decreases to zero very fast. The graph $(6.2)$, which refers to the different dose of $BCG$ at the same period of BCG pulsing, is on a rise in cycles and does not decrease. The explanation for this behavior is clear because during the treatment the interaction of $BCG$ and $CTL$ grows and goes down due to $BCG$'s periodicity.

Figure (7) presents the solution profile of the variable $E_{T}$ which describes the Tumor-effector $CTL$s interactions. The dynamics of $E_{T}$ are comprised of their migration rate, inactivation rate and death rate. The behavior of this variable is analogous to variables that periodically cycle through growth and decline. $E_{T}$in Figures (7.2-4) behaves almost similarly. The different graph is the first one $(7.2-4)$ which refers to same amount of $BCG$ at the same time. The amplitude, as well as the wavelength of this graph, are smaller and denser than other graphs. Then drop fast rather than cyclic to zero as the others graphs.

Figure (8) presents the dynamics of the variable interleukin ($IL-2$). $IL-2$ is a lymphokine that induces the proliferation of responsive $T$ cells. In addition, it acts on some $B$ cells, via receptor-specific binding as a growth factor and antibody production stimulant. The main aim of $IL-2$ is the importance in stimulation of the adaptive immune system. In the equation of $IL-2$ we insert the function $\mathcal{J}(t)$ that indicate different combination of $BCG$ doses and periods of BCG pulsing. The natural sources of $IL-2$ are the immune cells, i.e., the activated $APC$s and $CTL$s. In this work, an external source have been added to the treatment using the function $\mathcal{J}(t)$. In our numerical simulaiton we used with the same coefficient $b_{k}$ (which indicates the amount of external source) for both $BCG$ and $IL-2$. $IL-2$ is injected into the bladder and modeled in the same manner as for $BCG$ instillations. Graphs $(7.1)$ and $(7.4)$ behave the same (except the wavelength), a constant increase and decrease throughout the treatment until a final reduction at the end of the process. Graphs $(7.2)$ and $(7.3)$ also behave the same with the similar amplitude and wavelength. $IL-2$ increases linearly during the treatment in periodically manner. The optimal solution for effective treatment is the growth of $IL-2$ in controlled form. I.e., growth and stabilizing at a certain stage. As we can see, there is no graph solutions for an optimal behave of $IL-2$. Because there are graphs that only rise and on the other hand there are graphs that only stabilize. Hence, by applying the $SPVF$ method we obtain the optimal linear combinations of the dynamical variables of the system for an optimal treatment.

Figure (9). $T_{i}$ is the infected (or "marked" by the activation of $APC$ ($A$), $BCG-APC$ ($A_{B}$), and $BCG-CTL$ ($E_{B}$)) tumor cells. At the beginning of the $BCG$ treatment, there are no infected tumor cells and hence $T_{i} = 0$. After that, the number of marked cells increases because the cancer cells are marked by the others dynamics variables of the system. On the other hand, the cells $T_{i}$ are eliminated due to the interaction with $BCG-CTL$ effector cells, and hence the graph of $T_{i}$ decreases. The solutions profiles of $(9.2)$-$(9.4)$ are optimal solutions since $T_{i}$ cells decrease to zero. The solution of $T_{i}$, which refers to same amount of $BCG$ at different times is not optimal, on the contrary, is worse for the patient. Because at some stage, although the tumor does not increase, it stabilizes, what can develop the disease and the recurrence of the tumor cells.

Figure (10) is the solution of uninfected tumor cells. At the beginning of the treatment ($t = 0$) the initial number of uninfected tumor cells in the urothelium before $BCG$ treatment taken to be $T_{u}(0) = 2\cdot 10^{6}$. After that $T_{u}$ decreases as expected to zero. The number of uninfected cells is decreasing and they become infected cells because of the interaction with the $BCG$, i.e., there is a transitions of cells from $T_{u}$ to $T_{i}$ by $BCG$ action/ interaction. The growth of $T_{u}$ is due to the natural tumor growth that proportional to $T_{u}^{2}$ and $T_{u}$ according to logistic tumor growth rate. According to ^[24], the number of $T_{u}$ cells is limited by the tumor cell maximal number $K = 10^{11}$. Our simulation results show that $T_{u} = 2\cdot 10^{6} < < k$. The decrease of $T_{u}$ is due to the tumor infection by bacteria as well as tumor elimination by immune system cells. The interactions and activations of the dynamical variables of the system such as $APC$ ($A$), $TAA-CTL$ effectors ($E_{T}$), $IL-2$ that dependent Michaelis-Menten term with a Michaelis parameter $g_{L}$ cause the decrease of $T_{u}$.

Figure (11) describes the dynamics of a transforming growth factor-beta which behaves as the variable $T_{u}$, i.e., growth very fast and decrease to zero.

$General$ $comment$: The model's solutions in the new coordinates are irrelevant and biologically insignificant. What matters is the exposing of the system hierarchy and hence the decomposition of the model into fast and slow subsystems. And hence, the recipient of the linear combination of the variables of the (old) system with appropriate coefficients.

As we can see from the above analysis of the results, one can not obtain the optimal treatment for the Bladder cancer (It is not known what combination of quantities with different or equal period of $BCG$ and $IL2$ pulsing). Therefore, according to the $SPVF$ method one can extract the optimal combination between the variables of the (old) system to achieve the maximal effect of the treatment and cause the elimination of the cancer cells.

3.   Analysis and results

In this section we apply the Singular Perturbed Vector Field ($SPVF$) method to the bladder cancer model (2.1)–(2.10). The algorithm of the $SPVF$ method can be found in paper ^[20].

3.1.   The application of the $SPVF$ method

According to $SPVF$ method the eigenvalues received by the algorithm of the $SPVF$

According to the algorithm of the $SPVF$, the maximum gap is $gap_{max_i} = \frac{\lambda_{2}}{\lambda_{3}} = 8.292\cdot 10^{10}$. This means that the original system of equations can be decomposed into fast and slow subsystems, where the fast directions of the system are in the directions of the eigenvectors $\vec{v}_{1}$ and $\vec{v}_{2}$ that corresponds to the eigenvalues $\lambda_{1}$ and $\lambda_{2}$.

The corresponds eigenvectors are

,

The above results show clearly the hierarchy of the given system in the new coordinates presentation. The fast direction of the system is in a parallel direction to the eigenvector $\vec{v}_{1}$. Then the less fast direction of the system is in a parallel direction to the eigenvector $\vec{v}_{2}$ and so on according to the descending order of eigenvalues since $\lambda_{1} > > \lambda_{2} > >, ..., > > \lambda_{10}$. The slow directions of the system are in parallel directions to the eigenvectors $\vec{v}_{3}$, $\vec{v}_{4}$, ..., $\vec{v}_{10}$ in descending order, where $\vec{v}_{10}$ is the slowest direction of the system.

The next step of the $SPVF$ method is transforming the model (2.1)–(2.10) using the above eigenvectors, hence using matrix in a form we write

where $t$ refers to the transpose operation, i, e., it means that each eigenvector is placed in a row in the above matrix.

For the convenience of the reader, we define the above system as matrix form hence, we denote the above matrix as $\mathcal{A}$, the vectors of the variables of the original system as $\vec{V} = \left(B, A, A_{B}, A_{T}, E_{B}, E_{T}, I_{2}, T_{i}, T_{u}, F_{\beta}\right)$, and the vector of the new variables as $\vec{U} = \left(x_{1}, x_{2}, x_{3}, x_{4}, x_{5}, x_{6}, x_{7}, x_{8}, x_{9}, x_{10}\right)$. Hence, the system (3.5) can be rewritten as

The next step is expressing the old variables of the system as a function of the new variables. For this aim, we multiply the set of equations (3.6) by the inverse matrix of $\mathcal{A}$

In order to write the system of the $ODE$'s model in the new coordinate, we take the derivative of the system (3.6) with respect to time, i.e.,

and substitute the expressions of the $RHS$ (right hand side) from the system (2.1)–(2.10) instead of $\frac{d\vec{V}}{dt}$ in (3.8). In details

where

is the $RHS$ of the system (2.1)–(2.10). Substitute Equation (3.7) into Equation (3.9) and receive the model in the new coordinates with the initial conditions as

3.2.   The decomposition of the model into fast and slow subsystems

The model (2.1)–(2.10) that is rewritten in new coordinates as (3.11), is the model with the required hierarchy, and hence can be decomposed to fast and slow subsystems as follows

where $0 < \varepsilon_{1}, \varepsilon_{2} < < 1$, and the initial conditions are

The mathematical, and mainly biological significance of this decomposition is that it enables us to reduce the model and to investigate the dynamics of the system only on the fast subsystem without losing information if we would solve the complete system that includes all the variables of the system. The reduced model decreases the run-time of calculations and in some cases, one can even obtain an analytic expression for different variables of the reduced system.

The next analysis is investigating the fast subsystem. We take into account only the variables that refer to the fast directions of the system i.e.,

where $\delta, \xi < < 1$. We take into account only three coordinates from the eigenvectors and neglect the rest because the other coordinates belong to the order of magnitude $< 10^{-7}$. The coordinates of the eigenvectors are actually the desire optimal combination, i.e., in the direction of $\vec{v}_{1}$ the dominant variable of the system is $A$ (this variable will receive the largest weight) since its coefficient larger compared to the coefficients of $B$ and $A_{B}$. The other two variables will approximately the same weight since $\frac{0.4192}{0.4192} = 1$. The coefficients of the above linear combination indicate the proportions (relations) of the quantities (or percentages) that should be taken for the considered treatment for each variable of the old variables. The same analysis can be applied for the variable $x_{2}$, which in the direction of the vector $\vec{v}_{2}$.

4.   Preliminaries to fast and slow systems, stability analysis

In this section, we present a brief introduction to the fast and slow systems and investigate the stability of the considered model using the $SPVF$ method.

4.1.   Introduction to fast and slow system

Once we apply the $SPVF$ method, we can decompose the system into fast and slow subsystems. This means that the investigation of the model can be concentrated on the reduced model, i.e., one can investigate the stability analysis of the fast subsystem without loosing the information of the original system.

Given a model that presented with a hidden hierarchy, by applying the $SPVF$ method one can present the model with an explicit hierarchy, i.e., $SPS$ form as

where $0 < \varepsilon < < 1$, $\vec{v} = (v_{1}, ..., v_{m_{fast}})\in\mathbb{R}^{m_{fast}}$ called fast variables which change fast compared to the slow ones $\vec{u} = (u_{1}, ..., u_{k_{slow}})\in\mathbb{R}^{k_{slow}}$ and $m_{fast}+k_{slow} = n$ ($n$ is the dimension of the system. The functions $\vec{F}$ and $\vec{H}$ are assumed to be $C^{\infty}$ of $\vec{v}$, $\vec{u}$ and $\varepsilon$ in $(V\times U)\times I$ where $V$ is an open subset of $\mathbb{R}^{m_{fast}}$, $U$ is an open subset of $\mathbb{R}^{k_{slow}}$ and $I$ is an open interval containing $\varepsilon = 0$.

Setting $\tau = \varepsilon t$, the above model is transformed to

Since $0 < \varepsilon < < 1$ which describes the difference in time scales, than system (4.1) and (4.2) is called a $slow$ and $fast$ $system$, respectively.

These systems are equivalent whenever $\varepsilon\neq 0$ and they are labeled $singular$ $perturbation$ $problems$ when $0 < \varepsilon < < 1$. The label $"singular"$ stems in part from the discontinuous limiting behavior in the system (4.2) as $\varepsilon\rightarrow0$. At this case, when $\varepsilon\rightarrow0$, system (4.1) leads to differential-algebraic system

called $reduced$ $slow$ $system$ which dimension decreases from $m_{fast}+k_{slow}$ to $k_{slow}$.

On the other hand, system (4.2) reduces to an $m_{fast}$-dimensional system called $reduced$ $fast$ $system$, with the variable $\vec{u}$ as a constant parameter:

By exploiting the decomposition into fast and slow reduced systems (4.4) and (4.3), the geometric approach reduced the full singularly perturbed system to separated lower-dimensional regular perturbation problems in the fast and slow regimes, respectively.

In this paper we apply the $SPVF$ method and present the considered model at the form of system (4.1).

According to the manifolds theory approximation, the slow manifold is obtained by setting $\varepsilon = 0$ in system (4.1)

In addition, we assumed that the solutions of (4.5) tend toward an equilibrium

which is the root of equation (4.5).

The solutions of the system (4.1) have a fast transition from the initial points (or initial conditions) of the system to a point of the slow manifold $M_{slow}$. Hence, at the system (4.1) a slow motion takes place on the slow manifold, according to the equation

This equation is called the $reduced$ $problem$ or the $reduced$ $model$. This description of system (4.1) was given by $Tikhonov$ ^[27].

Follow the above theory of fast and slow system, the equilibrium points are obtained by solving the fast sub-system for the fast variables while setting all the slow variables as constant (the slow variables are "frozen" at their initial values and hence one can take $\vec{u} = \vec{u}(0)$ as the constants);

Here, we denote with the star notation the equilibrium points which indicate the fast variables at their equilibrium values. Substitute the equilibrium points of the fast variables into the slow subsystem and solve this subsystem to find the equilibrium points of the slow variables, i.e.,

In order to analyze the stability of the above equilibrium points, one should examine the sign of the real part of the eigenvalues of the following Jacobian matrix at each equilibrium point of (4.8)–(4.9)

The eigenvalues with negative real parts are locally stable.

Another method to find out the equilibrium points and their stability is to concentrate only on the fast subsystem, while the rest of the slow variables remaining constant. The advantage of this method is that we obtain more simple algebraic expressions and we reduce the running time of computer work. At this method, we solve only the fast subsystem (4.8) and the stability of the equilibrium points of the fast variables is given by analyzing the Jacobian matrix of the fast subsystem

which is a $m_{fast}\times m_{fast}$ matrix of partial derivatives of the fast subsystem with respect to the fast variables $\vec{v}$ while the slow variables are taken as constants.

This means that the information about the stability of the model can be extracted from the analysis of the fast manifold.

At the next section we apply the first method.

4.2.   Stability analysis of the model

In this section, we investigate the stability of the model that is the decomposed model into fast and slow subsystems.

The following steps are implements for finding the fixed points in the new coordinates and determining their stability.

$1.$ substitute the slow variables as a constant into the fast subsystem (one can take the initial condition of the slow variables as the constants).

$2.$ setting all fast variables derivatives to zero i.e., solve the fast subsystem (with slow variables as constants) and find the equilibria points of the fast variables $(x_{1}^{*}, x_{2}^{*})$

(here, we denoted the fast equilibrium points with star notation). Here, we have two equations with two unknown variables $x_{1}^{*}, x_{2}^{*}$.

$3.$ substitute the equilibrium points from step $2$ ($x_{1}^{*}, x_{2}^{*}$) into the slow subsystem, solve the slow subsystem and find the equilibrium points of the slow variables $(x_{3}^{*}, ..., x_{10}^{*})$

Here, we have eight equations with eight unknown variables $x_{3}^{*}, ..., x_{10}^{*}$.

$4.$ substitute the equilibrium points from steps $2$ and $3$ at the Jacobian matrix of the full system (the model at the new coordinates).

$5.$ compute the eigenvalues of the Jacobian matrix of the system (the model at the new coordinates) for each set of equilibrium point (the stable points are those with a negative real part of the eigenvalues).

$6.$ transform only the equilibrium points that are stable from steps $2$ and $3$ to the original coordinates using the inverse matrix of the eigenvectors, i.e., compute

where $i$ indicateS for different stable equilibrium points.

It should be noted here that the equilibrium points we have received in the new coordinate system are functions of the parameters of the original model and a function of the initial conditions of the new model, i.e.,

$Comment$ 1: the initial condition as well as the equilibrium points $\left\{x_{1}^{*}, x_{2}^{*}, x_{3}^{*}, ..., x_{10}^{*}\right\}$ are biologically meaningless as there is no biological significance to the variables in the new coordinates. Hence, in order to understand the biological meaning of the above results, we should transform these points to the original coordinates using the inverse matrix of the eigenvectors $\mathcal{A}^{-1}$. 2: Since the eigenvalues are invariant under change of coordinates, hence the stability equilibrium points will be remain stable equilibrium points at the original coordinates.

By applying the above steps $1-6$, we have found the following stability equilibrium point express in the original coordinates as a function of the parameters of the original model and the initial conditions at the new coordinates (For convenience, we denoted the initial condition $x_{i}(0)$ as $x_{i}^{0}$)

Here, we neglect elements of an order smaller than $10^{-4}$.

As we have mentioned before, this equilibrium point is stable since the eigenvalues are invariant under change of coordinates.

In order to complete writing the above stable equilibrium points in the original system, we write the initial conditions $\vec{U}$ as a function of the original initial conditions $\vec{V}$ using the matrix $\mathcal{A}$ as $\vec{U}(0) = \mathcal{A}\vec{V}(0)$ or in an explicit form

By substitute equation (4.17) into equation (4.16), we receive the equilibrium point written completely in the original coordinates of the model. Our analysis for the equilibrium point includes different values of the function $\mathcal{J}(t_{0}; \alpha, \beta)$ i.e., different combinations of dosages and periods of BCG pulsing.

Here, in table 1–4, we present the stable equilibrium point of the dynamical variable $T_{u}^{*}$ for different values of $b_{k}$ at different periods of BCG pulsing (after $30$ days, $44$ days, $58$ days and $65$ days).

Every table presents the stable equilibrium points expressed in the original variables of the bladder cancer model for different values of the Gamma functions, i.e., for different combinations of dosages and times. For example, 2 presents the stable equilibrium point for the function $\mathcal{J}^{(2)}$ which relates to the sub-graph $2$ at Figure 1 that indicates a different doses of $BCG$ at the permanent period of BCG pulsing: after $30$, $44$ $58$ and $65$ days. According to these results, we can see that this combination of treatment is optimal since it causes the variable $T_{u}$ (cancer) to decrease to zero very fast compared to the other combinations of the function $\mathcal{J}$, and this result is consistent with the results obtained from the numerical simulations presented in the graph of $T_{ut}$.

5.   Conclusions

In this study, we improved the mathematical model for the treatment of bladder cancer using the probability density function in the form of a gamma function to describe the introduction of BCG and IL-2. Gamma function allows the delivery of treatment at different periods of pulsing and at different doses. The mathematical model includes ten nonlinear ordinary differential equations of first order. The model describes the dynamics of the combination between $BCG$ and $IL-2$ influence on the bladder tumor. The model describes the dynamics variability of the variables of the system involved in the interaction between the variables of the immune system and the variables of treatment in a way that the hierarchy is not explicit. We solved the model using numerical simulations for varying amounts of $BCG$ and at varying period of BCG pulsing. Outcome following $BCG$ therapy depends on the dose and period of pulsing ^[28].

In order to obtain an optimal solution, we implemented the algorithm of the $SPVF$ method. This method transforms the model, using eigenvectors, to a model presented in new coordinates in which the system can be splitted into a fast subsystem and a slow subsystem according to the hierarchy eigenvalues of the system. This hierarchy enables one to study the model using various asymptotic methods without losing any mathematical or biological information from the original system.

After rewriting the model in the new coordinate using the eigenvectors received from the $SPVF$ method, we exposed the hierarchy of the new model and received two fast variables and eight slow variables. This decomposition enabled us to apply the stability analysis to the fast subsystem. We found one stable equilibrium point express by the parameters of the original system and the initial conditions in the new coordinates. In order to obtain a biological meaning of these equilibrium point, we inversed and transformed the stable equilibrium points using the inverse matrix of the considered eigenvectors. And since the eigenvalues are invariant under change of coordinates, the stable equilibrium points remain stable at the original coordinates. After, applying this procedure, we received the equilibrium points expressed by the parameters of the original model and the initial conditions of the original model. The equilibrium points are also depend on the function $\mathcal{J}$ which takes into account different combination of dosages and periods of the pulsing therapy. We have found that the optimal combination of different dosages and periods of $BCG$ pulsing at $30$, $44$, $58$ and $65$ days after the start of the treatment.

Due to the $SPVF$ method, we have explicitly received the equilibrium points, since this method enables us to reduce the dimension of the original model. We compared the results of the equilibrium point to the results we received from the model and found that they corresponded to our numerical simulations as one can see from the graphs.

Acknowledgments

This work was supported by Jerusalem College of Technology.

Conflict of interest

The authors declare that there is no Conflict of interest.

Appendix

$Gamma$ $distribution$

The general formula for the probability density function of the gamma distribution is:

for $x$ $\geq \mu$, $\alpha, \beta>0$, where $\alpha$ is the shape parameter, $\mu$ is the location parameter, $\beta$ is the scale parameter, and $\Gamma$ is the gamma function which has the formula

with the following properties:

$1.$ $\Gamma(1) = \int_{0}^{\infty}e^{-x}dx = 1$

$2.$ $\Gamma(\alpha+1) = \alpha\Gamma(\alpha)$

$3.$ $\frac{\Gamma(\alpha)}{\lambda^{\alpha}} = \int_{0}^{\infty}x^{\alpha-1}\cdot e^{-\lambda x}dx$

$4.$ $\Gamma(n) = (n-1)!$, $n = 1, 2, ...$

$5.$ $\Gamma\left(\frac{1}{2}\right) = \sqrt{\pi}$.

The case where $\mu = 0$ and $\beta = 1$ is called the standard gamma distribution and has the form of

for $x$ $\geq 0, \alpha>0$.

$Parameters$

$\mu_{A} = 0.038$ APC half life $[days^{-1}]$

$\mu_{A_{1}} = 0.138$ Activated APC half life $[days^{-1}]$

$\mu_{E_{1}} = 0.19$ Effector cells mortality rate W/O $IL-2$ $[days^{-1}]$

$\mu_{E_{2}} = 0.034$ Effector cells mortality rate $IL-2$ $[days^{-1}]$

$\mu_{B} = 0.1$ $BCG$ half life $[days^{-1}]$

$p_{1} = 1.25\cdot 10^{-4}$ The rate of $BCG$ binding with APC $[cells^{-1}][days^{-1}]$

$p_{2} = 0.028\cdot 10^{-6}$ Infection rate of tumor cells by $BCG$ $[cells^{-1}][days^{-1}]$

$p_{3} = 1.03\cdot 10^{-10}$ Rate of $E$ deactivation after binding with infected tumor cells $[cells^{-1}][days^{-1}]$

$p_{4} = 1.1\cdot 10^{-6}$ Rate of destruction of infected tumor cells by effector cells $[cells^{-1}][days^{-1}]$ $\lambda = 0.01\cdot 10^{-6}$ Production rate of $TAA-APC$ $[days^{-1}]$

$\beta_{B} = 1.45\cdot 10^{8}$ Recruitment rate of effector cells in response to signals released by $BCG$-infected and activated $APC$ $[cells^{-1}][days^{-1}][I_{2}^{-1}]$

$\beta_{T} = 1.51\cdot 10^{6}$ Recruitment rate of effector cells in response to signals released by $TAA$-infected and activated $APC$ $[cells^{-1}][days^{-1}][I_{2}^{-1}]$

$\gamma = 4700$ Initial $APC$ cell numbers $[cells^{-1}][days^{-1}]$

$\eta = 2.8\cdot 10^{-6}$ Rate of recruited additional resting $APC$s $[cells^{-1}][days^{-1}]$

$r = 0.0048-0.0085$ Tumor growth rate $[days^{-1}]$

$\beta = 0.034$ Migration rate of $TAA-APC$ and bacteria activation $APC$ to the lymph node $[cells^{-1}][days^{-1}]$

$\alpha = 3.7\cdot 10^{-6}$ Efficacy of an effector cell on tumor cells $[cells^{-1}][days^{-1}]$

g = 10$^{13}$ Michaelis-Menten constant for $BCG$ activated $CTL$s and for $TAA-CTL$s $[cells]$

$g_{T} = 5.2\cdot 10^{3}$ Michaelis-Menten constant for tumor cells $[cells]$

$k = 10^{11}$ Maximal tumor cell population $[cells]$

$q_{1} = 0.007$ Rate of $IL-2$ production $[cells^{-1}][days^{-1}]$

$q_{2} = 1.2\cdot 10^{-3}$ The production of $IL-2$ used for differentiation of effector cells $IU$ production $[cells^{-1}][days^{-1}]$

$\mu_{I_{2}} = 11.5$ Degradation rate $[days^{-1}]$

$\theta = 0.01$ 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/cell^{-1}]$

$a_{T, \beta} = 0.69$ Michaelis-Menten saturation dynamics. The dependence of $F_{\beta}$ is decreasing from $0$ to $a_{T, \beta}$ $[dimensionless]$

$b_{T, \beta} = 10^{4}$ Michaelis constant $[pg]$

$\mu_{\beta} = 166.32$ The constant rate, account for degradation of $F_{\beta}$ $[days^{-1}]$

$g_{L} = 10^{4}$ Michaelis-Menten constant for $IL-2$ $[cells]$