Correlated random walks of human embryonic stem cell in-vitro

We perform a detailed analysis of the migratory motion of human embryonic stem cells in two-dimensions, both when isolated and in close proximity to another cell, recorded with time-lapse microscopic imaging. We show that isolated cells tend to perform an unusual locally anisotropic walk, moving backwards and forwards along a preferred local direction correlated over a timescale of around 50 minutes and aligned with the axis of the cell elongation. Increasing elongation of the cell shape is associated with increased instantaneous migration speed. We also show that two cells in close proximity tend to move in the same direction, with the average separation of 70 um or less and the correlation length of around 25 um, a typical cell diameter. These results can be used as a basis for the mathematical modelling of the formation of clonal hESC colonies.


Introduction
There are many different types of active and spontaneous cell motion, e.g., swimming, gliding, crawling and swarming, detected in both prokaryotic and eukaryotic cells [1,2]. The favour of one mechanism over another depends on the environment and the balance of achieved displacement and energy expenditure. Cell motility and migration is essential in many biological processes including the development, morphogenesis and regeneration of multicellular organisms, wound healing, tissue repair and angiogenesis [3,4,5,6,7,8]. Anomalous cell migration can cause developmental abnormalities, tumour growth, neuronal migration disorders and the progression of metastatic cancer [9,10,11].
Unconstrained cell migration on a plane in vitro can often be described as a twodimensional random walk [12]. The simplest random walk, the Brownian motion, is uncorrelated (the current direction of movement is independent of the last) and unbiased (the direction of each step is random). Correlated random walks (CRWs) involve a directional bias; there is a preference for the direction of the next step to be related to that in the previous step. It is common for cells in the absence of external biases to migrate as CRWs: the migration of amoeboids [13], mammary epithelial cells [14] and mouse fibroblasts [15] have all been modelled as CRWs.
Adaptations in cell morphology facilitate migration. Some eukaryotic cells achieve motion through the coordinated and cyclic reorganisation of the actin cytoskeleton, which determines their speed, direction and trajectory [16]. Several types of protrusive pseudopodia structures have been characterised, which mainly differ in the organisation of actin [17]. Analysis of the formation of pseudopods has shown that cells extending pseudopodia which then split into two to allow a change of direction exhibit strong persistence and small turning angles [18,19].
Understanding the form of cell trajectories provides important insights into diverse cell motility modes and helps to design and interpret experiments. For example, understanding the role of cell migration in metastatic cancer has led to new treatments which modify signalling pathways and alter cell morphology to reduce cell motility [20,21]. A thorough understanding of the mechanisms underlying cell migration will not only deepen our understanding of many integral biological processes but also facilitate the development of therapies for treating migration-related disorders.
In this work we analyse the migration of human embryonic stem cells (hESCs) in vitro on a homogenous two-dimensional matrix. Due to the promises of clinical applications of hESCs and hiPSCs (human induced pluripotent stem cells) and the discovery of new engineered substrates for cell growth, data presented in this paper is of prime importance [22,23,24]. Surface-engineered substrates provide an attractive cell culture platform for the production of clinically relevant factor-free reprogrammed cells from patient tissue samples and facilitate the definition of standardised scale-up methods for disease modelling and cell therapeutic applications [25]. Because the clinical application of stem cells may require as many as 10 1 0 cells per patient and disease modelling efforts typically require more than 10 6 cells to make a single differentiated cell type, robust methods of producing cells under conditions that accelerate proliferation could be particularly valuable [26,27]. Feeder-free systems represent key progress in simplifying hESC/hiPSC production but most of these systems (synthetic polymers, peptide-modified surfaces, embryonic extra-cellular matrix (ECM) laminin isoforms, fibronectin from ECM with a small molecule mixture, and various vitronectin proteins) provide only modest gains in scaling-up hESC/hiPSC production because they still require seeding at a suitably high cell density and passaging through multicellular clumps. Even for the defined systems that support clonal growth [28,29,30], mass production of synthetic polymers, recombinant proteins, or small molecule mixtures may be a challenge, particularly when considering the number of cells needed for disease modelling and clinical application. These matters highlight the need to understand the factors that may facilitate clonal expansion of hESCs and hiPSCs for clinical needs. Current efforts are focused on optimising differentiation protocols in order to generate homogenous populations of cells of interest, hence an understanding of the features that charcaterise the starting cell population is essential for informing these protocols.
Unfortunately, the motion and dynamics of single and pairs of hESCs/hiPSCs has received limited attention. In culture, hESCs are anchorage-dependent and migrate through actin cytoskeleton reorganisation [31]. The main structures that define the leading edge on a migrating hESC are referred to as pseuodpodia. Motility is an intrinsic property of hESCs, and they perform an unbiased random walk when they are farther than 150 µm apart with cells closer to one another exhibiting coordinated motion [32]. Our previous work [33] investigated how the kinematics of single and pairs of hESCs impact colony formation. We performed statistical analysis on cell mobility characteristics (speed, directionality, distance travelled and diffusivity) from the timelapse imaging. We demonstrated that single and pairs of hESCs migrate as a diffusive random walk for at least 7 hours of evolution. We showed that for the cell pairs mutual interactions of closely positioned cells strongly affect the migration, and we identify two distinct behavioural regimes for cells resulting from a division. Also, the cell pair as a whole is shown to undergo a random walk with characteristic diffusivity [33].
Here we focus on the migration of single and pairs of hESCs by examining more subtle and yet significant aspects of migration as a further step towards understanding cell group formation from a single cell. We consider how the direction of the motion is related to the cell morphology and analyse how the separation of cells affects their coordinated movements.

