Accelerated tumor invasion under non-isotropic cell dispersal in glioblastomas

Glioblastomas (GLBs) are the most common variety of brain tumors and are highly malignant (median survival <15 months [1, 2]). They are known to display diffuse boundaries (unlike most solid tumors) and are very difficult to treat. Surgery often ends in recurrence, i.e. the regrowth and invasion from surviving cancerous cells. Such a pattern is seldom modified either by chemotherapy or by cytotoxic drugs. Over the last decade, biomedical efforts have been running in parallel with mathematical and computational approaches. Many of these efforts have considered the nonlinear character of invasion phenomena by cancer cells and its implications for tumor spreading [3, 4]. The problem has been addressed from a physics perspective approximation [5–9] where GLBs can be roughly described as a growing, isotropic mass of tumor cells replicating under limited resources. This can be approached through a model based on a Fisher equation, namely

$$\begin{equation} {\frac{\partial \rho }{\partial t}}=g\rho \left( 1-{\frac{\rho }{K}}\right) +\nabla \cdot (D(\mathbf{r})\nabla \rho ), \end{equation} \tag{ 1 }$$

where ρ(r,t) is the local density of cells at a given point r and time t, whereas g and K are the growth rate and maximal population density, respectively. The diffusion coefficient is a spatially dependent one, i.e. D = D(r), and allows us to take into account the heterogeneous structure of brain tissue [8]. This model is able to provide a good picture of the spatial tumor dynamics, and its outcome can be compared with clinical data involving the spatial and temporal dynamics of tumors. In figure 1(a) we see a typical example of a simulated GLB. The bright area corresponds to the dense domain visible through brain scans, whereas the black domain corresponds to the infiltrated tumor zone, not detectable in brain scans but naturally expected to be present if a decaying diffusion front is involved. If we ignore spatial heterogeneity, assuming a constant diffusion coefficient, the key prediction from equation (1) is Fisher's propagation speed,

$$\begin{equation} c_{\rm Fisher}=2\sqrt{gD}, \end{equation} \tag{ 2 }$$

for a radially traveling front (see figure 1(b)). Alternatively, Fisher's model predicts that the front width is proportional to $\sqrt {D/g}$ [10] and this gives a measure of the size of the infiltrated tumor region [10]. Fisher's model is very important in order to estimate the relative fraction of infiltrative versus total tumor mass [10], GLB invasion speeds [8], cancer recurrence [11] and tumor shape and spreading [12]. Fisher's model assumes that cells behave as random walkers, but several aspects of cancer cell population dynamics depart from simple diffusion [13]. In particular, it has been shown that some GLB cell lines display cell movements biased toward the local invasion direction (i.e. the outward radial direction in the case of a circular tumor) [14, 15]. Such a biased pattern should have a predictable impact on invasion speed. In order to provide an analytic estimate of this impact, here we present several increasingly realistic models of GLB spread with an outward dispersal preference by invasive cells.

Zoom In Zoom Out Reset image size

We want to describe the dynamics of a growing tumor. Below we will introduce the necessary dispersal parameters, which have been estimated from cell trajectories recorded experimentally in two dimensions [14, 15]. Accordingly, we consider a two-dimensional model and assume the tumor has a circular shape, with growing radius centered at position (x,y) = (0,0). For simplicity we consider only invasive cells, i.e. the tumor core is not described by our model. Therefore, the invasion speed considered is that of the outer rim of the invasive region [14]. This approach allows us to focus our attention on the effect we are interested in, i.e. that of the biased dispersal of invasive cells.

Let ρ(x,y,t) stand for the number of invasive cells per unit area at position (x,y) and time t. If cell positions are measured at times separated by the time interval T, we define the dispersal kernel ϕ(Δ_x,Δ_y) as the probability per unit area that a cell that was at (x − Δ_x,y − Δ_y,t) moves to (x,y,t). Let R[p(x,y,t)] stand for the number of new cells due to net reproduction over a time interval T. From these definitions, it follows that [16]

