Modeling the roles of cohesotaxis, cell-intercalation, and tissue geometry in collective cell migration of Xenopus mesendoderm

The collective migration of cells is a feature of many biological processes including wound healing, tissue regeneration, cancer invasion, and embryo and tissue morphogenesis (Friedl and Gilmour, 2009; Li et al., 2013; Mayor and Etienne-Manneville, 2016; Scarpa and Mayor, 2016; Spatarelu et al., 2019; Yang et al., 2019). Collectively migrating groups of cells are distinguished from motile single cells by the coordinated polarization and directed migration of distinct leader- and follower-cell populations connected via cadherin-based cell–cell adhesions (Gov, 2007; Omelchenko et al., 2003; Qin et al., 2021; Rørth, 2012; Shellard and Mayor, 2019). Leader cells exhibit polarized protrusive behavior in the direction of the collectively migrating front. In many cases, follower cells are pulled along as a result of substrate traction forces generated by leader cells. Leader–follower behavior is developed and reinforced largely through differential Rac- and Rho-mediated signaling cascades in these two populations (De Pascalis et al., 2018; Khalil and de Rooij, 2019; Sonavane et al., 2017; Weber et al., 2012; Zhou et al., 2022). A more complete understanding of the cellular and subcellular mechanisms underlying the regulation of collective cell migration in vitro remains a focus of active investigation. In addition, there is a need to identify the critical factors governing collective cell movements in vivo, where more complex tissue geometries (e.g. circularly arranged tissue) and increased length scales play roles (Davidson et al., 2002; Keller and Shook, 2008; Keller et al., 2003; Shook and Keller, 2008).

During gastrulation in the amphibian Xenopus laevis, a coextensive population of mesoderm and endoderm, collectively termed ‘mesendoderm’, involutes and spreads across the underside of the overlying blastocoel roof (BCR) as a circular mantle of tissue. The advancing 360° leading edge of the mesendoderm mantle ultimately meets and fuses to form a layer of contiguous tissue beneath the BCR by the end of gastrulation (Keller, 1976, 1991; Keller and Jansa, 1992; Keller and Tibbetts, 1989). The assembly and arrangement of these tissues are regulated through multiple mechanisms, including the initiation of mesendoderm involution, the coordination of single cell protrusive activity and collective cell movement, and the directionality and velocity of migration. Leading edge mesendoderm migrates across the BCR using integrin-dependent adhesions to fibronectin (FN) fibrils, which are assembled along the inner surface of the BCR (Davidson et al., 2002, 2004; Rozario et al., 2009; Smith et al., 1990; Sonavane et al., 2017; Winklbauer, 1990). Several mechanisms have been reported to influence this migration. These include a chemotactic gradient of platelet derived growth factor (PDGF) in the BCR (Ataliotis et al., 1995; Nagel et al., 2004; Scarpa and Mayor, 2016; Symes and Mercola, 1996) and a force-dependent process termed cohesotaxis (Weber et al., 2012). Both of these mechanisms result in polarized protrusions in the forward direction of travel. Cohesotaxis involves the tugging of cells on one another as a result of integrin-dependent traction forces generated at the front and balanced by rearward cell–cell adhesions and the recruitment of a mechanosensitive cadherin-keratin complex to sites of cell–cell contact (Weber et al., 2012). Under these conditions, Rac activity and protrusions are inhibited at the rear of leading row cells but increased at the cell front (De Pascalis et al., 2018; Helfand et al., 2011; Sonavane et al., 2017). A similar mechanism, termed ‘plithotaxis’, involves the orientation of collective cell migratory behavior in culture along the axis of maximal principal stress (Roca-Cusachs et al., 2013; Tambe et al., 2011). In cohesotaxis, the mechanical coordination of cells within the collective results in high traction stresses concentrated primarily at the leading edge (Sonavane et al., 2017). These forces are balanced by inter- and intra-cellular stresses in following rows. Additionally, Xenopus mesendoderm cells often rearrange relative to one another during closure of the tissue mantle. These cell intercalation movements can occur in both radial and medio-lateral directions (Davidson et al., 2002) and may affect the speed of mantle closure. Ultimately, factors that drive tissue deformations required for vertebrate morphogenesis such as mesendodermal mantle closure include: (1) individual cell protrusive and migratory behaviors that are self-driven; (2) cell protrusive and migratory behaviors that are influenced by connections to neighboring cells (e.g. cohesotaxis and intercalation); and (3) tissue-level geometries (e.g. a circular multi-cell layered mantle of tissue) that provide physical boundary constraints.

