Dynamics of a delayed reaction-diffusion predator-prey model with the effect of the toxins

1.   Introduction

The most severe brain diseases are glioblastoma and grade 4 astrocytoma. These two classes were identified as glioblastoma before 2021, but due to advances in molecular genetics, they have been recently distinguished by the presence or absence of an isocitrate dehydrogenase (IDH) mutation ^[1]. Both of these brain cancers are extremely difficult to treat, with most treatments leading to tumor recurrence ^[2,3]. Fairly recently, Osswald et al. ^[2] have performed mice experiments where they took grade Ⅱ–Ⅳ gliomas (IDH wild type) from human brain tumors and implanted them into mice. By observing the growth of glioma cells in mice, Osswald et al. ^[2] discovered an important structure which they named "tumor microtubes". Tumor microtubes (TMs) describe thin, long protrusions extending from the cell body of a cancer cell. The tumor microtubes contain cell cytoskeleton and can reach up to 5 times the diameter of a cell ^[2]. In 2017, Weil et al. ^[3] performed microsurgery experiments on mice, where they removed a cylindrical lesion from the tissue, to better understand the role of TM with respect to cancer regrowth and response to treatments. In this paper, we develop a mathematical model that can describe the tumor growth after the following treatments: surgical resection, surgical resection paired with targeted therapy, and surgical resection paired with an anti-inflammatory treatment in mice ^[3]. We discuss these treatments in more detail below. We find that wound-healing mechanisms such as growth signaling and the proliferative advantage inferred from the TMs are likely the most important contributors to tumor regrowth, while the TM-induced anisotropy helps to accelerate this process.

1.1.   Tumor microtubes

Tumor microtubes (TMs) can extend up to 500 $\mu$m from the cell body and some TMs can even exceed this length. Further, the number of TMs tends to increase as the tumor progresses ^[2]. For example, on day 20 (postimplantation) most glioblastoma cells have one TM whereas on day 60 (postimplantation) most cancer cells have at least 4 TMs ^[2]. Moreover, TMs have been linked to aid tumor growth, promote spatial spread, and have been shown to allow communication between cells ^[2,3,4]. Communication between cancer cells is possible due to TMs ability to link cancer cells to each other through connexin 43 gap junctions and the gap junctions facilitate cell communication by allowing calcium wave propagation ^[3,4]. Moreover, it has been identified that growth associated protein-43 (GAP-43) is essential for TM formation and TM network functionality ^[2,3,5]. TMs have also been linked to play a role in treatment resistance, where it was shown that the communication network between cancer cells plays a significant role in aiding the repair of the tumor after treatments such as surgery, radiotherapy, or temozolomide chemotherapy ^[3].

In 2017, Weil et al. ^[3] performed detailed mice experiments in order to test glioblastoma treatments. Mice had one of the primary glioblastoma cell lines (S24 or T269) implanted into their brain. The main difference between the S24 and T269 strands is that in vitro, S24 tends to form more TMs ^[3]. After sufficient growth of the tumor inside the mouse brain, Weil et at. ^[3] performed a surgical treatment where they resected a cylindrical volume from the tumor using a syringe with a diameter of approximately 300 $\mu$m. Weil et al. ^[3] examined the behavior of the tumor for 28 days after surgery, and found that the tumor regrows faster and denser inside the lesion. In addition to surgery, Weil et al. ^[3] tested other treatments combined with surgery. One of the treatments performed alongside surgery in ^[3] can be broadly categorized as targeted therapy, which involved specifically inactivating GAP-43 using small hairpin RNA (shGAP-43) to prevent normal formation of the TMs. Weil et al. ^[3] found that this was a successful treatment to delay the regrowth of the tumor in the lesion in comparison to the control treatment (surgery). Furthermore, an anti-inflammatory treatment with dexamethasone (DEX) was tested in ^[3] which was administered daily for 14 days, starting on the day of surgery. Weil et al. ^[3] found that this treatment only allowed for a transient benefit, where the tumor regrew slower on day 7 after surgery in comparison to the control treatment (surgery), but then on day 14 there was no significant difference between the control treatment (surgery) and the combination treatment (surgery combined with DEX).

TMs are able to aid with cell proliferation ^[2] and Weil et al. ^[3] showed that after surgery, nearby cells actually extend their TMs toward the lesion along which cell nuclei can travel. These nuclei can then begin repopulating the lesion. Weil et al. ^[3] also showed that glioblastoma cells can actually move toward the lesion, where on day 5 after surgery, more than half of the glioblastoma cells move toward the lesion. In addition, surgery creates a wound and a wound-healing response is triggered. Wound healing is classically a four-stage process with stages that may overlap: 1) hemostasis, 2) inflammation, 3) cell proliferation, and 4) remodeling ^[6,7]. Inflammation during a wound-healing process has been shown to have the potential to be tumor promoting, since immune inflammatory cells release growth factors and promote cancer cell proliferation and cancer spread ^[8,9]. Moreover, an injury may activate microglia and macrophages which aid with cancer cell survival, proliferation, and invasion ^[10,11]. Hence, in addition to TMs, wound-healing mechanisms may also be promoting tumor growth. We will use our mathematical model to investigate the significance of these processes related to the cancer regrowth data from ^[3].

1.2.   Modeling of brain tumors

