Simulation of the growth and differentiation of stem cells on a heterogeneous scaffold

We present Monte Carlo simulations illustrating that scaffolds with specially designed heterogeneities can be used as a tool for governing the growth and differentiation of stem cells.

This tool is predicted to be efficient provided that the size of the heterogeneities exceeds (but not appreciably) the cellular size and the rate of cell diffusion is relatively fast compared to that of cell division.


Depending on the source, stem cells have the potential to form one, many or all cell types of an organism (the diversity and abilities to evolve in different directions decrease from embryonic stem cells to “adult” stem cells) and accordingly can in principle be employed to build blocks for every tissue we comprise.1

Using this potential and bringing stem-cell therapies to the clinic demand a better understanding of the factors that maintain cells in a multipotent state or drive them to create differentiated cells.1,2

Specifically, there is a need to improve the control of differentiation of cells and to promote desirable division via external signals, in order to reach suitable spatial cellular arrangements.2–4

An additional challenge is to extend and sophisticate the tools which can be employed to govern the tissue evolution.

One of the opportunities here is to use specially designed heterogeneous scaffolds, on which the original cell culture(s) are deposited and eventually develop into a tissue5 (for other strategies, see, e.g., ref. 6).

To illustrate this concept, we present the first Monte Carlo (MC) simulations of the growth and differentiation of stem cells on a scaffold with μm-sized heterogeneities (such scaffolds can be fabricated by using, e.g., various lithographic techniques; for recent examples of experimental studies of behaviour of cells on heterogeneous surfaces, see ref. 7).

Physically, the cellular growth is complicated by cell–cell adhesion resulting usually in aggregation of similar cells and/or segregation of different cells.

Such phenomena can be treated by using a lattice approximation.

To our knowledge, the first lattice MC simulations of cell sorting were executed by Graner and Glazier.8

At present, there are various versions of lattice MC models focused primarily on the 2D cell segregation and/or growth (see, e.g., papers refs. 9–14 and references therein).

In particular, the proliferation and differentiation of stem cells on an uniform surface was analysed in refs. 13 and 14.

Our aim here is to simulate the evolution from a single stem cell into various cell types, on a heterogeneous surface, where the cell division/differentiation occurring according to certain rules is influenced by the surface.


In our present coarse-grained lattice model, cells are located on a 2D L × L square lattice.

(Alternatively, one could use, e.g., a triangular lattice.

It would not change the conclusions, because the segregation and/or island-growth kinetics are not too sensitive to the type of a lattice.) Each lattice site is either vacant or occupied by one cell (for the multi-site models making it possible to mimic changes of the cell shape, see refs. 8 and 11).

The symbol i (or j) indicating a cell type runs from 1 to 3.

Cells of type 1 are assumed to be totipotent (in our context, the terms “totipotent” and “pluripotent” are used interchangeably), i.e., they are able to give rise to all cell types.

Cells of types 2 and 3 represent differentiated cells.

Specifically, the latter cells are considered to be able to produce only cells of its own kind (e.g., the division of a cell of type 2 results in the formation of two cells of this type).

The cell size is assumed to be slightly larger than the site size but smaller than the distance between the next-nearest-neighbour (nnn) sites.

With this prescription, the interactions between nearest-neighbour (nn) and nnn cells, εijnn and εijnnn, should be repulsive and attractive due to deformation and adhesion, respectively.

These interactions, introduced to take into account softness, deformation, and adhesion of cells, are of course “effective”.

(In reality, the cell adhesion may be related, e.g., to interactions via membrane proteins.15

Explicit inclusion of such details is beyond our coarse-grained description.)

Cell diffusion is a complex process including deformation of the cell shape.16

In our coarse-grained treatment, we skip the details and realize this process via cell jumps to vacant nn sites.

A jump is allowed provided that in the final state (after the jump) the cell has at least one contact with other cells, i.e., after a jump at least one of the nn or nnn sites should be occupied (this condition is introduced in order to suppress rare events of escape of single cells from the central area where they are located; its role in the simulations is minor, because the probability of such events is low anyway).

The normalized jump probability is defined by the standard Metropolis rule.

In the case under consideration, the Metropolis algorithm is appropriate, because it correctly describes the tendency of cells to increase the number of adhesive contacts.