$$\begin{eqnarray} &&\fl \rho (x,y,t+T)-\rho (x,y,t)=\int_{-\infty }^{+\infty }\int_{-\infty }^{+\infty }\rho (x-\Delta _{x},y-\Delta _{y},t)\;\phi (\Delta _{x},\Delta _{y})\;\mathrm{d}\Delta _{x}\;\mathrm{d}\Delta _{y}\nonumber\\ &&-\rho (x,y,t)+R\left[\rho (x,y,t)\right] , \end{eqnarray} \tag{ 3 }$$

where the first two terms on the right-hand side correspond to cells arriving and leaving an area centered at (x,y), respectively.

This approach is well known [14, 17], but it will be useful to analyze its consequences before generalizing it (in section 3). Performing a Taylor expansion up to second order in space and first order in time, we obtain from equation (3)

$$\begin{equation} \frac{\partial \rho }{\partial t}=\sum_{\mu =x,y}\left( -U_{\mu }\frac{\partial \rho }{\partial \mu }+D_{\mu }\frac{\partial ^{2}\rho }{\partial \mu ^{2}}\right) +D_{xy}\frac{\partial ^{2}\rho }{\partial x\textrm{ }\partial y}+F(\rho ), \end{equation} \tag{ 4 }$$

where we have defined the terms U_μ,D_μ (with μ = x,y) and D_xy as follows:

$$\begin{equation} U_{\mu }(x,y)=\frac{1}{T}\int_{-\infty }^{+\infty }\int_{-\infty }^{+\infty }\Delta _{\mu }\;\phi (\Delta _{x},\Delta _{y})\ \mathrm{d}\Delta _{x}\;\mathrm{d}\Delta _{y}=\frac{\langle\Delta _{\mu }\rangle}{T}, \end{equation} \tag{ 5 }$$

$$\begin{equation} D_{\mu }(x,y)=\frac{1}{2T}\int_{-\infty }^{+\infty }\int_{-\infty }^{+\infty }\Delta _{\mu }^{2}\;\phi (\Delta _{x},\Delta _{y})\ \mathrm{d}\Delta _{x}\;\mathrm{d}\Delta _{y}=\frac{\langle\Delta _{\mu }^{2}\rangle}{2T}, \end{equation} \tag{ 6 }$$

$$\begin{equation} D_{xy}(x,y)=\frac{1}{T}\int_{-\infty }^{+\infty }\int_{-\infty }^{+\infty }\Delta _{x}\;\Delta _{y}\textrm{ }\phi (\Delta _{x},\Delta _{y})\ \mathrm{d}\Delta _{x}\;\mathrm{d}\Delta _{y}=\frac{\langle\Delta _{x}\;\Delta _{y}\rangle}{T}. \end{equation} \tag{ 7 }$$

$F\left ( \rho (x,y,t)\right )$ in equation (4) is the time derivative of ρ(x,y,t) due to net reproduction (population growth), i.e. [16]

$$\begin{equation} R\left[ \rho (x,y,t)\right] =T\!F\left( \rho (x,y,t)\right) +\frac{T^{2}}{2!}\frac{\partial }{\partial t}\left. F\left( \rho (x,y,t)\right) \right\vert _{g}+\cdots, \end{equation} \tag{ 8 }$$

where the subindex g stands for population growth [18]. For GLB tumors a logistic growth function is appropriate, i.e. [8]

$$\begin{equation} F\left( \rho (x,y,t)\right) =g\rho (x,y,t)\left( 1-\frac{\rho (x,y,t)}{K}\right) . \end{equation} \tag{ 9 }$$

Front speeds in two dimensions can be found most easily by considering a small area on the outer rim of the invasive region (moving attached to the front), and choosing the x-axis along the local invasion direction so that (within this small area) the front curvature becomes very small for large values of time, i.e. the y-dependence of ρ(x,y,t) can be neglected [19]. Within such a local area, the reaction–diffusion–advection (RDA) equation (4) with (9) simplifies to $\frac {\partial \rho }{\partial t}=-U_{x}\frac {\partial \rho }{\partial x}+D_{x}\frac {\partial ^{2}\rho }{\partial x^{2}}+g$ $\rho (x,y,t)\left ( 1-\frac {\rho (x,y,t)}{K}\right )$. Assuming as usual [19] that the minimum speed is that of the front (this is called the marginal stability criterion [20]), the speed of front solutions to this RDA equation is given by the well-known result (see e.g. [17, p 5])

