Evaluation of Classical Mathematical Models of Tumor Growth Using an On-Lattice Agent-Based Monte Carlo Model

Purpose: To analyze the capabilities of different classical mathematical models to describe the growth of multicellular spheroids simulated with an on-lattice agent-based Monte Carlo model that has already been validated. Methods: The exponential, Gompertz, logistic, potential, and Bertalanffy models have been fitted in different situations to volume data generated with a Monte Carlo agent-based model that simulates the spheroid growth. Two samples of pseudo-data, obtained by assuming different variability in the simulation parameters, were considered. The mathematical models were fitted to the whole growth curves and also to parts of them, thus permitting to analyze the predictive power (both prospective and retrospective) of the models. Results: The consideration of the data obtained with a larger variability of the simulation parameters increases the width of the χ2 distributions obtained in the fits. The Gompertz model provided the best fits to the whole growth curves, yielding an average value of the χ2 per degree of freedom of 3.2, an order of magnitude smaller than those found for the other models. Gompertz and Bertalanffy models gave a similar retrospective prediction capability. In what refers to prospective prediction power, the Gompertz model showed by far the best performance. Conclusions: The classical mathematical models that have been analyzed show poor prediction capabilities to reproduce the MTS growth data not used to fit them. Within these poor results, the Gompertz model proves to be the one that better describes the growth data simulated. The simulation of the growth of tumors or multicellular spheroids permits to have follow-up periods longer than in the usual experimental studies and with a much larger number of samples: this has permitted performing the type of analysis presented here.


Introduction
A multicellular tumor spheroid (MTS) is a spheroid-shaped cellular aggregate that grows in suspension under controlled culture conditions [1,2]. Though actual tumors are more complex than any in vitro biological model, MTS can be considered as an in vitro 3D cell structure with a complexity similar in some aspects to that of in vivo tumors or micro-metastasis with similar sizes [3][4][5][6][7]. As a consequence, MTS are useful for tumor biology studies, as well as for biomedical and clinical applications [8][9][10][11]. In this regard, it is worth noting that the growth kinetics of MTSs is similar to that of early in vivo avascular tumors, showing the physio-pathology of tumor micro-regions and metastases, and other characteristics such as metabolic alteration, altered gene and protein expression, necrotic nucleus, resistance to therapy, and appearance of hypoxic regions at similar distances from nutrients [1,10,12]. In MTSs, nutrients and oxygen penetrate their structure only by diffusion and this means that at distances between 50 µm and 250 µm from the surface, the cells are in a hypoxic, or even necrotic, state due to nutrient and/or oxygen deficiency [13].
MTSs also allow the study of how their proliferation kinetics is affected by treatment strategies such as chemotherapy, hormone therapy and/or radiotherapy [14][15][16][17][18][19][20]. More specific phenomena such as low-dose hypersensitivity [21] or the effects of different radiotherapy fractionation strategies [22] have also been tested in this experimental model.
Despite their complex nature, the kinetics of tumor and MTS growth seems to follow relatively simple laws that can be expressed through a great variety of mathematical models. The exponential, the Gompertz, the logistic, the Bertalanffy and the potential models are among the most widely used [23,24].
All these mathematical models include several free parameters that are chosen by fitting the corresponding functions to in vivo tumor [24,25] or in vitro MTS [23,26,27] data. The fitted models have been applied to determine the dosage design and scheduling of cancer drugs [28][29][30], or to predict the cytotoxic or anti-angiogenic effects of drugs on tumor growth [31,32]. Another common application in clinical practice is the evaluation of the modifications that need to be included in the treatment plan when delays or interruptions occur in radiotherapy. As far as we know, only the exponential model has been (and is still being) used for this purpose and we are not aware of any other of the aforementioned models being considered for this task [33].
One of the characteristics of these mathematical models, which is crucial for the aforementioned applications, is their predictive capability; that is, how they describe the tumor or MTS growth at times other than those used in the fitting procedure; however, this characteristic has not been studied in depth, may be due to the fact that, usually, only a few growth curves, in which the amount of volume data available is typically very small, are analyzed. For example, Benzekry et al. [24] studied the evolution of lung and breast cancers xenografted in in vivo animal models, testing the capabilities of various growth mathematical models. Their data series showed a follow-up of ∼15 days with ∼10 observations; they fitted the model using the first 3-9 pieces of data and compared their prediction with the experimental observation found between 1 and 7 days later. As another example, we can mention the work of Murphy et al. [34] who used a single growth curve including 14 data points and spanning 114 days. The situation is similar in many other works (see, e.g., Refs. [35][36][37][38]. Though the information obtained in this way is useful, the predictive capability of the models cannot be tested in depth due to the limited amount of data available for each single growth curve and the small number of such curves in the sample. Tumor growth can also be simulated by means of computational models of different types [39,40] that require fulfilling a number of conditions (see, e.g., Ref. [41] and references therein). These models have been used to study the effects of radiation on in vitro and in vivo tumor cells, predicting the efficacy of radiotherapy treatments [42][43][44][45], to study certain chemotherapy or surgical treatments for cancer in pre-clinical, clinical, biological, and molecular studies [46][47][48][49], or to investigate how to personalize cancer treatments [50].
Among the different computational approaches considered to describe tumor growth, it is worth pointing out the "agent-based models" (ABM) [51,52]. These models deal with a discrete population of cells whose states are characterized by vector variables that involve the position of the cell, its current cycle phase, details about how it interacts with the environment, etc. In these models, the agents can be genes, proteins, metabolites, or the cells themselves, and they interact with each other according to predefined rules that depend on their spatial situation and the instant of the system evolution that is being analyzed. Both space and time are considered as discrete variables, no formal spatial grids or time synchronies are required, at least a priori, and all the processes involved may include randomness.
Recently, an on-lattice ABM of MTS growth that uses simple mathematical tools to describe the cell evolution has been developed [53]. The model includes Monte Carlo techniques to take into account the variability observed among the cells forming the spheroids and takes into account the basic features that characterize the MTS behavior when irradiated, allowing for addressing the capabilities of different therapeutic strategies for solid tumors. This model has been developed and validated in such a way that the simulated spheroids show the same growth behavior as the experimental ones, reproducing characteristics such as the shapes of the different parts of the spheroids, the relation between the radius of the necrotic core and that of the whole MTS, or dependence of the number of proliferative plus hypoxic cells with the MTS volume.
An advantage of almost all computational models is their ability to adequately describe the tumor evolution in time intervals longer than those usually observed in experiments carried out with in vivo or in vitro tumor models. This is precisely the characteristic that we want to exploit in this work to determine the predictive capabilities of various classical tumor growth models. By fitting the models to a part of the growth curves generated by our ABM, the ability of these classical mathematical models to reproduce the past and future growth of the virtual tumors can be analyzed. The aim of this work is to perform a quantitative evaluation of the predictive power of the exponential, Gompertz, logistic, potential, and Bertalanffy models.
The results obtained in this analysis may help to establish the potential prognostic capabilities of the simple mathematical models in the clinic, addressing their use in applications such as those mentioned above.

Classical Mathematical Models of Tumor Growth
As said above, the purpose of this work is to analyze the predictive capabilities of various classical mathematical models that have been employed to describe MTS (and tumor) growth. In all the classical mathematical models of the tumor growth kinetics considered in the present work, the variable described is the total volume of the tumor, V, as a function of time, t. Furthermore, this total volume is supposed to be proportional to the total number of cells in the tumor. To simplify the labeling, E, G, L, P, and B have been used to refer to the exponential, Gompertz, logistic, potential, and Bertalanffy models that are described in what follows.
The exponential model assumes that all cells within the MTS have a cell cycle with a constant duration [54]. The same growth behavior occurs if it is supposed that the cell proliferation is due to a constant fraction of the MTS volume [25] or if the duration of the cell cycle is a random variable following an exponential distribution. The mathematical expression describing this model is: where V 0 , the volume at t = 0, and E are the model parameters.
The observation of non-constant doubling times in some tumors led to considering more elaborated models. One of the most widespread and accepted is the Gompertz model, which has been used to describe population growth in many different branches of knowledge [25,35,37]. The Gompertz model shows a sigmoid shape, i.e., a rising curve with an inflection point that converges asymptotically to a maximum volume. This qualitatively reproduces the growth slowdown observed experimentally [37,55] and is consistent with general growth patterns of organs and organisms. The essential characteristic of the Gompertz model is that it shows a relative growth rate that reduces exponentially when volume increases. One of the main criticisms of the Gompertz model is that the growth relative rate may become arbitrarily large (or, equivalently, that the volume doubling time may reach arbitrarily small values) for small tumor volumes. In the present work, the following expression has been used for this model: In this model, V 0 , which gives the MTS volume at t = 0, and A and a are the model parameters. The logistic model also shows a sigmoid shape [23,25]. It is based on the equilibrium of metabolic processes and is defined by assuming that the relative growth rate reduces linearly with the volume. In this work, the following analytical expression has been considered: In this case, the model parameters are V 0 ≡ V(t = 0), B, and b. Von Bertalanffy [56] and others [57] proposed deriving the general laws governing the organic growth from basic principles involving the energy. They established that the net growth rate should be the result of the balance of synthesis and destruction, noting that metabolic rates often follow the law of alometry (i.e., they scale with a power of the total size of the system analyzed), while catabolic rates are proportional to that total volume. The Bertalanffy model has already been successfully applied to describe tumor growth [58,59]. The mathematical function corresponding to this model is: Here, the model parameters are the volume at t = 0, V 0 , C, and c. The potential model was derived from the Bertalanffy model [56] and does not exhibit a clear saturation phase. Tumor growth is proportional to the number of proliferating cells, under the assumption of a constant cell cycle duration. It is a model that produces a tumor growth with a decreasing growth fraction and, therefore, a decreasing relative growth rate, which also provides a description in terms of a geometric characteristic of the proliferative tissue. The time dependence of the volume in this model is given by: The model parameters are V 0 , the volume at t = 0, P, and p.

On-Lattice Monte Carlo Agent-Based Model
In order to investigate their predictive capabilities, the classical mathematical models described in the previous section have been fitted to data obtained with an on-lattice ABM recently developed and tested [53], which has been used as a reference. In this ABM, cells are situated at the vertices of a cubic grid and are organized in layers forming regions characterized by one of the three possible states considered: proliferative (k = p), quiescent or hypoxic (k = h), or necrotic (k = n).
Proliferative cells may divide according to a probability binomial distribution, p div = p max g(d; d 0 , κ), where the nutrient gradient inside the MTS, g, depends on the depth of the cell inside the MTS, d, and the two parameters d 0 and κ. After each iteration, which corresponds to a complete cellular cycle, simple evolution rules, linked to the distance of cells to nutrients and O 2 , are used to simulate how these cells change from one state to another, updating their status in the aggregate. Exfoliation of the peripheral cells is also considered according to a binomial distribution with probability p exfol .
Region thicknesses are calculated as T k = Tn k , where T is a scaling factor that takes into account the cell size andn k is the average of the cell number in each of the three possible states, in the three orthogonal directions. With these thicknesses, MTS volumes are calculated as: Here, the f k are packing factors that take into account the intercellular space and the cell compression in each MTS region.
The feasibility of the model was tested by comparing simulated growth curves with the experimental data measured for an MTS sample of the MCF-7 cell line. The parameters p max = 0.37, p exfol = 0.01, and T = 15 µm, as well as the number of proliferative and hypoxic layers, L p = 6 and L h = 3, were chosen from the experiment [22,26,60]. The values of the packing factors, f p = 0.9, f h = 1.0, and f n = 1.1 were selected according to experimental evidence [61]. Finally, d 0 = 6, κ = 0.8, were fine-tuned to reproduce the experimental results obtained in [53] for the MCF-7 line.
The shapes of the experimental and simulated MTS are very similar, with the necrotic nucleus and the hypoxic and proliferating cell regions well described. The ABM also reproduces the linear relation experimentally found between the radius of the necrotic nucleus and that of the whole MTS: both the slope (b exp = 1.078 ± 0.096 and b sim = 1.006 ± 0.002) and the independent term (a exp = −151 ± 44 µm and a sim = −142.5 ± 1.6 µm) of the respective linear regressions coincide within the uncertainties (given with a coverage factor 2). Furthermore, the dependence of the number of hypoxic plus proliferative cells with the total volume of the MTS shows the same behavior as the experimental data [60], agreeing within uncertainties with them (see reference [53]).
Two samples with 1000 virtual MTS each have been generated with the ABM. These are the data to which the mathematical models have been fitted in order to test their predictive capabilities. In sample #1, the model parameters are left fixed at the values indicated above for all simulated tumors. In this way, only the variability linked to the Monte Carlo randomness of the different processes involved in the simulation has been taken into account. The variability observed in the MTS volumes measured in the experiment performed in [53] is within the statistical variability found in this sample #1 of simulated MTS. However, it has been realized that, in order to be fully described, some of the experimental details require an additional variability in the simulation. This is accounted for in sample #2, where the model parameters for each virtual MTS were sampled by assuming Gaussian distributions centered at the values considered in the first simulation and with widths equal to 10% of the respective centroids. Figure 1 shows the volume vs. time of growth plots corresponding to these two samples. Some of the simulated MTS are shown with solid lines. The whole region covered by the two samples are plotted with light colors. The large increase in the volume ranges of sample #2 produced by the extra 10% of variability in the simulation parameters is apparent. In addition, the crosses between the growth curves in sample #2 occur more often than in sample #1, thus reproducing the experimental observations better.
The initial configuration was the same for all simulated spheroids: seven cells occupying the center of the grid and the six near neighbor positions. Each simulated MTS has been followed for 60 days. Thus, sets of volumes {V i , i = 1, 2, . . . , 60} were obtained for each spheroid.
In [53], the experimental follow-up of the MTS gave us volume data each 2-3 days between 2 and 16 days after the cell sowing. By considering simulated data between days 1 and 60, we spanned a time range much larger than the experimental one, thus solving the drawbacks mentioned in the Introduction. It is worth mentioning that the simulations can be extended to larger time periods if necessary. The 60-day data available for each simulated MTS are enough to perform the analysis we have done.

Statistical Methods
The various mathematical models described in the previous section have been fitted to the volume vs. time data set obtained with the on-lattice Monte Carlo ABM for each one of the simulated MTS. These fits have been carried out by using the Levenberg-Marquardt method [62], in which the model parameters are adequately changed until the minimum of the function is reached. In the previous equation, t i indicates the time values at which the MTS volumes V i are obtained in the simulation, σ i labels the corresponding volume uncertainty provided by the simulation, and V(t i ) is the estimation of the MTS volume obtained with the specific mathematical model. The results of the fitting procedure include the fitting parameters that provide the best fit of the model to simulated data and their uncertainties that are calculated following standard procedures [62].
To measure the goodness of the fit, the figure-of-merit χ has been used; it is defined as: Here, ν = N − n par indicates the number of degrees of freedom in the fit that is calculated as the difference between the number of data in the sample to which the model is fitted, N, and the number of model parameters, n par . In this case, n par = 2 for the E-model and n par = 3 for the other models. The quantity χ represents the χ 2 per degree of freedom, and it is worth noting that good fits can be assumed if χ = 1 [62]. A first analysis consisted in evaluating the capabilities of the mathematical models to reproduce the behavior of the volume of the simulated MTS in the whole growth period. However, the spheroids showed significant volume changes in the first growth days, and the data were too noisy. To avoid numerical problems in the fitting procedure, the mathematical models were fitted to V i data in the interval i = . In this case, the figure-of-merit measuring the goodness of the fits carried out was labeled χ all .
To complete the analysis carried out, a comparison between the figure-of-merit χ reg , obtained when the various models were fitted to different regions of the MTS growth data, was carried out. Fits to subsets of volume data {V i , i = k, k + 1 . . . k + 19}, with k varying between 6 and 41, were carried out.
It is worth pointing out that, as the predictive capabilities of the fitted models are checked with data from the spheroid itself used to perform the fit, no type of normalization to common initial volume is necessary.
For all fits done, the distributions of the χ values obtained for the two samples of simulated MTS considered were calculated. In addition, the corresponding average values, χ, as well as the 95% confidence interval (CI), were obtained for each of these distributions. The 95% CI is given as [l : u], with l and u the lower and upper limits of the interval. Figure 2 shows the results of the fits to the volume vs. time simulated data obtained with the ABM for three different MTS of sample #1. Panel (a) corresponds to the MTS #278, which is the spheroid that is best fitted using the G-model. The results of this fit are shown with a dotted black line. The fits obtained with L-, P-, and B-models are shown with solid green, dashed blue, and dashed-dotted red lines, respectively. The E-model is analyzed independently due to its particular characteristics.

Results and Discussion
In the inset, the differences between these fits and the simulated volumes are plotted with open black squares, solid green circles, solid blue squares, and open red circles for the G-, L-, P-, and B-models. Panel (b) shows the results corresponding to MTS #179, which is the spheroid that is best described with the L-model. Finally, the results shown in panel (c) correspond to the MTS #138; this is the spheroid best described by both the P-and the B-models. In the table shown in Figure 2, the χ all values were obtained in all the fits of this figure.
G-, P-, and B-models provide a reasonable description of the data in all cases. However, it is evident that the L-model is not able to produce a good fit of all data simultaneously, even in the case of the MTS #179 (see panel (b)) that is the spheroid that best fits this model. This is even clearer if one looks at the χ all data shown in the table embedded in Figure 2. The values obtained for the G-model are an order of magnitude smaller than those found in the fits of the P-and B-models, while those found for the L-model are between ∼40 and ∼100 times larger.
It is interesting to discuss here about the exponential model. Figure 3 shows the results obtained when the E-model is fitted to the data of the MTS #1 (simulated in sample #1). The first point to be noted is that this model cannot give a reasonable description of the volume data. As it can be seen in the inset (where a semi-logarithmic plot is shown), these data do not show a constant slope as the E-model does. As a consequence, only partial sets of data can be described. In this case, the fits shown are those obtained when three different data intervals,  day,  day, and  day, are chosen.
Though the data in a given time period may be nicely reproduced, the E-model fails completely to describe the volume values outside that period. The fact that experimental data of tumor or MTS growth are usually available for relatively short time periods (due to the obvious difficulties inherent in the experimental monitoring of samples) may be one of the main reasons why this model has been widely used to describe tumor growth. is the spheroid, among all in sample #1, whose data best fit the potential and Bertalanffy models; the MTS #179 (b) is the spheroid whose data best describe the logistic model, and the MTS #248 (c) is the spheroid whose data best fit the Gompertz model. The values of χ all obtained in each case are also shown. MTS #1    shows the corresponding semi-logarithmic plot to emphasize the MTS volume behavior. Figure 4 shows the distributions of the χ all values obtained in the first analysis carried out. In this case, the mathematical models have been fitted to the whole sets of volume data obtained in the simulations. The fits to data in the sample #1 (those without variability in the ABM parameters) produced the distributions shown with dark colors, while those plotted with light colors (and outlined in black) correspond to the sample #2 in which a variability in the values of the ABM parameters was included in the simulations. A summary of these results can be found in the third column of Table 1 where the average valuesχ all and the 95% CI of the distributions shown in Figure 4 are given for the four models and the two MTS samples.
In the light of these results, the following facts deserve a comment: 1. For all models, the distributions corresponding to sample #2 were wider than those of sample #1. This is seen in Figure 4 and also in Table 1. The sizes of the 95% CI quoted in the table for sample #2 are larger than those of sample #1 by factors between ∼2 for the G-model and ∼5 for the P-model. In addition, theχ all found for sample #2 is larger than the corresponding for sample #1, the increase ranging between 0.2% and 17% for the G-and B-models, respectively. The variability in the AMB parameters that was considered in the simulation generating sample #2 produces more different growth shapes of the MTS and this results in better (those with smaller χ 2 all values) and worst (those with larger χ 2 all ) fits occuring independently of the specific mathematical model considered.

2.
In general, the best fits to the whole growth curves were provided by the G-model. As shown in Table 1, theχ all obtained for this model are an order of magnitude smaller than those corresponding to the other models considered. In the case of sample#1, the 95% CI obtained for the G-model is clearly below those found for the other models. For sample #2, there is a slight overlap between the CI of the G-and P-models.

3.
It is noticeable the large χ all values obtained for the L-model. Theχ all values are above 170 in both samples, and the respective 95% CI begin at 139 and 56, respectively. This again points out the fact already discussed in connection to Figure 2 about the difficulties of this model to describe the whole MTS growth.

4.
A curiosity of these results is that all the distributions are skewed to the right (to high χ all ) except that found for the P-model in the case of sample #1 where an almost Gaussian behavior is observed.  The predictive capabilities of the models were analyzed by looking at how well the models, fitted to reproduce a part of the growth curves, described the remaining volume data. The χ distributions obtained are shown in Figures 5 and 6 and summarized in Table 1 in terms of the corresponding values ofχ and of the 95% CI. In the two figures, the distributions shown in the left panels are those found in the fits, while those calculated for the predicted data are plotted in the right panels. Again, the distributions shown with light colors correspond to sample #2 while those with dark colors were obtained for sample #1.
The results obtained deserve several comments: 1.
The best fits were now much better than those found when all volume data were included. The maximumχ fit wasχ fit pro = 3.75, which corresponds to the L-model. These relatively low values are due to the fact that the number of data fitted is now 21 instead of the 55 considered before. It is worth noting that also the low limits of the 95% CI of χ fit were much lower than in the case of χ all . In all cases, the values obtained are close to 1, with the only exception of the 1.97 found for the L-model in the prospective prediction (see Table 1). This means that, a priori, all the mathematical models could produce a nice description of a more or less short data section of a given spheroid. However, as pointed out in Figure 4, not all the models can describe the whole growth curve of some of the simulated MTSs included in the samples. In this respect, it is important to point out that, in practice, these mathematical models are fitted to a small number of experimental volume data and for few spheroids or tumors in the sample, and the model that better fits the existing data may not be the one that best describes the subsequent growth [34].

2.
Contrary to what was observed in Figure 4, the χ distributions obtained for sample #2 are rather similar to those obtained for sample #1. This is corroborated by the results quoted in Table 1 and may indicate that the main ingredient controlling the results of the fit is the number of points to be fitted. An exception to this seems to occur for L-and P-models in the case of the prospective fits (see Figure 6c,e). The possible effect due the region of data fitted is discussed below.

3.
As can be checked by comparing the right panels in Figures 5 and 6, the models are more efficient in the retrospective predictions than in the prospective predictions.
As it can be seen, the distributions obtained in the latter case extend to values that are at least twice as large as those resulting for the former. The worst situation is that of the B-model in which the ratioχ pr pro /χ pr ret ∼ 13. This ratio is ∼2 for G-and L-models and ∼5 for the P-model. The retrospective predictive power of the models studied may be of interest in the clinic when the response of tumors to treatments is analyzed. In fact, the evolution of the regression of the irradiated tumors could be similar to looking back at the tumor growth and what has been observed in the present analysis could provide significant hints in that respect. Obviously, treatments may perturb the growth trends and more specific conclusions could be drawn when irradiation is included in the ABM simulation algorithm. Work in this direction is in progress.

4.
The best results in what refer to retrospective predictions were provided by the Bmodel for which the averageχ pr ret ∼ 23 (see Table 1). A slightly larger value, ∼30, was found for the G-model, whereas, for the other two models, much higher values are achieved, in particular for the L-model whereχ pr ret ∼ 480.

5.
The G-model seems to be the most robust in what refers to the prospective predictive capability. In this case,χ pr pro < 60 (slightly larger for sample #2 than for sample #1). For P-and B-models, this average was ∼300, while, for the L-model, values larger than 1000 were found. 6.
In any case, the values obtained for both χ pr pro and χ pr ret indicate that the predictive capabilities of the models considered here are discrete.  To finish the analysis of the performance of the various mathematical models with regard to the description of the simulated MTS, the goodness of the fits when they are done in different regions of the growth curves was studied. Figure 7 shows the average values of χ reg and the corresponding 95% CI obtained for the four mathematical models as a function of the initial day of the data region fitted (which include 20 consecutive data).
Solid circles and uncertainty bars correspond to sample #1, whereas solid curves and light colored areas stand for sample #2.
The following points must be noted: 1. Similar results were found for G-and B-models that showed average values tending to 1 with increasing t k . The former appeared to be slightly robust in this sense, showing a less noising behavior for t k < 20 day.

2.
For L-and P-models, a non-negligible dependence on the data region fitted occurred.
It is worth pointing out the fluctuating behavior with t k observed in both cases.

3.
This also results in the differences between the two samples studied being larger for L-and P-models than for G-and B-models.
The analysis that has been done may help to understand how tumors, particularly small tumors that could remain as surgical residues or micro-metastases, evolve before or during radiotherapy, or to investigate the effects provoked by the delays and interruptions that often occur throughout the treatment. In the case of delays in starting radiotherapy, the prospective predictive power of the models is of great interest. For interruptions, both the capacity for retrospection, if the trend of the growth is preserved, and prospection are relevant. These problems can be analyzed with the ABM model and, afterwards, it would be possible to see if the classical models are able to describe what happens in these cases. The E-model deserves a particular comment. As it has been shown, this model can reproduce the growth curves only when the data considered for the fit are limited to a short time period. On the other hand, its predictive capabilities are very poor. However, in the clinical applications of tumor radiobiology in radiotherapy, only the exponential growth model is used in practice. A relevant example of this is the management of interruptions and delays occurring in fractionated radiotherapy, an important clinical problem due to its effect on the treatment effectiveness [33]. In this case, the exponential growth model is used to account for the increase in tumor cells when the total treatment time is extended and to estimate the effects of this prolongation. As the sole sophistication, the so-called onset time is included. This time attempts to describe, somewhat crudely, the apparent period of slow growth of the tumors that is followed by an increase in the proliferation rate [63] and that occurs as a consequence of the reduction in the tumor cells produced by radiation after a few weeks from the beginning of the treatment [64].
It is evident that MTSs are much simpler systems than in vivo tumors. However, although MTSs are able to reproduce the early stages of the growth of the actual tumors, they do not include several of the elements that most influence the growth of the latter. The results of this work show that the E-model, even for such a simple "approach", is clearly inefficient for predicting tumor proliferation, and that, on the contrary, other models, such as the G-model, constitute a much more adequate scenario. The G-model also includes, in a simple way, the variation of the growth rate with the tumor size and, therefore, can account for the phenomenon of accelerated repopulation. Our results indicate clearly that it should be necessary to abandon the E-model in favor of any other, maybe the G-model, in clinical applications in order to produce more realistic estimations of the quantities calculated.

Conclusions
In this work, a detailed analysis of the capabilities of the Gompertz, logistic, potential, and Bertalanffy models to describe MTS growth curves was carried out. Spheroids simulated with an ABM recently developed were used as "pseudo-experimental" data. Two samples of MTS in which the ABM parameters were assumed to have a different variability were considered.
The best description of the whole growth data was provided by the Gompertz model that showed a value ofχ all an order of magnitude smaller than the other models.
In what refers to the predictive power, Gompertz and Bertalanffy models provided the best retrospective predictions, and the Gompertz model was the one giving the best prospective predictions.
Logistic and potential models had difficulties in fitting the data, even when restricted subsets were considered.
In principle, it is possible to find a certain (reduced) number of spheroids such that one of the models can be fitted to correctly describe a specific part of their growth curves. However, this fitted model does not describe adequately either the remaining volume data of the spheroids included in the fit or the data corresponding to other spheroids in the sample. Usually, experiments of MTS or tumor growth involve quite a few individuals that are followed during a relatively short period and this may explain the "success" of some mathematical models in describing the corresponding growth curves.
In summary, Gompertz model appeared to be the best option to describe the MTS growth data and showed to have enough flexibility to adapt its predictions to the volume growth data variability observed in the simulated spheroids. This indicates that this model should be considered in the clinical applications instead of the exponential model that is currently used. Funding: This work has been partially supported by the Spanish Ministerio de Ciencia y Competitividad (FPA2015-67694-P, PID2019-104888GB-I00), the European Regional Development Fund (ERDF), and the Junta de Andalucía (FQM0387, P18-RT-3237).

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data of the spheroids simulated with the on-lattice agent-based Monte Carlo model are available upon request from the authors.