Alignment of cell division axes in directed epithelial cell migration

Cell division is an essential dynamic event in tissue remodeling during wound healing, cancer and embryogenesis. In collective migration, tensile stresses affect cell shape and polarity, hence, the orientation of the cell division axis is expected to depend on cellular flow patterns. Here, we study the degree of orientation of cell division axes in migrating and resting epithelial cell sheets. We use microstructured channels to create a defined scenario of directed cell invasion and compare this situation to resting but proliferating cell monolayers. In experiments, we find a strong alignment of the axis due to directed flow while resting sheets show very weak global order, but local flow gradients still correlate strongly with the cell division axis. We compare experimental results with a previously published mesoscopic particle based simulation model. Most of the observed effects are reproduced by the simulations.


Introduction
The positioning and orientation of the cell division axis is a key element during the development of multicellular organisms and to maintain the functionality of the tissue. Many examples of pattern formation have been found in the developing embryo depending on emerging anisotropies in cell division, such as the imaginal wing disc of Drosophila [1]. These observations imply that the understanding of cell division orientation is fundamental. Single cell experiments have shown that the orientation of the cell division plane is determined by many factors such as cell shape, cell polarity and the spatial distribution of the extracellular matrix [2][3][4]. The orientation of cell division within tissues, however, is less well studied. Madin-Darby canine kidney (MDCK) cells are a widely used model system to study growing motile tissues [5][6][7][8]. This epithelial cell line exhibits strong cell-cell junctions and grows as a polarized and continuous monolayer in culture [9,10]. Cell division usually occurs within the plane of the epithelium, producing two daughter cells that remain within the cell sheet [11]. Cells continue to divide at confluent level with decreasing division rates until a maximum density is reached [7]. At the same time, cell motility is reduced with increasing cell density, displaying a glasslike transition in cellular behavior [12]. In general, directed cell migration can be initiated by introducing a model wound in the cell sheet, creating a new, free surface [8]. Cells spontaneously start to migrate and proliferate in order to close the gap, displaying striking features like the formation of fingers [13] and leader-cells [8]. During collective migration of cell sheets, forces are exerted on the substrate, leading to the development of a global tensile stress within the layer [14]. Individual cells react to generated stresses by changing their shape to elongated architectures or rearranging within the tissue in order to locally reduce the stresses. The orientation of the mitotic spindle is predicted to self-organize during cell division due to mechanical forces. As early as 1884, Oscar Hertwig proposed a simple concept, known as the 'long axis rule', implying that the axis of the spindle lies in the direction of the longest axis of the cell, cutting this axis transversely [15]. As a consequence, the possible correlation between emerging forces during collective migration and cell division orientation is an interesting issue. Poujade et al investigated the orientation of cell divisions in the central region of expanding, stripe shaped MDCK monolayers, finding no preferred direction of the division axis [8]. At the rim of the expanding culture, however, finger formation promotes strong correlation of cell division axis and polarity with the main finger direction as proven by Reffay et al [6]. A majority of the motility phenomena were explained recently by a simple mechanical model of growth [16] combined with a local orientation mechanism of cellular motility that does not require any alignment force from the neighboring cells, or chemical gradients to determine the direction of cell migration. This simple model was able to reproduce the glassy transition for confluent monolayers, finger formation at invading fronts, and the tension in expanding cell sheets [17].
In this work, we study the orientation of cell division axes in flowing and resting epithelial monolayers and compare the experimental results to a computational model. We use microstructured channels to create a directional cell flow of migrating and proliferating cells. The distribution and orientation of the cell division events as a function of distance to the leading edge and local flow properties are analyzed. We demonstrate that cell divisions are significantly oriented by the emergent cell flow in contrast to the scenario of cells in resting, confluent monolayers within the same channels. Computer simulations based on the algorithm developed by Basan et al [17] are able to reproduce the bulk properties of the experimental results, without further adaptation of the parameters. Near boundaries, however, they display strong deviations. Both the experimental data as well as simulations suggest a tight coupling between the direction of collective cell migration and the orientation of the cell division axis.

