Cellular Automaton Simulation of Tumour Growth – Equivocal Relationships between Simulation Parameters and Morphologic Pattern Features

Objective: To develop an interpretation procedure which estimates simulation parameters (tumour cell motility, tumour cell adhesion, autocrine and paracrine growth control, stroma destruction) of simulated patterns solely based on morphometric features of the morphologic pattern. Methods: A cellular automaton computer simulation program was developed which produces morphologic patterns by growth of a seed of tumour cells. At the beginning of each simulation run certain simulation parameters are assigned to the tumour cells. After the run has been completed, the resulting pattern is evaluated by a set of morphometric features. Simulation parameters and resulting morphometric features of 27,800 simulations were stored in a database and were used for the evaluation of potential relationships. Results: Correlation analysis showed highly significant correlations between morphometric features on the one hand and the preset simulation parameters (tumour cell motility, tumour cell adhesion, autocrine and paracrine growth control, stroma destruction) on the other. Correlation coefficients, however, varied from 0.72 to 0.99. When only one simulation parameter varied while all others were kept constant, morphometric features yielded a highly reliable estimate of the particular simulation parameter. When variability was extended to 4 simulation parameters, morphometric features were less effective in estimating the setting of the parameters. Though in all patterns tested several possible simulation parameter constellations could be ruled out, morphometric features were usually compatible with more than one set of simulation parameters thus preventing a straightforward interpretation. Conclusions: Though simulation parameters significantly and reproducibly influence the resulting morphologic pattern as characterized by morphometric features, estimates of the simulation parameters based on morphometric features yield equivocal results.


Introduction
Histologic diagnosis and prognostic assessment of neoplasms is based on the experience that morphologic features are related to functional properties [9]. Both subjective qualitative evaluation and objective image analysis procedures rely on this basic assumption. Relatively little, however, has been done in order to elucidate the mechanisms which relate function to morphology in the context of tumour growth [21]. Knowledge about this relationship would provide two advantages: First, it would help us to understand why certain morphologic features are associated with outcome and would probably draw our attention to additional morphologic features which have not been taken into account so far. Second, functional aspects of tumour cell behaviour are usually assessed in highly simplified in vitro assays which do not necessarily reflect biologic behaviour within complex tumour tissues. The possibility to extract knowledge about the functional behaviour of tumour components from morphologic evaluation of the complex tumour tissue would provide additional insights in tumour biology in addition to those obtained from arteficial in vitro systems.
There are two possible ways to investigate relationships between function and morphology: On the one hand, morphologic features of a tumour with certain biological properties of the tumour cells can be assessed. This approach is particularly worthwhile in animal models of tumour growth and metastasis, where basic biological properties of the tumour cells can be assessed in vitro [7,24]. On the other hand, a theoretical approach is possible using computer simulation programs based on the cellular automata concept [5,22]. In these models, certain simulation parameters representing functional properties can be assigned to the tumour cell population, and the pattern resulting from continuous growth of this tumour cell population can be observed and can be described by quantitative pattern features [14].
Previous studies using computer simulations have shown that variations of simulation parameters, such as growth control, motility, adhesiveness and degradation are accompanied by variations in the resulting morphologic pattern [14,18]. There is also evidence that quantitative morphologic information collected from the final pattern may facilitate probabilistic estimates of the preset simulation parameters of the particular simulation run [18]. Application of this approach to experimental and clinical settings has demonstrated that these estimates derived from pattern analysis are consistent with functional properties of the tumour cells in vitro [13,16,23] and may be related to outcome in malignant tumours in humans [17][18][19][20].
The present study addresses the question whether an unequivocal relationship between simulation parameters and quantitative features of the resulting pattern exists, and whether estimates based on pattern features are able to reflect reliably functional parameters of the tumour cell population.