One should however bear in mind that in our context the thermal energy, kBT, used in the Metropolis rule, is effective as well as the cell–cell interactions (for this aspect of the simulations, see ref. 10).

The mechanism of cell division is governed by the spatial constraints and cell–cell signaling17 including messenger transport between cells, signal transmission via the cell membrane, gene regulation, etc.

In our model, this complex process is mimicked by introducing prescribed probabilities of cell division, which are different for different arrangements of cells and depend on the underlying surface to which the cells are attached.

In particular, the division of a given cell is allowed provided that (i) all the nn sites are vacant (the cells located in nn sites are deformed and accordingly the probability of their division is assumed to be negligibly small) and (ii) the newly-born cell has at least one nn vacant site.

The signaling may occur via ligand diffusion between cells18 or via the so-called gap–junction communication.19

In the former case, the cells may be far from each other and the ligand concentration is usually fairly low.

In the latter case, the cells should be in close contact.

In open systems (e.g., on surfaces), ligand diffusion away from cells may easily result in decrease of the ligand concentration and accordingly the role of the long-range signaling may be less significant.

In our present simulations, we take into account only short-range signaling in order to keep the model as simple as possible (in principle, the long-range signaling can be incorporated as well as described in ref. 14, but this is beyond our present goals).

In particular, the range of cell–cell signaling is considered to be limited to next-nearest neighbours.

Thus, the division probability may depend on the occupation of nnn sites (nn sites are empty; see above).

This approximation is reasonable at the initial stages of cellular growth when the number of cells is not too large (at the later stages, the growth may be influenced by long-range inhibitory interactions among the cells20 and/or limitations in the supply of oxygen and nutrients21).

Specifically, we use the simplest rules including spontaneous division of stem and differentiated cells and induced division of stem cells.

Spontaneous division of a cell of type i, occurring with probability pisp, is considered to result in the formation of two cells of type i.

The rules for induced division of a stem cell, due to communication with adjacent cells, may be rather complex.

To be specific, we assume that division of a stem cell, related to communication with a nnn cell of type i, results in generation of an additional cell of type i + 1 with probability pi+1ind.

For example, a stem (type-1) cell communicating with a type-2 cell can give rise to a type-3 cell.

(To validate this rule, we may refer, e.g., to differentiation of neural stem cells.22

In reality, the probabilities pisp and piind depend not only on the cell type but also on the growth factors added into the nutrition solution for the cell culture.) During trials of the induced division of stem cells, one of the nnn site is chosen at random and then the process is performed with the prescribed probability depending on the occupation of this site.

The scaffold surface may be fabricated of “conventional” materials with heterogeneities fabricated on the μm scale.

Another promising approach is based on the use of selectively-functionalized supported lipid bilayers, offering efficient reduction of nonspecific cell and protein binding, in combination with appropriate attachment sites.23

To mimic such situations, the lattice sites representing the scaffold surface are assumed to form patches.

Specifically, the patches of type 1, 2 and 3 are considered to interact more strongly with cells of type 1, 2 and 3, respectively.

The cells are allowed to diffuse on all the patches, but as soon as a cell of type i reaches the patch formed of material i, its movement becomes limited by this patch, i.e., jumps back to the other patches are no longer possible (this is a prerequisite for efficient cell segregation).

To reduce the number of the model parameters, the rates of diffusion of different cells on different patches are set equal (this rule is sufficient for our present goals but if necessary it can easily be relaxed).

To characterize the relative rates of cell division/differentiation and diffusion, we use the dimensionless parameter pdiv.

The rates of these processes are considered to be proportional to pdiv and 1 − pdiv, respectively.

With this specification, our MC algorithm consists of sequential trials to realize diffusion or division events.

A site is chosen at random.

If the site is vacant, the trial ends.

Otherwise, a cell located in the site tries to diffuse if ρ > pdiv or to divide if ρ < pdiv, where ρ (0 ≤ ρ ≤ 1) is a random number.

In both case, one of the nn site is selected at random and if the latter site is vacant the processes are performed with the prescribed probabilities by using the rules described above.

All the MC runs were started with a single stem cell located in the center of the clean lattice (with L = 150) at the patch formed of material 1.

During the simulations, we used the reflective (no-flux) boundary conditions.