Results
In order to study the influence of directed cell growth on the orientation of cell division axes, we confine expanding epithelial cell sheets in microchannels. Using micromolding in capillaries, the cell-repellent polymer poly(ethylene glycol)-dimethacrylate (PEG-DMA) is shaped by a structured polydimethylsiloxane (PDMS) mold and cross-linked by UV irradiation [18]. MDCK cells are grown to confluence in front of the mold. After removal of the PDMS, the cells start to migrate into the three-dimensional PEG-DMA channels of different widths ranging from 200 μm-300 μm, corresponding to approximately 15-20 cell diameters (see figure 1). The expansion of the cell sheets is followed by time lapse microscopy over 10 h-35 h. We use stably transfected MDCK cells expressing mcherry labeled 2B histones, resulting in bright and located fluorescence of cell nuclei. Thus, the cell density can be identified by automated counting of cell nuclei. The local velocity field is determined by particle image velocimetry (PIV) analysis.
We compare the experimental results to particle based computer simulations. The model, as described in [17], is used without any parameter modifications. In short, cells are represented by two point particles, which interact with particles of other cells like sticky, soft colloids. The particles of one cell repel each other and cell division is initiated when a critical distance is reached. Additionally, cells can be in a non-motile or motile state. By switching into the motile state with a constant rate, cells pull in a random direction with a constant motility force. The rate of turning non-motile is dependent on the alignment of the motility force with the cell velocity (alignment results in lower rate than non-alignment).

Flow and density profiles in invading cell sheets
We study the orientation of cell division axes in confined, growing MDCK monolayers experimentally. First, we characterize the flow profile of the invading sheets for different channel widths (see figure 1). The displacement of the leading edge is monitored by the analysis of the averaged intensity profile of the brightfield images. The leading cell front advances at a constant speed of about μ ± − (22 5) m h 1 , independent of the channel width and invasion length as shown in [19]. Additionally, the average velocity of the invading cell sheet is measured to be much higher near the front compared to far away as shown in figure 2(a). The flow component in y-direction of the channel is negligible. Perpendicular to the channel walls, a flow profile is observed which is neither parabolic nor a simple plug flow ( figure 2(b)). The local velocity is reduced by a factor of three from the center of the channel towards the walls, referring to finite slip boundary conditions. During the expansion of the cell colony, a density gradient develops along the channels with increasing cell densities from about 1 cell per μ 1000 m 2 at the front to 2.5 cells per μ 1000 m 2 roughly μ 400 m behind the front (see figure 2(d)). The simulations show a constant invasion velocity, independent of channel width and penetration length, in agreement with the experimental results. The velocity profiles along the channels display the same trends as observed in the experimental data, however, simulations yield an almost perfect plug flow profile across the channels. The density profile is also consistent between simulations and experiments (see figure 2(d)). The similarity in the density profile shape as well as the constant invasion velocity suggest these quantities as units of measure to rescale the simulation results to physical units (for detailed information see Materials and methods).
In order to get good statistics, we averaged over the different experiments of channel sizes between μ 200 m and μ 300 m. We chose a bin size of μ 50 m ( μ 42 m) along the x-or s-axis for the experiments (simulations) and averaged over the whole channel width, several individual experiments and time (60-210 frames per experiment). Cross channel profiles, like the velocity profile in figure 2(b), are not averaged over the channel width but instead binned in y-direction with a size of μ 8 m ( μ 8.4 m).