Material and methods
Computer simulation program. A cellular automata program (CATS = Cellular Automata Tissue Simulation) was developed using Visual Basic (Microsoft), running on an IBM-compatible image analysis system (KS 400, Kontron, Muenchen-Eching, FRG). CATS is a further development of the concept of the SMN5 program used in previous applications [16,18,19]. Simulation takes place within a grid of 100 × 100 cells, which can either be occupied by tumour cells or stroma elements. Each simulation starts with a seed of tumour cells in predefined distribution. Prior to each run, simulation parameters are defined corresponding to biological properties of tumour tissues. The cells of the initial seed and all offsprings are considered as a homogeneous population, i.e., that all cells in the simulation run show identical functional properties. In particular, the following simulation parameters were taken into account: Overall motility: Tumour cells were allowed to move. Overall motility was defined by the parameter MOTILITY denoting the average number of movements in relation to a single cell division. A single step was defined as a change of a cell between neighbouring positions in a randomly selected direction. MOTILITY = 0 means that the tumour cells do not show any active movement at all. MOTILITY = 10, for example, means that for any single cell division of a tumour cell another tumour cell carries out an average of 10 single step movements. Since there is no absolute timing in this simulation procedure, the degree of motility has always to be considered in relation to cell division.
Clustering of motility: It makes a difference whether tumour cell division is followed by small movements and again by division, or whether tumour cells move across a long distance followed by a period of repeated cell divisions. MOTILCLUSTER = 1 denotes that about each cell division is followed by movement, MOTILCLUSTER = 5, for example, means that after an average of 5 cell divisions a movement occurs with a number of steps defined by MOTILITY · 5.
Cohesion: Cohesion affects the place that is selected for a daughter cell when tumour cell division has taken place, and controls whether tumour-tumour cohesion or tumour-stroma adhesion is preferred. COHESION = 0 means that the daughter cell is deposited to a stromal place selected at random. CO-HESION > 0 means that the daughter cell is deposited on a stroma place which yields a maximum of tumour-tumour contacts and a minimum of tumour-stroma contacts. In contrast, COHESION < 0 means that a place is selected which yields a minimum of tumour-tumour contacts and a maximum of tumour-stroma contacts. Absolute values of COHESION between 0 and 1 denote the probability that an individual cell division takes cohesive properties into account. COHESION = −0.5, for example, means that about half of all cell divisions will choose a place which maximizes tumour-stroma contact, while the other cell divisions select places at random.
Autocrine and paracrine growth stimulation: Without autocrine or paracrine growth stimulation (GROWTHFACTOR = 0), each tumour cell has the same probability for cell division. The release of an autocrine growth factor by tumour cells (GROWTHFACTOR > 0) leads to an increased probability of cell division for tumour cells surrounded by other tumour cells, while the release of a paracrine growth factor by stroma elements (GROWTHFACTOR < 0) results in an increased probability of cell division in tumour cells surrounded by stroma elements. The dependence of probability on the surrounding has previously been described [11,15].
Stroma destruction: When a tumour cell divides, a stroma element is occupied by the daughter cell. The control parameter STROMADESTRUCT defines whether this stroma element is destroyed (i.e., replaced by a tumour cell) or whether this stroma element is pushed aside by the tumour cell. STRO-MADESTRUCT denotes the probability that in a particular cell division action the stroma element is destroyed, with STROMADESTRUCT = 1 indicating that the stroma element is always destroyed, and STROMADESTRUCT = 0 indicating that the stroma element is never destroyed but always pushed ahead. This "pushing ahead" takes place in a particular direction depending on the relative position of the tumour cell proliferation site and the site of deposition of the daughter cell. In this direction all following tumour cells and stroma elements are pushed forward for one grid place until the margin of the simulation matrix is reached. The direction of this "pushing forward" is defined by the relative position of the daughter cell to its maternal cell and can be one of 8 directions in the rectangular matrix. All elements are pushed forward for one grid place along a straight line oriented in the particular direction. At each pushing step, however, the involved elements are allowed to rearrange within a 3 × 3 matrix in order to keep the number of tumour-tumour, tumour-stroma and stroma-stroma interfaces constant. In brief, the number of tumour-tumour and tumour-stroma contacts of a particular element is recorded before this element is pushed forward. After pushing forward has been performed, the actual number of tumourtumour and tumour-stroma contacts is recorded. If these values differ from those before movement, the particular element is allowed to switch place with this one of the eight nearest neighbours which minimizes the alteration of tumour-tumour and tumour-stroma contacts. If this difference cannot be set to zero in any of the 8 neighbouring places, the residual difference is recorded and taken into account when the next element is pushed forward. This procedure has turned out to be necessary to maintain pattern integrity. Otherwise simple linear "pushing ahead" would result in complete disintegration of tumour cell clusters.
Sets of computer simulations. Each simulation was started with a seed of 20 tumour cells randomly scattered along the top row of the simulation matrix. This type of seed was chosen in analogy to the growth properties of cutaneous melanoma which starts with atypical cells lined along the dermoepidermal junction. Each simulation was terminated after 5000 steps of tumour cell division. In order to evaluate the influence of a single simulation parameter on the resulting pattern, 1500 computer simulations were carried out with GROWTHFACT = −3.00, MOTILITY = 10, MOTILCLUSTER = 5 and STROMADESTRUCT = 0, while COHESION was allowed to vary between -1 and +1. For statistical analysis, these simulations were divided into a learning set of 1000 simulations and a test set of 500 simulations. For assessing the influence of several variable simulation parameters, only MOTILCLUS-TER = 5 was kept constant, while COHESION varied from −1 to +1, GROWTHFACT from −3 to +3, MOTILITY from 0 to 20, and STROMADESTRUCT from 0 to 1. A total of 26300 patterns was simulated, with 24300 patterns serving as a learning set and 2000 patterns being used as a test set. Each simulation was terminated after 5000 steps of tumour cell division.