A variety of tissue explant preparations has been used to investigate the collective cell movements of Xenopus mesendoderm (Fig. 1A). These include the dorsal marginal zone (DMZ) explant, the toroid or ‘donut’ explant composed of the bulk of the mesendodermal mantle, and the ‘in the round’ (ITR) explant prepared with four juxtaposed DMZ explants to approximate the toroidal geometry of the mesendoderm (Davidson et al., 2002). Despite the utility of these explant preparations for studies of collective cell migratory behaviors, how cell and tissue geometry contribute to mesendoderm mantle closure remains difficult to address via experiment alone. Dissection and removal of this tissue from the embryo disturbs features of its native geometry and association with neighboring embryonic tissues (e.g. the BCR and ventral mesoderm and endoderm). Moreover, individual cell rearrangements and intercalation behaviors at tissue-length scales (e.g. in the embryo and/or some explants) are either difficult or currently not possible to image live, or perturb in isolation in the in vitro setting.

The three-dimensional tissue geometry and migration behaviors that characterize mesendoderm mantle closure on the BCR arise in large part from the individual actions of leader and follower cells, and can be considered emergent properties of the tissue. Computational methods such as agent-based modeling are often applied to investigate emergent behaviors that may otherwise be resistant to experimental observation and manipulation (An et al., 2009; Chan et al., 2010; Glen et al., 2019; Yu and Bagheri, 2020). The Cellular-Potts model, also known as the Glazier-Graner-Hogeweg model, is a mathematical framework that allows for agent-based modeling while effectively capturing the inherent stochasticity of cell and tissue properties, including differential adhesion between cells and tissues, and forces exerted by cells upon one another and their extracellular environments. Additionally, the Cellular-Potts model framework has been successfully applied to investigate multiple collective cell movements of interest, such as convergent-extension in Xenopus laevis, making it a promising choice for investigating tissue emergent behaviors that drive morphogenesis (Belmonte et al., 2016; Graner and Glazier, 1992; Swat et al., 2012; Zajac et al., 2003).

We present a novel computational model in the Cellular-Potts framework of a DMZ explant of collectively migrating mesendoderm. This model can be used as a methodological examination of collective cell migration, and we apply it here to investigate mesendoderm movements in the Xenopus laevis gastrula. The model replicates the behaviors of a single DMZ biological explant with leader and follower cell agents. We use the model to simulate more complex geometries representative of those in the embryo (Fig. 1B), and we conduct sensitivity analyses and run in silico experiments to investigate how collective cell migration dynamics emerge as a function of cell intercalation, tissue geometry, and cohesotaxis. Our model suggests how cohesotaxis, cell intercalation, and circular tissue geometry each contribute to affect cell migration speed during morphogenesis.

Data from Davidson et al., 2002, compared the leading-edge migration speeds of ‘linear’ DMZ explants with toroid-shaped ‘donut’ explants. The migration speeds of the toroid explants, which are representative of the tissue geometry in situ were over double that of the linear DMZ explants (107 μm/h versus 237 μm/h, Fig. 1A). However, arranging four linear DMZ explants in a juxtaposed configuration to create the ITR configuration after collision, resembles the circular toroid explant (Fig. 1A), and exhibits migration speed similar to that of the toroid (236 μm/h).

We hypothesize that the circular geometry of the explant tissue in the ITR and toroid configurations directs migratory behavior towards the cell-free opening in the center. Further, the circular tissue geometry enforces a physical boundary on collective cell migration that enables the speed of mesendoderm mantle closure to surpass the rate of collective cell migration when the cells are unbounded on either side, as they are in the single DMZ explant. To address this, we constructed a computational model of the single DMZ explant, which is displayed in Fig. 1B. Cells were represented as leader and follower cell agents, two populations of cells with distinct behaviors as described for many examples of collective cell migration. We constructed the in silico representation of the DMZ explant using leader and follower cell agents with behaviors described in Materials and Methods. Timelapse images of live DMZ explants (representative example: Movie 1) were used to guide computational model construction. Representative images of the computational model with a schematic indicating critical model parameters are displayed in Fig. 1B, and video of the DMZ model migration is included in Movie 2.

As described in Materials and Methods, six key parameters were built into the model to control migratory behavior in the computational explant (see schematic in Fig. 1B). These were queried to determine which parameters had the greatest effect upon migration speed of the explant (i.e. those parameters most sensitive for migration speed). We performed a one-dimensional sensitivity analysis to measure the effect of each of these parameters on forward migration speed. The analysis was performed by sequentially increasing and decreasing each parameter, while holding all other parameters constant. Parameters were varied in increments of ±10% until ±50% variation from calibrated baseline levels was reached, with 10 simulated replicates performed for each parameter variation, and the percent change in migration speed from reference was observed. The parameters varied included: stiffness of cell–cell connections between cells in the tissue (λ_Tissue), cell-substrate connections for leader cell agents (λ_Lamellipodia) and follower cell agents (λ_{FollowerSubstrate}), the probability of forming new connections for cells to their substrate (ζ_Lamellipodia for leader cell agents, and ζ_Substrate for follower cell agents), and the distance away from a leader cell that a link is formed to reproduce a lamellipodial extension (χ_Lamellipodia), as shown in Fig. 2 and diagrammed in Fig. 1B. An exhaustive list of parameters and their descriptions are included in Materials and Methods.