According to the data from the World Health Organization (WHO), the incidence of brain cancer is about 4 out of 100,000 people worldwide (www.who.int). Treatments vary widely depending on the stage, size, and location of the cancer, and age and health of the patient ^[12,13]. Most patients would undergo surgery, followed by radiation treatment, but many other treatments such as chemotherapies, targeted therapies, and immunotherapies, are possible too ^[12].

There have been a number of mathematical models proposed to model the growth and spread of glioblastomas. One can use agent-based models to model glioblastoma dynamics to better capture single cell behavior. For example, Gao et al. ^[14] have used a cellular Potts model to capture cell structure, adhesion, and motility as well as the cell heterogeneity of glioblastoma tumors. Reaction diffusion equations have been especially popular to model the growth dynamics for glioblastoma tumors, as the diffusion term can be used to capture the movement of cells and the source term (in the form of exponential or logistic growth) can be used to capture the growth dynamics of glioblastomas ^[15,16]. The simulations from these models have been shown to have good agreement with MRI data of gliomas ^[15,16], where the model in ^[16] can adapt to different cell movement speeds along the white or gray matter tracts in the brain and the model proposed in ^[15] can use Bayesian personalization of model parameters to better fit to patient data. These models can also be adapted to include anatomical boundaries ^[17]. However, the models in ^[15,16] have only employed isotropic diffusion of the glioma cells, meaning that cells have equal probability to move in any given direction. It has been observed that glioma cells tend to move along white matter tracts in the brain, showing that glioma cells can have preferential movement ^[18]. To include preferential movement, one can use a mechanistic model where anisotropic diffusion of cancer cells arises as a byproduct of various mechanistic forces. For example, Colombo et al. ^[19] proposed a fully mechanistic model which incorporated patient specific anisotropy from diffusion tensor imaging (DTI) data and accounted for nutrient availability for the cancer cells. Alternatively, a fully anisotropic version of the isotropic models described above can be formally derived in the form of transport equations ^[20,21]. The anisotropic movement in the transport equation model is captured by the diffusion tensor, which can be formulated to include patient-specific DTI by specifying the orientational distribution of tissue fibers ^[21,22,23].

Recently, a transport equation model was proposed by Hillen et al. ^[24] which incorporated not only the dynamics of glioma cells, but also the TMs. This model includes TM invasion (where the preferential movement of the tips in the anisotropic environment is accounted for), TMs ability to retract and create new TMs, nuclei movement along the TMs, nuclei maturation into glioma cells, and glioma cell division. By comparing the simulated tumor growth to the data found in Osswald et al.'s ^[2] mouse experiments, they validated the model where they showed that the model is able to replicate the trends in the data. Moreover, they showed that under certain assumptions their model can reduce to key models found in previous studies, for example, the simple reaction diffusion model proposed by ^[16].

In this work, we use one of the models in ^[24] for the cancer cells and add a compartment for the healthy cells to describe the experiments of Weil et al. ^[3]. In Figure 1, we illustrate the basic ingredients for our tumor microtube model. Glioma cells are connected via TMs and they extend TMs into the lesion region to guide the oriented movement. Healthy cells are also present and they compete with the cancer cells for space as well as resources. Healthy cells also invade the lesion region at a lower rate. We will use the model to investigate the role of directed TM invasion on cancer regrowth inside the lesion. We expected that the oriented movement guidance of the TM was essential to repair the lesion in a short time. We find, however, that the dominating mechanism is likely an increased growth rate in the lesion caused by increased growth signaling from the wound-healing response as well as the proliferative advantage from the TMs. The TM orientation is a secondary contributing factor.

Figure 1.  Sketch of the basic microtube-driven cancer invasion process. Cancer cells are connected by tumor microtubes and extend microtubes into the lesion region. Healthy cells compete for space and resources and grow into the lesion region at a lower rate.

DownLoad: Full-Size Img PowerPoint

While developing our model, we make use of a radial symmetry assumption. This allows us to use an explicit form of anisotropy as developed in Bica et al. ^[25]. In general, the anisotropy term requires complex calculations of second-order moments of spherical distributions (see ^[26]). The radial case allows us to explicitly formulate anisotropy through a radial anisotropy parameter $\alpha(r)$. In doing this, we combine some theoretical results of ^[24] and ^[25], and apply them to a given situation of glioma lesion experiments of Weil et al. ^[3]. Our approach shows that our modeling is useful for this process. Moreover, we use the model to consider treatments and use it to explain success or failure of the treatments used in ^[3]. In future research, we can use the model to test other treatments and possibly suggest new experiments.

1.3.   Paper outline

In Section 2, we combine the more theoretical approaches of ^[24] and ^[25] to develop a model specifically for the experiments of Weil et al. ^[3], which we then parametrize in Section 2.3. In Section 3, we present the simulation results and relate them to the experiments. We compare isotropic versus anisotropic diffusion and consider smooth or sharp transitions of the parameter functions at the boundary of the lesion. We also perform a sensitivity analysis. Section 4 includes the two treatments that were used in the experiments. The first treatment is targeted therapy that reduces the formation of tumor microtubes paired with surgery. The second treatment is an immunosuppressant treatment with dexamethasone paired with surgery. Our model confirms the finding that inhibition of TM formation is a more effective treatment than immune suppression. We finish with a discussion in Section 5 where we relate our results to current research, explain their significance, and suggest avenues of further research.

2.   The model

In our modeling, we assume that the biological domain is roughly radially symmetric. This enables us to use results from Bica et al. ^[25] to explicitly derive the anisotropic diffusion term in a radial symmetric case.

