Patient-specific parameter estimates of glioblastoma multiforme growth dynamics from a model with explicit birth and death rates

1.   Introduction

Glioblastoma multiforme (GBM) is a highly agammaessive primary brain cancer, with median survival time from diagnosis on the order of 15 months; long-term survival is extremely rare ^[1]. Such rapid progression is promoted by highly proliferative and diffusely invasive cancer cells, which makes complete surgical removal impossible. Magnetic resonance imaging (MRI) is conventionally used to identify the location and characteristics of the tumor pre-operatively, to guide surgery, and to monitor and track progression and treatment response. Perioperatively, MRI is used to guide the resection of the tumor mass, to assess post-operatively the volume of tumor resected, and to target other adjunct treatment such as radiation therapy.

GBMs morphologically typically appear (at least at initial diagnosis) as roughly spherical but highly heterogeneous masses that often exhibit a (crudely speaking) three-layer structure. Within the tumor there is usually extensive cell necrosis, often accompanied by tumor cells, and a cystic component as well. An outer region, which typically appears as contrast-enhancing on T1-weighted gadolinium contrast-enhanced MRI, is cytologically typified by proliferating cells that then infiltrate into surrounding brain tissue. The surrounding brain tissue is generally seen to be edematous on T2-weighted or T2-FLAIR MRI and at surgery, due to vasogenic edema. Prior statistical analyses have found that edema is a prognostic indicator of patient survival ^[2], but the relationship is complex and appears to be mediated by the expression of vascular endothelial growth factor and the activity of related angiogenic genes ^[3] and various autocrine factors ^[4].

The standard of care for GBM patients was largely established by the 2005 clinical trial by Stupp et al. ^[5]. It comprises maximal surgical resection of the primary tumor, followed by six weeks of radiation to the gross tumor volume, plus a 2–3 cm margin, with concomitant oral temozolamide (TMZ), and 6–12 months of maintenance TMZ chemotherapy ^[6]. Maximal surgical resection appears to offer some survival benefit. Nevertheless, the absolute survival benefit of even the most effective therapy is typically on the order of months. The highly infiltrative nature of GBMs makes recurrence nearly inevitable, even with maximal resection and agammaessive adjuvant therapy, although individual tumors vary in their degree of invasiveness.

Given the grim situation, mathematical modeling has been proposed as a method to better understand the biophysical rules underlying GBM growth, with the ultimate goal to provide more effective therapy. Mathematical models have been widely applied to a variety of cancers and to cancer treatment in general ^[7], and GBM is the focus of many such works (see ^[8] for a review). A popular class of cancer models takes the form of a system of reaction-diffusion equations. In many cases ^[9,10,11], such systems generate a traveling-wave solution, with the traveling-wave speed of great interest, as it is an indicator of how fast the cancer progresses.

Variants of the Fisher-Kolmogorov equation ^[12], originally introduced in the 1930s, were first suggested (to our knowledge) as models for GBM growth by J. D. Murray and coworkers. The Fisher-Kolmogorov model is given by

where $c(x, t)$ is the cancer cell density at location $x$ and time $t$, $D$ is a diffusion coefficient, $\rho$ is the intrinsic tumor cell growth rate, and $K$ is the local carrying capacity. (Variants include a linear version that replaces the logistic growth term with a simple exponential growth rate, $\rho c$.) Murray and coworkers have used them to explore the effect of chemotherapy ^[13], to quantify patients' survival as a function of the extent of surgical resection ^[14], and to estimate the time of tumor initiation ^[15].

Other authors have suggested that the net growth and diffusion parameters of model (1.1) may be estimated by image differencing when two sequential, pre-treatment patient MR series are available ^[16,17,18]. Such a procedure is problematic, however, because changes in the tumor in images taken a few days or weeks apart tend to be small and are convolved with image co-registration errors ^[19]. Some patients may be treated with steroids following initial diagnosis to reduce tumor-related edema and resulting neurological symptoms, which may alter the brain geometry and imaging appearance of the tumor at subsequent times ^[20,21].