Our sensitivity analysis revealed that migration speed was significantly affected by the parameters χ_Lamellipodia, λ_Tissue, and λ_Lamellipodia according to one-way ANOVA. We fitted a line of best fit through the means of these data, the slope of which corresponds to sensitivity of the output to change in the parameter. Migration speed in the DMZ Explant model was most sensitive to the two parameters χ_Lamellipodia, and λ_Lamellipodia. These two parameters are specific to the leading-edge agents (i.e. cells in the leading row), suggesting that the actions of the leading edge cells have a greater effect on the migration speed of the DMZ explant than behaviors of the follower cells. Increasing and decreasing λ_Lamellipodia was expected to significantly affect migration speed because it represents the magnitude of force exerted by a leading edge cell's lamellipodial protrusion. However, the explanation for why altering the parameter χ_Lamellipodia, affects migration speed is not as obvious. While this parameter does not affect the force exerted upon leader cell agents, when this parameter is increased the Cellular-Potts link object is created farther away from a cell agent. In the simulation, leader cell agents migrate forward until they reach their lamellipodia link target, at which point they stop migrating. Link formation is a stochastic process, such that a longer link allows for more time for a new link to form, and thus fewer interruptions in migration. This may allow the DMZ explant to achieve a greater migration speed due to decreased interruptions in migration on a single-cell basis.

To investigate the effect of a migratory bias on the migration speed of the computational explant, we encoded a mechanism by which leading edge cells are more likely to form lamellipodial protrusions in the direction of migration rather than randomly towards areas of free space. A directional and persistent migratory bias of leader cells along their substrate can be explained by cohesotaxis, a mechanism by which migratory cells in a tissue coordinate mechanically resulting in directed polarized protrusive activity at the leading edge in the direction of migration where cells remain in contact with their neighbors (Weber et al., 2012). To represent cohesotaxis in the model of the DMZ explant, we created a rule for the leader cells such that they preferentially extend lamellipodia away from sites of end-to-end cell–cell adhesion resulting in lamellipodial extension in the direction of migration rather than towards any free edge, thus representing polarized protrusive activity in leader cell agents (Fig. 3A). Leading edge cells select a voxel through which a new link's direction is established preferentially from those with the lowest cumulative distance to cell–cell borders in the explant, resulting in creation of lamellipodia away from cell–cell connections in the tissue and towards the direction of migration to maintain persistent migration. The voxel is selected with a probability as displayed in Fig. 3B, where parameter κ approximates a kurtosis to modulate the extent to which migratory bias is applied during the simulation. When κ is smaller, there is less migratory bias in the direction of migration, and when κ is larger, there is a greater migratory bias. For further explanation of this parameter, please see Materials and Methods.

In silico experiments were performed to investigate the effects of the migratory bias on the migration speed of the explant. The cohesotaxis parameter κ was tuned from a value of −6, representing zero migratory bias, to +6, representing maximal migratory bias. The model was run for n=10 simulations for each parameter value κ (Fig. 3C). Imposing a stronger directional migratory bias by increasing the cohesotaxis parameter resulted in a significant increase in migration speed (P<0.05 by independent sample t-test).

We next investigated whether the ITR geometry would increase migration speed. The ITR configuration was simulated by placing four DMZ explants opposed to one another in the simulation environment, just as was performed with biological explants (Movie 3). Migration speed in the ITR configuration was significantly greater than the DMZ explant with cohesotaxis parameter, κ, equal to 6 (P<0.05, independent samples t-test, Fig. 4B).

To determine whether the ITR geometry provided a directional migratory bias to affect migration speed, we varied the cohesotaxis parameter, κ, throughout the range of 6 to −6 (Fig. 3). There were no significant differences between the measured migration speeds across the range of κ values (n=10, P>0.05 by one-way ANOVA, Fig. 4C). The ITR geometry results in the same migration speed regardless of the cohesotaxis parameter. Thus, we suggest that a similar directional bias on cell migration emerges from the ITR geometry. The ITR geometry restricts the direction of migration to the free area bounded by leading edge mesendoderm biasing the direction of migration entirely towards the center of the free area. This contrasts with a single DMZ explant, where the free area in front of the leading edge remains unbounded resulting in a ‘fanning out’ of the leading edge, which slows migration as can be observed in biological explant behavior (i.e. both convex morphology and speed of migration). The migrational bias that emerges from the tissue geometry may be a factor assisting mesendoderm mantle closure in vivo.

