PDE based geometry model for BCG immunotherapy of bladder cancer

1. Introduction and related work

Bladder cancer (BC) is the 10th most common form of cancer worldwide, with an estimated 549,000 new cases and 200,000 deaths. The highest incidence occurs in industrialized and developed areas, such as Europe, North America, and Australia (Bray et al., 2018). The primary cause of about half of bladder cancer cases is occupational exposure to chemicals in industrial areas processing paints, metals, dyes and petroleum products. Tobacco smoking and environmental carcinogens are another risk factor for bladder cancer (Bunimovich-Mendrazitsky et al., 2015a). 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., 2014).

Treatment of non-invasive BC has not advanced significantly over the past few decades following the treatment protocol suggested by Morales et al. (1976) that involves weekly instillations of Bacillus Calmette–Guérin (BCG) (Morales et al., 1976). BCG, an attenuated non-pathogenic strain of Mycobacterium bovis that was originally used as a vaccine against tuberculosis (TB), is a type of immunotherapy used to treat non-invasive bladder cancer (Herr et al., 1988; Simon et al., 2008; Redelman-Sidi et al., 2014). BCG immunotherapy has proven superior to chemotherapy in reducing the rate of tumor relapse (Wei, 2016). Although Lamm and others have found that BCG even reduces the progression of the disease (Lamm, 2006), it is necessary to understand why the standard BCG treatment protocol is not effective for non-responding or relapsing patients. The BCG treatment protocol has yet to be optimized specifically for those patients who do not achieve remission from treatment according to the standard scheme.

Mathematical modeling shown to be a useful tool in oncology, allowing to investigate both the disease and possible treatments (Bhattacharya et al., 2020; Jordão and Tavares, 2017; Hornberg et al., 2006). Several attempts have been made to develop the model for BCG treatment of BC as a response of the immune system to introduced bacteria into the bladder by means of Ordinary Differential Equations (ODE) in order to find the optimal treatment protocol (Shaikhet and Bunimovich-Mendrazitsky, 2018; Bunimovich-Mendrazitsky and Goltser, 2011; Bunimovich-Mendrazitsky et al., 2015b). In addition, several attempts were made to describe the cell dynamics taking into account biological interactions in the physical space based on partial differential equations (PDE) (Lazebnik et al., 2020; Matzavinos et al., 2004; Eikenberry et al., 2009; Fridman and Kao, 2014). One of them is the model is investigated by Fridman et al. (Fridman and Kao, 2014) which describes the case where the geometrical configuration is a sphere. The model describes the cell population dynamics including the immune system cells, cancer cells, healthy cells and disease-infected cells. Their model (Fridman and Kao, 2014) analyzes both logistic and exponential growth of tumor cells.

Another model describing BCG immunotherapy treatment dynamics that takes into consideration an approximation of the bladder’s geometry using PDE investigated by Lazebnik et al. (2020). Their model (Lazebnik et al., 2020) assumed continuous BCG instillation and logistic growth of tumor cells inside the bladder. In (Lazebnik et al., 2020) the changes in treatment protocol were studied by considering a sphere-ring approximation to the bladder’s geometry. Moreover, in this model a diffusion dynamics was added for all cell populations as BCG, tumor and effector cells. The model from (Lazebnik et al., 2020) suffered from numerical instability for long treatment time: cancer cell population shows divergence to infinity at the 35th day of the treatment.

In the research of Guzev et al. (2019) a BCG and interleukin 2 (IL-2) combined therapy was examined and presented validation of this protocol for BC patients (Guzev et al., 2019). However, the model from (Guzev et al., 2019) is lacking the geometrical understanding of the dynamics of the biological system and diffusion of cell population during therapy.

There are medical and biological investigations that show the prognostic significance of tumor location on survival outcomes in patients with BC (Grabnar et al., 2006; Weiner et al., 2019). Therefore, in order to obtain the best treatment protocol, it is important to consider the geometry of the bladder and the location of the polyps.

In this research, we propose a model that tackles the two main inefficiencies of the Guzev et al. (2019) and Lazebnik et al. (2020) models. This paper is organized as follows: in Section 2, we introduce our mathematical model of BCG treatment with the approximation of the geometrical configuration of the bladder. Afterwards, the numerical calculation used to solve the PDE and solution stability analysis. In Section 3, we present the model’s bifurcation points and their clinical impact. In Section 4, we offer a method to find the optimal time treatment protocol given the paitent’s start condition. In Section 5, we investigate the influence of the tumor depth and distribution in the urothelium on the optimal treatment protocol. In Section 6, we discuss the main clinical results arising from the model.

2. Mathematical modeling extension

Our aim in the development of this mathematical model is to allow clinical professionals to perform better and more personalized treatment in the scope of BCG immunotherapy of bladder cancer. The advantage of using a mathematical model is the ability to examine the biological system in relatively simple settings while producing clinical results which may be used on real patients. The mathematical model we present is a system of 10 s-ordered, nonlinear PDEs representing cell populations dynamic.

2.1. Model definition

The system of Eq. (1-14) represents the treatment of bladder cancer with BCG and IL-2 as the dynamics of cell populations (the values of the parameters that we used are shown in Table 3).

In Eq. (1), ^∂B(t,r) is the dynamical rate of BCG cell population distri∂t bution over time. It is affected by the following five terms. First, a quantity b of BCG instilled into the bladder every τ days. As the instillation of the BCG is modeled by a shifted Dirac delta function δ(t ˘mτ),m ∈ {0, 1,…,N^˘1}, the m_th dose raises B(t, r) by b units at t = mr. Second, the elimination of BCG by antigen presenting cells (APCs) according to the rate coefficient − p₁. Third, the BCG tumor cell groth at a rate coefficient − p₂. Fourth, the bacteria cell death with rate coefficient μ_β. Finally, the diffusion of the BCG cell population inside the bladder’s geometry is assumed to spread at rate coefficient D₁.

In Eq. (2) ^∂A(t,r) _∂t is the dynamic of nonactivated APCs. It is affected by the following six terms. First, the normal influx of APCs to the tumor at a constant rate λ. Second, the recruitment of APCs due to bacterial infection at a rate coefficient τ. Third, the activation of APCs by BCG at the rate coefficient − p₁. Fourth, the natural cell death at the rate coefficient, − τ_A. Fifth, the two-stage elimination of tumor cells, first by effector CTL activity on 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 APCs, such as macrophages, which in turn eliminate uninfected tumor cells, according to the rate − p₃. Finally, a diffusion of the APCs cell population inside the bladder’s geometry is assumed to spread at rate coefficient D₂.

