A MODEL OF THERMOTHERAPY TREATMENT FOR BLADDER CANCER

. In this work, we investigate chemo- thermotherapy, a recently clinically-approved post-surgery treatment of non muscle invasive urothelial bladder carcinoma. We developed a mathematical model and numerically sim-ulated the physical processes related to this treatment. The model is based on the conductive Maxwell’s equations used to simulate the therapy administration and Convection-Diﬀusion equation for incompressible ﬂuid to study heat propagation through the bladder tissue. The model parameters correspond to the data provided by the thermotherapy device manufacturer. We base our computational domain on a CT image of a human bladder. Our numerical simulations can be applied to further research on the eﬀects of chemo- thermotherapy on bladder and surrounding tissues and for treatment person-alization in order to maximize the eﬀect of the therapy while avoiding burning of the bladder.


Introduction.
1.1. Bladder Cancer. Urothelial Bladder Carcinoma or Bladder Cancer (BC) represents an increasing health problem worldwide. It is estimated that around 400,000 new cases are diagnosed annually and 150,000 people die directly from BC each year. The highest incidence of BC occurs in the industrialized and developed countries in Europe, North America, and Northern Africa. According to the current statistics urinary bladder carcinoma is the fourth most common new cancer in men and ninth in women [15]. BC is a sporadic disease and rarely inherited genetically, around 70% of the cases could be related to exposure to carcinogens, especially aromatic amines. In the past occupational exposure was common, however today most of the cases are related to cigarette smoking. The majority (∼ 85%) of new patients are to be diagnosed with the low grade non muscle invasive BC (NMIBC), and following local resection with or without complementary intravesical treatment they will need a lifelong follow-up to rule out tumor recurrence and progression [7,11].
NMIBC is an incurable disease and it may be considered a chronic disease, due to its natural history and tendency to recur in most of the patients. Due to the risk of progression many of the patients with NMIBC are offered radical cystectomy. A successful addition of a treatment modality is another step in order to preserve the native bladder despite aggressive NMIBC. The primary goal of such a treatment, especially for those at an increased risk, is to prevent recurrences and progression