Model (1.1) can yield a dense tumor core with an advancing front, but it cannot capture the heterogeneity between live and necrotic tumor cells, as it assumes that all cells are equally viable. While several modeling efforts have taken into account various proliferating, migrating, and necrotic cell components (e.g., ^[22,23]), they are too complicated to be reliably parameterized by the limited number of patient MRI series in typical clinical cases. The motivation for this work is to extend model (1.1) to include necrotic cells in a simplified way, such that patient-specific model parameters can be estimated from suitable measurements of MR images acquired at a single time point.

T1-weighted MRI sequences of GBM often show a partially necrotic core surrounded by a bright enhancing rim that correlates with high blood vessel density and, presumably, with rapid cell proliferation. Neurosurgical and biopsy studies indicate that this core and rim are usually surrounded by a large expanse of edema, which is best visualized on T2-weighted MRI and has been found to correspond with a component of diffusely invasive GBM cells ^[24]. By approximating the tumor as a sphere, we may be able to identify three idealized digital marks from imaging: necrotic radius, enhancing radius, and what we shall call the "T2" or "maximum" radius. We hypothesize that a relatively simple mathematical model framework can capture all these three digital marks and yield insights into the relative contributions of cellular proliferation, motility, and necrosis to the observed image features.

The next section describes our model and its assumptions. We demonstrate that the model has a traveling-wave solution and present the approximate wave profile. We describe a simple procedure to estimate patient-specific parameters by fitting the approximate wave profile to a tumor profile derived from patient MRIs. The identifiability of the model parameters is also discussed. We apply this parameter estimation procedure to obtain the key model parameters (consisting of the rate of cancer cell proliferation, death, and diffusion) for several patients.

2.   Model and method

2.1.   Model description

Our proposed model of the growth of GBM is a system of reaction-diffusion equations:

where

and $p(x, t)$ and $q(x, t)$ represent the proliferating and quiescent cell densities at time $t$ and location $x$, respectively; quiescent cells are functionally equivalent to necrotic cells in this framework. We assume that the flux of total population due to migration is $-D\nabla(p+q)$, where $D$ is a constant diffusion coefficient. It is further assumed that the proportion of the total flux contributed by each cell type equals its proportion of the total population. This form of diffusion was used in ^[25] to account for the key property of contact inhibition in cancer cell movement, with the underlying assumption that the two cell populations move together with equal motility, unaffected by necrotic cells. This type of model has successfully captured the structure of a growing tumor.

The per capita birth rate is $\tilde{g}(w)$; proliferating cells become quiescent at the per capita rate $\tilde{\delta}(w)$, where $w$ (Eq. 2.2) represents the availability of space or some generic nutrient, which we will call growth factor henceforth. We have scaled the maximum cell density to be $1$. In our model, necrosis is not explicitly included but can be regarded as being lumped into $q$. Insofar as quiescent cells cannot become proliferative, $\tilde{\delta}(w)$ can be viewed as a functional death rate. Our motivation in keeping the model framework relatively simple is to be able to estimate model parameters directly from clinical MRI imaging that is sparse in time.

To make the model biologically reasonable, we impose the following constraints on $g(w)$ and $\delta(w)$:

That is, birth (death) should increase (decrease) with the availability of the growth factor; there is more birth than death at maximum values of the growth factor; and there is only death with no growth in the absence of growth factor. It is also assumed that the death rate is negligible at maximum values of the growth factor. With these assumptions, we observe numerically that with suitable initial conditions, the solution of (2.1) stays positive and is bounded ($p+q\le 1$) for all $t$.