Methods
We follow the methods used to prepare and plate, and then image and track hESCs described in our previous work [33]. In brief, hESCs (WiCell, Madison WI) were plated at a density of 1500 cells/cm 2 onto 6-well plates pre-coated with Matrigel R Basement Membrane Matrix (Corning Inc.), in the presence of mTeSR TM 1 media (STEMCELL Technologies). ROCKi (10µM, Chemdea) was present for the first hours after plating, and removed before time-lapse imaging.
After 1 hour, the plates were imaged with time-lapse microscopy (Nikon Eclipse Ti-E microscope) images taken every 15 minutes over 66 hours at a resolution of 0.62 µm/pixel. From these images, we selected 26 single hESCs and 50 pairs of hESCs. Single (isolated) hESCs are defined as those that initially have no neighbour within a 150 µm radius; interactions of hESCs are negligible beyond this distance [32]. The lineage trees for these cells are provided in Ref. [33]. We define the time variable t as zero at the start of the image recording. The pairs of hESCs are those where the separation of two cells is less than 150 µm from each other and more than that from other cells. The cells either exist as pairs at the start of the imaging, or form a pair when a single isolated cell divides.
Each cell in our analysis was manually tracked throughout its motion, and its position in each image frame was defined as the location of its geometrical centre by eye, or 'centre of mass' if the mass within the cell density is considered constant. For the single cell considered in Section 3.2, the cell boundary and geometrical centre was tracked using ImageJ [34,35]. Comparison of this to the previous coordinates taken by eye showed no significant difference. Tracking of a single cell ceased when the cell died; cell pairs were tracked until one of them died or divided. We did not follow cell triples even when they were formed by division of a cell in a pair. Formation of a pair from convergence of two unrelated cells is rare since the individual random walks lead, on average, to the divergence of cell trajectories provided sufficient space is available.
The instantaneous velocity of a cell was obtained from its displacement between two consecutive frames. Circular statistics calculations were performed as described in Ref. [36] using Matlab and its Circular Statistics Toolbox (Directional Statistics) [37].

Results
Figure 1(a) shows images of one of the cells during its migration; its full trajectory is shown in Figure 1(b). This cell is elongated in the instantaneous direction of motion, with a pseudopodia protrusion leading its next movement. The relation between motion and morphology is discussed in Section 3.2. The single cell shape can vary between approximately circular, with diameter of around 20 µm, to more elongated with length of up to 70 µm. In comparison, hESCs in colonies tend to be circular and considerably smaller, with diameters typically about 10 µm [32,38].
The cells can, and often do, change their direction of motion by up to π. An example is shown in Figure 1(a). The cell moves in the direction of its persistent pseudopodia protrusion, before contracting and moving in the direction of a new pseudopodia, resulting in a change of direction by approximately π. The whole manoeuver in this example takes about 6 hours.
The lineage trees for the 26 single cells can be found in Ref. [33]. Death rates are low, with only two cells dying before dividing. The remaining cells have divided by t = 20 h, with division occurring at a mean time interval of t d = 7 h. The median speed of the cells is 16 µm/h, with the average of 23 µm/h and with no noticeable differences in the migration behaviour between single cells which eventually die or divide.