2.1.   Radially symmetric anisotropic diffusion

Before introducing our model, we first discuss the spatial model that we will use to describe anisotropic (and also isotropic) movement. The domain that we will work with is a two-dimensional circle, which we call $\Omega$. In this circle, the cells have the option to move along the radial lines, perpendicular to them, or simply diffuse (not have any relationship to the radial lines). This movement is sketched in Figure 2, where $\alpha(r)$ is a function depending on the radius, which denotes the degree of movement along the radial lines. If $\alpha(r) > 0$, then the cell moves preferentially along the radial lines and if $\alpha(r) < 0$, then the cell moves preferentially perpendicularly to the radial lines as shown in Figure 2. In the limit $\alpha \rightarrow 1$, cells move fully parallel to the radial lines and in the limit $\alpha \rightarrow -1$, cells move fully perpendicular to the radial lines. If $\alpha(r) = 0$, then the movement is isotropic (i.e., cells diffuse with no correlation to the radial lines). As in Bica et al. ^[25], we assume that $|\alpha(r)|\leq1$.

Figure 2.  Plot of cell movement along radial lines. The orange oval denotes a cell. If $\alpha(r) > 0$, then the cell moves preferentially parallel to the radial line and if $\alpha(r) < 0$, then the cell moves preferentially perpendicular to the radial lines. The lesion is denoted by the filled-in circle filled with blue, where $R_0$ is its radius. The outer circle has radius $R$.

DownLoad: Full-Size Img PowerPoint

To derive the model for the directed movement shown in Figure 2, Bica et al. ^[25] started with the fully anisotropic diffusion equation proposed by Hillen et al. ^[27,28], which is given by

with $u = u(t, x)$ denoting the cell density, $t\in\mathbb{R}$ denoting time, $x\in \Omega$ denoting space, $u_t$ denoting the rate of change of the cell density with respect to time, and $\mathbb{D}(x)$ is a diffusion tensor. Note that this model is not in the standard form of a Fickian diffusion since here we consider actively moving particles and the Fokker-Planck-type diffusion (2.1) is a more appropriate formulation ^[27,28,29].

We transform to polar coordinates $(r, \phi)$ by assuming no-flux boundary conditions at $r = R$:

and assuming that the solutions are bounded at $r = 0$.

The diffusion tensor used by Bica et al. ^[25] was derived by Hillen et al. ^[26] from the von Mises distribution, and the formula for the diffusion tensor is given by

in the planar polar coordinate system $(r, \phi)$. In (2.2), $\alpha (r)$ is the function describing the degree to which cells move along the radial lines as discussed earlier. Substituting (2.2) into (2.1), Bica et al. ^[25] derived the following model:

where $r\in(0, R]$ with $R$ being the maximum radius, and $\Phi[r, \alpha, u]$ denotes the linear, second order differential operator that depends on $r$ and $\alpha(r)$. Note that in the isotropic case $\alpha = 0$, this operator simplifies to the standard Laplacian in polar coordinates in the radially symmetric case with diffusion coefficient 1/2:

The model (2.3) is a one-dimensional PDE model which is applicable to situations in Figure 2 and can be easily solved numerically.

2.2.   Our model

To describe the experiments of Weil et al. ^[3], we inform the choice of $\alpha(r)$ from the presence of tumor microtubes (TMs). We assume, as seen in the experiments of Weil et al. ^[3], that TMs extend in the direction of the lesion, i.e., in the radial direction to the center of the lesion. As cell movement is guided by the TM, this corresponds to a choice of $\alpha(r) > 0$. The value of $\alpha$ then describes the sensitivity of movement to the TM orientation. Later we use 0.6 as a reasonable choice.

We extend the spatial model for two competing populations, cancer cells and healthy cells. Let $u(t, r)\geq 0$ be the density of the glioma cells and $v(t, r)\geq 0$ be the density of healthy cells, where $r > 0$. The model we propose is as follows:

where $\Phi$ is given in (2.3).

We consider the above model on a disk of radius $R > 0$ with boundary conditions

The initial conditions will be defined later once we study specific experiments.

In (2.4), $\alpha(r)$ is a function describing the degree to which cells move along the radial lines as discussed earlier. The parameters $d_u$ and $d_v$ are the respective diffusion coefficients of cancer cells and healthy cells. We use a logistic Lotka-Volterra competition term for growth and interactions where $g_u$ and $g_v$ are the respective growth rates of cancer cells and healthy cells, and $K$ is the carrying capacity ^[30]. Notice that in the second equation for the healthy cells, the spatial term is obtained by setting $\alpha(r) = 0$, since healthy cells do not extend TMs. Essentially, in Model (2.4), we are assuming that glioblastoma cell movement is guided by the radial lines and the healthy cells simply move by diffusing. Both the cancer and healthy cells grow logistically where both cell types compete for the limited space and nutrients which are capped at the carrying capacity. Note that in Model (2.4), we incorporate the wound-healing mechanisms such as increased growth signaling and increased carrying capacity, directly into the model parameters.

2.3.   Parameterization