Orientation of the cell division axis
The question we asked is to what extent cell division events are affected by the emergent flow in the expanding cell sheets. The axes of individual cell divisions are determined by identifying the position of the two daughter cells directly after division. The angle between the division axis and the x-direction of the channel is denoted with θ as described in figure 1(c). The orientation of the cell division is represented by a unit vector d. Due to mirror symmetry, it is best described by the nematic order tensor which is often used in the context of liquid crystals. The indices i and j represent the spatial directions. The order of one orientation relative to another, however, can be described by a simple scalar quantity: We find a significant alignment of the cell division axis in the bulk of the invading cell sheet, displaying order in the range of = S 0.2 x d to 0.5 (see figure 3(b)). The bulk of the invading cell sheet is defined as the part of the layer that is not affected by border instabilities at the front. Petitjean et al [13] reported a velocity correlation length in fingers of approximately μ 200 m, thus we assume that cells μ 200 m behind the front act as bulk cells (see also figures 2(a) and (d)). Towards the leading edge, however, the order parameter drops to negative values, indicating an orientation of the axis perpendicular to the channel. This phenomenon of division axis alignment perpendicular to the expansion direction (i.e. in the plane of the front), can also be seen in the videos of the experiments (see supporting information). Cells tend to be more elongated at the front and, furthermore, show alignment with it, which could be due to surface tension. Computer simulations are able to quantitatively reproduce the orientation of cell division axes. We find deviations close to the boundaries of the cell sheet (at the leading edge and the channel walls) between experiment and model, however, if divisions close to boundaries are not taken into account, agreement is quantitative (see figure 3(b), note that order is dimensionless, and thus without adjustable parameters or units).
It is known from literature [14,20,21] that a tensile stress is built up in advancing cell layers. As no direct stress measurements are accessible from the experimental data, we focus instead on the strain rate tensor In fluids (and most often in complex fluids as well), the stress tensor is directly proportional to the strain rate tensor. This tensor gives the deformation rates in all directions at each point and is symmetric by construction. The eigenvector λ, belonging to the larger eigenvalue λ 1 , gives the direction of largest extensile flow (or least contractile flow) and will be called 'main axis'. Due to the similarity of a growing cell to a force dipole, we define the difference between the largest and smallest eigenvalue λ λ = − p 1 2 as the dipole strength. It reflects the asymmetry of the flow gradient. The sum of both eigenvalues equals the divergence of the flow. Figure 2(a) shows a positive divergence behind the front, which is thus indicative of an expanding tissue (instead of a cell inflow from the channel entrance). The main axis of the strain rate tensor is an indication for the main axis of stress via the coupling between strain rate and stress tensor. We find no correlation between the overall average main axis of the strain rate tensor and the x-axis, while at the points of division the main axis is weakly correlated to the x-axis (see figure 3(c)). Similar to the division orientation, λ shows perpendicular alignment directly at the front. The simulations agree well in the bulk, also showing isotropic order. However, near the front, they display a conflicting behavior compared to the experimental data.

Boundary effects on cell division orientation
The orientation of the cell division axis shows a strong alignment with the main direction of the flow, however, we can not exclude an extensive effect of the channels wall influencing cell division events. In order to separate influences on the cell division axis generated by cell flow from boundary effects, we generated a closer-to-equilibrium scenario of non-flowing tissues within identical PEG-DMA channels. Confined cell monolayers are prepared by uniform cell seeding over the whole microstructures. Cells sediment and adhere solely inside the threedimensional channels due to the cell-repellent properties of PEG-DMA (figure 4(a)). A confluent and resting MDCK monolayer is formed and analyzed in terms of cell division events, cell density effects and the local velocity fields.
The order parameter S x d indicates a nearly isotropic distribution of the cell division axis orientation for the bulk of the channels ( figure 4(b)). Cell division events close to the PEG-DMA walls, however, show almost perfect order and are strongly aligned with the confinement. The profile shape is best described by an exponential and a decay length of ξ μ = ± (6.7 1.3) m. As the decay length shows, the alignment declines over roughly one cell layer and is well described by the assumption of total randomness in cell division orientation under constraint. In order to verify this assumption and compare it to the experimental data, we randomly draw center of mass positions and division angles from uniform distributions with a fixed daughter distance of μ = l 16 m. We define a division axis as the line connecting the positions of two daughter cells. All events, where the resulting division axis penetrates the wall, are discarded. Therefore, the range of possible orientations at distances smaller than the average distance between daughter cells decreases, i.e. leading to an increased alignment of the division axis with the wall. The resulting divisional order is shown in figure 4(b) and demonstrates a very good agreement with the experimental data. The observed alignment with the channel walls is also present in the case of invading cell sheets, although less significant as the bulk already shows a high degree of order (see supplements).
The particle based simulations reproduce both the random orientation of cell division axes in the bulk of the cell sheet for the confluently plated cells (figure see 4(b)) and the highly aligned division orientation for the invading cell sheets (see supplements), but, contrary to the experimental observations, show unoriented cell divisions close to the channel walls in both cases. Based on the observations in the resting cell monolayer, we can conclude that the observed alignment of the cell division axis in flowing cell sheets is mainly caused by the global orientation of cell flow.

