[3] Cell-Invasion : 2D Monte-Carlo simulation

We have implemented a Monte-Carlo simulation program for studying tumor cell invasion into a half space of complex bio-material (as a part of a new project described in a recent post and its followups). In the following, the tumor cells are treated as active agents in a multi agent system.

• Space is discretized as a 2D rectangular grid of sites (x,z). The positive z-direction corresponds to the invasion direction, pointing into the material.

• Each site can be occupied by at most one agent. Multiple occupation is forbidden (in a way that is analogous to the Hubbard model used in the theory of metal insulator transitions).

• Each site is assigned a potential energy U(z,x), reflecting the difficulty for an agent to occupy this site. It would be possible, in principle, to investigate potential landscapes of arbitrary complexity. Simple choices would be a completely flat potential, or a binary valued potential without spatial correlations, i.e. binary white noise. The two possible values of the local potential could be zero (no obstacle) and infinity (impenetrable obstacle). And so on.
  

• There is a certain probability per unit time that an agent moves from its present site to a neighboring site (This leads formally to a kind of "hopping rate", although the agents are thought to move continuously in reality. It is the discrete, "pixelized" detection of their position that implies discrete hopping steps). The maximum hopping rate is R0. However, when the agent moves from a site of low to a site of high potential energy, the hopping rate is reduced by a Boltzmann factor with some effective temperature (Of course, tumor cell invasion is not driven by thermal fluctuations, but by active motor forces. Nevertheless, the cell migration rate will slow down when the cell tries to move into obstacles, i.e. when it enters high energy sites in the potential landscape. This effect can be conveniently modelled by the Boltzmann factor.)

• For the horizontal (z) and vertical (x) edges of the simulation area, boundary conditions can be chosen between hard walls and periodic in our implementation.

• Many agents can be in the system simultaneously. They can interact with each other (in order to describe processes such as binding of tumor cells etc.). They can also change the potential landscape (describing, for instance, stigmergic effects such as the protolysis of collagen by chemicals emitted by the tumor cells).

• The Monte-Carlo algorithm is "exact": In each cycle, a complete list of all possible processes (all possible hops of all agents, etc.) is generated and the rates of these processes (depending on the current system state) are computed, as well as the sum of all rates, which determines the average waiting time dtMean until the next process happens. The MC algoritm then chooses a random, exponentially distributed waiting time according to dtMean. It also chooses one of the processes from the list, with the correct relative probabilities, according to the rates. The process chosen is executed, the system state is updated, the relevant variables are recorded, and the next cycle starts.

As a simple example, let us consider a MC simulation on a grid of 1000(x) times 100(z) sites. In the x-direction, we chose periodic boundaries, but in the z-direction the boundaries are hard (This is realized by placing an almost infinitely high potential wall just outside and along the west and east edge of the potential landscape). Initially, 1000 cells are placed onto the west surface (z=0), forming a complete monolayer. The cells then start to diffuse randomly (This is simply achieved by a flat potential landscape U(z,x)=const. inside the whole simulation area. Then the hopping rates are everywhere homogeneous and isotropic). The picture below shows a snapshot of the distribution of agents over the simulation area after 100 time units (black circles) and after 500 time units (red squares). One can see how the agents are invading the half space.