Time was measured in MC steps (MCS).

One MCS is defined as L × L MC trials.

The duration of the runs was selected so that in the end there were about 3000 cells.

Results of simulations

In reality, the typical cell cycle time varies from an hour in an early embryo to about 20 h in adult stem cells.

The rate of diffusion strongly depends on the cell-substrate interaction.

Thus, the ratio of the rates of these processes is expected to be in a wide range.

At present, quantitative data on this subject are however scarce.

To get a realistic value of pdiv, we scrutinized diffusion, proliferation and differentiation of mouse neural stem cells on various substrate-coated plates24 (these observations were possible as a by-product of experimental studies of neural stem cells in our laboratory).

Diffusion was observed to be relatively fast compared to division.

Specifically, it appears to be reasonable to put pdiv = 0.01 or somewhat lower.

In our present simulations, cell division is chosen to be two or three orders of magnitude slower compared to diffusion, i.e., we employ pdiv = 0.01 or 0.001.

The energetic parameters are selected so that cells of different types tend to segregate and so that this process is not too slow.

Practically, this means that the ratio of the effective cell–cell interactions and kBT should not be too large.

In particular, we use εijnn/kBT = 1, ε11nnn/kBT = −1, εiinnn/kBT = −2 for i ≥ 2, and εijnnn/kBT = −1 for ij.

For spontaneous and induced division of cells, we employ p1sp = 0.1, p2sp = p3sp = 1, p2ind = p3ind = 0.1, and piind = 0 for i ≥ 4.

As a reference case, we have calculated [Figs. 1(a) and 2(a)] the cellular growth with pdiv = 0.01 on the uniform surface.

In this case, the growth of cells of type 2 is rather rapid.

In the very beginning, the birth of these cells is preferable compared to that of type-3 cells due to cell–cell communication between the stem cells.

Later on, the type-2 cells surround stem cells and effectively suppress generation of cells of type 3 due to the spatial constraints (usually there are no or only a few cells of the latter type).

Introducing heterogeneities does not change qualitatively the growth kinetics and lattice patterns [see, e.g., the results shown in Figs. 1(b) and 2(b) for a three-strip scaffold with a narrow (two-side) central strip].

Physically, this growth mode seems to be connected with insufficiently high rate of diffusion of cells (note that this process is partly suppressed due to attractive cell–cell interactions and spatial constraints).

Specifically, diffusion is able to maintain a cell arrangement which is locally close to equilibrium but not able to result in global redistribution of cells.

To increase the relative rate of cell diffusion, we have reduced pdiv down to 0.001.

For the uniform surface or in the case of relatively wide heterogeneities [e.g., for a three-strip scaffold with a six-site (or wider) width of the central strip], this modification has not changed the growth mode.

In particular, as a rule, the cells of type 2 again dominate [Figs. 3(a) and 4(a)].

On the scaffolds with narrow heterogeneities [Figs. 3(b–d) and 4(b–d)], the growth kinetics and lattice patterns are however found to be qualitatively different compared to those presented earlier.

Specifically, the increase of the relative rate of cell diffusion has resulted in effective segregation of the cells.

Under such circumstances, the birth of cells of type 3 by cells of type 1 is not suppressed.

For this reason, the population of cells of type 3 is appreciable (the numbers of cells of types 2 and 3 are comparable).

For pdiv = 0.001 this effect holds (not shown) with decreasing kBT and/or increasing the absolute values of the energetic parameters (if, e.g., εijnnn/kBT = −3 for i ≥ 2).

With increasing kBT (if, e.g., εiinnn/kBT = −1 for i ≥ 2), the cell motility increases, the rate of cell division increases as well (primarily for cells of type 2), but the driving force for segregation obviously decreases.

The former two factors do not compensate the latter one, and the cells of type 2 dominate both for pdiv = 0.01 and 0.001 (not shown), because as a rule the cells of type 1 are surrounded by cells of type 2 and the birth of cells of type 3 by cells is suppressed.


In summary, we have illustrated that scaffolds with specially designed heterogeneities can be used as a tool for governing the growth and differentiation of stem cells.

This tool is predicted to be effective provided that the size of the heterogeneities exceeds (but not appreciably) the cellular size and the rate of cell diffusion is relatively fast compared to that of cell division.