Study of Tumor Growth under Hyperthermia Condition

Abstract

The new concept of keeping primary tumor under control in situ to suppress distant foci sheds light on the treatment of metastatic tumor. Hyperthermia is considered as one of the means for controlling tumor growth. To simulate the tumor growth, a continuum mathematical model has been introduced. The newest understanding of the Warburg effect on the cellular metabolism and diffusion of the nutrients in the tissue has been taken into consideration. The numerical results are compared with the in vivo experimental data by fitting the tumor cell doubling time/tumor cell growth rate under different thermal conditions. Both the tumor growth curve and corresponding average glucose concentration have been predicted. The numerical results have quantitatively illustrated the controlling effect on tumor growth under hyperthermia condition in the initial stage.

1. Introduction

Cancer is the second major cause of human death in the world, and its mortality rate is growing every year [1]. Treatments include surgery, radiotherapy, chemotherapy, and gene therapy. Thermal therapy has also been intended to locally destroy tumor cells or enhance the body defense against tumor cells. However, recurrent rate of malignant tumor is still high [2], and the efficacy of the existing therapeutic means is yet to be improved. A new concept has been proposed recently that the primary tumor suppresses distal foci [3, 4]. This sheds new light on tumor treatment. Keeping the primary tumor in situ but restricting its size might enable the host to impede the development of distal foci and progression of metastasis.

For tumor growth, there are three distinct stages: avascular, vascular, and metastatic/invade stage. Mathematical models have been developed to perform parametric studies on factors influencing tumor growth or to evaluate the outcome of tumor treatment modalities [5, 6]. Model-based numerical studies would enable one to extrapolate more spatial and temporal information from the experimental findings and to make predictions [7]. Laird [8] first found that the tumor growth data-fitted Gompertz function could be used to simulate the entire growth curve, which was defined as an empirical model. Hu and Ruan [9] studied the suppression effect of immunity system on tumor growth by merging the Gompertz function into a cellular automaton model. Other mathematical models based on certain biological assumptions have also been attempted to predict tumor growth curve using fundamental physics, such as mass/energy conservation. Greenspan [10] introduced surface tension into the diffusion model developed by Burton [11]. Tumor growth/inhibition factors [12, 13], cell adhesions [14, 15], angiogenesis [16, 17] and invasion [18, 19] were further considered to describe tumor growth at different stages.

Models focusing on the avascular stage [20–27] have been well studied and could be easily applied to in vitro experiment. Ward and King [23, 24] and Casciari et al. [28] proposed a continuum mathematical model focusing on how nutrients’ concentration affects tumor growth. These models typically consist of reaction-diffusion equations. Forbes [29] further incorporated energy metabolism (ATP production rate) into the growth model. However, most of these models have not taken the Warburg effect into consideration, which fundamentally differentiates the tumor cell metabolism from that of the normal cells.

In 1930, Warburg (1930) proposed that tumor cells preferentially underwent glycolysis when consuming glucose even under aerobic conditions. Unregulated glucose uptake and lactic acid production have been found in tumor cells as compared to normal cells [30, 31]. It indicates that tumor cells obtain energy to maintain their viability primarily relying on anaerobic metabolism. This phenomenon was termed as “the Warburg effect.” Anaerobic glycolysis consumes one molecule of glucose to produce 2 molecules of ATP as compared with oxidative phosphorylation which can produce 38 molecules of ATP [31–40]. Although the latter is much more efficient in glucose utilization, the rate of anaerobic glycolysis is much faster than aerobic metabolism. Therefore, the inefficient metabolism pathway might still supply enough energy for tumor cells to maintain their activities and differentiate at the cost of unreasonable consumption of glucose. The mechanisms causing the Warburg effect have been explained by gene mutation [38], signaling pathway alternations, possible defects in mitochondria [36, 41], and microenvironment deterioration (hypoxia or fluctuation of oxygen) [34, 37, 42]. Heiden et al. [32] have reported that biomass synthesis in tumor cells plays a role in the Warburg effect. Furthermore, he has determined nutrition utilizations in tumor cells: 85% of glucose converting to lactate in cytoplasm, 5% reacting in mitochondria, and 10% synthesizing biomass. As the metabolic activities greatly influence the growth of tumor, it is necessary to include this unique metabolic mode of tumor in mathematical models.