ITR experiments using live explants were repeated to confirm the findings reported in Davidson et al., 2002, and to further investigate whether the migratory behavior of DMZ explants was affected by tissue geometry. The migration speed of a single experimental DMZ explant (N=1, n=4) was compared with the measured speed of five sets of four DMZ explants each placed ITR (N=2, n=10) (Fig. 4E). A movie of a representative ITR experiment is included in Movie 4. A zoomed-in view of one corner highlighting the interface between two colliding explants in the ITR configuration both in vitro and in silico is included in Movie 5. Before the adjacent in vitro DMZ explants made contact with one other at their respective edges, there were no significant differences in migration speeds of the DMZ explants as compared with a single explant in vitro. However, after the adjacent DMZ explants made contact with one another at their respective edges, there was a significant increase in migration speed (Fig. 4E).

Davidson and colleagues (Davidson et al., 2002) previously reported intercalation behaviors in the DMZ explant, with cells leaving and entering the substrate level, as well as leaving the leading edge. These behaviors were absent in the model as described through Fig. 4. To investigate the role of these cellular rearrangement behaviors, we added an additional rule to our model where cell-to-cell link objects that orient cells in fixed locations relative to one another were allowed to break probabilistically based on the parameter ζ_Tissue, as described in Materials and Methods. This behavior allowed for passive cellular rearrangements to occur within the in silico tissue, resulting in a different appearance during ITR tissue closure in Fig. 5A as compared with Fig. 4A due to cell intercalation events within the tissue. A representative simulation time course is shown in Movie 6. While the simulations without intercalation (Fig. 4A) retained their rigid bi-layer structure, the simulations with intercalation behaviors (Fig. 5A) were observed to flatten to a monolayer over the time course of tissue closure.

Along with the observation that cellular rearrangement via radial intercalation occurs during the ITR closure, previous studies suggested that intercalation may contribute to increased migration speed in the ITR and toroid mesendoderm geometries (Davidson et al., 2002). To test this hypothesis, we performed an in silico experiment to compare the ITR model with a ITR model that includes intercalation behavior (parameter ζ_Tissue=0.1). A significant increase in explant migration speed (n=10, P<0.05, via independent sample t-test) was observed with the inclusion of this parameter (Fig. 5B). The probability of the behavior that allowed for cell rearrangement to occur was then progressively reduced, from the ζ_Tissue value of 0.1 where cellular rearrangement is observed, to ζ_Tissue=0.0001, where cellular rearrangement events become rare and migration speed of closure was approximately equivalent to the model without intercalation. We observed a monotonic decrease in the migration speed with decreasing probability of intercalation behaviors in the tissue (Fig. 5C). To validate this prediction in vitro, we sought to increase cell intercalation by reducing Ca⁺⁺ present in the DFA media in order to ‘loosen-up’ cadherin-based cell–cell adhesions, allowing for greater cell–cell rearrangement within the tissue. We measured mesendoderm explant migration rates on FN in standard DFA media (1.0 mM Ca++) versus DFA media containing 0.8 mM Ca++ (i.e. a 20% reduction in Ca++), to weaken cadherin-based cell–cell adhesion. Explants migrated significantly faster in 0.8 mM Ca++ DFA than in standard DFA, over the first hour of migration (Table S1), consistent with in silico model predictions from Fig. 5B,C. The migration rate in reduced Ca++ DFA is highest over the first hour (Table S1), which may reflect the limited number of mesendodermal cell layers (two to three layers) that are available to intercalate onto the fibronectin substrate, such that the frequency of intercalation events decreases over time.

Cellular rearrangement resulted in a flattening of the in silico explant, as cells from the layer above shifted down onto the substrate level. We quantified the number of cells on the substrate level for different values of parameter ζ_Tissue and observed an increase in the number of cells at the substrate level over the time course of the ITR simulation, as well as in a single DMZ explant simulation, with larger values of ζ_Tissue (Fig. 5D). This confirms that the top-down intercalation behavior occurs and with greater frequency in both the ITR model and in the single DMZ model with increasing ζ_Tissue.

We then confirmed by experiment the model prediction that cell intercalation towards the substrate occurs in the DMZ explant. Fluorescent Dextran-labeled mesendoderm cells were added to the top of an unlabeled DMZ explant migrating on a fibronectin substrate. Over the 30-min time course of the experiment, labeled cells intercalated in a downward or radial direction toward the substrate to become visible in the bottom layer of cells (Fig. 5E, Movie 7). In control DMZ explants lacking any additional labeled cells, cells initially residing in layers above the bottom-most layer of the DMZ explant gradually appeared in the bottom layer over the time course of the experiment, providing further evidence of intercalation towards the substrate (Movie 8).