We can estimate only up to three parameters based on the necrotic, enhancing, and maximum radii to be measured from MRI images. Therefore, we place a few more restrictions on $\tilde{g}(w)$ and $\tilde{\delta}(w)$ to simplify the estimation of model parameters and to ensure their identifiability. We assume that the proliferation rate at maximum growth factor is $\rho$ and that the death rate at zero growth factor is $k$ and incorporate these parameters into $\tilde{g}$ and $\tilde{\delta}$, respectively; that is, $\tilde{g}(w = 1;\rho) = \rho$ and $\tilde{\delta}(w = 0;k) = k$. For reasons that will become clear later, we pick a functional form that can be written as $\tilde{g}(w; \rho) = \rho g(w)$ and $\tilde{\delta}(w; k) = k\delta(w)$. Some examples include the cumulative distribution function of the beta distribution family (cf. the left pane of Figure 4). These additional assumptions impose little impact on the generality of our model. The benefit of including them will become clear in Section 2.3.

2.2.   Approximate wave profile

In most biological applications of reaction-diffusion models, solutions take the form of traveling waves. MRI images of GBM cancer growth suggest that we can approximate the evolution of the tumor by a traveling-wave solution of its growth model. To uniquely identify and accurately approximate GBM growth model parameters, it is highly desirable to obtain some analytic approximation of the traveling wave, to enable computational matching of the image wave profile and the approximate model wave profile. For this purpose, we consider one spatial dimension, which suffices insofar as the tumor is approximately spherical, and, at the time of diagnosis, its radius is large enough so that radial effects are negligible. With these assumptions, model (2.1) takes the form

We nondimensionlize the system using the characteristic length $\sqrt{D/k}$ and the characteristic time $1/k$ so that $x = \sqrt{D/k}\, \hat{x}$ and $t = \hat{t}/k$, which leads to

where $\hat{\rho} = \rho/k$. We seek a traveling wave solution of the form $p(\xi) = p(\hat{x}-c\hat{t})$, $q(\xi) = q(\hat{x}-c\hat{t})$, where $c$ is the wave speed. Substituting these into (2.5) gives

Linearizing at the wave head, i.e., substituting the ansatz $p = Ae^{-r\xi}$ and $q = Be^{-r\xi}$ into (2.6), gives $(r^{2}-cr+\rho)A = 0$. For a biologically realistic wave front, we expect $A > 0$, $B > 0$, and $r > 0$. This requires that $c^2 > 4\rho$, which implies that the minimum speed of the wave is $c_{\mbox{ min}} = 2\sqrt{\hat{\rho}}$. It is numerically verified that the minimum speed is exactly the asymptotic speed, i.e., $c = c_{\mbox{ min}}$.

To obtain an approximate wave profile, we adopt a method first used by Canosa ^[26]. We rescale the wave coordinate as $z = -\xi/c$, which leads to

Assuming that $1/c^{2}$ is small, we neglect each first term of (2.7). Writing the resulting system in terms of $p$ and $w$, we obtain the reduced system

which is amenable to phase-plane analysis. The approximate wave solution corresponds to a trajectory that leaves $(0, 1)$ and ends at $(0, w^{*})$, with $w^{*}\in[0, 1)$ (see Figure 1). (In the appendix, we show that such a trajectory exists, given the assumptions (2.3).) Dividing (2.8a) by (2.8b) yields

Upon integration, we obtain $p$ as a function of $w$, which we will use in the next section.

2.3.   Parameter estimation

From clinical MRI data, we may derive three idealized radii: $R_{0}$, $R_{1}$, and $R_{2}$, representing respectively the radius of the inner necrotic core, the radius to the edge of the contrast-enhancing rim, and the radius to the outer edge of tumor-associated edema. Such data have been extracted from a series of anonymized patient MRI data consisting of T1-contrast enhanced and T2-weighted MRIs at initial diagnosis. Using the publicly available MATLAB software package, Statistical Parametric Mapping 12 (SPM 12) ^[27], MRIs are initially registered to a standard brain space, and then, using Slicer 3D ^[28] software, the total necrotic core volumes, enhancing rim volumes, and tumor-associated edema volumes are determined from semi-manual tumor segmentation. Finally, these volumes are converted to radii assuming a spherical tumor geometry. The width of the proliferating rim, denoted as $L_{1}$, and the width of the edematous rim, denoted as $L_{2}$, can be calculated as $L_{1} = R_{1}-R_{0}$ and $L_{2} = R_{2}-R_{1}$, as demonstrated visually in Figures 2 and 3.