Although thermal treatment has been applied in clinical applications for many years, most of them were used as short-term treatments. There are three classes of the treatment strategies [43–45]: with a mild temperature at 40~41°C for 6~72 hours until the thermal dose is equivalent to 5 minutes at 43°C; a moderate temperature at 42~45°C for 15~60 minutes; and high temperature >50°C for 4~6 minutes. Both mild and moderate temperature treatments (40~45°C) termed as hyperthermia treatment impair mammalian cells by protein denaturation and membrane damage and could cause cell death in situ [46].

In the present study, a mathematical model of tumor growth has been built by combining ATP production rate and the mechanism of the Warburg effect. It is validated by the tumor growth measurements in situ and further applied to study the hyperthermia effect on tumor metabolism over a period of time.

2. Theory

2.1. Experimental Study of Tumor Growth

Animal study of tumor metabolisms under long-term hyperthermia has been studied [47]. In the study, tumor cells were injected into the back of Balb/c mice around 6~8 weeks. Mice were grouped randomly into the control and treatment groups with 6 in each group. In the treatment group, the tumor region of the mice was heated via a 15 mm-diameter circular heating pad for 28 days. On each day, it was heated for 12 hours with an interval of 12 hours. The supplying power of the heating pad was controlled to maintain the tumor surface tissue temperature at 39°C. Tumor sizes of both the control and treatment groups were measured and recorded at 3, 7, 10, 14, 21, 28 days, respectively. The experimental results are shown in Figures 2–5 with the permission of the authors.

2.2. Model Development

At the avascular stage of a solid tumor, tumor cells proliferate and tumor grows like a spheroid without restriction. In our model, the tumor is assumed to start from one single tumor cell (as illustrated in Figure 1), and to grow into a homogenous spheroid over a period of time. This is an assumption used by many models [7, 10, 14, 16, 23, 24, 29] and acceptable when the diameter of tumor is less than 2~4 mm prior to microvascular development. In addition, only living tumor cells are supposed to take the space of the spheroid. The proliferation of tumor cells is related to local concentration of nutrients. In this model, nutrients are simplified to the main source of energy (glucose) only. Glucose diffuses passively into tumor tissue from the outer rim of the tumor, where its concentration remains at a constant over the initial tumor growth period.

Figure 1 

Schematic of tumor growth. The tumor starts with one single cell (represented by the yellow circle), and the cells differentiate continuously. When the infinitesimal unit of the space volume is filled up with the tumor cells, they grow into the next volume unit.

Figure 2 

Tumor growth curve of the experimental measurements versus the numerical simulation of the control group.

Figure 3 

Tumor growth curve of the experimental measurements versus the numerical simulation of the treatment group (39°C) (𝑛=6).

Figure 4 

Average glucose concentration of the experimental measurements versus the numerical simulation of the control group.

Figure 5 

Average glucose concentration of the experimental measurement versus the numerical simulation of the treatment group.

Usually, the glucose is metabolized through two different pathways: aerobic and anaerobic, depending on the cell status and physiological conditions. According to the recent understanding of the Warburg effect [32, 38], in tumor metabolism, about 5% of glucose undergoes aerobic pathway, and 85% takes the anaerobic pathway to produce ATP. The rest of 10% glucose is utilized for biomass synthesis necessary for cell divisions.