$$\begin{equation} c_{\rm RDA}=2\sqrt{gD_{x}}+\frac{\langle \Delta _{x}\rangle }{T}. \end{equation} \tag{ 10 }$$

For isotropic dispersal, $\langle \Delta _{x}\rangle =0,$ D_x = D_y = D ≡ 〈Δ²_x + Δ²_y〉/(4T) [16] and we recover Fisher's speed (2). However, Stein et al [14, 15] observed a clear bias toward the GLB radial invasion direction by tracking many individual cells in vitro, measuring their motion in the radial direction (which corresponds to the x-direction in our model above), and averaging over many cells. In this way, they obtained [15] $\frac {\langle \Delta _{x}^{2}\rangle }{T}=1.8\times 10^{-4}\,{\mathrm { cm}}^{2}\,{\mathrm { day}}^{-1}$ or $D_{x}=\frac {\langle \Delta _{x}^{2}\rangle }{2T}=0.9\times 10^{-4}\,{\mathrm { cm}}^{2}\,{\mathrm { day}}^{-1}$ (see equation (6)) for the radial diffusion coefficient, and $U_{x}=\frac {\langle \Delta _{x}\rangle }{T} =0.012\,{\mathrm { cm}}\,{\mathrm { day}}^{-1}$ for the asymptotic advective speed. Additional measurements indicate that g < 0.3 day⁻¹ for the proliferation rate [14], so in figure 2(a) we have considered the range 0.001 < g < 0.3 day⁻¹. The curve labeled c_RDA in figure 2(a) is the total GLB invasion speed according to equation (10). The horizontal dashed line in figure 2(a) is the term due to biased dispersal, i.e. the last term in equation (10). The curve labeled c_Fisher in figure 2(a) is the first term in equation (10). Figure 2(a) also includes, as a hatched area, the experimental speed range in vitro according to the data by Stein and co-workers, namely 0.0067 < c < 0.0133 cm day⁻¹ (see supporting information, available from stacks.iop.org/NJP/15/055001/mmedia). The most important point from figure 2(a) is that, for all values of g in the range from [14], the bias contribution to the invasion speed, i.e. $\frac {\langle \Delta _{x}\rangle }{T}$, is larger than the speed computed according to the usual formula, namely c_Fisher [8]. Therefore, the bias effect cannot be neglected. Finally, although the value $\frac {\langle \Delta _{x}\rangle }{T}=0.012\,{\mathrm { cm}}\,{\mathrm { day}}^{-1}$ used above is realistic (because it has been determined from experiments), the dispersal bias effect on cell trajectories may depend on the cell line considered [14], i.e. different cell lines will display different values of $\frac {\langle \Delta _{x}\rangle }{T}$. Thus in figure 2(b) we plot the classical (Fisher) speed c_Fisher and the biased-model speed c_RDA given by equation (10) as a function of the radial bias parameter $\frac {\langle \Delta _{x}\rangle }{T}$ (for the value g = 0.1 day⁻¹, which lies in the range 0 < g < 0.3 day⁻¹ [14]). We observe that, according to the RDA model, within the experimentally observed speed range (hatched area in figure 2(b)), the bias effect (dotted line) is between 90 and 180%. This confirms that the effect of biased cell dispersal should be taken into account to understand GLB invasion rates.

Zoom In Zoom Out Reset image size

The RDA model predicts a bias effect of about 100% or higher for realistic parameter values (figure 2(b)). But in that model, the bias effect is taken into account by the U_μ-terms in equation (4), which are in turn first-order terms of a Taylor expansion of equation (3). If first-order terms introduce such a substantial effect, then higher-order terms are also likely important. Therefore, in order to make accurate predictions, we have to find the speed of front solutions of equation (3) without performing any Taylor expansion (i.e. we deal with an integro-difference rather than with an RDA equation). Using equation (9) and linearizing yields $\left . \frac {\partial F}{\partial t}\right \vert _{g}=\frac {\mathrm {d}F}{\mathrm {d}\rho }\left . \frac {\partial \rho }{\partial t}\right \vert _{g}=\frac {\mathrm {d}F}{\mathrm {d}\rho }F\simeq g^{2}\rho$ [18] and similarly $\left . \frac {\partial ^{n}F}{\partial t^{n}}\right \vert _{g}\simeq g^{n+1}\rho$. Thus equation (8) yields $R [ \rho (x,y,t)] \simeq ( \mathrm {e}^{gT}-1) \rho .$ We use this in equation (3) and look as usual [19] for constant-shaped solutions with the form $\rho (x,t)=\rho _{0}\exp [-\lambda (x-ct)]$ with c > 0 and λ > 0 (the y-dependence of ρ(x,y,t) can be neglected locally if we choose the local x-axis along the front propagation direction [19]; see the text below equation (9)). This yields