Cell division orientation and local flow fields
So far, the divisional order has only been examined on large scales. Cells, however, have no known mechanism of sensing large scale variables like their distance from the front of the sheet. Thus, the order should only depend on local variables. In this section, we study the dependence of division axis alignment on measurable local variables. Due to broken Galilean invariance, symmetry allows for the orientation of the division axis by the velocity vector which would not be possible without an underlying substrate. Note, that while the degree of alignment with the x-axis increases with distance from the front, the xvelocity decreases. This contradicts assumptions of flow aligning divisions, as opposed to flow gradients. In fact, the average order between the local velocity vector and the division axis is only ± 0.16 0.014 for the flow experiments, and ± 0.06 0.014 for the no-flow experiments (errors are standard deviation of mean). We interpret the (small) positive correlation in the flow experiments as a common cause correlation as both the local velocity vector and the division axis are on average parallel to the x-axis. As discussed above, it is known that mechanical stresses orient cellular divisions. Moreover, this is a requirement in theoretical works explaining fluidization of tissues due to cellular division [22], which suggest to treat tissues as fluids on time scales greater than the division time. It thus stands to reason that cellular divisions should be oriented by local flow gradients, which we investigate in the remainder of the section. To study correlations between the local flow gradients and division orientation, we calculate the average order of divisions in direction of the flow gradients: where ϕ is the angle between the division axis, and the eigenvector corresponding to the largest eigenvalue of the strain rate tensor (see figure 5(b), inset). We find a clear correlation with an average value of = ± λ S 0.29 0.02 d (and variance σ = 0.67) between divisions and the strain rate tensor in the invading cell sheets. Indeed, the degree of alignment scales linearly with the difference of the largest and smallest eigenvalue of the strain rate tensor, the dipolar strength, which quantifies the asymmetry of the flow gradient. In confluently plated, resting cells, the degree of alignment to the local flow gradient is also present but smaller . Note, that in both cases, flow and confluency, the alignment between flow gradient and division axis is twice as large as the alignment with the local flow velocity.

Conclusion
In this article, we have studied the influence of an emerging flow in collective cell migration on the orientation of cell division axes. Expanding epithelial cell sheets were confined in microstructured PEG-DMA channels and flow characteristics were determined by PIV analysis.
The tracking of single cell division events was facilitated by stable nuclei staining of the MDCK cells, which additionally provides the basis for automated investigation of the cell density. We correlated flow fields, division events and cell density and managed to extract many quantitative features of tissue migration. The cell monolayers invaded the channels at a constant front speed with a well defined flow profile. A density gradient, with low densities at the leading edge and higher ones in the bulk, developed. We found a strong alignment of the orientation of the cell division axis with the direction of the observed flow.
In order to distinguish phenomena caused by the cell flow from boundary effects of the channel walls, we generated resting cell sheets within the same PEG-DMA microchannels. In this scenario, cell division axes were randomly oriented, supporting our suggestion that the emergent flow in migrating monolayers was the main influence on the direction of cell division.
Analyzing local quantities revealed that the orientation of the cell division axis correlates best with local flow gradients, as opposed to the channel axis or velocity direction. In fact, we found the average order between the division axis and the main axis of the strain rate tensor to be twice as large as the average order between the division axis and the local velocity direction. This is consistent with stresses orienting cell division, and the stress being proportional to the flow gradient.
Experimental results were compared to previously published, particle based computer simulations [16,17,22]. The model was originally conceived to study rheology and growth mechanics of three-dimensional tissues. This required an intrinsic and oriented feedback between growth and mechanics, which was achieved by representing single cells by two particles that repel each other in order to grow. This simple mechanical feedback was sufficient to describe growth of tissue spheroids [23,24] and expected viscoelastic properties [22]. Later, the model was extended by a simple motility algorithm, in order to understand particular properties of growing motile tissues: tension in growing colonies, glass-like arrest with density, swirling motion in bulk and fingering instabilities at the leading edge. These simulations were never designed to study divisional orientation. Thus, it is striking how many phenomena are captured well by the simulations. In particular, they displayed a similar decrease in cell density towards the leading edge, comparable emergent flow profiles, a weak correlation of the main axis of the strain rate tensor with the x-axis as well as a strong correlation with the division axis and strong correlations between the main axis of the strain rate tensor and the division axis. However, comparing quantities near boundaries (channel walls or leading edge of invasion front) revealed strong deviations and indicate points, where the model needs to be extended. As real cells have a complex inner structure and are highly deformable, we speculate that representing cells by more than two particles in the simulations is mandatory to resolve the dynamics close to interfaces.
The presented approach of measuring cellular flow and cell division events can be easily adopted to several adherent cell lines in order to prove the universality of our reported results. Naturally, the orientation mechanism should also work in three dimensions. While the simulations are inherently three-dimensional, motility experiments are much more challenging in three dimensions.
The orientation of cell division axes in flow fields could be of significant relevance for growing tissues in wound healing and development. It seems natural that an expanding tissue will orient individual cell divisions in the direction of expansion in order to reduce the principal stress. It remains to be seen if in vivo experiments show similar oriented cell divisions in cellular flows.