Through anaerobic pathway, glucose first degrades into pyruvate, and pyruvate converts into lactate by lactate dehydrogenase in the cytoplasm. While in the aerobic metabolism, pyruvate will further react with oxygen, and produce water and carbon dioxide inside mitochondria. It is clear that not only the metabolic site and end products are different, but also the amount of energy produced differs. Through aerobic metabolism, one mole of glucose consumes 6 moles of oxygen and produces 38 moles of ATP, while in anaerobic metabolism only 2 moles of ATP per mole glucose could be produced without oxygen. These reactions are simplified and presented by the following formulae

Aerobic respiration: C₆H₁₂O₆+6O₂+38ADP⟶6CO₂+6H₂O+38ATP(1)

Anaerobic respiration: C₆H₁₂O₆⟶2C₃H₆O₃+2ATP(2)

The proliferation rate of tumor cells is supposed to rely on the ATP production rate, which is determined by the metabolic reactions: 𝑄_ATP=38𝑄_𝑔,AR+𝑄_lac,(3) where 𝑄_𝑔,AR is the glucose assumption rate due to aerobic respiration, 𝑄_lac is the lactate production rate.

Following the hypothesis of Heiden [32], the glucose consumption rate due to aerobic metabolism and the lactate production rate could be derived from the uptake rate of glucose by tumor cells: 𝑄𝑄_𝑔,ARlac=5%𝑄_𝑔,=85%×2×𝑄_𝑔,(4) where 𝑄_𝑔 is the cellular glucose uptake rate.

Then the ATP production rate 𝑄_ATP could be calculated: 𝑄_ATP=5%×38×𝑄_𝑔+85%×2×𝑄_𝑔.(5)

The cellular uptake rate of glucose depends on both the extracellular level of glucose concentration and the ability of cell uptake. The reaction is assumed to be enzyme reaction and followed Michaelis-Menten kinetics’ equation [10, 25]: 𝑄_𝑔=𝑄_𝑔max𝑐_𝑔𝑐_𝑔+𝐾_𝑔,(6) where 𝑄_gmax is the maximum glucose uptake rate per tumor cell, 𝑐_𝑔 is the glucose concentration in tumor tissue, and 𝐾_𝑔 is saturate concentration.

The concentrations of nutrients (both glucose and lactate) are assumed to be functions of tumor spheroid radius that are changing with the rate of cell proliferation [7]: 𝑑𝑐_𝑔𝑑𝑡=𝐷^𝑔∇²𝑐_𝑔−𝑄_𝑔𝑙,(7) where 𝐷_𝑔 is the diffusion coefficient of glucose and reported to be 1.05 * 10⁻⁶ cm²/s for EMT6/R0 spheroid, and 𝑙 is the number of living cell per unit volume.

The concentration of glucose at the boundary is decided by the average value of glucose concentration in blood—5 mM, and the tumor is assumed symmetric, thus the boundary conditions for (5) are as follows

Boundary condition: 𝜕𝑐_𝑔𝜕𝑡||||_𝑟=0=0,𝑐_𝑔||_𝑟=𝑅=5mM,(8) where 𝑅 is the radius of solid tumor.

The distribution of lactate in tumor is determined by 𝑑𝑐_lac𝑑𝑡=𝐷^lac∇²𝑐_lac−𝑄_lac𝑙,(9) where 𝑐_lac is the concentration of lactate in tumor tissue, and 𝐷_lac is the effective diffusion coefficient of glucose and lactate. The boundary conditions are: 𝜕𝑐_lac𝜕𝑡||||_𝑟=0=0,𝑐_lac||_𝑟=𝑅=0mM.(10)

The diffusion coefficient of lactate could be determined from the glucose diffusivity [25]: 𝐷_lac=𝐷_𝑔𝑀𝑀𝑊𝑊^𝑔_lac^3/4,(11) where 𝑀𝑊_𝑔 and 𝑀𝑊_lac are molecular weights of glucose and lactate, respectively.