In this study, we present an agent-based computational model of migrating layers of mesendoderm cells. Our model simulates multiple aspects of collectively migrating cells, including leader and follower cell dynamics, mechanical coordination of cells within the collective in the form of cohesotaxis, and stochastic cell–cell adhesion dynamics that allow for intercalation behaviors to emerge in the model.

This model allows us to investigate the respective roles of cohesotaxis, circular in vivo tissue geometry, and intercalation in driving tissue closure in the mesendoderm mantle. The migration speed of a single DMZ explant was sensitive to the directional bias provided by the computational representation of cohesotaxis that we developed and implemented in this work. At the demonstrated ranges, migration rates in our model were sensitive to the stiffness of cell–cell adhesion, leading edge cell–substrate adhesion and lamellipodial extension distance, but robust to cell–substrate link reformation rates and the stiffness of follower cell–substrate adhesion. Future work may explore these parameters over greater ranges. Cohesotaxis allowed cells to migrate synchronously, increasing the measured speed of migration. The method with which we encode cohesotaxis in our agent-based model is directionally agnostic, meaning it works no matter which direction the explant is migrating. This allowed us to perform the ITR experiments in silico, as four different explants would be migrating towards each other. Additionally, to our knowledge, this method of mathematically encoding cohesotaxis behavior has not yet been described in the literature. When the model was placed ITR, it exhibited migration speed that exceeded the migration speed of a single DMZ explant model even when synchronous migration through cohesotaxis was added. Additionally, this finding was true for all tested values of the cohesotaxis parameter. This suggests that the circular geometry of the tissue also allows for synchronous cell migration to occur, likely because it restricts the area towards which leader cells may migrate, resulting in cooperative cell migration. Finally, adding additional rules to allow for cell rearrangement behaviors in the explant resulted in a further increase in migration speed as well as the observation that cells intercalated radially (i.e. in the direction of the substrate). This radial intercalation is a model prediction that was validated biologically with the introduction of single labeled mesendoderm cells onto the top of an unlabeled DMZ explant, as well as in control explants. An increase in migration speed in response to intercalation was not reported in earlier studies (Davidson et al., 2002).

Additional validation of model predictions was carried out experimentally by lowering extracellular Ca++ concentration in order to reduce cadherin adhesion and promote radial intercalation. Migration rates are increased in DMZ explants under these conditions as predicted by the model. This change in rate is consistent with an increase in top-down (radial) intercalation brought about by a shift in the balance between cell–cell and cell–matrix adhesion in 0.8 mM Ca++. The likely consequence of additional cells moving to substrate level is more rapid spreading of the explant, thus contributing to the increased rate of migration observed and predicted by the model. Although Ca++ levels in these experiments were reduced by only 20% (standard DFA media contains 1.0 mM Ca++), even modest reductions in Ca++ concentration might affect cellular processes in addition to cadherin-based cell–cell adhesion. However, the 0.8 mM Ca++ used in these experiments is well within the physiological range for interstitial Ca++ levels at gastrulation in Xenopus. Gillespie, (1983) reported that blastocoel fluid contains 0.5 mM Ca++ at gastrula stages when multiple tissues including the mesendoderm begin to move. Interestingly, interstitial Ca++ concentrations exceed 1.0 mM prior to and after gastrulation (Gillespie, 1983), suggesting that a developmentally regulated drop in [Ca++] may trigger the onset of gastrulation movements by ‘loosening-up’ cells to promote cellular rearrangement and tissue morphogenesis, as first suggested by Holtfreter, (1944). We and others have proposed previously that there is a decrease in cell–cell adhesion at gastrulation (Brieher and Gumbiner, 1994) coincident with an increase in cell-ECM adhesion (e.g. Ramos et al., 1996). Regulation of these stage-dependent differences in Ca++ concentration may provide a mechanism to promote the dramatic cell and tissue movements that occur at gastrulation. These observations contribute important context for our CPM model and its experimental validation. Moreover, while we altered intercellular adhesion through changes in cell–cell link parameters to produce intercalation in this study, future work could explore replicating intercalation through altered interfacial tension between discrete cells in the CPM framework.