In Eq. (3) ^∂AT(t,r) is the tumor-Ag-activated APC (TAA-APC) dynamic. ∂t It is affected by the following five terms. First, the APCs which were activated by tumor antigen. Second, the tumor-Ag-activated APCs cells which destroy the uninfected tumor cells, with a rate coefficient λ. This term is multiplied by an IL-2-dependent parameter with a saturation constant g_I, to propose that in the absence of IL-2, A_T production ceases, while in the presence of external IL-2, the production term is close to 1. Third, the migration of TAA-APC to the draining lymphoid tissues at a rate of coefficient − β₁. Fourth, the natural death of TAA-APC at a rate coefficient μ_A1. Finally, the diffusion factor of the TAA-APC cells population in the bladder geometry is assumed to be D₃.

In Eq. (4) ^∂AB(t,r) is the dynamic of BCG-activated APCs. It is affected ∂t by the following four terms. First, the number of nonactivated APCs as well as BCG bacteria, with rate coefficient p₁. Second, the migration of the infected, activated APCs to the draining lymphoid tissues, at rate coefficient β₁. Third, the death of activated APCs at rate coefficient A₁. Finally, the diffusion factor of the BCG-activated APCs cells population in the bladder geometry is assumed at coefficient D₄.

In Eq. (5) ^∂EB(t,r) is the dynamic of effector CTLs that react with BCG ∂t infection. It is affected by the following four terms. First, the migration term is proportional to A_B and IL-2, with a maximal rate coefficient β_B. This rate is brought to saturation by large numbers of A_B, using a Michaelis–Menten saturation function, with Michaelis parameter g. Second, the decrease in CTLS population size from the inactivation of effector CTLs via their encounter with infected tumor cells (T_i) at a success rate coefficient − p₃. Third, the decrease in the BCG-effector CTL (E_B) population size from the (E_B) cells’ natural death rate μ_E. Finally, the diffusion factor of the effector CTLs that react with BCG infection cell population in the bladder geometry is assumed at rate coefficient D₅.

In Eq. (6) ^∂ET(t,r) is the dynamic of effector cells reacting with tumor ∂t Ag. It is affected by the following four terms. First, the migration element is proportional to A_T and IL-2 with a maximal rate coefficient β_T. This rate is brought to saturation by large numbers of A_T using a Michaelis–Menten saturation function, with Michaelis parameter g. Second, the inactivation of effector CTLs via their encounter with uninfected tumor cells (T_u), at a success rate coefficient − p₃. Third, the E_T natural death rate, with a rate coefficient μ_E. Finally, the diffusion factor of the effector cells reacting with tumor Ag cell population in the bladder geometry is assumed at rate coefficient D₆.

maximal growth rate coefficient (r), which is limited by the maximal tumor cell number (K). Second, bacterial infection, which is characterized by a coefficient rate of p₂. Third, capturing and elimination of T_u cells by APC cells (A), which were activated by tumor-Ag at a rate coefficient λ and to the activity of TAA-CTL effectors, (E_T), which destroy uninfected tumor cells, (T_u), at a rate coefficient α. The dependence in the equation of T_u on F_β is decreasing from 1 to a_Tβ with Michaelis constant e_T,β. And then there is a multiplication of those terms by an I₂-dependent Michaelis–Menten term, with Michaelis parameter g_I, to propose that in the absence of I₂, 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 _I2+gI by I_{2 Tu+gT}, to denote the inversely proportional reduction in effector cell acg_T tivity rate, such that when T_u = 0 the term is equal to 1 and when

In Eq. (7) ^∂I2(t,r) is the IL-2 dynamic. It is affected by the following five ∂t processes, with all processes assuming equal expression at a constant rate coefficient q₁. They reflect the IL-2 external source (i₂), which is injected into the bladder every θ time units. First, I₂ is consumed by APCs and CTLs. They assume that the rate of consumption is similar for both types of cells and denote its coefficient by q₂. The consumption depends on I₂ and is limited in a Michaelis–Menten fashion, with the Michaelis constant g_I. Second, introduction of − μ_I2 , the I₂ degradation rate coefficient. Finally, the diffusion factor of the IL-2 cell population in the bladder geometry is assumed at rate coefficient D₇.

In Eq. (8) ^∂Ti(t,r) is the dynamic of infected tumor cells, and it depends ∂t on three mechanisms. The first mechanism corresponds only to the rate of bacterial infection of uninfected tumor cells, (T_u), according to rate coefficient p₂. The second is the elimination of infected tumor cells (T_i) by their interaction with BCG-CTL effector cells (E_B), at a rate coefficient − p₄. Finally, the diffusion factor of the infected tumor cell population in the bladder geometry is assumed at rate coefficient D₈.

( ) lim_Tu→∞ ^gT = 0. Finally, the diffusion factor of the uninfected T_u+g_T tumor cell population in the bladder geometry is assumed at rate coefficient D₉. ∫ _R dF_β(t) (10) = α_βT T_u(t, r)dr − μ_βF_β(t) dt r₀

In Eq. (10) ^dFβ(t) _dt is the dynamic of a transforming growth factor-beta, proportional to the tumor cell population T_u with α_βT as a proportion coefficient and is destroyed at a rate of μ_β proportional to F_β. The dynamics of F_β are not geometry dependent but the dynamics of T_u(t, r) are geometry dependent. Therefore, the ∫^R _r0 T_u(t, r)dr provides all the uninfected tumor cells in the bladder’s geometry.

We assume the bladder’s geometry satisfies Eq. (11) as an approximation to the bladder’s geometrical configuration:

In Eq. (11), the variables x, y, z are the Cartesian coordinate system, r₀ and R are the radius of the internal and external spheres of the geometrical configuration, respectively. The bladder’s geometry is approximated using a perfect ring-sphere while the real human bladder is more like a ring-ellipsoid with three tunnels (Guzev et al., 2019). Fig. 1 visualizes the geometry of the system.

