Stochastic growth pattern of untreated human glioblastomas predicts the survival time for patients

Glioblastomas are highly malignant brain tumors. Knowledge of growth rates and growth patterns is useful for understanding tumor biology and planning treatment logistics. Based on untreated human glioblastoma data collected in Trondheim, Norway, we first fit the average growth to a Gompertz curve, then find a best fitted white noise term for the growth rate variance. Combining these two fits, we obtain a new type of Gompertz diffusion dynamics, which is a stochastic differential equation (SDE). Newly collected untreated human glioblastoma data in Seattle, US, re-verify our model. Instead of growth curves predicted by deterministic models, our SDE model predicts a band with a center curve as the tumor size average and its width as the tumor size variance over time. Given the glioblastoma size in a patient, our model can predict the patient survival time with a prescribed probability. The survival time is approximately a normal random variable with simple formulas for its mean and variance in terms of tumor sizes. Our model can be applied to studies of tumor treatments. As a demonstration, we numerically investigate different protocols of surgical resection using our model and provide possible theoretical strategies.

has been applied to several types of solid tumors 7 . However, the Gompertz growth model often exhibits discrepancies between clinical or experimental data and theoretical predictions 8 . These discrepancies may be due to internal environmental fluctuations and variation among patients. To consider such environmental fluctuations and variation in individual patients, some tumor growth models incorporate stochastic processes. For example, Albano & Giorno 9 and Lo 10 considered some stochastic models for tumor growth.
In a recent clinical study, Stensjøen et al. reported a data set of untreated human glioblastoma in vivo 3 . They calculated the tumor growth rates for each patient, and found considerable variation among individual patients.
In this study, we propose a Gompertz diffusion growth model that is the best fit to the data set of 94 untreated glioblastoma patients collected in the aforementioned study conducted by Stensjøen et al. 3 in Trondheim, Norway. Our model is a stochastic differential equation, where white noise captures the variants of glioblastoma growth among patients. Our model is re-verified by newly collected data sets from Seattle, US. By computer simulations, we find that the patient survival time is a normally distributed random variable and its mean and variance are explicit functions of current tumor volumes. Using our model, we can predict the survival time for each patient based on a magnetic resonance imaging (MRI) scan with a prescribed probability before treatments.
As proof of principle, we conduct a theoretical study of surgical resection using our mathematical model. The extent of resection [11][12][13][14][15] and repeated resections 16,17 have been important subjects in neurosurgery clinics and research. We will study a variety of resection protocols including extent and repeated resections. In particular, we will study the following questions. (1) How long could a patient survive without treatment given the initial diagnosed tumor size? (2) How long could a patient survive if the patient accepts surgical resection immediately after diagnosis of glioblastoma? (3) How long could a patient survive if the patient accepts surgical resection after some period of time after diagnosis of glioblastoma? (4) If a second surgical resection is performed, how long could a patient survive? Our computational study will answer these questions in detail and provide some suggestions which may help surgeons to determine optimal treatment strategies. Computationally, we confirm three established conclusions: patients with small glioblastomas have a longer survival time while large glioblastomas have a short survival time 11,12 ; the patient will have survival benefits even for 50% of the extent of resection 13,14 ; repeated resections will increase patient survival time 16,17 . Quantifying the extent of resection, we obtain the following theoretical conclusions. For large glioblastomas, the more a surgical resection can cut off from the tumor, the longer the survival time for the patient. For small glioblastomas, if we wait until the glioblastoma grows to a large tumor and then perform surgical resection, the patient survival time will increase; furthermore, if we wait again until the glioblastoma grows to an even larger tumor and then perform a second surgical resection, the patient survival time will increase further. However, these conclusions are purely theoretical outcomes which do not consider many medical factors, and therefore may not be immediately applicable in the clinics as many other aspects of the tumor need to be taken into account.