$$\begin{equation} \mathrm{e}^{\lambda cT}-\mathrm{e}^{gT}=\int_{-\infty }^{+\infty }\int_{-\infty }^{+\infty }\mathrm{e}^{\lambda \Delta _{x}}\;\phi (\Delta _{x},\Delta _{y})\;\mathrm{d}\Delta _{x}\;\mathrm{d}\Delta _{y}-1. \end{equation} \tag{ 11 }$$

In order to solve this equation for the front speed c > 0, we need an expression for the dispersal kernel $\phi (\Delta _{x},\Delta _{y})$. We consider three models.

3.1. Dirac-delta model

In this simple model, all cells move the same distance d along the radial direction (recall that this corresponds to the local x-direction in our model), with probability $p>\frac {1}{2}$ to move outwards from the tumor core (i.e. with Δ_x = d > 0) and probability 1 − p to move inwards (Δ_x = −d < 0). We also assume that the jump displacements along the x (radial) and y directions are independent of each other. Then, the dispersal kernel can be written as the product of two probability distributions,

$$\begin{equation} \phi (\Delta _{x},\Delta _{y})=\psi (\Delta _{x})\textrm{ }\chi (\Delta _{y}), \end{equation} \tag{ 12 }$$

where $\int _{-\infty }^{+\infty }\chi (\Delta _{y})\;\mathrm {d}\Delta _{y}=1$ and

$$\begin{equation} \psi (\Delta _{x})=p\textrm{ }\delta (\Delta _{x}-d)+(1-p)\textrm{ }\delta (\Delta _{x}+d), \end{equation} \tag{ 13 }$$

with $\delta (\Delta _{x}-d)$ the Dirac delta centered at $\Delta _{x}=d$. Then equation (11) becomes

$$\begin{equation} \mathrm{e}^{\lambda cT}-\mathrm{e}^{gT}=p\textrm{ }\mathrm{e}^{\lambda d}+(1-p)\textrm{ }\mathrm{e}^{-\lambda d}-1. \end{equation} \tag{ 14 }$$

From equation (14) we see that, even for this very simple kernel, we cannot obtain an explicit equation for the GLB invasion speed c. Thus we use numerical methods (see supporting information, available from stacks.iop.org/NJP/15/055001/mmedia) to obtain the GLB speed as a function of the dispersal bias probability p (curve c₁ in figure 3). In order to compare the Dirac-delta model c₁ to the less accurate RDA model (section 2), in figure 3 we have also included the speed c_RDA from equation (10) by taking into account that, according to equations (5) and (12)–(13),

$$\begin{equation} \langle \Delta _{x}\rangle \textrm{ }=(2p-1)d. \end{equation} \tag{ 15 }$$

We see from figure 3 that the second-order approximation (c_RDA, from equation (10)) substantially overestimates the speed c₁ (and the observed range) if the dispersal bias is relevant (p > 0.8), by up to 30%. This shows that, in general, it is necessary to take into account higher-order terms, as done by the model in this subsection (c₁ in figure 3, from equation (14)). We also see in figure 3 that the model in this subsection (c₁) is consistent with the experimental range (hatched region) for almost the full range of values of the dispersal bias probability p. Finally, as an example we can consider the experimental value $\langle \Delta _{x}\rangle$ =0.012 cm day⁻¹ [15], already used in section 2. According to equation (15), this corresponds to p = 0.95, which is very close to the maximum possible value p = 1 for the bias probability p. For this value p = 0.95, according to figure 3 the higher-order effect (relative to the RDA approximation) is 26%. Moreover, the higher-order prediction c₁ is within the experimental range (hatched region) whereas the RDA prediction is inconsistent with the experimental range. This confirms the need to take higher-order biased-dispersal terms into account. In other words, for some real GLBs the RDA approximation (section 2) is not appropriate.