Measurement features.
Prior to the measurements, the area of interest was determined by applying a binary close operation with a 3×3 operator applied for five times, with the tumour cells used as "objects". Within this area, basic morphometric parameters (number of tumour cells, number of stroma elements, tumour cell density, surface density of the tumour-stroma border [25]) were calculated. Furthermore, mean and standard deviation of chord length of tumour as well as stroma components was assessed in X and Y directions [1]. City block distance mapping was performed in the tumour and in the stroma component, and the relative proportions of tumour and stroma components in the distance classes 1-10 were recorded [12]. Clustering of tumour cells was determined using 2×2, 4×4, 8×8 and 16×16 square boxes [3]. Affinity of tumour-tumour and tumour-stroma elements was calculated [6], and the fractal dimensionality of the tumour-stroma interface was assessed by a box-counting method [4]. Acronyms and definitions of the set of pattern measurement features are provided in Table 1.
Interpretation procedure. Two different types of interpretation were applied. At first, multivariate stepwise linear regression analysis [2] was applied in order to evaluate the relationship of the set of measurement features with each of the simulation parameters. There is, however, one major drawback to this procedure: It is possible that morphologically similar patterns may evolve in completely different settings of the simulation parameters, and this fact would not become appearent in linear regression analysis. Therefore, a second approach was chosen. For this particular purpose, each simulation parameter was classified by dividing the whole range into discrete categories. Cohesion, for example, was classified into values ranging from −1.00 to −0.33 (class 1), from −0.32 to +0.33 (class 2), and from +0.34 to +1.00 (class 3). In case of 4 simulation parameters with three categories each, each simulation run of the training set belongs to one out of 81 distinct groups. For each group, the minimum and maximum value of each measurement feature found in the training set was determined and stored in a database. Subsequently, each pattern of the test set was evaluated of being compatible with one or more of the 81 groups. A pattern was regarded as being compatible with a particular group, when each of the pattern features was within the range of the features defined for this group. By this approach, a list of "compatible" groups was determined for each pattern of the test set. In plain language, each of the groups found to be "compatible" represents a simulation parameter setting capable of producing a morphologically similar pattern.
Statistical evaluation. The "goodness of fit" in linear regression analysis was expressed by r 2 (squared regression coefficient) values. In the grouping procedure, the number of compatible groups found, the Table 1 Pattern interpretation in simulated tumour patterns using cellular automata. Acronyms and definition of pattern measurement features. All measurements are carried out within the "area of interest" defined by a closing operation of the tumour cells with a 3 × 3 operator applied for three times AFFINITY denotes the affinity of tumour cells to stroma elements calculated on the frequency of tumour-tumour, tumour-stroma and stroma-stroma contacts [6] FRACTDIM fractal dimensionality of the tumour-stroma interface calculated by box counting [4] frequency of the correct group found, the correlation of the real simulation parameter setting with the median of the compatible range, and the number of "inconclusive" patterns, where no particular parameter setting could be ruled out, was assessed.