CHRISTOPH SADÉE AND EUGENE KASHDAN
Intravesical BCG (Bacillus Calmette-Guérin) immunotherapy is considered the gold standard adjuvant treatment for the prevention of BC recurrence, although very effective BCG has numerous side effects some serious and some even life threatening. Unfortunately, all other topical treatments developed to reduce the threat of side effects were not effective enough indicating the need for alternative measures [18].
1.2. Chemo-thermotherapy. The usual procedure for NMIBC is local removal of the tumor with a cystoscope and follow up chemotherapy. Unfortunately as shown in the study [9], the recurrence of bladder cancer after local removal and chemotherapy is rather high. Several devices have been invented to further reduce this recurrence rate by heating the chemotherapy drug solution, which is circulated within the bladder cavity or by heating the bladder wall itself, therefore increasing the diffusivity of the chemotherapy drug into the bladder tissue and/ or inducing hyperthermia in the bladder wall.
One of the most successful chemo-thermotherapy treatments that has been tested in vivo since mid-1990s is a microwave induced chemo-hyperthermia (C-HT) [18]. The C-HT device developed by the compnay called Synergo© [3] has recently passed clinical trials and has been cleared for clinical use in the US and in Europe. However, it is not often used by physicians due to its unpredictability. The optimal treatment schedule for C-HT is unknown, and current protocols have been based largely on existing chemotherapy trials.
During C-HT treatments the tissue surrounding bladder can be heated up to 44 o − 45 o C, which is the maximal tolerated temperature that may induce biological response without tissue damage [18]. Serious burn injuries to the bladder were not reported in the publications based on the analysis of the clinical trials of the Synergo© device, but posterior wall thermal reaction (PWTR) was commonly seen [17]. At cystoscopy, PWTR appears as a small, superficial, darkly discolored patch surrounded by hyperemia in the location where the heating device contacted the bladder wall during treatment. The overall rate of PWTR was 40.2% [12,10] and most cases were mild, however, few cases of severe and prolonged (asymptomatic) PTWR were reported [12].
On the other hand, it was reported [12] that no significant difference in the extent of the tumor cell killing was observed compared to the regular chemotherapy when the tissue was heated below 37.8 o C. Intravesical administration of chemotherapy delivers the drug to the bladder tissue by diffusion. Although diffusion is very common in biological systems, it is also a very slow process (on a macroscopic scale).
The most commonly administered chemotherapy drug for NMIBC treatment is Mitomycin C (MMC), which has low toxicity and is considered safe. In regular body temperature, MMC penetration to the bladder wall is time dependent. In 8.3 min, MMC will penetrate to a depth of 1 mm and to achieve penetration of 1 cm at least 14 hours of exposure are needed [22]. Looking at the anatomy of the bladder wall (see, for instance, [19]) one can notice that the width of the mucosal layer is between 50 − 200µm, the lamina propria width is between 200 − 1000µm, hence the detrusor muscle is located around 200−1200µm from the bladder lumen. In [22], the authors showed that there is a considerable drop of MMC concentration relative to the bladder cavity. In accordance with the clinical protocol, the duration of routine instillation of MMC patient are instructed to hold the drug in their bladder for 60-90 minutes which is insufficient for adequate full bladder width penetration. The diffusion of MMC into the bladder wall is also impaired by other factors such as dilution. MMC urine concentration is diluted to about 1:10 at the end of 1 hour instillation by constantly produced urine [6].
In [12], the authors reported on cell survival of four human bladder cancer cell lines that were exposed to several concentrations of MMC heated to 37 o C or 43 o C for 60 min. They found that thermotherapy alone did not increase cell death rate; however, combined C-HT with MMC improved the drug effect in all cell lines. The authors hypothesized that the synergistic effect may be due to higher intracellular MMC-uptake and increased activation of MMC.
In this work, we provide, first to our knowledge, detailed analysis of major physical processes behind chemo-thermotherapy. Based on the data available from the device manufacturer and imaging of the patient bladder we simulate chemothermotherapy treatment with an emphasis on potential bladder overheating and a view on the C-HT protocol optimization and treatment personalization.
The manuscript is organized as follows: In section two we describe the C-HT procedure and its physical principles, including parameters and characteristics of the the Synergo© device. Section three is dedicated to the mathematical model of the treatment based on the coupled Maxwell's and Convection-Diffusion equations. In Section four we present the results of our numerical simulations. Conclusions and future work are outlined in the Section 5. The details of the numerical schemes used in our simulations are summarized in the Appendix.
2.1. Human bladder. The CT image of a human bladder (Fig. 1) shows its almost symmetric form. The urinary bladder is located within the pelvis and it consists of the following layers: bladder lumen, urothelium, lamina propria, muscle and fat [16,19] as sketched in Fig. 2. The urothelium is a highly specialized layer of epithelial cells lining the bladder. It has to maintain a tight barrier against urea and other toxins, while accommodating large changes in bladder volume. The epithelial cells of the urothelium form part of an integrated network, which plays a major role in bladder sensory system. It is responsible for local circulation of blood and removal of pathogens. Under normal conditions the urothelium has a very low cell turnover rate. However, when the urothelium is damaged it is able to repair itself rapidly [16,23].
Majority of NMIBC (∼ 95%) originate in urothelium and it suffers a major damage following tumor removal. With high probability, the urothelial cells damaged during the surgery and the post-surgery treatment will be replaced by the cancerous cells. New polyps will continue to grow, and thus the BC recurrence is a very common follow-up scenario.
Therefore the administration of C-HT should reflect the current state of the urothelium, and the protocol of the treatment should be optimized to avoid additional damage to this tissue.