While our model is designed to provide insight into mechanisms driving collective cell migration during mesendoderm mantle closure in Xenopus gastrulation, the cell agent behaviors are not specific or particular to this system. Therefore, our model may be extended to other systems known to exhibit collective cell migration, such as cell crawling mechanisms involved in neural crest migration or epithelial gap closure for wound healing or cancer metastasis (Shellard et al., 2018; Spatarelu et al., 2019; Staddon et al., 2018). Our model is specified in the accessible open-source CompuCell3D modeling framework and the behaviors of cohesotaxis and intercalation are specified in Python, allowing users to modify, extend, or re-use model mechanisms in a modular fashion (Sego et al., 2020). Many computational models of cell migration focus on single-cell migration in detail (Fortuna et al., 2020; Kumar et al., 2018; Scianna et al., 2013) or collective migration in two-dimensional epithelial cell monolayers (Harrison et al., 2011; Khataee et al., 2020; Nguyen Edalgo et al., 2019; Pan et al., 2021; Staddon et al., 2018; Tetley et al., 2019; Zhao et al., 2017). However, agent-based models built to understand the emergent role of complex three-dimensional tissue geometries have yet to be widely employed. Three-dimensional models are better suited to studying tissue morphogenesis because they can capture the mechanics and behaviors of multiple layers of cells.

Our initial model recapitulates many general behaviors prevalent in collectively migrating mesendoderm tissue. At present, it offers limited insight into environmental drivers of migration, such as chemotaxis through PDGF or other biochemical gradients, haptotaxis, or durotaxis. Additionally, the simulated cells in the model's current formulation have limited capacity to deform and are unable to demonstrate dramatic cell shape elongation known to be present in mesendoderm mantle closure (Sonavane et al., 2017; Weber et al., 2012). These limitations will be addressed in future model iterations by applying the agent behaviors developed in this work to other modeling methodologies, such as vertex model or finite element methods which are able to represent forces and stresses more explicitly present in the tissue during morphogenesis. In addition, the current CPM can be extended to capture the more representative geometry of the toroid explant, which has been investigated in prior experimental work (Davidson et al., 2002; Sonavane et al., 2017). Further interrogation of this model in the context of representative in vivo tissue geometry will be pursued in future work. The CPM model adapted to represent more complex geometry would allow for the testing of additional hypotheses for mesendoderm mantle closure in vivo, such as the effects of tissue surface tension throughout the mesendoderm (Shook et al., 2022).

In conclusion, we have developed and interrogated a model of collectively migrating cells in a tissue to elucidate the roles of cohesotaxis, circular geometry, and intercalation in modulating tissue migration speed. We describe how cohesotaxis is an important mechanism for enabling cell migration, however the circular geometry of the tissue results in greater migration speed of the explant relative to cohesotaxis. Thus, while cohesotaxis is important for the polarized protrusive behavior of migratory cells, it does not affect the speed of migration in situations of constrained tissue geometry such as ITR. Top-down (radial) cell-intercalation behavior provides an additional mechanism contributing to migrating mesendoderm but is not required for mantle closure. Our model predicts previously undescribed contributions of cell intercalation to collective cell migration. Three-dimensional computational models such as the one developed in this study represent a promising tool for the thorough investigation of interdependent mechanisms of tissue morphogenesis. This work has described the development, validation, and application of a novel Cellular–Potts model (CPM) of collective cell migration to investigate the respective roles of cohesotaxis, cell intercalation, and tissue geometry during tissue morphogenesis of the Xenopus laevis gastrula.

Poisson parameter ζ was calibrated for leader and follower cell agents to allow the model to reflect observed mesendoderm cell–substrate attachment and single-cell movement behavior replicated from video data collected for the development of the model. Representative video of a live DMZ explant with leader and follower cells is included in Movie 1. Parameter values selected to reproduce observed biological explant behavior are displayed in Table 1. Units for CPM parameters are given when possible, however exact units for certain parameters in the CPM methodology are often difficult to determine given the method's reliance on statistical minimization of the abstract effective energy functional to reproduce stochastic cellular behaviors.

The CPM temperature parameter affects the likelihood of accepting voxel copy attempt if system effective energy is increased, whereas voxel copy attempts that decrease system effective energy are always accepted. is the change in system effective energy after a voxel copy attempt.

The geometry of a biological DMZ explant was approximated by simulating a tissue eight cells in width, four cells in length, and two cells in height where each cell occupies 5×5×5 voxels to represent cells approximately 30×30×30 μm in size. The explant is placed upon a flat substrate for the simulated cells to form mechanical link attachments in order to represent cell–substrate adhesion.

During model initialization and generation of explant geometry, eight cells at an assumed leading edge of the explant model are assigned to be leader cell agents and the remaining cells are assigned to be follower cell agents. Cells do not change type during the simulation. Leader cell agents extend mechanical link objects of constant tension from their cell centroids out to the substrate in the direction of a randomly chosen voxel along the cell border with its medium to establish an initial forward direction for the explant to begin migration. Follower cells bordering the substrate similarly extend Hookean spring-like mechanical link objects to their substrate to simulate their cell–substrate adhesion behaviors. All cells establish Hookean spring-like mechanical links to their neighbors at the start of the simulation to recreate cell–cell adhesion.