We assume that contrast-enhancing regions of T1-weighted images correspond to high densities of proliferating tumor cells and that edematous regions on T2-weighted imaging correspond to low densities. We denote the respective detection thresholds for T1 and T2 imaging as $a_1p_{\mbox{ max}}$ and $a_2 p_{\mbox{ max}}$, where $0 < a_2 < a_1 < 1$ and $p_{\mbox{ max}} = \max_{z} p(z)$, i.e., the maximum density of proliferating cells given by the traveling-wave solution.

Often only a single MRI series is available before surgery, although in some cases, a diagnostic MRI followed some days or weeks later by a pre-surgery MRI may be available. In the latter case, the image-derived wave velocity $V$ is the change in tumor radius divided by the length of the time interval.

From our approximate wave profile, we can compute the corresponding quantities to match with MR images (cf. Figure 3). The wave-solution based approximation of the width of the proliferating rim, denoted $\ell_1$, and the width of the edematous rim, denoted $\ell_2$, (in dimensional form) are computed as

respectively, where $w_{1}^{\pm}$ and $w_{2}$ satisfy, respectively, $p(w_{1}^{\pm}) = a_{1}p_{\mbox{ max}}$ and $p(w_{2}) = a_{2}p_{\mbox{ max}}$. Here $p(w)$ is obtained by integrating Eq. (2.9) (see appendix for details). Additionally, the model-derived wave speed $c = 2\sqrt{\rho D}$ can be matched with the image-derived speed $V$. Thus we have three nonlinear equations

from which we hope to find the parameters $D$, $\rho$, and $k$. Given our assumptions, we can simply take the ratio of (2.10a) and (2.10b), which gives

insofar as the integrals are functions of $\hat{\rho}$. Equation (2.12) can be solved for $\hat{\rho}$ analytically in special cases or numerically in general. The monotonicity of $f(\hat{\rho})$ is important for the identifiability of parameters. Once we find $\hat{\rho}$, i.e., the ratio $\rho/k$, all parameters can be found by back substitution.

The above method requires two MR scans taken at two consecutive times prior to surgery to obtain an image-derived estimate of wave speed. If no second MR series is available, then tumor age may be estimated by the tumor radius divided by the wave speed. However, the estimate depends on which radius ($R_{1}$ or $R_{2}$) is used, because the tumor grows exponentially at first and linearly later on ^[7]. This initial exponential growth stage needs to be taken into account as a correction to the aforementioned tumor age estimation. Suppose that for $0\le t\le t^*$, quiescence is negligible and the proliferating cancer cells grow exponentially from a point source of density $p_0$, and that for $t > t^*$, the tumor grows as a traveling wave with speed $2\sqrt{\rho D}$. By equating the two age estimates, we obtain

where

where $R_1$ and $R_2$ are respectively the T1 and T2 radii at $t = t^*$ (see details of $R_1^*$ and $R_2^*$ in the appendix).

Replacing the last equation in (2.11) with (2.13), we again have three equations. To solve them for the unknown parameters, we first take (2.10a) over (2.10b) as before to obtain (2.12). It can then be solved for the ratio $\hat{\rho} = \rho/k$. Substituting this expression back to either $\ell_1 = L_1$ or $\ell_2 = L_2$ gives the raito $\rho/D$. Finally, expressing $D$ and $k$ in terms of $\rho$, (2.13) can be solved for $\rho$, and $D$ and $k$ follow.