We take our domain $\Omega$ to be a circle with a radius of $R = 600 \ \mu$m. The domain $\Omega$ will have two distinct regions: inside the lesion (including the lesion boundary) and the region outside the lesion (see Figure 2). We call the lesion and its boundary $\Omega_{in}$ and the region outside the lesion $\Omega_{out}$. The radius for $\Omega_{in}$ is taken to be $R_0 = 150 \ \mu$m since the diameter of the lesion in Weil et al. [3, Figure 1] is approximately $300 \ \mu$m. We assume that the interior of the lesion $\Omega_{in}$ is homogeneous tissue. $\Omega_{out}$ contains blood vessels and other tissue structures, but they are not the focus of this study and we assume they are essentially uniformly distributed.

Due to the TMs and wound-healing mechanisms, the cells growing in the lesion area will have an advantage, therefore requiring us to distinguish between parameters in $\Omega_{in}$ and in $\Omega_{out}$. Hence, to summarize the parameters, we represent them in a tuple $(\cdot \; \cdot)$ in Table 1, where the first entry is the parameter inside $\Omega_{in}$ and the second entry is inside $\Omega_{out}$. To make it clear which parameter we are discussing, we use the notation $(\underline{p}, \bar{p})$ where $\underline{p}$ denotes the parameter inside $\Omega_{in}$ and $\bar{p}$ denotes the parameter in $\Omega_{out}$.

Table 1.  Summary of the parameters used. The tuple denotes parameters within and outside the lesion. The first value contains the parameter inside the lesion, and the second value contains the parameter outside the lesion. In the Reference row, "a'' stands for assumed, and "e'' stands for estimated.

Parameter	Cancer cell	Healthy cell	Cancer cell	Healthy cell	Carrying
	diffusion	diffusion	growth rate	growth rate	capacity
Units	$\mu$m$^2$/h	$\mu$m$^2$/h	1/h	1/h	cells/$\mu$m$^2$
(inside;outside)	$(\underline{d_u};\bar{d_u})$	$(\underline{d_v};\bar{d_v})$	$(\underline{g_u};\bar{g_u})$	$(\underline{g_v};\bar{g_v})$	$(\underline{K};\bar{K})$
Values	$(6.4;6.4)$	$(3.2;1.6)$	$(0.018;0.006)$	$(0.0012;0.0006)$	$(0.01,0.0018)$
Reference	([24];[24])	(a; a)	(a; e)	(a; [24])	(e; e)

---

 | Show Table

DownLoad: CSV

In Table 1, we choose to set the diffusion coefficients $\underline{d_u} = \bar{d_u} = \frac{s_p^2}{2\mu_p}$ as was done in ^[24] where the factor of $1/2$ comes from (2.3). This estimate is based on a typical cell speed of $s_p = 4.8 \ \mu$m/h and a turning rate of $\mu_p = 1.8$/h, which are parameters taken from ^[24]. For $d_v$, we assume that inside $\Omega_{in}$, $\underline{d_u} = \underline{d_v}/2$, and inside $\Omega_{out}$, $\bar{d_u} = \bar{d_v}/4$, because healthy cells move slower than cancer cells ^[8], and since there is more space available in $\Omega_{in}$, the cells can move faster. To obtain $g_u$ for cancer cells, we do a rough estimate based on the figures given in Weil et al. ^[3] which we summarize here. From Figure 1 in ^[3] (also shown in Figure 5, bottom two rows), we see that the rough doubling time in $\Omega_{out}$ is 5 days based on images of day 3 and day 5 in ^[3]. Then using the formula $\text{growth rate} = ln(2)/\text{doubling time}$ ^[31], we get a rough growth rate of $\bar{g_u} = 0.006$/h in $\Omega_{out}$. To get the growth rate for cancer cells in $\Omega_{in}$, we triple the $\bar{g_u}$ from $\Omega_{out}$ as cancer cells grow faster inside the lesion due to the wound-healing and TMs. For the healthy cell growth rate, we assume that $\bar{g_v}$ is the same as the bulk growth rate of a glioma tumor, hence we set it to be $\bar{g_v} = 0.0006$/h ^[24]. Since the cells in $\Omega_{in}$ have a growth advantage due to the wound healing response, we double the growth rate in $\Omega_{in}$ and obtain that $\underline{g_v} = 2\bar{g_v}$. We also estimate the carrying capacity from Figure 1 in ^[3]. For $\bar{K}$ in $\Omega_{out}$, based on day 0 in Figure 1 in ^[3], we can roughly count how many cancer cells fit into a $150 \ \mu$m $\times$ $150 \ \mu$m square. We then take that estimate as $\bar{K}$ for $\Omega_{out}$ for both healthy and cancer cells. For $\underline{K}$ in $\Omega_{in}$, we note that a glioma cell is approximately $10 \ \mu$m from Figure 2 in ^[3] and estimate that 225 cells fit into a $150 \ \mu$m $\times$ $150 \ \mu$m square. From there we obtain $\underline{K}$ in $\Omega_{in}$. Note that we assume a higher carrying capacity in $\Omega_{in}$ due to wound-healing signaling and angiogenesis, and hence cells in $\Omega_{in}$ have more favorable growth conditions ^[32].

For the function $\alpha(r)$ we use the hyperbolic tangent function given by

