Absence of convection in solid tumors caused by raised interstitial fluid pressure severely limits success of chemotherapy—a numerical study in cancers

Abbreviations: TFIP: Tumor interstitial fluid pressure; DTC: Disseminated tumor cells

1.   Introduction

There is a remarkable difference between the treatment responses of generalized lymphomas and leukemias on the one hand and solid neoplasms, i.e. cancer, on the other. While the former can be successfully treated with chemotherapy and/or antibody treatment, the latter cannot. This treatment failure in carcinomas has been attributed to the high tumor interstitial fluid pressure (TIFP) observed in most cancers, which is absent in lymphomas and leukemias. TIFP is caused by a fluid imbalance within the tumor mass. In particular, functioning lymphatic vessels are absent within the tumor so that interstitial fluid, normally transported away by the lymphatic vessels, remains in the local area. Two forms of fluid transport within tissues can be distinguished, namely diffusion and convection. Diffusion is the random movement of molecules and particles from a higher concentration region to one of lower concentration which is driven by a concentration gradient. In contrast, convection is the directed flow of particles within the interstitial fluid from the blood-tissue exchange microvessels to the lymphatic drainage system, which is pressure-related. The latter is absent in cancers so transport processes are more or less limited to diffusion ^[1,2]. For both macromolecules such as antibodies ^[3,4] and small molecules such as doxorubicin ^[5] or cisplatin ^[6], xenograft experiments showed that anti-cancer drugs only penetrate about 100 µm around blood vessels, and therefore target only a minority of cancer cells. As these drugs are transported in the same way as oxygen and nutrients, their access to the tumor cells via the existing blood vessels is accordingly limited. Due to the lack of oxygen and nutrients, a fraction of the oxygenated proliferating cells, those that are more than 100 µm from a blood vessel, switch into a hypoxic and later into a necrotic state.

Due to this low penetration depth, cancers consist mainly of hypoxic and necrotic areas; this leads to treatment resistance because hypoxia enhances chemo- and/or radioresistance of cancer cells ^[7]. After the removal of oxygenated cells by chemotherapy or radiotherapy, formerly hypoxic cells are resupplied with oxygen and nutrients, leading to an oxygenated state. This process is called reoxygenation, and it also rescues hypoxic cells from becoming necrotic. The process of reoxygenation is utilized during radiotherapy, where small radiation fractions are applied at short intervals, leading to radiosensitization of previously hypoxic tumor areas ^[8]. Saggar and Tannock investigated chemotherapy-induced reoxygenation using hypoxia markers to quantify the number of cells that are reoxygenated after therapy in MCF-7 and PC-3 xenograft mouse tumor models ^[9]. They showed that after 24 hours, approximately 90% of the hypoxic cells were already reoxygenated in xenograft primary tumors of both cell lines.

Mathematical models are widely used to describe the increase in tumor volume over time, the spread of metastases, and the effect of different treatment strategies ^{[10,11,12,13,14]}. In order to create biologically meaningful translational predictions (e.g., optimal treatment doses), these models have to be calibrated with experimental data, validated with untrained data, and evaluated on predictive performance for known treatments ^[15]. Mathematical models can also be useful for developing a better understanding of biological mechanisms, which can stimulate new experimental approaches. A variety of mathematical models exist that can describe tumor growth including Gompertz ^[16,17], exponential ^[18], von Bertalanffy ^[19], or the generalized logistic model ^[20,21]. These models only consider the change in tumor volume over time, not the behavior of individual fractions within the tumor. Ribba et al. published a mathematical model that successfully describes tumor growth and tumor fractions over time based on longitudinal data and histological biomarkers ^[10].

In order to better understand the limited success of chemotherapy in solid neoplasms, we modeled the effect of limited distribution of chemotherapeutic drugs into the primary tumor by extending this mathematical model ^[10], adding a pharmacokinetic-pharmacodynamic model ^[22], and including the reoxygenation of hypoxic cells. In a second step, the effect of chemotherapy on disseminated tumor cells (DTCs) and multicellular metastases was studied. DTCs and small metastases are of particular relevance as they consist only of tissue that is well supplied with nutrients and oxygen and have not built up increased TIFP. Only once metastases have grown to a certain size does TIFP build up, leading to the formation of hypoxic and necrotic tumor fractions. Changes were identified in the overall effectiveness of chemotherapy, which decreased with the progression of the disease. TIFP also has adverse effects on the treatment of metastases and is, therefore, a critical factor that must be considered in developing treatment strategies for malignant tumors.

2.   Materials and methods

2.1.   Longitudinal tumor data and histological biomarkers

To investigate the limited convection of chemotherapeutic drugs into the primary tumor, we used experimental data, including tumor size and histological biomarkers, which describe two fractions of the tumor over time: hypoxic and necrotic tissue ^[10]. In this study 6 × 10⁶ HT29 cells were subcutaneously injected into 15 female athymic nude mice. Tumor size was measured at repeated intervals using a caliper. The percentage of hypoxic and necrotic tissue was determined using immunohistochemistry. The mean values of tumor size (mm) and hypoxic and necrotic fractions (%) for each measured time point were obtained using WebPlotDigitizer ^[23]. The oxygenated fraction of the tumor was deduced from the hypoxic and necrotic fractions.