Results
Experiments with one variable simulation parameter. In the set of simulations where all simulation parameters except COHESION were kept constant, clear differences in the resulting morphologic pattern depending on the value of COHESION were evident (Fig. 1). Multivariate linear regression analysis showed r 2 = 0.99, indicating that the degree of cohesiveness is almost completely expressed in the resulting pattern and is also reflected by the pattern measurement features. When the simulation parameter COHESION was transformed into a categorical variable, the quality of the parameter estimation clearly depended on the circumstances of the estimation procedure (Tables 2-4). When the number of categories was increased stepwise from 2 to 10, there is a definite increase in the quality of grouping of the test set: The number of simulations that appeared incompatible with any group decreased, the range of compatible values of COHESION decreased, and correlation between the median of the compatible range and the real value of COHESION on the other increased (Table 2). When the number of simulations in each group of the learning set was increased from 5 to 200, the percentage of patterns that appeared incompatible showed a dramatic decrease. The quality of the assessment of the patterns where a compatible group  was found, however, obviously did not depend on the number of simulations per group (Table 3). When the number of measuring features was increased by starting with a randomly selected feature, and adding two randomly selected at a time (Table 4), correlation of the estimate and the real value of COHESION   showed an increase at the beginning, but did not change considerably whether 5 or more features were used, indicating that the pattern measurement features are redundant. Experiments with multiple variable simulation parameters. Multivariate linear regression analysis showed a range of r 2 values ranging from 0.89 for GROWTHFACTOR to only 0.52 for STROMADE-STRUCT (Table 5). Obviously the interaction of various simulation parameters obscures the pattern changes which are due to changes of a single simulation parameter. The grouping procedure based on a learning set of 24,300 simulations divided into 3 4 = 81 groups revealed that many patterns were compatible with a wide range of simulation parameter settings. Mean number of compatible groups was 16 ± 10 out of 81 possible groups. When each of the 4 variable simulation parameters was evaluated separately, there was a considerable proportion of patterns which would have been compatible with the Table 6 Pattern interpretation in simulated tumour patterns using cellular automata. Variation of multiple parameters. Example of a case compatible with the whole range of MOTILITY and STROMADESTRUCT. Though no value of MOTILITY or STROMADESTRUCT could be ruled out, still certain constellations of MOTILITY and STROMADESTRUCT appeared to be incompatible with this particular pattern. In plain language, the particular pattern seems to be compatible with either low motility and a low degree of stroma destruction, or a high degree of motility and stroma destruction, but is incompatible with high motility together with low destruction or low motility with high destruction.  Table 6. Though the pattern appears to be compatible with the whole range of MOTILITY and TUMOURDESTRUCT, the combination of high MOTILITY with low TUMOURDESTRUCT -which would yield a less compact pattern -and the combination of low MOTILITY with high TUMOURDESTRUCT -which would not show individually scattered cells beyond the tumour bulk -can be ruled out.

MOTILITY
whole range of the simulation parameter: This was true for 9.6% of patterns of the test set as far as GROWTHFACTOR was concerned, but 78.4% appeared to be compatible with the whole range of CO-HESION (Table 5). When individual patterns, where the whole range of two simulation parameters was possible, were examined, it turned out, that nevertheless not all constellations of these two parameters were compatible, but that certain constellations could be ruled out (Table 6; Fig. 2). When for each simulation parameter the median of the compatible range was taken as an estimate of this parameter, correlation with the real value was slightly weaker than the correlation obtained by multivariate linear regression analysis (Table 5).