Each timestep, leader cell agents delete and form new cell–substrate links to create lamellipodial extension and retraction behavior with a frequency governed by their defined probability (Eqn. 3). New links are established in a direction decided by sampling randomly from a location along the leader cell agent bordering the medium, that is, the free edge of a leader cell agent. Should the leader cell be surrounded by other cells such that it does not have a free edge as can be seen in representative screenshots of the ITR experiments (Figs 4A and 5A), leader cells do not extend links representative of lamellipodia, thus behaving effectively as follower cells. Follower cell agent cell–substrate adhesion similarly deletes and re-forms according to its defined probability during each computational timestep. A computational timestep was selected to represent 5 s of real-world time. Representative video of a DMZ explant migration along a simulated substrate is demonstrated in Movie 2.

To validate the model, the CPM model was applied to reproduce an experiment previously performed in Davidson et al., 2002. The model was simulated to migrate for approximately 17 min to reach a steady state of migratory behavior before stopping all further leading edge lamellipodial link objects from being formed. The remaining links at this timepoint and beyond were allowed to disappear as a function of their probabilistic lifetime in the model, which resulted in a retraction event comparable to an experiment introducing a function-blocking monocolonal antibody (mAB) to fibronectin preventing any further lamellipodial connections (Fig. 6A-D). The baseline or reference parameterization that was chosen initially to match migratory behavior from movies of the biological explant was able to reproduce comparable mean retraction distances to those reported in Davidson et al., 2002 from n=10 simulations (Fig. 6E). This demonstrates that the chosen reference parameterization and resulting model agent-based behaviors are able to reproduce a similar directional tension known to be present in the biological explant and that these mechanical properties produce model events on a spatial scale that can be compared to experimental conditions.

Cohesotaxis is a mechanism by which adherent cells coordinate mechanically to establish a persistent direction of migratory behavior (Weber et al., 2012). We have encoded a bias to the direction to which leader cell agents in our computational model migrate to represent cohesotaxis along with a tunable parameter to modulate the extent to which this bias is applied. To allow for investigation of the ITR geometry, our method of simulating cohesotaxis must be agnostic to the direction the tissue faces. To implement this mechanism, when leader cell agents form new cell–substrate mechanical links to represent lamellipodial cycling behavior during model evaluation, leader cell agents must select preferentially from free-edge voxels in the direction most opposed to sites of cell–cell adhesion rather than selecting a random voxel from the list of eligible free-edge voxels. Free-edge voxels are defined as voxels belonging to a cell agent that exist in contact with medium as opposed to in contact with other cells. Therefore, during formation of a new mechanical link object, the list of eligible free-edge voxels is sorted by cumulative distance to all cell border voxels. The free edge voxel through which a new link's direction is established is selected preferentially from those in the direction of migration, or most opposite from cell–cell contacts. This behavior was created in the model by selecting with greater weight given to voxels with the lowest cumulative distance to the list of cell border voxels. The voxel is selected by the probability defined in Fig. 7B, where cohesotaxis parameter κ approximates a kurtosis of the probability distribution function to tune the extent to which the migratory bias is applied during the simulation.

Values of f(x) for each function were normalized by the sum of the function within the viewing window to convert the vector f(x) into a vector of probability weights that sum to 1 such that a list of an arbitrary number of sorted voxel objects could be selected using the vector of probabilities and subsequently used to encode agent-based behaviors in the CompuCell3D simulation software. Normalized probability functions for different values of κ are shown in Fig. 7B. The direction of a link object representing lamellipodia is defined by the selecting free-edge voxels through this algorithm, whereas the distance away from the selected voxel to which a link object is created is defined by parameter χ­­_Lamellipodia.

We simulated the ITR configuration by placing four DMZ explants orthogonal to one another in the same simulation space. We did not change other parameters intrinsic to the model or CPM framework to recreate the ITR configuration, and allowed the single DMZ explant model behaviors drive the ITR simulation.

Cell rearrangement behaviors that result in cell intercalation were enabled in the model by deleting and re-forming mechanical link objects that represent cell–cell adhesion within the tissue. When mechanical link objects are removed and re-formed, the CPM Potts algorithm causes cells to shift and rearrange in the tissue encouraged by link formation to new cell neighbors. Link objects were defined to have a lifetime modeled as a Poisson process as described by Eqn 3, where Poisson parameter ζ_Tissue represents a parameter to allow for variable frequency to which cell–cell adhesion links break and re-form. During model simulation, if a cell–cell adhesion link has been determined to break based on evaluating this probability, then the cell will form a new cell–cell adhesion link with a random neighbor to which it does not currently have a cell–cell adhesion link until a maximum of four neighbors is reached. This allows cell intercalation behaviors to emerge within the model. The maximum number of neighbors remains a tunable parameter in the model code, however we have observed that parameter values of 3 or lower, or 6 and greater, result in unrealistic tissue behavior. This value allows the tissue to be fluid enough to rearrange but connected enough to prevent it from falling apart at a range of ζ_Tissue values.