Experimental procedure
In this study, we analyzed the flow behavior and orientation of the cell division axis in confined MDCK cell sheets invading PEG-DMA microchannels. During the experimental procedure, we prepared the microstructured surfaces with standard cell culture dishes as substrates, followed by cell seeding (for better illustration, see figure 6). The PEG-DMA channels were shaped by the PDMS mold as a perfect negative replica and the cell-repellent polymer was cross-linked by radical polymerization, forming a three-dimensional, dense network. The microstructures were uncovered by the removal of the PDMS stencil, however, to protect the microstructured area from initial cell attachment, the PDMS mold remained on top of the structures during cell seeding. The cells were grown to a confluent level around the covered structures forming a proper edge at the border of the PDMS. After removal of the PDMS, the substrate consisted of free-standing PEG-DMA walls that guided the expanding cell sheets as both chemical and physical barriers. The free surface presented to the cells was the untreated petri-dish. The channels exhibited a height of approximately μ 15 m and a length of several millimeters, however, we only monitored the cell sheet expansion to a penetration depth of up to μ 700 m. In order to ensure the correct ascription of the observed effects to the emergent flow in collectively migrating MDCK cells as opposed to effects of the PEG-DMA confinement, we created resting cell sheets within the same kind of microchannels. In this case, the PDMS mold was removed before cell seeding to bare the channels to sedimenting cells.

Preparation of PEG-DMA microchannels.
Three-dimensional PEG-DMA channels were prepared according to the protocol described previously [18]. Briefly, a PDMS mold that features channels with 200 μm, 250 μm and 300 μm widths and approximately 15 μm height was created by curing the prepolymer against a master prepared by conventional soft lithography. Therefore, the prepolymer solution (Sylgard 184, Dow Corning, USA) was carefully mixed with the curing agent in a 10 : 1 ratio (w/w) and degassed for 15 min. The mixture was deposited on top of the master followed by an additional degassing step for 15 min to avoid bubble formation. Next, the polymer was cured at°70 C for 4 h. The molded stamp was peeled off and cut in a way so that the channels were protected by a small barrier of PDMS on one end of the structure since the other end of the channels were open. The PDMS was activated by argon plasma (Diener Electronic, Nagold, Germany) for 30 s. Immediately after the treatment, the stamp was placed upside down in close contact with a hydrophilic ibidi μ-dish (ibidi, Martinsried, Germany). A solution of PEG-DMA (Mn = 550) containing 2% (v/v) of the photoinitiator 2-hydroxy-2-methylpropiophenone (both purchased from Sigma-Aldrich, Germany) was freshly prepared. A small drop was deposited in front of the open ends of the channels and the structures were filled spontaneously by capillary force induced flow. The process was visually controlled under a microscope to ensure complete filling of the channels and the formation of proper edges. Once the PEG-DMA filled out the whole channel, it was polymerized by the use of an UV-ozone cleaning system (UVOH 150 LAB, FHR, Ottendorf, Germany) with a wavelength ranging from 185 nm-546 nm and a power of >50 mW − cm 2 for 10 min. The structures were stored in a drying oven (Binder GmbH, Tuttlingen, Germany) overnight at°50 C. Immediately before cell seeding, samples were sterilized with 80% ethanol for 5 min, washed three times with PBS and covered with cell medium while the PDMS mold still remained on top of the structures to protect the channels from cell adhesion.

4.1.2.
Experiments with global flow. Cells were trypsinized and centrifuged at 1000 rpm for 3 min. The cell pellet was resuspended in medium and cells were seeded in front of the PDMS stamp. Then, cells were grown to confluence overnight under culture conditions. The PDMS stamp was removed to bare the protected PEG-DMA channels and the sample was washed with fresh medium to eliminate dead cells. The cell layer started to expand towards the microchannels. Samples were stored under culture condition until the cell front reached the entrance of the channels (typically after 2 d). For time-lapse measurements, the medium is replaced with CO 2 independent Leibovitzʼs L15 medium, supplemented with 10% fetal bovine serum (FBS).