which is one of the functions for $\alpha(r)$ proposed in ^[25]. We choose it as it allows for a smooth and gradual transition between regions. We set $\rho = 150$ because between the regions $r\in(0,150]$ ($\Omega_{in}$) and $r\in(150,600]$ ($\Omega_{out}$), the degree of anisotropy changes. We reason that the closer the cells are to the lesion, the more radially they move as evidence from ^[3]. In ^[3], they note that nearby cancer cells tend to move toward the lesion. Hence, we set $\beta = 0.3$ which makes the upper bound of $\alpha(r)$ to be 0.6, meaning that most cells inside and near $\Omega_{in}$ move parallel to the radial lines. As the cells are further away from the lesion, we assume they mainly move via diffusion, so $\alpha (r) \rightarrow 0$ as $r\rightarrow R$. We consider two choices of $\gamma$; for a smooth transition we choose $\gamma = 0.06$ and for sharp transition $\gamma = 1$ in $\alpha(r)$ (2.6). We plot the two possible choices of $\alpha(r)$ in Figure 3. We choose the maximal value of $\alpha$ to be 0.6 to express small variations of the TM from the radial direction. Note that for the isotropic case, we have $\alpha(r) = 0$.

Figure 3.  Plot of the smooth and sharp transitions for $\alpha(r)$ from (2.6) where $r$ has the units $\mu$m. The gray shaded region indicates the lesion in these experiments. We use $\beta = 0.3$ and $\rho = 150$ for both curves. For the smooth transition, we use $\gamma = 0.06$, and for the sharp transition, we use $\gamma = 1$ in (2.6).

DownLoad: Full-Size Img PowerPoint

For model (2.4), we use the following initial conditions:

meaning that there are no cells in the lesion, and outside the lesion, the cell density is at carrying capacity. We assume that healthy cells make up $90\%$ while cancer cells make up $10\%$ of the carrying capacity. For the numerical implementation, we simply assume the parameters inside and outside the lesion have different values. We also tried a smooth transition for the parameters, similar to the choice of $\alpha(r)$ as the hyperbolic tangent, but we observed no difference in the numerics.

3.   Results

3.1.   Simulation of the model

To simulate Model (2.4) with the initial conditions (2.7) and boundary conditions (2.5), we use MATLAB's pdepe solver which is suitable for solving parabolic and elliptic PDEs. The parameters are shown in Table 1 and $\alpha(r)$ is as defined in (2.6) applicable to the "smooth" case. The time mesh is taken to be for 0 to 672 hours with a step size of 0.2 and the radial interval is taken to be 0 to 600 with a step size of 0.2. We then obtain the 1D solutions for Days 0, 1, 3, 7, 14, 21, and 28, which are illustrated in Figure 4. The solutions in Figure 4 were then rotated radially to obtain the visualization in $\Omega$ as shown in the first two rows of Figure 5.

Figure 4.  Plots of 1D solutions of Model (2.4) at particular days which illustrate the cancer cell density $u$ (cells/$\mu$m$^2$) with respect to the radius $r$ ($\mu$m). The right graph shows a close-up view of the curves on the left. The shaded region represents $\Omega_{in}$. Model (2.4) was simulated with theparameters given in Table 1, initial conditions (2.7), and boundary conditions (2.5).

DownLoad: Full-Size Img PowerPoint

Figure 5.  Comparison of the simulated tumor density (first two rows) with experimental results (last two rows). The first two rows show plots of 2D solutions of Model (2.4) at particular days where each plot illustrates the cancer cell density $u$ on a circle $\Omega$. The radius of $\Omega$ is 600 $\mu$m and $x, y$ on the axes have units $\mu$m. Model (2.4) was simulated with the parameters given in Table 1, initial conditions (2.7), and boundary conditions (2.5). The last two rows show the experimental images of glioblastoma growth from Weil et al. ^[3]. The dashed circle is the lesion site, in green are the glioblastoma cells, and in blue are the blood vessels (used with permission).

DownLoad: Full-Size Img PowerPoint

In Figure 4, we see that the cancer density $u$ invades the lesion region $\Omega_{in}$ as time progresses and by Day 28, $u$ almost grows to the carrying capacity $\underline{K} = 0.01$ cells/$\mu$m$^2$ which is also illustrated in Figure 5. We chose to visualize these particular days in order to see if we could match the trends seen in Weil et al. [3, Figure 1] (bottom two rows of Figure 5). Indeed, with these chosen parameters we are able to replicate the general growth speed and general shape of the tumor. In Figure 5, from the bottom two rows, it can be seen visually that on day 7 the tumor has started regrowing visually matching the density on day 0. Here, we see a similar effect where on day 7 (in the top two rows), the density inside the lesion starts to match the density outside the lesion. In our simulations, we see that the tumor has grown significantly inside the lesion on day 14, which matches the observations in ^[3]. We note that on this day, the tumor starts growing outside the lesioned area as well. On day 21, the tumor has more than doubled in size, which can also be seen in ^[3]. The tumor keeps growing in size on day 28, almost reaching the carrying capacity $\underline{K} = 0.01$ cells/$\mu$m$^2$ inside the lesion. On days 14, 21, and 28, we see from Figure 4 that the cancer $u$ not only grows inside the lesion, but outside the lesion as well, where $u$ progressively invades the nearby space of the lesioned area. Evidence of this is also shown in ^[3], where the tumor begins growing outside the lesion as well, which is particularly evident on day 21 and onward.

From the above observations, we believe that Model (2.4) does well in replicating the time scale of regrowth as seen in ^[3]. However, due to a lack of quantitative measures, such as the $u$ density per $\mu$m$^2$, we cannot validate the model further.

3.2.   Isotropic and anisotropic diffusion

