Control of cell colony growth by contact inhibition

We investigate the dynamics of a colony of crawling, proliferating cells with a minimal, mechanical cell model. The cells consist of two disks, modelling the cell body and a pseudopod, connected by a finite extensible spring. The cells exhibit locomotion due to a linear coupling of the motility force to the cell extension. With a simple mechanism for contact inhibition of proliferation, we find the typical regimes of colony growth, with exponential growth at short times turning into sub-exponential growth at long times. In the latter regime the colony boundary moves outwards with a constant speed. We identify simple scaling relations for both regimes and the crossover between them. We find that the shape of the cells (the ratio between the two disk's radii) determines the efficiency with which cells orient themselves away from the colony. The better the cells are at aligning themselves away from the colony, the faster the colony expands.

contact inhibition is a cell property that limits the migration and proliferation of cells in crowded environments. Here we investigate the growth dynamics of a cell colony composed of migrating and proliferating cells on a substrate using a minimal model that incorporates the mechanisms of contact inhibition of locomotion and proliferation. We find two distinct regimes. At early times, when contact inhibition is weak, the colony grows exponentially in time, fully characterised by the proliferation rate. At long times, the colony boundary moves at a constant speed, determined only by the migration speed of a single cell and independent of the proliferation rate. further, the model demonstrates how cell-cell alignment speeds up colony growth. our model illuminates how simple local mechanical interactions give rise to contact inhibition, and from this, how cell colony growth is self-organised and controlled on a local level.
Cells move collectively and proliferate [1][2][3][4][5][6] as the embryo develops during morphogenesis 7 , as cancer spreads or as wounds close [8][9][10] . The way in which migration and proliferation interact with each other is complex 11 . Essential for the regulation of these processes is contact inhibition of locomotion (CIL), which describes the tendency of cells to stop migration or change direction when coming into contact with other cells [12][13][14][15][16] . CIL has been shown to enable cells to collectively migrate 14 , follow chemical gradients 17,18 , and disperse. It is now believed that CIL helps the control of collective tissue migration 16,17,19,20 , tissue growth 21,22 , morphogenesis 14 , wound healing and tumour development 23 .
Potentially distinct from CIL is contact inhibition of proliferation (CIP) which refers to the suppression of cell divisions in dense regions of tissues [24][25][26] , which in turn regulates their growth. There is evidence that CIP does not require direct cell contact [27][28][29][30] and as a consequence the effect is also called "density dependent inhibition of cell growth" 31 .
Modelling approaches for cell migration are manifold 32 , ranging from strongly idealised models for single cells crawling on substrates [33][34][35] to cells in confluent tissues 36 to continuum theories [37][38][39] . The collective behaviour of cells is under intense study and many questions about contact inhibition are still open 16,26 .
In order to reduce the complexity of the systems under study, it is valuable to investigate well controlled model systems. One such model system deals with the crawling and proliferation of a monolayer of cells seeded onto a substrate. How a few cells develop into an extended colony has been investigated in a recent experiment 21 , as well as in simulations, e.g. 22,[40][41][42][43][44][45][46][47][48] . The colony in said experiment consists of of Madin-Darby canine kidney epithelial (MDCK) cells 49 and exhibits two growth regimes. In the beginning the colony's area grows exponentially with time, followed by sub-exponential growth. Similar growth stages are also found in tumour growth dynamics, e.g. 50,51 In the former regime, cells are highly mobile and divide frequently, while in the latter, both the motion and proliferation of the cells becomes suppressed, linking the transition to contact inhibition. Such a transition is readily found in simulations if the proliferation rate of cells is locally coupled to density or available space 40,41,44 , or stress 37 .
Previously, we developed a minimal, mechanical model for crawling cells in which contact inhibition of locomotion arises from the internal mechanics of the cells 52 . Our model focuses purely on the cell mechanics, since mechanical forces inside of and between cells are now understood to be of crucial importance for the understanding of cell dynamics 9,21,53-58 . Our model differs from previous models in that it incorporates a coupling of motility to cell polarisation, and asymmetric cell shapes into a minimal model; both of which influence mechanical interactions between cells. We found that the model naturally exhibits a range of realistic behaviours, including the emergence of collective migration, without including an explicit alignment term. Our model is thus a natural candidate to investigate colony growth.
In this paper, we extend our model to include cell proliferation in such a way that the motility and division dynamics are entirely governed by the internal dynamics of the cells. The cells cycle between a motile phase and a division state. In the division state, the cells attempt to proliferate by making space for two cells; otherwise, the cell division is aborted. This naturally gives rise to a form of contact inhibition of proliferation. Similar approaches have been used e.g. by 40,44 .
We find that our model reproduces the typical colony dynamics; with exponential growth at short times turning into sub-exponential growth with a constant boundary speed at long times. Coinciding with this transition, the average cell speed decreases strongly, because of CIL occurring in the inside of the colony. As a result of contact inhibition, cells close to the boundary have higher speeds and proliferation rates. By varying the proliferation rate over a wide range, we then identify simple scaling relations for both regimes and the crossover between them. We had previously found that cell shape has a strong effect on cell collisions and that cells with large fronts align and coherently migrate 52 . In this work, we now demonstrate that cell with large fronts orient themselves away from the colony, which enhances the speed of colony expansion.

Model
We build on a model for crawling cells 52 , and include cell division. Each cell consists of two disks with distinct roles, see Fig. 1. One models a pseudopod and is at the front of the cell (index f ), the other disk represents the cell body and is at the back of the cell (index b). The positions of the disks r f and r b define the distance between the elements and the orientation of the cell → = → − → r r r bf f b . The dynamics of the cells are coarse-grained over the typical idealised crawling cycle 59 . The substrate exerts a drag force ζ − → v i i on the disks, with → v i being the velocity of the disks ( ∈ i f b , ) and ζ i being the respective drag coefficients. For simplicity, we set ζ 1 = ζ 2 = ζ. We employ the commonly made assumption that the friction with the substrate is the dominant contribution to friction in the system, e.g. by 60 , and neglect intracellular friction and friction between different cells.
It is now understood that the shape of cells can be highly variable and has a strong influence on their collective behaviour 36,[61][62][63][64][65][66] . In our model, the cell shape is determined by the configuration of the two disks with diameters σ f and σ b , and can be understood as representing a statistical average over the real cell shape. Such a coarse graining of the highly variable shape of real cells is a promising approach for cells whose shape does not deviate too much from the average, e.g. for epithelial cells. In our previous work, we showed that the present model exhibits an alignment transition in coherent migration as a function of cell shape 52 . Cells with a small front disk (as compared to the back disk) are bad at moving apart after colliding, while cells with a large front disk push each other out of their path, which over time leads to global alignment. In this sense, cells with a large front disk in our model exhibit another form of contact inhibition. For this work, we mainly investigate cells with a shape ratio of σ σ = .
/ 079 b f which exhibits realistic alignment and migration behaviour 52 . In our previous work, interactions between cells were purely repulsive. Here, we introduce new, adhesive interactions between cells. The adhesiveness of the potential is characterised by its well depth ε well in respect to the potential height ε core . The force acting on a disk α by other cells is denoted α F cc , (α ∈ b f [ , ]). For details, see Methods.
cell life cycle. Each cell switches independently between two states, a motile state and a division state. The duration of the motile state is determined when switching to it from the division state, by drawing a random number from a normal distribution with mean T mot and standard deviation T mot /2. With this distribution we model natural fluctuations in cell cycle duration. The duration of the division state is held constant at T div . On average, the duration of the whole cell cycle is thus = + T T T div m ot . In the motile state, the cell behaves as in ref. 52 . The front disk exerts a migration force F mig in the direction of the orientation of the particle, while the back disk is passive. The migration force is given by bf mig with motility strength m, see Fig. 1(a,d). The connection between the disks is modelled as a finitely extensible nonlinear elastic (FENE) spring 67 , which gives rise to a contracting force between the disks bf bf bf fene max 2 with coupling parameter κ determining the strength of the contraction, and R max the maximum distance between the two disks.
In the division state, cells attempt to make space for two daughter cells. Cells only divide if at the end of the division state the cell extension reaches a division threshold R div . Coupling cell proliferation to cell area is an idealisation of the observation that larger cells divide more frequently in experiments 21,68 and models contact inhibition of proliferation.
Cells elongate by having both disks enact the same migration force F div with opposite sign. For convenience we set = F F div m ig . The contracting force between the disks remains unchanged. Cells do not migrate in the division state. To make space for the new cells, we increase the size of the smaller of the two disks linearly until it matches the size of the larger disk at the end of the state. After a successful division we construct two cells in place of the original cell's disks. We randomly displace the new cell's disks to randomise the orientation of the new cells. If < r R bf div at the end of the division state, the cell division is aborted, the cell contracts again, and the migration state is entered.

equations of motion.
For each of the cells, we now have two coupled non-linear equations of motion, with χ = 1 in the division state and χ = 0 in the motile state. Apart from the randomised positions of daughter cells' disks, our model does not include random forces. This is a reasonable assumption when collisions (and cell division) dominate the dynamics 69,70 . In the migration state, the cell is only motile when its disks have some separation, > r 0 bf , and thus when its shape deviates from a circle. This kind of coupling of motility and deformation is typical in migratory cells 71 .
If the cell is in the motile state for long enough, it then enters a steady state with constant extension r bf ss and constant speed v ss in which the forces acting on the cell balance 52 .

Results
colony growth. At first, we simulated cell colonies of non-adhesive cells to determine what results the simplest possible interactions yield in our model, see Fig. 2a) and Video S1 (Supplementary Materials). At early times, the colony grows exponentially, but eventually crosses over into sub-exponential growth. In the exponential regime, all cell division attempts are successful. Since cells attempt to double with a rate of = + T T T mig d iv and we always start with one cell at = t 0, the number of cells grows as In the experiment 21 , the sub-exponential growth is characterised by the boundary of the colony moving outwards at a constant speed. If R t ( ) is the radius of the approximately circular colony, then the outwards speed of the boundary can be extracted from the area The area of the colony can be calculated from the areas of all the individual cells, for which we use Eq. (11). At long times, the speed of the boundary saturates in our simulation to a constant speed as well, see Fig. 2b). We find the speed to be ≈ .
v v 0 2 B ss . The speed of cells allows quantifying the activity of the colony over time. In the exponential regime, the average cell speed is constant and then in the sub-exponential regime decreases over time, eventually dropping below the boundary speed, see Fig. 2b). The transition in the average cell speed occurs at the same time as the transition of the boundary speed. All of this is qualitatively similar to what is observed in the experiment.
Adhesive cells exhibit similar growth dynamics, see Fig. 2 and Video S3 (Supplementary Materials). The slope of the exponential regime is the same, with all divisions being successful at early times. The average cell speed in the exponential regime remains unchanged as well. However, the exponential regime only extends until the colony consists of tens of cells. The transition to a constant boundary speed takes much longer, as does the slowing of the average cell speed. On average, we find cells to be faster, but the colony to expand slower.

