Using a continuum model to decipher the mechanics of embryonic tissue spreading from time-lapse image sequences: An approximate Bayesian computation approach

Advanced imaging techniques generate large datasets capable of describing the structure and kinematics of tissue spreading in embryonic development, wound healing, and the progression of many diseases. These datasets can be integrated with mathematical models to infer biomechanical properties of the system, typically identifying an optimal set of parameters for an individual experiment. However, these methods offer little information on the robustness of the fit and are generally ill-suited for statistical tests of multiple experiments. To overcome this limitation and enable efficient use of large datasets in a rigorous experimental design, we use the approximate Bayesian computation rejection algorithm to construct probability density distributions that estimate model parameters for a defined theoretical model and set of experimental data. Here, we demonstrate this method with a 2D Eulerian continuum mechanical model of spreading embryonic tissue. The model is tightly integrated with quantitative image analysis of different sized embryonic tissue explants spreading on extracellular matrix (ECM) and is regulated by a small set of parameters including forces on the free edge, tissue stiffness, strength of cell-ECM adhesions, and active cell shape changes. We find statistically significant trends in key parameters that vary with initial size of the explant, e.g., for larger explants cell-ECM adhesion forces are weaker and free edge forces are stronger. Furthermore, we demonstrate that estimated parameters for one explant can be used to predict the behavior of other similarly sized explants. These predictive methods can be used to guide further experiments to better understand how collective cell migration is regulated during development.


Introduction
A fundamental understanding of the physical mechanisms driving morphogenesis during embryogenesis, during disease progression, and within biosynthetic engineered structures is progressing rapidly through a powerful combination of advanced microscopy tools and detailed biophysical models. With advances in microscopy, very large time-lapse datasets can be generated over the course of a few hours. Such time-lapse data contains a wealth of information describing both the structures involved in morphogenesis and their kinematics; this information can serve as the input for computational models that allow us to explore the biological and biophysical principles of morphogenesis and predict the behavior of cells and tissues under a variety of perturbations. This presents many technical challenges, from collecting quality images suitable for machine vision tools to developing computational models that implement appropriate biophysics and physiological rules at useful spatial and temporal scales. To maximize the benefits of these approaches, experimental designs must be able to integrate image data with model simulations in a way that allows robust statistical assessment of model fitness or failure.
We seek to design studies that integrate time-lapse image data with computational models in order to reveal differences in mechanical processes as experimental conditions are varied. A standard method of estimating physically relevant model parameters is through deterministic approaches such as the gradient descent method. Deterministic approaches identify a single set of optimal model parameters that offer a best fit to experimental data [32]. However, deterministic approaches do not quantify uncertainty in either the source data or the mathematical model, and thus are ill-suited for testing the significance of parameter fits between two experimental groups. Thus, we implement a Bayesian method that attempts to define a probability density distribution for the estimated model parameters. In particular, we implement the approximate Bayesian computation (ABC) rejection algorithm [33,34] to estimate biomechanical properties of spreading embryonic tissue via a mathematical model of collective cell migration. ABC has also been used in studying complex biological problems such as in population genetics, ecology, epidemiology, systems biology, and cell biology [33][34][35][36].
While our main goal is to develop a methodology for combining image data and computational modeling, here we apply this methodology to study the collective migration of tissue explants of various initial sizes. Epiboly is a key tissue movement during gastrulation in Xenopus laevis embryos, and it describes collective cell movements of the ectoderm of the animal cap region as is spreads and maintains coverage over the exterior of the embryo [37]. We make use of an in vitro model of epiboly where ectoderm from the late blastula stage animal cap region of the Xenopus embryo is microsurgically dissected and allowed to adhere to a fibronectin coated substrate. Dissected tissues, e.g. explants, spread over 10 to 24 hours in a manner similar to that observed of the tissue in vivo, suggesting that this spreading is at least partially driven by active processes in the ectoderm. Deep cells throughout the ectoderm are thought to contribute to epiboly through adhesion to and traction on a fibronectin coated substrate [38]. Cells at the margin of the explant are thought to generate outwardly directed traction forces that aid spreading. Additionally, as the tissue spreads, the cells of the ectoderm decrease in height but it is not known if this activity is active or merely due to passive strain in the tissue. Our goal in integrating time-lapse image data with a computational model is to determine the extent to which cell shape change, traction forces, adhesion, and tissue elasticity contribute to the mechanics of tissue spreading.
The physical mechanisms of Xenopus embryonic tissue spreading are investigated here by tracking both the movement of the tissue edge and internal tissue deformations. When cell nuclei can be tracked, internal deformations of spreading tissues may be tracked as local tissue density by counting cells in a subregion and then dividing by the area of that subregion [39][40][41]. Since tissues are generally opaque in the Xenopus ectoderm, we instead adopted an alternative automated method of extracting the local change in texture density between a pair of images in time-lapse sequences. This method involves calculating the strain, or local deformation, between a pair of images using elastic registration [42]. Elastic registration, a form of differential image correlation [43], fits well with an elastic continuum model of tissue spreading, and it can also serve as a reasonable approximation for more general models of cell migration, such as ones that incorporate constitutive properties such as viscosity or plasticity [44].
In this paper, we adapt an established continuum mechanical mathematical model from Arciero et al. [11] for collective cell migration in which the model parameters correspond to physical quantities.
We identify differences in physically meaningful parameters that regulate tissue spreading rates and deformations that depend on the initial explant size. Thousands of simulations of each experiment reveals multiple parameter regimes that lead to similar results. Single parameter sets identified using traditional deterministic methods may or may not represent these broader regimes and would be illsuited for statistical tests. To identify parameters that modulate the role of initial explant size on spreading, we use approximate Bayesian computation (ABC) to determine the probability density distribution for each parameter and identify statistically significant differences using Tukey's multiple comparisons test.

