From pathway to population – a multiscale model of juxtacrine EGFR-MAPK signalling
© Walker et al. 2008
Received: 03 June 2008
Accepted: 26 November 2008
Published: 26 November 2008
Skip to main content
© Walker et al. 2008
Received: 03 June 2008
Accepted: 26 November 2008
Published: 26 November 2008
Most mathematical models of biochemical pathways consider either signalling events that take place within a single cell in isolation, or an 'average' cell which is considered to be representative of a cell population. Likewise, experimental measurements are often averaged over populations consisting of hundreds of thousands of cells. This approach ignores the fact that even within a genetically-homogeneous population, local conditions may influence cell signalling and result in phenotypic heterogeneity.
We have developed a multi-scale computational model that accounts for emergent heterogeneity arising from the influences of intercellular signalling on individual cells within a population. Our approach was to develop an ODE model of juxtacrine EGFR-ligand activation of the MAPK intracellular pathway and to couple this to an agent-based representation of individual cells in an expanding epithelial cell culture population. This multi-scale, multi-paradigm approach has enabled us to simulate Extracellular signal-regulated kinase (Erk) activation in a population of cells and to examine the consequences of interpretation at a single cell or population-based level using virtual assays.
A model consisting of a single pair of interacting agents predicted very different Erk activation (phosphorylation) profiles, depending on the formation rate and stability of intercellular contacts, with the slow formation of stable contacts resulting in low but sustained activation of Erk, and transient contacts resulting in a transient Erk signal. Extension of this model to a population consisting of hundreds to thousands of interacting virtual cells revealed that the activated Erk profile measured across the entire cell population was very different and may appear to contradict individual cell findings, reflecting heterogeneity in population density across the culture. This prediction was supported by immunolabelling of an epithelial cell population grown in vitro, which confirmed heterogeneity of Erk activation.
These results illustrate that mean experimental data obtained from analysing entire cell populations is an oversimplification, and should not be extrapolated to deduce the signal:response paradigm of individual cells. This multi-scale, multi-paradigm approach to biological simulation provides an important conceptual tool in addressing how information may be integrated over multiple scales to predict the behaviour of a biological system.
Experimentally, cell populations are usually considered to be homogeneous. Techniques such as Western Blotting inherently make the assumption that variations between individual cells are not significant. Furthermore, most computational models of signalling pathways relate to events in a 'typical' cell. Whilst adequate for most purposes, these models do not address the issue of heterogeneity – particularly in the context of how variations in the micro-environment influence gene and protein expression of individual cells.
This issue of local heterogeneity in microenvironment, and the potential influence on cell signalling and fate has yet to be studied extensively. However, recent publications have suggested that such factors may play a critical role in determining cell fate, with implications for clinically important phenomena such as stem cell lineage fate in MLL-AF9 leukaemia . The architecture of intercellular contacts within beta cell islets has been shown to result in heterogenous production of insulin . In this paper, we describe how we have used multi-scale computational simulations, supported by experimental data, to explore how local differences in intercellular contact may influence intracellular signalling, with implications for individual cell response and population heterogeneity.
Cell behaviour can be modelled by representing individual cells as computational entities or software agents. Unlike conventional modelling of intracellular signalling pathways, agent modelling allows the response of signalling to be interpreted at the cellular level. For example, we can incorporate the hypothetical rule that a particular threshold amount or concentration of a product is required for an individual cell to progress through the cell cycle. This method allows the study of emergent behaviour of a system (e.g. tissue growth or wound healing) as the outcome of the interaction of the individual components (the cells). We have previously developed Epitheliome, an agent-based model of epithelial cell populations. Software agents represent individual cells and iteratively change their state (e.g. cell cycle stage, location, shape) according to a number of pre-programmed rules representing biological behaviour, such as proliferation, intercellular adhesion, migration and apoptosis. We have used this simple rule-based model to simulate monolayer growth  and wound healing  in normal human urothelial cell cultures and have adapted it to explore stratification and differentiation in normal and transformed keratinocytes . More recently, we have combined this model with a mathematical representation of autocrine ligand release, diffusion and binding . The biological basis for this work has been experimentally-generated observations in cultures of normal human urothelial (NHU) cells where, in the absence of exogenous growth factors, proliferation is driven through autocrine production of EGFR-binding factors . In the NHU cell culture system, inhibition either of EGFR, or of the downstream MAPK/Erk pathway components, results in dephosphorylation of Erk and an inhibition of proliferation that is reversible upon removal of the inhibitor . Thus, although it is clear that multiple, interacting signalling pathways influence cell cycle progression and proliferation ([8–10]), there is compelling evidence for a causal link between EGFR-Erk signalling pathway and proliferation, at least in these cells. Autocrine growth mechanisms may operate through the release of soluble ligands, or through juxtacrine signalling mediated via cell:cell interactions and we have used a modelling approach to examine the implications of these modes of communication. Our approach of using an individual-based paradigm naturally lends itself to examining heterogeneity, either by varying the internal parameters (memory state) of each agent, or by observing the emergent behaviour of agents in differing micro-environments.
The initial realisation of the Epitheliome model can be considered as a phenomenological representation of actual cell behaviour, as the rules governing state transitions are based purely on observations made experimentally, many of which were extracted from the scientific literature. For instance, agents progress around the cell cycle consisting of phases representing G1, S, G2 and M, traversing a single checkpoint in G1. As described in , transition through this checkpoint is dependent on cell:cell bonding and cell spreading, with the underpinning rule set mimicking the phenomenon of contact inhibition of growth. Results produced by this model successfully predicted the reduction of growth rate observed experimentally when epithelial cell cultures were switched to culture conditions that promoted formation of intercellular adherens junctions, through the calcium-dependent homotypic binding of E-cadherin . However, a consistent enhancement of the in vitro growth rate in low density populations grown in physiological calcium conditions was not predicted computationally, and this divergence of in vitro and in virtuo systems led us to investigate potential mechanisms for enhanced growth mediated via intercellular contact, and in particular, juxtacrine signalling via the Epidermal Growth Factor Receptor (EGFR).
EGFR is known to be a critical mediator of growth control in many cell types, including urothelium  and there is growing experimental evidence that EGFR ligands are biologically active prior to cleavage from the membrane . In order to investigate a possible role for contact-mediated juxtacrine signalling in epithelial cell cultures, we have utilised a mathematical ordinary differential equation (ODE) model of membrane-bound EGFR-ligand interaction, receptor dimerisation/activation and downstream signalling via the Ras-Raf-MAPK pathway. The end point of this model pathway is the activation (diphosphorylation) of the cytoplasmic protein, Erk. It is known that activated Erk (Erk-PP) translocates to the nucleus where it can activate transcription of key cell cycle regulators, such as cyclin D1. The temporal characteristics of Erk activation have previously been shown to be closely related to cell fate, with sustained, moderate Erk activation correlated with cell cycle progression [12, 13].
The primary criteria for selecting a suitable model of intracellular signalling for integration into our existing modelling framework were: 1) ability of the model to simulate the key features of the output (temporal changes in Erk-PP expression) with respect to input (receptor activation), 2) minimum complexity (to introduce parsimony by restricting the number of ODEs required for future multi-agent and hence multi-contact simulations), and 3) inclusion of key stages from the entire pathway, as opposed to a limited subset, to allow future mapping of signal to cellular response.
On the basis of these criteria, the model chosen as a starting point was that of , which was originally used to examine the temporal pattern of Erk activation in response to soluble EGF or NGF binding to surface receptors. This model was selected due to its relative simplicity compared to other more detailed descriptions, although this does not impact on its ability to predict key features of the pathway, notably the high amplitude transient peak in activated Erk, followed by a lower, but sustained activation for periods of several hours. The receptor-ligand dynamics were slightly simplified and the equations used for this part of the model were those implemented in . A schematic of the kinetic model is shown in figure 1a.
EGFR signalling and trafficking.
Intercellular contact dynamics.
σ = 3.5/60
The first three terms on the right hand side of the equations (1) and (2) in table 1 represent changes in the numbers of species arising from receptor-ligand interactions. The fourth term represents receptor internalisation or ligand cleavage, the fifth represents receptor/ligand synthesis and the last terms describe recruitment of new receptors or ligands into expanding contact areas. The constant k 1 is the receptor-ligand association constant; k -1 the dissociation constant; k int the rate of internalisation of unoccupied receptors, k Rsyn and k Lsyn , the rates of de novo receptor and ligand synthesis respectively per cell; SA denotes the total surface area of the cell; σ is the rate constant for E-cadherin mediated contact growth; R 0, L 0 and S A0 are the total receptor or ligand numbers and surface area on the cell not associated with intercellular contacts at any given time (given by equations 7–10 in table 2); and Ω AB is a flag indicating active contact growth, where Ω AB = 1 if the contact is defined to be expanding in size, and = 0 otherwise. Diffusion of EGFR and ligands in and out of the contact area is ignored. This is justified by the observation that typical diffusivities are extremely small (of the order of 10-10 cm2/second), and that movement of surface EGF receptors is restricted by the cytoskeletal network .
The number of occupied, dimerised and phosphorylated (active) receptors associated with contact AB on cell A will then be given by equation (5) in table 1. Similar equations will exist for EGF receptor-associated species on cell B associated with ligands on cell A.
Intracellular pathway model. Rate equations relating to intracellular EGFR-MAPK pathway.
k5 = 12
K5 = 6 × 103
V6= 3.0 × 105
K6 = 6 × 103
K7 = 2 × 10-3
k-7 = 3.81
K8 = 1.6;
K8 = 6 × 105
V9 = 75;
K9 = 2.0 × 104
k10 = 1.63 × 10-2
k10 = 10
k11 = 15
k12 = 5.0 × 10-3
k-12 = 60
k13 = 7.2 × 102
k14 = 1.2 × 10-3
k-14 = 3.0
k15 = 27
V16 = 9.7 × 104
K16 = 6 × 103
k17 = 50;
K17 = 9 × 103
V18 = 9.2 × 105
K18 = 6 × 105
k19 = 50;
K19 = 9 × 103
V20 = 9.2 × 105
K20 = 6 × 105
k21 = 8.3;
K21 = 9 × 104
V22 = 2.0 × 105
K22 = 6 × 105
k23 = 8.3;
K23 = 9 × 104
V24 = 4.0 × 105;
K24 = 6 × 105
The intracellular pathway associated with each cell is thus represented by 19 ODEs (equations 11–29, table 3), with individual steps represented by standard mass action kinetics, or, where appropriate, Michaelis-Menten kinetics . This model was solved using the Mathworks Matlab ode23tb solver.
The agent (rule-based) and signalling (ODE-based) models are run alternately. Each agent model iteration represents a total period of 30 minutes and the total solution period of the following call to the ODE model represents the same 30 minute period. The data passed between the models on a per contact and per agent basis is summarised in figure 3. It is assumed that agent sizes and positions are static during the solution of the signalling model (though existing contacts can grow in length). The creation or destruction of contacts is determined entirely by the agent model, so the minimum length of time that a transient (non E-cadherin stabilised) contact exists is 30 minutes. This is reasonable, given that a detailed study of high-resolution time-lapse images of urothelial cells grown in low calcium conditions (0.09 mM) suggested that most transient intercellular contacts persist between 10 and 60 minutes (e.g. additional files 1 and 2). At the end of every simulation, the updated profile of Erk-PP (equation 29) is returned to the memory of each agent. Numbers of all intracellular molecular species (11–29) are stored for every agent and surface species (1–5) for every contact, in order to provide the initial conditions for the next set of calculations. The identities of cells sharing any contact are also stored. In the event that a pre-existing contact is registered as broken following the subsequent agent iteration (a frequent occurrence for non-E-cadherin-mediated cell contacts), the receptor and ligand terms in equations (1) and (2), and equivalent equations for cell B, are set to zero and the ODE solution deals with residual intracellular signalling only.
The model described above was used to study differences in the activated EGFR-mediated conversion of Erk into its active diphosphorylated form (active Erk or Erk-PP) associated with the differing nature of contacts seen in low (0.09 mM) and near physiological (2 mM) calcium environments. This is considered to be the primary 'output' of the model, due to its role in mediating cell behaviour [12, 13]. Simulations were initially carried out for the simplified case of a single pair of agents sharing a transient or E-cadherin-mediated contact. This concept was then extended to a much larger model consisting of hundreds to thousands of agents forming contacts on a stochastic basis.
The basis for this component of the computational model is the prototype Epitheliome model, running on the Mathworks Matlab platform, as described in detail in [3, 4]. Briefly, individual cells are represented by software agents (instances of class objects in object-oriented Matlab). Each agent contains memory parameters representing its current state (e.g. location, size and cell cycle position). The simulation is run iteratively with each agent updating its memory parameters ('state') according to a number of pre-programmed rules representing biological behaviour such as proliferation, intercellular adhesion, migration and apoptosis. Parameters are scaled such that each model iteration represents 30 minutes in real time. A numerical algorithm is employed after each agent update to shift the cells to minimise edge overlap arising due to cell growth, division and migration. We refer the reader to our earlier publication for a comprehensive description of basic model operation and the details of this numerical algorithm .
In order to focus on modelling cellular interaction, and in particular, biochemical signalling processes that might be influenced by direct cellular contact, it was necessary to increase the complexity of rule sets governing migratory and adhesive behaviour. These changes are described below:
The cadherins are a family of transmembrane proteins that mediate intercellular adhesion via interactions between homotypic proteins on opposing cell surfaces. Epithelial tissues express the cadherin subgroup, E-cadherin, which is critical for forming initial adherens contacts between cells in developing tissues or cell cultures, allowing more established structures such as tight junctions and desmosomes to develop. In order to assume a functional conformation, E-cadherin requires the presence of extracellular calcium ions. It has been demonstrated that intercellular contacts form in cultures of normal human urothelial (NHU) cells at a critical calcium concentration ≥ 1 mM . A sigmoidal relationship between bonding activity and exogenous calcium, with the inflection point located at 1 mM was measured in doublets of Chinese hamster ovary cells .
where P = amount of functional E-cadherin, A = lower asymptote (= 0.1), C = upper asymptote (= 0.57), M = point of max growth (= 1 mM), B = growth rate (= 1), T determines relationship between max growth and asymptotes (= 1) and X is the exogenous calcium concentration in mM. The relationship between functional E-cadherin amount, P and exogenous calcium thus closely mirrors the relationship described by  between binding activity and calcium concentration.
The pair of cells will form an initial E-cadherin mediated contact if a pseudo-random number is less than BP.
Previously, we have modelled intercellular bonds simply as being either present or absent. However, in order to explore juxtacrine signalling, it was necessary to incorporate the concept of a finite bond length. On forming an initial contact, E-cadherin molecules become aggregated, leading to the 'zippering' of adjacent cell membranes over a number of hours . We devised a rule dictating the increase in bond length with time from initial contact formation, based on a quantitative interpretation of these results:
for (0.5 < t < 12)
BL = 2 + (σ *t)
BL = 15;
where BL = bond length and t = time in hours and σ = 3.5 μm/hour. As this linear expansion of contacts occurs on a similar timescale as the intracellular signalling dynamics, this 'rule' is incorporated into the signalling model (equation 6, table 2).
Close inspection of time-lapse microscopy images of NHU cells cultured in low calcium medium (0.09 mM) revealed that although significantly fewer stable contacts are formed than in physiological calcium, cells tended to remain in contact with one another for periods of up to one hour, before migrating away to form new contacts (see additional files 1 and 2). During this period of transient contact, it is feasible that juxtacrine signals may be exchanged between cells and therefore, we have chosen to include the concept of transient cell contact in our agent model. These could represent weak interactions mediated by inter-membrane surface tension, as implied by immunofluorescence microscopy (see , and additional files 1 and 2). We have included the rule:
for cells with edge separation, sep, less than 0.1 μm
if sep < 0
clen = max_r;
clen = (0.1-sep)*max_r/0.1;
where clen is the length of the non-specific contact length and max_r is 1.5 times the larger of the two cell radii The form of this relationship is somewhat arbitrary, but the range of contact lengths produced gives a reasonable agreement with that observed in microscopy images of growing NHU cell cultures (, additional file 1). Note that cell agents in low calcium are also permitted to form E-Cadherin mediated contacts, but according to the rules described above, these occur with a much lower frequency than in higher calcium environments.
In our earlier agent models (e.g. [3, 6]), only cells with zero intercellular contacts were able to migrate. Introduction of the new bond probability rules and the concept of transient contacts described above resulted in very limited migration even in low calcium conditions; this situation did not mirror that observed in real NHU cell cultures. Hence, new rules were introduced that determined the ability of an agent to migrate to be the outcome of the difference between the active migration force (reported to be approximately 20 nN ) and intercellular adhesion force resulting from stable and transient cell contacts (reported in the former case to be approximately 100 nN after 30 minutes, and 200 nN after one hour of intercellular contact for normal E-cadherin expression levels, and linearly related to the square of the amount of E-cadherin present at the cell surface ). In the case of transient contacts, surface tension between non-adherent cells was calculated to be approximately 0.32 nN per micrometre of boundary in contact . This leads to the following formulation for total adhesion force acting on a cell:AF = (0.32 * clen) + (P 1 * P 2(100 * b1 + 200 * b2))
where clen is the length of cell boundary involved in non specific contacts, b1 and b2 are the number of E-cadherin mediated bonds that have existed for less and more than one hour respectively, and P1 and P2 are the levels of functional E-cadherin expression (between 0 and 1).
The model rule implemented was that a total intercellular adhesion force less than 20 nN would result in migration of a particular cell, but for larger values the cell was not permitted to migrate. For normal levels of functional E-cadherin expression, this effectively means that cells with one or more E-cadherin contacts were prevented from migration (as was the case in the earlier versions of our model), whereas those with contacts defined as transient were permitted to migrate only if the total boundary length involved is less than 62.5 μm. The maximum migration velocity was estimated from time-lapse images of NHU cells to be in the region of 80 μm/hour.
A potential flaw of many computational models is hyper-sensitivity to one or more parameters. This can be particularly problematic when the parameters are difficult to measure experimentally and are thus estimated or assumed. This is a critical issue in multi-scale models such as that described here, where sensitivity to parameters at one model scale could potentially propagate through into higher scales, often with unforeseen effects.
In order to assess the sensitivity of the results described above to the model parameters, a sensitivity analysis was carried out as follows. Each of the 50 parameters in the signalling pathway model was varied systematically by first dividing, then multiplying the value of each parameter by a factor of two. The reason for selecting a relatively large perturbation factor was that the kinetic parameters for the upstream pathway model were derived from radio-labelling studies utilising the soluble form of an EGFR binding ligand (e.g. ) and parameters for membrane-bound receptor-ligand interactions might be expected to vary quite significantly.
Extracellular calcium is known to be a key mediator of intercellular contacts, with concentrations greater than 1 mM resulting in stable adherens junctions mediated by the homotypic interaction of E-Cadherin molecules expressed on the membranes of adjacent cells . Additionally, cells may form short-term or transient contacts that are independent of the exogenous calcium concentration (for example, see additional files 1 and 2).
Although the results from the simple two cell model suggest that Erk activation associated with an individual transient contact is relatively short lived (figure 4), the summation over relatively large number of cells 'cultured' in low calcium medium, continuously making and breaking contacts, leads to continuously elevated Erk-PP. By contrast, in physiological calcium, stable contacts are formed quickly. The consequence of this is that averaged levels of Erk activation measured at fixed time points in cell populations exposed to the two different environments both appear very similar (figure 9), even though the activation levels in individual cells can vary significantly. These results suggest that in populations where the cell density is relatively high, but growth is not contact-inhibited (ie pre-confluent cultures with several hundred cells per mm2), the activated Erk signal averaged over the population would be higher in low calcium (e.g. figure 9b), than in physiological calcium environments.
In simulations of very low cell density (figure 9a) intercellular contacts are scarce, and Erk activation is always higher in physiological calcium where any intercellular contacts that occur are maintained for longer durations. For a starting density of 200 cells/mm2, the Erk-PP signal in the low calcium simulations exceeds that in the physiological calcium simulations after 12 hours (figure 9b), and this relative elevation persists until at least 24 hours. For seeding densities of 500 and 700 cells/mm2 (figure 9c–d), the population levels quickly approach confluence. In these conditions, space constraints limit cell migration and the opportunity for cells to form transient contacts, and the Erk-PP profiles are similar in both low and physiological calcium environments.
A qualitative, second level sensitivity analysis was conducted for a starting cell density of 200 cells/mm2. Five kinetic parameters pertaining to ligand/receptor interactions, five downstream parameters with relatively high sensitivity coefficients and the E-cadherin contact growth rate parameter were systematically varied by a factor of two. The cell culture simulation was run for a period representing 25 hours. Although the amplitudes of the virtual Erk activation level were altered, the prediction of comparable or higher activation levels in the low calcium, relative to the physiological calcium simulations continued to hold true in every case (results not shown).
Ideally, we would validate our single cell pathway model by tracking Erk activation in cells making different types of contacts under controlled conditions. Unfortunately, this is not possible due to difficulties in visualising protein activation levels and intercellular contact formation in real time. However, our modelling results predict that an otherwise homogeneous population of non-confluent cells sampled from a small spatial region would at any moment in time have a unique pattern of intercellular contacts and hence would be expected to show heterogeneity in Erk activation if examined by immunofluorescence microscopy.
Figure 6a–b shows a snap-shot of the distribution of Erk-PP in urothelial cells maintained in low and physiological calcium visualised by immunofluorescence microscopy. All cells show some degree of Erk activation (indicated by the presence of nuclear Erk), which can be attributed to pathway activation by soluble EGF present in the culture medium. However, in both images there is clear evidence of heterogeneity, with the cells indicated showing low, moderate and high levels of Erk activation, respectively, which suggests that there is a secondary mode of stimulation in these cells.
Our model of juxtacrine signalling between individual cells has yielded two interesting results. At the level of a pair of cells forming a single intercellular contact, the model suggests that the temporal characteristics of activated Erk are highly dependent on the rate of formation and stability of the contact. Specifically, our results suggest that transient contacts result in only transient Erk activation, irrespective of the spatial extent of the contact, whereas more slowly formed, stable contacts stabilised by the interaction of E-Cadherin result in a small but sustained Erk activation, due to the engagement of new receptors and ligands. Secondly, when incorporated into a multi-agent simulation and used to derive an averaged Erk-PP value for the population (analogous to a Western blotting experiment), seemingly contradictory results can be obtained. Specifically, for low density populations (<400 cells/cm2), the averaged Erk activation level in low calcium environments (where cells are migratory and most contacts are transient) is equal, or even higher than in physiological calcium environments (where most contacts develop slowly, then remain stable).
Our model predictions relating to Erk activation at the single cell level has potential implications for our understanding of how proliferation through juxtacrine mechanisms is regulated within cell populations. The temporal characteristics of Erk activation – the Erk-PP "signature" – is recognised as a critical factor in determining cell fate. Thus, an initial peak in Erk-PP followed by sustained activation for several hours has been shown to be associated with cell cycle progression, whereas a transient peak followed by a rapid decline is insufficient to induce cycle progression ([12, 13]). This suggests that by providing a sustained basal level of Erk activation, juxtacrine signalling via EGFR could provide a mechanism for the increased growth rates that we have previously observed in NHU cells cultured in physiological versus low calcium environments . A positive relationship between cell contact and proliferation has also been observed in endothelial and smooth muscle cells  and in kidney epithelial cells .
The negative association of cell-cell contact and proliferation (contact inhibition of growth) is well accepted. In our agent-based model, this is currently represented by an arbitrary rule that a cell exits the mitotic cycle when it has ≥ n bonds with neighbouring cells. At this stage we do not know which other intercellular mechanisms are involved in positive and negative growth modulation, but one candidate is β-catenin, which is a nuclear transcription factor inactivated by E-cadherin sequestration. Future computational and experimental exploration of intracellular processes involving β-catenin that builds on the work of others  should indicate whether this is a plausible explanation.
We have not yet fully developed rules for the agent level simulation to relate the proliferation of cells (or more specifically, their decision to pass the G1-S checkpoint) to the Erk-PP profile. This is a natural extension of the model and would allow us to explore how intercellular interactions impact on intracellular molecular state and ultimately on normal and abnormal tissue growth.
By necessity, the model presented here is a simplification of the real intracellular pathway, which may include up to 19,000 distinct molecular species (). We have so far neglected possible membrane diffusion of receptors and ligands, potential targeting of newly synthesised EGFR and ligand to, or clustering at, stable E-cadherin-mediated contacts. We have also ignored the presence of soluble autocrine or exogenous ligands in the medium, which might provide some level of pathway activation independently of any juxtacrine effect (as demonstrated by immunofluorescence in figure 6, where exogenous EGF was present in the growth medium) and other interacting intracellular cascades that may also be downstream of the EGFR (e.g. the P13K-Akt pathway ). It is possible that membrane-bound and soluble forms of a ligand may mediate differential effects on EGFR activation. However, sensitivity testing suggests that our basic finding that sustained Erk activation is associated with gradually-formed, stable contacts is not affected by a perturbation in kinetic parameters.
The multi-agent simulations have yielded results that have potential implications for the interpretation of experimental data. The frequency of intercellular contact in a cell population is related to cell density – the greater the number of cells in a given surface area (or volume in 3D culture) the greater the likelihood of contact formation. We have observed that even in low calcium environments, cells can form E-cadherin-mediated adherens junctions in medium to high density cultures/simulations, although these are substantially less common than in physiological calcium (figures 7, 8,  and additional files 3 and 4). Simulated cell cultures with initial cell densities of up to the maximum simulated value of 700 cells/mm2 are sub-confluent and remain within the exponential phase of population growth (a 1 mm × 1 mm model reaches visual confluence at approximately 1000 cells). This observation, along with model predictions of the relative frequency and perimeter length associated with E-cadherin mediated cell contact in the different environments (as shown in figure 7), indicates that at high density, cells in low calcium are still able to migrate and make contacts that are predominantly transient in nature (as shown in additional file 2). This suggests that the elevated population-based Erk-PP signal in low calcium simulations (figure 9) represents the summation of multiple transient signals over relatively large cell numbers, rather than an alteration of the nature of intercellular contacts at increased cell density.
These results provide an important insight into the potential pitfalls of interpreting experimental data that is averaged over relatively large cell numbers. Practical limitations mean that by necessity, data is often collected from cell populations and used to draw inferences on individual cell behaviour or intracellular state that cannot be easily observed in the laboratory. A common of example of the latter in a cell biology context is the interpretation of Western blotting data obtained from cultured cell populations. In this case, amounts of a particular protein (active, inactive or total) averaged across a large cell population are measured and used to draw conclusions relating to expression levels in individual cells in the population. We have presented the simulated levels of activated Erk across our entire agent population in the form of a 'virtual' Western blot (figure 9). These model predictions highlight the fact that experimental data obtained from populations of cells should be viewed with caution, particularly in terms of inferring behaviour on the level of the individual cell. The methodology of extracting and analysing our simulated data is directly analogous to the process that an experimental cell biologist would routinely use in the laboratory and highlights the versatility of the agent-based paradigm and its natural relationship to the experimental world.
Our multi-scale model has demonstrated how a) changes in the extrinsic environment affect intercellular contact formation and impact upon intracellular molecular state and b) how the combination of many instances of these states can impact on the measurements that are made on a population level. Clearly, our model is at present simplistic, as it describes a direct causality between EGFR activation, Erk status and proliferation, when in reality, cross-talk from other signalling pathways may modulate at all levels. Nevertheless, the example used to illustrate our novel approach to modelling heterogeneity within populations is biologically-relevant, and this single example of extrinsic heterogeneity could be extended to examine other nutrients, growth factors or cytokines in the environment. The second obvious extension of this idea would be to investigate the concept of intrinsic heterogeneity, where individual cells in a population differ in terms of their gene expression and numbers or concentrations of particular receptors, ligands or intracellular species. An example of this could be the situation where an E-cadherinnull cell arises by mutation and the fate of the population subset is determined relative to its interactions with the rest of the population. This could have important potential implications for understanding how proliferation controls are lost in the event of malignant transformation and how such malignant cells behave in comparison to their normal counterparts during carcinogenesis.
By combining a deterministic, ODE-based model of a signalling pathway with a non-deterministic agent-based representation of a multi-cellular population, we have developed a simulation environment that can be used to compare how emergent heterogeneity in the cell population can impact on system level behaviour. Our model predicts that the heterogeneity of cell contacts, arising as a result of the inherently stochastic processes of plating cells in culture and the subsequent random migratory behaviour of these cells, provides a basis for heterogeneity in biological response within a 'homogeneous' population. In this paradigm, factors that modulate the random nature of the cell contact may have a major role in influencing the behaviour or phenotype of emergent subpopulations. We have also demonstrated how an agent-based computational approach has a role as an iterative, predictive tool by generating 'virtual assays' to bring new insight into the interpretation of cell biological experiments.
The authors would like to thank RCUK, who provide the Academic Fellowship held by DW. NTG carried out experimental work as a postdoctoral research fellow funded by the Engineering and Physical Sciences Research Council (EPSRC) on grant GR/S62338/01. JS holds a research chair supported by York Against Cancer. Gemma Hill produced the time lapse movie in additional file 2.
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.