In Eq. (9) ^∂Tu(t,r) is the dynamic of uninfected tumor cells. It depends ∂t on four processes. First, the natural tumor growth characterized by a The boundary condition is based on two surfaces, the inner and outer sphere, respectively. The inner sphere boundary condition is assumed according to Eq. (12) known to be exactly b and decreases over time according to the system dynamics. The cells population associative to the immune system (e.g. A, A_T, A_β, E_β, E_T, I₂ ) is assumed to be equal to zero as the immune system does not allocate resources to the area until the BCG is injected. T_i is assumed to be equal to zero as well because cancer cells were not able to be infected by BCG before any BCG is injected into the system. T_u is assumed to be equally distributed inside each sphere approximating the cancer tumor. The inner sphere boundary condition is given to be:

The boundary condition of the external sphere is unknown. It is assumed that the natural cell population spread over time satisfies diffusion equations. Therefore, one can find the boundary condition of the external sphere by reverse engineering the values that best satisfy the known start conditions and internal boundary sphere conditions. Given the inner sphere boundary condition from Eq. (12) and the start condition from Eq. (14), algorithm (1) returns the outer sphere boundary condition.

Specifically, we used finite-elements in a 2d matrix where the geometry is represented by a tensor with 10⁶ values and the step in time (Δt) is 2 h (12 steps each day). Unless otherwise stated, it is assumed that the cancer polyps’ sizes are equal in the beginning of the treatment (t₀). In addition, the cancer cells population is equally distributed in the urothelium’s geometry.

The ^dFβ(t0) equation is the outcome of solving the linear ordinary Eq. dt (13) for F_β:

The initial condition is assumed to be:

B(r, t₀) = 0, A(r, t₀) = a, A_T(r, t₀) = 0, A_B(r, t₀) = 0,

where a > 0 is the natural influx of APC cells, n > 0 the number of polyps at the beginning of the treatment, Sp(α, θ, R) is a sphere with radius R and origin in angles (α, θ) on the (xy, xz) plain, respectively.

2.2. Numerical solution

To obtain a better understanding of how different parameter values influence the system dynamics, in this section we illustrate the behavior of the system using numerical analysis. To carry out the numerical simulations of the tumor-immune model, we used the parameter values from Table 3.

Eq. (1-10) are PDEs, second order, nonlinear, from R² to R¹⁰, where R² is the space of both time (marked by t) and radial distance from the center of the bladder’s geometry configuration (marked by r) and R¹⁰ is the population distribution of all nine populations (marked by B(t,r), A(t, r), A_T(t,r), A_B(t,r), E_B(t,r), E_T(t,r), I₂(t,r), T_i(t,r), T_u(t,r)) and the value of growth factor-beta F_β(t)). Galerkin-Petrov’s method is suitable for approximating the solution for such a system of equations (Skeel and Berzins, 1990).

In Eqs. (1)–(9) the leading (second order) factor is a diffusion dynamics and in Eq. (10) an ODE so the equations are elliptic. Therefore, Galerkin-Petrov’s method takes the form:

All the numerical calculations in this paper have been performed with Matlab software (version 2019b) using the pdepe function (Skeel and Berzins, 1990). However, the pdepe function has been modified to use specificity Eq. (15), provided with the system dynamics in Eq. (1-10), inner sphere boundary conditions from Eq. (12), and start conditions from Eq. (14).

The overall treatment has been divided into several components, following τ-long treatments as a result of the injection of BCG Σ^{N− 1} _m=0bδ(t − mτ) and IL-2 Σ^{N− 1} _m=0i₂δ(t − m) in Eqs. (1) and (7), respectively. The start condition of each τ-long treatment, except the first one, has been updated according to Eq. (16). Similarly, the boundary condition has been updated according to Eq. (17).

B(r, t_τ⋅i+Δt) = B(r, t_τ⋅i) + b, A(r, t_τ⋅i+Δt) = A(r, t_τ⋅i),

where i ∈ N is the i_th day of the overall treatment. The values of the model’s parameters are shown in Table 3 (see Appendix). The solutions of the system (1–14) with the parameter assumptions and values we used are shown in Fig. 2. Fig. 2 shows the cell population sizes over time of B(t,r), A(t,r), A_B(t,r), E_B(t,r), E_T(t,r), and T_i(t, r) where the x-axis in all nine graphs represents the time (in days) that has passed from the beginning of the treatment and the y-axis is the size of each cell population size, respectively. Fig. 3 shows the cell population of I₂(t,r), A_T(t, r), and T_u(t,r), where part of each graph is presented in a different scale.

The BCG (B) population reaches an upper limit, as shown in Fig. 2a. The maximum values of BCG (on days {7i}⁹ _i=1) occur on the days of the BCG injection, and during the week they decrease to a level of 2⋅ 10¹⁰. Fig. 3a shows the harmonic behavior of IL-2 (I₂), where i₂ is introduced every {7i}⁹ _i=1 day, which decreases to around 750 in the same day (as shown in Fig. 3a). In Fig. 3b nonactivated APC (A) cells significantly decrease until the third week, where there is converge to a harmonic oscillation between 5 and 10 (as shown in Fig. 3b).

In Fig. 2b it was shown that the population of APC activated by BCG (A_b) cells converges to the upper limit value of 1.4⋅10⁶. The same phenomenon occurs in CTL effector cells infected by BCG, as shown in Fig. 2d. Effector cells reacted to the Ag tumor grow in the first four weeks and then decrease. In addition, on the BCG injection days the cell population sharply increases, as shown in Fig. 2e. In Fig. 2c TAA-APC (A_T) cells growth during first four weeks. After the fourth week, A_T decreases over time, except for a local increase on the days when BCG is administered.

In Fig. 2f, the population of cancer cells infected by BCG increases during the first week. After the second injection of BCG into the bladder, the T_i population decreases over time with the local maximums on the BCG injection days, as shown in Fig. 2f. Additionally, in Fig. 3fc uninfected cancer cell populations (T_u) decrease over time, and in the first week this sharp decrease occurs with a constant rate (̃ 0.7). Each week, the population decreases in one factor of magnitude (̃ 10^{− 1}) as shown in Fig. 3c.

2.3. Solution stability

Lyapunov’s stability analysis method cannot be used for the system (1–14) because it does not satisfy the needed conditions (Buis, 1968). The system does not diverge to infinity on a representative case as presented in Fig. 2. Therefore, it is possible to analyze its stability for a given set of parameters using the system’s image space (Guzev et al., 2019).