2.2.   Metastasis formation and dormancy behavior

Since no metastases were examined in the original study, we used data from an HT29 xenograft mouse model ^[24], which provided information about the number and size of the metastases at the end of the experiment. In this study 1 × 10⁶ HT29 cells were subcutaneously injected into pore-forming protein and recombination activating gene 2 double knockout (pfp/rag2) mice and in recombination activating gene 2 only knockout (rag2) mice. In pfp/rag2 mice typical metastases (~ 789) consisted of approximately 10–100 tumor cells. In rag2 mice, only disseminated tumor cells were generally found (~ 200). Computer simulation of the metastatic spread proved that metastases stayed dormant for at least 30 days after reaching a size between 10 and 100 cells in the pfp/rag2 group. In contrast, the metastatic cells in the rag2 group entered an initial dormant phase, in which they remained for 30 days before starting to proliferate.

2.3.   Reoxygenation of hypoxic tumor cells

Investigations of human breast (MCF-7) and prostate cancer (PC-3) primary xenograft tumors have shown that chemotherapy induces reoxygenation and repopulation of hypoxic cells ^[9]. In this study 5 × 10⁶ MCF-7 and 2 × 10⁶ PC-3 cells were subcutaneously injected into female and male athymic nude mice, respectively. The tumor hypoxia marker pimonidazole was given concurrently with doxorubicin. At 24, 48, 72, 96 and 120 hours following initial treatment, a second tumor hypoxia marker EF5 was injected to quantify reoxygenated tissue. The reoxygenation process rescues hypoxic cells from conversion to a necrotic state by resupplying oxygen and nutrients (Figure 1). Approximately 90% of the hypoxic cells were converted to an oxygenated state after 24 hours. After 120 hours, around 100% of formerly hypoxic cells were reoxygenated. The number of reoxygenated cells should depend on the number of eliminated oxygenated cells. For example, if only 20% of the tissue is destroyed during therapy, it is implausible that all hypoxic cells will be converted back into an oxygenated state if the number of cells in the hypoxic tissue is higher than the initial number of cells in the oxygenated tissue. We assume that rescued hypoxic cells are able to proliferate after being converted to an oxygenated state.

Figure 1.  Schematic representation of the effect of reoxygenation of the hypoxic fraction after treatment. The drug destroys a specific number of tumor cells in the oxygenated fraction of the tumor. Depending on the number of cells in the oxygenated fraction that are destroyed, some (or all) of the hypoxic cells are resupplied with nutrients. The cell numbers in the hypoxic fraction decrease. P: Oxygenated tissue, Q: Hypoxic tissue, N: Necrotic tissue.

DownLoad: Full-Size Img PowerPoint

2.4.   Mathematical model

An existing tumor growth model was coupled with models of the pharmacokinetics and pharmacodynamics of a cell-cycle non-specific drug such as cisplatin (an alkylating agent), reoxygenation of hypoxic tissue after treatment, and metastasis formation and dormancy behavior (Figure 2).

Figure 2.  Schematic description of the full mathematical model. The mathematical model consists of two main parts: The pharmacokinetic/pharmacodynamic (PK/PD) part describing drug behavior and the part modeling tumor growth (Tumor/Metastasis Growth ^[10]). Each compartment of the primary tumor/metastasis contains cancer cell species (oxygenated, hypoxic, necrotic) that characterize the amount or concentration of an entity, e.g., the number of oxygenated cells (P). Solid black arrows describe the interaction between compartments. For example, the drug in the central compartment moves into the peripheral compartment, where the rate constant k₁₂ regulates the amount of the drug. Gray dotted arrows indicate where the quantity of the cancer cell species is specified by the amounts or concentrations of other cancer cell species and not by rates. Black dashed arrows describe dependencies, e.g., the maximal tumor size is limited by the carrying capacity K, which depends on the rate of increase in the carrying capacity. A complete description of the model is given in the main text.

DownLoad: Full-Size Img PowerPoint

2.4.1.   Mathematical model of tumor growth

The mathematical model of tumor growth published by Ribba et al. ^[10] describes the growth of cancer and its tumor fractions over time. The tumor consists of three main compartments: Oxygenated (P), hypoxic (Q), and necrotic tissue (N). P^* represents the total tumor size. The oxygenated tissue follows a generalized logistic equation, where the maximum size is limited by a carrying capacity (K) that describes the maximum tumor size that can be reached with the amount of physical space and nutrients available. Therefore, the carrying capacity depends on the rate of increase in carrying capacity for a specific cell line (b) induced by angiogenic processes. Since we assume that only oxygenated cells are able to proliferate, λ_P is the growth rate constant for the oxygenated tissue. The parameter s regulates the slowdown of the logistic growth of oxygenated tissues and $\alpha$ reflects the effect of hypoxic stress. Cell transfer from oxygenated tissue (P) to hypoxic tissue (Q) and hypoxic tissue (Q) to necrotic tissue (N) is modeled via the rate constants k_PQ and k_QN, respectively. The system equations are given by the following:

2.4.2.   Pharmacokinetic-pharmacodynamic model

We modeled the pharmacokinetics (PK) of a cell-cycle non-specific drug using a two-compartment model. The system is divided into two kinetically different regions: Central (C_D) and peripheral (P_D) compartment ^[22]. The central compartment consists of the plasma and tissues where the drug is homogeneously distributed right after the administration. In contrast, drug distribution in the tissue of the peripheral compartment is slower. The intravenous bolus administration is modeled via the addition of the dose to the central drug compartment. The parameters k₁₂ and k₂₁ are rate constants characterizing drug movement between the central and peripheral compartment. k_e represents the rate constant for drug elimination (clearance) out of the central compartment, which leads to a drug concentration c(t) that is following an exponential decay function (as illustrated in the "PK/PD" sub-model in Figure 2). The corresponding model equations are given by:

To include the effect (PD) of limited diffusion of chemotherapeutic agents into the primary tumor, the drug only affects cells in the oxygenated layer (P), because drugs cannot diffuse into the hypoxic layer due to the limited penetration depth of 100 µm around blood vessels ^[5,6]. The parameter k₂ represents the drug potency. Treatment was applied on days 10, 20 and 30.

2.4.3.   Reoxygenation of hypoxic cells

The process of reoxygenation is modeled by the fraction of the hypoxic cells (k_QPQ) that is transferred to the oxygenated compartment after treatment, where k_QP is the constant transfer rate. This approach is similar to models of the reoxygenation process after chemotherapy or radiotherapy treatment ^[25,26,27]. Since approximately 90% of hypoxic cells were already reoxygenated after 24 hours in a previous study ^[9], we assume that the reoxygenation process begins immediately after chemotherapy. Cell transfer from the hypoxic to the oxygenated compartment remains active until the number of cells in the oxygenated layer reaches the number of cells before treatment. This decision is based on experiments showing that chemotherapy does not lead to better blood supply through angiogenesis, but to a reduction in blood vessel density, which is still detectable several days after the therapy ended ^[28]. Therefore, the number of cells in the oxygenated layer during the reoxygenation process cannot be higher than the number before treatment.

2.4.4.   Spread and growth of metastases

The primary tumor spreads single malignant cells at a colonization rate $\mu \left({P}^{*}\right)$, which depends on the per day probability m that any individual cell will spread and establish a distant metastasis, the number of cells in the primary tumor P*, and the fractal dimension of blood vessels $\delta$ ^[29]. A fractal dimension of 2 implies that the tumor is superficially supplied with blood vessels, while a value of 3 means a completely supplied tumor.

A cumulative trapezoidal numerical integration of the colonization rate was performed to derive the time of development of each metastasis from the integrated values. Each time the integrated values reach the next integer value, a new metastasis is created.

The growth of metastases starts from a single cell and follows the growth pattern of the primary tumor. Metastases also consist of oxygenated, hypoxic, and necrotic tissue. Since nutrients are only transported up to 100 µm away from the blood vessels, we assume that the formation of a hypoxic and necrotic layer only begins when the metastasis reaches a threshold diameter of 200 µm with an assumed superficial blood supply. Therefore, the initial values of the rate constants k_PQ and k_QN were set to 0 until the size of the metastasis reached this threshold.

2.4.5.   Modeling late dormancy

Malignant cells emitted from the primary tumor can stay dormant for a certain length of time or enter a late dormancy state after they reach a specific size ^[24]. In this article we focus on the case of metastases entering a late dormancy phase of 30 days. With the chosen treatment strategy, it makes no difference whether 1 or 100 cells are affected since chemotherapy is parameterized in such a way that all treatable cells will be destroyed. Furthermore, all metastases are in a dormant phase at the selected treatment times (days 10, 20 and 30). The growth parameter λ_P and transfer rate constants (k_PQ, k_QN) are set to a value of 0 during the time span of dormancy (30 days). Based on previous experiments on metastasis formation, the number of cells at which metastases enter a late dormant state is calculated randomly for each metastasis based on a uniform distribution between 10 and 100 cells ^[24].

2.4.6.   Combined model

The full mathematical model including tumor growth, chemotherapy treatment, reoxygenation of hypoxic cells, and dormancy of metastases is defined by:

2.5.   Model calibration

Experimental data was used to determine model parameters and rate constants in order to calibrate the model.

2.5.1.   Primary tumor growth