2.2.
Principles of the C-HT administering device. The Synergo© device model SB-TS 101 [1] approved for clinical use is a computer controlled, intravesical irrigation system combined with an energy delivering unit. The scheme of the device is presented in Fig. 3. It has a 915 MHz microwave applicator contained within a catheter which is induced into the bladder through the urethra. A specialized urinary catheter is used during thermotherapy treatment and it transmits microwaves with the aim to heat the bladder to 42 o C ± 2 o C causing hyperthermia in the bladder wall. The catheters position is checked with ultrasound and then secured in the bladder with a water inflatable balloon as seen in figure 4. There are a total of six sensors (thermocouples), measuring the superficial temperature of the bladder wall [3], these are the only means to determine the temperature within the bladder by the device. Chemotherapy drug Mytomicin C (MMC) dissolved in distilled water is continuously pumped out and into the bladder through a small outlet and inlet in the catheter. An external degassing system removes any air bubbles within the fluid and an external cooling system controls the temperature. The treatment is applied for up to one hour (usually 2 consecutive 30 min sessions) while the patient has local anaesthesia.  3. Mathematical model. In this Section, we present a theoretical model of the SB-TS 101 device and its effect -the induced hyperthermia to the bladder wall. We choose the 2D slice of bladder as our computational domain. The rationale behind this choice will be discussed in more detail in Section 4. We also consider bladder to be quasi-symmetric (see Fig. 4) with the symmetry axis z running along the catheter. This assumption allows us to simplify the model while preserving its validity and it also makes our approach applicable for future "real-time" treatment adjustment and optimization.
The microwave applicator is modeled as a dipole antenna radiating electromagnetic waves in all azimuthal directions perpendicular to its axis, surrounded by a plastic casing with negligible absorption. The electric and magnetic field just outside of the casing is then given in spherical coordinates in accordance with [20]: where I 0 is electric current, L is a length of the antenna, r case is a distance between antenna and plastic casing, c is a speed of light and ω represents the angular frequency (equal to 2π · 915 MHz). Here E and H correspond to the components of electric and magnetic fields respectively. The waves radiated by the dipole antenna form the 3D shape of a torus. To incorporate them in our 2D model, a slice of the torus was taken where the magnitude E θ reaches its maximum. This occurs at θ = π /2, which is at the center of the dipole.
To analyze the propagation of the electric field outside the catheter casing one has to consider loss of the electric field strength due to absorption. The propagation of microwaves through lossy media is given by Maxwell's curl equations that we need to solve numerically: The electromagnetic field pattern consists of a transverse electromagnetic (TEM) wave. That could be separated into TE (transverse electric) and TM (transverse magnetic) modes. In accordance with the principles of microwave heating, in our model we consider a TM-mode, which involves electric field E z propagating along z-axis (catheter) and magnetic field perpendicular to it. Here ε is a real part of the complex dielectric permittivity and σ is electric conductivity.
In a 2D TM cylindrical representation the equations (5)-(6) become: Equation (5) does not contain a magnetic field loss term since any human tissue has a negligible magnetic field absorption coefficient σ m . The second term on the RHS of equation (6) is the electric loss term, showing how much is absorbed by the different tissues, depending on their electric conductivity. The full dielectric permittivity is a complex number given by ε = ε + iε where ε is the loss factor, which in turn is related to the electric conductivity: ε = σ /ω. All parameters are space dependent due to different dielectric properties of bladder tissue, water and fat. The electric conductivity is assumed to be a constant for each tissue type, neglecting any changes due to temperature variations. Equation (7) -(9) are numerically solved using the time-domain finite difference scheme [21]. The details of the discretization are summarized in the Appendix.
The conversion of the microwave energy into heat in tissue and the diffusion and convection of heat through the bladder cavity and tissues are modeled using the Convection-Diffusion equation for an incompressible fluid (chemotherapy drug solution inside the bladder): where u is a temperature, v is a velocity vector and D = kt /ρcp is a thermal diffusion coefficient that depends on heat capacity c p , thermal conductivity k t and material density ρ. The term S in equation (10) stands for the source and it corresponds to the heat delivered to the tissue by microwaves: Here: ω = 2πf is angular frequency (f = 915 MHz) ε is dielectric loss factor and |E| stays for the amplitude of electric field. The source term connects electromagnetic and fluid dynamics components of the model. It is obtained by integrating the Poynting vector, giving the power flow through a closed surface and dividing by the specific heat capacity and thermal conductivity (see [13] for the derivation).
For computational purposes, equation 10 is expressed in 2D polar coordinates (r, φ): The equation is solved numerically using the Crank-Nicolson finite difference scheme as summarized in the Appendix.