3.   Results

Monotonicity is crucial for parameter identifiability, so we first investigate the monotonicity of $f$ for some specific choices of $g(w)$ and $\delta(w)$. Given the restrictions described in Section 2.1, the cumulative density function (CDF) of the Beta distribution family suits our purposes. Therefore, we let $g(w) = B(w; \alpha_g, \beta_g)$ and $\delta(w) = 1-B(w; \alpha_{\delta}, \beta_{\delta})$, where $B(w; \alpha, \beta)$ is the CDF of the beta distribution with shape parameters $\alpha$ and $\beta$. By varying $\alpha$ and $\beta$, we can get linear, sigmoidal, and concave up/down curves (see the left pane of Figure 4). Our framework is robust to those choices, that is, the monotonicity of $f(\hat\rho$), defined by Eq. (2.12), is preserved (right pane of Figure 4). Sigmoidal-shaped growth and death functions ($g$ and $\delta$, respectively) may provide biologically realistic response functions to limited growth factors (most enzymatic reaction rates have sigmoidal shapes with respect to reactant concentration). Given this family of functions, we now consider the question of estimating patient-specific tumor growth and death rates from MR imaging.

We parameterize our model with patient data in which there is only one MRI scan before surgery. In Table 1, we summarize the image-derived tumor radii and the corresponding parameters estimated by the method introduced in the previous section. The parameters $a_1 = 0.9$ and $a_2 = 0.1$ are adapted from values found in the literature ^[16], while $p_0 = 0.02$ and $t^* = 60$ days are hypothetical values. The parameters vary considerably among individual patients.

We compare our approximate quantities to those obtained from the numerical solution of the model. As shown in Figure 5, the approximated results match well with the numerical results except for some discrepancy for $L_{2}$ when $\hat{\rho}$ is small. This result is not a surprise, because the approximation assumes that $c = 2\sqrt{\hat{\rho}}$ is large. Moreover, the numerical approximation of $L_{2}$ is prone to errors due to the fixed grid size and large rate of change around the threshold of $L_{2}$. Overall, we believe that our approximation is accurate for the parameter ranges estimated from the image data.

4.   Discussion

In this work, we have extended the Fisher-Kolmogorov reaction-diffusion model of GBM growth, Eq. (1.1), to explicitly separate the cancer cell birth and death (or quiescence) processes, which are described in terms of generic functions that depend upon an implicit nutrient or growth factor. We specify the birth and death processes, $g(w)$ and $\delta(w)$, respectively, by the cumulative distribution function of a beta distribution, each uniquely specified by a single parameter, $\rho$ and $k$. Thus, along with the diffusion coefficient, $D$, our model describes cancer growth via three parameters, $D$, $\rho$, and $k$, and yields a tumor morphology (in one dimension) consisting of a necrotic core, a high-density rim, and an outer low density rim, which we may correlate to three radii, $R_0$, $R_1$, and $R_2$, that can be estimated from a single patient MR image.

We have demonstrated that our reaction-diffusion system has a traveling-wave solution, which is common in such systems. Studies on this topic date back to the Fisher's work in the 1930s on the spread of advantageous genes ^[12]. Rigorous proof of the existence of a traveling-wave solution in a reaction-diffusion system often leads to phase-space analysis such as the one on the diffusive Lotka-Volterra equations ^[29]. Although in general a rigorous proof of a traveling-wave solution is a daunting task, our reduced system is amenable to phase-plane analysis, and the orbit that represents the traveling-wave solution can be identified (see appendix).

Via traveling-wave analysis, we have developed a method to estimate $D$, $\rho$, and $k$ from as few as a single magnetic resonance image, based on certain growth assumptions. We have estimated these parameters for six patient test cases, as shown in Table 1. Because of the sparsity of imaging data for a typical patient, parameter identifiability in this case is provided by the monotonicity of the function $f(\hat\rho)$ (as seen in left pane of Figure 4); our approach differs from more common statistical practices ^[30] that are appropriate when more data are available.