Measurements of non-flowing cell sheets.
To prepare confluent cell layers within the PEG-DMA microchannels, cells were seeded all over the uncovered structures. In this case, the PDMS mold was removed after storing the samples overnight in the drying oven. The bare channels were sterilized analogously to the protected ones with 80% ethanol for 5 min, followed by three washing steps with PBS. After cell trypsination, centrifugation and resuspension, cells were seeded all over the whole petri-dish. The non-fouling properties of PEG-DMA prevented cell adhesion on top of the channel walls, and consequentially, cells adhered exclusively inside the channels. The sample was stored under culture condition until confluence was reached. For time-lapse measurements, the medium is replaced with CO 2 independent Leibovitzʼs L15 medium, supplemented with 10% FBS.

Cell culture and microscopy
Doubly transfected MDCK cells with lifeact (EGFP) and H2B (mcherry) were kindly provided by Dr Wedlich-Söldner of the University of Münster. Cells were routinely cultured in Dulbeccos Modified Eagle Medium, supplemented with 10% FBS, 20 mmol L-Glutamine, high-glucose level (4.5 g l −1 and 110 mgl −1 pyruvate). Cells were kept at°37 C in a humidified atmosphere, 5% CO 2 level. Cell migration within the channels was followed by time-lapse microscopy by using an epifluorescence motorized Nikon Eclipse Ti microscope equipped with an ibidi heating stage (ibidi GmbH, Martinsried, Germany) and a CCD camera (model Clara E from Andor Technology). Image acquisitions were controlled through μ-Manager open source software. Brightfield and fluorescence images were taken at 10 min intervals using a 10× objective. Cells were maintained at°37 C using a temperaturecontrolled mounting frame. . The same applies to the angle ϕ between d and the eigenvector of the greatest eigenvalue of the velocity gradient tensor λ. These angles were then used to calculate the according order parameters S x d and λ S d .

