Essential operating principles for tumor spheroid growth
© Engelberg et al; licensee BioMed Central Ltd. 2008
Received: 15 July 2008
Accepted: 23 December 2008
Published: 23 December 2008
Our objective was to discover in silico axioms that are plausible representations of the operating principles realized during characteristic growth of EMT6/Ro mouse mammary tumor spheroids in culture. To reach that objective we engineered and iteratively falsified an agent-based analogue of EMT6 spheroid growth. EMT6 spheroids display consistent and predictable growth characteristics, implying that individual cell behaviors are tightly controlled and regulated. An approach to understanding how individual cell behaviors contribute to system behaviors is to discover a set of principles that enable abstract agents to exhibit closely analogous behaviors using only information available in an agent's immediate environment. We listed key attributes of EMT6 spheroid growth, which became our behavioral targets. Included were the development of a necrotic core surrounded by quiescent and proliferating cells, and growth data at two distinct levels of nutrient.
We then created an analogue made up of quasi-autonomous software agents and an abstract environment in which they could operate. The system was designed so that upon execution it could mimic EMT6 cells forming spheroids in culture. Each agent used an identical set of axiomatic operating principles. In sequence, we used the list of targeted attributes to falsify and revise these axioms, until the analogue exhibited behaviors and attributes that were within prespecified ranges of those targeted, thereby achieving a level of validation.
The finalized analogue required nine axioms. We posit that the validated analogue's operating principles are reasonable representations of those utilized by EMT6/Ro cells during tumor spheroid development.
Extensive study of EMT6/Ro (hereafter EMT6) multicellular tumor spheroids grown in culture has provided useful insight into important aspects of tumor growth and tumor cell culture models. The behaviors of EMT6 cells in culture fall reliably within narrow ranges, as if cell behavior and thus the underlying mechanisms are tightly choreographed. Those actions can be thought of as being constrained and guided by a set of genetically specified biological operating principles. Can we discover and attribute a small, robust set of operating principles that combine to create the system level phenomena that characterize EMT6 growth in vitro? How can we represent and challenge those operating principles? What organization of the subcellular molecular biology enables the operating principles to emerge, and be sustained at the cellular level? Before addressing the last question, we need answers to the first two, which has been the objective of this project.
Efforts to model tumor spheroid growth characteristics (see , two recent reviews [3, 4], and references therein) have been extensive, informative, and successful. However, no one has proposed a cohesive set of cell level operating principles. Only recently has it become feasible to design and instantiate quasi-autonomous, cell mimetic analogues, [5–7] capable of exhibiting a rich phenotype of their own. The focus of most modeling and simulation efforts has not been in that direction. It has been primarily to provide precise, mostly mathematical descriptions of growth dynamics in terms of measured biochemical and physical factors combined with detailed descriptions of essential cell processes. The resulting models have been successful in explaining the time course and limits of spheroid growth in terms of nutrient depletion , increased acidity near the spheroid's center , and the dynamics of tumor spheroid metabolism . Jiang et al. combined these features into a comprehensive model that separately considered each cell and spanned three levels of mechanistic resolution . Other modeling efforts such as [11, 12] have used hybrid mathematical and individual based approaches that have shown initially promising qualitative results.
Cells consume resources, change state, proliferate, lose adhesion, die, shed, and move.
Cells proliferate throughout the duration of growth of the EMT6 spheroid.
Cells behave autonomously and locally.
The EMT6 spheroid develops an inner necrotic core, a middle quiescent layer, and an outer proliferating layer.
The EMT6 spheroid initially grows exponentially, then linearly, and then stabilizes.
The EMT6 spheroid has different growth characteristics at different levels of nutrient.
Necrosis onset occurs when the EMT6 spheroid has an area of roughly 0.2 mm2 at high nutrient and 0.02 mm2 at low nutrient.
The viable rim has a width of roughly 240 μm at high nutrient and 60 μm at low nutrient.
The measured initial doubling times are roughly 22 hours at high nutrient and 26 hours at low nutrient.
The mean error percentage between EMT6 spheroid and SMS growth is within 15% at high and low nutrient levels.
A cohesive set of operating principles (as distinct from isolated principles) can provide a framework into which more detailed, subcellular and molecular level information can be connected directly to system level phenotype. The plan was to work backward from a targeted set of in vitro observations of EMT6 cell and spheroid phenomena to a plausible set of analogue AXIOMS, which would be necessary and sufficient to generate in silico counterparts of the targeted phenomena. With that vision, this project has been motivated by three expectations: 1) Understanding hypothesized mechanisms in vitro would be facilitated by successfully building and studying analogous mechanisms in silico. 2) Achieving and refining validated analogues would offer a scientific, experimental approach to discovering and studying cohesive sets of operating principles. 3) Knowledge of axiomatic operating principles would facilitate exploration of their biological counterparts. This article reports on the design and implementation of the analogue, and the results of its execution.
Parameter names, values, units and sources.
In silico value
In vitro value
Proliferating NUTRIENT critical level (proNut)
3.0 × 10-3
3.0 × 10-19 mol/μm3
Quiescent NUTRIENT critical level (quiNut)
8.0 × 10-4
8.0 × 10-20 mol/μm3
Proliferating CELL'S NUTRIENT uptake (proConsumeRate)
5.0 × 10-4/SEC
5.0 × 10-17 mol/(cell s)
Quiescent CELL'S NUTRIENT uptake (quiConsumeRate)
1.0 × 10-4/SEC
1.0 × 10-17 mol/(cell s)
Delay before dead CELL is removed (removeDelay)
3.6 × 104 SEC
1.8 × 104 s
Movement bias (moveEmptyBias)
Delay between CELL creation events (prolifDelay)
Proliferation bias (proBias)
NUTRIENT diffusivity (diffusionRate)
Initial NUTRIENT concentration (initialVal)
0.165 or 0.008
16.5 mM or 0.8 mM
1 grid space
Average cell cycle**
~4.24 × 104 SEC
~4.24 × 104 s
Average time of removal after cell death**
1.8 × 104 SEC
1.8 × 104 s
In vitro source
NUTRIENT > proNut
Switch to PROLIFERATING state
Cell quiescence is regulated by the glucose and oxygen supply .
Switch to NECROTIC state
Cell death is regulated by the glucose and oxygen supply .
quiNut <NUTRIENT <proNut
Switch to QUIESCENT state
Cell quiescence is regulated by the glucose and oxygen supply .
State = PROLIFERATING or QUIESCENT
Consume NUTRIENT equal to proConsumeRate or quiConsumeRate
Cells consume oxygen and glucose at varied levels .
State = NECROTIC; removeCounter < 0
Necrotic cells eventually break up and are consumed .
Inside CELL adjacent to empty space
Move into empty space
Cells move and mix with other cells within the spheroid .
Outside CELL adjacent to empty space
Move into empty space with prob. pm
Cells move and mix with other cells within the spheroid .
Outside CELL with 0 neighbors
Randomly move in space
Cells can be shed from the exterior of the spheroid .
State = PROLIFERATING; prolifCounter < 0; CELL has empty neighbors
Create new CELL with prob. pb
Cells create new cells within the SMS, causing it to increase in size .
In silico growth curves matched in vitro growth curves
For the parameter values listed in Table 2, SMS growth curves were quantitatively similar to those of EMT6/Ro spheroids for both high and low nutrient conditions. CELLS within an SMS proliferated initially at an exponential pace. Growth then slowed and became linear because only CELLS near the outer SMS rim could reproduce. The increase in cross-sectional area was linear until CELLULAR NECROSIS began. Thereafter, SMS growth rate began decreasing toward zero. A stable size was reached when CELL creation was balanced by CELL removal. Plots of SMS cross-sectional area over time (Fig. 3) closely mirrored EMT6 spheroid growth  for both high and low levels of NUTRIENT. For simplicity, as discussed under Methods, we conflated measured concentrations of glucose and oxygen, along with the other in vitro nutrients, and represented the entire collection using the factor NUTRIENT. High NUTRIENT level mapped to 16.5 mM glucose and 0.28 mM oxygen. Low NUTRIENT level mapped to 0.8 mM glucose and 0.07 mM oxygen. The reported coefficient of variation of mean cross-sectional EMT6 spheroid area between multiple in vitro experiments was roughly 29% at 7 days, increasing over time . Given that, and the fact that EMT6 spheroids increase their size by many orders of magnitude during growth, we judged that having simulated values within 15% of referent values would be reasonable, and made that a targeted attribute (Table 1). The mean percent error between in silico and in vitro data was 12% for high and 8% for low NUTRIENT, which was within the targeted 15% range. The only parameter changed between the two conditions was initialVal, the level of NUTRIENT present outside the SMS during the simulation. For both conditions, CELLS used the same Table 3 set of operating principles. Under Methods, we describe that only minimal tuning of the indicated subset of the parameters in Table 2 was needed to achieve these matching growth characteristics.
In silico doubling times were similar to in vitro doubling times
Comparison of in vitro and in silico growth characteristics.
Initial doubling time
Viable rim width
Necrosis onset time
Necrosis onset size
21.4 HOURS/19.2 HOURS‡
1.61 mm2 ‡/1.46 mm2‡
40.0 HOURS/19.2 HOURS‡
16.5 mM glucose & 0.28 mM oxygen
23.0 hours†/21.6 hours‡
3.25 mm2 †/2.79 mm2‡
0.8 mM glucose & 0.07 mM oxygen
17.0 hours†/26.4 hours‡
0.221 mm2 †/0.0725 mm2‡
Measured viable rim widths were similar
Viable cell rim widths, eighth in Table 1, have been characterized, and were used to further validate SMS attributes. The data in Table 4 show that VIABLE SMS rim widths were close to in vitro values. Because the coefficient of variation of EMT6 spheroid areas was at least 29%, the corresponding value for radius was roughly 15%. We specified that any SMS radius within 15% of a referent radius would be acceptably similar because that radius would be experimentally indistinguishable from a repeat EMT6 experiment, had one been preformed. We observed a mean SMS rim width under the high NUTRIENT condition that mapped to 245 μm, compared to 240 μm in vitro, a difference of 2%. At low NUTRIENT, mean SMS rim width mapped to 62 μm, compared to 60 μm in vitro, a 3.3% difference. Although we performed some tuning of the critical levels required to remain in the PROLIFERATING or QUIESCENT state, these similarities are still noteworthy. They reinforce the likelihood that the principles of operation used by SMS CELLS may map to a corresponding set of operating principles used by EMT6 cells.
NECROSIS onset and final saturation size were similar
Under high NUTRIENT, measures of SMS diameters at NECROSIS onset, listed in Table 4, achieved the targeted similarity measure. They were within 15% of those observed by Freyer and Sutherland . Following SMS execution under high NUTRIENT, mean diameter at which the SMS first underwent NECROSIS mapped to roughly 530 μm, compared to 516 μm for the EMT6 spheroids under comparable conditions, a 2.7% difference. Under the low NUTRIENT condition, SMS underwent NECROSIS when the system reached a diameter corresponding to approximately 180 μm, compared to 152 μm for EMT6 spheroids, an 18.4% difference.
The maximum sizes attained by SMS were not similar to those predicted by Freyer and Sutherland , (Table 4). Freyer and Sutherland did not measure maximum sizes, but instead inferred them by fitting data to the Gompertz equation and then using the fitted equation to predict an expected maximum size. The data fit were averages of results from experiments on different batches of EMT6 spheroids. Because the SMS were being compared to data from a single experiment, we fit the Gompertz equation to that referent data (Fig. 3). Table 4 shows that the new result did not differ significantly from the one originally reported at high nutrient levels: the moderate in silico-in vitro discrepancy remained. However, the maximum fit size was smaller at low nutrient concentrations, down from 0.221 mm2 to 0.0725 mm2. The maximum size reached by the SMS at low nutrient mapped to 0.0645 mm2, a difference of 11%. That was judged acceptably similar to our Gompertz equation fit. A reasonable conjecture for the discrepancy at high nutrient levels is that the set of operating principles used by cells in maturing EMT6 spheroids were somewhat different than the set used earlier, during spheroid expansion. Our goal was to seek one set of SMS operating principles that would enable validation for both high and low NUTRIENT conditions. It would be straightforward to relax that requirement and achieve improved similarity at high NUTRIENT levels.
SMS shape and stability were controlled by proBias
SMS long-term shape changes lead to instability
Varying parameters changed growth curve and SMS shape
Effects of increasing parameters on in silico measures.
Viable rim width
Necrotic core size
As demonstrated by Fig. 7, when moveEmptyBias was set to zero, SMS grew linearly at a high rate in low NUTRIENT conditions and failed to saturate. Fissures appeared, which caused the SMS to destabilize and grow chaotically. Increasing moveEmptyBias by as little as 0.25 resulted in almost complete SMS saturation. Further increasing moveEmptyBias did not significantly affect growth rates or stability, but small changes were evident at both low and high NUTRIENT levels. While moveEmptyBias does not directly map to an in vitro quantity, these results indicate that there may be threshold values for shape maintenance mechanisms, below which an EMT6 spheroid would generally become unstable.
SMS events and mechanisms were not intended to be exact replicas of the actual physical or chemical events ongoing in vitro during EMT6 spheroid growth. Nor were predictions of specific events part of the intended SMS use. Rather, the intent behind our method has been, given a set of EMT6 spheroid attributes, to discover SMS computational mechanisms that might map logically and intuitively to in vitro counterparts. This has been accomplished by exploring the inverse map from phenomena to mechanism. The primary functional unit of an SMS – a CELL – does map 1:1 to an EMT6 cell. Because an EMT6 cell is autonomous, we designed SMS CELLS to be quasi-autonomous. SMS CELLS currently have no internal components. As atomic software objects, they needed operating principles to function. Most of the principles that cause an EMT6 cell to act in a particular way when faced with specific circumstances in culture were unknown. Consequently, we needed to discover and implement operating principles that each SMS would use, evaluate those mechanisms through simulation and observation, and modify them based on the results. Following [5, 7], CELL operating principles were formulated as AXIOMS. Their specifications were tightly guided by available knowledge of EMT6 behaviors in culture [17–20]. By iteratively following the diagram in Fig. 1, we narrowed and refined early candidate AXIOMS to nine. These AXIOMS were refined further so that measures of SMS growth characteristics would match prespecified, iteratively expanded, targeted sets of EMT6 spheroid growth characteristics according to specific similarity measures. Having achieved that objective, we suggest that the resulting SMS operating principles (Table 3) can stand as an abstract representation of EMT6 operating principles under comparable growth conditions. We posit that the larger the targeted set of EMT6 attributes satisfactorily matched, the more realistic the mapping between SMS and EMT6 operating principles.
It is significant that within a simulation cycle one CELL can apply more than one AXIOM. This reflects the complexity inherent in even the simplest interpretation of a biological system. The amount of nutrients or growth factors in the environment, for example, can be independent of whether a cell is surrounded by other cells or isolated.
Use of AXIOM 8 in combination with the others resulted in an extreme degree of contact inhibition: CELLS that were surrounded by other CELLS did not create new CELLS. That was a purposeful simplification. Nevertheless, the targeted attributes were achieved. The evidence indicates that some cell proliferation does occur throughout EMT6 spheroids , but that the frequency decreases dramatically with distance from the surface. If those observations were to be added to the list of targeted attributes, it would falsify the current SMS. Validation against that expanded attribute set would require increasing SMS complexity, possibly revising, as well as extending the list of AXIOMS. Relative to the current SMS, the fraction of CELL creation events occurring at the surface would be reduced and counterbalanced by division events occurring elsewhere. Because SMS components are quasi-autonomous, when the current set of targeted attributes is expanded one at a time, it is relatively straightforward to revise an SMS to match each new, expanded set. We achieved the targeted attributes using an SMS CELL that exists in three states. When the attribute list is expanded (even to include pathological attributes of drug treatments), it is straightforward to add new CELL states that possess different axiomatic operating principles.
The current set of abstract, axiomatic operating principles is believed to be the source of the discrepancy between in silico and in vitro growth at high NUTRIENT (Fig. 3). The SMS can be parameterized so that simulated growth more closely matches the higher NUTRIENT data (not shown), but at the expense of achieving a much poorer match to the low NUTRIENT data. Note that the differences in growth properties at low and high concentrations of oxygen and glucose are more extreme for the referent data than is seen with other available sets of growth data, such as the data used by . Achieving a tighter match would require adding more detail.
Whereas Freyer and Sutherland described the inhibitory actions of a tumor extract on proliferating cells , they did not separate the components to identify the source of inhibition. LaRue et al.  observed cyclin-dependent kinase inhibitors that are associated with cell-cycle arrest, but they did not demonstrate a causative role. Researchers have speculated that a factor in spheroid growth stabilization may be cell inhibition caused by some material being released from necrotic cells [14, 19]. We did not include such an attribute among those targeted, in part because it had not been confirmed. Nevertheless, the current SMS successfully produced stable spheroids without the production and action of such a factor, effectively establishing that one is not required for growth stabilization at biologically realistic SMS parameter settings. Of course, we cannot conclude from this in silico evidence that a necrotic inhibitor is absent in vitro. It is instead evidence that EMT6 spheroid growth stabilization need not require the presence of such an inhibitor. It is also useful to contrast the modeling approach used here with that used by [21, 22]. Longo et al., having achieved some degree of satisfaction about the mechanisms implemented, focused on replication and prediction of referent results from particular AXIOMS in an exploration of the forward map from generator to phenomenon. Our approach focused on discovering appropriate AXIOMS, such as the need for a potential inhibitor or a particular arrangement of neighboring components, and which were necessary and/or sufficient. We relied on falsification to select from the plausible generators.
CELLS that experienced a high STRESS were likely to move to reduce STRESS, while CELLS experiencing low STRESS were likely to proliferate and create more CELLS. As shown in Fig. 5, some large SMS destabilized during long-term growth. We determined that this behavior was caused by the probabilistic, local nature of the STRESS based movement and proliferation algorithms. At small SMS sizes, all deviations from the minimum-STRESS, convex curvatures are corrected by the movement and proliferation algorithms within a small number of simulation cycles. For much larger sizes, however, local curvature can be within the variability of the STRESS algorithm, yet the shape that emerges can be non-circular and irregular. That is because all AXIOM preconditions used only local information. When SMS are very large, the surface adjacent to every surface CELL can be relatively flat (the CELLS are experiencing low STRESS), yet the overall SMS can be non-circular. If needed, the effect could be minimized in several ways, all of which would require increasing SMS complexity. The simplest for the current SMS design would be to enable sharing information about each CELL'S current STRESS with a larger cluster of neighbors.
Although other models have not explicitly controlled spheroid shape, they have nevertheless done so implicitly. For example, by placing an adhesion term in their models, Schaller et al. and Jiang et al. caused CELLS to cling together, thus minimizing surface irregularities [2, 6]. In fact, Schallar et al. noticed differences in overall shape when they used different values for the adhesion parameter. Anderson et al. found that changing the EXTRACELLULAR MATRIX structure in a simulated model of tumor invasion produced dramatic differences in tumor morphology . Our analogue did not initially contain a mechanism to control SMS shape, but we found that the analogue could not mimic the targeted attributes without one. Although the SMS did not explicitly define and implement cellular adhesion like [2, 6], stress based movement and proliferation produced a similar effect.
We presented an idea: under the conditions of EMT6 spheroid growth in culture, molecular cell biology manifests at the cell level in what can described as a small set of operating principles that are responsible for the characteristic in vitro phenotypic attributes. We anticipated needing to identify and understand the operating principles in order to better understand how specific, detailed subcellular events may be linked to attributes of systemic EMT6 spheroid growth. Our method and approach are diagrammed in Fig. 1. We designed, refined, and tuned quasi-autonomous software components that, upon execution, formed abstract SMS analogues. We showed that measures of SMS behaviors during simulated growth were similar to available wet-lab data using a quantitative similarity measure. We submit that SMS mechanisms, with emphasis on the explicit AXIOMS, may stand as a plausible, abstract hypothesis for what was observed during those EMT6 spheroid growth experiments.
A future challenge will be to build a parallel system in which each (or some) atomic CELL component and its operating principles are replaced with a composite CELL object containing a set of interacting components intended to map to modular components within EMT6 cells. During INTRACELLULAR interaction, specific internal components would each use a portion of the same local environment information to act on other internal components such that actions are essentially identical to the current SMS events. The resulting growth characteristics would be indistinguishable from those described herein. The two systems could be iteratively advanced in parallel as new information and data were added to the set of targeted attributes. Using cross-model validation in that way is expected to provide a systematic strategy to answer the third of the three questions posed in the Introduction. What organization of modular and molecular biological details enables operating principles to emerge, and be sustained at the cell level?
In vitro system: historical context
Because cancer is such a complex and heterogeneous disease, researchers develop and study model systems. One is the in vitro EMT6/Ro multicellular tumor spheroid system. Freyer and Sutherland used the system to study avascular cancers in the 1970s and 1980s [13, 19, 24–27]. Their initial goal was to create a system that would allow many EMT6 spheroids to be grown in the same flask under identical, controlled conditions. Study of that model was expected to improve our understanding of how early stage cancer forms and improve our ability to treat it when the cells have not reached total genomic instability and still have much in common with normal cells .
In order to obtain adequate numbers of spheroids for measuring growth curves, experiments employed spinner flasks containing hundreds of spheroids. The cultures were initiated in monolayer and then grown in dishes until small spheroids were present (95–100 hours). These spheroids (usually 400–600 cells) were sorted and transferred to flasks, which contained a solution of glucose-free Eagle's Basal Medium, Fetal Bovine Serum, and varied concentrations of glucose. Oxygen was bubbled through the flask, and glucose was replenished roughly every 10–14 hours . Most early experiments were designed to characterize the system and its behavior. Eventually, however, many researchers shifted to using the system as a tool rather than studying the system itself. The seminal studies describing the behavior and characteristics of EMT6 spheroids were primarily completed by 1992. Important work continued nevertheless. There was an effort  to identify a potential necrotic inhibitor, and confocal microscopy was used to assess growth fractions .
Targeted attributes describe in vitro EMT6 spheroid growth
Although there is variation in EMT6 spheroids' growth, it reliably follows the same well-defined pattern . EMT6 spheroids initially grow exponentially without constraint from nutrient or other cells. This gives way to linear growth as cells become quiescent due to nutrient depletion within the spheroids. Eventually the spheroids begin stabilizing, both in volume and cell number, though cells continue to reproduce on the outer edge. They develop a concentric layered structure: an outer layer of actively proliferating cells, a middle layer of quiescent cells, and a core of necrotic cells and cellular debris. The material released by dying cells is thought to inhibit cellular proliferation , but it is not clear if this material actually affects cells in the rim and thus spheroid growth rates. During growth, an EMT6 spheroid maintains a generally spherical shape, but neither shape nor relative EMT6 spheroid stability have been quantified, which makes shape validation difficult. Additionally, though the width of the viable rim has been measured for different cell types, only one group has attempted to quantify the ratio of proliferating and quiescent cells . Wartenburg and Acker performed these measurements on human glioma spheroids. They do form spheroids with concentric layers, but have quantitative growth characteristics (such as initial doubling time and saturation size) that are different from those of EMT6 spheroids. We elected to falsify and validate our experimental results using data from Freyer and Sutherland , because they performed the most thorough and complete analysis of EMT6 spheroid growth, including measuring growth curves, the width of the viable rim, and the size at the approximate time of necrosis onset.
Analogue construction within an agent-based paradigm
SMS construction used agent-based methods available in the simulation toolkit MASON . This framework was used for data generation, scheduling, and visualization. The SMS is an example of a class of simulation models referred to as executable biology [31, 32]. Our simulation shares some similarities with  and is closely related to , though the system under study and the simulation framework are distinct. Executable biological analogues are poorly suited for precise prediction, but are ideally suited for testing hypotheses about mechanisms. The basic method requires building mechanisms at the functional unit level closest to the targeted phenomena. Here, that unit is the cell. An SMS is comprised of quasi-autonomous agents. Each maps to an EMT6 cell. The initial limit for SMS resolution was the CELL. If achieving the list of targeted attributes required doing so, the resolution of the SMS could be increased. CELLS interact with each other and their environment during each simulation cycle within a two-dimensional, hexagonal grid. Earlier versions used a square grid, but in addition to requiring a higher order implementation of discrete diffusion, it also generated artifacts that are not present with the hexagonal grid. We tested different orders of discrete diffusion to verify that artifacts were not caused by the diffusion algorithm. CELL actions are mandated by AXIOMS[5, 7]: when a specified condition is met, a specified action occurs. Together, these AXIOMS are a CELL's operating principles. A goal has been to improve the variety of SMS attributes that are similar to corresponding EMT6 spheroid attributes: the expectation being that with increasing phenotypic similarity, the higher the likelihood SMS AXIOMS will map to corresponding EMT6 operating principles (Fig. 1). As the list of targeted attributes expands, analogue resolution can be adjusted as needed.
A second grid, adjacent to the CELL'S grid contains a diffusible substance called NUTRIENT. NUTRIENT adjacent to each CELL is detected by and available to that CELL for consumption. During each simulation cycle, each CELL uses the AXIOMS in Table 3 to select actions based on how its local environment has changed since the last simulation cycle. Examples of actions include move, change state, create new CELL, DIE, and shed. AXIOMS are implemented by algorithms that utilize the parameter values listed in Table 2. Where appropriate, parameter values intentionally mirror values measured in vitro. The remaining parameters were tuned using an iterative process: change parameter value, execute, evaluate relative to referent observables, cogitate, change again, etc.
Tuning parameter values improves the analogue's ability to survive falsification
Once a set of targeted attributes had been specified, parameter tuning and SMS validation became closely linked. When seeking fundamental necessary and sufficient in silico mechanisms, we incremented the complexity upwards. We started with the simplest possible system and used an iterative falsification process, beginning with the first of the targeted attributes listed in Table 1. That iterative refinement method, of which parsimony is a factor, has been used successfully in addressing other simulation goals . While exploring early AXIOM specifications and the in silico conditions needed to achieve the first attribute, we mostly ignored our larger knowledge of EMT6 spheroid biology. At that stage, the analogue had one and only one goal: achieve the targeted attribute. Once that was achieved, that early SMS was valid for that one targeted attribute. As shown in Table 3, during the process of achieving Attributes 1 and 2 individual AXIOMS were qualitatively validated against their in vitro counterparts. For instance, we verified that individual CELLS did not create new CELLS more frequently than is observed in vitro. In this context, an analogue was considered valid if it exhibited attributes that matched the targeted set according to some prespecified similarity measure. We then added a new attribute, such as no. 2 in Table 1, to the targeted list. Doing so often (but not always) immediately falsified that SMS, which was the case with the addition of attribute no. 2. To revise the construct and form a new, more valid SMS with the expanded set of targeted attributes, we found it essential to introduce a volume loss mechanism and a mechanism to stabilize SMS growth and shape: we added CELL shedding and STRESS states along with AXIOMS to manage those new states. The new AXIOMS necessitated adding new parameters: moveEmptyBias and proBias. Following a period of iterative refinement, these additional mechanisms enabled the SMS to survive our attempts to falsify it with the expanded attribute list. We executed that same protocol for each of the other attributes in Table 1. Following each expansion of the attribute list, we reconsidered all AXIOMS, revising and merging parsimoniously where needed. We initially coarse-tuned parameter values and subsequently fine-tuned them. We continued that process for all attributes listed in Table 1, until the SMS was able to produce the matching growth curves in Fig. 3. That same iterative refinement method can be used to further improve SMS behaviors, and – presumably – bring SMS principles of operation into closer alignment with those of EMT6 cells. Our explicit process of iterative falsification contrasts to most prior work in executable biology, including [21, 33], which describe the completed models and predictions, but do not list the attributes targeted for reproduction or the order in which they were achieved. We believe added transparency will allow others to build on the work described here.
Some in silico parameters, such as the diameter of a CELL, mapped directly to measured observations of in vitro EMT6 quantities. These are noted by their source in literature within Table 2. One exception was the mean value of the in silico interval between when a CELL entered the NECROTIC state and when it was removed (creating an empty space). A value of five HOURS (18,000 SECONDS) was used. Doing so required three assumptions. The first was that the experimental setup used by Harris et al. to obtain these measurements did not contribute excessively to the measured apoptosis duration . The second was acceptance of the authors' assumptions about apoptosis: apoptosis begins when apoptotic morphology was observed and ended when the cell began to fragment. The final and most significant was that we could map these values to EMT6 cells undergoing necrotic cell death induced by nutrient depletion.
The actual number of simulation cycles that elapsed from when a CELL became NECROTIC and when it disappeared depended on the value of removeCounter, a pseudo-random number (PRN) drawn from a uniform distribution over the interval [0-removeDelay). RemoveCounter was decremented each cycle that a CELL was in the NECROTIC state, resulting in its removal when the value reached 0. Setting removeDelay to 36,000 SECONDS (SEC), resulted in a mean removeCounter = 18,000 SEC, which mapped directly to the reported mean duration of apoptosis .
In order to achieve the targeted attributes, it was sometimes necessary to select parameters that mapped to values that were toward the extreme end of an observed, referent range. For instance, in order to avoid excessive NUTRIENT consumption resulting in premature appearance of NECROSIS, CELLS consumed NUTRIENT at a rate of 5.0 × 10-17 MOL/CELL/SEC. Observed glucose consumption rates were between 5.5 × 10-17 and 36.0 × 10-17 mol/cell/s .
Once a subset of parameter values had been set to map to in vitro counterparts, the remaining parameter values were tuned empirically so that the similarity between SMS and in vitro attributes achieved a specified measure of similarity. Previous agent-based simulation projects demonstrated that the empirical tuning approach is an effective strategy for locating biologically relevant regions of an analogue's parameter space [5, 7, 35]. Initially, parameter values were varied extensively to discover ranges for which qualitative SMS behavior could be mapped to a corresponding biologically plausible behavior. For instance, if prolifNut (the value that must be exceeded for a CELL to remain in the PROLIFERATING state) was higher than initialVal, proliferation did not occur. Similarly, we found that moveEmptyBias had to be higher than 0.5 to prevent fissure formation and eventual SMS destabilization. Following empirical tuning to the in vitro doubling times of 18 to 24 hours , we selected a value of 2.25 for proBias and 800 SEC for prolifDelay. Each CELL had its own individual prolifCounter that specified the number of SEC that must pass before it attempted to create a new CELL. We calculated prolifCounter using the method in Walker et al. : prolifCounter = prolifDelay/2 + RG, where RG was a pseudo-random number drawn from a Gaussian distribution having mean = prolifDelay/2 and standard deviation = prolifDelay/10 . Consequently, the average prolifCounter value was roughly equal to prolifDelay. Once parameter ranges were identified that achieved the targeted measure of similarity, each parameter was adjusted in sequence over a narrow range, and the consequences for SMS properties were recorded. Values that brought simulated behaviors closer to targeted values were retained. The parameter values obtained following that protocol are identified in Table 2. Note that with the possible exception of the critical NUTRIENT levels, none of these tuned parameter values map directly to measurable in vitro counterparts, and it would be problematic to obtain such values through experimentation.
Measuring in silico and in vitro values
In vitro doubling time is the time required for an average EMT6 cluster to grow from 600 to 1,200 cells . These numbers corresponded to an SMS expanding from 8.6 × 10-3 mm2 to 1.35 × 10-2 mm2, which we used to determine SMS doubling time. Both high and low NUTRIENT VIABLE rim values were calculated by averaging the VIABLE rim width at NECROSIS onset and at the end of the simulation. Individual values of VIABLE rim width were found by counting the number of CELLS between the SMS center and the edge in three directions: from right of center, above the center, and diagonally left of the center. These three values were averaged to obtain the final width. Because the initial occurrence of NECROSIS was not necessarily stable (during early growth NECROSIS could appear and later vanish), we estimated the time of NECROSIS onset by moving backward; we specified it to be the latest time at which no NECROTIC CELLS were present. The method is similar to the one used by wet-lab researchers to estimate necrosis onset . As done with EMT6 spheroids, the maximum SMS size reached was estimated by fitting the Gompertz equation to growth data and taking the maximum size predicted by the equation.
Growth rates in vitro were based on mean measures of spheroid diameter . We converted those values into cross-sectional area in order to compare them with our in silico results. To establish a measure of SMS cross-sectional area we followed the method used in vitro, adding to the explicit phenomenological mapping between the SMS and its referent. We assumed SMS are roughly circular, which our observations demonstrated was the case. We first calculated the X and Y extents defined as follows: largest X (East-West or left-right) and Y (North-South or up-down) differences between CELLS at the edge of the SMS, ignoring detached, isolated CELLS. Those two values were averaged to obtain the measure of SMS diameter used to calculate area. This measurement adequately described the trends in SMS growth and remained quite close to the actual area occupied by all CELLS.
Analogue environment on a hexagonal grid
The width of each grid location mapped to 10 μm. Each location was either empty or held a single CELL. A second, identical sized hexagonal grid was overlaid on the first. It contained the NUTRIENT consumed by CELLS, and its value was specified using a floating-point value from 0 to 1. The NUTRIENT within the system mapped primarily to glucose, but other medium components, such as diffusible growth factors, were conflated into the referent. At this early stage, the targeted attributes selected did not include a requirement for glucose and oxygen being separate components. The single NUTRIENT factor was deemed sufficient. NUTRIENT diffusion used a discretization of the two-dimensional continuous diffusion equation du/dt = D∇2u, where D is the diffusion constant and u is the amount of diffusible material. NUTRIENT moved from high to low density areas. In a given simulation cycle, each location in NUTRIENT space calculated a new value based on the values of itself and its neighbors during the previous cycle using a method adapted for the hexagonal grid from . The new value was u new = u (1 - λ) + λ (u NE + u SE + u S + u SW + u NW + u N )/6, where λ is the discretized diffusion value and u NE , u SE , etc. are the NUTRIENT values at the neighboring locations.
CELLULAR actions, such as NUTRIENT consumption and creating new CELLS (proliferation), occur on larger timescales than diffusion, so for convenience, the NUTRIENT space underwent ten steps for each SEC in CELL space. Having multiple time scales allowed for model accuracy without wasting computation time. Physical and temporal resolutions were purposefully mapped to specific values obtained from the EMT6 system in order to improve simulation realism and help ensure that observed SMS behaviors were not artifacts of unrealistic scaling. The desired time step (Δt) for diffusion within the analogue was related to the unit distance (Δx), D, and λ, such that Δt = 3Δx2λ/8D . For the purposes of the simulation, λ = 0.28 was chosen in order to allow a Δt of 0.1 SEC when Δx is 10 μm and D is 105 μm2/SEC.
The NUTRIENT space was replenished by an algorithm that detects which empty locations lie outside the SMS and which lie inside. Empty locations outside the SMS were replenished to initialVal every ten SEC. We replenished during a simulation more frequently than was done in vitro in order to simulate the effect of stirring within EMT6 cultures: diffusion is not responsible for moving glucose and oxygen toward the EMT6 spheroid during growth, only through it . In order to calculate which depleted spaces resided inside and which were outside, the replenishment algorithm completes multiple passes from the top left corner of the grid to the bottom right and back. On each pass any empty space that is adjacent to an outside empty space is labeled as outside the SMS, sequentially replenishing deeper and deeper fissures. The algorithm used a multiple-pass approach rather than a recursive approach to avoid memory overflow errors. This algorithm was capable of replenishing fissures that extend relatively deep into the SMS.
Analogue is local based
Cells within biological systems evaluate their surroundings through direct interaction with their environment. Most information is transmitted through diffusible signals (which can travel long distances but require contact with a receptor to be recognized), cell environment, or cell-cell interactions. In order to mimic that important biological reality and to preserve a clean separation between mechanism and phenomena we added a key targeted attribute to the list: CELLS must use only local mechanisms (the third attribute in Table 1). Each CELL can query the level of NUTRIENT in its local neighborhood along with the characteristics of each neighboring location. In order to locally control SMS surface irregularities and prevent fissure formation the STRESS based movement and proliferation algorithms were developed. They only required CELLS to query their immediate neighbors. Requiring that all mechanisms must be a consequence of local events would falsify some existing individual based models of tumor spheroid growth, such as , where an artificial gradient toward the center of the simulated tumor spheroid is created that uses global knowledge.
Initial analogue state
Execution begins with a single CELL that grows unrestricted to 50 CELLS. At that stage, the AXIOMS are implemented. This simplification creates a slight, negligible shift in the growth curves. It should be noted that the initial growth rate of EMT6 spheroids was not clear because EMT6 cells were first grown in a monolayer, and then transferred to cultures dishes until cluster size reached about 400 cells. Only then were they placed in the spinner flasks .
CELLS in context
CELLS follow axiomatic operating principles that determine state change, movement, proliferation, and resource consumption. The flow chart in Additional file 1, Fig. S4 demonstrates the full range of actions CELLS can take. CELLS exist in three states: PROLIFERATING, QUIESCENT, and NECROTIC. CELL state is determined by the amount of NUTRIENT to which CELLS are exposed. PROLIFERATING and QUIESCENT CELLS consume NUTRIENT equal to proConsumeRate and quiConsumeRate, while NECROTIC CELLS do not consume NUTRIENT.
Within a simulation cycle, when a CELL finds itself adjacent to an empty space, it will move into that space, simulating random cell movement and churning. That action has the net effect of causing spaces vacated by NECROTIC CELLS to move randomly within the SMS, eventually merging with the external space. The process is illustrated in Additional file 1, Fig. S5. When a CELL is on the outer edge, it will stochastically determine if it moves into an adjacent internal space. The probability of doing so is adjusted based on the CELL'S STRESS and is biased by the value of moveEmptyBias. The higher the STRESS, the greater the likelihood the CELL will move into the empty space, smoothing the local SMS edge. An outside CELL adjacent to an interior space will move into an adjacent space if PRN < pm. The PRN is drawn from [0,1), and pm is specified by an empirically derived exponential function of moveEmptyBias and STRESS, as detailed in Additional file 1. At constant STRESS, increasing moveEmptyBias decreases the likelihood of movement, while if moveEmptyBias is constant, increased STRESS will increase the likelihood of movement.
PROLIFERATING CELLS decrement prolifCounter during each cycle. When this value drops below zero, the CELL will have an opportunity to create a new CELL. If the CELL is adjacent to empty spaces, it will select one randomly, and then be given an opportunity to place a daughter CELL (a copy) at that location. A CELL given an opportunity to proliferate will do so if PRN < pb. The PRN is drawn from [0,1), and pb is specified by an empirically derived exponential function of proBias and STRESS, as detailed in Additional file 1. At constant STRESS, the likelihood of proliferation will decrease as proBias is increased, and at constant proBias increasing STRESS causes qualified CELLS to be less likely to proliferate.
After an attempt at creating a new CELL, prolifCounter is reset regardless of whether or not the attempt was successful. A daughter CELL has the same parameter values as the parent, except for prolifCounter and removeCounter, which are set to unique random values. For simplicity, we specify that CELLS are subject to contact inhibition: only CELLS adjacent to empty space can create new CELLS. Although it is not clear to what extent contact inhibition occurs in vitro, LaRue et al.  observed that only the outer two or three cell layers proliferate at the same rate as exponentially growing cells. Proliferation beyond the SMS surface was not necessary to achieve the targeted attributes (i.e., to survive falsification with the current set of targeted attributes).
METABOLISM requires a single source of NUTRIENT
An SMS differs significantly from the individual based models of Chignola et al.  and Schaller et al. , especially in its simple representation of metabolism. The main similarities are that SMS use a diffusible NUTRIENT and CELLS DIE when insufficient NUTRIENT is available. There was no need to represent a particular type of metabolism (aerobic or anaerobic), only that CELLS consume NUTRIENT equal to proConsume or quiConsume and change state based on NUTRIENT level. We achieved the targeted attributes without being forced to add additional METABOLIC complexity. Because we achieved those attributes using a simple representation, we can achieve the same behaviors using a more complicated representation of metabolism.
Fissure formation is related to STRESS
Early SMS versions had no means to control shape, either directly or indirectly, leading to fissure formation similar to that seen in . SMS fissures were induced by diffusion-limited aggregation (DLA). DLA is a phenomenon that occurs when objects move randomly in space until they encounter and adhere to each other, forming structures with crystalline appearance . SMS fissures form as the empty spaces move about because of CELL movement, and adhere when another space is encountered. Empty spaces, even though they are not actively moving objects, are subject to DLA rules because they effectively walk randomly through the SMS. Fissures form as spaces connect to each other. The inner extreme of a fissure is the closest outside spaces to the NECROTIC core, where spaces are generated. In order to prevent fissure formation, we developed the STRESS-based proliferation algorithm. It helps prevent fissure development. CELLS at the inner extreme of any fissure will have low STRESS values, leading to preferential proliferation at that location.
Shedding of cells from the SMS surface
An occasional CELL will become isolated near the SMS surface because of normal AXIOM operation. In order to prevent their local accumulation, we implemented an algorithm that simulates shedding and the consequences of shear force caused by stirring. Any CELL that has no CELL neighbors will move randomly, selecting one of its six immediate neighbors using a uniform distribution, stopping when it encounters another CELL. Most CELLS reattach to the SMS or form small clusters. An occasional CELL will move far enough to exit the grid; it is then removed from the simulation. The number of CELLS shed currently is significantly smaller than that reported in . Shedding was not deemed a sufficiently important attribute to target at this stage. Should it be targeted, it will be straightforward to add a shedding AXIOM and then adjust parameter values to reestablish matching behaviors. As an example, a similar modification was carried out between an earlier version of the analogue  and the current one. Cells in  consumed the same quantity of NUTRIENT regardless if they were in the QUIESCENT or PROLIFERATING state, but in the current analogue the amount consumed is different.
We thank Teddy Lam, Sunwoo Park, Sean HJ Kim, Li Yan, Jonathan Tang, and Shahab Sheikh-Bahaei for useful discussion and for providing valuable feedback. Virginia Platt's assistance and support was invaluable for the completion of this research. We also thank Roberto Chignola and James Freyer for assistance in answering queries regarding tumor spheroid biology. Previous SMS versions were presented at the 2007 International Conference on Computer Applications in Industry and Engineering (CAINE 2007), San Francisco, CA, November 7–9, 2007 and the 2007 International Conference on Complexity in Acute Illness (ICCAI), Long Beach, CA, October 5–7, 2007. This work was supported in part by a Fellowship to JAE provided by the CDH Research Foundation.
- Webster's Third New International Dictionary, Unabridged. 2002, Springfield, MA: Merriam-Webster
- Jiang Y, Pjesivac-Grbovic J, Cantrell C, Freyer J: A multiscale model for avascular tumor growth. Biophys J. 2005, 89: 3884-3894. 10.1529/biophysj.105.060640PubMed CentralView ArticlePubMedGoogle Scholar
- Roose T, Chapman SJ, Maini PK: Mathematical models of avascular cancer. SIAM Review. 2007, 49: 179-208. 10.1137/S0036144504446291.View ArticleGoogle Scholar
- Araujo RP, McElwain DLS: A history of the study of solid tumour growth: the contribution of mathematical modelling. Bull Math Biol. 2004, 66: 1039-1091. 10.1016/j.bulm.2003.11.002View ArticlePubMedGoogle Scholar
- Grant MR, Mostov KE, Tlsty TD, Hunt CA: Simulating properties of in vitro epithelial cell morphogenesis. PLoS Comput Biol. 2006, 2 (10): e129- 10.1371/journal.pcbi.0020129PubMed CentralView ArticlePubMedGoogle Scholar
- Schaller G, Meyer-Hermann M: Multicellular tumor spheroid in an off-lattice Voronoi-Delaunay cell model. Phys Rev E Stat Nonlin Soft Matter Phys. 2005, 71: 051910-View ArticlePubMedGoogle Scholar
- Tang J, Ley K, Hunt CA: Dynamics of in silico leukocyte rolling, activation, and adhesion. BMC Systems Biology. 2007, 1: 14- 10.1186/1752-0509-1-14PubMed CentralView ArticlePubMedGoogle Scholar
- Burton AC: Rate of growth of solid tumours as a problem of diffusion. Growth. 1966, 30: 157-176.PubMedGoogle Scholar
- Casciari JJ, Sotirchos SV, Sutherland RM: Mathematical modelling of microenvironment and growth in EMT6/Ro multicellular tumour spheroids. Cell Prolif. 1992, 25: 1-22. 10.1111/j.1365-2184.1992.tb01433.xView ArticlePubMedGoogle Scholar
- Chignola R, Milotti E: A phenomenological approach to the simulation of metabolism and proliferation dynamics of large tumour cell populations. Phys Biol. 2005, 2: 8-22. 10.1088/1478-3967/2/1/002View ArticlePubMedGoogle Scholar
- Kim Y, Stolarska M, Othmer HG: A hybrid model for tumor spheroid growth in vitro I: theoretical development and early results. Math Models Methods Appl Sci. 2007, 17: 1773-1798. 10.1142/S0218202507002479.View ArticleGoogle Scholar
- Galle J, Loeffler M, Drasdo D: Modeling the effect of deregulated proliferation and apoptosis on the growth dynamics of epithelial cell populations in vitro. Biophys J. 2005, 88: 62-75. 10.1529/biophysj.104.041459PubMed CentralView ArticlePubMedGoogle Scholar
- Freyer JP, Sutherland RM: Regulation of growth saturation and development of necrosis in EMT6/Ro multicellular spheroids by the glucose and oxygen supply. Cancer Res. 1986, 46: 3504-3512.PubMedGoogle Scholar
- LaRue KEA, Khalil M, Freyer JP: Microenvironmental regulation of proliferation in multicellular spheroids is mediated through differential expression of cyclin-dependent kinase inhibitors. Cancer Res. 2004, 64: 1621-1631. 10.1158/0008-5472.CAN-2902-2View ArticlePubMedGoogle Scholar
- Gompertz B: On the Nature of the Function Expressive of the Law of Human Mortality, and on a New Mode of Determining the Value of Life Contingencies. Philosophical Transactions of the Royal Society of London. 1825, 115: 513-583. 10.1098/rstl.1825.0026.View ArticleGoogle Scholar
- Chignola R, Schenetti A, Andrighetto G, Chiesa E, Foroni R, Sartoris S, Tridente G, Liberati D: Forecasting the growth of multicell tumour spheroids: implications for the dynamic growth of solid tumours. Cell Prolif. 2000, 33: 219-229. 10.1046/j.1365-2184.2000.00174.xView ArticlePubMedGoogle Scholar
- Casciari JJ, Sotirchos SV, Sutherland RM: Variations in tumor cell growth rates and metabolism with oxygen concentration, glucose concentration, and extracellular pH. J Cell Physiol. 1992, 151: V386-394. 10.1002/jcp.1041510220.View ArticleGoogle Scholar
- Dorie MJ, Kallman RF, Rapachhietta DF, Van Antwerp D, Huang YR: Migration and internalization of cells and polystyrene microspheres in tumor cell spheroids. Exp Cell Res. 1982, 141: 201-209. 10.1016/0014-4827(82)90082-9View ArticlePubMedGoogle Scholar
- Freyer JP: Role of necrosis in regulating the growth saturation of multicellular spheroids. Cancer Res. 1988, 48: 2432-2439.PubMedGoogle Scholar
- Landry J, Freyer JP, Sutherland RM: Shedding of mitotic cells from the surface of multicell spheroids during growth. J Cell Physiol. 1981, 106: 23-32. 10.1002/jcp.1041060104View ArticlePubMedGoogle Scholar
- Longo D, Peirce SM, Skalak TC, Davidson L, Marsden M, Dzamba B, DeSimone DW: Multicellular computer simulation of morphogenesis: blastocoel roof thinning and matrix assembly in Xenopus laevis. Dev Biol. 2004, 271: 210-222. 10.1016/j.ydbio.2004.03.021View ArticlePubMedGoogle Scholar
- Delsanto PP, Condat CA, Pugno N, Gliozzi AS, Griffa M: A multilevel approach to cancer growth modeling. J Theor Biol. 2008, 250: 16-24. 10.1016/j.jtbi.2007.09.023View ArticlePubMedGoogle Scholar
- Anderson ARA, Weaver AM, Cummings PT, Quaranta V: Tumor morphology and phenotypic evolution driven by selective pressure from the microenvironment. Cell. 2006, 127: 905-915. 10.1016/j.cell.2006.09.042View ArticlePubMedGoogle Scholar
- Sutherland RM, Durand RE: Radiation response of multicell spheroids – an in vitro tumour model. Curr Top Radiat Res Q. 1976, 11: 87-139.PubMedGoogle Scholar
- Freyer JP, Sutherland RM: Selective dissociation and characterization of cells from different regions of multicell tumor spheroids. Cancer Res. 1980, 40: 3956-3965.PubMedGoogle Scholar
- Freyer JP, Sutherland RM: A reduction in the in situ rates of oxygen and glucose consumption of cells in EMT6/Ro spheroids during growth. J Cell Physiol. 1985, 124: 516-524. 10.1002/jcp.1041240323View ArticlePubMedGoogle Scholar
- Freyer JP, Sutherland RM: Proliferative and clonogenic heterogeneity of cells from EMT6/Ro multicellular spheroids induced by the glucose and oxygen supply. Cancer Res. 1986, 46: 3513-3520.PubMedGoogle Scholar
- Leaf C: Why we're losing the war on cancer (and how to win it). Fortune. 2004, 149: 76-82. 84–6, 88PubMedGoogle Scholar
- Wartenberg M, Acker H: Quantitative recording of vitality patterns in living multicellular spheroids by confocal microscopy. Micron. 1995, 26: 395-404. 10.1016/0968-4328(95)00009-7View ArticlePubMedGoogle Scholar
- Luke S, Balan GC, Panait LA, Paus S: MASON: a Java multiagent simulation library. Proceedings of Agent. Edited by: Sallach DL, Macal C. 2003, 49-64. http://cs.gmu.edu/~eclab/projects/mason/papers/Agent2003Presentation.pdfGoogle Scholar
- Fisher J, Henzinger TA: Executable cell biology. Nat Biotech. 2007, 25: 1239-1249. 10.1038/nbt1356.View ArticleGoogle Scholar
- Hunt CA, Ropella GEP, Park S, Engelberg JA: Dichotomies between computational and mathematical models. Nat Biotech. 2008, 26: 737-738. 10.1038/nbt0708-737.View ArticleGoogle Scholar
- Christley S, Alber MS, Newman SA: Patterns of mesenchymal condensation in a multiscale, discrete stochastic model. PLoS Comp Biol. 2007, 3 (4): e76-10.1371/journal.pcbi.0030076.View ArticleGoogle Scholar
- Yan L, Ropella GEP, Park S, Roberts MS, Hunt CA: Modeling and simulation of hepatic drug disposition using a physiologically based, multi-agent in silico liver. Pharm Res. 2008, 25: 1023-1036. 10.1007/s11095-007-9494-yView ArticlePubMedGoogle Scholar
- Harris LK, Keogh RJ, Wareing M, Baker PN, Cartwright JE, Aplin JD, Whitley GSJ: Invasive trophoblasts stimulate vascular smooth muscle cell apoptosis by a fas ligand-dependent mechanism. Am J Pathol. 2006, 169: 1863-1874. 10.2353/ajpath.2006.060265PubMed CentralView ArticlePubMedGoogle Scholar
- Walker DC, Southgate J, Hill G, Holcombe M, Hose DR, Wood SM, Neil SM, Smallwood RH: The epitheliome: agent-based modelling of the social behaviour of cells. Biosystems. 2004, 76: 89-100. 10.1016/j.biosystems.2004.05.025View ArticlePubMedGoogle Scholar
- Ermentrout GB, Edelstein-Keshet L: Cellular automata approaches to biological modeling. J Theor Biol. 1993, 160: 97-133. 10.1006/jtbi.1993.1007View ArticlePubMedGoogle Scholar
- Dormann S, Deutsch A: Modeling of self-organized avascular tumor growth with a hybrid cellular automaton. In Silico Biol. 2002, 2: 393-406.PubMedGoogle Scholar
- Rejniak KA: A single-cell approach in modeling the dynamics of tumor microregions. Math Biosc Eng. 2005, 2: 643-655.View ArticleGoogle Scholar
- Witten TA, Sander LM: Diffusion-limited aggregation, a kinetic critical phenomenon. Phys Rev Lett. 1981, 47: 1400-1403. 10.1103/PhysRevLett.47.1400.View ArticleGoogle Scholar
- Engelberg JA, Hunt CA: A local mechanism generates saturation in an in silico model of in vitro multicellular tumor spheroid growth. J Crit Care. 2007, 22: 331-10.1016/j.jcrc.2007.10.002.View ArticleGoogle Scholar