Radial analysis.
To understand the growth of the colonies in more detail, we look at the spatial distribution of the following key quantities: cell density, cell divisions, and cell speed. For this we analyse one exemplary simulation with non-adhesive and adhesive cells. Non-adhesive cells form colonies with a diffuse boundary, with some cells even escaping, Fig. 3a), whereas adhesive cells form a denser colony, with no cells escaping the bulk, see Fig. 3b).
For the non-adhesive cells, we find that in the exponential regime ( < t 150) the density is quite low but later quickly reaches a high value in the bulk of the colony, Fig. 3c). The boundary of the colony remains diffuse with a low density and constant width over the whole simulation. The colony of adhesive cells is different at early times, see Fig. 3d), because it is dense and cohesive from the beginning. At long times, the adhesive colony appears qualitatively quite similar to the non-adhesive colony, with the same bulk density, albeit with a sharper boundary.
We then calculated the radial distribution of successful cell divisions N r t ( , ) div , see Methods. In the exponential regime, all attempted divisions are successful, and thus are distributed homogeneously over the colony, Fig. 3e,f). For non-adhesive cells, divisions at later times mostly occur in a ring of constant thickness at the boundary of the colony, where there is more space available to the cells, see Fig. 3a,e). This is a result of our cell division mechanism which naturally gives rise to contact inhibition of proliferation. Cells in the bulk cannot make the necessary space for the cell division to occur, except for the colony centre, where space becomes available as cells migrate outwards. For adhesive cells we find a similar pattern, see Fig. 3f), where a ring of increased proliferation probability can still be discerned. However, the probability of a division occurring at the boundary is considerably reduced, because the local density tends to be higher as compared to the non-adhesive colony. In addition, cell divisions are more frequent in the bulk of the colony. The reason for this becomes clearer with an analysis of the cell speeds.
Non-adhesive cells are mostly mobile at the boundary of the colony, with motion being strongly suppressed by CIL in the colony bulk, see Fig. 3a). The cells on the boundary are on average pointing away from the colony centre, but there is considerable local variation due to contact inhibition of locomotion after collisions and noise introduced by cell divisions. Cells can momentarily obtain high speeds after a successful cell division when they move away from the other daughter cell.
In the exponential regime, all cells are mobile, but in the sub-exponential regime, only cells at the boundary exhibit speeds close to the single-cell steady state speed v ss , Fig. 3g). Motion in the bulk is suppressed strongly by contact inhibition of locomotion. In comparison, adhesive cells are slower on the border but more mobile in the bulk and more aligned with each other, see Fig. 3b,h). We attribute this to the cells at the boundary being held back by the cells at their back and in turn the cells of the bulk being pulled outwards by the boundary cells. This is commonly called the "tug of war" between cells 9,57,72 . The tug of war leads to a different bulk structure between the colonies, even though the densities are similar, Fig. 3a,b). There is more free space in the non-adhesive colony and cells are more compressed due to contact inhibition of locomotion. The adhesive colony, on the other hand, is fully cohesive, with cells always being at contact and pulling on each other and therefore, on average, more extended. This is reflected in the higher cell speeds in the bulk of the adhesive colony. As a result, it becomes easier for the adhesive cells in the bulk to reach the division threshold. In conclusion, we find that cell-cell adhesion considerably alters the colony structure on a local level, while leaving the qualitative colony dynamics with the two growth regimes unchanged. www.nature.com/scientificreports www.nature.com/scientificreports/ Variation in the cell cycle. We illustrate in the following that at short times, the colony dynamics is entirely determined by the proliferation rate, whereas at long times, the dynamics are determined by the migration properties of the cells. For this purpose we vary the durations of the migration and division states of the cells, Fig. 4a). For comparison, we also simulated the situation where the cells are not migrating at all, but divide regularly.
Since the colony area grows as at early times, we find that all simulations collapse onto one asymptote when time is rescaled by = + T T T mot d iv , see Fig. 4b). Therefore, while contact inhibition of proliferation is weak, the dynamics are purely determined by the rate of cell divisions 1/T. As a consequence, it also does not play a role whether the cells are migratory or not.
However, this simple rescaling cannot collapse the data at long times, as the time at which growth becomes sub-exponential varies between the simulations. To calculate those crossover times, we make use of the www.nature.com/scientificreports www.nature.com/scientificreports/ observation that the boundary speed of the colonies is limited. In the exponential regime, the boundary speed increases as, see Eqs. (3) and (4), The crossover time ⁎ t is reached when the colony reaches the terminal boundary velocity which we assume to be proportional to the steady state speed of a single cells, = ⁎ v Cv B ss with some factor C. Then we obtain From the data, we find that = .
C 0 25 holds generally for migrating cells, see Fig. 4d). With this value for C, the crossover times correctly mark the transitions to sub-exponential growth for all simulations, see Fig. 4b). From Eq. (6), we see that the crossover is most strongly influenced by the duration of the cell cycle, followed by the crawling speed of the cells. At the crossover, the colony is of size = At long times,  ⁎ t t , the radius of the colony grows with constant speed Cv ss , so we find that approximately it must hold that ss Rewritten for the colony area, we have ss This long-time scaling is exposed in Fig. 4c). Most of the data collapse onto a master function that is close to Eq. (10). The two simulations that deviate from this are the two cases in which the cells do not actively migrate and thus necessarily violate the scaling.
We therefore find that if the cells are actively migrating, the colony boundary moves at a constant speed determined by cell motility, and that if the cells are not migrating, long time growth becomes severely suppressed. Or stated differently, regardless of the proliferation rate, at long times the colony expands as fast as the cells on the boundary are able to migrate away from the colony centre.
Notably, the boundary speed is smaller than the steady state speed of a solitary cell v ss by the factor C. This has two reasons: (1) the cells on the boundary adhere to the rest of the colony and are pulled back by them (2)  . This parameter choice therefore presents our best match to the experiment. For = T 64 mot , we find that the crossover is reached at a colony size of ~1000 cells after ≈10 full cell cycles, see Fig. 4b). This compares well with the 5-6 days, i.e. ≈7 full cell cycles, that the cells undergo in the experiment until the crossover occurs at a similar colony size, considering that our colonies originate from a single cell while the experiment starts with multiple cells.
Influence of cell shape. We investigate the interplay of the alignment transition with colony growth, see Fig. 5. The exponential growth regime is found to be independent of shape, see Fig. 5a), while at long times colonies with large-front cells tend to expand more rapidly. Colonies for σ σ = .
/ 067 b f are about 3 times larger than for / 126 b f at the end of a simulation. For all times, cells with larger fronts tend to be faster, see Fig. 5b). This is because cells with larger fronts are better at aligning away from the colony, as the average orientation of the cells in respect to the colony ⋅r r bf shows, Fig. 5c). Cells with small fronts tend to be slightly aligned towards the colony centre. If the cells have large fronts, they tend to be aligned away from the colony, with alignment getting more pronounced towards the boundary of the colony. In the growth of a colony, the alignment of the cells due to shape acts as a mechanism for orienting the cells away from the colony, which is what is expected from cells exhibiting contact inhibition of locomotion 16 . As a result, the colonies of cells with large fronts expand much more quickly in the sub-exponential regime, see Fig. 5d) with the boundary speed being about twice as large for σ σ = .
/ 067 b f than for σ σ = . neighbour information. The colony can be further characterised by a measurement of the local structure.
We calculated the coordination number n, i.e. the number of nearest neighbours, for each cell. For each simulation snapshot, we calculate a Voronoi tessellation of the cell centres → + → r r Fig. 6(a). To cut off the large Voronoi cells on the boundary of the colony, we intersect the Voronoi cells with circles of radius σ + . R ( 15 )/2 b max which represents the maximum interaction radius from the cell centre. With this corrected Voronoi tessellation, we calculate n for all cells and then calculate the histogram p n ( ), see Fig. 6(b-d). The data is shown for multiple times relative to ⁎ T for each T mot . Additionally, we show a histogram at a late stage in the simulation. The results are mostly independent of T mot presented in this way.
As time progresses, cells with 6 neighbours become quickly the most prevalent, while having 3 or 4 neighbours becomes less likely. Cells with 5 neighbours lightly decline, and 7 neighboured cells become more likely. This is quite similar to the MDCK experiment 21 and other cell types 73 . The notable difference is that our simulations overestimate configurations with 6 neighbours at long times. This is directly visible in simulation snapshots as crystallisation of the tissue. This could have many reasons such as a lack of polydispersity of cell sizes, which would hinder crystallisation, or could be caused by the strength of contact inhibition of proliferation present in the model which leads to arrest of the cells.