The growth of tumor size in experiments is usually measured as diameter in mm. For computer models the number of cells is required. The diameter of the tumor (mm) was converted to volume (mm³) and multiplied by 10⁶ using the well-established conversion rule (1 mm³ ~ 1 × 10⁶ cells) to determine the number of cells ^[30]. The number of hypoxic and necrotic cells was calculated from the percentage of tumor fractions at each measuring time point. Least-squares minimization was performed using fmincon from the MATLAB Optimization Toolbox (MATLAB R2019b, The MathWorks Inc., Natick, USA) to obtain all model parameters based on experimental tumor data. The number of initial tumor cells in the oxygenated compartment (P) was set to 6 × 10⁵ cells during the parameter estimation process because our previous study showed that only 10% of the subcutaneously injected tumor cells survive the first 24 hours ^[31]. Furthermore, fixation of the initial tumor cell number to a well-known value reduces the number of degrees of freedom and improves the parameter estimation results.

2.5.2.   Pharmacokinetic-pharmacodynamic model

Due to the absence of pharmacokinetic and pharmacodynamic data for the experimental data analyzed here, we used determined parameters from Simeoni et al. ^[32], which describe the drug concentration of the chemotherapy agent paclitaxel over time. We adjusted the drug potency parameter k₂ until all cells in the oxygenated fraction of the tumor and metastases were eliminated (in the absence of reoxygenation) to mimic the most effective treatment. Since the cells in the oxygenated layer will be entirely destroyed, the limiting effect of the interstitial fluid pressure on the efficacy of the chemotherapeutic drug can be studied most clearly.

2.5.3.   Reoxygenation

Previous experimental data from human MCF-7 and PC-3 primary xenograft tumors in mice were used to determine the reoxygenation rate k_QP of hypoxic cells ^[9]. These data provide the percentage of reoxygenated cells 24, 48, 72, 96 and 120 hours after the application of chemotherapy. They show that 90 and 100% of hypoxic cells were reoxygenated by 24 and 120 hours, respectively. These percentages represent the mean values from four mice. Based on these percentages, growth data for the reoxygenation process after the application of chemotherapy was generated. For example, 1 × 10⁵ cells were killed by chemotherapy in the oxygenated layer. After 24 hours, 9 × 10⁴ cells (90%) are already reoxygenated. After 120 hours, the oxygenated layer has regained its pre-treatment cell count. The generated data were used to estimate the reoxygenation rate constant k_QP using lsqnonlin (trust-region reflective algorithm). The number of cells that are reoxygenated can change depending on the drug concentration in the system. Since we assume that reoxygenation starts immediately after the therapy, reoxygenated cells can also be killed. The reoxygenation process remains active until the number of cells has reached the cell count in the oxygenated layer before treatment.

2.5.4.   Spread of metastases

We assumed a superficial blood supply (fractal dimension α = 2) for the primary tumor due to the fast growth behavior that is observed in typical xenograft experiments. Based on experimental data ^[24], we adjusted the per day probability m for each cell to spread and establish a distant metastasis until the metastasis count reached a value of 789.

2.6.   Simulation procedure

The mathematical model was implemented in MATLAB SimBiology^®-Command line (MATLAB R2019b Update 4, 9.7.0.1296695, The MathWorks Inc., Natick, USA). We used the SUNDIALS solver (Suite of Nonlinear and Differential/Algebraic Equation Solvers) to solve the model's ordinary differential equations. The simulation time was set to 50 days. The first step was the parametrization of the primary tumor, the spread of metastases, and reoxygenation based on the model calibration results. After this simulation, the metastasis creation time was calculated based on the colonization rate $\mu \left({P}^{*}\right)$. The last step was the simulation of each metastasis. This procedure allowed us to examine the course of the individual fractions of each metastasis for further evaluations.

3.   Results

3.1.   Calculations and observations for primary tumor fit well

The tumor growth model was calibrated with experimental data, including longitudinal tumor data and histological biomarkers (hypoxic and necrotic fractions). The initial number of surviving engrafted tumor cells after subcutaneous injection (P) was set to 6 × 10⁵ cells. The predictions were in good agreement with the observed values from experimental data (tumor cells R² = 0.95, hypoxic cells R² = 0.82, necrotic cells R² = 0.97). The parameterized model was able to describe the size of the tumor and the dynamics of its hypoxic and necrotic fractions based on the number of cells (Figure 3). The estimated parameters are displayed in Table 1.

Figure 3.  Model predictions vs observed values and simulation results of primary tumor growth including 95% confidence intervals. Panels A–C show the individual model predictions of tumor cells (A), hypoxic cells (B), and necrotic cells (C) compared to the observed values from experimental data. Panels D–E show model simulations with a 95% confidence interval. Tumor size (A, D), hypoxic (B, E) and necrotic (C, F) fractions.

DownLoad: Full-Size Img PowerPoint

Table 1.  Parameter description and estimated values.