Zoom In Zoom Out Reset image size

3.2. Gaussian model

The Dirac-delta model above assumes that all cells move approximately the same distance. A more realistic model, widely used in microbiology [21], assumes that the dispersal probability is a Gaussian. For biased dispersal,

$$\begin{equation} \psi (\Delta _{x})=\left\{ \begin{array}{@{}l@{}} \frac{p}{\alpha }\exp \left[ -\left( \frac{\Delta _{x}}{\alpha }\right) ^{2}\right] \qquad \;{\rm if}\quad\Delta _{x}\geqslant 0, \\ \frac{(1-p)}{\beta }\exp \left[ -\left( \frac{\Delta _{x}}{\beta }\right) ^{2}\right] \quad {\rm if} \quad\Delta _{x}<0. \end{array} \right. \end{equation} \tag{ 16 }$$

which used in equations (11) and (6) yields

$$\begin{equation} \mathrm{e}^{\lambda cT}-\mathrm{e}^{gT}=p\,\mathrm{e}^{\frac{\lambda ^{2}\alpha ^{2}}{4}}\mathrm{fcer}\left[ -\frac{\lambda \alpha }{2}\right] +(1-p)\,\mathrm{e}^{\frac{\lambda ^{2}\beta ^{2}}{4}}\mathrm{fcer}\left[ \frac{\lambda \beta }{2}\right] -1, \end{equation} \tag{ 17 }$$

where fcer[] is the complementary error function. In this model, we can avoid an unphysical probability jump by requiring the continuity of the kernel (16) at $\Delta _{x}=0.$ This implies that β = (1 − p)α/p. Using this condition and the kernel (16) in equation (6),

$$\begin{equation} \alpha =\sqrt{\frac{4D_{x}T}{p+\frac{(1-p)^{3}}{p^{2}}}}. \end{equation} \tag{ 18 }$$

In figure 3 we include the speed predicted by this model (c₂), obtained by computing the characteristic distance α from equation (18) for the value of p considered (and the experimental values D_x = 0.9 × 10⁻⁴ cm² day and T = 1 day in section 2), and finding out the minimum speed c for which equation (11) has a real solution (see supporting information, available from stacks.iop.org/NJP/15/055001/mmedia).

3.3. Laplacian model

A decreasing exponential with increasing distance (Laplacian kernel) is reasonable according to some experimental data for glioma cells (see e.g. figure 2B in [22]),

$$\begin{equation} \psi (\Delta _{x})=\left\{\begin{array}{@{}l@{}} \frac{p}{\alpha }\exp \left[ -\frac{\Delta _{x}}{\alpha }\right] \quad {\rm if}\quad \Delta _{x}\geqslant 0, \\ \frac{(1-p)}{\beta }\exp \left[ \frac{\Delta _{x}}{\beta }\right] \quad {\rm if}\quad \Delta _{x}<0. \end{array} \right. \end{equation} \tag{ 19 }$$

Combined with the kernel (19) and equation (11), this leads to

$$\begin{equation} \mathrm{e}^{\lambda cT}-\mathrm{e}^{gT}=\frac{p}{1-\lambda \alpha }+\frac{1-p}{1+\lambda \beta }-1. \end{equation} \tag{ 20 }$$

Continuity at $\Delta _{x}=0$ implies again that β = (1 − p)α/p, which used in equations (19) and (6) yields

$$\begin{equation} \alpha =\sqrt{\frac{D_{x}T}{p+\frac{(1-p)^{3}}{p^{2}}}}. \end{equation} \tag{ 21 }$$

We compute the speed by following the same approach as for model 2 above (and requiring that $\lambda <\frac {1}{\alpha }$, so that the integral in equation (11) is finite). In figure 3 we observe that the results for the Laplacian kernel (c₃) are almost the same as for the Gaussian kernel (c₂) and rather similar to those from the Dirac-delta kernel (c₁). The three models are widely consistent with the experimental range of GLB invasion (shaded region). Thus, interestingly, the kernel shape does not seem to affect our conclusions. Figure 3 also confirms the breakdown of the RDA model, which leads to large errors for real GLBs, e.g. 26% (see model 1 above). This shows the need for higher-order biased models of GLB invasion (as first reported in this work).

We have shown that (i) for some GLBs displaying biased diffusion, it is necessary to apply higher-order models (beyond second order, i.e. beyond the RDA approximation), and (ii) such models are consistent with experimental data. These conclusions have been reached by estimating invasion speeds and parameter values from experimental work by other authors. This raises some interesting points, which can motivate future research. For example, in [15] an outward bias in the motion of cells was observed. However, the bias was stronger during the first day ($U_{x}=\frac {\langle \Delta _{x}\rangle }{T}=0.048\,{\mathrm { cm}}\,{\mathrm { day}}^{-1}$) than during the second day ($U_{x}=\frac {\langle \Delta _{x}\rangle }{T}=0.012\,{\mathrm { cm}}\,{\mathrm { day}}^{-1}$). In the absence of experimental data for subsequent days, it is reasonable to use the value for the second day (as done in this paper), but it would be interesting to determine experimentally whether this value remains constant later on or not (GLB speeds were measured during 6 days in [14, figure 2A]). In turn, this point has implications concerning other parameter values, as follows. In case the bias (as estimated from the value of $U_{x}=\frac {\langle \Delta _{x}\rangle }{T}$) did not reach a constant value after several days, this would imply non-constant values for the dispersal bias probability p and/or the characteristic dispersal distance (d for the Dirac-delta model, see equation (15); α for the Gaussian and Laplacian models). Clearly, it would be very useful to extend cell tracking experiments [15] to longer times, and estimate all of these parameter values during many days. In case some parameter did not reach a fairly constant asymptotic value, analytical approaches would probably fail and numerical simulations would be the only way to compute the GLB invasion speed.

Since different cell lines display different invasion speeds [14], it would be very interesting to determine experimentally if the observed range of the GLB invasion speed (hatched area in figure 3) actually depends on the value of p, and especially if the average speed increases with increasing values of p (as predicted by the biased (non-Fisher) models in figure 3). Unfortunately, GLB cell-tracking experiments have not been performed yet for different cell lines [15], so it is not possible to determine the value of p for cell lines with different invasion speeds [14]. Thus, further experimental work will be necessary before such open questions can be tackled.

We have shown that, for realistic parameter values, the speed of GLB invasion is substantially affected by the fact that the motion of cells is not, in general, equally probable in all directions but can be biased (i.e. more probable) in directions closer to that of the invasion front. Indeed, such a bias has been observed experimentally [14, 15]. We have derived GLB invasion speeds when this effect is present, and found them to be consistent with experiments. Our analytical results make it possible to predict the importance of the bias effect very easily from the values of three relevant parameters (g, D_x and $\frac {\langle \Delta _{x}\rangle }{T}$), without the need to perform numerical simulations. The effect of biased dispersal is substantial, and the RDA approximation seriously overestimates GLB invasion speeds.

Future extensions should consider the non-circular shape of real tumors [23], the shedding of invasive cells by the tumor core [14], the effects of radiotherapy and surgery [23], the impact of potential treatment affecting cell migration (considered as targets for future therapies [24]), the heterogeneous nature of brain structure (particularly in terms of gray and white matter [8]) and the role of stem cells [25] as the originators of gliomas.

In future models, biased dispersal will surely lead to wider predicted infiltrative areas (not detectable by present imaging techniques), and this could improve predictions for the optimal level of tumor resection in order to decrease the probabilities of tumor recurrence after surgery [12].

Finally, it is important to note that simulations of GLB invasion are nowadays currently performed for individual patients (in order to predict their survival times and the effects of several therapies) by making use of parameter values measured for the individual patient considered [12]. Our work allows us to predict the effect of biased dispersal in terms of the values of three relevant parameters (g, D_x and $\frac {\langle \Delta _{x}\rangle }{T}$) that can be measured for each patient. This bias effect is currently not taken into account in such patient-specific simulations [12] but can drastically change their predictions and, therefore, lead to potentially crucial clinical implications.

This work was supported by MEC (FIS-2012-31307 grant, JF), the Generalitat de Catalunya (2009-SGR-374 grant, JF), the Marcelino Botín Foundation (RVS) and the Santa Fe Institute (RVS).