The recent advances in the mathematical modelling of human pluripotent stem cells

Human pluripotent stem cells hold great promise for developments in regenerative medicine and drug design. The mathematical modelling of stem cells and their properties is necessary to understand and quantify key behaviours and develop non-invasive prognostic modelling tools to assist in the optimisation of laboratory experiments. Here, the recent advances in the mathematical modelling of hPSCs are discussed, including cell kinematics, cell proliferation and colony formation, and pluripotency and differentiation.


Introduction
Human pluripotent stem cells (hPSCs) have the ability to self-renew indefinitely through repeated divisions (mitosis) and can differentiate into any bodily cell type (the pluripotency property). The latter property underpins their promising clinical applications in drug discovery, cell-based therapies and personalised medicine [1,2]. Amongst others, cardiomyocytes [3], pancreatic cells [4] and corneal cells [5] have all been successfully created from hPSCs. In the lab, hPSCs are grown in mono-layer colonies of up to thousands of cells ( Fig. 1) from which they can be directed for specific experiments or therapies, or expanded to produce further hPSC colonies. They occur either as human embryonic stem cells (hESCs) derived from the early embryo, or human induced pluripotent stem cells (hiPSCs) which are derived by the genetic reprogramming of differentiated cells [6]. The latter approach, which received the 2012 Nobel Prize in Medicine or Physiology for its discovery, offer patient-specific hPSCs without the ethical issues associated with hESCs.
Emerging biomedical technologies require the efficient, large-scale production of hPSCs [7]. Furthermore, applications of hPSCs in the clinic require great control over the pluripotency, clonality (the proportion of identical cells that share a common ancestry) and differentiation trajectories in-vitro. However, the existing procedures for large scale experiments remain inefficient and expensive due to low cloning efficiencies of 1% to 27% (the percentage of single cells seeded that form a clone) [8,9]. Understanding factors which promote the efficient generation and satisfactory control of hPSC colonies (and their derivatives) is a key challenge.
Mathematical and computational modelling allows the identification of generic behaviours, providing a framework for rigorous characterisation, prediction of observations, and a deeper understanding of the under-lying natural processes. The application of mathematics to biology [10] has led to many significant achievements in medicine and epidemiology (for example, predicting the spread of 'mad cow' disease [11,12] and influenza [13]), evolutionary biology [14] and cellular biology (descriptions of chemotaxis [15] and predicting cancer tumour growth [16]). Similarly, mathematical models are a powerful tool to further our understanding of hPSC behaviours and optimise crucial experiments.
The first mathematical model of stem cells, a stochastic model of cell fate decisions [17], has since been extended to include many other aspects of cell behaviour [18][19][20][21][22]. In particular, when such mathematical models are rigorously underpinned and validated on experimental observations, the reciprocal benefit for experimentation can be profound: an example is the development of an experimentally rained model of hiPSC programming, which led in turn to strategies for marked improvements in reprogramming efficiency [23].
Coherent mathematical models of hPSC properties may provide non-invasive prognostic modelling tools to assist in the optimisation of laboratory experiments for the efficient generation of hPSC colonies. Statistical analysis of experimental data allows the quantification of stem cell behaviour which can then inform the development of these models. Here we shall discuss recent advances in the mathematical modelling of hPSCs and their impact.
This review focuses mostly on hESCs, with some limited discussion of hiPSCs. We first outline some of the key properties of hPSCs before focussing on recent developments in mathematical models of the key properties: • Section 2: Key biological properties of hPSCs • Section 3: Cell kinematics. The movement of cells alone, in relation to one another and within hPSC colonies. • Section 4: Colony growth. Models capturing cell proliferation, with and without a spatial component. • Section 5: Cell pluripotency. Pluripotency regulation models, both intra-cellular and at the colony scale.
Finally, in Sect. 6 we provide a summary of the models discussed, their impact on biological experiments and the next steps for model development.