Single cells: correlated random walk
First we seek to test for a bias in the direction of the single cell movements. We measured the turning angle, that is, the change in direction of the cell from one time frame to the next, denoted θ and illustrated in Figure 2(a). As well as the turning angle with respect to the earlier direction of motion, we also considered the angle φ between the cell displacement and the global frame that does not change with time. Figure 2(b) shows the polar histogram of θ for 26 single cells, while Figure 2(c) presents the corresponding linear histogram. It is evident that the distribution has maxima at θ = 0 and θ = π: the cell preferentially moves directly forwards or directly backwards with a roughly equal frequency between the two directions. The bias is robust, remaining even if small steps (< 7µm) are removed from the dataset. The mean axis of movement, shown in Figure 2(b), is approximately along the θ = 0 or θ = π (with the standard deviation of σ θ = 0.19). In this manner, the motion represents a quasi-one-dimensional random walk. Both the χ 2 and V tests reject the null hypothesis that the probability density of the turning angle θ is uniform at the 99.5% confidence level. The probability density distribution can be approximated by p = a + b cos(2θ) with a = 0.16, b = 0.04 and the R 2 value of 0.55. This fit suggests a symmetric spread of the distribution about θ = 0 and θ = π.
The distribution of the turning angles has a distinct temporal pattern. Figure 3 shows the polar histograms of θ at early (0-5 h), intermediate (5-10 h) and late times (10-18 h). At early times, the distribution is slightly biased towards θ = π, indicating a weak dominance of the back-and-forth motion over a systematic forward motion. However, this effect is weak and the distribution is approximately uniform over angles. This is consistent with our previous observations that the motion of hESCs is close to an isotropic random walk at early times [33]. By late times, however, the distribution is strongly biased towards θ = 0, that is, persistent forward motion. What we see on average for all times is a mixture of persistent and back-and-forth motions. This feature can be characterised with the temporal autocorrelation function, C θ (τ ), for twohourly moving averages of the angle θ. For each cell, C θ (τ ) is calculated as the circular correlation for θ with itself, delayed by a time lag of τ . The average autocorrelation over all single cells, C θ (τ ), with least-squares fitting C θ (τ ) = e −τ /τc , τ c = 0.8, is shown in Figure 4(a). We see a temporal correation in θ, with an average correlation time of τ c = 0.8 h.
We find no significant correlation in the global direction of movement, φ, when individual steps are considered. We attribute this to the dominance of the backand-forth motion over short periods of time. However, considering two-hourly moving   averages we again find a systematic trend. The average autocorrelation over all single cells, C φ (τ ), is shown in Figure 4(b). We see a temporal correlation in φ, with leastsquares fitting C φ (τ ) = e τ /τc , τ c = 0.7 shown in Figure 4 and hence a correlation time of τ c = 0.7 h, similar to that found considering the change in direction θ. Notably there is significant anticorrelation in φ for the time interval [2 h< δt < 5 h], in agreement with the biased nature of the random walk. For comparison, unbiased random walk has no such anticorrelation.

Single cells: direction of motion and cell elongation
It is evident, from the images in Figure 1 in particular, that the direction of motion appears to be aligned with the elongation axis of the cell structure including its pseudopodia. A further example is shown in Figure 5. This is unsurprising as cell branching and elongation has been shown to be involved in cell motion and directional persistence, although it has not been fully quantified [39].
To analyse quantitatively the alignment of the direction of motion and the elongation of the cell we measure the alignment angle of the cell, α, with respect to a global reference frame. Consider R(α), the vector from the geometric centre to the boundary of the cell and R max corresponding to the maximum magnitude of R. The alignment angle α is defined as the angle between R max and the horizontal, as shown in Figure 5(b). The polar histograms of α and the direction of travel on the plate, φ, both in the same global reference frame, are shown in Figure 5(c). Their mean values are α = 0.79 ± 0.34 and φ = 0.72 ± 0.35. The difference is insignificant as the Watson-Williams and Kuiper's tests provide no evidence to reject the null hypothesis that α and φ are from the same distribution at the 99% confidence level.
The speed of migration, v, and the measure of elongation of the cell, R max /R min , where R max = max|R| and R min = min|R|, are shown as functions of time in Figure 6. The hourly moving averages of R max /R min and the cell speed v have a Pearson correlation coefficient of 0.53 suggesting a slight positive correlation between the elongation of the cell and its speed. Hourly moving averages of α and φ are shown in Figure 6(c). This shows that directed movement is in the direction of the pseudopodia and suggests that the cell moves faster when it is more elongated.