Assuming that all viable space in tumor tissue is filled with living and tightly packed tumor cells (see Figure 1), the number of living cells per unit tumor tissue volume is 𝑙=1𝑉_𝐿,(12) where 𝑉_𝐿 is the volume of a living tumor cell. It is obtained through in vitro measurement. Diameters of more than 150 suspended 4T1 tumor cells are measured and averaged to get the volume of a single cell. Its average radius is 13.41 μm with a standard deviation of 2.21 μm. With this radius, the tumor cell volume was calculated and the result is listed in Table 1.

Table 1 

List of variables used in modeling.

By solving the above listed equations together, the ATP production rate and concentrations of glucose and lactate inside the tumor tissue could be obtained. The relationship between ATP production rate and cell growth rate is proposed and well parameterized by Forbes et al. [29]: 𝑘_𝑙=𝐴𝑄^ATP𝑄_ATP+𝐾_ATP,(13) where 𝑘_𝑙 is cell growth rate of tumor cells, A is the maximum cell growth rate, and 𝐾_ATP is saturation ATP production rate.

In this model, tumor growth is considered as tumor cells moving outward while all space occupied. In other words, tumor growth is mainly due to living tumor cells’ proliferation, then the moving velocity of the tumor spheroid rim is given as 𝑣𝑡,𝑅_(𝑡)=∫^𝑅(𝑡)0𝑘_𝑙𝑙𝑉_𝐿𝑟²𝑑𝑟𝑅^2(𝑡),(14) where 𝑟 is radial distant from center of tumor spheroid, and 𝑅_(𝑡) is the radius of tumor spheroid at time 𝑡. The radius of tumor spheroid at time 𝑡+𝑑𝑡 is 𝑅_{(𝑡+𝑑𝑡)}=^{𝑡𝑡+𝑑𝑡}𝑣𝑡,𝑅_(𝑡)𝑑𝑡+𝑅_(𝑡).(15)

The average concentration of the metabolites is calculated from the integration of the substances’ distribution throughout tumor tissue divided by the tumor tissue volume as 𝑐_(𝑡)=∫^𝑅(𝑡)0𝑐(𝑟)⋅4𝜋𝑟²𝑑𝑟(4/3)𝜋𝑅^3(𝑡).(16)

By substituting the dimensionless quantities and parameters as defined below into the above equations, the dimensionless equations could be obtained.

The dimensionless variables are defined as 𝑐^∗𝑔=𝑐^𝑔𝑐_𝑔,0,𝑐^∗_lac=𝑐^lac𝑐_𝑔,0,𝑟^∗=𝑟𝑅^∗(𝑡),𝑡^∗=𝑡1/𝐴,(17) where 𝑐_𝑔,0=5mM.

Thus, the dimensionless equations are 𝑄^∗_ATP=85%∗2∗𝑄^∗𝑔+5%∗38∗𝑄^∗𝑔𝑄^∗𝑔=𝑄^𝑔𝑄_𝑔max=𝑐^∗𝑔𝑐^∗𝑔+𝐾_𝑔/𝑐_𝑔,0𝑄^∗_lac=𝑄^lac𝑄_𝑔max=𝑐^∗lac𝑐^∗_lac+𝐾_lac/𝑐_𝑔,0𝑑𝑐^∗𝑔𝑑𝑡^∗=1𝐴∇²𝑐^∗𝑔−𝑅(𝑡)²𝑄_𝑔max𝐷𝑉_𝐿𝑐_𝑔,0𝑄^∗𝑔𝑐_∗𝑔||_𝑟∗=1=1𝑑𝑐^∗𝑔||_𝑟∗=0𝑑𝑡^∗𝑑𝑐_∗lac=0𝑑𝑡^∗=1𝐴∇²𝑐^∗_lac−𝑅(𝑡)²𝑄_𝑔max𝐷𝑉_𝐿𝑐_𝑔,0𝑄^∗_lac𝑐_∗lac||_𝑟∗=1=0𝑑𝑐^∗lac||_𝑟∗=0𝑑𝑡^∗𝑘=0_∗𝑙=𝑘^𝑙𝐴=𝑄^∗_ATP+𝑄𝐾^∗ATP_ATP/𝑄_𝑔max𝑅_(𝑡)=^𝑡0𝑣^∗𝑅_(𝑡)𝐴𝑑𝑡^∗1/𝐴+𝑅⁰.(18)