Key biological properties of hPSCs
The satisfactory understanding and control of hPSC evolution remains elusive due to their complex behaviour over multiple scales: the intra-cellular scale (processes happening within cells), the cellular or micro-environment scale (the environmental effects on individual cells) and the colony scale (collective cell behaviours throughout colonies), as illustrated in Fig. 2. Advances in imaging and molecular profiling (classification based on gene expression) have identified the core processes within the evolving colony [8,[24][25][26]. Here we outline some of these key biological properties across these scales and their relevance for mathematical modelling.

Intra-cellular scale
The key intra-cellular behaviours integral to hPSC modelling are the cell cycle and pluripotency regulation. The cell cycle is the timed series of events controlling DNA replication and resulting in a cell division. The phases of the cell cycle are: G1 (growth phase), S (synthesis phase in which DNA is replicated), G2 (further growth) and M (mitosis, the cell division). The G1 phase is shortened for hPSCs, leading to more rapid proliferation than for somatic cells [27]. The maintenance of pluripotency depends on the stable inter-regulation of pluripotency transcription factors (PTFs) [28], mainly by the genes OCT4, SOX2 and NANOG [29]. Fluctuations of the PTF abundances are believed to cause the variation in pluripotency in different subpopulations [28]. Destabilisation and the interaction of these PTFs with chemical signalling pathways triggers Scales of hPSC behaviour: a intra-cellular scale e.g., cell cycle, division and inheritance of pluripotency factors. b Cell micro-environment e.g., interaction with other cells, the medium and substrate. c Colony-scale phenomena e.g., patterning of differentiated cells differentiation, the departure from the pluripotent state [28,30] towards specific cell fates [31]. The cell cycle also affects pluripotency and cell fate [32] and vice versa [29,33,34]. Moreover, recent work suggests that the PTFs are inherited asymmetrically as a cell divides [35], biasing the fate of the daughter cells and contributing to colony heterogeneity.

Micro-environment
As in the embryo, the local environment of the cell is key to its in-vitro evolution. One of the leading environmental factors affecting hPSCs is the substrate on which they are grown. Substrates may either consist of a layer of mouse or human 'feeder' cells or a protein substrate, with the latter growing in popularity for clinical application since they avoid the risk of genetic contamination. The substrate influences pluripotency [36] and mobility [37] through its growth factors and adhesion forces. Low cell motility improves clonality by suppressing cross-contamination of colonies [38], although its role in colony heterogeneity is yet to be established.
As well as the substrate, cell-cell interactions are also important. hPSCs benefit from being in colonies where they exhibit higher viability and pluripotency [39]. hPSCs apparently sense each other up to a distance of around 150 μm (of order 5 cell diameters) [40,41]. Meanwhile, as the colony grows and becomes denser, the mutual mechanical pressure of the hPSCs can affect the cell cycle [42].

Colony scale
Perhaps most intriguing, yet least understood, are behaviours that emerge on a colony scale. The promotion of pluripotency in larger colonies [43,44] shows that single cells are influenced by the whole colony. Indeed, it has been suggested that pluripotency is a collective statistical property of cells [45], rather than a well-defined property of individual cells.
Further colony-scale effects are evident in the spatial patterning of the cell fates after differentiation. Mechanical forces and chemical signals operating over distances larger than the cell separation influences single-cell genetic expression to form bands of differentiated cells [46] (illustrated in Fig. 2); these structures are enhanced under imposed boundaries, emphasizing the role of mechanical forces [47,48]. With further understanding, mechanical effects and boundaries could be harnessed to engineer specific desired differentiated cells [49].
Incorporating these complex behaviours over multiple scales into mathematical models is challenging. A key goal is to develop coherent models which capture the individual cell behaviours, e.g., cell kinematics and the inter-cellular maintenance of pluripotency, and lead to the observed collective effects on the colony scale, e.g., collective migration and the spatial patterning of pluripotency and differentiation.

Cell kinematics
Motility is an intrinsic property of hPSCs; they can increase their migratory activity under certain conditions [50]. Their migration is achieved through adaptations in cell morphology via the reorganisation of the actin cytoskeleton to form a leading edge pseudopodia [51]. Unregulated cell migration in-vitro can cause clonality loss as the cell population grows which is undesirable when a genetically identical clonal population is required [52,53]. Furthermore, anomalous cell migration has been linked to deviations in the undifferentiated state of hiPSCs [54]. A thorough understanding of the migration of hESCs is needed to optimise in-vitro clonality and facilitate the development of therapies for migration related disorders. Here we discuss the kinematics of isolated cells and their pairs as well as cell migration within colonies.

Kinematics of isolated cells and pairs
hPSCs are often seeded at low density to preserve the clonal purity of the emerging colonies. Migration of individual cells between the incipient colonies can result in clonality loss. It is important therefore to quantify the migration of individual cells upon a growth plate.
The unconstrained motion of cells on a 2D plane can often be described as a 2D random walk, the simplest being Brownian motion [55,56]. Random walks can be biased by an external source giving preference to movement in a particular direction (a biased random walk or BRW). A correlated random walk (CRW) involves a correlation in the direction of the next step in relation to the previous step, i.e., persistence, where the next step is more likely to be in the direction of the previous step, or antipersistence, where the next step is more likely to be in the opposite direction. CRWs often occur in cell kinematics in the absence of external biases [57][58][59].
The diffusive nature of a random walk can be quantified by considering the mean square displacement (MSD) of cell trajectories. The MSD is a measure of the trajectory of a particle from its starting position over time, ⟨r 2 ⟩ = ⟨( − ) 2 ⟩ , where (t) is the position of the particle, 0 is the initial position at t = 0 and angular brackets denote the average taken over all trajectories. For a typical diffusive particle, the MSD increases linearly with time, ⟨r 2 ⟩ ∝ Dt , where D is the diffusion coefficient. The root mean square displacement is given by ⟨r 2 ⟩ 1∕2 = √ 2Dt , from which D can be calculated. If ⟨r 2 ⟩ ∝ Dt , with < 1 the motion is sub-diffusive or super-diffusive with > 1.
The nature of individual cell movements has been observed through direct experiments with hPSCs (in particular hESCs) and analysed within the random walk framework [40,41,60]. The movement of single hESCs has been described as an isotropic random walk when the cells are in isolation, i.e., more than approximately 150 μm away from any neighbouring cells. As the separation distance decreases the cell movements become more directed towards each other, with motility-induced re-aggregation occurring in 70% of instances when the distance been two hESCs is less than 6.4 μm [40]. A minority of isolated single cells exhibit super-diffusive behaviour, contributing heavily to the motility related clonality loss [8,40,60]. Example experimental trajectories for cells exhibiting typical diffusive behaviour and super-diffusive migration are shown in Fig. 3. These results show that individual cell movement influences hPSC clonality, although the biological causes of the distinct diffusive behaviours remains to be explored with further experiments. They can provide additional guidance for improvement of clonogenic assays in the analysis of hPSC self-renewal [40] and be used to identify timescales for motility-driven cross-contamination between colonies which is of practical use when producing high clonality colonies.
Our study containing further experimental analysis of hESCs [41] has shown evidence of correlated random walks of individual isolated stem cells. Single hESCs (more than 150 μm away from any neighbouring cells, as in [40]) tend to perform a locally anisotropic walk, moving backwards and forwards along a preferred local direction correlated over a time scale of around 50 min, becoming more persistent over time. The motion is also aligned with the axis of cell elongation (Fig. 3) which could suggest an attempt to locate other neighbouring cells. Further experiments should quantify how the presence of multiple neighbours affects this anisotropic movement.
Our study also found that pairs of hESCs in close proximity tend to move in the same direction, with the average separation of 70 μm or less and a correlation length (the length scale of communication) of around 25 μm. Often the pairs of cells remained connected by their pseudopodia, even at larger distances (> 100 μm) when they exhibited independent movements. For the correlated pairs, it is not known whether the movement correlation is facilitated by the physical connection or the coordination is due to cell-cell chemical contact alone.
This quantification of cell motility allows the direct comparison of cell movement between cell types and under different experimental conditions. For example, the addition of a CellTracer (a common biological marker allowing the tracking of cell generations) results in significantly reduced migratory behaviour for individual cells [60].
There is evidence that cell migration in 3D does not follow a persistent random walk and new models will need to be developed to accurately describe this motion [61]. These experimental results further inform the development of individual based models for cell migration as a random walk and can be integrated into more complex models of cell movement within colonies in-vitro.

Colony kinematics
Stem cells also exhibit motion as part of larger groups and colonies. The coordinated migration of large numbers of hPSCs in-vivo is essential in tissue generation [62] and wound healing [63]. The modelling of such larger groups and colonies of hPSCs is more complex, as both collective and individual behavioural effects are involved [64].
Popular agent-based models have been developed to incorporate these results into colony models, but the challenges still remain to fully capture the experimental behaviours, especially collective aspects and cell migration in 3D. These agent-based migration models are often combined with models of colony growth and proliferation [65,66].  [41] hPSCs show coordinated intra-colony movements which cease upon differentiation [67]. Cell movement speed varies within colonies, with higher average speed at the periphery and lower in the central region [54]. Recently, a two-dimensional individual-based stochastic model was developed of cell migration, cell-cell connections and cell-substrate connections and captures well these experimental observations [65]. The model introduces the energies of cell-cell and cell-substrate connections. Any energy released by breaking and forming these connections allows cell migration to one of the eight directions on a square lattice. The direction of movement is determined at random based on a probability related to the cell's energy and a spatial weighting which favours a side rather than a diagonal direction (as described in [68]). Cell proliferation and quiescence (the reversible state of a cell in which it does not divide) are also included. The observation of the spatial difference in average movement rates is difficult to explain by experimental results alone. This computational model suggests that cell division is a leading factor in the increased mobility at the colony edges and will be useful for studying further behaviours of hiPSCs and improving bio-processing experiments.
So far we have considered the movement of cells in 2D, analogous to the common experimental practice of growing cells on a flat substrate. However, culturing cells in 3D is becoming prevalent in order to provide a more realistic representation of the in-vivo behaviours of hPSCs [69] and for modelling in-vitro engineering of tissues on 3D scaffolds [70].
There are recently developed models which provide a good starting point for a 3D simulation, such as the Physi-Cell model, originally developed for cancer cells but transferable to other cell types including hPSCs [71]. The model implements cell movement by defining a persistence time, a migration speed and a migration bias, allowing for a range of cell motions from Brownian to deterministic. Movement due to the mechanical interactions between neighbouring cells is also included.
Modelling 3D cell movement on a discrete lattice is widely used, e.g., for mesenchymal stem cell tissue differentiation [72] and cancer stem cell driven tumour growth [73]. Some models allow many lattice nodes per cell as in the Potts model [74]. There is also a range of 3D agentbased continuous models where cell movement is not restricted to a grid but a cell can move continuously in any direction as illustrated in 2D in Fig. 4 [75,76]. Here the movement is described using forces or potentials with positions obtained from differential equations of motion for each cell. In centre based models (CBM), each cell is represented by a simple geometrical object, such as a circle, whereas in vertex models a cell is defined by a number of connected nodes [77]. These models will be discussed in more detail in Sect. 4.
There are also models which focus on the cells' changing morphology. For example, a model has been developed for mesenchymal stem cells which includes the random formation, elongation and retraction of pseudopodia, resulting in dragging forces which lead to cell movement [66]. This model can quantitatively reproduce the spatiotemporal organisation of cells and emphasises the importance of cell-cell interactions in tissue formation. However, the model of Ref. [66] shows more ballistic and accelerated dynamics than experimental results [78]. How much of this discrepancy is due to differences in cell type and culturing conditions should be investigated further to clarify the model's applicability to different experiments.
Informed by experimental results, these general cellular computational models could be adapted to describe the 3D movement of hPSCs both in culture and in-vivo.

Colony growth
Colonies of hPSCs are formed by repeated mitosis in which two genetically identical daughter cells are produced from the division of the mother cell. The cell cycle is the sequence of events that occur in a cell in preparation for the division as described in Sect. 2.1. The simplest mathematical models incorporate cell proliferation probabilistically, with the division time for each cell drawn at random from a suitable probability distribution [65]. Others go a step further by moving cells through each cell cycle phase according to timings based on experimental data [79] or as cell volume increases [66]. Sometimes divisions do not occur; this probabilistic nature of self-renewal can be incorporated when the end of the cell cycle is reached [80]. There are also more complex models which describe the relationship between inter-cellular processes based on growth factors (proteins that regulate cell growth) [76] and Fig. 4 The migration of cells can be modelled either on a lattice or in continuous space more sophisticated mathematical models describing the cell cycle in terms of limit cycles [81].
The doubling time of stem cells number varies and can be affected by various environmental and chemical factors, including cell density and the colony maturity [8,[82][83][84]. Models of colony growth can be dynamical-system type models that address the time evolution of the colony size, or spatial models which track individual cells and the growing colony in space and time.

Population dynamics models
Population models have been used to understand the process by which blood cells are formed [22], cancer tumours grow [85] and the impact of hPSC colony growth on clonality [86]. Early population dynamics models for stem cells were based on stochastic birth-death processes [17] involving systems of ordinary differential equations [87]. One of the most popular models for hPSCs includes two populations of dividing and non-dividing cells, with a term for accounting for cell loss through death or differentiation (often referred to as the Deasy model, which is a development of the Sherley model to include cell loss) [88,89]. The evolving number of cells over time N(t) is obtained as where N 0 is the initial number of cells, is the mitotic fraction, D t is the cell division time, and M is the number of lost cells.
More recently, hyperbolastic growth models (a new class of parameter model for self-limited growth behaviours [90]) have been introduced for both adult and embryonic stem cells [91]. These growth models provide more flexibility in the growth rate as the population reaches its carrying capacity and have been demonstrated to capture experimental data well [90,91]. The population in this case is governed by a non-linear differential equation with the initial condition N(0) = N 0 , and the parameters L (representing the limiting value, or carrying capacity of the population), (the intrinsic growth rate), (a dimensionless allometric constant) and (additional term allowing for the variation in the growth rate). This model can be used to describe both proliferation and cell death rates more accurately than Eq. (1) [91] and helps identify when the growth of cells becomes self-limiting, a biological problem currently not fully understood. (1) Our most recent work develops a population model of the growth for hESC colonies based on experimental data [86]. We analysed the evolution of the colony populations and found that the distribution of colony sizes was multimodal, corresponding to colonies formed from a single cell and colonies formed from pairs of cells as shown in Fig. 5. This importantly shows inherent differences in the biological behaviours of cells with different numbers of neighbours. The colony populations can be described using a stochastic exponential growth model, with the growth rates of colonies emerging from single cell and cell pairs being drawn from normal distributions: with A = 0.039 and 2 A = 0.006 2 , B = 0.043 , 2 B = 0.002 2 , = 0.77 and = 0.23 inferred from the fitting to the experimental data shown in Fig. 5. The growth rate for colonies emerging from pairs of cells is greater than for colonies founded by single cells. This means that colonies that have grown from cell pairs are larger not only due to the initial condition but also because their proliferation rate is larger. This is consistent with observations that hPSCs proliferate more effectively when in close proximity to other cells [39,92]. This difference is important when the clonality of a colony needs to be assessed non-invasively, e.g., from its size.
Upon collection of further experimental data, the model can be expanded to describe colony growth from larger groups of founder cells. It is expected that the growth rate for colonies will increase with number of starting cells before reaching its peak and this should be quantified. These growth rates are also expected to vary under different experimental conditions. The model can be used to predict hPSC colony growth and to calculate the time scales over which colony size The growth rate probability distributions for both populations. Adapted from [86] no longer predicts the number of founding cells based on their seeding density. Up to this critical time, colony size can be used as a non-invasive marker of clonal single founder cell colonies. This model can also be used to simulate colony growth in space which is discussed in the next section.

Spatial modelling
Colony growth can also be modelled spatially and, as with cell migration, the models can either be set on a regular or irregular lattice or in continuous space. Each cell can be modelled individually in an agent-based model, or for large numbers of cells where agent-based models become computationally challenging, using continuum models. A thorough summary of these different model types, along with their advantages and disadvantages with a view to tissue mechanics is provided in [77]. Here the recent attempts to model hPSC colonies using a variety of these techniques will be discussed.
Our multi-population model, Eq. (3), can be implemented to explore the impact of colony growth on clonality [86]. Generating homogeneous populations of clonal cells is of great importance [52,53] as clonally derived stem cell lines maintain pluripotency and proliferative potential for prolonged periods [93]. To achieve this, cross-contamination and merger of colonies (illustrated in Fig. 6a) should be avoided.
Assuming that, initially, the cells are randomly scattered in a growth area with a particular seeding density (the average number of cells per unit area), each cell (or group of cells) proliferates according to Eq. (3). Each colony is then approximated by a circle, with a certain position in space (the geometric centre of the founding cells) and a radius based on the population size and an assumed cell area of 250 μm 2 [94]. The time at which a colony begins to merge with its neighbour, , is the time at which the perfect clonality is lost as illustrated in Fig. 6. Using the simulated values of and a least squares fitting leads to the equation = (−0.007 ± 0.0001)n 0 + (102 ± 3) with R 2 = 0.99 , in hours and n 0 in cells/cm 2 . We are therefore able to estimate the time taken for the first colony merge to occur from the simplified version of the fitting equation, where n 0 is the initial seeding density of cells before their attachment to the substrate in cells/cm 2 and is produced in hours. These results can be used to achieve the best outcome for homogeneous colony growth in-vitro by choosing the optimal cell seeding density.
Other spatial models consider each individual cell's position in space. Common vertex based models for adult stem cell proliferation use Voronoi tessellation to describe cell position and areas. The colony area is divided so that the area occupied by a cell is obtained by tracing straight lines between the position of a cell and all its neighbours and drawing a perpendicular line in the middle as shown in Fig. 7a. These lines form a convex polyhedron called the Voronoi cell. The Voronoi cells are not uniform in shape and their number of sides varies. The tessellation can be constructed from experimental images using the cell centroid or cell nuclei positions, as shown in Fig. 7b [94]. Voronoi tessellation has been used to model adult stem cells in intestinal crypts in 2D [95,96] and is now being transferred to hESCs. The model uses an agent-based approximation in which each cell is represented as a Voronoi tessellation of the space [96,97]. The domain grows according to the pressure flow due to mitotic divisions in the colony. The dynamics between the cells are described by an elastic potential acting on each cell i as with k v and k c elastic constants, i the area of each cell, 0 the equilibrium area and i the initial positions of the cells, which do not necessarily correspond to the centroids denoted with 0i . The first term in the right hand side of Fig. 6 a An example of two colonies merging from experimental images. The two colonies, shown in blue and orange are beginning to merge at 5 days after seeding. The scale bar represents 100 μm. b Diagram illustrating initially seeded cells and the colonies at time , the first time at which the two growing colonies touch each other from a simulation of the cell seeding model. The orange cells are classed as a pair and grow accordingly faster. From [86] Eq. (5) tends to enforce uniform cell size and the second one gives the shape of the cells. Since the forces are conservative, applying the gradient operator to Eq. (5) and adding a drag force, the total force acting on each cell is obtained. The boundary of the colony is modelled using 'ghost cells' whose only function is to bound the domain. Figure 7c shows a simulated colony undergoing a cell division. Cells in the middle of the colony experience a higher pressure and show mitotic arrest, i.e. they do not divide.
Spatially modelling each individual cell in a colony in this way raises an important question about the physical process involved in cell division: how does the colony rearrange to make space for new cells? In Voronoi tessellation models [96,97] the cells re-accommodate themselves according to the potential from the neighbouring cells or the crypt walls. In most square or hexagonal lattice-based models, one daughter cell is placed in the same position as the mother cell while the other is put in a neighbouring position, chosen at random [98], isotropic mitosis. If there is no free position available next to the dividing cell, the neighbouring cells are re-arranged into other available free spaces stochastically until there is a free space next to the dividing cell [65] or, if this is impossible, mitosis is suppressed (quiescence) [72,99]. Further experimental timelapse image data is needed to clarify exactly how the new cells are placed in real colonies.
Proliferation also depends on spatial and environmental factors. There is evidence that high cell density reduces cell proliferation [42], which has been captured in a model showing preferential cell division at the colony edge [65]. Self-organisation of cells has also been observed, where the newly divided (smallest) cells cluster together in patches, separated from larger cells at the final stages of the cell cycle [94]. This segregation by cell size allows the interchange of neighbours as the colony grows and could directly influence cell-to-cell interactions and community effects.
All the models mentioned in this section (with the exception of [72]) consider the spatial formation of cells in 2D. The general 3D cellular PhysiCell [71] model describes cells with a volume which varies with the cell cycle, with daughter cells having half the volume of their parent cell and are placed accordingly by their parent cell position. A combination of this model and the more hPSC specific spatial models could be adapted to describe the structure of 3D colonies.
Spatial models of hPSCs become increasingly complex with colony size, and it is difficult to successfully incorporate many properties of colony growth along with any collective migratory effects. The question of how colonies re-arrange upon cell divisions requires more experimental investigation to elucidate the best models. The development of these models has already had an impact in understanding the growth of cancer tumours [100] and wound healing [101].

Cell pluripotency
Pluripotency is the defining characteristic of stem cells, often referred to as a cell's 'stemness' . It is hPSCs pluripotency that gives them the capability of differentiating into any type of specialised cell in the human body. However, hPSCs can undergo spontaneous differentiation which is undesirable for further experimental applications. Mathematical models of pluripotency are deepening our understanding of how pluripotency is regulated, leading to the optimisation and control of pluripotency in the laboratory.
The decision of a stem cell to remain pluripotent or to differentiate into a particular specialised cell is known as its fate decision. It is not known when a cell makes this decision. Even clonal cells under the same conditions make different fate decisions and it remains unclear how much fate choice is lead by inherited factors versus environmental Fig. 7 a Voronoi diagram illustrating how colony area is split into tessellated cells. b The Voronoi tessellation obtained from the centroid positions of cells in an experimental microscopic image [94]. c Voronoi tessellation to simulate a proliferating hESC colony. The cells divide and give rise to two daughter cells under suitable conditions, see highlighted cells outlined in yellow. The colour bar shows the elastic field in Eq. (5) factors and intracellular signalling [102]. There are several thorough reviews of the computational models of cell fate decisions [103][104][105]. Here we focus on the regulation of pluripotency and spatial patterning within colonies.
Biomedical and clinical applications of hPSC colonies demand tight control of colony pluripotency and homogeneity [43], yet this remains challenging. At a single-cell level, pluripotency is inherently stochastic; indeed, it has been proposed that pluripotency is only defined statistically within a population [45]. Cells are regulated by their local environment [54,106], notably their beneficial interactions with neighbours [44,46]. Colonies exhibit heterogeneous subpopulations of cells with differing levels of PTF expression [28,30] suggesting a play-off between disruptive single-cell and regulatory community effects. Such heterogeneity is undesirable, biasing evolution the trajectories and leading to spatially disordered differentiation [47]. Here we will consider intra-cellular models of pluripotency based on PTFs, and the spatial organisation of pluripotency at the colony level.

Fluctuating PTFs
The positive-feedback regulation between PTFs (the transcription factors which regulate pluripotency, see Sect. 2.1) was first described as a first order differential equation model using the Hill equations [107]. However, the parameters of such a model are difficult to estimate accurately [108]. More recently, PTFs have been modelled through branching processes [109]. A thorough review of the models of pluripotency is available [18], along with a review of computational modelling of the fate control of mouse embryonic stem cells, with many models transferable to hPSCs [105].
Recent experimental work has investigated how the PTFs vary over time, and how maternal PTFs are transmitted and distributed between the daughter cells [35]. The OCT4 abundance in the cells was tracked over time before and after the addition of an agent which induces differentiation (BMP4). The cell fates were also recorded. The OCT4 values over time for all cells, organised by cell fate (pluripotent, unknown or differentiated), are shown in Fig. 8a.
We are currently working on modelling the trends and fluctuations in pluripotency over time based on the experimental OCT4 data in [35]. First we quantified the nature of the persistence of the OCT4 time series. The Hurst exponent, H is a measure of the the long-term memory of a time series, with H = 0.5 corresponding to Brownian motion, 0 < H < 0.5 anti-persistence (a preference to change the direction of the last step) and 0.5 < H < 1 persistence (a preference to continue the trend of the last step). The mean Hurst exponent for the OCT4 data is 0.36, signifying anti-persistence and importantly suggesting self-regulation of pluripotency. We are exploring stochastic modelling techniques, particularly fractional Brownian motion to capture the anti-persistence and the stochastic logistic equation to model the evolutions of the cells. Both of these models are well established, however their application to modelling pluripotency is novel.
As the general OCT4 levels is inherited after cell division, pluripotency levels are most similar among closely related cells even when a reasonable level of randomness is allowed for [35]. The analysis in [35] also shows that OCT4 is not always equally allocated between daughter cells upon cell division with the split being sometimes asymmetric, as shown in Fig. 8b. Models of pluripotency inheritance should take into account this variation in the splitting ratio upon cell division. This study also suggests that a cell's decision to differentiate is largely determined before the differentiation stimulus is added and can be predicted by a cell's pre-existing OCT4 signalling patterns. These results imply that the choice between developmental cell fates can be largely predetermined at the time of cell birth through inheritance of a pluripotency factor [35]. Note that although the cell pre-stimulation status can influence the stimulation efficacy, the addition of BMP4 still favours the shift towards a differentiated phenotype.
It is worth noting that here only OCT4 is considered due to the availability of experimental data. For future work  [35]. b The OCT4 splitting ratio between daughter cells before and after BMP4 addition. Figure from [35] similar experiments should be conducted for NANOG and SOX2 to investigate the relationships between the three PTFs and any differences in their dynamics. This would then allow the development of a system of coupled equations to describe the PTF behaviour.
These results highlight the important properties for models of hPSC pluripotency to capture at the individual cell level: the stochastic inheritance of PTFs, the anti-persistence or self-regulation of pluripotency and the predetermined cell fate decision. Suitable models can then be developed to not only represent the behaviour on an individual cell scale, but also the colony scale.

Spatial organisation
Pluripotency also shows spatial variation on the colony scale. Preliminary experiments monitoring the OCT4 levels in colonies grown from single cells at 72 h post seeding show that pluripotency is clustered, with highly pluripotent cells grouped together, as shown in Fig. 9.
The differentiation of hPSCs also shows distinctive spatial patterning [46,47]. Experiments monitoring the pluripotency marker SOX2 and the differentiation marker AP2 have shown that differentiation occurs preferentially at the colony periphery in a band of constant width, independent of colony size, as illustrated schematically in Fig. 2c and shown in Fig. 10 [46]. These differentiated cells originate from the outer third of the colony, and remain at the edge. This provides important information for modelling the spatial patterning of the pluripotent state.
This within-colony spatial patterning behaviour of the differentiation has been captured by a mechanical bidomain model [110], a continuum model first developed to describe the elastic behaviour of the cardiac tissue [111]. The model predicts that differentiation and traction forces occur within a few length constants of the colonies edge, consistent with the experimental results for differentiation in hPSCs [46,47]. The model assumes that differences in displacement are responsible for any mechanotransduction (chemical processes through which cells sense and respond to mechanical stimuli) and describes both the intra and extra-cellular spaces in colonies with relationships between stress, strain and pressure forces. The basic equation for the difference between the intra and extra-cellular displacements for changing distance from the colony centre r, u r and w r respectively as where T is a uniform stress caused by the growth and crowding of cells, is the shear modulus, is a length constant and R is the colony radius. This model shows that if the difference between the intra-cellular and extracellular displacements drives the differentiation, then differentiation is confined to the edge of the colony. This model could be further developed to include more complicated geometries as currently the colony is assumed to be circular to allow analytical solutions to the model equations. Furthermore, it is worth investigating whether the cell growth represented by the tension T is a function of u r − w r alone, as observations for hESCs suggest distinct actin organization and greater myosin activity near the colony edge, implying that T could be non-uniform [46].
Further experiments are needed to collect data on the pluripotency of cells across colonies. Analysis of the data using techniques common in spatial statistics will allow Fig. 9 A microscopy image of a hESC colony at 72 h after seeding, alongside a colour-coded version of the same colony quantifying the level of expression of OCT4. Red represents the highest pluripotency with blue representing the lowest. Scale bar represents 50 μm Fig. 10 a Phase (top) and immunostaining images (bottom) of hESC colonies before and after BMP4 addition. b Analysis of expression of a pluripotency marker SOX2 and differentiation marker (AP2 ) 3 days after BMP4 treatment. Fluorescent intensity is plotted as a function of distance from the colony edge and normalized to the maximum intensity of each colony [n = 20 colonies, p < 0.0001 and represents statistics for AP2 (green) and SOX2 (red) levels between distance 35 μm and 175 μm from the edge using a twotailed paired t test]. Error bars represent standard deviations from the mean. Adapted from [46] the continued development of pluripotency models on the colony scale.

Discussion
Mathematical and computational models of hPSC growth are essential in formulating non-invasive predictive tools.
Although we have focussed on hPSCs here, it is worth noting that similar models are used to describe the reprogramming of somatic cells into iPSCs, which is still a lowyield process with the underlying processes of cell fate decision uncharacterised [6]. As the reprogramming is a stochastic process, most mathematical models in this area probabilistic [23]. A model describing cell types as a set of hierarchically related dynamical attractors representing cell cycles has lead to the identifications of two mechanisms for reprogramming in a two-level hierarchy: cyclespecific perturbations and a noise-induced switching [21]. These reprogramming protocols make specific predictions concerning reprogramming dynamics which are broadly in line with experimental findings. Another reprogramming model using a two-type continuous-time Markov process with a constant reprogramming rate has revealed two different modes of cellular reprogramming dynamics: TF expression alone leads to heterogeneous reprogramming while TFs plus certain other factors homogenise the dynamics [112].
Here we have discussed some key properties of hPSCs: cell kinematics, cell proliferation and cell pluripotency. However, there are other important factors which could be included in modelling, e.g., environmental factors, cell-cell signalling, intra-cellular properties and collective migration. Models isolating a few of these key properties have often captured experimental results well. For example, focussed migration models have lead to a greater understanding of the behaviour of isolated cells [40,41,60] and the movement of cells within colonies [65,66]. It is worth noting here that since it is not known what causes homogeneity in the motility characteristics of individual cells, the analysis is often statistical, considering average properties of the population as a whole.
There are many population models for colony proliferation, taking into account cell divisions and deaths, providing a distinct computational advantage over more complex spatio-temporal models. Models of colony growth have been used to investigate the impact of colony expansion on clonality [86], cell regeneration within intestinal crypts [95,96] and tumour growth [100].
Many current efforts focus on modelling cell pluripotency and cell fate, as applications of hPSCs require greater control over pluripotency and differentiation trajectories. The stochastic nature of pluripotency at the single cell level [45], along with regulatory community effects leads to heterogeneous sub-populations across colonies [28,30]. Recent studies of the fluctuations of PTFs throughout colonies [35] and spatial patterning of differentiation [46,47] are being used to inform the development of models of pluripotency and cell fate.
Developing comprehensive models of hPSCs remains challenging, due to their many complex properties across multiple scales, and not yet characterised collective behaviour effects. It is also difficult to match parameters with experimental observations. Model refinement should be based on a two-way interaction with experiments; model parameters should be informed by experimental results, and models should influence experimental design. Such models have already helped provide an insight into tissue formation, wound healing, tumour growth and the reprogramming of iPSCs and will no doubt continue to do so as these models progress.
Acknowledgements We acknowledge financial support from Newcastle University, and European Community (IMI-STEMBANCC, IMI-EBISC, ERC #614620), NC3R NC/CO16206/1) and BBSRC UK (BB/I020209/1) for providing financial support for this work. IN acknowledges the grant from the Russian Government Program for the recruitment of the leading scientists into Russian Institution of Higher Education 14.w03.31.0029. SOF thanks the National Council for Science and Technology (CONACYT), Mexico, for the scholarship CVU-174695. AS acknowledges partial financial support of the Leverhulme Trust (Grant RPG-2014-427).

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.