Calculation of cell density.
For cell density determination, positions of all stained cell nuclei (x i ) were captured using the 'find maxima' function of ImageJ with a manually adjusted threshold for every movie. The relative error was determined to be below 6% by manual cell counting. In this two-dimensional system of a cell monolayer, the density ρ was defined as number of cells per area. Thus, the density was calculated from the cell nuclei positions x i for any area A as follows: 3. Detection of front displacement. In order to identify the leading edge of the growing cell sheet, we used the total intensity as a function of x (i.e. intensities of all pixels with the same y-coordinate are summed up). These spectra showed a sharp drop in intensity, corresponding to the position of the front. A custom Matlab script (The MathWorks, Natick, MA) was used to detect this sharp drop by calculating a relative 'coarse-grained' derivative. It turns out that the front position is usually well described by this derivative being above a certain threshold. In order to calculate this quantity, we first summed the intensity values of ten adjacent pixels and subtracted it from the sum over the ten adjacent intensity values shifted one pixel to the left from the original interval (meaning from the cell-free surface in the channel towards the advancing cell front). This difference was divided by the first of the two sums, and the quotient compared to a chosen threshold value. This process was repeated, shifting the intensity value intervals for both sums by one pixel to the left each time until the quotient exceeds the threshold for the first time, marking the point we defined as the front. The result was compared to the original movie and outliers were identified visually and deleted. The emerging gap was closed by linear interpolation. 4.3.4. PIV analysis. Velocity fields were mapped by PIV analysis using the MatPIV software package for Matlab (J Kristian Sveen: http://folk.uio.no/jks/matpiv/, GNU general public license), with some slight, custom modifications. The size of the interrogation window was 32 × 32 pixels, i.e. 21 × 21 μm, with 62.5% overlap in a single iteration. The resulting velocity vectors were filtered with a signal-to-noise ratio filter, a global histogram operator, a peak height filter and a local filter (all included in the MatPIV package) to produce smoother vector fields. Vectors removed by the filtering process were replaced by linear interpolation from the surrounding vectors, as long as there were at least five neighbors remaining that were not excluded by the previous step. Requiring this minimal number of neighbors prevented the analysis from interpolating into the region ahead of the leading edge of the cell sheet. 4.3.5. Calculation of velocity field derivations. Spatial derivatives were obtained by using a plain fitting method. A plane was fitted through all v x that lie within a certain radius r c of a point x i . The gradient of the resulting plane in x and y direction then corresponds to ∂ ∂ v x x and ∂ ∂ v y x , respectively. The same procedure was run for v y to get ∂ ∂ v x y and ∂ ∂ v y y . For a chosen radius r c that only includes the four nearest neighbors, the derivatives exactly matched the ones obtained by a finite differences method. By changing r c , the area over which the derivatives were averaged could easily be adjusted. Furthermore, this method allowed the definition of a second radius < r r s c of excluded points, in order to omit the velocities of the cell itself in the gradient calculation. The plane fitting was done with a linear least squares algorithm. This transformed the problem into a linear system of equations, which was solved by applying the LU decomposition algorithm implemented in the GNU Scientific Library. The choice of r c and r s , within reasonable limits, did not change the qualitative behavior of the results. In this work, we set r s = 11 μm and r c = 23 μm. This choice approximately calculates the derivatives on the scale of the cells neighbors.
From the derivatives, the strain rate tensor E ij (see equation 4) was calculated. The eigenvalues and corresponding eigenvectors of this symmetric 2 × 2 matrix were then obtained through the characteristic polynomial. Furthermore, with the unit vectors n parallel to d and ⊥ n perpendicular to d, the strain rate tensor E ij was used to calculate ∂ =

Computer simulations
Growth of motile tissues is modeled as presented in [17]. In brief, cells are represented by two point-particles, that repel each other with a growth force B. When a critical size is reached, the cells divide with a rate k d . The apoptosis rate was set to zero since experimental data showed an insignificant number of cell deaths. Mechanical interactions between the cells is achieved through a soft volume exclusion and an adhesion force, which basically resembles soft sticky spheres. Energy dissipation and random noise is introduced by a DPD thermostat that acts between the particles of one cell as well as between different cells. Cell motility is introduced as a simple two state model, in which each cell can be either in a motile or non-motile state. A nonmotile cell changes state with a rate k wake , in which case a uniformly distributed random motility direction is chosen. The transition rate from the motile into a non-motile state depends on whether the motility direction is aligned with the cell velocity ( + k ) or not ( > − + k k ). Motile cells pull against a background friction with a predefined motility force in the motility direction.
In order to resemble the experimental conditions as close as possible, we confined the cells in a rectangular box with a channel connected at one side (see figure 1(b)). The rectangular box acts as a reservoir and is initially covered by a confluent cell layer. We chose the reservoir dimensions as three times the channel width perpendicular to the channel and one times the channel width parallel to it. To approximate the finite slip boundaries in the experiments, reflective boundary conditions were used. The cells of the confluent simulations were initially distributed randomly over the whole channel with a predefined density of 2.3 cells per μ 1000 m . 2 In order to obtain the velocity gradient, we directly used the individual velocity of each particle as input to the plane fitting algorithm instead of calculating an average velocity on a grid, following the PIV analysis. 4.4.1. Units. A prominent problem of mesoscale simulations is their relation to real units. Computer simulations rely on contracting time and length scales, in order to make them feasible [25]. This naturally leads to discrepancies if one tries to compare all involved time and length scales. Thus, taking the model cells as one cell length for example is overinterpreting the model. Instead, one has to carefully choose how to relate simulation with experimental results. In order to reproduce mesoscopic continuum properties, we choose our rescaling quantities on a mesoscopic level. The unit of length is chosen such that the density of cells in the bulk matches. This gives us a length scale to convert simulation lengths to real lengths by μ = l 42 m 0 . For the unit of time we use the the fact that both experiments and simulations lead to a robust constant invasion velocity, independent of channel width and invasion length. With the length scale already defined, this gives us a timescale of = t 0.11 h 0 to convert simulation times into real time units.
The parameters used in the simulations are the same as in [17].

Error estimation
Unless otherwise noted, error bars are the standard deviation of mean, where averages over one channel are considered as one experiment. For example, the density profile in figure 2 is calculated for different distances from the front for several channels. First, the cell densities are determined and related to the actual front position in every frame. We obtain a mean density profile for each channel by averaging the densities of corresponding distances from the front over time. Subsequently, each channel's profile is considered as one measurement from which we obtain the mean and standard error of mean. We have no indication of systematic error of any kind.