In this subsection, we study how the degree of anisotropy affects the model. As discussed before, setting $\alpha(r) = 0$ and its consecutive derivatives to zero in (2.3) yield a model with isotropic diffusion. We will study if this causes any significant changes to the dynamics that are seen in Figure 5. Further, it is of interest to examine how the degree of smoothness in the $\alpha(r)$ affects the dynamics. We will study this by simulating Model (2.4) for $\alpha(r)$ given in (2.6) for each case of smooth and sharp transitions. In Figure 6, we see a comparison of what occurs to the cancer cell density $u$ for different choices of $\alpha(r)$ at particular days 3, 14, and 28. For each case, we use the parameters in Table 1.

Figure 6.  Plots of 1D solutions of Model (2.4) at particular days which illustrate the cancer cell density $u$ (cells/$\mu$m$^2$) for each $\alpha(r)$ choice with respect to the radius $r$ ($\mu$m). The shaded region illustrates $\Omega_{in}$. The second row shows a close-up view of the curves in the above row to illustrate the transition region. Model (2.4) was simulated with the parameters given in Table 1, initial conditions (2.7), and boundary conditions (2.5).

DownLoad: Full-Size Img PowerPoint

The blue curves in Figure 6 refer to the sharp $\alpha(r)$ case, the pink curves refer to the smooth $\alpha(r)$ case, and the black curves refer to the isotropic case ($\alpha(r)\equiv 0$). We see that the day 3 dynamics for $u$ are very similar in each case, with only slight differences. The sharp $\alpha(r)$ case has sharper transitions between the cancer density within the lesion and the cancer cell density outside the lesion whereas the transition for the smooth $\alpha(r)$ case is smooth between the regions. We see that the anisotropic cases (smooth and sharp $\alpha(r)$) show that the cancer cell density invades further into the lesion initially than in the isotropic case (seen by the leading edge on the left). The isotropic case shows steady and slower invasion into the lesion, which is also seen by the peak in the black curve on day 14. On day 14, we see that the isotropic cases have completely invaded the interior of the lesion while the isotropic case is still growing. On day 28, there is little difference between the cases. Further, for each case (isotropic and anisotropic), we see that the cancer cells invade the healthy tissue outside the lesion at very similar speeds. This overgrowth was also observed in the experiments in ^[3].

For the sharp case, we observe in the pink profiles in Figure 6 that a local peak forms near the boundary of the lesion at $r = 150$. This peak would yield a ring around the bulk tumor if we were to plot these cases in 2D. The presence of this peak is not surprising, as it was shown that the anisotropic diffusion equation (2.1), with sharp transitions in the diffusion coefficient, can form local peaks at such transition layers ^[33]. Such a ring is not observed in the experiments of ^[3], indicating that the smooth transition for $\alpha(r)$ is a more appropriate assumption.

3.3.   Sensitivity

In addition, we study the dynamics for each case (smooth, sharp, and isotropic $\alpha(r)$) if the parameters in Table 1 vary. The results are summarized in Table 2. We choose to examine the most interesting variations, which are those that may affect the growth of $u$. For all three cases, if the carrying capacity $\underline{K}$ is lowered but does not decrease beyond $\bar{K}$, then there is no effect on the speed of growth of $u$. Note that this trend still holds if $\underline{K}$ is increased. If the growth rate of cancer $\underline{g_u}$ is lowered but not lower than $\bar{g_v}$, then $u$ grows slower, as expected. Similarly, if $\underline{g_u}$ is increased then $u$ grows faster in all cases. If we set the growth rates to be the same for both the cancer cell and healthy cell densities, then $u$ will not grow. In this situation, the healthy cell density $v$ is the dominant density in both regions ($\Omega_{in}$ and $\Omega_{out}$). Finally, increasing the diffusion of $v$, does not impact the growth of $u$ significantly. Note that for this study, the changes in parameters were significant, and were not small perturbations that we will explore when performing the sensitivity analysis in the following section. We did this analysis to see if there were any significant changes in the dynamics between the cases when varying parameters. After performing this analysis, we see that the three cases (smooth, sharp, and isotropic $\alpha(r)$) have very similar overall dynamics, where the anisotropy accelerates the regrowth initially.

Table 2.  Summary of trends from varying parameters for the three cases of $\alpha(r)$ (smooth, sharp, isotropic). Note that when we write $(\underline{p} \downarrow \; \bar{p})$, we do not decrease the parameter past $\bar{p}$, that is, $\underline{p} \downarrow \ > \bar{p}$.