Basically, the main goal is to find a treatment protocol resulting in a tumor-free equilibrium given a patient condition in the beginning of the treatment. We define a successful treatment as a treatment resulting in tumor-free equilibrium (T_u(t^*, r) = 0) and unsuccessful treatment otherwise, where t^* is the time at the end of the treatment.

In our model, four parameters affect the success of treatment: First, the initial cancer cell population size T_u(t₀). Second, the amount of BCG b injected over the course of the treatment. Third, the overall time of the treatment t in days; Fourth, the amount of IL-2 i₂ injected over the course of the treatment. Based on these, it is possible to define a four dimensional space to investigate the influence of T_u(0), b, t, and i₂ on the success of the treatment protocol. Determining which initial condition of the patient and which treatment protocol leads to successful treatment can be performed using the solution stability method described in (Lazebnik et al., 2020).

For our analysis we neglect IL-2 from the parameters and, therefore, are left with a three dimensional space defined by T_u(0), b, and t. We define a function S : R³→Z₂ such that S(T_u(t₀), b, t) ∈{0, 1} where 1 represents successful treatment and 0 represents unsuccessful treatment. For any vector v ∈(T_u(t₀),b,t), it is satisfied that S(v) is the solution for system (1–14) with the parameters from Table 3. The binary classification is determined according to a set of factors C = {c_i}⁹ _i=0. Where ∀i ∈ {0, …, 9} : c_i ∈ R⁺ are thresholds of the nine population sizes B(t), A(t, r), A_T(t), A_B(t), E_B(t), E_T(t), I₂(t), T_i(t), and T_u(t), respectively. It is assumed that there are lower and upper boundaries for each one of the parameters (T_u(t₀),b,t). All three parameters are lower-bounded by 0 as they cannot be negative. The upper boundary for T_u(t₀) is the number of all cells in the bladder. The amount of b injected is bounded by the bladder’s volume and the treatment time is bounded because medical treatment cannot be provided for eternity.

This results in space P⊂R³. P is a compact parameter’s set because it is complete (as a sub-set of R³) and bounded. We assume that the image of function S|_P is continuous and can be restored from discrete sampling.

Fig. 4 presents the space S|_P which has been sampled 16000 times, T_u(0) sampled 40 times ranging from 0 to 2.5⋅10⁸ in equal steps, b marked as BCG sampled 40 times ranging from 0 to 2.14⋅ 10⁶ in equal steps, and the treatment time in weeks t sampled 10 times ranging from 0 to 10 in equal steps. A treatment is considered successful if it satisfies T_u = 0 ∧BCG < 10⁸. The sampled space size is the largest that a modern personal computer was able to calculate in 8 h.

All the pixels in Fig. 4 have been computed successfully (no diversions at any T_u(0), b, t). Using these values, it is easy to see that by picking k = 1 then function S satisfies the Lipschitz continuity condition

d₁(S(v₁), S(v₂)) ≤k⋅d₃(v₁, v₂), where d_i is the Euclidean distance function in Rⁱ. Recall, we assume that S continues and can be restored from a discrete sample. Therefore, for some value ε > 0, there exists a unique solution to the initial value problem S(v) according to the Picard–Lindelöf theorem (Coddington and Levinson, 1955). Therefore, the system (1–14) is numerically stable on the sampled sub-space P.

3. Bifurcation analysis

Bifurcations in the cell population of A_B, A_T, T_i, T_u arise from changes in T_u(0), r, and b. Such bifurcations indicate various clinical results for different treatment protocols and allow to drawing the line between successful and unsuccessful treatment protocols.

The bifurcation numerically emerges in the cases where the sensitivity analysis (Section 7.2) shows behaviorally-different dynamics for different values of the parameter in question. We define two functions P₁ and P₂ behaviorally-different if there is no t_i ∈{0, …, t_max} such that

where P₁, P₂ are a polynomial interpolation of a cell population that originated from two different values and x is a manually picked threshold. Fig. 5 is the result of picking x = 1.

The motivation for Eq. (18) is to find an interval [t_i, t_max] such that the difference between curvatures of two functions is small enough (x) in time. If two functions satisfy this condition, then, there is some point in time t_i in the treatment protocol where the sum of differences in the changes of the cell populations are significant enough and the original functions should present an entirely different behavior.

Fig. 5a and b presents the bifurcation in the dynamics of T_u and A_b, respectively. In the case where T_u(0) = 2⋅10⁶ then the population of T_u is decreasing over time and A_b converge to 1.5⋅10⁶ after three weeks. On the other hand, where T_u(0) = 2.7⋅10⁸ the population of T_u increases over time and A_b has harmonic behavior oscillating between 7.5⋅10⁵ and 1.75⋅ 10⁶.

This bifurcation is associated with the fact that a small enough amount of cancer cells in the beginning of the treatment leads to relatively small amplitude in the BCG-activated APC A_B cell population size, after converging to 1.3⋅10⁶. The small amplitude reflects the overall immune system’s response; the A_B population is slightly affected by different amounts of injected BCG b as shown in Fig. 10. Similarly, a large amount (2.7 ⋅10⁸) of cancer cells in the beginning of the treatment leads to more sporadic behavior of the immune system.

Fig. 5c presents the bifurcation in the dynamics of T_i from the changes in the natural tumor growth r. In the case where r = 2.73⋅10^{− 1} the population of T_i increases in the first seven weeks and then decreases. On the other hand, where r = 2.92⋅10^{− 1} the population of T_i monotonically increases over time.

Fig. 5d and e presents the bifurcation in the dynamics of T_u and A_t from changes in the BCG installations b, respectively. In the case where b = 10⁵ the population of both A_t and T_u monotonically increases over time. On the other hand, where b = 2.5⋅10⁵ the populations of both T_u and A_t decrease to zero.

4. Time optimal treatment protocol

The question is whether it is possible to find the optimal treatment protocol in accordance with the initial conditions of the patient. By finding the equation describing the border between successful and unsuccessful treatment protocols in continuous settings, it can be determined if it is feasible to predict if a treatment will result in tumor-free equilibrium or not. In space P this is a phase transformation between unsuccessful and successful treatment in time as shown in Fig. 4.