Single DMZ in silico experiments were performed by allowing the DMZ model to run for a period of 2 h while sampling leading edge centroid locations approximately every 17 mins (every 200 timesteps). The average distance traversed by the center four leading edge cells of the explant were then calculated. Migration speed was determined to be the average speed of these cells of the explant in the direction of migration over the 2 h.

Closure speed in the ITR experiments was determined by measuring the free area towards which the explants migrated during simulation of closure of the mesendoderm mantle. The area of the cell-free space in the center of the ITR configuration was measured and the rate of change of the length of one side (⁠⁠) was used to determine the migration speed of the tissue.

Gametes were obtained from adult Xenopus laevis, eggs fertilized and embryos reared to gastrula stages using standard methods (Sive et al., 2000). Animal housing and care protocols were approved by the Animal Care and Use Committee (ACUC) at the University of Virginia.

DMZs were dissected from stage 10.5 embryos and placed on glass bottomed dishes coated overnight with 10 μg/ml Fibronectin (Sigma F1141) at 4°C. For ITR measurements four DMZs were arranged orthogonally as described in Davidson et al., 2002. Images were acquired every 5 min for 3 h. The distance migrated by the free edge over 40 min was calculated before and after adjacent DMZs contacted one another. Single DMZs were allowed to migrate for at least 1 h then imaged for 30 min. The average distance migrated by four leading-edge cells per DMZ was calculated.

DMZs ITR were fixed for 10 min with 3.7% formaldehyde, 0.1% glutaraldehyde in 20 mMTris 150 mM NaCL pH7.4 containing 0.1% Tween20. The explants were then incubated overnight with a 1:100 dilution of Alexa488 Acti-stain (Cytoskeleton PHDG1) and imaged using a Nikon AX confocal microscope and a 10× objective lens.

DMZ explants were dissected from stage 11 embryos and placed on glass bottomed dishes coated with FN. Mesendoderm tissue was dissected from embryos that were injected with Alexa488 dextran and dissociated in calcium and magnesium free MBS to obtain single cells. After 30 min of migration dissociated cells were ‘sprinkled’ on the migrating DMZ explants. Beginning 90 mins after the cells were added, confocal z stacks (13×2 μm sections) were obtained at 1 min intervals with a Nikon AX confocal microscope.

For data displayed in Table S1, ME was explanted from embryos at stage 10.5 to 11.5 and cultured in DFA (Sater et al., 1993) or DFA made with 0.8 mM Ca++ and Mg++ (i.e. 20% less Ca++ and Mg++ than standard DFA with 1.0 mM Ca++ and Mg++). Petri dish surfaces were coated with 10 or 20 mg/ml FN for 30' at 37°C.

Images were captured every 2 min on a Zeiss AxioZoom dissecting microscope. One to three points across the leading edge of the migrating ME were tracked every 10 mins using the ObjectJ plugin in ImageJ (Schindelin et al., 2012). Displacement of each point over time for the first 1 and 2 h of migration was determined, and the point with maximum velocity was used for the analysis of each explant. Explants which showed a great deal of slippage (tendency to snap back) or failure to migrate over the first 2 h were not included in data analysis. Migration rates were subjected to two-way analysis of variance (ANOVA), to account for differences in rate between clutches, differences in stages, and FN amount.

For pairwise comparisons, unpaired, two-tailed, Student's t-tests were used to determine P values. For multiple-group comparisons, we used one-way ANOVA to determine whether multiple groups had a common mean. We did not perform post-hoc testing following an ANOVA because the presence or absence of a difference in the means of the multiple groups were sufficient for our analysis. P values were as follows: P>0.05 (not significant, ns), *P<0.05, ***P<0.001. All error bars represent standard deviation (s.d.) except Fig. 6, which shows standard error of the mean (s.e.m) to remain consistent with data reported in Davidson et al., 2002.

We acknowledge Anita Impagliazzo for graphical illustrations. Fig. 6 panels A and B reprinted Davidson, L. A., Hoffstrom, B. G., Keller, R. and DeSimone, D. W. (2002). Mesendoderm extension and mantle closure in Xenopus laevis gastrulation: combined roles for integrin α5β1, fibronectin, and tissue geometry. Dev. Biol. 242, 109-129, p. 117, with permission from Elsevier.