Disagammaegating the net cell proliferation into birth and death processes not only aids in relating (simplified) tumor appearance on MRI to the model parameters, but it may provide useful valuable information for personalized treatment design, insofar as chemotherapy and radiotherapy target proliferating cells. Moreover, the structural information is also potentially useful, as drug dosages might be selected to ensure penetration through the width of the proliferating rim. Research along this general line has been conducted using model (1.1) ^[31].

Our analysis uses a more complex description of motility than simple diffusion. The diffusion term in Eq. (2.4) belongs to a more general category called cross diffusion ^[32]. It represents the phenomenon in which the gradient in the concentration of one species causes a flux of another species. The type of cross diffusion considered in this paper has been studied in a more general and theoretical context ^[33]. In a modeling study of avascular tumour growth ^[25], the authors justified the adoption of a proportion-based cross diffusion in a tumor-growth model by recognizing that tumor cell migration is "contact inhibited": the presence of one type of cell halts the movement of the other. This type of cross diffusion can cause the solution to become negative ^[32], but we have not encounted this difficulty so far in our numerical simulations. We conjecture that, for suitable initial conditions, solutions remain positive, which will be a topic for future study. Other types of density-dependent diffusion have been considered in modeling GBM migration ^[34].

Despite the existence of many possible diffusion terms, the exact form of diffusion does not matter in the one-dimensional analysis, because the second derivatives are dropped in model (2.7). However, The diffusion coefficient does play an important role in the linearized wave head, where it affects the wave speed, and in the characteristic length where its square root scales the space. The scale-invariant part of the wave profile is mostly determined by the exact forms of the birth and death functions.

The major contribution of this paper is two novel methods to make patient-specific estimates of a three-parameter model of GBM grwoth from the limited MRI data that is typically available in clinical settings. It is possible to estimate the model parameters from a single pre-surgery image. Improved estimates may be possible when images are acquired at multiple time points.

As cancers progress, the underlying parameters describing their growth are unlikely to be static. Data assimilation refers to basic method of updating parameters as more data is acquired, and there has been research on applying full-fledged data assimilation to cancer modeling ^[35,36]. The methods described in this paper can be incorporated as part of a future data assimilation system.

Acknowledgments

The work is mainly supported by a grant from Arizona Biomedical Research Commission and by funds from the Newsome Chair in Neurosurgery Research held by Mark C. Preul. It is also partially supported by NIH grant 5R01GM131405-02 to Yang Kuang. The authors would like to thank Dr. Leslie Baxter and Dr. Leland Hu at the Barrow Neurological Institute for providing MR images.

Conflict of interest

The author declares there is no conflict of interest.

Appendix

Traveling-wave solutions

In the following, we rigorously establish the existence of traveling-wave solutions in system (2.8). We first show that the solutions of system (2.8) with positive initial values are non-negative and bounded.

Lemma 1. Assume that $g(w) = wG(w)$, where $G(w)$ is a bounded function. Then the solutions of system (2.8) with positive initial values are positive and bounded.

Proof: If $x' = x\, f(t, x)$ and $f$ is a bounded function, then $x(t) = x(t_0)\exp \left(\int^t_{t_0} f(s, x(s))\, ds\right)$, which is positive whenever $x(t_0) > 0$. Since

we see that $(p+w)$ is bounded, which implies that both $p$ and $w$ are also bounded.

For any $w^*\in [0, 1]$, the point $(0, w^*)$ is an equilibrium of (2.8).

Theorem 1. System (2.8) admits positive traveling-wave solutions that correspond to heteroclinic orbits connecting the steady state $(0, 1)$ to another steady state, $(0, w^*)$.

Proof: We have performed a detailed phase-plane analysis of (2.8) to show the existence of a trajectory that starts from $(1, 0)$ and ends at $(w^{*}, 0)$, where $w^*\in [0, 1)$. (See Figure 1.) First we notice that

