Cellular automata and integrodifferential equation models for cell renewal in mosaic tissues

Abstract

Mosaic tissues are composed of two or more genetically distinct cell types. They occur naturally, and are also a useful experimental method for exploring tissue growth and maintenance. By marking the different cell types, one can study the patterns formed by proliferation, renewal and migration. Here, we present mathematical modelling suggesting that small changes in the type of interaction that cells have with their local cellular environment can lead to very different outcomes for the composition of mosaics. In cell renewal, proliferation of each cell type may depend linearly or nonlinearly on the local proportion of cells of that type, and these two possibilities produce very different patterns. We study two variations of a cellular automaton model based on simple rules for renewal. We then propose an integrodifferential equation model, and again consider two different forms of cellular interaction. The results of the continuous and cellular automata models are qualitatively the same, and we observe that changes in local environment interaction affect the dynamics for both. Furthermore, we demonstrate that the models reproduce some of the patterns seen in actual mosaic tissues. In particular, our results suggest that the differing patterns seen in organ parenchymas may be driven purely by the process of cell replacement under different interaction scenarios.

Figure 1.

The updating process for the asynchronously updated cellular automata. (a) Description of the process for the CA based on the majority conceptual model, while (b) describes that of the CA based on the single-cell conceptual model.

Figure 2.

Solution of the cellular automata (CA) as described above, at times t = 0, 10⁷, 10⁸ and 10⁹. Here, we use mixed initial conditions by randomly assigning to each grid square the value of 0 or 1. In (a–d) left to right, we consider the majority conceptual model, and see cells quickly forming large agglomerations, before they either die out completely or dominate the grid. Out of 1000 runs, we see quick domination by a single cell type across the entire domain in about 30% of cases, while the remaining 70% end with both species present as seen in (e–h). In (i–l), a solution to the single-cell conceptual model is shown, demonstrating persistence of both cell lines over time. For (a–d) and (e–h), cells are updated according to the process described in figure 1a, while the CA in (i–l) is updated according to the process described in figure 1b. All CAs are carried out on a grid of size 256 × 256, and use an eight-neighbour Moore neighbourhood. On the boundary and at the corners of the domain, only cells in the domain are considered, and an average is taken over that reduced number of cells in the majority model, whilst in the single-cell model a cell that picks a neighbour outside the domain does not change state.

Figure 3.

Schematic of renewal function in the continuous model. (a) A smooth continuous approximation to a step function, symmetric about ½, representing the locally biased model. (b) Linear f, representing the locally unbiased model.

Figure 4.

The homogeneous mixed initial conditions used for the two-dimensional model. The two cells are evenly mixed across the domain. The domain is a square with sides of length 10 dimensionless space units.

Figure 5.

A solution of the two-dimensional locally biased model equations (2.1) with homogeneous initial conditions. We plot the proportion of cells of type A in space at dimensionless times t = 0, 10, 30 and 250. (a–d) The domain quickly evolves to an all white domain of A cells; (e–h) quick dominance by B (black). (i–l) Dominance by neither cell type, which held for long times (solutions were found to be stable in runs up to t = 10⁸; not shown). Out of 20 runs, A dominated 10 times, B 4 times, and neither 6 times. We begin with initial conditions of a = 0.5 + 0.02 × c where c is chosen randomly between 0 and 1 at each numerical grid point. In (a–d) and (e–h) the function f is given by f(I) = 0.5 tanh(tan(Iπ − 0.5π)) + 0.5, a continuous approximation to a step function, while f = I in (i–l). The dimensionless parameter values are R = 1.0, α = 1.0. The domain is of size 10 dimensionless space units. We set absolute error tolerance in the rowmap scheme to 10^{− 6}.

Figure 6.

Time evolution plots of the two-dimensional locally unbiased model equations (2.1). The domain is initially an even mix of the cell populations A and B, as in figure 4, and is of size 10 dimensionless space units. We plot the proportion of cell type A across space in the x direction for y = 2 at various times t, until t = 100. We see the proportion of A spreading homogeneously across the domain over time until A = 0.5 everywhere. All parameter values and numerical details are as in figure 5, although with a linear f as stipulated by this model (see figure 3b).

Figure 7.

A solution of the two-dimensional locally biased model equations (2.1) with mixed homogeneous initial conditions. We plot the proportion of cells of type A in space at dimensionless time t = 50. We see the formation of stable stripes. The domain is 40 dimensionless space units wide, and 4 high. All other numerical details are as in figure 5.

Figure 8.

More initial conditions (i.c.s) used for the two-dimensional model. The domain is again a square with sides of length 10 dimensionless space units. (a) Split i.c.s, (b) curved i.c.s and (c) island i.c.s.

Figure 9.

Solutions to the two-dimensional locally biased model equations (2.1) with curved initial conditions as in figure 8b. We plot the density of cells in space at (a) t = 10, (b) t = 20 and (c) t = 1000. Neither cell type dominates over long times. All numerical details and parameter values are as in figure 5.

Figure 10.

Solutions to the two-dimensional locally biased model equations (2.1) with island initial conditions as in figure 8c. We plot the density of cells in space at (a) t = 1, (b) t = 2 and (c) t = 5. The dominant cell type rapidly engulfs the ‘island’, demonstrating that we will not see spotted patterns. All numerical details and parameter values are as in figure 5.