Tumor growth in vivo is actually a complex process, which involves many influencing factors such as gene mutation, immune system, tumor cell mechanism effect, tumor angiogenesis, metabolic waste, and tumor microenvironment. It is difficult to include all these factors into one single model. The present model is built based on energy production only and all other influencing factors are lumped into the maximum cell growth rate parameter 𝐴. The behavior of tumor growth is determined by the maximum cell growth rate 𝐴 and it is determined by fitting to the animal experiment data. The other parameters used are all listed in Tables 1 and 2. The constants (𝐷_𝑔, 𝑄_𝑔max) listed in Table 2 are taken from EMT6/Ro tumor, whose biophysical constants have been well studied. EMT6 is a mouse breast cancer cell line that grows in the Balb/C strain. 4T1 cells are assumed to be similar to EMT6 biophysically in the present modeling.

Table 2 

List of constants used in modeling (cell type—EMT6/Ro).

To solve the equations listed above, discrete algorithm is built. The one-dimensional tumor space is divided into 1000 intervals evenly. The time increment Δt is set to be small enough to guarantee the solver stability. The diffusion-reaction equations are then differentiated using the finite difference method. The tumor radius at certain time 𝑡 is used as an input into the differentiated equation, and the forward elimination and backward substitution method is used to solve these differentiated equations to obtain the glucose concentration distribution, the corresponding ATP production rate, tumor growth rate, the velocity of the moving boundary of tumor, and the new tumor radius at a given time 𝑡. The updated tumor radius is then used as an input into the equations iteratively, until the difference between two iterations is small enough. For the next time step (𝑡+Δ𝑡), the radius of tumor is derived from the original radius at time 𝑡 and radius increment during Δ𝑡. This process is repeated until the final convergence reached in simulation.

The maximum cell growth rate parameter 𝐴 is fitted from the animal experimental data. The reciprocal of 𝐴 ranged from 0~50 h, with an interval of 0.2 h. Each 𝐴 is then substituted to the mathematical model to obtain a simulated curve, which is compared with the animal experimental data using the coefficient of determination (R2) as a criteria. The optimal 𝐴 is chosen in correspondance to the maximum achieved R2 value.

3. Results and Discussion

By fitting the numerical model with the experimentally measured tumor growth data, the parameter 𝐴 for tumor growth with and without treatment was obtained. The fitted curve and the experimental results were shown in Figures 2 and 3, respectively.

For the control group, the tumor growth curve during the first 28 days was well captured by the energy-based model developed in this study (Figure 2). The fitted maximum cell growth rate parameter 𝐴 was (18 h)⁻¹. However, after 28 days, the tumor growth slowed down and became nonlinear, which could no longer be described by the current mathematical model. It implied that alternations in tumor metabolism and some impedimental mechanisms might appear which led to cell death when tumor radius reached about 4~5 mm in diameter. These factors were not considered in the present model.

Tumor growth in mice under long-term mild hyperthermia treatment was also fitted and shown in Figure 3. It was obvious that the growth rate of tumor was much slower and the fitted maximum cell growth rate 𝐴 was (24 h)⁻¹. The model predicted the general trend of tumor growth well in the first 15 days. The oscillation occurred in the growth curve that was not observed in the control group could be attributed to the thermal interference, although the causing mechanism is yet to be further explored. Moreover, it was clear that the model overestimated the tumor growth under the hyperthermia condition after 15 days. This indicates that more complex effects on the system might be triggered by hyperthermia after a period of time. It could be related to tumor angiogenesis, the activities of enzymes in different pathways associated with tumor cell metabolism and proliferation and so forth, which cannot be simply modeled based on energy consumption.