Discussion
We investigated the dynamics of a colony of crawling, proliferating cells with a minimal, mechanical cell model. With a simple mechanism for contact inhibition of proliferation (CIP), we find the typical regimes of colony growth, with exponential growth at short times turning into sub-exponential growth at long times. The latter regime is characterised by the colony boundary moving outwards with a constant speed. We identify simple scaling relations for both regimes and the crossover between them.
We find that the crossover colony size ⁎ A between the two regimes is a marker of contact inhibition of proliferation and that the boundary speed v B in the sub-exponential growth regime expresses the strength of contact inhibition of locomotion. From Eq. (8), we see that ⁎ A depends on the steady state speed of the cells v ss and the duration of the cell cycle T . The faster the cells can migrate and the longer the cell cycle is, i.e. the longer the time between division attempts, the larger the colony can become before CIP sets in. The long-time behaviour of v B only depends on the speed at which cells are able to travel, and thus primarily measures contact inhibition of locomotion (CIL), see Eq. (10). This is corroborated by the investigation of cell shape on colony dynamics. There, the crossover ⁎ A is unchanged by cell shape, but large-front cells are better at aligning away from the colony, and as a result, tend to be more mobile and tend to expand the colony faster at long times. It remains to be seen how close the observed faster expansion for large-front cells in the present model is to contact enhancement of locomotion 74 .
The mechanisms of CIL and CIP as exhibited by our model are idealised compared to the experimental situation. While in the experiment the bulk of cells is still mobile, our cells arrest more quickly and tend to crystallise. Also, while the mechanism of contact inhibition of locomotion -which arises from the motile force being coupled to the cell polarisation as well as from the shape of the cell itself -controls alignment of the cells away from the colony and colony growth, our cells are unable to expand the colony at their full crawling speed. We also neglected the polydispersity of cells sizes and noise.