Parameter	Description	Value
P0	Initial tumor size	6.00e5 cells*
K0	Initial carrying capacity	1.1193e9 cells
λP	Growth rate for the oxygenated tissue	1.8420 day-1
kPQ	Transfer rate from oxygenated to hypoxic tissue compartment	0.3595 day-1
kQN	Transfer rate from hypoxic to necrotic tissue compartment	0.1141 day-1
kQP	Transfer rate from hypoxic to oxygenated tissue compartment	2.2216 day-1
b	Rate of increase of carrying capacity	0.5811 day-1
m	Per day probability to spread and establish a distant metastasis	2.732e-8 (cell day)-1
δ	Fractal dimension of blood vessels	2
α	Effect of the hypoxic stress on the growth	0.1* [10]
dose	Drug dose	20 arbitrary unit*
k12	Transfer rate from central to peripheral compartment (PK/PD)	0.0115 day-1 [32]
k21	Transfer rate from peripheral to central compartment (PK/PD)	0.0616 day-1 [32]
ke	Rate constant for drug elimination (clearance)	3 day-1 *
k2	Drug potency (effectivity)	2.1 *
*Known parameters were set to a fixed value during the parameter estimation procedure, which increased the accuracy of the estimates and reduced the number of degrees of freedom.

---

 | Show Table

DownLoad: CSV

3.2.   Percentage of oxygenated cells in a tumor decreases rapidly with time

Some time after subcutaneous injection of HT29 tumor cells into mice, the engrafted surviving cells form a local tumor nodule and develop hypoxic and necrotic fractions depending on the blood vessel geometry. Experimental tumor data, including histological biomarkers of hypoxic and necrotic tissue, were used to estimate model parameters and rate constants. Figure 4 shows the simulation output of untreated tumor growth including the percentage of tumor fractions over time. According to the model output, the percentage of oxygenated cells decreased significantly over time compared to the tumor size. This finding implies that the tumor consists mainly of hypoxic and necrotic tissue, which is in good agreement with the experimental data ^[10].

Figure 4.  Dynamics of untreated tumor growth and tumor fractions over time. Panel A shows the change in size of the tumor over time. Panel B shows percentages of tumor fractions over time. The number of cells in the oxygenated tissue increased slowly (A), while the percentage decreased in relation to the size of the tumor (B). Hypoxic and necrotic tissue compartments represent the most substantial parts of the tumor. Panel C shows the tumor structure in layer form over time.

DownLoad: Full-Size Img PowerPoint

3.3.   Effectiveness of chemotherapy on tumors decreases rapidly over time

Chemotherapeutic drugs only penetrate about 100 µm around blood vessels, and therefore only target cancer cells within this range ^[5,6]. The percentage of oxygenated cells in a tumor after chemotherapy could be understood as an indicator of the success of the therapy since hypoxic (and necrotic) cells are chemotherapy resistant ^[7]. One reason for this resistance is their increased distance from blood vessels, so that no chemotherapeutic drugs can reach these cells ^[5,6]. Furthermore, this distance also causes nutrient and oxygen deficiency, which slows down the cell cycle and, as a result, the proliferation rate ^[33]. Therefore, cell-cycle specific drugs are less effective in these cells because they are only active in a particular phase of the cell cycle. We investigated the effectivity of cell-cycle non-specific chemotherapy, which can treat cells in all phases of the cell cycle. Since cell-cycle non-specific drugs destroy more cells than cell-cycle specific drugs, we can focus on the limited distribution of the chemotherapeutic drug.

The drug potency parameter k2 for cell-cycle non-specific chemotherapy was manually adjusted until the number of cells in the oxygenated layer was smaller than one cell layer after treatment in the absence of reoxygenation. A value of k2 = 2.1 was determined for all three treatment time points (days 10, 20 and 30). The number of cells in the oxygenated layer is affected by the occurrence of reoxygenation of hypoxic cells, which is assumed to start immediately after cells are destroyed in the oxygenated layer. This approach is based on our assumption that more cells are destroyed if a sufficient concentration of the chemotherapeutic drug is still present, depending on the PK/PD model's parameterization. Figure 5 shows the results of model simulations of chemotherapy treatment at three different treatment times (days 10, 20 and 30, panels A–F) as well as the percentage distribution of the tumor fractions before and after treatment (panels G–I). Since the formation of the necrotic portion of the tumor is not reversible and treatment of these cells has no effect, the oxygenated and hypoxic fractions represent the potentially treatable part of the tumor. Over time, the treatable part of the tumor became smaller and smaller, and hence a more substantial portion of the tumor could be destroyed with early treatment. However, the fact that penetration of the drug is limited to 100 µm around blood vessels prevents treatment of the hypoxic fraction, even for small tumors. Panels A–C show the course of tumor growth for three different treatment times (days 10, 20 and 30). As shown in panel D (blue dotted lines), the oxygenated cells were killed in the first days after the application of chemotherapy. The decreasing concentration of the drug, combined with reoxygenation of the hypoxic cells and the proliferative behavior of the oxygenated cells, led to recovery of the oxygenated tissue. The hypoxic layer also decreased in size due to the reoxygenation process, as hypoxic cells switched into an oxygenated state (panel E). Since reoxygenation prevented a proportion of hypoxic cells from entering the necrotic state, the growth of the necrotic layer slowed down (panel F). The values for the oxygenated, hypoxic, and necrotic fractions were determined when the drug concentration (central drug) reached a value below 1% of its peak concentration (day 2 post-treatment) in order to compare the tumor fraction distribution before and after therapy. This drug concentration is considered to have a negligible therapeutic effect ^[34,35]. If treatment commenced on day 10, the proportion of potentially targetable tumor cells was higher than if it starts on day 20 or 30 (panels G–I). The percentage of necrotic tissue before and after treatment provides information about the therapeutic efficacy. The necrotic portion did not grow abruptly after treatment but did represent a higher percentage of the tumor volume since the oxygenated tissue had been killed and a fraction of the hypoxic cells had been reoxygenated. The difference between the necrotic tissue fraction before and after therapy was 49% (G), 47% (H) and 37% (I) on treatment days 10, 20 and 30, respectively. These results imply that the timing of chemotherapy treatment is critical in the presence of interstitial fluid pressure, which leads to the formation of hypoxic and necrotic tumor fractions.