Considering the initial state of the patient as the size of the cancer cell population at the beginning of treatment (T_u(t₀)), the treatment protocol as the amount of BCG injected (b) and the frequency of its introduction into the bladder (m), the optimal treatment ensures a minimum duration (if it exists) so that the treatment is successful. This phase transformation can be defined by a border function BF : R²→R such that BF(T_u(0),b)→t. Given (T_u(0),b, m = 7), any time t that satisfies t >= t_min will result in a successful treatment and t < t_min otherwise, where t_min = BF(T_u(0), b).

It is possible to approximate the border function using the border pixels and then using the least mean square (LMS) method (Björck, 1996) to approximate the border function itself. First, a pixel (i,j,k) will be defined as a border pixel iff it satisfies

where {0, 27} are the cases where all the values of 3 × 3 × 3 window centered in the pixel (i,j,k) are the same (either 0 or 1). Fig. 4 shows the border pixels in green. This method is inspired by computer vision threshold based edge detection algorithms (Yellasiri et al., 2011).

Second, to use the LMS method one needs to define the family function approximating the function. The function family has been chosen to balance between the accuracy of the sampled data on the one hand and simplicity of usage on the other (Shanock et al., 2010). The border function is obtained with a coefficient of determination R² = 0.85, using the LMS method. Therefore, it is safe to claim that function f(b, T_u(0)) is well fitting the data and presents a good approximation for the border function, which is easy and stable to compute.

f (b, T_u(0)) = a₁ + a₂b + a₃T_u(0) + a₄bT_u(0) + a₅b² + a₆T_u(0)²

For example, consider a patient who has a polyp in the urinary bladder of size T_u(t₀) = 5⋅10⁶ cells, and receives BCG treatment b = 2⋅ 10⁶ once a week. By setting these values in the Eq. (19), the minimal treatment time, 10 weeks, will be possible for successful results.

5. Treatment protocol based on initial tumor distribution

Weiner et al. (2019) have shown that the localization of a cancer polyp in the bladder affects the dynamics of cells inside the bladder as a result of differences in the biological reaction to treatment. The model proposed by Grabnar et al. (2006) is a mathematical model that describes the interaction of different layers of the urinary bladder from a biological point of view. They presented the variable drug concentration due to urine formation and voiding diffusion in the bladder tissue with parameters from in vitro experiments.

The model proposed in Eq. (1 - 14) takes into consideration the bladder’s geometry (Eq. (11)) and the diffusion dynamics of cell population biological reactions. Specifically, adding the diffusion dynamics to Eqs. (1)–(9) and the border condition from Eq. (12) and Algorithm 1. The model allows to fine-tune the prediction according to the tumor’s depth in the urothelium at the beginning of the treatment which in turn improves the accuracy of the model and therefore more accurately predicts the treatment result.

Cancer polyps depth and distribution inside the bladder’s geometry can be approximated using multiple ring-sphere shaped layers of the urothelium, centralized in the center of bladder. According to this approximation, it is possible to define a two-parameter space which describes non-isomorphic instances of tumor depths in the bladder’s geometry at time t₀ of the treatment. The normalized size of cancer cell ( ) population ^Tu(r,t) and the distribution of the population in the eighth ||T_u(r,t)|| layer of the urothelium.

Figs. 6 and 7 are derived using Eq. (15) as well, at the beginning of the treatment t₀, the geometry of the bladder as presented in Fig. 1, has been divided into eight separated geometries (layers). Each bladder tissue layer’s geometry satisfies the condition:

where i ∈[0, …, 7] is the index of the layer. Each geometry has been represented by a two-dimensional array (grid), where cancer cells have been allocated to a layer, with the values of the array assigned to be the amount of the uninfected cancer cells.

It is possible to analyze the differences in the dynamics between the layers of the urothelium, allowing us to better understand the differences in the dynamic in each tissue layer. Fig. 6 shows the dynamics of the system where the cancer cells T_u are equally distributed at a single layer of the urothelium. The population of BCG infected cancer cells T_i grows faster in the first 2 weeks as the cancer found only in a deeper layer but is reduced at a faster rate after the second week, as shown in Fig. 6a. Similarly, when the cancer cells are found deeper in the urothelium at the beginning of the treatment than the population of cancer cells T_u decays at a slower rate, as shown in Fig. 6b. In addition, the reaction of the immune system is relatively small, as shown in Fig. 6c and d.

	1 layer	2 layers	3 layers	4 layers	5 layers	6 layers
SP [m3t]	1.90 ⋅ 107	1.63 ⋅107	1.36 ⋅107	0.88 ⋅107	0.70 ⋅ 107	0.54 ⋅107
AP [m3t]	1.157 ⋅ 109	1.155 ⋅109	1.159 ⋅ 109	1.158 ⋅109	1.159 ⋅109	1.157 ⋅ 109

All cell population sizes for the different cases converge after 10 weeks, which can be associated with the fact that after enough time the diffusion of the cells’ populations are spread all over the geometry of the bladder and from that point operating as instance response, making the location insignificant, as shown in Fig. 6.

It is possible to divide the geometry to the eight layers of the urothelium to examine the influence of the initial cancer polyp depth on the layer level. Fig. 7 shows the dynamics of each individual layer of the system’s geometry where the cancer cells are equally distributed in the first (most shallow) layer of the urothelium.

The cancer starts to spread to other layers where layers closer to the polyps are affected faster. Nevertheless, in the first week of the treatment, the overall population of uninfected cancer cells (T_u) decreases, as shown in Fig. 7b. As a result, the BCG-infected cell population (T_i) shows a different increment rate in different layers until the end of the second week, where the changes between the different layers become negligible. Moreover, the first and eighth layers show a slightly different behavior in comparision to the other six layers. This phenomenon can be explained by the fact that these are the border layers of the model and have only one neighbor layer to spread to other layers which have two neighbor layers to spread to, as shown in Fig. 7a where the first and eighth layers show faster increase in the first three weeks of the treatment, reaching a higher maximum (1.2⋅10⁶) compared to the other six layers (1⋅ 10⁶). The populations of TAA-APC A_t and CTL reacting to tumor Ag. E_t cells show similar behavior, as shown in Fig. 7c and d, respectively. In both Figs. 6 and 7, the cell populations B, I₂, A, A_b, and E_B are not shown as the changes are neglected relative to the dynamics shown in Fig. 2.