Results
In order to explore the mechanical processes that are involved in the collective migration of developmental tissues, we microsurgical isolated ectodermal tissue from the animal cap region of gastrulating Xenopus embryos. During gastrulation, this tissue spreads to cover the exterior of the embryo (Fig 1A). Once isolated and placed on a fibronectin-coated substrate, tissue explants spread essentially symmetrically outward (Fig 1B). Fibronectin is the principle extracellular matrix ligand for animal cap cells and the number of binding sites is likely to control cell-substrate adhesion. The initial area of the explant affects the rate at which it spreads with larger explants spreading faster than smaller ones (Fig 1C-D), suggesting that the initial area of an explant affects the mechanics of the tissue. Since spreading rates may be affected by either the strength of adhesion, lamellipodia formation, or tissue resistance to spreading from either its elastic modulus or the amount of material available, we use parameter estimation and statistical techniques to analyze predictions from simulations of a mathematical model to examine trends in estimated parameter values as they relate to the explants' initial areas. Experimental data such as tissue boundaries and strain are extracted from time-lapse sequences of images to integrate with the mathematical model (Fig 1E; Materials and Methods). (A) Schematic of animal cap ectoderm of a Xenopus laevis embryo during the late blastula/midgastrula stage and at the tailbud stage. At the late blastula stage, a single layer of epithelial cells (blue) lay on top of a deep multilayer of mesenchymal cells (red), which is on top of the blastocoel, a fluid filled cavity in the upper half of the embryo. The red bars indicate where an animal cap explant is cut from the embryo. By the tailbud stage, the ectoderm spreads to cover the entire embryo. The underlying mesenchymal multilayer has undergone intercalation to become a single layer while the cells in the epithelial layer stretch. (B) Schematic of animal cap ectoderm removed from the embryo during the late blastula through mid-gastrula stages and cultured on fibronectin adsorbed glass initially (green), and then after time has passed. As the tissue spreads, deep mesenchymal cells radially intercalate into the epithelial cell layer (purple). Dashed box region shown enlarged below. Note: Experimental images are taken from above, not from the side, so only the pigment-containing apical surface of the epithelial layer is visible. (C) Area of explants over course of the experiment for representative explants in the model building set (blue) and test set (red). Explants shown here are the same explants analyzed in Figure 4. (D) Average change in area over time (ΔA/Δt), calculated as the slope of the regression line of area measurements over time (as in panel C), for Xenopus animal cap explants of different initial sizes. The explants are grouped by initial area into subgroups Ia, Ib, II, and III, indicated by the dashed lines. The blue filled-in circles correspond with the explants included in our model building set and the red circles correspond with the explants included in our test set. The length of time was 10 hours. (E) Time-lapse images of explant spreading are shown in the first row, followed by the boundary identification results of the Level Sets plugin. Then the strains as defined in Eq. 21 are shown and the density ratios (Eq. 16) are shown in the last row. See the Supporting Information Videos for time-lapse sequences of these still frames.
Our experimental data set contained 41 explants from two clutches of eggs. These explants were microsurgically shaped so the initial areas varied eighteen-fold ( Fig 1C). We grouped these explants into four subgroups based on their initial area, where groups Ia and Ib are comprised of the smallest explants, group II has medium sized explants, and group III has large sized explants. We chose eighteen representative explants for our model building set to run the parameter estimation analysis on, leaving the remaining explants available to test the predictive power of the method. See Fig 1C and Supporting Information S1 Table for the naming convention used for the explants and how it corresponds to the initial area of the explant (the smallest explant is named Explant #1).

Mathematical Model
To study Xenopus ectoderm explant spreading, we modify the two-dimensional (2D) cell migration model of Arciero et al. [11], which is derived from mechanics principles and uses a small number of parameters to capture many of the mechanisms that are important in 2D tissue spreading, in particular, adhesion between the tissue and substrate, elastic coupling of cells within the tissue, and forces generated by lamellipodia. This model also reflects the limited complexity of Xenopus explant mechanics and movements that we are able to capture through image analysis.
The model of Arciero et al. [11] was originally applied to small intestinal epithelial cell layers (IEC-6 cells) and Madin-Darby canine kidney (MDCK) epithelial cell layers. Minor modifications to the model are required to apply to explant spreading since the biological mechanisms of material growth in the visible cell layer are different. For example, epithelial cell layers in MDCK sheets are one cell thick and exhibit 2D volumetric growth through cell proliferation. By contrast, the Xenopus animal cap does not undergo 3D volumetric grow through cell proliferation [45] and is a composite tissue consisting of a single epithelial layer underlain by two or more layers of deep mesenchymal cells. Spreading in the animal cap can undergo apparent 2D volumetric growth as cells move from one layer to another. 2D volume (e.g. area) can change as superficial cells actively change shape or as deep cells move into the top layer from deeper within the tissue through radial intercalation ( Fig 1B). (Note: Xenopus embryos do not change total volume from fertilization to feeding tadpole stages.) Since the time-lapse images of explant migration are taken normal to their plane of movement, we assume the explant tissue is a single layer that can be tracked using only the cells in the visible epithelial layer of the tissue. Like the model for MDCK movement we also assume that the explant is uniformly thick, but note that changing this condition would be a useful future extension.
In the model of Arciero et al. [11], a 2D single cell layer is represented by an elastic continuum capable of deformation, motion, and proliferation of cells within the layer. Since the multicellular ectoderm of Xenopus laevis animal caps is continuous and does not break or shear, a continuum model of tissue migration is appropriate [46]. The tissue is represented as a 2D compressible inviscid fluid and the motion of cells is described in spatial (Eulerian) coordinates. The parameters in this model correspond to physical properties in the tissue, and thus we can deduce the material properties of the tissue by estimating model parameter values that best represent the behavior of an experiment. The governing equation of the model of Arciero et al. [11] is where variable ρ(x,t) describes the tissue density as a function of position x=(x,y) and time t, k and b are parameters described in Table 1, q is a growth term, and Ω t is the domain that describes the extent of the tissue at time t (Fig 2A). See the Supporting Information S2 File or Arciero et al. [11] for further details on the derivation of Eq. 1 and the boundary conditions below. Visible material growth rate due to radial intercalation or active cell shape changes in the epithelial layer. Given the presence of lamellipodia on the edge of the ectoderm, i.e. tissue boundary ∂Ω 1 t , we assume there is a constant force per unit length F (see Table 1 and Fig 2A-B) exerted outward at the tissue boundary that is equal in magnitude to that of the force of the cells in the interior. To express this boundary condition in mathematical notation, a function describing the forces within the tissue is necessary. The constitutive relation chosen for the pressure within the tissue by Arciero et al. [11] was where ρ unstressed is a parameter described in Table 1. Thus with Eq. 2 and setting p = −F at the boundary [11,47], we obtain the boundary condition The constitutive relation in Eq. 2 is also the cause for the appearance of the Laplacian in Eq. 1, hence the Laplacian does not arise from any underlying diffusion process or Brownian motion and k/b should not be interpreted as a diffusion constant (see Supporting Information S2 File or Arciero et al. [11]).

An additional Stefan moving boundary condition is given by
where v(x,t) is the velocity of the layer and n 1 (x,t) is the outward unit normal to the tissue boundary (see Supporting Information S2 File or Arciero et al. [11]).
On the edge of the computational domain ∂Ω 2 , we assume that there is no flux of cells, i.e. cells are unable to move beyond this boundary, and so we have the Neumann boundary condition where n 2 (x,t) is the outward unit normal to the edge of the computational domain. We note that in our simulations, the computational domain ∂Ω 2 is large enough so that the tissue boundary ∂Ω 1 t is never close to ∂Ω 2 .
By segmenting cells in confocal images of the epithelial layer of a representative animal cap explant labeled with a GFP-membrane tag, we measured the tissue density in the epithelial layer at the initial imaging time point to be 0.0047 cells/µm 2 , so we take the initial condition to be In the original model of Arciero et al. [11], the growth term q was included in Eq. 1 to model the volumetric proliferative growth which is typical for spreading sheets of cultured cells, such as IEC-6 or MDCK cells, and taken to be logistic growth. However, the cleavages in early Xenopus embryonic tissues are reductive, so migration of explants does not depend on volumetric growth provided by cell proliferation [45]. Therefore, we modify q in Eq. 1 to represent material added through other means, such as active cell shape change or radial intercalation in the epithelial layer.
We observe that since there is a finite number of cells in the initial explant and since tissue mass is conserved, there cannot be unconstrained material growth in the visible epithelial layer.
Furthermore, we estimate that there is at most a 50% increase in number of cells on the surface of the initial explant from either active shape changes or intercalated cells from the lower mesenchymal layer, which is a reasonable upper bound based on the estimate that no more than 20% of the surface cells in the embryo come from the deep mesenchymal layer [38]. Thus, after testing various possible growth functions for qualitative fits to explants' change in area over time, we use a mass-limited densitydependent logistic growth function to describe this addition of visible material into the surface layer, where α is a parameter described in Table 1, ! ρ is the limiting density assumed to correspond to a compressed tissue such that ! ρ = ρ unstressed e F/k [11], m(t) is the cumulative amount of mass that has been added to the visible surface layer through time t, and m 0 is the initial amount of mass.
Details of the numerical method implemented to solve the model equations may be found in Arciero et al. [11], and the output of a toy simulation is shown in Fig 2C.

Parameter Estimation Using Approximate Bayesian Computation (ABC)
Equipped with a mathematical model of collective cell migration in Xenopus ectoderm during gastrulation to tailbud stages and kinematic data extracted from time-lapse images, we aim to infer material properties of the tissue that are not directly measurable as well as predict tissue behavior in new situations via estimating model parameters. We use a Bayesian framework, the approximate Bayesian computation (ABC) rejection method [33,34], to determine probability distributions for the model parameters for an individual tissue explant; see Fig 2E for a flow chart depicting this process.
The results for the individual explants from our model building set of explants are then compared to determine whether there is a difference in parameter values for explants of various sizes (Fig 1C). A set of explants is reserved for our test set to examine the predictive power of the mathematical model.
Estimating parameters for our model in a deterministic manner appeared infeasible as preliminary parameter space exploration hinted that our model had parameter structural nonidentifiability, or functionally related model parameters [48]. While parameter identifiability has been explored extensively in the case of ordinary differential equations, which involves studying issues such as whether different parameter sets give rise to different model predictions, it has not been studied as much in partial differential equation models [49]. Thus, due to our preliminary parameter space exploration analysis and aware that other partial differential equation models that are similar to our collective cell migration model have shown non-identifiability relations between parameter values [50], we use ABC rejection because it allows for quantification of the uncertainty of our parameter estimates. Furthermore, since we aim to compare multiple data sets, being able to quantify uncertainty in a single data set permits us to analyze the then more complex problem of comparing estimated parameter values among multiple data sets.
To use ABC rejection, the first step is to specify a prior distribution for each of the parameters.
This prior distribution is our best educated guess as to what the landscape of the parameter value looks like. The parameter values to be estimated are F, k, b, α, and ρ unstressed (Table 1) where index n denotes the computational points counted along the edge and N is the total number of these points. Summing over all time points gives the distances error term For the density ratios, at each pixel we approximate the experimental density ratio, denoted as ξ exp,j , from Eq. 16 (Materials and Methods Section) with the deformation gradient Eq. 18 for the j th (corresponding with ρ r ) and (j+δ) th (corresponding with ρ s ) images in the time-lapse sequence, where δ denotes the frame increment between the pair of images. To limit the impact of inelastic mechanical events such as cell rearrangement or cell division, e.g. cases where tissue position becomes discontinuous, we limit the time between the pairs of registered images. δ is chosen to be the smallest frame increment in which changes between images can be detected, and in our case, δ=5 (Supporting Information S1 File). Since the size of the matrix ξ exp,j may not be the same size as the computational grid (chosen for computational efficiency), we first linearly interpolate the experimental density ratio values at the computational grid nodes using the MATLAB (The MathWorks, Natick, MA) function interp2. The computational density ratio is denoted ξ comp,j , which is computed at each grid node by taking the ratio of the density from the j th time point (corresponding with ρ r ) and the (j+δ) th time point (corresponding with ρ s ). In other words, to compute ξ comp,j where ρ(x,y,t) is the computational density for each element ℓ 1 , ℓ 2 ( ) in the grid, The square root of the average of the squares of the differences between the experimental and computational density ratios is where ξ comp,j and ξ exp,j are evaluated at grid node ℓ 1 , ℓ 2 ( ) , and the size of the computational grid is L 1 ×L 2 . Summing over all time points (minus frame increment δ) gives the density ratios error term We take the total error to be where we use weight w = 1000 to ensure that the distances error term z d and density ratios error term z p are approximately the same order. If the simulated explant migrates outside the computational domain, the error is taken to be "NaN" in our numerical code to signify an extremely poor fit since the entire explant remains within the domain during the time-lapse sequence. Furthermore, if the density in the center of the explant increases in the numerical simulation, the error is also taken to be "NaN" since explants thin out in our experiments.
We collected 10,000 parameter sets that result in total error (Eq. 13) less than or equal to a  Model parameters describing explant spreading reveals an inverse relationship between F/k and k/b, which implies non-identifiability, and less distinct relationships between the other pairs of parameters, which may imply identifiability or weak non-identifiability ( Fig 3A). This relationship is further indicated in Supporting Information S1 Figure, which colors the results of the triangle plot in Fig 3A by the value of the error (Eq. 13) instead of frequency. To verify the association between the different parameters, we calculated the correlation coefficient for each pair of parameters for each explant (Fig 3B). We found that F/k and k/b were strongly negatively correlated, pair k/b and α and pair α and ρ unstressed were weakly correlated, and the other parameter pairs were moderately correlated.
We also calculated 95% confidence intervals for the one-dimensional projections of the posterior distributions for each parameter for each explant; see Supporting Information S1 Table. We observe that there appears to be trends in the mean parameter value versus the initial area of the explant, but to verify this we ran ANOVA tests. Using 10,000 parameter sets per explant and grouping them into regions Ia, Ib, II, and III (in other words, the parameter sets for Explants #1-5 were combined into one set, the parameter sets for Explants #6-12 were combined into one set, etc.), we determined that, for all the parameters F/k, k/b, α, and ρ unstressed , the p-value < 10 -29 . This implies that there is statistically significant evidence that not all of the average parameter values were the same for the four groups of explants. We repeated this analysis without grouping the explants into regions but instead individually comparing each, and found that the p-values were even smaller.
Additionally, we ran Tukey's multiple comparisons test to compare the average parameter values of the eighteen explants, grouped into the four regions and individually, and to determine which showed statistically significant differences. Figure 3C shows which is a realistic expectation since explants are isolated from similar locations, then these results imply that for larger explants the forces produced by the lamellipodia (F) are stronger and the adhesion forces (b) are weaker. A similar relationship between the lamellipodia force and the size of the initial explant was also observed in primary mouse keratinocyte cell colonies [51,52], where traction stresses increase with colony radius. In addition, from the inverse relationship between F/k and k/b in the triangle plot (Fig 3A), we predict that if adhesion is increased, then the measured traction force will also increase.

Prediction
To

Discussion
Time-lapse imaging of tissue explant migration generates complicated datasets that are information-rich. Here, we present a formal methodology for extracting mechanical mechanisms of tissue spreading by integrating complex image information with a mathematical model of cell migration. This approach identifies key differences between experimental groups of differently sized embryonic tissue explants. In order to directly relate the experimental kinematics of spreading of multiple tissue explants to mechanical properties of the tissue, we focused on one mathematical model, the previously-developed Eulerian mechanical model of Arciero et al. [11], that represents forces involved in migration within a 2D elastic continuum. In this modeling framework, extended here to represent Xenopus animal cap epiboly, we can directly couple experimentally derived kinematic data to model parameters that correspond to physical properties in the tissue.
Since cell nuclei positions were not recorded in our time-lapse images, we mapped strain across the tissue and used changes in strain to estimate local changes in tissue density. While in the current study we used variations in pigment to calculate strain and time-dependent changes in density, any landmark-rich image, for example, random textures of sub-cellular organelles, dye, or fluorescence, could be used to obtain density ratios. Furthermore, our assumption that the tissues are elastic for the image registration and the mathematical model, could be relaxed to accommodate more realistic constitutive models such as viscoelasticity or viscoplasticity [53][54][55]. Since strain mapping techniques do not depend on free edges, the strain mapping techniques presented here can also be extended beyond tissues in two dimensions to analyze more complex structures, organs, and whole organisms.
Previous studies of tissue spreading that used methods developed by Arciero et al. [11] and Mi et al. [14] analyzed two different cultured cell lines, rat intestinal epithelial cells (IEC-6) and Madin-Darby canine kidney (MDCK) epithelial cells, and estimated model parameters for individual experiments. While these studies report one parameter set that matches a single set of experimental data, we show here that while there may be a single "optimal" parameter set for each set of experimental data, due to the inverse relation between model parameter ratios F/k and k/b, there is a wide range of "close to optimal" parameter sets. Thus, to compare between different sets of experimental data, we needed to implement a robust statistical analysis capable of accounting for the uncertainty in the model parameter estimates.
The approximate Bayesian computation (ABC) rejection method was utilized to collect multiple parameter sets that, when used in a computer simulation of the mathematical model, resulted in sufficiently good fits to the experimental data. Fit was determined via an error function that calculated the difference between the model generated boundary and strain field and image-based measurements of tissue boundary and density changes throughout the simulation, rather than just matching a single end time point. Hence, we used as much as possible of the available image dataset to estimate physically-relevant parameter values. Statistical analysis then allowed us to examine the difference in means in multiple parameter set distributions for the multiple experiments of explants of varying sizes. Using this methodology, we found that there are interdependencies between parameters F/k and k/b. We also detected trends that for larger explants, adhesion forces (b) are weaker and forces produced by lamellipodia (F) are stronger if we assume that the stretching modulus k is fixed among explants since they are comprised of the same tissue. We used the same prior uniform distribution for each of the explants (Fig 2D) and observed that in general, more simulations needed to be run to collect 10,000 parameter sets that resulted in an error (Eq. 13) under the required threshold for smaller explants compared to larger explants (Supporting Information S5 Figure). The parameter values that fit the experimental data the best for smaller explants appear to be more on the extremes of the prior uniform distribution than for larger explants.
Large computational expenses limited the number of parameter sets gathered for each explant to 10,000 for the ABC method as well as the size of our model building set to 18 of 41 available explants in our experimental data set. However, we observed that the parameter values found were robust among similarly sized explants (Fig 4). Indeed, for explants with similar average change in area over time (ΔA/Δt), simulations with a parameter set that resulted in a small error (Eq. 13) for an explant in the model building set also resulted in a small error for a similarly sized explant in the test set. Thus we would expect that the trends we discovered would be the same if we analyzed more explants in the future, taking into account embryo to embryo variation in stiffness [56], clutch to clutch variation [57], and possibly environmental variations such as room temperature that might affect rates and deformation maps of tissue spreading [58]. We also propose adding the velocity or the change in area over time as another input to the error function Eq. 13 as it appears to have an effect on estimating parameters (Fig 1C and 4). Gathering more parameter sets through the ABC method will also allow for the probability distributions of the parameters to become more precise. Further investigation of the estimated set of model parameters would require the ability to measure one of the parameters directly from experiments, reducing the total number of parameters.
The probability distributions provide ranges for parameters F/k and k/b that can be compared to direct measurements or may be manipulated experimentally (Fig 4).

Embryos and sample preparation
Eggs were collected from Xenopus laevis frogs and fertilized in vitro using standard methods [45]. After fertilization eggs were dejellied and cultured in 1/3 X Modified Barth's Saline (MBS) [60].
Embryos used to analyze intercalation were transferred to a 3% ficoll solution in 1X MBS and were injected (Harvard Apparatus) with approximately 1ng mRNA encoding a membrane tagged GFP (mem-GFP [61]) at the two cell stage, after which they were returned to 1/3 X MBS. Embryos were raised to approximately Stage 10 [37] and were transferred to explant culture medium (Danalchik's for Amy; DFA) [62]. Embryos were devitelinized with forceps and the animal cap ectoderm was microsurgically removed using hair tools (Fig 1A-B). Tissue explants were cultured on either a petri dish (Fisher) or in a glass-bottomed chamber coated with human plasma fibronectin (Roche Molecular Biochemicals) in DFA supplemented with antibiotic/antimycotic (Sigma The embryonic tissue explants generally spread symmetrically outward over their fibronectincoated substrate in our experiments. Any variation in spreading was likely due to unseen substrate defects or anisotropy in fibronectin substrate density but did not affect the overall spreading rate of explants. As expected, reducing the fibronectin concentration below 1.0 µg/mL reduced the spreading rate (Supporting Information S4 Figure).

Imaging
Stereoscope images of tissue explants spreading were collected using a CCD camera (Scion Corporation) mounted on a stereoscope (Stemi 2000; Zeiss) ( Fig 1D; first row). Multiple explants could be imaged in the same dish using a motorized XY stage (Marzhauser and Ludl) controlled through microscope automation and image acquisition software (µ-Manager [63]). In order to determine cell density, high resolution confocal images of membrane tagged GFP in live samples were collected using an inverted compound microscope equipped with a laser-scanning confocal scanhead (SP5; Leica Microsystems). In some cases, a stack of images over a z-range of 50 to 100 µm were collected and combined into a single image using maximum intensity projection prior to analysis.

Image Analysis
Subsequent image analysis of time-lapse sequences were performed using ImageJ ( Tissue edge boundaries were identified from single stereoscope images using the Level Sets plugin in Fiji [65], which is an image segmentation technique based on partial differential equations ( Fig 1D; second row). The error is expected to be less than 5 pixels, which corresponded to no more than 2 cell widths (approximately 40 µm). To calculate the average areal spreading rate, the area within this boundary was determined at each time point. These areas were then fit using linear regression. The slope of the resulting regression line was defined as the average change in area over time (ΔA/Δt) ( Fig   1C).

Strain Mapping Method
Time-lapse sequences provide a rich set of information on the movement of tissues. In addition to edge detection, digital image correlation can be used to track individual features such as cell nuclei, cell membranes, or even textural features such as pigment granule distribution that do not change significantly between sequential images [42,67]. By tracking these features, a "translation" map or Eulerian description of flow can be produced showing where parcels of tissue move over time [43].
This map can then give an estimate of local density change.
Elastic registration between two images in a time-lapse sequence separated by a fixed amount of time represents the kinematic transformation between the two images and makes minimal assumptions about the mechanical properties of the materials involved. The usefulness of the registration, i.e. that it reflects linear deformations in the ectoderm, relies on additional assumptions based on the biology of Xenopus embryos. Since cell division, rearrangement, and radial intercalation events are typically slow, these processes do not disrupt the pigment patterns observed at the apical face of the ectoderm.
Since we are examining a short-term analysis of strain in time, we assume that the Xenopus tissue is in a quasi-static equilibrium and mass is added very slowly between each pair of images, and thus we may assume that there is local conservation of mass. Letting dV r be the infinitesimal volume in the reference material coordinates and dV s be the corresponding infinitesimal volume in the spatial coordinates, the change of volume due to deformation [68] can be written as where F is the deformation gradient defined as F mn (X, t) = ∂x m (X, t) ∂X n , where X=(X 1 ,X 2 ) is the (x,y)coordinate positions in the first still image in pixels and x=(x 1 ,x 2 ) is the (x,y)-coordinate positions in the second still image.
Eq. 14 implies that the Jacobian det F measures the ratio of the volumes. The conservation of mass equation ρ s dV s = ρ r dV r , thus implies that, locally, Hence, if we can calculate the deformation gradient F, we can obtain an estimate of the local change in density.
To perform digital image correlation on image pairs of pigmented ectoderm tissue and calculate strain, or local deformation, we use a texture mapping strategy with the ImageJ software plugin bUnwarpJ, which is used for elastic and consistent image registration ( [69]; available for download at http://fiji.sc/BUnwarpJ). By using differential image correlation techniques [70], we calculate deformation of the tissue via the translation map (Fig 1D; third-sixth rows). For details on the strain mapping method discretization to obtain numerical estimates of the x-strain ε xx , y-strain ε yy , xy-strain ε xy , yx-strain ε yx (which is equal to ε xy ), and displacement gradient ∇u, please see the Supporting Information S1 File.
Since the gradient of the displacement vector is defined as where I is the identity matrix, using the discretization of the displacement gradient ∇u(i,j) we numerically approximate the deformation gradient at each pixel (i,j) as Taking the determinant for each i and j, Eq. 18 gives a numerical approximation of the density ratio in Eq. 16 at each pixel ( Fig 1D; last row).

Supporting Information
Supporting Files S1 File. Strain mapping method discretization.
Here we describe the implementation of the strain mapping method to calculate the deformation of a tissue via estimates of the x-strain ε xx , y-strain ε yy , xy-strain ε xy , yx-strain ε yx , and displacement gradient ∇u between two images in a time-lapse sequence.
For images of width n pixels and height m pixels, the properties of each pixel can be represented by an entry in an m×n matrix. For each pair of images from the time-lapse sequence, let X=(X 1 ,X 2 ) be the (x,y)-coordinate positions in the first still image in pixels and x=(x 1 ,x 2 ) be the (x,y)coordinate positions in the second still image. The top left corner of an image is the origin, and x increases from left to right and y increases from top to bottom. The entries of X are therefore defined by X 1 (i,j)=j−1 and X 2 (i,j)=i−1 for i=1,2,…,n, and j=1,2,…,m.
We mask the pair of images so that registration occurs only in the actual location of cells. We then initialize bUnwarpJ to calculate the coefficients of the cubic B-spline map β that defines the transformation (X 1 , X 2 ) β ! → ! (x 1 , x 2 ) [69]. Initializing bUnwarpJ again, we apply β to X by converting the transformation to "raw" data, which reports the elastically-mapped position of each pixel in the first image to the "fitted" position in the second image. We obtain x, the mapped position of each pixel from the first image to its position in the second image, in pixels. Note that pixels outside of the mask will be mapped as well, but we will remove this extraneous data after we have obtained the strains. Next, we calculate the displacement vector u by u 1 (i, j) = x 1 (i, j) − X 1 (i, j), i = 1, 2,..., n, j = 1, 2,..., m, u 2 (i, j) = x 2 (i, j) − X 2 (i, j), i = 1, 2,..., n, j = 1, 2,..., m.
The engineering, or Cauchy, strain is defined as where ΔL is the change in length of the tissue, L 0 is the original length, and L is the current length. The displacement vector u is converted into x-strain, y-strain, xy-strain, and yx-strain by