Pairs of cells
Wadkin et al. [33] considered the movement of cell pairs (two cells within 150 µm of each other at the start of imaging) as a whole and found that the motion of their geometric centre is approximated by an isotropic random walk for up to around 7 hours of their evolution, albeit with reduced motility compared to that of single cells [33]. The diffusivity is reduced from 80 µm 2 /h for single cells, to 60 µm 2 /h for pairs. In this section we look in greater detail at the dynamics of pairs of hESCs, in particular the correlations between the individual motions of a pair's cells.
For the cell pairs in the experiment, the mean separation at time t, r(t) = (δx(t)) 2 + (δy(t)) 2 , where δx and δy are the distances between two cells in the x and y directions respectively, varies with time as shown in Figure 7. By performing a least-squares fit of the functional form r = A − Be −t/C , for parameters A, B and C we obtain the line r = (68 ± 0.6) − (37 ± 3)e −t/(2±0.03) . The asymptotic nature of r indicates an optimal separation of pairs at around 70 µm.
To quantify the coordination between the movements of the two cells in a pair, we measure the smaller angle between their velocities, 0 < ψ < π, illustrated in Figure 8(a). If the cells travel in the same direction on the plate, then ψ = 0, and if they travel in opposite directions ψ = π; note that ψ = π does not distinguish between the two cells moving exactly towards each other or exactly apart. The polar histogram of ψ for all the pairs is shown in Figure 8(b), with the corresponding linear histogram in Figure 8(c). There is a bias in the distribution towards ψ = 0, confirmed by the χ 2 test which rejects the null hypothesis that the distribution is uniform at the 95 % level, i.e., there is a significant preference towards pair cells moving in the same direction. Example microscopy images of a pair that move in this way are shown in Figure 9.
Binning ψ according to the separation distance, r, between two cells shows that this bias primarily occurs at small separations as shown in Figure 10. The χ 2 test provides evidence to reject that each of the histograms in Figure 10 is uniform at the 95 % level. However, a measure of the skew is shown in the first moment, i.e., the arithmetic mean, ψ (as opposed to the circular mean). For a uniform distribution between 0 and π the arithmetic mean would be ψ = π/2 or 90 • . For the ψ distributions for r < 20 µm, between 20-50 µm, between 50-100 µm and r > 100 µm the arithmetic mean values are respectively, ψ =73 • , 79 • , 89 • and 88 • , indicating there is bias towards ψ = 0 at smaller separations. Pearson's moment coefficient of skewness, γ = E[(ψ − ψ) 3 ]/σ 3 ψ , also provides a measure of the asymmetry in the distributions. For a perfectly symmetrical distribution γ = 0, while for a distribution skewed towards lower values γ > 0 and for skew towards higher values γ < 0. For ψ where r < 20 µm γ = 0.40, for 20 < r < 50 µm γ = 0.24, for 50 < r < 100 µm γ = 0.04 and for r > 100 µm γ = −0.02, showing reducing skewness towards ψ = 0. The Kolmogorov-Smirnov test provides no evidence to reject the null hypothesis that the distributions for r < 20 µm and 20 < r < 50 µm are the same. Similarly for 50 < r < 100 µm and r > 100 µm. However the test rejects the null hypothesis that the two smaller separation distributions are the same as the two larger separation distributions. Calculating ψ with separations binned more frequently shows the length at which the movement is correlated. By performing a least-squares fit of the form ψ = 90(1 − e −(r+r 0 )/m ), for parameters r 0 and m, we obtain the line ψ = 90(1 − e −(r+23.0)/25.9) ) with an R 2 value of 0.6, shown in Figure 11. The characteristic length of the decay is therefore 26 µm, a typical cell diameter, suggesting   the pairs only exhibit correlated motion while they are in contact.
We also analysed the motion of the cells in the pair via the pair correlation function. This function was not found to be sensitive to the correlations between the cell motions, and was unable to distinguish the cell motions from from IRWs. This analysis is presented in the Appendix.