Figure 5.  Simulation results of tumor growth under chemotherapy on different days. Panels A–C show the growth curves of an untreated tumor (solid lines) and the treated tumor (dashed lines). Treatment was applied on day 10 (A), day 20 (B), and day 30 (C). Panels D–F show the effects of different treatment times on the oxygenated (D), hypoxic (E) and necrotic layers (F). Panels G–I compare the percentage distribution of the individual tumor fractions before and two days after treatment for tumors treated at day 10 (G), 20 (H) and 30 (I), respectively.

DownLoad: Full-Size Img PowerPoint

The application of chemotherapy can possibly disrupt the blood vessel geometry ^[28]. As an effect of the altered geometry, cell transfer from the hypoxic to the oxygenated compartment after chemotherapy might lead to a reduced number of cells after treatment only due to the smaller amount of cells that are in contact with blood vessels. To study this effect, we performed the same simulation with a modified reoxygenation behavior (Figure S1). In this scenario, cells were transferred from the hypoxic to the oxygenated compartment until the number of cells in the oxygenated compartment reached 80% of the cell count before treatment. This percentage represents a reasonable assumption in our view since no experimental data is available. Similar results were obtained in comparison to the simulation with a reoxygenation behavior in which the number of cells is completely restored during the reoxygenation process indicating that the percentage of the oxygenated fraction is of lesser importance.

3.4.   Interstitial fluid pressure protects metastases from chemotherapy

Inspired by the simulation results for the effects of chemotherapy on the size of the primary tumor, we applied the same procedure to understand its impact on the spread of malignant cells and metastases. The colonization rate µ(P^*) was parametrized in such a way that an untreated primary tumor spread and established 789 distant metastases, which was the observed mean number from previous HT29 xenograft experiments ^[24]. We assumed a superficial blood supply (fractal dimension of 2) due to the rapid growth of the primary tumor. In some of our previous research, the size of metastases at the end of the experiment ranged from 10 to 100 cells, which implied a late dormancy after metastases reached a specific size ^[24]. Two main scenarios were considered, namely metastases with and without dormancy behavior. Chemotherapy was applied on day 30 in both scenarios with a simulation time of 50 days. In addition, effects of the resection of the primary tumor two days after chemotherapy treatment commenced were simulated for both scenarios.

Figure 6 shows different simulation scenarios to investigate the influence of chemotherapy on metastasis formation. The untreated primary tumor established 789 distant metastases. In contrast, the treated primary tumor, which received chemotherapy on day 30, established 683 distant metastases or disseminated tumor cells due to a reduced spreading behavior because the reduction in the size of the primary tumor ${P}^{*}$ affected the colonization rate $\mu \left({P}^{*}\right)$. Chemotherapy had no substantial effect on metastases that had already built up interstitial fluid pressure. Only 56 metastases that were smaller than 200 µm in size at the time of commencement of chemotherapy were destroyed (blue bar). In the scenario where metastases switched into a late dormancy state (Figure 7B, C), chemotherapy destroyed 266 metastases. The remaining 417 metastases were disseminated after application of chemotherapy and could not be treated. With an additional resection two days after chemotherapy, 761 distant metastases or disseminated tumor cells were eliminated.

Figure 6.  Effect of chemotherapy on the number and size distribution of metastases in different scenarios. Panel A shows the number of metastases after different treatment strategies and in varied growth conditions. Panel B shows the size distribution of metastases for chemotherapy when there was no dormancy of the metastases, leading to build-up of interstitial fluid pressure (ChT). Panel C shows the size distribution of metastases for chemotherapy when metastases showed dormancy behavior (ChT + dormancy). The chemotherapy was applied on day 30. The simulation time was set to 50 days. Gray bars represent reference values (Ref.) for an untreated and a tumor treated with chemotherapy to compare treatment effectivity.