parameters.
Introducing cell-cell adhesions made it necessary to modify the model's parameters with respect to ref. 52 . We needed to shorten the maximum extension of the cells so that the disks could not be pulled apart far enough to leave a space in-between. The changed parameters also yielded cells that more readily expand from their symmetric, circular state. We confirmed that modifying the parameters did not qualitatively change the dynamics discussed in ref. 52 . All of the results reported here are for σ = . . The time step of the simulations is . ⋅ − 2 1 10 3 . The simulations are initialised with a single cell in one of the states randomly. We average our results over 10 individual simulation runs with different seeds for the pseudo-random number generator used for the calculation of the state durations and randomisation of cell orientation after a cell division. cell area. The area A of a cell with a back particle of diameter σ b , a front particle of diameter σ f and a distance r bf between the particles is given by between the disks, potential minima can form due to superposition of the elements' potentials at certain positions along the boundary of the cells. Those can be either behind the smaller cell element or on the sides of the cell, see Fig. 7a). This can lead to alignment of cells purely due to the potential, which we want to avoid. Additionally, cell elements can interact with other cells' elements with which they are not in contact and even through the other cell element, since the well depth of the potential is comparable to the range of the potential, see Fig. 7b). This can cause cells which are in contact to contract each other, see Fig. 7c). In contrast, real cell interactions occur due to direct contact of cell membranes. Therefore, to accurately model cell-cell interactions we must find a way to only have surface-surface interactions as well.
We implement surface-surface interaction by letting two cells only have up to one adhesive interaction, between the two elements whose surfaces are closest, see Fig. 7(d,e). To keep the model similar to our previous paper, we still let all the cell elements repel each other.
Finally, using the Lennard-Jones potential had one additional drawback. The width of the potential well scales with the elements' radii, which leads to a discontinuity when only the closest surface is considered adhesive and Scientific RepoRtS | (2020) 10:6713 | https://doi.org/10.1038/s41598-020-62913-z www.nature.com/scientificreports www.nature.com/scientificreports/ the elements have different sizes, see Fig. 7f). This was not an issue to our previous work in which the cutoff was set to make the forces purely repulsive, which makes the interaction very short ranged.
We solved these issues by switching to a soft-core potential, in which the width of the well can be set independently from the element radius, which yields a potential landscape such as the one shown in Fig. 7g). We first define a cubic helper function t y ( ) to serve as a continuous step. The function describes a continuous and monotonous step between = t(0) 0 and ξ = t( ) 1, in which the end points are saddle points , the force and the potential between two cell disks at separation → r are given by www.nature.com/scientificreports www.nature.com/scientificreports/ See also Fig. 7(h,i). The minimum of the potential is at σ = r , identifying σ as the characteristic length scale of the potential, see Fig. 7h). ε well is the energy in the minimum, i.e. φ σ ε = − ( ) well . The distance between the two inflection points of the potential (the extremes of the force) is given by ξ, which is therefore a measure of the width of the potential well. We require ξ σ ≤ ≤ 0 to enforce that ε core is the energy for = r 0, i.e. φ ε = (0) core , and that = F(0) 0. The natural cutoff of the potential is σ ξ = + r cut . Consider now two cells with elements α and β, respectively (α, β ∈[b, f]). The separation between any pair of elements α β ( , ) is given by → = → − → . With this soft-core potential, the potential well around a cell is now of constant width and of shorter range, see Fig. 7g). For all our simulations we chose a width of the potential well of ξ = .