Discussion
In culture, hESCs are anchorage-dependent: they adhere to the surface and sense external cues by extending lamellipodia and filopodia, referred to in a general way as pseudopodia. For directed movement in response to external factors, cells acquire a defined front-rear polarity extending a protrusive structure at the leading edge before subsequently moving the cell body, and retracting the trailing edge [40]. The integration of negative and positive chemical feedback loops accounts for the oscillatory behaviour of pseudopodia, i.e. cycles of protrusion and retraction which result in cell movement [2]. Observations of single cell movement in two-dimensions cultures, in the absence of external cues, indicate a production of pseudopodia structures in random directions, a behaviour observed in other cell types [41].
Our results are summarised in Figure 12. The relative angle of movement, θ, characterises the dynamics of random walks further to the mean-square displacement [42]. Our results show that isolated single cells migrate in an unusual uni-directional walk, moving backwards and forwards along a preferred local axis, with cells becoming more persistent over time. Hence, the longest lived isolated cells show the strongest directional persistence. Broadly, there are a wide range of example cells that exhibit a preferential turning angle; those that can be modelled as a correlated random walk as previously discussed, e.g., [13,14,15]. There are also examples of a biomodal preference for turning angle, similar to the one we see for single hESCs [43,44]. The bias in the Figure 12. (a) Single hESCs preferentially move along their elongation axis, at speed higher for a stronger elongation. (b) Cells separated by 70 µm or less move in a coordinated manner, whereas a wider separation implies independent biased random walk [33].
walk is further shown in the temporal correlation in both the change in direction, and the direction of movement with a correlation time of around 0.8 h. The microscopy images in Figure 1 show the elongated morphology of the single cells, with movement in the direction of the leading pseudopodia, leading to this motion along a local axis.
These single cells demonstrate random migratory patterns, travel large distances and do not result in colony formation. Isolated cells seeded at low density display directional migration towards neighbours [38]. Perhaps in the absence of neighbours, as in this experiment, the cells employ the uni-directional walk along the local axis in an attempt to locate neighbours. Our quantitative analysis of a directed cell trajectory confirms the axis of cell motion is aligned with the elongation axis of the cell. Increased elongation is also linked to increased speed, corresponding to previous results suggesting that persistence in direction of motion is linked to increased speed as a universal rule for all types of cells [45].
An understanding of the migration of single hESCs is integral to colony growth at low-density platings. Their directed, super-diffusive migration can facilitate colony expansion at low-density platings by the finding and joining neighbours, however this re-aggregration is undesirable in experiments which require colonies originating from a single cell to achieve a homogenous clonal population [32,46].
For pairs of hESCs their separation over time increases exponentially before approaching an asymptote at a distance of 70 µm. This shows that, on average, 70 µm is the optimal separation for pairs of cells. There is a preference for the cells to move in the same direction as each other on the plate at small separations (< 70 µm). At these small separations it can be seen from the microscopy imaging that the cells are physically connected by their pseudopodia, as in Figure 9. This coordinated movement could be due to an external stimulus, but the connection of the cell bodies facilitates this motion. At separations greater than ≈ 70 µm the motion of each cell in a pair appears uncorrelated. Often there is still a connection between the cell bodies at these distances, but the cells move in independent directions whilst maintaining the connection, and as an isotropic random walk when considered as a whole entity [33]. Neighbouring cells are integral to colony formation as cell survival and cell divisions are highly correlated with the number of neighbouring cells [38].
Another ramification would be an exploration of the effects of stem cell markers, such as NANOG, OCT and KLF, on cell migration. These factors have been shown to affect the migration, invasion and colony formation of various cancer stem cells [47,48,49]. Effects of pluripotency markers on the migration and motility of single hESCs have not been explored. hESCs with NANOG overexpression form colonies efficiently even at very low seeding densities. Cell motility and colony formation affected by stem cell markers are subjects of our future work.
Further experiments need to verify the robustness of these results under different culture conditions. This additional information on low density plated cells will assist in the development of agent-based models, combining the motion of diffusive and superdiffusive cells with their biological states and cell-cell interactions. where circumflex denotes a unit vector, r 12 = |r 12 |, δ(r − r 12 )=1 if r < r 12 < r + δr and zero otherwise, and δr is the width of a bin. A positive correlation indicates that the cells tend to approach one another, whereas C(r) < 0 indicates that they systematically move apart. The cells in pairs with C(r) ≈ 0 move with little or no coordination. The pair correlation for all 50 pairs considered together is approximately zero due to the averaging of positive and negative correlations, see Figure A1. However, we can assess the average degree of correlation (positive or negative) by considering the magnitude of the correlation, |C(r)|. The absolute value of the correlation for all pairs, calculated by taking |v i ·r ij | in Eq. A.1 and is within errors to the equivalent for a random isotropic walk for both cells in the pair. A comparison of θ (the angle of movement for each individual cell), C(r) and |C(r)| for the experimental data and for a simulated IRW for both cells is shown in Figure A1. For an IRW with no correlation between cells in a pair, the expected value of |C(r)| is 2/π, resulting from E[| cos(θ)|] = 2/π. Figure A1. (a) The angle of movement relative to the last, θ, (b) the correlation C(r), (c) the absolute correlation |C(r)| and (d)v ·r for i) the experimental pairs and ii) a simulated IRW for two cells. |C(r)| is theoretically constant at 2/π for an uncorrelated IRW pair.