DownLoad: Full-Size Img PowerPoint

Figure 7.  Synthetic parameter samples and tumor growth response. Panels A–D show the histogram and distribution of drug clearance k_e and drug potency k₂ varying around its parameter value (gray dash-dot line). Panels E–H show simulation outputs of tumor growth behavior for tumor size (E), oxygenated cells (F), hypoxic cells (G), and necrotic cells (H). The black line indicates the model output with the estimated values whereas gray dashed lines represent model outputs with synthetic samples.

DownLoad: Full-Size Img PowerPoint

In addition, we studied blood vessel disruption after chemotherapy to analyze its influence on spreading behavior. As in the utilized model of vascular tumor growth ^[10], the disruption of blood vessels was realized with the formula $\frac{dK}{dt} = b{P}^{*}-\xi K$, where $\xi$ was set to an arbitrary value of 0.1 during treatment. Since no data about the alteration of the blood vessel geometry is available for the used experiment, the changed fractal dimension of blood vessels after the chemotherapy was estimated from our previous studies ^[28]. To obtain a biologically realistic value we adjusted the fractal dimension after chemotherapy to the mean value of 1.58. Supplementary Figure S2 represent the same simulation scenarios to investigate the influence of chemotherapy on metastasis formation. In general, the disruption of blood vessels has a positive effect on the spreading behavior of malignant cells, since the primary tumor spread fewer malignant cells. Although surgical removal of the primary tumor stops the tumor's spread, the resistance of the metastases due to the increased interstitial fluid pressure is the most prominent obstacle in the treatment of cancer.

As with the primary tumor, the timing of the initiation of chemotherapy plays an important role in the treatment of metastases, because in small metastases interstitial fluid pressure has not been built up. Interstitial fluid pressure protects metastases from chemotherapy if they are not in a dormant state or have grown larger than 200 µm in diameter, leading to the development of hypoxic and necrotic fractions. Furthermore, the results show that one cycle of chemotherapy is not sufficient to kill all metastases if the primary tumor is not removed by surgery.

3.5.   Drug clearance is the most sensitive parameter for tumor growth behavior

Chemotherapy profiles from real published experiments were used ^[32]. To understand the limitations of the present approach, we generated 25 synthetic parameter sets to simulate different treatment effects on a primary tumor that received chemotherapy on day 30, as previous results (Figure 5) showed that this treatment time point was the most ineffective of the whole series of time points used. The effects on the number and size of distant metastases in the absence and presence of interstitial fluid pressure were also investigated. We performed a local sensitivity analysis (SimBiology^® Model Analyzer) of the model to investigate which PK/PD parameters had a significant influence on the growth behavior of the primary tumor or metastases (P^*, P, Q and N). Since the anti-cancer drug influenced the size of the primary tumor and the number of metastases, the tumor growth response was simulated for each of the synthetic parameter sets. The local sensitivity analysis of the combined model showed that drug clearance k_e (µ = 3.0) and drug potency k₂ (µ = 2.1) are the two most sensitive parameters that affect tumor growth behavior. This is especially true for the oxygenated fraction P since the treatment is only effective in this tissue which is limited due to the low penetration depth when diffusion takes place, but convection is absent. The full sensitivity analysis for the combined model is displayed in the supplementary material Figure S3 which increases the understanding of the relationships between input and output parameters in our model for further investigation. We drew parameter values from a joint probability distribution with a standard deviation of θ = 0.20 for both parameters to cover a wide parameter range. To select a subset of samples, random sampling with a rank correlation matrix, where the matrix was an identity matrix, was performed. Simulations were then performed for each parameter set.

Figure 7 shows the distribution of each parameter and simulation outputs for the primary tumor based on the generated synthetic samples. The influence of the different therapies led to an overall reduction in tumor size (P^*; range −22 to +8%) at the end of the experiment (simulation time 50 days, treatment on day 30). Compared with the drug potency k₂, drug clearance k_e had a considerably greater influence, as reduced clearance meant the drug was present within the system and affecting tumor cells for a prolonged period of time. Furthermore, depending on the length of time the drug was present in the system, the oxygenated tissue recovered faster or slower (Figure 7F), which influenced the overall reduction in tumor size. Drug potency k₂ can also have a substantial impact but this is limited to the oxygenated compartment due to the low diffusion depth. The various tumor size reductions led to different spreading behaviors as well as different size distributions of those metastases that did not switch into a dormant phase. No differences were found in metastases that switched into a late dormant state after reaching a size of 10 to 100 cells. This result indicates that all metastases present at the time of therapy are killed in the absence of interstitial fluid pressure.

4.   Discussions

4.1.   Mathematical modeling of metastatic spread