Computational experiments.
4.1. Assumptions and parameters. We consider a 2D slice of the bladder shown in Fig. 1. The model equations are derived for the full 3D case, however, we can exploit quasi-symmetric form of bladder and reduce the computations. This also makes our model well suitable for embedding in the treatment optimization algorithm and enables its use in a "real-time" treatment adjustment.
Our computational domain represents a bladder with the urothhelium fully recovered from the surgery, which is a very natural assumption taking into account that urothelium has a very rapid damage repair rate and considering the requirements of the clinical protocol to start a treatment several weeks after the surgery.
The dielectric properties of urothelium and surrounding tissues (conductivity and permittivity) are temerature dependent. However, this dependence is negligible for the range of the temperatures we are working in. The chemotherapy drug MMC is dispersed in water. However, its concentration is too small to affect the conductivity. The dielectric permittivity of tissues is considered a constant, in particular, due to the missing experimental information. Our computational experiments (not shown here) have also confirmed that slight variation in the values of muscle and fat permittivity do not affect the results of numerical simulations.
The manufacturer of the system did not provide details about the antenna used in SB-TS 101 except its resonant frequency equal to 915 Mhz, so we considered a standard dipole antenna in our simulations. The operational current I 0 applied to the antenna and required for equations (1)-(2) was not provided as well. A current was chosen in such a way that the surface temperature of the bladder wall would have an equilibrium temperature of 42 o C±2 o C with the applied cooling and in the ideal case of deionized water with MMC in the bladder cavity as recommended by the manufacturer. A Gaussian pulse was chosen as a source based on [1].
In order to test our model two numerical simulations with different electric permittivity values of the solution in the bladder cavity are carried out. One with the electric permittivity values for deionized water ε = 75 − 1i and one for salty water with 10 parts per thousand (ppt) ε = 75 − 35i. Both values where assumed to be constant over the applied temperature change, taken at regular body temperature from [2]. The initial temperature was set equal to 36.6 o C and all other temperatures were computed in relation to it. The cooling applied during the simulation is set to be at the initial temperature as well.
The physical size of the urothelium is estimated using the image in Fig. 1 that has a measurement scale. The chosen spatial discretization was 100 r-steps and 100 θ-steps as shown in Fig. 5 where the different colors correspond to different values of electric permittivity, blue for fat tissue, orange for bladder tissue and yellow for the solution.
The time step in Yee algorithm is chosen below the CFL condition to ensure numerical stability. The boundary conditions implemented in both Yee and Crank-Nicolson schemes are discussed in the Appendix. The results of both simulations show that the system has reached equilibrium after 10min. Hence we are showing the temperature distribution at any later time.

4.2.
Simulations. Our first computational experiment shows a situation when the bladder is emptied from urine and filled with the deionized water plus chemotherapy drug following the clinical protocol. The isoclines in Fig. 6 correspond to the areas of the constant temperatures. In this case the heat induced by antenna propagates beyond the bladder wall and forms the areas of the high temperature that might potentially lead to the local burning of bladder and surrounding tissues or another damage due to the overheating. The maximum temperature exceeds 10 degrees above the initial (body temperature) and reaches 47 o C near the upper right corner of the bladder lumen. Note that the thermocouples embedded in the device would only measure a temperature of ≈ 41 o C at the surface of the bladder wall, hence not registering the local hot spot in the top right corner.
Another simulation considers the bladder cavity filled with slightly salty water. This situation might correspond to the low amount of urine remaining within the Temperature difference Figure 6. C-HT treatment using deionized Water bladder and also that urine will constantly accumulate within the bladder cavity. The results are shown in Fig. 7. One can observe from Fig. 7 that the outcome of the simulation is drastically different compared to Fig. 6: microwaves are getting absorbed by the water and the bladder wall is heated homogeneously as a result of convection. The temperatures within the bladder seems to get very high, but during a real treatment that would not be the case since the cooling and power output of the antenna would be adjusted to such a level that the superficial bladder wall would only be at 42 o C ±2 o C. Hence avoiding burning but also only allowing superficial temperature penetration. This is qualitatively the same as heating the

5.
Conclusions. In this work we developed a mathematical model and presented the first quantitative analysis of the chemo-thermotherapy. Our results show that the current recommended administration of C-HT can lead to local hot spots and burning some areas of the bladder and surrounding tissue as seen in Fig. 6. Incorrect administration of C-HT with residual urine or a slightly salty solution within the bladder cavity, will only give ineffective superficial heating of the bladder wall, not resulting in the desired effect of deep heat penetration. As mentioned previously this is qualitatively equivalent to heating the fluid externally and the circulating it within the bladder. To avoid ineffective heating and local hot spots one should continuously check for the salinity of the solution. Small amounts of salinity (e.g. 0.1 ppt) might be beneficial in damping excess heat received by some areas of the bladder and more evenly distributing the energy within the bladder and surrounding tissue.
Our computatinal experiments essentially show that the treatment is dependent on the loss factor ε of the fluid within the bladder. Another application of our work is a personalized approach to the chemo-thermotherapy. We have shown that the effect of C-HT and the protocol of administration can be analyzed based on patient-specific CT images of the bladder and the patients individual core temperature. The model can be adjusted to be based on Ultrasound images, which would allow for a continuous monitoring of the bladder while the treatment is performed.
Thermotherapy is a rather old form of cancer therapy, but due to its limitation in clinical applications and its unpredictability has not spread as widely as one would expect. Only marginal success was reported in clinical trials involving thermotherapy alone. However, with the new types of thermal delivery systems such as the SB-TS 101 device and more accurate models, advances can be made for more precise heat and drug delivery and better control mechanisms [18]. Successful administration of chemo-thermotherapy to treat BC, which has a high potential of recurrence, gives hope that these new types of therapeutic approaches may one day be extended to the treatment of a broader spectrum of tumors.