From clinical trials, it is known that cancer cells usually spread across multiple layers of the urothelium (Weiner et al., 2019). To examine the influence of the number of tumor polyps on the system, multiple layers with tumors can be initialized. Table 1 shows the differences in the dynamics of the uninfected cancer cell T_u population size given a combination of {k}⁶ ₀ different layers where the cancer cells are equally distributed. SP is defined as the difference between the treatment protocol with the biggest population size and smallest population size [m³t]. This metric allows for measurement of the difference between the different treatment protocols of cancer cell population that originated in multiple layers of the urothelium, calculated as follows

where T^α _u(t, r) and T^β _u(t, r) are two different cancer cell populations dynamics from different initial conditions.

AP is defined as the average cancer population size for all the possible combinations. The values are calculated over the first four weeks of the treatment as follows

where P is the set of all the possible isomorphic permutations of k layers from the eight layers of the urothelium.

The differences between the case with the biggest population size of uninfected cancer cells at time t_max and with the smallest population size are reduced as the cancer is initialized in more layers. In addition, the average case for different amounts of initial layers that cancer cells are similarly distributed for all {k}⁶ ₀ layers, is shown in Table 1.

Based on Table 1 and Fig. 7b, it can be noticed that the treatment starts to be efficient after the BCG arrives at the deepest layer where the cancer polyp is located. In Fig. 7b, the amount of uninfected cancer cells T_u is almost equal after the first week across all the layers.

We argue that a treatment protocol should aim to arrive at the point where the BCG arrives at the deepest layer where the cancer polyp is located in the dynamics as early as possible. From Eq. (1) the diffusion of the BCG influences the distribution of the BCG in the geometry. The BCG diffusion is

Both r and D₁ are properties of the bladder and a way to influence them does not exist in the scope of this treatment. On the other hand, the size of the population B(r, t) can be changed in the treatment protocol by introducing into the system an increased amount of BCG. A larger amount of BCG results in faster spread of the BCG in the geometry of the bladder. Fig. 6b shows that the size of T_u in the whole geometry for the same amount of injected BCG b.

Layer	1st (baseline)	2nd	3rd	4th	5th	6th	7th	8th
BCG (b⋅ 106)	1.07	1.16	1.48	1.91	2.49	3.12	3.88	5.04

The motivation is to reduce the initial spread of the cancer cell population (T_u(0)) inside the geometry of the bladder relatively quickly. In doing so, the point can be reached where the system’s dynamic operates as an instant response and the geometry can be neglected (Lazebnik et al., 2020). Such an approach may be used to personalize the treatment protocol according to the patent’s initial spread of the cancer cell population at the beginning of the treatment, aiming to arrive at the point in the treatment where the treatment protocol can be replaced with a generic one with the best results.

Recall, we define a treatment protocol by the four parameter T_u(0), b, t, and i₂. We define two treatment protocols TP₁, TP₂ which differ only in the injected amount of BCG b to be T_u − similar after some time t^* if and only if the T_u(t^*) resulting from treatment protocol TP₁ (marked as T_u(t^*)|_TP1) and TP₂ (marked as T_u(t^*)|_TP2) satisfies

where k > 1 ∈ R. The motivation of the definition is to declare that the population of uninfected cancer cells T_u(t, r) in the whole geometry for two protocols is in one level of magnitude defined by some factor k.

Table 2 shows the amount of BCG that is needed in the first week such that a treatment protocol (TP) will be T_u − similar at time t^* = t₇ to the treatment protocol that satisfies T_u(0) = 1⋅10⁶, k = 10 and b = 1.07⋅10⁶, where t₇ is the seventh day of the treatment.

6. Conclusions

Mathematical modeling has already been shown to be a useful tool for studying the mechanism of tumor growth and response to therapy (Bunimovich-Mendrazitsky et al., 2015a; Bunimovich-Mendrazitsky and Goltser, 2011; Matzavinos et al., 2004; Bunimovich-Mendrazitsky et al., 2019). Models which better represent the biological and clinical dynamics and complexity can provide a better understanding of the system, resulting in more accurate prediction of an outcome of a treatment and determination of better therapeutic protocols (Shaikhet and Bunimovich-Mendrazitsky, 2018; Guzev et al., 2019). Specifically, models based on population analysis are a common way of describing such systems (Bunimovich-Mendrazitsky et al, 2015a; Kirschner and Panetta, 1998).

Based on the proposed model, the bifurcation for T_u, A_b, T_i, T_u, and A_t resulted in changes in T_u(0), r, and b as shown in Fig. 5. We argue that these bifurcations are the border line between a successful and unsuccessful treatment and the found values can be used to assist in clinical decisions based on the proposed treatment protocol.

In addition, a time optimal treatment protocol has been proposed using the stable matrix defined by the parameters that affect the success of a treatment based on the system’s image space (Guzev et al., 2019). A formula that, given the patient’s cancer cell population size at the beginning of the treatment T_u(0) and the amount of BCG installations b, returns the minimal treatment time in days for the treatment to be successful as shown in Eq. (19).

Furthermore, we show that the proposed model takes into consideration the initial distribution of the cancer cells over the geometry of the bladder and as such can provide more customized treatment by providing tumor polyps depth in the urothelium. Fig. 6 shows that a cancer tumor that originates in a shallower layer of the urothelium can be treated with less aggressive treatment either in treatment duration or injection of BCG b. Layer specific treatment is insignificant and in the case where the tumor polyps are spread out, it is easier to instead treat the case where the tumor is localized in a deep layer of the urothelium, as shown in Fig. 7 and Table 1.

Moreover, Table 2 shows the amount of BCG that is needed to be introduced in the first week of the treatment such that the model can neglect the geometry of the bladder without meaningful loss of accuracy. Providing a personal BCG injection protocol at the beginning of

Appendix

7.1 Computationally Parameter’s Values

treatment gives the best results.

These results are important for BCG immunotherapy, which modulates the healing effect. Understanding the key processes in tumor-immune interactions will be crucial for the development of effective treatments, for setting goals, as well as for optimizing dosage and schedule. The model presented here takes a step towards demonstrating the effect of the depth of the cancer and its spread, and further extensions of the model will be used to learn how to manage the treatment protocol for the successful elimination of bladder cancer.