The mathematical tumor growth model presented here and calibrated with longitudinal tumor data and histological biomarkers successfully describes the development of hypoxic and necrotic fractions reported by Ribba et al. ^[10]. Since no metastases were investigated in the original experiment, we used data from other HT29 xenograft experiments to calibrate the model for the primary tumor to spread and establish distant metastasis ^[24]. Previous research showed that chemotherapy modifies the geometry and density of the blood vessels ^[28]. This change was still detectable 15 days post-treatment compared to the control group, and led to a reduction in the number of malignant cells that were able to spread to distant sites. Even if the correct number of metastases is not known for the underlying experimental data, chemotherapy will have a similar effect on the number of metastases because the increased interstitial fluid pressure limits the distribution of the drug. The simulation results clearly showed the impact of interstitial fluid pressure inside the tumor, which results in a barrier to oxygen, nutrients, and drugs.

4.2.   Build-up of interstitial fluid pressure in metastases

In cancers, oxygen and nutrients penetrate only about 100 µm around the blood vessels in the presence of increased interstitial fluid pressure. This increased interstitial fluid pressure results in the formation of hypoxic and necrotic fractions once the tumor has reached a specific size. In contrast to large primary tumors, micrometastases are too small for interstitial fluid pressure to build up and the blood supply can reach all cells. This is the equivalent of a fractal dimension of 3 for the blood vessels. Micrometastases consist only of oxygenated cells until their size reaches a certain threshold when the fluid pressure starts to build up. The threshold value for the formation of hypoxic and necrotic layers may vary with different blood vessel geometries. Further investigations of the dynamics of hypoxic and necrotic tissues are important in order to better understand this process. This variation in the threshold may have a negligible effect when metastases grow as rapidly as the primary tumor, leading to a superficial blood supply after several days. Furthermore, metastases may switch into a late dormant state for at least 30 days after reaching a size between 10 and 100 cells ^[24]. However, this threshold may be important for simulation scenarios that extend far beyond a period of about 50 days.

4.3.   Chemotherapy resistance of dormant cells

In a previous HT29 xenograft study, the metastases present at the end of the experiment consisted of 10–100 tumor cells in pfp/rag2 mice and only disseminated tumor cells were generally found in rag2 mice. This difference in the number of metastases was caused by a late dormancy phase of at least 30 days when cells entered a reversible cell cycle arrest in the pfp/rag2 group. The emitted cells in the rag2 group underwent an initial dormancy phase before they were able to proliferate. Cell-cycle specific chemotherapy agents such as doxorubicin that target actively dividing cells cannot kill cells that are in a dormant state because they are not cycling ^[36]. In contrast, cell-cycle non-specific drugs kill cells in any phase of the cell cycle. In order to be able to affect both types of cell simultaneously, up to five agents could be combined, including cell-cycle specific and cell-cycle non-specific chemotherapy agents (polychemotherapy) ^[37]. The mathematical model must be modified to model cycling and resting cells in the oxygenated layer to perform simulations of cell-cycle specific treatments ^[38]. Since cell-cycle non-specific chemotherapeutic agents such as cisplatin were used in our simulations, the simulated treatment was able to kill dormant cells as well. Therefore, simulation results obtaining using cell-cycle specific chemotherapy may show reduced treatment efficacy compared with cell-cycle non-specific chemotherapy.

5.   Conclusions

In this article, a mathematical model is presented to model tumor growth and the dynamics of oxygenated, hypoxic and necrotic tumor fractions. This model was combined with a pharmacokinetic-pharmacodynamic model that describes the behavior of a chemotherapeutic drug over time. We used this combined model as a theoretical tool to investigate the effects of limited distribution of chemotherapy into the primary tumor and metastases as a result of interstitial fluid pressure, which means chemotherapy drugs can only penetrate about 100 µm around blood vessels. The number of cells in the oxygenated tissue increased slowly over time. However, the size of the oxygenated tissue fraction as a percentage of the total tumor size decreased significantly over time, thereby reducing the effect of chemotherapy on the entire tumor. This is due to the limited distribution of the drug, which can only reach oxygenated cells as oxygen and drug diffuse together to the cancer cells. Therefore, early treatment is recommended, especially before surgery, to reduce the base number of oxygenated and hypoxic tumor cells and the vessel density to change the spreading behavior of malignant cells. Similar results were obtained in the treatment of metastases, where interstitial fluid pressure also limits the effect of chemotherapy. The timing of the initiation of chemotherapy is an important factor in the success of the treatment, especially when interstitial fluid pressure is increased. Priority should be given to overcoming the convection barrier to achieve a higher penetration depth for chemotherapeutic drugs ^[2]. This approach would increase the chances of a cure for solid tumors as a more substantial proportion of the tumor cells would be reached by chemotherapy.

Acknowledgements

Bertin Hoffmann gratefully acknowledges funding by the European Social Fund through a scholarship from the State Graduate Funding Program (Landesgraduiertenförderung) of the University of Rostock, Germany. The funder did not have any role in the design and conduct of the study, the analysis and interpretation of the data, the decision to publish, or the preparation of the manuscript.

We thank Jenny Cattermole for critical reading and copyediting.

Conflict of interest

All authors declare no conflicts of interest in this paper.