Discussion
Our study shows that tumour growth can be simulated using a computer model based on cellular automata concepts. Furthermore, objective evaluation of the resulting morphologic pattern is possible using several pattern measurement features. With this simulation model and the pattern measurement procedures, theoretical relationships between functional conditions determining tumour growth and the resulting pattern can be examined.
The simulation parameters of the simulation model were initially designed in analogy to basic mechanisms of tumour growth as known from in vitro and in vivo studies [8,10]. Pattern measurement features were those previously described for morphometric purposes [1,3,4,6,12,25].
The highly significant correlation between simulation parameters and pattern measurement features strongly indicates that the functional conditions of tumour growth reproducibly affect morphology of the evolving tumour. This significance was found when only one simulation parameter was altered, and also when 4 simulation parameters were simultaneously altered. With respect to the "goodness of fit", however, a clear-cut relationship was only evident in the simulation set where only one parameter has been altered, while in the simulation set where 4 parameters were altered the relationship was less clear. The attempt to estimate the preset simulation conditions based only on the knowledge of the final pattern harbours a considerable degree of uncertainty, as expressed by low r 2 values in multivariate regression analysis.
These results indicate that various simulation parameters interact when influencing the resulting pattern, and that it is often impossible to get an unequivocal idea of the underlying growth conditions when only the final pattern features are known. This finding raises the suspicion that certain morphologic patterns may be compatible with more than one possibility of simulation parameter settings.
In order to elucidate this problem more clearly, the learning set was divided into groups by classifying the simulation parameters. For each group, the range of values for each pattern measurement feature was determined. For each pattern of the test set all groups of the learning set were tested for "compatibility", i.e., whether all feature values of the test set pattern were within the range of the particular group. It turned out that in fact most patterns were compatible with more than one group of parameter settings, and that in many patterns no possibility of a particular parameter could be ruled out. Nevertheless, even in these cases only a subset of possible parameter combinations was compatible.
These data show that an unequivocal relationship between simulated growth mechanisms on the one hand and morphologic patterning on the other is only present when a single simulation parameter is changed. When several parameters change, only a stochastic relationship with a considerable degree of uncertainty is found. This relationship, however, is still of high statistical significance. Obviously in the context of the computer simulation program used, simulation parameter settings and morphology are statistically related and some hints as to the simulation parameters can be estimated when only the pattern features are known. Multivariate linear regression analysis is better than the grouping procedure in calculating an average parameter estimate, but the grouping procedure honestly discloses the limits of a functional interpretation of the static pattern.
When the procedure of estimating functional properties on the base of pattern measurement features is applied to real tumours, further limitations have to be taken into account: There are the limitations inherent to the simulation model, as, for example, the assumptions of a discrete, two-dimensional space, isotropy, and homogeneity of the cell populations. Vascularity has not been taken into account, since too many parameters (e.g., pre-existent vasculature, growth properties, branching, formation of anastomoses) would have to be arbitrarily set. Additionally, estimates which are derived from the knowledge obtained in the simulation sets may be influenced by the parameter distribution of these sets. Since certain parameter constellations may be more frequent in real tumours and other parameter constellations which were included in the learning sets may be completely absent in real tumours, estimates which appear "likely" in the simulation set may be rather misleading in real tumours. Therefore every application of such a procedure to real tumours can only be viewed and interpreted when the exact conditions of the simulation procedure and the learning set, from which this knowledge was derived, are disclosed and the associated limitations are kept in mind. The uncertainty may be narrowed down when additional knowledge about a real tumour exists, as, for example, the degree of proliferation or motility of the tumour cells in vitro.
In conclusion, the study shows strong statistical relationships between functional simulation parameters and resulting morphologic patterns, but also a considerable uncertainty when an attempt is made to estimate the simulation parameters solely based on the knowledge of the morphologic pattern. When applied to real tumours, cautious interpretation is mandatory.