$\alpha(r)$	$(\underline{d_u}; \bar{d_u})$	$(\underline{d_v}; \bar{d_v})$	$(\underline{g_u}; \bar{g_u})$	$(\underline{g_v}; \bar{g_v})$	$(\underline{K}; \bar{K})$	Result in $\Omega_{in}$
Isotropic	-	-	-	-	$(\underline{K}\downarrow \; \bar{K})$	$u$ grows at the same pace
Isotropic	-	-	$(\underline{g_u} \downarrow \; \bar{g_u})$	-	-	$u$ grows much slower
Isotropic	-	-	$(\bar{g_v}; \bar{g_v})$	$(\bar{g_v}; \bar{g_v})$	-	$u$ does not grow
Isotropic	-	$(\underline{d_v} \uparrow \; \bar{d_v} \uparrow$)	-	-	-	$u$ grows at the same pace
Smooth	-	-	-	-	$(\underline{K}\downarrow \; \bar{K})$	$u$ grows at the same pace
Smooth	-	-	$(\underline{g_u} \downarrow \; \bar{g_u})$	-	-	$u$ grows much slower
Smooth	-	-	$(\bar{g_v}; \bar{g_v})$	$(\bar{g_v}; \bar{g_v})$	-	$u$ does not grow
Smooth	-	$(\underline{d_v} \uparrow \; \bar{d_v} \uparrow$)	-	-	-	$u$ grows at the same pace
Sharp	-	-	-	-	$(\underline{K}\downarrow \; \bar{K})$	$u$ grows at the same pace
Sharp	-	-	$(\underline{g_u} \downarrow \; \bar{g_u})$	-	-	$u$ grows much slower
Sharp	-	-	$(\bar{g_v}; \bar{g_v})$	$(\bar{g_v}; \bar{g_v})$	-	$u$ does not grow
Sharp	-	$(\underline{d_v} \uparrow \; \bar{d_v} \uparrow$)	-	-	-	$u$ grows at the same pace

---

 | Show Table

DownLoad: CSV

3.4.   Sensitivity analysis

We perform a sensitivity analysis following the procedure outlined in ^[34]. Note that we adapt the sensitivity analysis for ODEs to PDEs by fixing a representative spatial point inside the lesion ($r = 20$). The sensitivity indices $S^u_p$ for the cancer cell density $u$ and given parameter $p$ is summarized in Figure 7. The sensitivity indices were calculated numerically for a fixed location inside the lesion $r = 20 \ \mu$m at time 28 days according to the relative change ^[34]:

where $p$ is the parameter of interest, $\triangle p$ is 1% of the $p$ value, $u_p(28, 20)$ is the solution with default parameters $p$ at location $r = 20$ and $t = 28$, and $u_{p+\triangle p}$ is the solution with a slightly increased parameter value. We see in Figure 7 that the parameters that have the greatest impact in affecting the cancer cell density are the cancer growth rate $\underline{g_u}$ and carrying capacity $\underline{K}$. This is not surprising as increasing the growth rate allows for the cancer cell density to grow faster, and increasing the carrying capacity allows for the cancer cell density to grow larger. Increasing the growth rate of the healthy cell density $\underline{g_u}$ and increasing the diffusion constant of the healthy cell density $\underline{d_u}$ have a minor negative effect on the growth of the cancer cell density. This is as expected, since faster growing healthy cell density can populate more space decreasing the cancer cell density, and increasing the diffusion constant allows for the healthy cell density to spread more, thereby decreasing the space available for the cancer cells. The sensitivity analysis reveals that increasing the degree of anisotropy (increasing $\beta$ in $\alpha(r)$) slightly decreases the growth of the cancer cell density. This means that moving via isotropic diffusion is slightly more advantageous in increasing cancer cell density. This is what we see in Figure 6, where the isotropic movement of cancer cell density has a more concentrated core. Negative sensitivity to $d_u$ can be explained by noting that higher diffusion will cause more cancer cells to move out from the lesion area $\Omega_{in}$ and those cells lose the additional growth benefit $\Omega_{in}$ provides.

Figure 7.  Summary of the sensitivity indices for the cancer cell density $u$. The parameters in the image refer to the parameters in $\Omega_{in}$, where $\beta$ is used in $\alpha(r)$ (2.6).

DownLoad: Full-Size Img PowerPoint

4.   Treatment

Here, we test two treatments in combination with surgery as was done in ^[3]. The first is a targeted therapy which was performed by injecting a small hairpin RNA (shGAP-43) to inhibit TM formation (by inhibiting GAP-43). We apply this treatment by reducing the growth rate of the cancer cells by 50% in the lesion, and we reduce the value of $\alpha(r)$ by 75% in $\Omega_{in}$ because targeted therapy targets the TMs reducing their ability to form, and therefore hindering directed motion and cancer cell proliferation within the lesion. We see the results of this treatment in Figure 8 (middle row). Note that the colorbar differs from that of Figure 5, in order to see the smaller densities easier. For the targeted therapy treatment, we see that the cancer cell density grows significantly slower than the control, which is expected since the growth rate is much smaller. The same effect was seen in Figure 3 in ^[3] and we are able to recreate the dynamics for this treatment.

Figure 8.  Illustration of different treatments where the colorbar shows the cancer cell density $u$ (cells/$\mu$m$^2$). In each figure, $x$ and $y$ have units $\mu$m. In row 1, Model (2.4) is simulated with the control parameters in Table 1. In row 2, targeted therapy with shGAP-43 is simulated (2.4), and in row 3, the anti-inflammation treatment is simulated by using (2.4) and the parameters outlined in Section 4. For days 3, 7, 14, the colorbar on the left applies and for day 28, the colorbar on the right applies in order to show changes in the densities more clearly.

DownLoad: Full-Size Img PowerPoint

The other treatment we test is the anti-inflammatory treatment with dexamethasone (DEX) ^[3]. DEX is a glucocorticoid (GC) which binds to the GC-receptor (GR) of the cells and deactivates NF-$\kappa$B ^[35,36,37]. NF-$\kappa$B is essential for the homeostasis of the immune system and lower levels of NF-$\kappa$B lead to reduced immune response ^[36]. However, the immune suppression is not permanent and several resistance mechanism are known ^[38,39]. For example, in leukemia, it was shown ^[38] that immune cells exposed to GC reduce the expression of GR over time, consequently becoming less sensitive to glucocorticoids.