Methods
In the study conducted by Stensjøen et al. 3 , 106 untreated human glioblastoma data sets were collected. In that study, there are tumor volumes from two MRI scans for each patient, and the time between the two MRIs. These measurements were taken before surgery. Specifically, for the i-th patient, we have the record ∆ ( V , t , V )  Fig. 1, there are two key observations on the growth dynamics: Figure 1. Scatter plot of SGR against tumor sizes of the first MRI scan: The horizontal axis presents the tumor size in the first MRI scan and the vertical axis presents the specific growth rate (SGR). The plot shows there is a negative correlation between the first scanned tumor volume and SGR; the variation of SGR decreases with respect to tumor size before the tumor volume reaches 100 mL.

Model set up. From
1. There is a significant negative correlation between the first scanned tumor volume and SGR. Stensjøen et al. also pointed this out. They carried out a statistical analysis on the distribution of SGR with different tumor sizes 3 . 2. The variation of SGR decreases with respect to tumor size before the tumor size reaches 100 mL.
Our model will incorporate these two characteristics. Let = V V(t) be the tumor size at time t. In most cases, phenomenological tumor growth laws are described by a single deterministic differential equation describes the tumor growth dynamics and v 0 is the initial tumor volume. Typically, the model (1) has the following simplifying assumptions 9,10 : 1. The tumor has one cell population, and there is no cell heterogeneity; 2. Tumor growth is spatially independent; 3. Tumor growth does not show explicit dependence on nutrients, host vasculature, or age.
The function g V ( ) can be expressed as the product of V and the intrinsic net growth rate r V ( ), thus It is well known that glioblastomas are composed of many cell types and are heterogeneous 19 . The heterogeneity in cell types and space may contribute to variation of intrinsic growth rate. We employ a stochastic process to express the tumor intrinsic growth rate. Let represents the deterministic intrinsic growth rate, which will be a decreasing function since growth rate decreases with tumor size, ξ t ( ) represents white noise, and σ V ( ) represents the strength of the white noise, and is a decreasing positive function. To parameterize the model, we need to undertake two steps of data fitting to determine α V ( ) and σ V ( ) from our data, separately.
Model fitting. In the literature, many studies report that Gompertz curve fits for tumor growth are better than others [20][21][22] . Therefore, we use a Gompertz curve to fit the deterministic intrinsic growth rate, i.e., , 0 are parameters. More concretely, a represents the intrinsic growth rate and b represents the carrying capacity. The individual conditions including nutrients, host vasculature, and age may contribute to the carrying capacity. Using pairs V SGR ( , ) a 0 009916 and = . b 121 6. The fitted curve for α V ( ) is shown in Fig. 2. Secondly, we fit the strength of the white noise σ V ( ). Since there is no commonly accepted functional form for that, we will find a function that is simple and fits our data well. To fit σ V ( ), we need pairs σ V ( , ) i i . For the sake of simplicity, we partition the tumor size of the first MRI scans, V i 1 , into 7 subgroups, D j , for =  j 1, , 7 based on the tumor size. To be concrete, let 10, 20, 30, 40, 50, 60 be volume separators for tumor size at the first MRI scan, , where V j 1 is the average value of D j , and D j is the total number of elements in D j . We dropped the group D 7 since it only has 7 data points in the range from 61.7 to 146.5. Thus, we have six ordered pairs of the form σ V ( , ). Using these six ordered pairs, we find c 0 0769 and = .
h 0 2241. The fit is shown in Fig. 3. , is used to fit the deterministic intrinsic growth rate, where a represents the intrinsic growth rate and b represents the carrying capacity. The best fitted values for = .
c 0 0769 and = . h 0 2241. Rewriting (3) in the usual form of a stochastic differential equation (SDE), we have where W t ( ) represents the standard Wiener process 24 . This is a stochastic differential equation, and is our mathematical model for glioblastoma growth.
Model verification with the original 94 data sets from trondheim, norway. Theoretically, for a given SDE and its initial condition, we can provide a confidence region of a given probability for the tumor size at any specific future time t. For our Gompertz type diffusion model (4), since there is no closed form solution, we do not have the exact transition probability distribution for predicting tumor sizes. However, we can numerically predict the tumor size with bounds given by a prescribed probability. For example, we can predict the tumor size at future time t after the first MRI scan with upper bound and lower bound given by a probability. Here we resort to numerical methods to verify our model as follows. Based on the data of the i-th patient, ∆ ( V , t , V ) i as an initial condition for the SDE (4) and then we simulate 500 solution sample paths on the time interval ∆t ( 0, ) i . From these 500 solution sample paths, we obtain an empirical distribution at time = ∆ t t i which provides the mean tumor size, µ i , and standard deviation, , 94, which means we give a 99% confidence interval. The results are shown in Fig. 4 and Table 1. Figure 4 shows our model  www.nature.com/scientificreports www.nature.com/scientificreports/ predictions of the tumor sizes at the time the second MRI scan was done against the real tumor sizes. Table 1 shows how the real tumor sizes will fall in the region of our model prediction if we give a 99% confidence interval. As the table shows, about 60+ percent of real tumor sizes fall into our prediction intervals.
In addition, Table 1 summarizes the results of predictions of the tumor size of the second MRI scan to a 99% confidence level.
Model re-verification with 3 data sets from Seattle, US. Our research group recently collected 3 untreated human glioblastoma data sets of tumor volumes at two different time points. The tumor volume was computed from three measurements in three directions that are perpendicular to each other as an ellipsoid. If we use our mathematical model with exact parameter values estimated from the original 94 data sets, two of three tumor volumes are in our model prediction range with probability 0.99 as shown in Table 2.
If we change the growth rate from the estimated value = .
a 0 009916 to = . a 0 015, all three tumor volumes are in our prediction range with 0.99 confidence level as shown in Table 3.
We do not have age records attached to the data sets. We also notice that the methods of measurements are different. The contrast-enhanced T1-weighted MRI scan was used in the data sets collected in Trondheim, while three-dimensional measurements with ellipsoid computation was used in the data sets collected in Seattle. Those methods may contribute some systematical differences in data sets 25 . In either case above, we consider our mathematical model is re-verified by these newly collected data sets.
Ethical approval. All procedures performed in studies involving human participants were in accordance with the ethical standards of the institutional and/or national research committee and with the 1964 Helsinki declaration and its later amendments or comparable ethical standards. They were approved by the Institutional Review Board of Fred Hutchinson Cancer Research Center. The data privacy regulations were followed.
Informed consent. Informed consent was obtained from all individual participants included in the study.

Results
Empirical distribution of survival times obtained by using numerical simulations. Since SDE (4) can predict the tumor size at future time t after the first MRI scan, we can obtain more insight on the glioblastoma growth dynamics based on this model. It is natural to explore the survival time for each patient. If we make the simplifying assumption that tumor size is the only reason for patient death, then the survival time for each patient can be expressed as the first time the tumor size reaches a specific value. We can use the first passage time random 0} is a diffusion process, we define the first passage time (waiting time) random variable, , v 0 is the tumor size for which we start our observations, v c is the tumor size for which we stop our observations. Since there is no closed form solution to SDE (4), we study the first passage time random variable numerically. For initial tumor size = are close to normal distributions. So, we assume that the distribution of the first passage time follows , we simulated 500 sample paths, thus we obtain 39 triples . Table 4 shows triple of    Table 4, Triples of  Table 4. Figure 7 shows the fitted surface σ v v ( , ) c 0 for the standard deviation of the waiting time, which depends on initial tumor size, v o , and tumor stopping size (stopping observation), v c , based on the simulated data in Table 4. Application of our model: theoretical surgery strategies. The extent of resection has been an important subject in neurosurgery research and practice. It has only recently been established that there is a statistical correlation between the extent of resection and survival [11][12][13][14][15]26 . However, almost all of the studies are retrospective and thus subject to numerous sources of bias and variation. Furthermore, there are only three categories for the extent of resection: stereotactic biopsy, subtotal resection, and gross total resection (radical or complete resection). In clinics, tumors in some locations may not be amenable to gross total resection due to proximity to  functional tissue. For example, the relation of the tumor to eloquent cortex (i.e. primary motor cortex or speech centers) is important since it directly corresponds to the proportion of the tumor that can be resected without permanent morbidity 11 . For elderly patients, the optimal extent of resection remains uncertain because of potential higher rates of mortality and morbidity associated with more extensive degrees of resection 12 . The choice of optimal extent of resection may depend on the tumor size and location, the patient's general and neurological status, and the neurosurgeon's experience 13 . Although there are some prospective studies 15,27,28 , a prospective study with quantifying extent of resection is necessary. Our SDE model is a good tool to study extent of resection prospectively. We can explore the survival time after surgical resection and may help surgeons to choose a best protocol for the treatment based on an individual patient's conditions, such as tumor size and location and possible resection proportion (quantified extent of resection). In the following, we will explore how various combinations of different diagnosed tumor sizes, tumor sizes when surgery is performed, different resection proportion, and repeated resections, will affect the patient survival time.
We assume that the glioblastoma growth pattern after resection is the same as before the resection. Specifically, the waiting time =  . In Table 5, we present the expectation and standard deviation (SD) of the survival time for a patient without treatment with an initial diagnosed tumor size v 0 , where we take the initial diagnosed tumor size from 1 mL to 50 mL for demonstration, and the unit of the survival time is days. It is easy and obvious to conclude that small glioblastomas have longer survival times.
The second question. For the second question, we use the random variable and c% is cut-off percentage. In Table 6, we present our computation results for the , T for a patient who accepts the surgical resection immediately after tumor diagnosis and the resection proportion is c%. We only present the cases where the diagnosed tumor size has volume from 1 mL to 50 mL and the resection proportions are percentages from 50 to 98. From our computation, we learn that a surgical resection does not increase survival times much for small tumors no matter what percentage of the tumor can be cut off, while for large tumors the more a surgical resection can cut off from the tumor, the longer will be the survival time for the patient. In the clinical study 14 , the authors reported that a high extent of resection improves overall survival, and even for 78% of resection extent, it also gives survival benefits, which agrees with our conclusions.

The third question.
For the third question, the tumor grows for some period of time from the initial diagnosed size v 0 to size v c and then a surgical resection is performed. The survival time is . Since the survival time is a sum of two random variables, we need to consider how these two random variables affect each other and their impact on the survival time. To include the relation between these two random variables, we use the correlation coefficient ρ. From the property that the sum of normally distributed random variables is still a normally distributed random variable, we know that In this computation, there are four parameters: the initial diagnosed tumor size v 0 , the tumor size v c when surgical resection is performed (after random time T 1 ), the resection proportion c%, and the correlation coefficient ρ. To demonstrate, we choose the initial diagnosed tumor size to be = v 1 mL 0 , then let the tumor grow to size v c for surgical resection with cut-off percentage c%. To compute the survival time, we choose three values for the correlation coefficient ρ for each theoretical treatment. Table 7 shows the expectation and standard deviation of the survival time for the parameter values we choose.
From our computation, we can conclude that for a small tumor, if we wait until the tumor grew to a big tumor and then perform surgical resection, the survival time would increase. If two survival times are positive linear related, the variance of the survival time will increase; if they are negative linear related, then the variance of the survival time will decrease. It may be reasonable to assume the survival times are positive linear related. However, the mean survival time keeps the same.
The fourth question. We consider a second surgical resection. Based on the procedure in the third question, after the first surgical resection for some period of time the tumor will grow to size v c 2 , then a second surgical resection is performed. After the second resection, if the random time for the tumor growth to death size is T 3 , the survival time for the patient is as follows  This is a normally distributed random variable with distribution µ σ ″ ′′ N( , ), w here .
In this computation, there are eight parameters: the initial diagnosed tumor size v 0 , the tumor size v c 1 when the first surgical resection is performed (after random time T 1 ), the first resection proportion c % 1 , the tumor size v c 2 when the second surgical resection is performed (after random time T 2 ), the second resection proportion c % 2 , and three correlation coefficients ρ ρ ρ , , 12 13 23 within T 1 , T 2 and T 3 . To demonstrate and simplify presentation, we choose the initial diagnosed tumor size v 0 =1 mL, the first and second resection proportion = c % 80% 1 and = c % 50% 2 , and give the triple ρ ρ ρ ( , , ) 12 13 23 three choices. For the tumor sizes v c 1 and v c 2 , we choose some volume from 2 mL to 50 mL and 15 ml to 90 mL respectively. Table 8 shows the expectation and standard deviation of the survival time for these parameter values.
From our computation, we may conclude that for a small tumor, if we wait until it grows to a large tumor and then perform a surgical resection with a high cut-off percentage, and wait again until the tumor grows to an even larger tumor then perform a second surgical resection, the patient survival time will significantly increase. We also conclude that correlation coefficients within waiting times for tumor growth do not have effect much on the survival time. This confirms the results reported in the studies 16,17 .
Our computational investigation first confirms the intuitively expected result that small glioblastomas have a longer survival time while large glioblastomas have a short survival time. We also confirm that the patient will have survival benefits even if the extent of resection is 70%. Our computation thirdly confirms that a second resection will increase patient survival time. In addition, we obtain the following conclusions. For large glioblastomas, the more a surgical resection can cut off from the tumor, the longer will be the patient survival time. For small glioblastomas, if we wait until the glioblastoma grows to a large tumor and then perform surgical resection, the patient survival time will increase; furthermore, if we wait again until the glioblastoma grows to an even larger tumor and then perform a second surgical resection, the patient survival time will increase further.  Assume tumor initial size is v 0 . After random time T 1 , the tumor reaches size v c1 , then resection with c 1 %; then after another random time T 2 , the tumor reaches size v c2 , the second resection with c 2 %; then the patient survives with some period of time T 3 . To demonstration, taking v 0 = 1 mL, c 1 % = 80%, c 2 % = 50%, and correlation coefficients ρ 12 , ρ 13 , ρ 23 among T 1 , T 2 and T 3 . For example, the total survival time will be 492 days if the first resection with 80% is perform ed when the tumor reaches 50 mL and the second resection with 50% is performed when the tumor reaches 90 mL.
It is clear that our theoretical study provides ample sources of thought for the surgical management of glioblastomas, and while our predictions provide results that are intuitively obvious, they quantify the possible benefits of resections. There are many factors in clinical practice which should be taken into consideration when we make decisions. For example, waiting for a glioblastoma in the brain to get bigger is not only dangerous in terms of letting cells penetrate deep into the tissue but also dangerous in terms of acute symptoms due to increased brain pressure. This requires us to consider parameters such as metastasis, extent of spatial heterogeneity, and invasiveness. Most importantly, the genetics of the tumors should be considered 29,30 .

Discussion
The clinical data on the growth rate of untreated glioblastomas show a large variability among different patients 3 . Several deterministic mathematical models have been proposed to describe glioma growth patterns. However, these mathematical models are unable to describe and predict the variance among patients (or experimental animals). Based on a data set of 94 glioblastoma patients, we find the best fitted mathematical model based on Gompertz growth with white noise. We use white noise to model the variability among patients and random factors in tumor cell environment. The strength coefficient of the white noise we fitted is of fraction form. Instead of growth curves predicted by deterministic models, we predict a band where the center curve is the mean of the solution (the tumor size), and its width is given by a prescribed probability that describes the variance of the tumor size deviating from the average over time. As mentioned in the Materials and Methods section, almost 2/3 of glioblastoma growth in the data set from Trondheim, Norway, falls into our predicted band. Newly collected human glioblastoma data from Seattle, US, also fall into our model prediction band after the growth rate is increased. Although our stochastic differential equation model captures some variability of the growth dynamics of glioblastomas, we cannot expect a model as simple as the one we propose to hold for all human glioblastomas, which have great heterogeneity among individual patients.
Since there is no closed form solution of our stochastic differential equation model, we numerically simulate our model for a large number of solution paths, and obtain some empirical results. One result is the empirical distribution of the survival time. We find it to be a normally distributed random variable, with mean and variance dependent on the tumor sizes. From our computational simulations, we numerically fit simple formulas to calculate the mean and variance in terms of tumor sizes. If the tumor size of a glioblastoma patient is found from an MRI scan, our empirical distribution of the survival time can predict how long this patient can live without treatment under a given probability, and can also calculate their average survival time and variance. This may give physicians and patients some information to prepare for a better life in recovery.
One interesting question is: what is the best time for surgical resection when a patient is diagnosed with a glioblastoma? This problem involves several parameters including tumor size, tumor location in the brain, physician's experience, possible extent of surgical resection estimated, etc. For example, for a patient with a glioblastoma, each physician has their own judgment of how much they can cut the tumor based on their experience. As a theoretical study, we consider the percentage of the tumor to be removed which is quantified extent of resection, the timing of surgical resections, and times of surgical resections. We confirm an obvious conclusion that small glioblastomas have a longer survival time while large glioblastomas have a short survival time. We also confirm that the patient will have survival benefits even if the extent of resection is low. Our computation shows that for large glioblastomas, the higher the percentage a surgical resection can cut off from the tumor, the longer will be the patient survival time. However, for small glioblastomas, a surgical resection does not increase survival time much no matter what percentage of the tumor can be cut off. If we wait until a small glioblastoma grows to a large tumor and then perform surgical resection, the patient can attain a longer survival time; furthermore, if we wait again until the glioblastoma grows to an even larger tumor and then perform a second surgical resection, the patient will get a significant longer survival time. We may conclude that repeated surgical resections will give more benefits to patients theoretically. In fact, the reoperation has been done in some clinics and experiments 16,17 . A recent systematic review and meta-analysis of contemporary literature suggests that repeated surgery at glioblastoma recurrence in select patients confers a significant, prognostic overall survival advantage independent of other prognostic factors 31 . These newer studies are significantly more likely to suggest greater benefit than are older studies, while larger prospective randomized controlled studies are needed to validate these findings 31 .
It is clear that our results are based on our simple mathematical model and its assumptions. These assumptions do not take all necessary factors into account and so our model predictions must be treated with great caution. For instance, we assume the tumor size determines patient survival. However, in many cases, the real driver of survival is the genetics of the tumors 29,30 which may outweigh the tumor size at diagnosis. In the theoretical resection, we also assume that resection does not alter growth dynamics. This may not be true in some cases since the neurological conditions may change after resection. The result that the patient can attain a longer survival time when a small glioblastoma grows to a large tumor and then a surgical resection is performed, may be explained by observing that the time for a small tumor to grow to a large tumor added to the time that the patient survives after the resection is longer than the time taken for a small tumor to grow to a small or medium tumor plus the time that the patient survives after the resection. However, glioblastomas are known to be spatially diffuse, so that by allowing a tumor to grow, we risk potentially fatal spread to other parts of the brain. It is clear that, in neurosurgery, the aim of resection is always total resection, but without impacting on patient's neurological function, and to minimize the risk of new postoperative deficits such as brain mapping/stimulation, navigation with diffusion tensor imaging. The next step in the type of analysis presented here would be to link the tumor growth rate with the genetic information and to develop a spatially extended model or adapt some of the models in the literature [32][33][34][35] so that one could apply some of the ideas from the present paper. However, that would involve a significant extension and is beyond the scope of the present paper.
In addition, we did not consider other treatments which can be easily incorporated into our model. Surgery is usually combined with other treatments, such as radiotherapy and chemotherapy. When relevant clinical or Scientific RepoRtS | (2020) 10:6642 | https://doi.org/10.1038/s41598-020-63394-w www.nature.com/scientificreports www.nature.com/scientificreports/ experimental data are available, we should be able to incorporate other treatments into our stochastic differential equation model or into a spatially extended stochastic differential equation model, and provide more insight into glioblastoma growth and treatments.

conclusion
Human untreated glioblastoma growth in vivo may follow a new type of Gompertz diffusion dynamics, which is described by a stochastic differential equation. Our mathematical model not only captures the average growth pattern and the variance among individual patients, but also can predict, in many cases, individual glioblastoma growth paths. The empirical distribution of survival times simulated from our mathematical model can be used to calculate a patient's survival time with a prescribed probability. We obtain empirical formulas to easily calculate the average survival time and its variance. As proof of principle, our mathematical model can be used to provide different protocols for performing surgical resections according to tumor sizes which will give patients long survival times. The conclusion should be interpreted with caution, owing to the number of simplifications of the mathematical modeling and small size of the data set.