Declarations of competing interest

None.

Table 3 describes the parameter values used in the calculation of the model not mentioned otherwise. All the values were taken from (Guzev et al., 2019). Parameters {D_i}⁹ _i=1 taken from (Lazebnik et al., 2020). D₁ satisfies the same conditions as Eq. (1) from (Lazebnik et al., 2020). D₂,D₃,D₄,D₅,D₆, and D₇ are equal to the diffusion factor of the effector cells E assuming the diffusion of the cell population related to the immune system is identical. D₈ and D₉ are identical to the diffusion coefficients of the BCG infected and uninfected cancer cells.

Parameter	Value	Source
βB	7.20⋅105	Guzev et al. (2019)
βT	7.50⋅103	Guzev et al. (2019)
Γ	4.70⋅103	Guzev et al. (2019)
R	8.50⋅10−3	Guzev et al. (2019)
B	1.07⋅106	Guzev et al. (2019)
Tu(0)	1.00⋅106	Guzev et al. (2019)
i2	1.00⋅106	Guzev et al. (2019)
i2	1.00⋅106	Guzev et al. (2019)
eTb	1.00⋅104	Guzev et al. (2019)
q1	7.00⋅10−3	Guzev et al. (2019)
q2	1.20⋅10−3	Guzev et al. (2019)
μA	3.80⋅10−3	Guzev et al. (2019)
μA1	4.00⋅10−1	Guzev et al. (2019)
μE1	1.90⋅10−1	Guzev et al. (2019)
μE2	3.40⋅10−3	Guzev et al. (2019)
А	3.70⋅10−6	Guzev et al. (2019)
αβT	1.38⋅10−4	Guzev et al. (2019)
αTβ	6.90⋅10−1	Guzev et al. (2019)
μB	1.50⋅10−1	Guzev et al. (2019)
μI2	1.15⋅101	Guzev et al. (2019)
μβ	1.66⋅102	Guzev et al. (2019)
gT	5.20⋅103	Guzev et al. (2019)
G	1.00⋅1013	Guzev et al. (2019)
gI	1.00⋅105	Guzev et al. (2019)
p1	1.25⋅10−4	Guzev et al. (2019)
p2	0.28⋅10−5	Guzev et al. (2019)
p3	1.03⋅10−10	Guzev et al. (2019)
p4	2.32⋅10−5	Guzev et al. (2019)
D1	1.00⋅10−4	Lazebnik et al. (2020)
D2	7.00⋅10−5	Lazebnik et al. (2020)
D3	7.00⋅10−5	Lazebnik et al. (2020)
D4	7.00⋅10−5	Lazebnik et al. (2020)
D5	7.00⋅10−5	Lazebnik et al. (2020)
D6	7.00⋅10−5	Lazebnik et al. (2020)
D7	7.00⋅10−5	Lazebnik et al. (2020)
D8	6.00⋅10−1	Lazebnik et al. (2020)
D9	6.00⋅10−1	Lazebnik et al. (2020)

7.2 Sensitivity Analysis