With the fitted parameter A, the distribution of glucose in tumor has also been obtained. The average concentration of glucose in the tumor was calculated and compared to the experimental results. Seen from Figures 4 and 5, the average concentrations of glucose tended to decrease which agreed well with the experimental results. There existed significant differences between experimental and numerical results, especially for the group under the long-term mild hyperthermia treatment. This was likely because in the present model only glucose consumption for cell proliferation was considered. There could be excessive supply of glucose from the conceivable blood vessels inside tumor and the decreased cellular uptake of glucose due to heat. A more sophisticated model taking these factors into consideration should be developed to accurately predict the glucose distribution in tumor tissue.

As there are no publications on the nutrition diffusion and cell uptake rate of 4T1 tumor cells, the data for EMT6/Ro tumor has been used. Although EMT6/Ro tumor cells are expected to possess similar properties, the dependence of these parameters on tumor species might also introduce some errors in the simulations. Therefore, parametric studies were performed and the influences on the glucose uptake rate and the diffusion coefficient were studied as shown in Figures 6 and 7. It was clear that the variation of tumor growth caused by 10% changes of either two factors was less than 5%.

Figure 6 

Parametric study of the tumor growth rate under different glucose uptake rates in tumor cells. Solid line stands for glucose uptake rate from reference [29]. Dot and dash lines represent 10% increase and decrease, respectively.

Figure 7 

Parametric study of the tumor growth rate under different glucose diffusion coefficients in tumor tissue. Solid line stands for glucose diffusion value from reference [29]. Dot and dash lines represent 10% increase and decrease, respectively.

In the present model, besides the Warburg effect, the influence of all other factors such as blood perfusion, immunity, and other growth or necrosis factors, has been lumped into a simple parameter 𝐴. Using tumor growth data from the experimental studies, the maximum tumor cell growth rate in vivo under the normal condition was fitted to be (18 h)⁻¹, and the linear tumor growth at the early stage was successfully modeled without consideration of cell death. The growth was found to be excessively inhibited by the long-term mild hyperthermia treatment. In fact, during the treatment, both the temperature gradient and inhomogeneity existed, and the higher temperature could impair tumor cells by denaturation or destruction of cellular membrane, cellular skeleton, and nucleus [48]. The long-term mild hyperthermia treatment might also affect the activities of enzymes in different pathways associated with tumor cell metabolism and proliferation, stimulate the immunology factors such as hsp70, which could arouse body system defense and eliminate tumor cells specifically [49]. Besides, Song’s study [50] revealed that under mild hyperthermia (41~42°C, 30 min) blood perfusion in tumor tissue could increase 1.5~2 folds as compared to that prior to the treatment. Tumor vasculature is very sensitive to heat, as it is loosely organized and usually lacks of smooth muscles [51–54]. Tumor cells would suffer from starvation due to the damage of the angiogenesis. All these influences on tumor growth are not linear and could not be accurately represented by a single parameter. This warrants further investigation into more detailed modeling of the long-term hyperthermia effect on tumor growth in the near future.

4. Conclusion

An energy-based model linking tumor growth with cell proliferation rate has been developed in this study to investigate the hyperthermia treatment effect. In the model, the new understanding of the Warburg Effect was for the first time taken into account for tumor cellular metabolism regardless of the concentration of oxygen. The maximum cell growth rate was used as an integrated variable responding to changes under different environments. Trends of initial tumor growth and changes of the average glucose concentration in tumor were successfully modeled. The comparison of the maximum tumor cell growth rate has revealed a slowdown of tumor growth under the long-term mild hyperthermia condition. To accurately predict the tissue glucose level and the corresponding metabolites, especially under the long-term hyperthermia, the model needs to be further developed.

Acknowledgments

This work has been supported by the National Natural Science Foundation of China (NSFC51076095, NSFC50725622) and Shanghai Municipal Science and Technology Commission (10QA1403900).