We employ the following sigmoidal curve to describe how the model parameters return to their base values as treatment resistance sets in. We define

where the sub-index $x$ stands for either $x = g_u, g_v,$ or $K$, i.e., the two growth rates and the carrying capacity. Here $p_x$ denotes the base value in Table 1 and $a_x$ is the value of the parameter under DEX treatment. The time $t^* = 3.5$ is set to describe the onset of resistance between 3 and 4 days, which was the best choice to describe the experimental observations. See Figure 9 for an illustration of (4.1).

Figure 9.  Plot of (4.1). The horizontal axis is time, which is measured in days, and the vertical axis is the parameter values. $a_x$ denotes the value of the parameter under DEX treatment and $p_x$ denotes the base value in Table 1.

DownLoad: Full-Size Img PowerPoint

In $s_K(t)$, we set $a_K = 0.0018 = \bar{K}$ while $p_K = 0.01 = \underline{K}$. For the growth rates $g_u$ and $g_v$, we assume that both are initially reduced by 50% inside the lesion. Hence $a_{gu} = \underline{g_u}/2$, $p = \underline{g_u}$ and $a = \underline{g_v}/2$, $p = \underline{g_v}$ with $g_u$ and $g_v$ from Table 1.

The model we use inside $\Omega_{in}$ during the DEX treatment is

whereas model (2.4) is unchanged in $\Omega_{out}$. Note that (4.1) could be replaced by any other sigmoid function with very similar results.

The simulations for the DEX treatment are shown in the third row of Figure 8. We can see that the cancer cell density is lower on day 7 than in the control, on day 14, it begins to catch up, and by day 28, it is at the same level as the control. In the corresponding Figure 4 in ^[3], they report that there is no significant difference between the control and the DEX treated tumor after 14 days.

5.   Discussion

In this study, we have proposed a mathematical model that applies to the mice experiments conducted by Weil et al. ^[3] after surgery of glioblastoma tumors. We find that the model is able to replicate the dynamics of glioma tumor growth and spread after surgery. Our model heavily simplifies the wound-healing mechanisms as well as TM dynamics, which serve to elevate the proliferation rates of cells and carrying capacity within the lesion site. The simulations show that the tumor grows back faster and denser within the lesion site matching the observations seen in Weil et al. ^[3]. This shows that the wound-healing processes and the proliferative advantage from the TMs are potentially the key to explaining the faster and denser regrowth that was observed experimentally.

Our model uses an anisotropic model proposed by Bica et al. ^[25] in order to account for the potential directed motion of brain cancer cells. We found that anisotropic movement speeds up the regrowth when compared to the isotropic case. However, with the available measurements, we find that isotropic as well as anisotropic diffusion can explain the data well.

We tested the combined treatments: surgery with targeted therapy (shGAP-43) and surgery with an anti-inflammatory treatment (DEX). Due to a lack of data, we kept the addition in the model simple where treatments predominately affect the growth and carrying capacity parameters. Our model was able to generally match the trends seen in the treatment experiments of Weil et al. ^[3], where the targeted therapy matches the experimental data well but the anti-inflammatory treatment is not quite able to match the speed of regrowth that is seen in the experiments. We find that the treatments that are able to significantly reduce the proliferation rate of cancer cells as well as the carrying capacity for a prolonged period of time are most effective at delaying the growth of the tumor. The targeted therapy treatment performs much better than the anti-inflammatory treatment as it is able to significantly reduce the proliferation rate of cancer cells by inhibiting the TMs whereas the anti-inflammatory treatment only provides a transient benefit. In future studies, it would be of interest to test other treatments numerically. For example, we could study the outcome if chemotherapy or radiotherapy treatment are combined with surgery or with a targeted therapy treatment. Adding chemotherapy and radiotherapy treatments into Model (2.4) is straightforward through extra treatment terms ^[40,41].

Model (2.4) can also be used to describe the initial growth of a tumor or tumor invasions as done in ^[24]. We note that Model (2.4) can be extended to include more aggressive competition between cell types, as follows:

where $c_1$ is the fitness rate of healthy cells and $c_2$ is the fitness rate of cancer cells. That is, if $c_1$ is higher than $c_2$, then healthy cells are more fit than the cancer cells and therefore, more easily grow in the domain (and vise versa). We analyzed Model (5.1) but we found that Model (2.4) is sufficient at explaining the results of Weil et al.'s ^[3] experiments.

Our study has several limitations. We assumed that the domain outside the lesion site is homogeneous, where blood vessels are ignored. It is possible that the cancer cells would prefer to grow near the blood vessels to have easier access to nutrients which could explain the asymmetrical tumors that are seen in the experimental images. Further, we did not include a full immune response into the model. In order to obtain biologically realistic parameters for the model extension, more data needs to be collected, in particular, the movement and activity of immune cells. If more data is available, more accurate effects of treatments could be incorporated into the model. Further, more complicated interactions between healthy and cancer cells could be studied, such as the competition for nutrients. We refrained from putting in extra processes in order to keep the model simple and to study the main processes.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

AS acknowledges the funding from the Alberta Graduate Excellence Scholarship. TH is supported through a discovery grant of the Natural Science and Engineering Research Council of Canada (NSERC), RGPIN-2023-04269.

Conflict of interest

All authors declare no conflicts of interest in this paper.