The sensitivity of the model to changes in several parameters has been explored. Using the sensitivity analysis it is possible to examine the lim itations and robustness of the model. The parameters in the model that have be explored are α, β, β_B, r, λ, b, and T_u(r, t₀) which are chosen because of their biological importance in the treatment protocol. It is known that each parameter has upper and lower boundary values in the scope of this treatment (Shaikhet and Bunimovich-Mendrazitsky, 2018). Figs. (8-12) present the different dynamics of the system for each one of the nine population sizes where the color of each plot defined the value of the considered parameter used in each calculation. The sample values for each parameter {

the time in days from the beginning of the treatment. The y-axis is the cell’s population size.

Fig. 8 shows the sensitivity of the system to parameter α ranging from 10³ to 9⋅10³ with step size 1.25⋅10³. Each color (and line style) represents a different sample. It is easy to see that there is no change in the population sizes as a function of α.

Fig. 9 shows the sensitivity of the system to parameter β ranging from 2.55⋅10^{− 2} to 4.25⋅10^{− 2} with step size 4.25⋅10^{− 3}. Sub figures A_t and A_b show that smaller β resulted in a stronger response of the immune system.

Fig. 10 shows the sensitivity of the system to parameter b ranging from 10⁵ to 10⁸ with step size 1.98⋅10⁶. Where b is low (blue line) then the tumor cells (T_i, T_u) are increasing and as a result the immune system cells A_t, E_t are increasing in the injection over time as well. It is easy to notice that a too little amount of b does not lead to tumour-free equilibrium. As b grows, more T_u cells are converted faster to T_i and the immune system’s cell populations grow but converge to some plateau. In addition, sub figures T_u and A_b show inherently different dynamics for different values of b which are further analyzed in the bifurcation section.

Fig. 11 shows the sensitivity of the system to parameter β_B ranging from 1.08⋅10⁴ to 1.81⋅10⁴ with step size 1.45⋅10³. Sub Figure A_t shows that A_t decreases as β_B increases while β_b has a minor effect on the other cell populations.

Fig. 12 shows the sensitivity of the system to parameter r ranging from 1⋅10^{− 3} to 5⋅10^{− 1} with step size 1.25⋅10^{− 2}. As r increases the needed time such that the cancer cell population size (T_i, T_u) decay rate is decreasing. Furthermore, for a large enough r the system does not converge to a tumor-free equilibrium, as shown in sub figures T_u and T_i. In addition, a small enough r results in a more stable immune system reaction A_b, A_t, A. Sub Figure T_u shows inherently different dynamics for different values of r which are further analyzed in the bifurcation section.

Article notes

Publication history

Received 21 August 2020 · Accepted 1 December 2020 · Published 10 December 2020

References

• Bhattacharya, S., Sah, P.P., Banerjee, A., Ray, S., 2020. Structural impact due to PPQEE deletion in multiple cancer associated protein - integrin V: an in silico exploration. ABiosystems 104216. link
• Björck, Å., 1996. Numerical Methods for Least Squares Problems. SIAM Journal on Scientific and Statistical Computing. Book OT51. link · link
• Bray, F., Ferlay, J., Soerjomataram, I., Siegel, R.L., Torre, L.A., Jemal, A., 2018. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA A Cancer J. Clin. 68 (6), 394–424. link · link
• Buis, R.G., 1968. Lyapunov stability for partial differential equations. NASA 1100. link · link
• Bunimovich-Mendrazitsky, S., Goltser, Y., 2011. Use of quasi-normal form to examine stability of tumor-free equilibrium in a mathematical model of BCG treatment of bladder cancer. Math. Biosci. Eng. 8, 529–547. link · link
• Bunimovich-Mendrazitsky, Pisarev, V., E. Kashdan, E., 2015a. Modeling and simulation of a low-grade urinary bladder carcinoma. Comput. Biol. Med. 58, 118–129. link · link
• Bunimovich-Mendrazitsky, S., Halachmi, S., Kronik, N., 2015b. Improving Bacillus Calmette Guérin (BCG) immunotherapy for bladder cancer by adding interleukin-2 (IL-2): a mathematical model. Math. Med. Biol. 159–188. link · link
• Bunimovich-Mendrazitsky, S., Kronik, N., Vainstein, V., 2019. Optimization of interferon-alpha and imatinib combination therapy for chronic myeloid leukemia: a modeling approach. Adv. Theor. Simul. 1800081. link · link
• Coddington, E.A., Levinson, N., 1955. Theory of Ordinary Differential Equations. New York McGraw-Hill. link · link
• Eikenberry, S., Thalhauser, C., Kuang, Y., 2009. Tumor-immune interaction, surgical treatment, and cancer recurrence in a mathematical model of melanoma. PLoS Comput. Biol., e1000362 link · link
• Eylert, M., Hounsome, L., Persad, R., Bahl, A., Jefferies, E., Verne, J., Mostafid, H., 2014. Falling bladder cancer incidence from 1990 to 2009 is not producing universal mortality improvements. J. Clin. Urol. 7, 90–98. link · link
• Fridman, A., Kao, C.Y., 2014. Mathematical Modeling of Biological Processs, Lecture Notes on Mathematical Modeling in the Life Sciences. Springer, Cham. link · link
• Grabnar, I., Bogataj, M., Belic, A., Logar, V., Karba, R., Mrhar, A., 2006. Kinetic model of drug distribution in the urinary bladder wall following intravesical instillation. Int. J. Pharma. 52–59. link · link
• Guzev, E., Halachmi, S., Bunimovich-Mendrazitsky, S., 2019. Additional extension of the mathematical model for BCG immunotherapy of bladder cancer and its validation by auxiliary tool. Int. J. Nonlinear Sci. Numer. Stimul. 20 (6), 675–689. link · link
• Herr, H.W., Laudone, V.P., Badalament, R.A., Oettgen, H.F., Sogani, P.C., Freedman, B. D., Melamed, M.R., Whitmore, W.F., 1988. Bacillus Calmette-Guérin therapy alters the progression of superficial bladder cancer. J. Clin. Oncol. 1450–1455. link · link
• Hornberg, J.J., Bruggeman, F.J., Westerhoff, H.V., Lankelma, J., 2006. Cancer: a systems biology disease. Biosystems 81–90. link
• Jordão, G., Tavares, J.N., 2017. Mathematical models in cancer therapy. Biosystems 12–23. link · link
• Kirschner, D., Panetta, J.C., 1998. Modeling immunotherapy of the tumor–immune interaction. J. Math. Biol. 37, 235–252. link · link
• Lamm, D.L., 2006. Improving patient outcomes: optimal BCG treatment regimen to prevent progression in superficial bladder cancer. Eur. Urol. Suppl. 5, 654–659. link · link
• Lazebnik, T., Yanetz, S., Bunimovich-Mendrazitsky, S., Haroni, N., 2020. Treatment of bladder cancer using BCG immunotherapy: PDE modeling. Partial Differ. Equ. doi:10.26351/FDE/26/3-4/5 · link
• Matzavinos, A., Chaplain, M.A., Kuznetsov, V.A., 2004. Mathematical Modelling of the Spatio-Temporal Response of Cytotoxic T-Lymphocytes to a Solid Tumour, Mathematical Medicine and Biology, pp. 1–34. doi:10.26351/FDE/26/3-4/5 · link
• Morales, A., Eidinger, D., Bruce, A.W., 1976. Intracavity Bacillus Calmette-Guérin in the treatment of superficial bladder tumors. J. Urol. 116, 180–183. link · link
• Paterson, D.L., Patel, A., 1998. Bactillus calmette-guerin (BCG) immunotherapy for bladder cancer: reivew of complications and their treatment. Aust. N. Z. J. Surg. 340–344. link · link
• Redelman-Sidi, G., Glickman, M., Bochner, B., 2014. The mechanism of action of BCG therapy for bladder cancer–a current perspective. Nat. Rev. Urol. 11, 153–162. link · link
• Shaikhet, L., Bunimovich-Mendrazitsky, S., 2018. Stability analysis of delayed immune response BCG infection in bladder cancer treatment model by stochastic perturbations. Comput. Math. Methods Med. doi:10.1155/2018/9653873 · link
• Shanock, L.R., Baran, B.E., Gentry, W.A., Pattison, S.C., Heggestad, E.D., 2010. Polynomial regression with response surface analysis: a powerful approach for examining moderation and overcoming limitations of difference scores. J. Bus. Psychol. 543–554. doi:10.1155/2018/9653873 · link
• Simon, M.P., O’Donnell, M.A., Griffith, T.S., 2008. Role of neutrophils in BCG immunotherapy for bladder cancer. Urol. Oncol.: Semin. Orig. Invest. 341–345. link · link
• Skeel, R.D., Berzins, M., 1990. A method for the spatial discretization of parabolic equations in one space variable. SIAM J. Sci. Stat. Comput. 11, 1–32. link · link
• Wei, H.C., 2016. Polynomial Regression with Response Surface Analysis: A Powerful Approach for Examining Moderation and Overcoming Limitations of Difference Scores, Discrete and Continuous Dynamical Systems - Series B., pp. 1279–1295. link · link
• Weiner, A.B., Desai, A.S., Meeks, J.J., 2019. Tumor location may predict adverse pathology and survival following definitive treatment for bladder cancer: a national cohort study. Eur. Urol. Oncol. 2 (Issue 3), 304–310. link · link
• Yellasiri, R., Poornima, B., Sridevi, T., 2011. Threshold based edge detection algorithm. Int. J. Eng. Technol. 3. link · link