Appendix.
Discretisation of Maxwell's equations. To analyze the propagation and absorption of microwaves through the bladder cavity and tissue, the Yee finite difference scheme was derived for equations (7) -(9) in accordance with [21]. The algorithm sets off each field component by half steps in space and the magnetic field from the electric field by half a time step.
The electric loss term in equation (15) causes a problem since its time index is located at a half step. This can be resolved by averaging over the previous and subsequent time step: Subbing into (15) and rearranging for E z | n+1 i,j yields: The inside boundary where the antenna is surrounded by water is described as a Dirichlet boundary. The magnitude of equation (2) is combined with a Gaussian pulse to simulate one oscillation of the wave.
The amplitude of the electromagnetic waves decays after it goes into the fat and muscle layers of bladder. Therefore its naturally to truncate the computational domain using absorbing boundary conditions [21]. In our simulations we used the 1st order Mur conditions on the outer boundary: Discretization of Convection-Diffusion equation. To simulate the diffusion of heat in the bladder wall and the cooling by convection a Crank Nicolson finite difference scheme was constructed for equation (10) [8]. In the scheme each spatial derivative is averaged over the n th and n+1 th time step, where as the time derivative is discretized as usual Note that in equation (20), τ was chosen as discrete time step label instead of ∆t as it is not the same time interval as in equations (13) - (19). Subbing equation (20), (21) and (22) into equation (12) and arranging all n + 1 terms on the LHS and all n terms on the RHS yields: where a = τ Equation (23) is implicit and can be written in the form: A and B contain all the coefficients from α i,j to ζ i,j that are in front of u n+1 and u n respectively. Equation (24) represents a linear system and since all terms on the RHS are known, it can be solved by matrix inversion: The inside boundary where the antenna is located and the solution is constantly pumped in and out of the bladder cavity is defined by a Dirichlet boundary conditions. Its temperature is given by the temperature of the fluid after it has been cooled externally and it is about to enter the bladder cavity again.
The outside boundary is slightly more complex. The bladder is not perfectly round hence a 2D polar finite difference scheme is not optimal. For once the finite difference scheme will make the bladder wall boundary edged, this becomes less important when making the grid cells sufficiently small. Another issue are that the outside boundary where the bladder stops and the surrounding tissue begins is located at different values of r and θ which would result in some coefficients in A being rearranged into B since their u n+ 1 values are known. To circumvent this, a type of Robin boundary condition was constructed.
We considered Fourier's Law where q = heat flux k = thermal conductivity of surrounding tissue (fat) A = surface area and adapted the thermal contact resistance formula from [14]: where u body is body temperature and s j is a position in bladder surrounding tissue where temperature equals body temperature. The distance between the heated bladder wall and the point where the tissue has cooled down to body temperature again is estimated to be 5cm. Combining equations (27) and (28) and rearranging for u N +1 yields: Subbing (29) into equation (23) at u n+1 gives (α i,j + γ i,j ) u n+1 N −1,j + (β i,j − 2h r L j γ i,j ) u n+1 N,j + δ i,j u n+1 N,j+1 + 2h r L j γ i,j u body = (−α i,j − γ i,j ) u n N −1,j + (−β i,j + 2 + 2h r L j γ i,j ) u n N,j − δ i,j u n+1 N,j+1 − 2h r L j γ i,j u body .
Simulation parameters. Thermal parameters are very hard to find for the individual layers of the bladder such as the urothelium or lamina propria. Instead general values for the bladder and fat tissue were used. All values originate from the same source [4] to allow for consistency in measurements. To obtain a proper velocity field one would have to solve fluid dynamics equations in the bladder and see how the velocity field will distribute. In practice, the velocity will become very low when approaching the bladder wall but will be at its peak when leaving and entering the bladder through the catheter. We consider water flowing outwards in all directions radially from the center. The average value of velocity that we have chosen (0.0005m/s) most closely resembles an overall heating of bladder wall up to 5.4 o C.
The dielectric parameters are based on the Cole-Cole dielectric model of bladder as presented in [5]: