A dynamic model for tumour growth and metastasis formation
 Volker Haustein^{1}Email author and
 Udo Schumacher^{1}
DOI: 10.1186/20439113211
© Haustein and Schumacher; licensee BioMed Central Ltd. 2012
Received: 7 January 2012
Accepted: 1 May 2012
Published: 5 July 2012
Abstract
A simple and fast computational model to describe the dynamics of tumour growth and metastasis formation is presented. The model is based on the calculation of successive generations of tumour cells and enables one to describe biologically important entities like tumour volume, time point of 1^{st} metastatic growth or number of metastatic colonies at a given time. The model entirely relies on the chronology of these successive events of the metastatic cascade. The simulation calculations were performed for two embedded growth models to describe the Gompertzian like growth behaviour of tumours. The initial training of the models was carried out using an analytical solution for the size distribution of metastases of a hepatocellular carcinoma. We then show the applicability of our models to clinical data from the Munich Cancer Registry. Growth and dissemination characteristics of metastatic cells originating from cells in the primary breast cancer can be modelled thus showing its ability to perform systematic analyses relevant for clinical breast cancer research and treatment. In particular, our calculations show that generally metastases formation has already been initiated before the primary can be detected clinically.
Keywords
Breast cancer Computational calculations Gompertzian growth function Tumour growth models Metastasis formationBackground
In the mathematically oriented medical literature different models are applied to describe the process of tumour growth and metastasis formation. Most of these models fall in one of the three following categories: The first ones are discrete models on the basis of single cell interactions which are then described by the aid of M^{te} Carlo simulations. The second ones are complex mathematical analyses of continuum models on the base of differential equations. A good overview of these approaches can be found in the articles of Ward and King [1, 2] and Roose, Chapman and Maini [3]. A third interesting alternate ansatz was developed by Iwata, Kawasaki and Shigesada [4, 5] which is in the following referred to as the IKSmodel. They model metastasis formation from the primary tumour and from metastases from metastases and give complex analytical solutions for the density respective the abundance of metastatic colonies depending on different growth functions of the primary tumour.
All the abovementioned methods have the disadvantage of complex reanalysis or the need for time consuming numerical recalculations when input functions or constraints are to be varied. Systematic investigations and the analysis of metastasis modulating events or treatment effects upon metastasis formation are limited due to the complexity or the computing power required.
In the following a mathematical model is presented which is based upon a series of successive generations of tumour development. This model enables a fast calculation of macroscopic relevant entities of the metastatic cascade. The entire programming was carried out in the C language using the graphical analysis package root, developed at CERN [6].
Results
The computational model
Metastasis formation is a complex process often referred to as a cascade as each step has to be performed in a certain order. It is initiated, when the first primary malignant cell starts to proliferate. If the developing primary tumour has reached a certain size, it sends out angiogenetic signals and blood vessels grow into the primary tumour. The future metastatic cell has to dissolve itself from the tumour mass by loosening of cell to cell contacts and has to degrade the basal lamina and the surrounding connective tissue. Having achieved this step in malignant progression, the future metastatic cell has to enter the bloodstream by migrating through the blood vessel endothelium. Once arrived in the circulation, the future metastatic cell has to survive in it and has to attach to the endothelium in the organ of the future metastasis. After attachment to the endothelial cell, the cell has to transmigrate through the endothelium and has to lodge in the stroma of the host organ. Presumably under the influence of local growth factors, the metastatic tumour cell has to proliferate in order to become a clinically detectable metastasis.
The characterized cascade can be effectively modelled by following this chronology of the events and making some realistic assumptions on the underlying distribution functions. This approach will be outlined in the following.
At each stage or generation of development a malignant cell inside a tumour has three possibilities: mitosis with doubling, apoptosis or migration into the next compartment where it becomes a potential metastatic cell. Each of these processes follows an exponential distribution with a characteristic constant λ_{a,m,d} = log(2)/T_{a,m,d}. With the restriction of no overlap in time, that implies that the 1^{st} started process will be executed, this results in a common exponential with λ_{G}=λ_{d}λ_{a}λ_{m} and a time per generation T_{G} = log(2.d)/λ_{G}. The fractions λ_{a,m,d}/λ_{G,} takes the values a,d and m and fulfil the constraint a + d + m = 1; the numbers are not necessarily constant over all considered generations. After n cycles this leads to (2.d)^{n} tumour cells. The number of potential metastatic cells is simply ∑(2.d)^{(n1)}·m. Either taking m(n) = m·δ^{n} or for calculation purposes more convenient leaving m constant and multiplying with a power of the actual number of cells, a metastasis formation process proportional to tumour volume V (δ=1), surface V^{2/3} or diameter V^{1/3} (δ<1) can be realized. Different interactions in the environment of the tumour or inside the lymphatic or blood vessel system will then lead to a finite life cycle of these disseminated cells either while being killed by the immune system respective by apoptosis or due to successful colonisation into the stroma of a peripheral organ. Again we assume an exponential distribution; now with the decay constant λ_{env}= λ_{k}+λ_{c}.
In continuation of the generation model with the time steps T_{G} we have to distinguish between cells that just enter the circulation and those that have already populated the blood or lymphatic system. The later group are surviving cells originating from former generations which had already entered circulation prior to the actual time step T_{G}. These cells will simply be successively reduced by a factor F = exp[−λ_{env}·T_{G}]. Accordingly the part (1F) will be eliminated from the blood system. In our model only the small fraction λ_{c}/λ_{env} of these cells will each colonise and develop metastases. The mean time point can be calculated by integrating the distribution function of such an exponential decay. The other group of cells, cells that just enter the circulation are subject to a different treatment. The process of creation and immediate elimination during the same time step T_{G} has to be accounted for. Number of surviving cells as well as mean time point and number of colonizing cells can be calculated by the combination of both the distribution function for dissemination into the blood stream and function for subsequent colonization of the stroma. Especially when the time scales for the life cycle T_{env} in the environment respective the generation time T_{G} differs significantly, this approach is necessary to calculate a more precise time of 1^{st} metastasis formation.
From the computational point of view simply a loop over N generations was generated, where each cycle generates the cell size of the primary tumour, number of disseminated cells in the blood or lymphatic vessels and the number of metastatic cells per T_{G} and in total at time n·T_{G}. Following the same strategy and by the use of recursion techniques the development in time of per T_{G} released metastatic cells and the process of secondary respective multiple metastases formation from metastases was calculated. If not indicated otherwise the further calculation is performed under the assumption, that metastases grow with the same speed and to the same maximal tumour size as the primary tumour.
Modelling tumour growth: the Gompertzfunction
In the following we will demonstrate the features of our straightforward strategy using the widely used Gompertzian growth function given by g(x)=μ⋅x⋅log(b/x). The parameter b is the asymptotic maximal reachable cell or tumour size and μ is the growth constant. Integration gives a tumour size of G(t) = b(1exp(μ∙t)). For metastasis formation a rate of the following form was taken by IKS: β(x) = γ·x^{α}. The parameter γ is simply the colonization coefficient and α stands for the fractal dimension of blood vessels infiltrating the tumour. In principle α denotes the fraction of tumour cells which participate in metastasis formation. For example α = 2/3 reflects a superficial angiogenesis of the tumour and dissemination occurs then notably only from the surface.
In the previous paragraph we proposed two models that describe a mitotic behaviour variable in time, but emanate from biologically totally different approaches. Both models show a Gompertzian like tumour growth and reproduce the metastasis formation of a given hepatocellular carcinoma. Before demonstrating the validity for breast cancer research we next show model systematic spread.
Systematic investigations
As can be deduced from the above given equation for G(t) the characteristic bending of the Gompertzian curve depends both on the maximal tumour size b and the rate parameter μ. On the other hand a constant μ translates in our model into a slightly b dependent initial T_{D}. For our calculations we choose a T_{D} of 5, 10 (comparable to the 9.8 days from IKS) and 20 days at a reference cell size of b = 10^{11}. To be comparable with the previous results dissemination from the primary occur proportional to V for the MS respective V^{2/3} for the GDRmodel. All rates and the lifetime of the tumour cell were chosen the same as in Figure 2, 3. A lifetime of 1 day is small or comparable to T_{D} and T_{G}; hence a realistic chance for colonisation into the stroma is given only for the time step being disseminated or the following one. A variation of the lifetime of a malignant cell within reasonable limits leads therefore simply to a scaling up of the combined rate of dissemination and colonisation. Our calculations confirm this assumption; systematic effects except those that can be seen by a variation of the combined rate were not regarded.
In the upper part of Figure 4 results are shown for the MSmodel for the three different developments of the tumour growth with time given by T_{D} = 5, 10 and 20 days (dasheddotted lines with black filled circles). In the lower part of Figure 4 the corresponding data for the GDRmodel are presented. To demonstrate that systematic differences between our models exist we include in both parts two further curves: The red lines with open circles represents a by one tenths reduced tumour growth rate, the blue lines and circles show the dissemination step following the next logical lower power of V^{m}, corresponding to the series volume, surface and diameter. For small maximal tumour sizes within the MSmodel, T^{1stM} takes a constant and depends only on the bending behaviour of the tumour growth curve with time. The graphs with the reduced rate and the lowered dependency to V^{2/3} underline this strong correlation. The value for T^{1stM} fits reasonable well with the time when the variation of the Gompertzfunction has reached its maximum, or when it is expressed with an equation when d/dt[dG/dt] = 0. Since the dissemination step is coupled to mitosis this corresponds to the timepoint when maximal metastases formation occurs. The probability for metastases formation is for small maximal cell sizes b in the percentage level and increases to 1 at around b = 10^{9} cells. It is obvious that the MSmodel would be able to describe a tumour entity which shows an extreme early but low rate of metastasis formation. In particular, in the next section we will argue that the MSmodel is a suitable candidate to describe a subgroup of breast cancer data. The monotonous decreasing data points above 10^{9} cells and the range b ≫ 5⋅10^{7} cells for the GDRmodel can be fitted with f(b) = −γ/μ⋅log(1log(β)/log(b)); a functional relation which is already known in a similar form from the combination of μ with T_{D}. The parameter γ is simply a scaling factor and β depends both on μ or T_{D} and the dissemination characteristics V^{m}. For the GDRmodel the comparison of the three graphs of T_{D} = 10 days but with different rates or dissemination characteristics emphasize again that a common asymptotic value for T^{1stM} will be reached if small tumour sizes are considered. Different to the MSmodel, the GDRmodel T^{1stM} does not correspond to the maximum of the variation of the Gompertzian growth. If the follow up time is sufficiently long, metastases formation will occur even if the primary does not change its size any longer as it has reached its maximum size. In contrast to the MSmodel where metastases start early and with low rate the metastases in the GDR model will colonize relative late but more frequent. This metastatic pattern reflects the fundamental differences between the two models. On one hand we assume a continuous prolongation of the tumour generation time T_{G}. This involves a naturally ageing of the cells with slower and slower running processes but with a regular and balanced sequence in mitosis and apoptosis. On the other hand we have highly active tumour cells; T_{G} remains constant but everything run with a high and lethal error rate. The fraction for doubling and apoptosis are shifted against each other which results in the decreasing tumour growth. Both models find their analogy in the biology of cells. It is known that a misbalance of anabolic and metabolic processes, reduced concentrations of enzymes or a failure in signal transduction is jointly responsible for an ageing of cells. Inadequate repair mechanism or missing stop signals in time during G_{0}phase of mitosis on the other hand leads to a slightly increasing of unformed and later apoptotic cells.
A clinical application in breast cancer
Probability and time of 1 ^{ st } metastasis formation
MSModel  GDRModel  

Dissemination from primary  Metastases growth characteristis  Single Cell size [10^{3} µm^{3}]  Metastasis Formation [%]  Mean time of 1^{st} colonizing cell [months]  Metastasis Formation [%]  Mean time of 1^{st} colonizing cell [months]  
Prim.  Meta.  Visible at pT2  Visible at pT4  Before pT1  Visible at pT2  Visible at pT4  Before pT1  
V  P  1  1  95.0  100.0  100.0  22.6  99.5  100.0  100.0  23.7 
2  1  93.3  100.0  100.0  22.1  97.0  100.0  100.0  23.3  
V^{2/3}  P  1  1  36.6  100.0  100.0  27.5  48.2  100.0  100.0  28.3 
2  1  32.3  100.0  100.0  27.5  40.5  100.0  100.0  28.1  
1  0.5  41.2  100.0  100.0  23.9  60.7  100.0  100.0  24.7  
A  1  1  10.6  87.9  100.0  16.2  5.5  22.7  71.6  13.4  
2  1  9.5  79.6  100.0  17.1  5.8  24.2  89.1  14.1  
1  0.5  10.3  90.5  100.0  14.5  6.0  25.2  77.7  11.9  
V^{1/3}  P  1  1  6.9  42.4  45.9  53.1  9.4  63.9  61.2  45.2 
2  1  6.6  44.4  59.2  51.2  8.4  58.2  76.6  45.2  
1  0.5  7.9  54.1  63.9  49.4  10.9  75.4  80.3  40.0  
A  1  1  4.0  17.3  98.1  29.5  2.5  6.1  15.2  22.7  
2  1  4.5  17.2  99.3  30.7  3.5  7.9  22.1  23.6  
1  0.5  5.0  20.2  99.8  25.8  3.5  8.7  19.9  19.6 
In the upper part of Figure 5 simulation results are shown for the MSmodel with dissemination proportional to V^{2/3} and V^{1/3}, the lower part shows the results for the GDRmodel, again with a V^{2/3} and V^{1/3} dependency. For both models a variation of the growth characteristics with an artificially slewed down growth behaviour of the metastatic cells are calculated and included in the plots. For the GDRmodel this growth behaviour was achieved by starting the metastatic growth process with the doubling rate of the time step the malignant cell has being disseminated. Taking the same generation time T_{G} and error rate for mitosis as for the primary tumour, a reduced initial doubling rate consequently leads to a reduced maximal colony size of the metastases. The modified growth function for the MSmodel was achieved by starting the metastatic growth just with the prolonged T_{G} of the generation the cell was disseminated. The ageing process is inherited. These colonies will reach the same maximal size as the primary but need for this growth significantly more time. Figure 5 clearly demonstrates that the metastases formation rate of 4.2% at time of initial diagnosis for pT2 patients and also the 21% for the pT4 patients can hardly be reached. Dependencies as V^{1/3} for the MSmodel or even lower ones for the GDRmodel are needed for the dissemination step, in order to achieve partly congruence between the clinical breast cancer data and our calculated probabilities.
Tumours of pT1 stage are in diameter just 3 times larger than our current clinical detection limit. To detect metastases of nearly equal size as the primary tumour indicates that the metastasis initiating cell must have been disseminated extremely early during tumour progression and even more importantly at a considerable rate of spread. Consequently the numbers of metastases will then increase nearly exponentially. Some mouse models [9] suggest an extremely early start of the dissemination process and are at least therefore in good agreement with our calculations. For reasons of clearness we recapitulate our assumptions: the whole data sample can be described by a single model and pT1 and pT2 mean age reflects the Gompertzian like tumour growth. Then it is mandatory, that either proportionalities ≤V^{1/3} should be taken into account to describe the low clinical probabilities for metastasis formation or the growth characteristics of metastatic cells should be different from that of the primordial stem cell initiating the primary tumour. Our models with the artificially reduced growth functions of the cells within metastases apparently come into the range of the clinical data. Due to the reduced maximal colony size in the variation of the GDRmodel metastases do not become large enough to be clinically detectable and the probability for the presence of metastases at stage pT2 is around 3%. At pT4 the probability falls below 10% where 21% was given for the MCR data. This indicates that the reduced growth as chosen for the GDR model somewhat underestimate the real growth characteristics. Anyhow it clearly demonstrates that a similar mechanism would be helpful to reproduce the data. More favourable is the situation for the MSmodel. Metastases that colonize distant organs at a later stage of malignant progression grow much slower than the primary tumour initially did. They remain hidden for a long time (“dormancy”). Under these assumptions we achieve a fair accordance with the data from the cancer registry, both pT2 and pT4stage probabilities are reproduced. The complete data set is summarized in Table 1, given are the mean values for the interval [7.5⋅10^{11}, 1.25⋅10^{12} cells] of the maximal reachable tumour size.
We have verified that a variation of the single cell size of both, primary and metastatic cells in the same direction do not lead to any noteworthy shift in the probabilities of metastasis formation. This is mainly due to the normalization at pT1stage. To show that our method is in general insensitive to the exact size of a tumour cell calculations where the quotient of metastatic and primary cell diameter varies by a factor ±2 are also given. Small systematic effects can only be seen for the standard growth behaviour of the metastatic cells. A reduction of the metastatic single cell size to half volume results in a relative 510% higher probability for primary metastasis formation at the pT2stage. This is a direct consequence of the reduced visibility during the normalisation procedure to 1.1% at pT1stage. On the other hand, doubling of the primary single cell size leads to a smoother bending of the tumour growth curve with a reduced growth rate. The increase in metastasis formation between pT1 and pT2 will therefore also be reduced. Again a 510% effect relative to the standard values can be seen.
In the following we will focus on those models in which metastases grow like the primary tumour. The time a primary tumour of the considered size (10^{12} cells) needs to reach the pT1 stage is single cell size dependent 49.4 month for cells of 2⋅10^{3} μm^{3} in volume, respective 54.9 month for 10^{3} μm^{3} cells. The mean time for the 1^{st} malignant cell to colonize into the stroma of the target organ lay around 23 month for the Vdependency and 28 and 45–53 months for the surface and diameter dependencies, respectively. All time distributions have a full width at half maximum of around 70%. These findings indicate that at least in half of the patients metastasis formation has taken place before the primary tumour became visible. The relative survival after 15 years with an over all metastasis formation rate was calculated to 77.6% for pT1 and 24.1% for pT4 [7], not differentiated by any treatment modalities. The MS as well as the GDRmodel with a V^{1/3} dependency for dissemination would best have the ability to explain the data if a reduction of metastases due to radio or chemotherapy is incorporated.
Discussion
We proposed a simple model of metastasis formation based on successive series of generations of tumour cells. With relative low computational power our model enables a fast insight into the growth and spreading behaviour of malignant tumours. The modelling itself is independent from the specific growth characteristics of a particular tumour. Here we concentrate on the Gompertzian growth and developed models rooted in the biological behaviour of malignant cells to describe such a growth function. Inside our framework we have demonstrated which mandatory implications can be deduced from the occurrence of metastases at a definite time. Especially the calculations based upon clinical data support the hypothesis that formation of metastases is a continuous and extreme early event during malignant progression. Our results are in good accordance with the analytical solution of Iwata et al. [4], who calculated their model according to a real clinical case. This accordance is remarkable, because we use a simple, straightforward simulation of successive generations of tumour cells whereas the IKSmodel is a complex solution for the development in time of the size distribution of metastases. Moreover we have demonstrated that our models should in principle be able to describe the breast cancer data of the Munich Cancer Registry as well. A combination of different Vdependencies for metastasis formation, a small but fast component that dominates the probability at pT1 stage and a slow V^{1/3} dependency for the further observed low numbers at pT2 and pT4 stages should be able to cover the whole range of growth and metastasis pattern. Additional and more detailed clinical data are however necessary before definite statements can be made.
Conclusions
A novel approach to simulate tumour growth and metastasis formation is presented. Within the framework of our model growth and dissemination characteristics of metastatic cells originating from cells in the primary tumour can be modelled. We adopted our model to clinical breast cancer data thus showing its ability to perform systematic analyses relevant for clinical breast cancer research and treatment. In particular, our calculations using these clinical data show that generally metastases formation has already been happened before the primary tumour can be detected with current clinical methods.
Abbreviations
 IKS = The authors Iwata:

Kawasaki and Shigesada
 MS:

Metabolic stagnation
 GDR:

Generation dependent rates
 MCR:

Munich Cancer Registry
 pT1:

pT2,pT3,pT4 = histopathological TNM classification of malignant tumours.
Declarations
Acknowledgements
We thank Prof. N. Shigesada for her comments while performing the numerical calculations of the analytical solution of the IKSmodel. We are especially indebted to Dr. J. Prahl for reviewing and deeply discussing our mathematical ansatz.
Authors’ Affiliations
References
 Ward JP, King JR: Mathematical modelling of avascular tumour growth. IMA J Math Appl Med Biol. 1997, 14: 3669.View ArticleGoogle Scholar
 Ward JP, King JR: Mathematical modelling of avascular tumour growth. II. Modelling growth saturation. IMA J Math Appl Med Biol. 1999, 16: 171211. 10.1093/imammb/16.2.171.View ArticlePubMedGoogle Scholar
 Roose T, Chapman SJ, Maini PK: Mathematical Models of Avascular Tumour Growth. SIAM Rev. 2005, 49: 179208.View ArticleGoogle Scholar
 Iwata K, Kawasaki K, Shigesada N: A Dynamical Model for the Growth and Size Distribution of Multiple Metastatic Tumours. J Theor Biol. 2000, 203: 177186. 10.1006/jtbi.2000.1075.View ArticlePubMedGoogle Scholar
 Struckmeier J: A Mathematical Investigation of a Dynamical Model for the Growth and Size Distribution of Multiple Metastatic Tumours. Hamburger Beiträge zur Angewandten Mathematik. 2002, 113.Google Scholar
 Brun R, Rademaker F: Root User’s Guide v5.26. 2011, http://root.cern.ch/drupal/content/usersguideGoogle Scholar
 Engel J, Eckel R, Kerr J, Schmidt M, Fürstenberger G, Richter R, Sauer H, Senn HJ, Hölzel D: The Process of Metastasisation for Breast Cancer. Eur J Cancer. 2003, 39: 17941806. 10.1016/S09598049(03)004222.View ArticlePubMedGoogle Scholar
 Hölzel D, Engel J, Schmidt M, Sauer H: Modell zur primären und sekundären Metastasierung beim Mammakarzinom und dessen klinische Bedeutung. Strahlenther Onkol. 2001, 1: 1024.View ArticleGoogle Scholar
 Schumacher U, Adam E: Lectin histochemical HPAbinding pattern of human breast and colon cancers is associated with metastases formation in severe combined immunodeficient mice. Histochem J. 1997, 29: 677684. 10.1023/A:1026404832394.View ArticlePubMedGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.