Since $g(w)$ and $\delta(w)$ are both positive functions and $g'(w) > 0$ and $\delta'(w) < 0$ on $w\in [0, 1]$, it follows that $d^{2}p/dw^{2} < 0$ for all $w\in [0, 1]$.

We show first that the trajectory starting from $(1, 0)$ will never cross the line $p = 1-w$. Because $g(1) = 1$, we have $\delta(1) = 0$, and therefore, $dp/dw\big\vert_{w = 1} = -1$. Since $d^2 p/dw^{2} < 0$, the slope of a trajectory with $w < 1$ will be greater than $-1$, which means it will not cross the line $p = 1-w$. Also, because the solution components stay positive, the trajectory starting from $(1, 0)$ will never cross the $p$-axis from right to left.

Because the functions $g$ and $\delta$ are monotone, there is a unique value $w^{\dagger}\in (0, 1)$ such that $dp/dw\big\vert_{w = w^{\dagger}} = 0$. Moreover, $dp/dw\big\vert_{1 > w > w^{\dagger}} > 0$ while $dp/dw\big\vert_{0 < w < w^{\dagger}} < 0$.

Insofar as $w$ is strictly decreasing and bounded from below by 0, there exists some $w^*\in(0, 1)$ such that $\lim_{z \rightarrow \infty} w(z) = w^*$. We claim that $w^* < w^{\dagger}$. Otherwise, $p$ is a non-decreasing function, which implies that the trajectory approaches a positive steady state $E^* = (p^*, w^*)$. However, the system (2.8) does not admit any positive steady states.

Let $a \in (w^*, w^{\dagger})$ and $b = \hat{\rho}g(a)-\delta(a) < 0$. Since $\lim_{z \rightarrow \infty}w(z) = w^*$, there is a $z^* > 0$ such that for $z \geq z^*$ and $w < a$. Therefore, for $z \geq z^*$, $dp/dz < bp$, which implies that $\lim_{z \rightarrow \infty} p(z) = 0$. Hence, system (2.8) admits positive traveling-wave solutions that correspond to heteroclinic orbits connecting the steady state $(0, 1)$ to another steady state, $(0, w^*)$.

Rim width

The trajectory in Figure 1 corresponds to a traveling-wave profile in the $z$ coordinate as shown on the left pane of Figure 3. Because of our choice of nondimensionlization, $x = -\left(2 \sqrt{D\rho}\, /k\right) z$, by definition we have the rim width in the dimensional form

where $p(z^{\pm}_1) = a_1 p_{\max}$. We have taken $p$ as a function of $z$ (cf. the left pane of Figure 3). We can also take $w$ as a function of $z$ to make a change of variable to the above integral, which gives Eq. (2.16). A similar argument applies to (2.17).

Derivation of $R_1^*$ and $R_2^*$

Tumor growth can be separated into two stages. First, the tumor cells grow exponentially until the cell density is high enough ($p_{\mbox{ max}}$) to form a stable wave profile. Afterward, the growth of tumor cells is described by (2.8). The radii $R_1$ and $R_2$ are observed from clinical MRI data and correspond respectively to the distance from the center of the tumor to the edge of the enhancing rim and to the edge of the edematous region on T2 imaging; they correspond to tumor dynamics before the stable wave profile has formed (cf. Figure 6).

During the exponential growth phase, say from $0 < t < t^*$, quiescence is negligible and the governing equation of tumor cell density is

where spherical symmetry is used, as we assume that the tumor is spherical when its radius is small. This linear equation has the Green's function

(see, e.g., ^[37] Page 314 and ^[38] Page 285). Suppose that the tumor starts at $t = 0$ as a point source with density $p_0$. It follows that at $t = t^*$, the position of the tumor front at density $p = a_i$ is given by Eq. (2.14).